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
1233 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1234 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1235 character(len=80) :: bxname
1236 character(len=2) :: licz1,licz2
1237 character(len=3) :: licz3,licz4
1238 character(len=5) :: ctemper
1241 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1242 real(kind=8) :: enepot(MaxStr)
1243 integer :: iperm(MaxStr)
1245 integer,dimension(0:nprocs) :: scount_
1248 if (me.eq.Master) then
1250 write (licz2,'(bz,i2.2)') islice
1251 if (nslice.eq.1) then
1252 if (.not.separate_parset) then
1253 bxname = prefix(:ilen(prefix))//".bx"
1255 write (licz3,'(bz,i3.3)') myparm
1256 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1259 if (.not.separate_parset) then
1260 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1262 write (licz3,'(bz,i3.3)') myparm
1263 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1264 "_slice_"//licz2//".bx"
1267 open (ientout,file=bxname,status="unknown",&
1268 form="unformatted",access="direct",recl=lenrec1)
1273 if (iparm.ne.iparmprint) exit
1274 call restore_parm(iparm)
1277 write (iout,*) "iparm",iparm," ib",ib
1279 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1280 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1287 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1289 !el old rescale weights
1291 ! if (rescale_mode.eq.1) then
1292 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1294 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1295 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1296 #elif defined(FUNCT)
1307 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1309 ! else if (rescale_mode.eq.2) then
1310 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1311 !#if defined(FUNCTH)
1312 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1313 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1314 !#elif defined(FUNCT)
1322 ! fT(l)=1.12692801104297249644d0/ &
1323 ! dlog(dexp(quotl)+dexp(-quotl))
1325 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1326 ! else if (rescale_mode.eq.0) then
1332 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1336 ! el end old rescale weihgts
1337 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1344 evdw=enetb(1,i,iparm)
1345 ! evdw_t=enetb(21,i,iparm)
1346 evdw_t=enetb(20,i,iparm)
1348 ! evdw2_14=enetb(17,i,iparm)
1349 evdw2_14=enetb(18,i,iparm)
1350 evdw2=enetb(2,i,iparm)+evdw2_14
1352 evdw2=enetb(2,i,iparm)
1356 ees=enetb(3,i,iparm)
1357 evdw1=enetb(16,i,iparm)
1359 ees=enetb(3,i,iparm)
1362 ecorr=enetb(4,i,iparm)
1363 ecorr5=enetb(5,i,iparm)
1364 ecorr6=enetb(6,i,iparm)
1365 eel_loc=enetb(7,i,iparm)
1366 eello_turn3=enetb(8,i,iparm)
1367 eello_turn4=enetb(9,i,iparm)
1368 eello_turn6=enetb(10,i,iparm)
1369 ebe=enetb(11,i,iparm)
1370 escloc=enetb(12,i,iparm)
1371 etors=enetb(13,i,iparm)
1372 etors_d=enetb(14,i,iparm)
1373 ehpb=enetb(15,i,iparm)
1375 estr=enetb(17,i,iparm)
1376 ! estr=enetb(18,i,iparm)
1377 ! esccor=enetb(19,i,iparm)
1378 esccor=enetb(21,i,iparm)
1379 ! edihcnstr=enetb(20,i,iparm)
1380 edihcnstr=enetb(19,i,iparm)
1381 ecationcation=enetb(42,i,iparm)
1382 ecation_prot=enetb(41,i,iparm)
1385 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1387 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1388 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1389 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1390 ! +ft(2)*wturn3*eello_turn3 &
1391 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1392 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1395 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1396 ! +ft(1)*welec*(ees+evdw1) &
1397 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1398 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1399 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1400 ! +ft(2)*wturn3*eello_turn3 &
1401 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1402 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1407 etot=wsc*evdw+wscp*evdw2+welec*ees &
1409 +wang*ebe+wtor*etors+wscloc*escloc &
1410 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1411 +wcorr6*ecorr6+wturn4*eello_turn4 &
1412 +wturn3*eello_turn3 &
1413 +wturn6*eello_turn6+wel_loc*eel_loc &
1414 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1415 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1417 etot=wsc*evdw+wscp*evdw2 &
1418 +welec*(ees+evdw1) &
1419 +wang*ebe+wtor*etors+wscloc*escloc &
1420 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1421 +wcorr6*ecorr6+wturn4*eello_turn4 &
1422 +wturn3*eello_turn3 &
1423 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1424 +wtor_d*etors_d+wsccor*esccor &
1425 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1430 beta_h(ib,iparm)*etot-entfac(i)
1433 write (iout,*) i,indstart(me)+i-1,ib,&
1434 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1435 -entfac(i),Fdimless(i)
1438 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1444 Fdimless_(i)=Fdimless(i)
1446 call MPI_Gatherv(Fdimless_(1),scount(me),&
1447 MPI_REAL,Fdimless(1),&
1448 scount_(0),idispl(0),MPI_REAL,Master,&
1451 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1452 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1453 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1455 call MPI_Gatherv(entfac(1),scount(me),&
1456 MPI_DOUBLE_PRECISION,entfac(1),&
1457 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1460 if (me.eq.Master) then
1462 write (iout,*) "The FDIMLESS array before sorting"
1464 write (iout,*) i,fdimless(i)
1471 call mysort1(ntot(islice),Fdimless,iperm)
1473 write (iout,*) "The FDIMLESS array after sorting"
1475 write (iout,*) i,iperm(i),fdimless(i)
1480 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1482 ! write (iout,*) "qfree",qfree
1485 do i=1,min0(ntot(islice),ensembles)
1486 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1488 write (iout,*) i,ib,beta_h(ib,iparm),&
1489 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1490 potE(iperm(i),iparm),&
1491 -entfac(iperm(i)),fdimless(i),sumprob
1493 if (sumprob.gt.0.99d0) goto 122
1499 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1501 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1506 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1509 if (iproc.ge.nprocs) then
1510 write (iout,*) "Fatal error: processor out of range",iproc
1515 close (ientout,status="delete")
1519 ik=ii-indstart(iproc)+1
1520 if (iproc.ne.Master) then
1521 if (me.eq.iproc) then
1523 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1524 " energy",potE(ik,iparm)
1526 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1527 Master,i,WHAM_COMM,IERROR)
1528 else if (me.eq.Master) then
1529 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1530 WHAM_COMM,STATUS,IERROR)
1532 else if (me.eq.Master) then
1533 enepot(i)=potE(ik,iparm)
1538 enepot(i)=potE(iperm(i),iparm)
1542 if (me.eq.Master) then
1544 write(licz3,'(bz,i3.3)') iparm
1545 write(licz2,'(bz,i2.2)') islice
1546 if (temper.lt.100.0d0) then
1547 write(ctemper,'(f3.0)') temper
1548 else if (temper.lt.1000.0) then
1549 write (ctemper,'(f4.0)') temper
1551 write (ctemper,'(f5.0)') temper
1553 if (nparmset.eq.1) then
1554 if (separate_parset) then
1555 write(licz4,'(bz,i3.3)') myparm
1556 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1558 pdbname=prefix(:ilen(prefix))
1561 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1563 if (nslice.eq.1) then
1564 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1565 ctemper(:ilen(ctemper))//"pdb"
1567 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1568 ctemper(:ilen(ctemper))//"pdb"
1570 open(ipdb,file=pdbname)
1572 read (ientout,rec=iperm(i)) &
1573 ((csingle(l,k),l=1,3),k=1,nres),&
1574 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1575 nss,(ihpb(k),jhpb(k),k=1,nss),&
1576 eini,efree,rmsdev,iscor
1590 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1600 close(ientout,status="delete")
1603 scount(i)=scount_(i)
1606 end subroutine make_ensembles
1607 !--------------------------------------------------------------------------
1608 subroutine mysort1(n, x, ipermut)
1610 integer :: i,j,imax,ipm,n
1611 real(kind=4) :: x(n)
1612 integer :: ipermut(n)
1613 real(kind=4) :: xtemp
1618 if (x(j).lt.xtemp) then
1626 ipermut(imax)=ipermut(i)
1630 end subroutine mysort1
1631 !--------------------------------------------------------------------------
1632 subroutine alloc_enecalc_arrays(iparm)
1635 use geometry_data, only:maxlob
1637 !---------------------------
1638 ! COMMON.ENERGIES form wham_data
1640 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1641 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1642 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1643 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1645 ! allocate ENECALC arrays
1646 !---------------------------
1649 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1650 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1651 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1652 aksc_all(maxbondterm,ntyp,iparm),&
1653 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1654 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1655 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1656 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1657 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1658 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1659 allocate(theta0_all(-ntyp:ntyp,iparm),&
1660 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1661 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1662 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1663 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1664 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1665 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1666 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1667 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1668 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1669 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1670 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1671 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1672 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1673 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1674 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1675 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1676 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1677 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1678 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1679 -nthetyp-1:nthetyp+1,iparm))
1680 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1681 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1682 -nthetyp-1:nthetyp+1,iparm))
1683 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1684 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1685 -nthetyp-1:nthetyp+1,iparm))
1686 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1687 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1688 -nthetyp-1:nthetyp+1,iparm))
1689 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1690 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1691 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1692 allocate(bsc_all(maxlob,ntyp,iparm))
1693 !(maxlob,ntyp,max_parm)
1694 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1695 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1696 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1697 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1698 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1699 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1700 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1701 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1702 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1703 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1704 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1705 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1706 -ntortyp:ntortyp,2,iparm))
1707 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1708 -ntortyp:ntortyp,2,iparm))
1709 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1710 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1711 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1712 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1713 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1714 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1715 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1716 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1717 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1718 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1719 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1720 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1721 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1722 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1723 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1724 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1725 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1726 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1727 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1728 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1729 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1730 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1731 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1732 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1733 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1734 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1735 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1736 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1737 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1738 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1739 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1740 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1741 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1742 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1743 -ntortyp:ntortyp,2,iparm))
1744 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1745 -ntortyp:ntortyp,2,iparm))
1746 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1747 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1748 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1749 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1750 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1751 nsingle_all(iparm),&
1752 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1753 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1755 end subroutine alloc_enecalc_arrays
1756 !--------------------------------------------------------------------------
1757 !--------------------------------------------------------------------------