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,&
1235 elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang, &
1239 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1240 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1241 character(len=80) :: bxname
1242 character(len=2) :: licz1,licz2
1243 character(len=3) :: licz3,licz4
1244 character(len=5) :: ctemper
1247 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1248 real(kind=8) :: enepot(MaxStr)
1249 integer :: iperm(MaxStr)
1251 integer,dimension(0:nprocs) :: scount_
1252 write(iout,*) "Begin make ensemble"
1254 if (me.eq.Master) then
1256 write (licz2,'(bz,i2.2)') islice
1257 if (nslice.eq.1) then
1258 if (.not.separate_parset) then
1259 bxname = prefix(:ilen(prefix))//".bx"
1261 write (licz3,'(bz,i3.3)') myparm
1262 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1265 if (.not.separate_parset) then
1266 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1268 write (licz3,'(bz,i3.3)') myparm
1269 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1270 "_slice_"//licz2//".bx"
1273 open (ientout,file=bxname,status="unknown",&
1274 form="unformatted",access="direct",recl=lenrec1)
1278 write(iout,*) "iparmprint",iparmprint,iparm
1280 if (iparm.ne.iparmprint) exit
1281 call restore_parm(iparm)
1284 write (iout,*) "iparm",iparm," ib",ib
1286 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1287 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1294 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1296 !el old rescale weights
1298 ! if (rescale_mode.eq.1) then
1299 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1301 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1302 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1303 #elif defined(FUNCT)
1314 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1316 ! else if (rescale_mode.eq.2) then
1317 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1318 !#if defined(FUNCTH)
1319 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1320 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1321 !#elif defined(FUNCT)
1329 ! fT(l)=1.12692801104297249644d0/ &
1330 ! dlog(dexp(quotl)+dexp(-quotl))
1332 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1333 ! else if (rescale_mode.eq.0) then
1339 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1343 ! el end old rescale weihgts
1344 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1351 evdw=enetb(1,i,iparm)
1352 ! evdw_t=enetb(21,i,iparm)
1353 evdw_t=enetb(20,i,iparm)
1355 ! evdw2_14=enetb(17,i,iparm)
1356 evdw2_14=enetb(18,i,iparm)
1357 evdw2=enetb(2,i,iparm)+evdw2_14
1359 evdw2=enetb(2,i,iparm)
1363 ees=enetb(3,i,iparm)
1364 evdw1=enetb(16,i,iparm)
1366 ees=enetb(3,i,iparm)
1369 ecorr=enetb(4,i,iparm)
1370 ecorr5=enetb(5,i,iparm)
1371 ecorr6=enetb(6,i,iparm)
1372 eel_loc=enetb(7,i,iparm)
1373 eello_turn3=enetb(8,i,iparm)
1374 eello_turn4=enetb(9,i,iparm)
1375 eello_turn6=enetb(10,i,iparm)
1376 ebe=enetb(11,i,iparm)
1377 escloc=enetb(12,i,iparm)
1378 etors=enetb(13,i,iparm)
1379 etors_d=enetb(14,i,iparm)
1380 ehpb=enetb(15,i,iparm)
1382 estr=enetb(17,i,iparm)
1383 ! estr=enetb(18,i,iparm)
1384 ! esccor=enetb(19,i,iparm)
1385 esccor=enetb(21,i,iparm)
1386 ! edihcnstr=enetb(20,i,iparm)
1387 edihcnstr=enetb(19,i,iparm)
1388 ecationcation=enetb(42,i,iparm)
1389 ecation_prot=enetb(41,i,iparm)
1390 evdwpp = enetb(26,i,iparm)
1391 eespp = enetb(27,i,iparm)
1392 evdwpsb = enetb(28,i,iparm)
1393 eelpsb = enetb(29,i,iparm)
1394 evdwsb = enetb(30,i,iparm)
1395 eelsb = enetb(31,i,iparm)
1396 estr_nucl = enetb(32,i,iparm)
1397 ebe_nucl = enetb(33,i,iparm)
1398 esbloc = enetb(34,i,iparm)
1399 etors_nucl = enetb(35,i,iparm)
1400 etors_d_nucl = enetb(36,i,iparm)
1401 ecorr_nucl = enetb(37,i,iparm)
1402 ecorr3_nucl = enetb(38,i,iparm)
1403 epeppho= enetb(49,i,iparm)
1404 escpho= enetb(48,i,iparm)
1405 epepbase= enetb(47,i,iparm)
1406 escbase= enetb(46,i,iparm)
1407 ecation_nucl=enetb(50,i,iparm)
1408 elipbond=enetb(52,i,iparm)
1409 elipang=enetb(53,i,iparm)
1410 eliplj=enetb(54,i,iparm)
1411 elipelec=enetb(55,i,iparm)
1412 ecat_prottran=enetb(56,i,iparm)
1413 ecation_protang=enetb(57,i,iparm)
1414 eliptran=enetb(22,i,iparm)
1421 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1423 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1424 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1425 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1426 ! +ft(2)*wturn3*eello_turn3 &
1427 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1428 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1431 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1432 ! +ft(1)*welec*(ees+evdw1) &
1433 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1434 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1435 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1436 ! +ft(2)*wturn3*eello_turn3 &
1437 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1438 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1443 etot=wsc*evdw+wscp*evdw2+welec*ees &
1445 +wang*ebe+wtor*etors+wscloc*escloc &
1446 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1447 +wcorr6*ecorr6+wturn4*eello_turn4 &
1448 +wturn3*eello_turn3 &
1449 +wturn6*eello_turn6+wel_loc*eel_loc &
1450 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1451 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1452 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1453 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1454 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1455 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1457 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1458 +wcatnucl*ecation_nucl&
1459 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang+&
1464 etot=wsc*evdw+wscp*evdw2 &
1465 +welec*(ees+evdw1) &
1466 +wang*ebe+wtor*etors+wscloc*escloc &
1467 +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
1468 +wcorr6*ecorr6+wturn4*eello_turn4 &
1469 +wturn3*eello_turn3 &
1470 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1471 +wtor_d*etors_d+wsccor*esccor &
1472 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1473 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1474 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1475 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1476 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1478 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1479 +wcatnucl*ecation_nucl&
1480 +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1489 beta_h(ib,iparm)*etot-entfac(i)
1492 write (iout,*) i,indstart(me)+i-1,ib,&
1493 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1494 -entfac(i),Fdimless(i)
1497 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1504 Fdimless_(i)=Fdimless(i)
1506 write(iout,*) "before gather Fdimless"
1507 write(iout,*) scount(me),scount_(0),idispl(0)
1508 write (iout,*) "added update of scount_"
1509 call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
1510 MPI_INTEGER, WHAM_COMM, IERROR)
1513 call MPI_Gatherv(Fdimless_(1),scount(me),&
1514 MPI_REAL,Fdimless(1),&
1515 scount_(0),idispl(0),MPI_REAL,Master,&
1517 write(iout,*) "after gather Fdimless"
1519 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1520 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1521 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1523 call MPI_Gatherv(entfac(1),scount(me),&
1524 MPI_DOUBLE_PRECISION,entfac(1),&
1525 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1528 if (me.eq.Master) then
1530 write (iout,*) "The FDIMLESS array before sorting"
1532 write (iout,*) i,fdimless(i)
1539 call mysort1(ntot(islice),Fdimless,iperm)
1541 write (iout,*) "The FDIMLESS array after sorting"
1543 write (iout,*) i,iperm(i),fdimless(i)
1548 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1550 ! write (iout,*) "qfree",qfree
1553 do i=1,min0(ntot(islice),ensembles)
1554 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1556 write (iout,*) i,ib,beta_h(ib,iparm),&
1557 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1558 potE(iperm(i),iparm),&
1559 -entfac(iperm(i)),fdimless(i),sumprob
1561 if (sumprob.gt.0.99d0) goto 122
1567 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1569 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1574 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1577 if (iproc.ge.nprocs) then
1578 write (iout,*) "Fatal error: processor out of range",iproc
1583 close (ientout,status="delete")
1587 ik=ii-indstart(iproc)+1
1588 if (iproc.ne.Master) then
1589 if (me.eq.iproc) then
1591 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1592 " energy",potE(ik,iparm)
1594 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1595 Master,i,WHAM_COMM,IERROR)
1596 else if (me.eq.Master) then
1597 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1598 WHAM_COMM,STATUS,IERROR)
1600 else if (me.eq.Master) then
1601 enepot(i)=potE(ik,iparm)
1606 enepot(i)=potE(iperm(i),iparm)
1609 write(iout,*) "DEBUG",me
1611 if (me.eq.Master) then
1613 write(licz3,'(bz,i3.3)') iparm
1614 write(licz2,'(bz,i2.2)') islice
1615 if (temper.lt.100.0d0) then
1616 write(ctemper,'(f3.0)') temper
1617 else if (temper.lt.1000.0) then
1618 write (ctemper,'(f4.0)') temper
1620 write (ctemper,'(f5.0)') temper
1622 if (nparmset.eq.1) then
1623 if (separate_parset) then
1624 write(licz4,'(bz,i3.3)') myparm
1625 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1627 pdbname=prefix(:ilen(prefix))
1630 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1632 if (nslice.eq.1) then
1633 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1634 ctemper(:ilen(ctemper))//"pdb"
1636 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1637 ctemper(:ilen(ctemper))//"pdb"
1639 open(ipdb,file=pdbname)
1641 read (ientout,rec=iperm(i)) &
1642 ((csingle(l,k),l=1,3),k=1,nres),&
1643 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1644 nss,(ihpb(k),jhpb(k),k=1,nss),&
1645 eini,efree,rmsdev,iscor
1659 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1669 close(ientout,status="delete")
1672 scount(i)=scount_(i)
1675 end subroutine make_ensembles
1676 !--------------------------------------------------------------------------
1677 subroutine mysort1(n, x, ipermut)
1679 integer :: i,j,imax,ipm,n
1680 real(kind=4) :: x(n)
1681 integer :: ipermut(n)
1682 real(kind=4) :: xtemp
1687 if (x(j).lt.xtemp) then
1695 ipermut(imax)=ipermut(i)
1699 end subroutine mysort1
1700 !--------------------------------------------------------------------------
1701 subroutine alloc_enecalc_arrays(iparm)
1704 use geometry_data, only:maxlob
1706 !---------------------------
1707 ! COMMON.ENERGIES form wham_data
1709 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1710 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1711 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1712 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1714 ! allocate ENECALC arrays
1715 !---------------------------
1718 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1719 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1720 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1721 aksc_all(maxbondterm,ntyp,iparm),&
1722 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1723 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1724 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1725 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1726 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1727 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1728 allocate(theta0_all(-ntyp:ntyp,iparm),&
1729 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1730 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1731 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1732 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1733 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1734 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1735 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1736 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1737 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1738 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1739 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1740 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1741 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1742 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1743 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1744 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1745 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1746 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1747 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1748 -nthetyp-1:nthetyp+1,iparm))
1749 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1750 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1751 -nthetyp-1:nthetyp+1,iparm))
1752 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1753 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1754 -nthetyp-1:nthetyp+1,iparm))
1755 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1756 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1757 -nthetyp-1:nthetyp+1,iparm))
1758 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1759 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1760 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1761 allocate(bsc_all(maxlob,ntyp,iparm))
1762 !(maxlob,ntyp,max_parm)
1763 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1764 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1765 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1766 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1767 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1768 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1769 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1770 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1771 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1772 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1773 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1774 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1775 -ntortyp:ntortyp,2,iparm))
1776 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1777 -ntortyp:ntortyp,2,iparm))
1778 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1779 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1780 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1781 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1782 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1783 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1784 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1785 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1786 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1787 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1788 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1789 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1790 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1791 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1792 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1793 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1794 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1795 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1796 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1797 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1798 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1799 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1800 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1801 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1802 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1803 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1804 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1805 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1806 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1807 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1808 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1809 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1810 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1811 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1812 -ntortyp:ntortyp,2,iparm))
1813 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1814 -ntortyp:ntortyp,2,iparm))
1815 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1816 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1817 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1818 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1819 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1820 nsingle_all(iparm),&
1821 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1822 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1824 end subroutine alloc_enecalc_arrays
1825 !--------------------------------------------------------------------------
1826 !--------------------------------------------------------------------------