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
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)
1409 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1411 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1412 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1413 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1414 ! +ft(2)*wturn3*eello_turn3 &
1415 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1416 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1419 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1420 ! +ft(1)*welec*(ees+evdw1) &
1421 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1422 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1423 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1424 ! +ft(2)*wturn3*eello_turn3 &
1425 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1426 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1431 etot=wsc*evdw+wscp*evdw2+welec*ees &
1433 +wang*ebe+wtor*etors+wscloc*escloc &
1434 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1435 +wcorr6*ecorr6+wturn4*eello_turn4 &
1436 +wturn3*eello_turn3 &
1437 +wturn6*eello_turn6+wel_loc*eel_loc &
1438 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1439 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1440 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1441 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1442 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1443 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1445 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1448 etot=wsc*evdw+wscp*evdw2 &
1449 +welec*(ees+evdw1) &
1450 +wang*ebe+wtor*etors+wscloc*escloc &
1451 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1452 +wcorr6*ecorr6+wturn4*eello_turn4 &
1453 +wturn3*eello_turn3 &
1454 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1455 +wtor_d*etors_d+wsccor*esccor &
1456 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1457 +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1458 +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1459 +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1460 +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
1462 +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
1468 beta_h(ib,iparm)*etot-entfac(i)
1471 write (iout,*) i,indstart(me)+i-1,ib,&
1472 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1473 -entfac(i),Fdimless(i)
1476 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1483 Fdimless_(i)=Fdimless(i)
1485 write(iout,*) "before gather Fdimless"
1486 write(iout,*) scount(me),scount_(0),idispl(0)
1487 call MPI_Gatherv(Fdimless_(1),scount(me),&
1488 MPI_REAL,Fdimless(1),&
1489 scount_(0),idispl(0),MPI_REAL,Master,&
1491 write(iout,*) "after gather Fdimless"
1493 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1494 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1495 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1497 call MPI_Gatherv(entfac(1),scount(me),&
1498 MPI_DOUBLE_PRECISION,entfac(1),&
1499 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1502 if (me.eq.Master) then
1504 write (iout,*) "The FDIMLESS array before sorting"
1506 write (iout,*) i,fdimless(i)
1513 call mysort1(ntot(islice),Fdimless,iperm)
1515 write (iout,*) "The FDIMLESS array after sorting"
1517 write (iout,*) i,iperm(i),fdimless(i)
1522 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1524 ! write (iout,*) "qfree",qfree
1527 do i=1,min0(ntot(islice),ensembles)
1528 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1530 write (iout,*) i,ib,beta_h(ib,iparm),&
1531 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1532 potE(iperm(i),iparm),&
1533 -entfac(iperm(i)),fdimless(i),sumprob
1535 if (sumprob.gt.0.99d0) goto 122
1541 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1543 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1548 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1551 if (iproc.ge.nprocs) then
1552 write (iout,*) "Fatal error: processor out of range",iproc
1557 close (ientout,status="delete")
1561 ik=ii-indstart(iproc)+1
1562 if (iproc.ne.Master) then
1563 if (me.eq.iproc) then
1565 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1566 " energy",potE(ik,iparm)
1568 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1569 Master,i,WHAM_COMM,IERROR)
1570 else if (me.eq.Master) then
1571 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1572 WHAM_COMM,STATUS,IERROR)
1574 else if (me.eq.Master) then
1575 enepot(i)=potE(ik,iparm)
1580 enepot(i)=potE(iperm(i),iparm)
1584 if (me.eq.Master) then
1586 write(licz3,'(bz,i3.3)') iparm
1587 write(licz2,'(bz,i2.2)') islice
1588 if (temper.lt.100.0d0) then
1589 write(ctemper,'(f3.0)') temper
1590 else if (temper.lt.1000.0) then
1591 write (ctemper,'(f4.0)') temper
1593 write (ctemper,'(f5.0)') temper
1595 if (nparmset.eq.1) then
1596 if (separate_parset) then
1597 write(licz4,'(bz,i3.3)') myparm
1598 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1600 pdbname=prefix(:ilen(prefix))
1603 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1605 if (nslice.eq.1) then
1606 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1607 ctemper(:ilen(ctemper))//"pdb"
1609 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1610 ctemper(:ilen(ctemper))//"pdb"
1612 open(ipdb,file=pdbname)
1614 read (ientout,rec=iperm(i)) &
1615 ((csingle(l,k),l=1,3),k=1,nres),&
1616 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1617 nss,(ihpb(k),jhpb(k),k=1,nss),&
1618 eini,efree,rmsdev,iscor
1632 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1642 close(ientout,status="delete")
1645 scount(i)=scount_(i)
1648 end subroutine make_ensembles
1649 !--------------------------------------------------------------------------
1650 subroutine mysort1(n, x, ipermut)
1652 integer :: i,j,imax,ipm,n
1653 real(kind=4) :: x(n)
1654 integer :: ipermut(n)
1655 real(kind=4) :: xtemp
1660 if (x(j).lt.xtemp) then
1668 ipermut(imax)=ipermut(i)
1672 end subroutine mysort1
1673 !--------------------------------------------------------------------------
1674 subroutine alloc_enecalc_arrays(iparm)
1677 use geometry_data, only:maxlob
1679 !---------------------------
1680 ! COMMON.ENERGIES form wham_data
1682 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1683 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1684 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1685 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1687 ! allocate ENECALC arrays
1688 !---------------------------
1691 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1692 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1693 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1694 aksc_all(maxbondterm,ntyp,iparm),&
1695 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1696 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1697 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1698 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1699 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1700 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1701 allocate(theta0_all(-ntyp:ntyp,iparm),&
1702 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1703 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1704 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1705 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1706 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1707 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1708 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1709 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1710 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1711 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1712 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1713 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1714 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1715 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1716 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1717 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1718 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1719 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1720 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1721 -nthetyp-1:nthetyp+1,iparm))
1722 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1723 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1724 -nthetyp-1:nthetyp+1,iparm))
1725 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1726 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1727 -nthetyp-1:nthetyp+1,iparm))
1728 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1729 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1730 -nthetyp-1:nthetyp+1,iparm))
1731 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1732 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1733 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1734 allocate(bsc_all(maxlob,ntyp,iparm))
1735 !(maxlob,ntyp,max_parm)
1736 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1737 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1738 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1739 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1740 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1741 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1742 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1743 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1744 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1745 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1746 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1747 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1748 -ntortyp:ntortyp,2,iparm))
1749 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1750 -ntortyp:ntortyp,2,iparm))
1751 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1752 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1753 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1754 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1755 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1756 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1757 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1758 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1759 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1760 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1761 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1762 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1763 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1764 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1765 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1766 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1767 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1768 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1769 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1770 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1771 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1772 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1773 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1774 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1775 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1776 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1777 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1778 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1779 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1780 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1781 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1782 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1783 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1784 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1785 -ntortyp:ntortyp,2,iparm))
1786 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1787 -ntortyp:ntortyp,2,iparm))
1788 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1789 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1790 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1791 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1792 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1793 nsingle_all(iparm),&
1794 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1795 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1797 end subroutine alloc_enecalc_arrays
1798 !--------------------------------------------------------------------------
1799 !--------------------------------------------------------------------------