2 !-----------------------------------------------------------------------------
5 use control_data, only: tor_mode
7 use geometry_data, only:nres,boxxsize,boxysize,boxzsize
9 ! use control_data, only:maxthetyp1
10 use energy, only:etotal,enerprint,rescale_weights
14 ! include "COMMON.MPI"
17 !-----------------------------------------------------------------------------
20 real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW
21 real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm)
22 real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,&
23 aksc_all,abond0_all !(maxbondterm,ntyp,max_parm)
24 real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm)
25 real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,&
26 bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
27 real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm)
28 real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm)
29 real(kind=8),dimension(:,:),allocatable :: theta0_all,&
30 sig0_all,sigc0_all !(-ntyp:ntyp,max_parm)
31 real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all
32 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
33 real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all
34 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
35 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,&
36 ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
37 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
38 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,&
39 ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3,
40 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
41 real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm)
42 real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm)
43 real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm)
44 real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm)
45 real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm)
46 real(kind=8),dimension(:,:,:,:),allocatable :: v0_all
47 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
48 real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,&
49 v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
50 real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,&
51 vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm)
52 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,&
53 v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
54 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all
55 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
56 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all
57 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
58 real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm)
59 real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,&
60 ee_all !(2,2,-maxtor:maxtor,max_parm)
61 real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,&
62 dtilde_all !(2,2,-maxtor:maxtor,max_parm)
63 real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm)
64 real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,&
65 ael6_all,ael3_all !(2,2,max_parm)
66 real(kind=8),dimension(:,:,:),allocatable :: aad_all,&
67 bad_all !(ntyp,2,max_parm)
68 real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,&
69 augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm)
70 real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm)
71 real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,&
72 akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm)
73 real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,&
74 v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
75 integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm)
76 integer,dimension(:,:,:,:),allocatable :: nlor_all,&
77 nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
78 integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,&
79 ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
80 integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm)
81 integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm)
82 integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,&
83 ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,&
84 nntheterm_all !(max_parm)
85 integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
86 !-----------------------------------------------------------------------------
89 !-----------------------------------------------------------------------------
91 !-----------------------------------------------------------------------------
92 subroutine enecalc(islice,*)
95 use control_data, only:indpdb
96 use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,&
97 vbld,rad2deg,dc_norm,dc,vbld_inv
98 use io_base, only:gyrate!,briefout
99 use geometry, only:int_from_cart1,returnbox
100 use io_wham, only:pdboutW
101 use io_database, only:opentmp
102 use conform_compar, only:qwolynes,rmsnat
103 ! include "DIMENSIONS"
104 ! include "DIMENSIONS.ZSCOPT"
105 ! include "DIMENSIONS.FREE"
109 ! include "COMMON.MPI"
111 ! include "COMMON.CHAIN"
112 ! include "COMMON.IOUNITS"
113 ! include "COMMON.PROTFILES"
114 ! include "COMMON.NAMES"
115 ! include "COMMON.VAR"
116 ! include "COMMON.SBRIDGE"
117 ! include "COMMON.GEO"
118 ! include "COMMON.FFIELD"
119 ! include "COMMON.ENEPS"
120 ! include "COMMON.LOCAL"
121 ! include "COMMON.WEIGHTS"
122 ! include "COMMON.INTERACT"
123 ! include "COMMON.FREE"
124 ! include "COMMON.ENERGIES"
125 ! include "COMMON.CONTROL"
126 ! include "COMMON.TORCNSTR"
129 integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
131 character(len=64) :: nazwa
132 character(len=80) :: bxname
133 character(len=3) :: liczba
134 !el real(kind=8) :: qwolynes
135 !el external qwolynes
136 integer :: errmsg_count,maxerrmsg_count=100
137 !el real(kind=8) :: rmsnat,gyrate
138 !el external rmsnat,gyrate
139 real(kind=8) :: tole=1.0d-1
140 integer i,itj,ii,iii,j,k,l,licz
141 integer ir,ib,ipar,iparm
143 real(kind=4) :: csingle(3,nres*2)
144 real(kind=8) :: energ
146 !el integer ilen,iroof
147 !el external ilen,iroof
148 real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini
149 !el real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini
150 real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2
152 integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm)
154 character(len=64) :: bprotfile_temp
157 integer,dimension(0:nprocs) :: scount_
158 !el real(kind=8) :: rmsnat
160 rescale_mode=rescale_modeW
162 call opentmp(islice,ientout,bprotfile_temp)
168 write (iout,*) "enecalc: nparmset ",nparmset
177 do i=indstart(me1),indend(me1)
178 write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1)
188 write(iout,*)"enecalc_ i ntot",i,ntot
190 read(ientout,rec=i,err=101) &
191 ((csingle(l,k),l=1,3),k=1,nres),&
192 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
193 nss,(ihpb(k),jhpb(k),k=1,nss),&
194 eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
196 !write(iout,*)"co wczytuje"
197 ! write(iout,*)((csingle(l,k),l=1,3),k=1,nres),&
198 ! ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
199 ! nss,(ihpb(k),jhpb(k),k=1,nss),&
200 ! eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
203 !write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3)
204 if (indpdb.gt.0) then
212 c(l,k+nres)=csingle(l,k+nres)
223 csingle(l,k+nres)=c(l,k+nres)
226 anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
227 q(nQ+1,iii+1)=rmsnat(iii+1)
229 q(nQ+2,iii+1)=gyrate(iii+1)
230 ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
231 ! fT=T0*beta_h(ib,ipar)*1.987D-3
232 ! ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
233 ! EL start old rescale
234 ! if (rescale_modeW.eq.1) then
235 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
237 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
238 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
239 !#elif defined(FUNCT)
249 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
251 ! else if (rescale_modeW.eq.2) then
252 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
254 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
255 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
256 !#elif defined(FUNCT)
264 ! fT(l)=1.12692801104297249644d0/ &
265 ! dlog(dexp(quotl)+dexp(-quotl))
267 ! else if (rescale_modeW.eq.0) then
272 ! write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
278 ! write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
279 ! & " kfac",kfac,"quot",quot," fT",fT
281 write(iout,*)"weights"
282 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
283 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
292 call int_from_cart1(.false.)
295 ! call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
298 write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
299 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
300 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
303 call restore_parm(iparm)
304 call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
306 write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
307 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
308 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
311 ! call etotal(energia(0),fT)
312 call etotal(energia(0))
313 !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
315 !write(iout,*)k,"c=",(c(l,k),l=1,3)
316 !write(iout,*)k,"dc=",(dc(l,k),l=1,3)
317 !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
320 !write(iout,*)k,"vbld=",vbld(k)
321 !write(iout,*)k,"vbld_inv=",vbld_inv(k)
324 !write(iout,*)"energia",(energia(j),j=0,n_ene)
325 !write(iout,*)"enerprint tuz po call etotal"
326 call enerprint(energia(0))
328 write (iout,*) "Conformation",i
329 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
330 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
331 ! call enerprint(energia(0),fT)
332 call enerprint(energia(0))
333 write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
334 write (iout,*) "ftors",ftors
335 !el call briefout(i,energia(0))
336 temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
337 write (iout,*) "temp", temp
338 call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
340 if (energia(0).ge.1.0d20) then
341 write (iout,*) "NaNs detected in some of the energy",&
342 " components for conformation",ii+1
343 write (iout,*) "The Cartesian geometry is:"
344 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
345 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
346 write (iout,*) "The internal geometry is:"
348 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
349 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
350 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
351 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
352 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
353 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
354 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
355 write (iout,*) "The components of the energy are:"
356 ! call enerprint(energia(0),fT)
357 call enerprint(energia(0))
359 "This conformation WILL NOT be added to the database."
364 if (ipar.eq.iparm) write (iout,*) i,iparm,&
365 1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
367 if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
368 dabs(eini-energia(0)).gt.tole) then
369 if (errmsg_count.le.maxerrmsg_count) then
370 write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
371 "Warning: energy differs remarkably from ",&
372 " the value read in: ",energia(0),eini," point",&
373 iii+1,indstart(me1)+iii," T",&
374 1.0d0/(1.987D-3*beta_h(ib,ipar))
376 call pdboutW(indstart(me1)+iii,&
377 1.0d0/(1.987D-3*beta_h(ib,ipar)),&
378 energia(0),eini,0.0d0,0.0d0)
379 write (iout,*) "The Cartesian geometry is:"
380 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
381 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
382 write (iout,*) "The internal geometry is:"
384 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
385 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
386 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
387 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
388 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
389 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
390 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
391 call enerprint(energia(0))
392 ! call enerprint(energia(0),fT)
393 errmsg_count=errmsg_count+1
394 if (errmsg_count.gt.maxerrmsg_count) &
395 write (iout,*) "Too many warning messages"
396 if (einicheck.gt.1) then
397 write (iout,*) "Calculation stopped."
400 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
407 potE(iii+1,iparm)=energia(0)
409 enetb(k,iii+1,iparm)=energia(k)
411 ! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
412 ! write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
414 write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
415 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
416 ! call enerprint(energia(0),fT)
419 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
420 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
421 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
422 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
423 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
424 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
425 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
426 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
427 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
428 write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
429 write (iout,'(f10.5,i10)') rmsdev,iscor
430 ! call enerprint(energia(0),fT)
431 call enerprint(energia(0))
432 call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
439 if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
440 write (ientout,rec=iii) &
441 ((csingle(l,k),l=1,3),k=1,nres),&
442 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
443 nss,(ihpb(k),jhpb(k),k=1,nss),&
444 potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
445 ! write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
447 if (separate_parset) then
448 snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
450 snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
452 ! write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
453 ! & " snk",snk_p(iR,ib,ipar)
455 snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
461 write (iout,*) "Me",me," scount",scount(me)
463 ! Master gathers updated numbers of conformations written by all procs.
464 call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
465 MPI_INTEGER, WHAM_COMM, IERROR)
469 indstart(i)=indend(i-1)+1
470 indend(i)=indstart(i)+scount_(i)-1
473 write (iout,*) "Revised conformation counts"
475 write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
476 "Processor",i," indstart",indstart(i),&
477 " indend",indend(i)," count",scount_(i)
480 call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
481 MaxR*MaxT_h*nParmSet,&
482 MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
488 stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
492 write (iout,*) "Revised SNK"
495 write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
496 iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
497 (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
498 write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
501 write (iout,'("Total",i10)') stot(islice)
507 101 write (iout,*) "Error in scratchfile."
511 end subroutine enecalc
512 !------------------------------------------------------------------------------
513 logical function conf_check(ii,iprint)
515 use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
516 use geometry, only:int_from_cart1
517 ! include "DIMENSIONS"
518 ! include "DIMENSIONS.ZSCOPT"
519 ! include "DIMENSIONS.FREE"
523 ! include "COMMON.MPI"
525 ! include "COMMON.CHAIN"
526 ! include "COMMON.IOUNITS"
527 ! include "COMMON.PROTFILES"
528 ! include "COMMON.NAMES"
529 ! include "COMMON.VAR"
530 ! include "COMMON.SBRIDGE"
531 ! include "COMMON.GEO"
532 ! include "COMMON.FFIELD"
533 ! include "COMMON.ENEPS"
534 ! include "COMMON.LOCAL"
535 ! include "COMMON.WEIGHTS"
536 ! include "COMMON.INTERACT"
537 ! include "COMMON.FREE"
538 ! include "COMMON.ENERGIES"
539 ! include "COMMON.CONTROL"
540 ! include "COMMON.TORCNSTR"
543 integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
545 integer :: j,k,l,ii,itj,iprint,mnum,mnum1
546 if (.not.check_conf) then
550 call int_from_cart1(.false.)
555 if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
556 (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
558 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
559 " for conformation",ii,mnum
560 if (iprint.gt.1) then
561 write (iout,*) "The Cartesian geometry is:"
562 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
563 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
564 write (iout,*) "The internal geometry is:"
565 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
566 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
567 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
568 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
569 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
570 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
572 if (iprint.gt.0) write (iout,*) &
573 "This conformation WILL NOT be added to the database."
582 if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
583 (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
585 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
586 " for conformation",ii
587 if (iprint.gt.1) then
588 write (iout,*) "The Cartesian geometry is:"
589 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
590 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
591 write (iout,*) "The internal geometry is:"
592 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
593 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
594 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
595 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
596 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
597 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
599 if (iprint.gt.0) write (iout,*) &
600 "This conformation WILL NOT be added to the database."
606 if (theta(j).le.0.0d0) then
608 write (iout,*) "Zero theta angle(s) in conformation",ii
609 if (iprint.gt.1) then
610 write (iout,*) "The Cartesian geometry is:"
611 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
612 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
613 write (iout,*) "The internal geometry is:"
614 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
615 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
616 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
617 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
618 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
619 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
621 if (iprint.gt.0) write (iout,*) &
622 "This conformation WILL NOT be added to the database."
626 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
629 ! write (iout,*) "conf_check passed",ii
631 end function conf_check
632 !-----------------------------------------------------------------------------
634 !-----------------------------------------------------------------------------
635 subroutine store_parm(iparm)
637 ! Store parameters of set IPARM
638 ! valence angles and the side chains and energy parameters.
641 ! include 'DIMENSIONS'
642 ! include 'DIMENSIONS.ZSCOPT'
643 ! include 'DIMENSIONS.FREE'
644 ! include 'COMMON.IOUNITS'
645 ! include 'COMMON.CHAIN'
646 ! include 'COMMON.INTERACT'
647 ! include 'COMMON.GEO'
648 ! include 'COMMON.LOCAL'
649 ! include 'COMMON.TORSION'
650 ! include 'COMMON.FFIELD'
651 ! include 'COMMON.NAMES'
652 ! include 'COMMON.SBRIDGE'
653 ! include 'COMMON.SCROT'
654 ! include 'COMMON.SCCOR'
655 ! include 'COMMON.ALLPARM'
656 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
658 call alloc_enecalc_arrays(iparm)
659 !el allocate(ww_all(n_ene,iparm))
663 ww_all(3,iparm)=welec
664 ww_all(4,iparm)=wcorr
665 ww_all(5,iparm)=wcorr5
666 ww_all(6,iparm)=wcorr6
667 ww_all(7,iparm)=wel_loc
668 ww_all(8,iparm)=wturn3
669 ww_all(9,iparm)=wturn4
670 ww_all(10,iparm)=wturn6
671 ww_all(11,iparm)=wang
672 ww_all(12,iparm)=wscloc
673 ww_all(13,iparm)=wtor
674 ww_all(14,iparm)=wtor_d
675 ww_all(15,iparm)=wstrain
676 ww_all(16,iparm)=wvdwpp
677 ww_all(17,iparm)=wbond
678 ww_all(19,iparm)=wsccor
679 ww_all(42,iparm)=wcatprot
680 ww_all(41,iparm)=wcatcat
681 ! Store bond parameters
682 vbldp0_all(iparm)=vbldp0
685 nbondterm_all(i,iparm)=nbondterm(i)
687 vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
688 aksc_all(j,i,iparm)=aksc(j,i)
689 abond0_all(j,i,iparm)=abond0(j,i)
692 ! Store bond angle parameters
695 a0thet_all(i,iparm)=a0thet(i)
699 athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
700 bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
705 polthet_all(j,i,iparm)=polthet(j,i)
708 gthet_all(j,i,iparm)=gthet(j,i)
710 theta0_all(i,iparm)=theta0(i)
711 sig0_all(i,iparm)=sig0(i)
712 sigc0_all(i,iparm)=sigc0(i)
715 IF (tor_mode.eq.0) THEN
716 nthetyp_all(iparm)=nthetyp
717 ntheterm_all(iparm)=ntheterm
718 ntheterm2_all(iparm)=ntheterm2
719 ntheterm3_all(iparm)=ntheterm3
720 nsingle_all(iparm)=nsingle
721 ndouble_all(iparm)=ndouble
722 nntheterm_all(iparm)=nntheterm
724 ithetyp_all(i,iparm)=ithetyp(i)
727 do i=-nthetyp-1,nthetyp+1
728 do j=-nthetyp-1,nthetyp+1
729 do k=-nthetyp-1,nthetyp+1
730 aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
732 aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
736 bbthet_all(m,l,i,j,k,iblock,iparm)= &
737 bbthet(m,l,i,j,k,iblock)
738 ccthet_all(m,l,i,j,k,iblock,iparm)= &
739 ccthet(m,l,i,j,k,iblock)
740 ddthet_all(m,l,i,j,k,iblock,iparm)= &
741 ddthet(m,l,i,j,k,iblock)
742 eethet_all(m,l,i,j,k,iblock,iparm)= &
743 eethet(m,l,i,j,k,iblock)
749 if (iblock.eq.1) then
750 ffthet_all1(mm,m,l,i,j,k,iparm)=&
751 ffthet(mm,m,l,i,j,k,iblock)
752 ggthet_all1(mm,m,l,i,j,k,iparm)=&
753 ggthet(mm,m,l,i,j,k,iblock)
755 ffthet_all2(mm,m,l,i,j,k,iparm)=&
756 ffthet(mm,m,l,i,j,k,iblock)
757 ggthet_all2(mm,m,l,i,j,k,iparm)=&
758 ggthet(mm,m,l,i,j,k,iblock)
768 write(iout,*) "Need storing parameters"
772 ! Store the sidechain rotamer parameters
775 !! write (iout,*) i,"storeparm1"
777 nlob_all(iii,iparm)=nlob(iii)
779 bsc_all(j,iii,iparm)=bsc(j,iii)
781 censc_all(k,j,i,iparm)=censc(k,j,i)
785 gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
793 sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
797 IF (tor_mode.eq.0) THEN
798 ! Store the torsional parameters
800 do i=-ntortyp+1,ntortyp-1
801 do j=-ntortyp+1,ntortyp-1
802 v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
803 nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
804 nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
805 do k=1,nterm(i,j,iblock)
806 v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
807 v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
809 do k=1,nlor(i,j,iblock)
810 vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
811 vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
812 vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
817 ! Store the double torsional parameters
819 do i=-ntortyp+1,ntortyp-1
820 do j=-ntortyp+1,ntortyp-1
821 do k=-ntortyp+1,ntortyp-1
822 ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
823 ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
824 do l=1,ntermd_1(i,j,k,iblock)
825 v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
826 v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
827 v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
828 v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
830 do l=1,ntermd_2(i,j,k,iblock)
831 do m=1,ntermd_2(i,j,k,iblock)
832 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
839 ! Store parameters of the cumulants
842 b1_all(j,i,iparm)=b1(j,i)
843 b1tilde_all(j,i,iparm)=b1tilde(j,i)
844 b2_all(j,i,iparm)=b2(j,i)
848 cc_all(k,j,i,iparm)=cc(k,j,i)
849 ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
850 dd_all(k,j,i,iparm)=dd(k,j,i)
851 dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
852 ee_all(k,j,i,iparm)=ee(k,j,i)
857 write(iout,*) "NEED storing parameters"
859 ! Store the parameters of electrostatic interactions
862 app_all(j,i,iparm)=app(j,i)
863 bpp_all(j,i,iparm)=bpp(j,i)
864 ael6_all(j,i,iparm)=ael6(j,i)
865 ael3_all(j,i,iparm)=ael3(j,i)
868 ! Store sidechain parameters
871 aa_all(j,i,iparm)=aa_aq(j,i)
872 bb_all(j,i,iparm)=bb_aq(j,i)
873 r0_all(j,i,iparm)=r0(j,i)
874 sigma_all(j,i,iparm)=sigma(j,i)
875 chi_all(j,i,iparm)=chi(j,i)
876 augm_all(j,i,iparm)=augm(j,i)
877 eps_all(j,i,iparm)=eps(j,i)
881 chip_all(i,iparm)=chip(i)
882 alp_all(i,iparm)=alp(i)
884 ! Store the SCp parameters
887 aad_all(i,j,iparm)=aad(i,j)
888 bad_all(i,j,iparm)=bad(i,j)
891 ! Store disulfide-bond parameters
900 ! Store SC-backbone correlation parameters
901 do i=-nsccortyp,nsccortyp
902 do j=-nsccortyp,nsccortyp
904 nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
908 do k=1,nterm_sccor(j,i)
909 v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
910 v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
915 write(iout,*)"end of store_parm"
917 end subroutine store_parm
918 !--------------------------------------------------------------------------
919 subroutine restore_parm(iparm)
921 ! Store parameters of set IPARM
922 ! valence angles and the side chains and energy parameters.
925 ! include 'DIMENSIONS'
926 ! include 'DIMENSIONS.ZSCOPT'
927 ! include 'DIMENSIONS.FREE'
928 ! include 'COMMON.IOUNITS'
929 ! include 'COMMON.CHAIN'
930 ! include 'COMMON.INTERACT'
931 ! include 'COMMON.GEO'
932 ! include 'COMMON.LOCAL'
933 ! include 'COMMON.TORSION'
934 ! include 'COMMON.FFIELD'
935 ! include 'COMMON.NAMES'
936 ! include 'COMMON.SBRIDGE'
937 ! include 'COMMON.SCROT'
938 ! include 'COMMON.SCCOR'
939 ! include 'COMMON.ALLPARM'
940 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
945 welec=ww_all(3,iparm)
946 wcorr=ww_all(4,iparm)
947 wcorr5=ww_all(5,iparm)
948 wcorr6=ww_all(6,iparm)
949 wel_loc=ww_all(7,iparm)
950 wturn3=ww_all(8,iparm)
951 wturn4=ww_all(9,iparm)
952 wturn6=ww_all(10,iparm)
953 wang=ww_all(11,iparm)
954 wscloc=ww_all(12,iparm)
955 wtor=ww_all(13,iparm)
956 wtor_d=ww_all(14,iparm)
957 wstrain=ww_all(15,iparm)
958 wvdwpp=ww_all(16,iparm)
959 wbond=ww_all(17,iparm)
960 wsccor=ww_all(19,iparm)
961 wcatprot=ww_all(42,iparm)
962 wcatcat=ww_all(41,iparm)
963 ! Restore bond parameters
964 vbldp0=vbldp0_all(iparm)
967 nbondterm(i)=nbondterm_all(i,iparm)
969 vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
970 aksc(j,i)=aksc_all(j,i,iparm)
971 abond0(j,i)=abond0_all(j,i,iparm)
974 ! Restore bond angle parameters
977 a0thet(i)=a0thet_all(i,iparm)
981 athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
982 bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
987 polthet(j,i)=polthet_all(j,i,iparm)
990 gthet(j,i)=gthet_all(j,i,iparm)
992 theta0(i)=theta0_all(i,iparm)
993 sig0(i)=sig0_all(i,iparm)
994 sigc0(i)=sigc0_all(i,iparm)
997 IF (tor_mode.eq.0) THEN
998 nthetyp=nthetyp_all(iparm)
999 ntheterm=ntheterm_all(iparm)
1000 ntheterm2=ntheterm2_all(iparm)
1001 ntheterm3=ntheterm3_all(iparm)
1002 nsingle=nsingle_all(iparm)
1003 ndouble=ndouble_all(iparm)
1004 nntheterm=nntheterm_all(iparm)
1006 ithetyp(i)=ithetyp_all(i,iparm)
1009 do i=-nthetyp-1,nthetyp+1
1010 do j=-nthetyp-1,nthetyp+1
1011 do k=-nthetyp-1,nthetyp+1
1012 aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
1014 aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
1018 bbthet(m,l,i,j,k,iblock)= &
1019 bbthet_all(m,l,i,j,k,iblock,iparm)
1020 ccthet(m,l,i,j,k,iblock)= &
1021 ccthet_all(m,l,i,j,k,iblock,iparm)
1022 ddthet(m,l,i,j,k,iblock)= &
1023 ddthet_all(m,l,i,j,k,iblock,iparm)
1024 eethet(m,l,i,j,k,iblock)= &
1025 eethet_all(m,l,i,j,k,iblock,iparm)
1031 if (iblock.eq.1) then
1032 ffthet(mm,m,l,i,j,k,iblock)= &
1033 ffthet_all1(mm,m,l,i,j,k,iparm)
1034 ggthet(mm,m,l,i,j,k,iblock)= &
1035 ggthet_all1(mm,m,l,i,j,k,iparm)
1037 ffthet(mm,m,l,i,j,k,iblock)= &
1038 ffthet_all2(mm,m,l,i,j,k,iparm)
1039 ggthet(mm,m,l,i,j,k,iblock)= &
1040 ggthet_all2(mm,m,l,i,j,k,iparm)
1050 write (iout,*) "Need storing parameters"
1053 ! Restore the sidechain rotamer parameters
1058 nlob(iii)=nlob_all(iii,iparm)
1060 bsc(j,iii)=bsc_all(j,iii,iparm)
1062 censc(k,j,i)=censc_all(k,j,i,iparm)
1066 gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1074 sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1078 IF (tor_mode.eq.0) THEN
1079 ! Restore the torsional parameters
1081 do i=-ntortyp+1,ntortyp-1
1082 do j=-ntortyp+1,ntortyp-1
1083 v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1084 nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1085 nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1086 do k=1,nterm(i,j,iblock)
1087 v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1088 v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1090 do k=1,nlor(i,j,iblock)
1091 vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1092 vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1093 vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1098 ! Restore the double torsional parameters
1100 do i=-ntortyp+1,ntortyp-1
1101 do j=-ntortyp+1,ntortyp-1
1102 do k=-ntortyp+1,ntortyp-1
1103 ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1104 ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1105 do l=1,ntermd_1(i,j,k,iblock)
1106 v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1107 v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1108 v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1109 v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1111 do l=1,ntermd_2(i,j,k,iblock)
1112 do m=1,ntermd_2(i,j,k,iblock)
1113 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1121 ! Restore parameters of the cumulants
1124 b1(j,i)=b1_all(j,i,iparm)
1125 b1tilde(j,i)=b1tilde_all(j,i,iparm)
1126 b2(j,i)=b2_all(j,i,iparm)
1130 cc(k,j,i)=cc_all(k,j,i,iparm)
1131 ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1132 dd(k,j,i)=dd_all(k,j,i,iparm)
1133 dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1134 ee(k,j,i)=ee_all(k,j,i,iparm)
1139 write(iout,*) "need storing parameters"
1141 ! Restore the parameters of electrostatic interactions
1144 app(j,i)=app_all(j,i,iparm)
1145 bpp(j,i)=bpp_all(j,i,iparm)
1146 ael6(j,i)=ael6_all(j,i,iparm)
1147 ael3(j,i)=ael3_all(j,i,iparm)
1150 ! Restore sidechain parameters
1153 aa_aq(j,i)=aa_all(j,i,iparm)
1154 bb_aq(j,i)=bb_all(j,i,iparm)
1155 r0(j,i)=r0_all(j,i,iparm)
1156 sigma(j,i)=sigma_all(j,i,iparm)
1157 chi(j,i)=chi_all(j,i,iparm)
1158 augm(j,i)=augm_all(j,i,iparm)
1159 eps(j,i)=eps_all(j,i,iparm)
1163 chip(i)=chip_all(i,iparm)
1164 alp(i)=alp_all(i,iparm)
1166 ! Restore the SCp parameters
1169 aad(i,j)=aad_all(i,j,iparm)
1170 bad(i,j)=bad_all(i,j,iparm)
1173 ! Restore disulfide-bond parameters
1175 d0cm=d0cm_all(iparm)
1176 akcm=akcm_all(iparm)
1177 akth=akth_all(iparm)
1178 akct=akct_all(iparm)
1179 v1ss=v1ss_all(iparm)
1180 v2ss=v2ss_all(iparm)
1181 v3ss=v3ss_all(iparm)
1182 ! Restore SC-backbone correlation parameters
1183 do i=-nsccortyp,nsccortyp
1184 do j=-nsccortyp,nsccortyp
1186 nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1188 do k=1,nterm_sccor(j,i)
1189 v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1190 v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1196 end subroutine restore_parm
1197 !--------------------------------------------------------------------------
1199 !--------------------------------------------------------------------------
1200 subroutine make_ensembles(islice,*)
1201 ! construct the conformational ensembles at REMD temperatures
1202 use geometry_data, only:c
1203 use io_base, only:ilen
1204 use io_wham, only:pdboutW
1205 use geometry, only:returnbox
1207 ! include "DIMENSIONS"
1208 ! include "DIMENSIONS.ZSCOPT"
1209 ! include "DIMENSIONS.FREE"
1212 ! include "COMMON.MPI"
1213 integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1215 ! include "COMMON.IOUNITS"
1216 ! include "COMMON.CONTROL"
1217 ! include "COMMON.FREE"
1218 ! include "COMMON.ENERGIES"
1219 ! include "COMMON.FFIELD"
1220 ! include "COMMON.INTERACT"
1221 ! include "COMMON.SBRIDGE"
1222 ! include "COMMON.CHAIN"
1223 ! include "COMMON.PROTFILES"
1224 ! include "COMMON.PROT"
1225 real(kind=4) :: csingle(3,nres*2)
1226 real(kind=8),dimension(6) :: fT,fTprim,fTbis
1227 real(kind=8) :: quot,quotl1,quotl,kfacl,&
1228 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1229 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1230 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1231 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
1232 ecation_prot, ecationcation,evdwpp,eespp,evdwpsb,eelpsb,evdwsb, &
1233 eelsb,estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
1234 ecorr_nucl,ecorr3_nucl,escbase, epepbase,escpho, epeppho,ecation_nucl
1236 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1237 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1238 character(len=80) :: bxname
1239 character(len=2) :: licz1,licz2
1240 character(len=3) :: licz3,licz4
1241 character(len=5) :: ctemper
1244 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1245 real(kind=8) :: enepot(MaxStr)
1246 integer :: iperm(MaxStr)
1248 integer,dimension(0:nprocs) :: scount_
1251 if (me.eq.Master) then
1253 write (licz2,'(bz,i2.2)') islice
1254 if (nslice.eq.1) then
1255 if (.not.separate_parset) then
1256 bxname = prefix(:ilen(prefix))//".bx"
1258 write (licz3,'(bz,i3.3)') myparm
1259 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1262 if (.not.separate_parset) then
1263 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1265 write (licz3,'(bz,i3.3)') myparm
1266 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1267 "_slice_"//licz2//".bx"
1270 open (ientout,file=bxname,status="unknown",&
1271 form="unformatted",access="direct",recl=lenrec1)
1276 if (iparm.ne.iparmprint) exit
1277 call restore_parm(iparm)
1280 write (iout,*) "iparm",iparm," ib",ib
1282 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1283 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1290 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1292 !el old rescale weights
1294 ! if (rescale_mode.eq.1) then
1295 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1297 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1298 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1299 #elif defined(FUNCT)
1310 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1312 ! else if (rescale_mode.eq.2) then
1313 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1314 !#if defined(FUNCTH)
1315 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1316 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1317 !#elif defined(FUNCT)
1325 ! fT(l)=1.12692801104297249644d0/ &
1326 ! dlog(dexp(quotl)+dexp(-quotl))
1328 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1329 ! else if (rescale_mode.eq.0) then
1335 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1339 ! el end old rescale weihgts
1340 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1347 evdw=enetb(1,i,iparm)
1348 ! evdw_t=enetb(21,i,iparm)
1349 evdw_t=enetb(20,i,iparm)
1351 ! evdw2_14=enetb(17,i,iparm)
1352 evdw2_14=enetb(18,i,iparm)
1353 evdw2=enetb(2,i,iparm)+evdw2_14
1355 evdw2=enetb(2,i,iparm)
1359 ees=enetb(3,i,iparm)
1360 evdw1=enetb(16,i,iparm)
1362 ees=enetb(3,i,iparm)
1365 ecorr=enetb(4,i,iparm)
1366 ecorr5=enetb(5,i,iparm)
1367 ecorr6=enetb(6,i,iparm)
1368 eel_loc=enetb(7,i,iparm)
1369 eello_turn3=enetb(8,i,iparm)
1370 eello_turn4=enetb(9,i,iparm)
1371 eello_turn6=enetb(10,i,iparm)
1372 ebe=enetb(11,i,iparm)
1373 escloc=enetb(12,i,iparm)
1374 etors=enetb(13,i,iparm)
1375 etors_d=enetb(14,i,iparm)
1376 ehpb=enetb(15,i,iparm)
1378 estr=enetb(17,i,iparm)
1379 ! estr=enetb(18,i,iparm)
1380 ! esccor=enetb(19,i,iparm)
1381 esccor=enetb(21,i,iparm)
1382 ! edihcnstr=enetb(20,i,iparm)
1383 edihcnstr=enetb(19,i,iparm)
1384 ecationcation=enetb(42,i,iparm)
1385 ecation_prot=enetb(41,i,iparm)
1386 evdwpp = enetb(26,i,iparm)
1387 eespp = enetb(27,i,iparm)
1388 evdwpsb = enetb(28,i,iparm)
1389 eelpsb = enetb(29,i,iparm)
1390 evdwsb = enetb(30,i,iparm)
1391 eelsb = enetb(31,i,iparm)
1392 estr_nucl = enetb(32,i,iparm)
1393 ebe_nucl = enetb(33,i,iparm)
1394 esbloc = enetb(34,i,iparm)
1395 etors_nucl = enetb(35,i,iparm)
1396 etors_d_nucl = enetb(36,i,iparm)
1397 ecorr_nucl = enetb(37,i,iparm)
1398 ecorr3_nucl = enetb(38,i,iparm)
1399 epeppho= enetb(49,i,iparm)
1400 escpho= enetb(48,i,iparm)
1401 epepbase= enetb(47,i,iparm)
1402 escbase= enetb(46,i,iparm)
1403 ecation_nucl=enetb(50,i,iparm)
1410 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1412 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1413 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1414 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1415 ! +ft(2)*wturn3*eello_turn3 &
1416 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1417 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1420 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1421 ! +ft(1)*welec*(ees+evdw1) &
1422 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1423 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1424 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1425 ! +ft(2)*wturn3*eello_turn3 &
1426 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1427 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1432 etot=wsc*evdw+wscp*evdw2+welec*ees &
1434 +wang*ebe+wtor*etors+wscloc*escloc &
1435 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1436 +wcorr6*ecorr6+wturn4*eello_turn4 &
1437 +wturn3*eello_turn3 &
1438 +wturn6*eello_turn6+wel_loc*eel_loc &
1439 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1440 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1441 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1442 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1443 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1444 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1446 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1447 +wcatnucl*ecation_nucl
1450 etot=wsc*evdw+wscp*evdw2 &
1451 +welec*(ees+evdw1) &
1452 +wang*ebe+wtor*etors+wscloc*escloc &
1453 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1454 +wcorr6*ecorr6+wturn4*eello_turn4 &
1455 +wturn3*eello_turn3 &
1456 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1457 +wtor_d*etors_d+wsccor*esccor &
1458 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1459 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1460 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1461 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1462 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1464 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1465 +wcatnucl*ecation_nucl
1471 beta_h(ib,iparm)*etot-entfac(i)
1474 write (iout,*) i,indstart(me)+i-1,ib,&
1475 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1476 -entfac(i),Fdimless(i)
1479 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1486 Fdimless_(i)=Fdimless(i)
1488 write(iout,*) "before gather Fdimless"
1489 write(iout,*) scount(me),scount_(0),idispl(0)
1490 write (iout,*) "added update of scount_"
1491 call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
1492 MPI_INTEGER, WHAM_COMM, IERROR)
1495 call MPI_Gatherv(Fdimless_(1),scount(me),&
1496 MPI_REAL,Fdimless(1),&
1497 scount_(0),idispl(0),MPI_REAL,Master,&
1499 write(iout,*) "after gather Fdimless"
1501 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1502 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1503 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1505 call MPI_Gatherv(entfac(1),scount(me),&
1506 MPI_DOUBLE_PRECISION,entfac(1),&
1507 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1510 if (me.eq.Master) then
1512 write (iout,*) "The FDIMLESS array before sorting"
1514 write (iout,*) i,fdimless(i)
1521 call mysort1(ntot(islice),Fdimless,iperm)
1523 write (iout,*) "The FDIMLESS array after sorting"
1525 write (iout,*) i,iperm(i),fdimless(i)
1530 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1532 ! write (iout,*) "qfree",qfree
1535 do i=1,min0(ntot(islice),ensembles)
1536 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1538 write (iout,*) i,ib,beta_h(ib,iparm),&
1539 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1540 potE(iperm(i),iparm),&
1541 -entfac(iperm(i)),fdimless(i),sumprob
1543 if (sumprob.gt.0.99d0) goto 122
1549 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1551 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1556 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1559 if (iproc.ge.nprocs) then
1560 write (iout,*) "Fatal error: processor out of range",iproc
1565 close (ientout,status="delete")
1569 ik=ii-indstart(iproc)+1
1570 if (iproc.ne.Master) then
1571 if (me.eq.iproc) then
1573 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1574 " energy",potE(ik,iparm)
1576 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1577 Master,i,WHAM_COMM,IERROR)
1578 else if (me.eq.Master) then
1579 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1580 WHAM_COMM,STATUS,IERROR)
1582 else if (me.eq.Master) then
1583 enepot(i)=potE(ik,iparm)
1588 enepot(i)=potE(iperm(i),iparm)
1592 if (me.eq.Master) then
1594 write(licz3,'(bz,i3.3)') iparm
1595 write(licz2,'(bz,i2.2)') islice
1596 if (temper.lt.100.0d0) then
1597 write(ctemper,'(f3.0)') temper
1598 else if (temper.lt.1000.0) then
1599 write (ctemper,'(f4.0)') temper
1601 write (ctemper,'(f5.0)') temper
1603 if (nparmset.eq.1) then
1604 if (separate_parset) then
1605 write(licz4,'(bz,i3.3)') myparm
1606 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1608 pdbname=prefix(:ilen(prefix))
1611 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1613 if (nslice.eq.1) then
1614 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1615 ctemper(:ilen(ctemper))//"pdb"
1617 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1618 ctemper(:ilen(ctemper))//"pdb"
1620 open(ipdb,file=pdbname)
1622 read (ientout,rec=iperm(i)) &
1623 ((csingle(l,k),l=1,3),k=1,nres),&
1624 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1625 nss,(ihpb(k),jhpb(k),k=1,nss),&
1626 eini,efree,rmsdev,iscor
1640 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1650 close(ientout,status="delete")
1653 scount(i)=scount_(i)
1656 end subroutine make_ensembles
1657 !--------------------------------------------------------------------------
1658 subroutine mysort1(n, x, ipermut)
1660 integer :: i,j,imax,ipm,n
1661 real(kind=4) :: x(n)
1662 integer :: ipermut(n)
1663 real(kind=4) :: xtemp
1668 if (x(j).lt.xtemp) then
1676 ipermut(imax)=ipermut(i)
1680 end subroutine mysort1
1681 !--------------------------------------------------------------------------
1682 subroutine alloc_enecalc_arrays(iparm)
1685 use geometry_data, only:maxlob
1687 !---------------------------
1688 ! COMMON.ENERGIES form wham_data
1690 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1691 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1692 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1693 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1695 ! allocate ENECALC arrays
1696 !---------------------------
1699 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1700 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1701 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1702 aksc_all(maxbondterm,ntyp,iparm),&
1703 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1704 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1705 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1706 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1707 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1708 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1709 allocate(theta0_all(-ntyp:ntyp,iparm),&
1710 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1711 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1712 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1713 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1714 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1715 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1716 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1717 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1718 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1719 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1720 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1721 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1722 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1723 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1724 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1725 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1726 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1727 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1728 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1729 -nthetyp-1:nthetyp+1,iparm))
1730 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1731 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1732 -nthetyp-1:nthetyp+1,iparm))
1733 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1734 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1735 -nthetyp-1:nthetyp+1,iparm))
1736 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1737 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1738 -nthetyp-1:nthetyp+1,iparm))
1739 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1740 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1741 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1742 allocate(bsc_all(maxlob,ntyp,iparm))
1743 !(maxlob,ntyp,max_parm)
1744 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1745 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1746 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1747 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1748 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1749 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1750 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1751 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1752 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1753 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1754 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1755 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1756 -ntortyp:ntortyp,2,iparm))
1757 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1758 -ntortyp:ntortyp,2,iparm))
1759 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1760 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1761 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1762 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1763 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1764 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1765 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1766 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1767 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1768 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1769 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1770 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1771 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1772 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1773 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1774 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1775 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1776 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1777 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1778 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1779 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1780 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1781 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1782 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1783 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1784 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1785 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1786 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1787 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1788 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1789 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1790 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1791 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1792 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1793 -ntortyp:ntortyp,2,iparm))
1794 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1795 -ntortyp:ntortyp,2,iparm))
1796 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1797 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1798 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1799 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1800 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1801 nsingle_all(iparm),&
1802 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1803 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1805 end subroutine alloc_enecalc_arrays
1806 !--------------------------------------------------------------------------
1807 !--------------------------------------------------------------------------