2 !-----------------------------------------------------------------------------
6 use geometry_data, only:nres,boxxsize,boxysize,boxzsize
8 ! use control_data, only:maxthetyp1
9 use energy, only:etotal,enerprint,rescale_weights
13 ! include "COMMON.MPI"
16 !-----------------------------------------------------------------------------
19 real(kind=8),dimension(:,:),allocatable :: ww_all !(max_ene,max_parm) ! max_eneW
20 real(kind=8),dimension(:),allocatable :: vbldp0_all,akp_all !(max_parm)
21 real(kind=8),dimension(:,:,:),allocatable :: vbldsc0_all,&
22 aksc_all,abond0_all !(maxbondterm,ntyp,max_parm)
23 real(kind=8),dimension(:,:),allocatable :: a0thet_all !(-ntyp:ntyp,max_parm)
24 real(kind=8),dimension(:,:,:,:,:),allocatable :: athet_all,&
25 bthet_all !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
26 real(kind=8),dimension(:,:,:),allocatable :: polthet_all !(0:3,-ntyp:ntyp,max_parm)
27 real(kind=8),dimension(:,:,:),allocatable :: gthet_all !(3,-ntyp:ntyp,max_parm)
28 real(kind=8),dimension(:,:),allocatable :: theta0_all,&
29 sig0_all,sigc0_all !(-ntyp:ntyp,max_parm)
30 real(kind=8),dimension(:,:,:,:,:),allocatable :: aa0thet_all
31 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
32 real(kind=8),dimension(:,:,:,:,:,:),allocatable :: aathet_all
33 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
34 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: bbthet_all,&
35 ccthet_all,ddthet_all,eethet_all !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
36 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
37 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet_all1,&
38 ggthet_all1,ffthet_all2,ggthet_all2 !(maxdouble,maxdouble,maxtheterm3,
39 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
40 real(kind=8),dimension(:,:),allocatable :: dsc_all,dsc0_all !(ntyp1,max_parm)
41 real(kind=8),dimension(:,:,:),allocatable :: bsc_all !(maxlob,ntyp,max_parm)
42 real(kind=8),dimension(:,:,:,:),allocatable :: censc_all !(3,maxlob,-ntyp:ntyp,max_parm)
43 real(kind=8),dimension(:,:,:,:,:),allocatable :: gaussc_all !(3,3,maxlob,-ntyp:ntyp,max_parm)
44 real(kind=8),dimension(:,:,:),allocatable :: sc_parmin_all !(65,ntyp,max_parm)
45 real(kind=8),dimension(:,:,:,:),allocatable :: v0_all
46 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
47 real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_all,&
48 v2_all !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
49 real(kind=8),dimension(:,:,:,:),allocatable :: vlor1_all,&
50 vlor2_all,vlor3_all !(maxlor,maxtor,maxtor,max_parm)
51 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v1c_all,&
52 v1s_all !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
53 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2c_all
54 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
55 real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: v2s_all
56 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
57 real(kind=8),dimension(:,:,:),allocatable :: b1_all,b2_all !(2,-maxtor:maxtor,max_parm)
58 real(kind=8),dimension(:,:,:,:),allocatable :: cc_all,dd_all,&
59 ee_all !(2,2,-maxtor:maxtor,max_parm)
60 real(kind=8),dimension(:,:,:,:),allocatable :: ctilde_all,&
61 dtilde_all !(2,2,-maxtor:maxtor,max_parm)
62 real(kind=8),dimension(:,:,:),allocatable :: b1tilde_all !(2,-maxtor:maxtor,max_parm)
63 real(kind=8),dimension(:,:,:),allocatable :: app_all,bpp_all,&
64 ael6_all,ael3_all !(2,2,max_parm)
65 real(kind=8),dimension(:,:,:),allocatable :: aad_all,&
66 bad_all !(ntyp,2,max_parm)
67 real(kind=8),dimension(:,:,:),allocatable :: aa_all,bb_all,&
68 augm_all,eps_all,sigma_all,r0_all,chi_all !(ntyp,ntyp,max_parm)
69 real(kind=8),dimension(:,:),allocatable :: chip_all,alp_all !(ntyp,max_parm)
70 real(kind=8),dimension(:),allocatable :: ebr_all,d0cm_all,&
71 akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all !(max_parm)
72 real(kind=8),dimension(:,:,:,:,:),allocatable :: v1sccor_all,&
73 v2sccor_all !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
74 integer,dimension(:,:),allocatable :: nlob_all !(ntyp1,max_parm)
75 integer,dimension(:,:,:,:),allocatable :: nlor_all,&
76 nterm_all !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
77 integer,dimension(:,:,:,:,:),allocatable :: ntermd1_all,&
78 ntermd2_all !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
79 integer,dimension(:,:),allocatable :: nbondterm_all !(ntyp,max_parm)
80 integer,dimension(:,:),allocatable :: ithetyp_all !(-ntyp1:ntyp1,max_parm)
81 integer,dimension(:),allocatable :: nthetyp_all,ntheterm_all,&
82 ntheterm2_all,ntheterm3_all,nsingle_all,ndouble_all,&
83 nntheterm_all !(max_parm)
84 integer,dimension(:,:,:),allocatable :: nterm_sccor_all !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
85 !-----------------------------------------------------------------------------
88 !-----------------------------------------------------------------------------
90 !-----------------------------------------------------------------------------
91 subroutine enecalc(islice,*)
94 use control_data, only:indpdb
95 use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,&
96 vbld,rad2deg,dc_norm,dc,vbld_inv
97 use io_base, only:gyrate!,briefout
98 use geometry, only:int_from_cart1,returnbox
99 use io_wham, only:pdboutW
100 use io_database, only:opentmp
101 use conform_compar, only:qwolynes,rmsnat
102 ! include "DIMENSIONS"
103 ! include "DIMENSIONS.ZSCOPT"
104 ! include "DIMENSIONS.FREE"
108 ! include "COMMON.MPI"
110 ! include "COMMON.CHAIN"
111 ! include "COMMON.IOUNITS"
112 ! include "COMMON.PROTFILES"
113 ! include "COMMON.NAMES"
114 ! include "COMMON.VAR"
115 ! include "COMMON.SBRIDGE"
116 ! include "COMMON.GEO"
117 ! include "COMMON.FFIELD"
118 ! include "COMMON.ENEPS"
119 ! include "COMMON.LOCAL"
120 ! include "COMMON.WEIGHTS"
121 ! include "COMMON.INTERACT"
122 ! include "COMMON.FREE"
123 ! include "COMMON.ENERGIES"
124 ! include "COMMON.CONTROL"
125 ! include "COMMON.TORCNSTR"
128 integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
130 character(len=64) :: nazwa
131 character(len=80) :: bxname
132 character(len=3) :: liczba
133 !el real(kind=8) :: qwolynes
134 !el external qwolynes
135 integer :: errmsg_count,maxerrmsg_count=100
136 !el real(kind=8) :: rmsnat,gyrate
137 !el external rmsnat,gyrate
138 real(kind=8) :: tole=1.0d-1
139 integer i,itj,ii,iii,j,k,l,licz
140 integer ir,ib,ipar,iparm
142 real(kind=4) :: csingle(3,nres*2)
143 real(kind=8) :: energ
145 !el integer ilen,iroof
146 !el external ilen,iroof
147 real(kind=8) :: energia(0:n_ene),rmsdev,efree,eini
148 !el real(kind=8) :: energia(0:max_ene),rmsdev,efree,eini
149 real(kind=8) :: fT(6),quot,quotl,kfacl,kfac=2.4d0,T0=3.0d2
151 integer :: snk_p(MaxR,MaxT_h,nParmSet)!Max_parm)
153 character(len=64) :: bprotfile_temp
156 integer,dimension(0:nprocs) :: scount_
157 !el real(kind=8) :: rmsnat
159 rescale_mode=rescale_modeW
161 call opentmp(islice,ientout,bprotfile_temp)
167 write (iout,*) "enecalc: nparmset ",nparmset
176 do i=indstart(me1),indend(me1)
177 write(iout,*)"enecalc_ i indstart",i,indstart(me1),indend(me1)
187 write(iout,*)"enecalc_ i ntot",i,ntot
189 read(ientout,rec=i,err=101) &
190 ((csingle(l,k),l=1,3),k=1,nres),&
191 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
192 nss,(ihpb(k),jhpb(k),k=1,nss),&
193 eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
195 !write(iout,*)"co wczytuje"
196 ! write(iout,*)((csingle(l,k),l=1,3),k=1,nres),&
197 ! ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
198 ! nss,(ihpb(k),jhpb(k),k=1,nss),&
199 ! eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
202 !write(iout,*)"ipar",ib,ipar,1.0d0/(beta_h(ib,ipar)*1.987D-3)
203 if (indpdb.gt.0) then
211 c(l,k+nres)=csingle(l,k+nres)
222 csingle(l,k+nres)=c(l,k+nres)
225 anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
226 q(nQ+1,iii+1)=rmsnat(iii+1)
228 q(nQ+2,iii+1)=gyrate(iii+1)
229 ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
230 ! fT=T0*beta_h(ib,ipar)*1.987D-3
231 ! ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
232 ! EL start old rescale
233 ! if (rescale_modeW.eq.1) then
234 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
236 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
237 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
238 !#elif defined(FUNCT)
248 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
250 ! else if (rescale_modeW.eq.2) then
251 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
253 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
254 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
255 !#elif defined(FUNCT)
263 ! fT(l)=1.12692801104297249644d0/ &
264 ! dlog(dexp(quotl)+dexp(-quotl))
266 ! else if (rescale_modeW.eq.0) then
271 ! write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
277 ! write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
278 ! & " kfac",kfac,"quot",quot," fT",fT
280 write(iout,*)"weights"
281 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
282 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
291 call int_from_cart1(.false.)
294 ! call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
297 write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
298 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
299 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
302 call restore_parm(iparm)
303 call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
305 write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
306 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
307 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
310 ! call etotal(energia(0),fT)
311 call etotal(energia(0))
312 !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
314 !write(iout,*)k,"c=",(c(l,k),l=1,3)
315 !write(iout,*)k,"dc=",(dc(l,k),l=1,3)
316 !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
319 !write(iout,*)k,"vbld=",vbld(k)
320 !write(iout,*)k,"vbld_inv=",vbld_inv(k)
323 !write(iout,*)"energia",(energia(j),j=0,n_ene)
324 !write(iout,*)"enerprint tuz po call etotal"
325 call enerprint(energia(0))
327 write (iout,*) "Conformation",i
328 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
329 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
330 ! call enerprint(energia(0),fT)
331 call enerprint(energia(0))
332 write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
333 write (iout,*) "ftors",ftors
334 !el call briefout(i,energia(0))
335 temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
336 write (iout,*) "temp", temp
337 call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
339 if (energia(0).ge.1.0d20) then
340 write (iout,*) "NaNs detected in some of the energy",&
341 " components for conformation",ii+1
342 write (iout,*) "The Cartesian geometry is:"
343 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
344 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
345 write (iout,*) "The internal geometry is:"
347 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
348 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
349 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
350 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
351 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
352 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
353 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
354 write (iout,*) "The components of the energy are:"
355 ! call enerprint(energia(0),fT)
356 call enerprint(energia(0))
358 "This conformation WILL NOT be added to the database."
363 if (ipar.eq.iparm) write (iout,*) i,iparm,&
364 1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
366 if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
367 dabs(eini-energia(0)).gt.tole) then
368 if (errmsg_count.le.maxerrmsg_count) then
369 write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
370 "Warning: energy differs remarkably from ",&
371 " the value read in: ",energia(0),eini," point",&
372 iii+1,indstart(me1)+iii," T",&
373 1.0d0/(1.987D-3*beta_h(ib,ipar))
375 call pdboutW(indstart(me1)+iii,&
376 1.0d0/(1.987D-3*beta_h(ib,ipar)),&
377 energia(0),eini,0.0d0,0.0d0)
378 write (iout,*) "The Cartesian geometry is:"
379 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
380 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
381 write (iout,*) "The internal geometry is:"
383 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
384 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
385 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
386 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
387 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
388 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
389 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
390 call enerprint(energia(0))
391 ! call enerprint(energia(0),fT)
392 errmsg_count=errmsg_count+1
393 if (errmsg_count.gt.maxerrmsg_count) &
394 write (iout,*) "Too many warning messages"
395 if (einicheck.gt.1) then
396 write (iout,*) "Calculation stopped."
399 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
406 potE(iii+1,iparm)=energia(0)
408 enetb(k,iii+1,iparm)=energia(k)
410 ! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
411 ! write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
413 write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
414 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
415 ! call enerprint(energia(0),fT)
418 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
419 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
420 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
421 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
422 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
423 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
424 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
425 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
426 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
427 write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
428 write (iout,'(f10.5,i10)') rmsdev,iscor
429 ! call enerprint(energia(0),fT)
430 call enerprint(energia(0))
431 call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
438 if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
439 write (ientout,rec=iii) &
440 ((csingle(l,k),l=1,3),k=1,nres),&
441 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
442 nss,(ihpb(k),jhpb(k),k=1,nss),&
443 potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
444 ! write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
446 if (separate_parset) then
447 snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
449 snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
451 ! write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
452 ! & " snk",snk_p(iR,ib,ipar)
454 snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
460 write (iout,*) "Me",me," scount",scount(me)
462 ! Master gathers updated numbers of conformations written by all procs.
463 call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
464 MPI_INTEGER, WHAM_COMM, IERROR)
468 indstart(i)=indend(i-1)+1
469 indend(i)=indstart(i)+scount_(i)-1
472 write (iout,*) "Revised conformation counts"
474 write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
475 "Processor",i," indstart",indstart(i),&
476 " indend",indend(i)," count",scount_(i)
479 call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
480 MaxR*MaxT_h*nParmSet,&
481 MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
487 stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
491 write (iout,*) "Revised SNK"
494 write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
495 iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
496 (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
497 write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
500 write (iout,'("Total",i10)') stot(islice)
506 101 write (iout,*) "Error in scratchfile."
510 end subroutine enecalc
511 !------------------------------------------------------------------------------
512 logical function conf_check(ii,iprint)
514 use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
515 use geometry, only:int_from_cart1
516 ! include "DIMENSIONS"
517 ! include "DIMENSIONS.ZSCOPT"
518 ! include "DIMENSIONS.FREE"
522 ! include "COMMON.MPI"
524 ! include "COMMON.CHAIN"
525 ! include "COMMON.IOUNITS"
526 ! include "COMMON.PROTFILES"
527 ! include "COMMON.NAMES"
528 ! include "COMMON.VAR"
529 ! include "COMMON.SBRIDGE"
530 ! include "COMMON.GEO"
531 ! include "COMMON.FFIELD"
532 ! include "COMMON.ENEPS"
533 ! include "COMMON.LOCAL"
534 ! include "COMMON.WEIGHTS"
535 ! include "COMMON.INTERACT"
536 ! include "COMMON.FREE"
537 ! include "COMMON.ENERGIES"
538 ! include "COMMON.CONTROL"
539 ! include "COMMON.TORCNSTR"
542 integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
544 integer :: j,k,l,ii,itj,iprint,mnum,mnum1
545 if (.not.check_conf) then
549 call int_from_cart1(.false.)
554 if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
555 (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
557 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
558 " for conformation",ii
559 if (iprint.gt.1) then
560 write (iout,*) "The Cartesian geometry is:"
561 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
562 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
563 write (iout,*) "The internal geometry is:"
564 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
565 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
566 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
567 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
568 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
569 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
571 if (iprint.gt.0) write (iout,*) &
572 "This conformation WILL NOT be added to the database."
581 if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
582 (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
584 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
585 " for conformation",ii
586 if (iprint.gt.1) then
587 write (iout,*) "The Cartesian geometry is:"
588 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
589 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
590 write (iout,*) "The internal geometry is:"
591 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
592 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
593 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
594 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
595 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
596 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
598 if (iprint.gt.0) write (iout,*) &
599 "This conformation WILL NOT be added to the database."
605 if (theta(j).le.0.0d0) then
607 write (iout,*) "Zero theta angle(s) in conformation",ii
608 if (iprint.gt.1) then
609 write (iout,*) "The Cartesian geometry is:"
610 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
611 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
612 write (iout,*) "The internal geometry is:"
613 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
614 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
615 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
616 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
617 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
618 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
620 if (iprint.gt.0) write (iout,*) &
621 "This conformation WILL NOT be added to the database."
625 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
628 ! write (iout,*) "conf_check passed",ii
630 end function conf_check
631 !-----------------------------------------------------------------------------
633 !-----------------------------------------------------------------------------
634 subroutine store_parm(iparm)
636 ! Store parameters of set IPARM
637 ! valence angles and the side chains and energy parameters.
640 ! include 'DIMENSIONS'
641 ! include 'DIMENSIONS.ZSCOPT'
642 ! include 'DIMENSIONS.FREE'
643 ! include 'COMMON.IOUNITS'
644 ! include 'COMMON.CHAIN'
645 ! include 'COMMON.INTERACT'
646 ! include 'COMMON.GEO'
647 ! include 'COMMON.LOCAL'
648 ! include 'COMMON.TORSION'
649 ! include 'COMMON.FFIELD'
650 ! include 'COMMON.NAMES'
651 ! include 'COMMON.SBRIDGE'
652 ! include 'COMMON.SCROT'
653 ! include 'COMMON.SCCOR'
654 ! include 'COMMON.ALLPARM'
655 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
657 call alloc_enecalc_arrays(iparm)
658 !el allocate(ww_all(n_ene,iparm))
662 ww_all(3,iparm)=welec
663 ww_all(4,iparm)=wcorr
664 ww_all(5,iparm)=wcorr5
665 ww_all(6,iparm)=wcorr6
666 ww_all(7,iparm)=wel_loc
667 ww_all(8,iparm)=wturn3
668 ww_all(9,iparm)=wturn4
669 ww_all(10,iparm)=wturn6
670 ww_all(11,iparm)=wang
671 ww_all(12,iparm)=wscloc
672 ww_all(13,iparm)=wtor
673 ww_all(14,iparm)=wtor_d
674 ww_all(15,iparm)=wstrain
675 ww_all(16,iparm)=wvdwpp
676 ww_all(17,iparm)=wbond
677 ww_all(19,iparm)=wsccor
678 ! Store bond parameters
679 vbldp0_all(iparm)=vbldp0
682 nbondterm_all(i,iparm)=nbondterm(i)
684 vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
685 aksc_all(j,i,iparm)=aksc(j,i)
686 abond0_all(j,i,iparm)=abond0(j,i)
689 ! Store bond angle parameters
692 a0thet_all(i,iparm)=a0thet(i)
696 athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
697 bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
702 polthet_all(j,i,iparm)=polthet(j,i)
705 gthet_all(j,i,iparm)=gthet(j,i)
707 theta0_all(i,iparm)=theta0(i)
708 sig0_all(i,iparm)=sig0(i)
709 sigc0_all(i,iparm)=sigc0(i)
712 nthetyp_all(iparm)=nthetyp
713 ntheterm_all(iparm)=ntheterm
714 ntheterm2_all(iparm)=ntheterm2
715 ntheterm3_all(iparm)=ntheterm3
716 nsingle_all(iparm)=nsingle
717 ndouble_all(iparm)=ndouble
718 nntheterm_all(iparm)=nntheterm
720 ithetyp_all(i,iparm)=ithetyp(i)
723 do i=-nthetyp-1,nthetyp+1
724 do j=-nthetyp-1,nthetyp+1
725 do k=-nthetyp-1,nthetyp+1
726 aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
728 aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
732 bbthet_all(m,l,i,j,k,iblock,iparm)= &
733 bbthet(m,l,i,j,k,iblock)
734 ccthet_all(m,l,i,j,k,iblock,iparm)= &
735 ccthet(m,l,i,j,k,iblock)
736 ddthet_all(m,l,i,j,k,iblock,iparm)= &
737 ddthet(m,l,i,j,k,iblock)
738 eethet_all(m,l,i,j,k,iblock,iparm)= &
739 eethet(m,l,i,j,k,iblock)
745 if (iblock.eq.1) then
746 ffthet_all1(mm,m,l,i,j,k,iparm)=&
747 ffthet(mm,m,l,i,j,k,iblock)
748 ggthet_all1(mm,m,l,i,j,k,iparm)=&
749 ggthet(mm,m,l,i,j,k,iblock)
751 ffthet_all2(mm,m,l,i,j,k,iparm)=&
752 ffthet(mm,m,l,i,j,k,iblock)
753 ggthet_all2(mm,m,l,i,j,k,iparm)=&
754 ggthet(mm,m,l,i,j,k,iblock)
765 ! Store the sidechain rotamer parameters
768 !! write (iout,*) i,"storeparm1"
770 nlob_all(iii,iparm)=nlob(iii)
772 bsc_all(j,iii,iparm)=bsc(j,iii)
774 censc_all(k,j,i,iparm)=censc(k,j,i)
778 gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
786 sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
790 ! Store the torsional parameters
792 do i=-ntortyp+1,ntortyp-1
793 do j=-ntortyp+1,ntortyp-1
794 v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
795 nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
796 nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
797 do k=1,nterm(i,j,iblock)
798 v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
799 v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
801 do k=1,nlor(i,j,iblock)
802 vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
803 vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
804 vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
809 ! Store the double torsional parameters
811 do i=-ntortyp+1,ntortyp-1
812 do j=-ntortyp+1,ntortyp-1
813 do k=-ntortyp+1,ntortyp-1
814 ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
815 ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
816 do l=1,ntermd_1(i,j,k,iblock)
817 v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
818 v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
819 v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
820 v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
822 do l=1,ntermd_2(i,j,k,iblock)
823 do m=1,ntermd_2(i,j,k,iblock)
824 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
831 ! Store parameters of the cumulants
832 do i=-nloctyp,nloctyp
834 b1_all(j,i,iparm)=b1(j,i)
835 b1tilde_all(j,i,iparm)=b1tilde(j,i)
836 b2_all(j,i,iparm)=b2(j,i)
840 cc_all(k,j,i,iparm)=cc(k,j,i)
841 ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
842 dd_all(k,j,i,iparm)=dd(k,j,i)
843 dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
844 ee_all(k,j,i,iparm)=ee(k,j,i)
848 ! Store the parameters of electrostatic interactions
851 app_all(j,i,iparm)=app(j,i)
852 bpp_all(j,i,iparm)=bpp(j,i)
853 ael6_all(j,i,iparm)=ael6(j,i)
854 ael3_all(j,i,iparm)=ael3(j,i)
857 ! Store sidechain parameters
860 aa_all(j,i,iparm)=aa_aq(j,i)
861 bb_all(j,i,iparm)=bb_aq(j,i)
862 r0_all(j,i,iparm)=r0(j,i)
863 sigma_all(j,i,iparm)=sigma(j,i)
864 chi_all(j,i,iparm)=chi(j,i)
865 augm_all(j,i,iparm)=augm(j,i)
866 eps_all(j,i,iparm)=eps(j,i)
870 chip_all(i,iparm)=chip(i)
871 alp_all(i,iparm)=alp(i)
873 ! Store the SCp parameters
876 aad_all(i,j,iparm)=aad(i,j)
877 bad_all(i,j,iparm)=bad(i,j)
880 ! Store disulfide-bond parameters
889 ! Store SC-backbone correlation parameters
890 do i=-nsccortyp,nsccortyp
891 do j=-nsccortyp,nsccortyp
893 nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
897 do k=1,nterm_sccor(j,i)
898 v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
899 v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
904 write(iout,*)"end of store_parm"
906 end subroutine store_parm
907 !--------------------------------------------------------------------------
908 subroutine restore_parm(iparm)
910 ! Store parameters of set IPARM
911 ! valence angles and the side chains and energy parameters.
914 ! include 'DIMENSIONS'
915 ! include 'DIMENSIONS.ZSCOPT'
916 ! include 'DIMENSIONS.FREE'
917 ! include 'COMMON.IOUNITS'
918 ! include 'COMMON.CHAIN'
919 ! include 'COMMON.INTERACT'
920 ! include 'COMMON.GEO'
921 ! include 'COMMON.LOCAL'
922 ! include 'COMMON.TORSION'
923 ! include 'COMMON.FFIELD'
924 ! include 'COMMON.NAMES'
925 ! include 'COMMON.SBRIDGE'
926 ! include 'COMMON.SCROT'
927 ! include 'COMMON.SCCOR'
928 ! include 'COMMON.ALLPARM'
929 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
934 welec=ww_all(3,iparm)
935 wcorr=ww_all(4,iparm)
936 wcorr5=ww_all(5,iparm)
937 wcorr6=ww_all(6,iparm)
938 wel_loc=ww_all(7,iparm)
939 wturn3=ww_all(8,iparm)
940 wturn4=ww_all(9,iparm)
941 wturn6=ww_all(10,iparm)
942 wang=ww_all(11,iparm)
943 wscloc=ww_all(12,iparm)
944 wtor=ww_all(13,iparm)
945 wtor_d=ww_all(14,iparm)
946 wstrain=ww_all(15,iparm)
947 wvdwpp=ww_all(16,iparm)
948 wbond=ww_all(17,iparm)
949 wsccor=ww_all(19,iparm)
950 ! Restore bond parameters
951 vbldp0=vbldp0_all(iparm)
954 nbondterm(i)=nbondterm_all(i,iparm)
956 vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
957 aksc(j,i)=aksc_all(j,i,iparm)
958 abond0(j,i)=abond0_all(j,i,iparm)
961 ! Restore bond angle parameters
964 a0thet(i)=a0thet_all(i,iparm)
968 athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
969 bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
974 polthet(j,i)=polthet_all(j,i,iparm)
977 gthet(j,i)=gthet_all(j,i,iparm)
979 theta0(i)=theta0_all(i,iparm)
980 sig0(i)=sig0_all(i,iparm)
981 sigc0(i)=sigc0_all(i,iparm)
984 nthetyp=nthetyp_all(iparm)
985 ntheterm=ntheterm_all(iparm)
986 ntheterm2=ntheterm2_all(iparm)
987 ntheterm3=ntheterm3_all(iparm)
988 nsingle=nsingle_all(iparm)
989 ndouble=ndouble_all(iparm)
990 nntheterm=nntheterm_all(iparm)
992 ithetyp(i)=ithetyp_all(i,iparm)
995 do i=-nthetyp-1,nthetyp+1
996 do j=-nthetyp-1,nthetyp+1
997 do k=-nthetyp-1,nthetyp+1
998 aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
1000 aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
1004 bbthet(m,l,i,j,k,iblock)= &
1005 bbthet_all(m,l,i,j,k,iblock,iparm)
1006 ccthet(m,l,i,j,k,iblock)= &
1007 ccthet_all(m,l,i,j,k,iblock,iparm)
1008 ddthet(m,l,i,j,k,iblock)= &
1009 ddthet_all(m,l,i,j,k,iblock,iparm)
1010 eethet(m,l,i,j,k,iblock)= &
1011 eethet_all(m,l,i,j,k,iblock,iparm)
1017 if (iblock.eq.1) then
1018 ffthet(mm,m,l,i,j,k,iblock)= &
1019 ffthet_all1(mm,m,l,i,j,k,iparm)
1020 ggthet(mm,m,l,i,j,k,iblock)= &
1021 ggthet_all1(mm,m,l,i,j,k,iparm)
1023 ffthet(mm,m,l,i,j,k,iblock)= &
1024 ffthet_all2(mm,m,l,i,j,k,iparm)
1025 ggthet(mm,m,l,i,j,k,iblock)= &
1026 ggthet_all2(mm,m,l,i,j,k,iparm)
1036 ! Restore the sidechain rotamer parameters
1041 nlob(iii)=nlob_all(iii,iparm)
1043 bsc(j,iii)=bsc_all(j,iii,iparm)
1045 censc(k,j,i)=censc_all(k,j,i,iparm)
1049 gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1057 sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1061 ! Restore the torsional parameters
1063 do i=-ntortyp+1,ntortyp-1
1064 do j=-ntortyp+1,ntortyp-1
1065 v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1066 nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1067 nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1068 do k=1,nterm(i,j,iblock)
1069 v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1070 v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1072 do k=1,nlor(i,j,iblock)
1073 vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1074 vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1075 vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1080 ! Restore the double torsional parameters
1082 do i=-ntortyp+1,ntortyp-1
1083 do j=-ntortyp+1,ntortyp-1
1084 do k=-ntortyp+1,ntortyp-1
1085 ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1086 ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1087 do l=1,ntermd_1(i,j,k,iblock)
1088 v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1089 v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1090 v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1091 v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1093 do l=1,ntermd_2(i,j,k,iblock)
1094 do m=1,ntermd_2(i,j,k,iblock)
1095 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1102 ! Restore parameters of the cumulants
1103 do i=-nloctyp,nloctyp
1105 b1(j,i)=b1_all(j,i,iparm)
1106 b1tilde(j,i)=b1tilde_all(j,i,iparm)
1107 b2(j,i)=b2_all(j,i,iparm)
1111 cc(k,j,i)=cc_all(k,j,i,iparm)
1112 ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1113 dd(k,j,i)=dd_all(k,j,i,iparm)
1114 dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1115 ee(k,j,i)=ee_all(k,j,i,iparm)
1119 ! Restore the parameters of electrostatic interactions
1122 app(j,i)=app_all(j,i,iparm)
1123 bpp(j,i)=bpp_all(j,i,iparm)
1124 ael6(j,i)=ael6_all(j,i,iparm)
1125 ael3(j,i)=ael3_all(j,i,iparm)
1128 ! Restore sidechain parameters
1131 aa_aq(j,i)=aa_all(j,i,iparm)
1132 bb_aq(j,i)=bb_all(j,i,iparm)
1133 r0(j,i)=r0_all(j,i,iparm)
1134 sigma(j,i)=sigma_all(j,i,iparm)
1135 chi(j,i)=chi_all(j,i,iparm)
1136 augm(j,i)=augm_all(j,i,iparm)
1137 eps(j,i)=eps_all(j,i,iparm)
1141 chip(i)=chip_all(i,iparm)
1142 alp(i)=alp_all(i,iparm)
1144 ! Restore the SCp parameters
1147 aad(i,j)=aad_all(i,j,iparm)
1148 bad(i,j)=bad_all(i,j,iparm)
1151 ! Restore disulfide-bond parameters
1153 d0cm=d0cm_all(iparm)
1154 akcm=akcm_all(iparm)
1155 akth=akth_all(iparm)
1156 akct=akct_all(iparm)
1157 v1ss=v1ss_all(iparm)
1158 v2ss=v2ss_all(iparm)
1159 v3ss=v3ss_all(iparm)
1160 ! Restore SC-backbone correlation parameters
1161 do i=-nsccortyp,nsccortyp
1162 do j=-nsccortyp,nsccortyp
1164 nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1166 do k=1,nterm_sccor(j,i)
1167 v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1168 v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1174 end subroutine restore_parm
1175 !--------------------------------------------------------------------------
1177 !--------------------------------------------------------------------------
1178 subroutine make_ensembles(islice,*)
1179 ! construct the conformational ensembles at REMD temperatures
1180 use geometry_data, only:c
1181 use io_base, only:ilen
1182 use io_wham, only:pdboutW
1183 use geometry, only:returnbox
1185 ! include "DIMENSIONS"
1186 ! include "DIMENSIONS.ZSCOPT"
1187 ! include "DIMENSIONS.FREE"
1190 ! include "COMMON.MPI"
1191 integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1193 ! include "COMMON.IOUNITS"
1194 ! include "COMMON.CONTROL"
1195 ! include "COMMON.FREE"
1196 ! include "COMMON.ENERGIES"
1197 ! include "COMMON.FFIELD"
1198 ! include "COMMON.INTERACT"
1199 ! include "COMMON.SBRIDGE"
1200 ! include "COMMON.CHAIN"
1201 ! include "COMMON.PROTFILES"
1202 ! include "COMMON.PROT"
1203 real(kind=4) :: csingle(3,nres*2)
1204 real(kind=8),dimension(6) :: fT,fTprim,fTbis
1205 real(kind=8) :: quot,quotl1,quotl,kfacl,&
1206 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1207 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1208 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1209 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
1210 ecation_prot, ecationcation
1211 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1212 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1213 character(len=80) :: bxname
1214 character(len=2) :: licz1,licz2
1215 character(len=3) :: licz3,licz4
1216 character(len=5) :: ctemper
1219 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1220 real(kind=8) :: enepot(MaxStr)
1221 integer :: iperm(MaxStr)
1223 integer,dimension(0:nprocs) :: scount_
1226 if (me.eq.Master) then
1228 write (licz2,'(bz,i2.2)') islice
1229 if (nslice.eq.1) then
1230 if (.not.separate_parset) then
1231 bxname = prefix(:ilen(prefix))//".bx"
1233 write (licz3,'(bz,i3.3)') myparm
1234 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1237 if (.not.separate_parset) then
1238 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1240 write (licz3,'(bz,i3.3)') myparm
1241 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1242 "_slice_"//licz2//".bx"
1245 open (ientout,file=bxname,status="unknown",&
1246 form="unformatted",access="direct",recl=lenrec1)
1251 if (iparm.ne.iparmprint) exit
1252 call restore_parm(iparm)
1255 write (iout,*) "iparm",iparm," ib",ib
1257 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1258 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1265 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1267 !el old rescale weights
1269 ! if (rescale_mode.eq.1) then
1270 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1272 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1273 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1274 #elif defined(FUNCT)
1285 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1287 ! else if (rescale_mode.eq.2) then
1288 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1289 !#if defined(FUNCTH)
1290 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1291 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1292 !#elif defined(FUNCT)
1300 ! fT(l)=1.12692801104297249644d0/ &
1301 ! dlog(dexp(quotl)+dexp(-quotl))
1303 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1304 ! else if (rescale_mode.eq.0) then
1310 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1314 ! el end old rescale weihgts
1315 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1322 evdw=enetb(1,i,iparm)
1323 ! evdw_t=enetb(21,i,iparm)
1324 evdw_t=enetb(20,i,iparm)
1326 ! evdw2_14=enetb(17,i,iparm)
1327 evdw2_14=enetb(18,i,iparm)
1328 evdw2=enetb(2,i,iparm)+evdw2_14
1330 evdw2=enetb(2,i,iparm)
1334 ees=enetb(3,i,iparm)
1335 evdw1=enetb(16,i,iparm)
1337 ees=enetb(3,i,iparm)
1340 ecorr=enetb(4,i,iparm)
1341 ecorr5=enetb(5,i,iparm)
1342 ecorr6=enetb(6,i,iparm)
1343 eel_loc=enetb(7,i,iparm)
1344 eello_turn3=enetb(8,i,iparm)
1345 eello_turn4=enetb(9,i,iparm)
1346 eello_turn6=enetb(10,i,iparm)
1347 ebe=enetb(11,i,iparm)
1348 escloc=enetb(12,i,iparm)
1349 etors=enetb(13,i,iparm)
1350 etors_d=enetb(14,i,iparm)
1351 ehpb=enetb(15,i,iparm)
1353 estr=enetb(17,i,iparm)
1354 ! estr=enetb(18,i,iparm)
1355 ! esccor=enetb(19,i,iparm)
1356 esccor=enetb(21,i,iparm)
1357 ! edihcnstr=enetb(20,i,iparm)
1358 edihcnstr=enetb(19,i,iparm)
1359 ecationcation=enetb(42,i,iparm)
1360 ecation_prot=enetb(41,i,iparm)
1363 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1365 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1366 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1367 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1368 ! +ft(2)*wturn3*eello_turn3 &
1369 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1370 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1373 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1374 ! +ft(1)*welec*(ees+evdw1) &
1375 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1376 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1377 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1378 ! +ft(2)*wturn3*eello_turn3 &
1379 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1380 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1385 etot=wsc*evdw+wscp*evdw2+welec*ees &
1387 +wang*ebe+wtor*etors+wscloc*escloc &
1388 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1389 +wcorr6*ecorr6+wturn4*eello_turn4 &
1390 +wturn3*eello_turn3 &
1391 +wturn6*eello_turn6+wel_loc*eel_loc &
1392 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1393 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1395 etot=wsc*evdw+wscp*evdw2 &
1396 +welec*(ees+evdw1) &
1397 +wang*ebe+wtor*etors+wscloc*escloc &
1398 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1399 +wcorr6*ecorr6+wturn4*eello_turn4 &
1400 +wturn3*eello_turn3 &
1401 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1402 +wtor_d*etors_d+wsccor*esccor &
1403 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1408 beta_h(ib,iparm)*etot-entfac(i)
1411 write (iout,*) i,indstart(me)+i-1,ib,&
1412 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1413 -entfac(i),Fdimless(i)
1416 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1422 Fdimless_(i)=Fdimless(i)
1424 call MPI_Gatherv(Fdimless_(1),scount(me),&
1425 MPI_REAL,Fdimless(1),&
1426 scount_(0),idispl(0),MPI_REAL,Master,&
1429 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1430 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1431 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1433 call MPI_Gatherv(entfac(1),scount(me),&
1434 MPI_DOUBLE_PRECISION,entfac(1),&
1435 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1438 if (me.eq.Master) then
1440 write (iout,*) "The FDIMLESS array before sorting"
1442 write (iout,*) i,fdimless(i)
1449 call mysort1(ntot(islice),Fdimless,iperm)
1451 write (iout,*) "The FDIMLESS array after sorting"
1453 write (iout,*) i,iperm(i),fdimless(i)
1458 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1460 ! write (iout,*) "qfree",qfree
1463 do i=1,min0(ntot(islice),ensembles)
1464 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1466 write (iout,*) i,ib,beta_h(ib,iparm),&
1467 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1468 potE(iperm(i),iparm),&
1469 -entfac(iperm(i)),fdimless(i),sumprob
1471 if (sumprob.gt.0.99d0) goto 122
1477 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1479 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1484 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1487 if (iproc.ge.nprocs) then
1488 write (iout,*) "Fatal error: processor out of range",iproc
1493 close (ientout,status="delete")
1497 ik=ii-indstart(iproc)+1
1498 if (iproc.ne.Master) then
1499 if (me.eq.iproc) then
1501 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1502 " energy",potE(ik,iparm)
1504 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1505 Master,i,WHAM_COMM,IERROR)
1506 else if (me.eq.Master) then
1507 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1508 WHAM_COMM,STATUS,IERROR)
1510 else if (me.eq.Master) then
1511 enepot(i)=potE(ik,iparm)
1516 enepot(i)=potE(iperm(i),iparm)
1520 if (me.eq.Master) then
1522 write(licz3,'(bz,i3.3)') iparm
1523 write(licz2,'(bz,i2.2)') islice
1524 if (temper.lt.100.0d0) then
1525 write(ctemper,'(f3.0)') temper
1526 else if (temper.lt.1000.0) then
1527 write (ctemper,'(f4.0)') temper
1529 write (ctemper,'(f5.0)') temper
1531 if (nparmset.eq.1) then
1532 if (separate_parset) then
1533 write(licz4,'(bz,i3.3)') myparm
1534 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1536 pdbname=prefix(:ilen(prefix))
1539 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1541 if (nslice.eq.1) then
1542 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1543 ctemper(:ilen(ctemper))//"pdb"
1545 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1546 ctemper(:ilen(ctemper))//"pdb"
1548 open(ipdb,file=pdbname)
1550 read (ientout,rec=iperm(i)) &
1551 ((csingle(l,k),l=1,3),k=1,nres),&
1552 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1553 nss,(ihpb(k),jhpb(k),k=1,nss),&
1554 eini,efree,rmsdev,iscor
1568 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1578 close(ientout,status="delete")
1581 scount(i)=scount_(i)
1584 end subroutine make_ensembles
1585 !--------------------------------------------------------------------------
1586 subroutine mysort1(n, x, ipermut)
1588 integer :: i,j,imax,ipm,n
1589 real(kind=4) :: x(n)
1590 integer :: ipermut(n)
1591 real(kind=4) :: xtemp
1596 if (x(j).lt.xtemp) then
1604 ipermut(imax)=ipermut(i)
1608 end subroutine mysort1
1609 !--------------------------------------------------------------------------
1610 subroutine alloc_enecalc_arrays(iparm)
1613 use geometry_data, only:maxlob
1615 !---------------------------
1616 ! COMMON.ENERGIES form wham_data
1618 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1619 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1620 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1621 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1623 ! allocate ENECALC arrays
1624 !---------------------------
1627 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1628 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1629 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1630 aksc_all(maxbondterm,ntyp,iparm),&
1631 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1632 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1633 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1634 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1635 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1636 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1637 allocate(theta0_all(-ntyp:ntyp,iparm),&
1638 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1639 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1640 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1641 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1642 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1643 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1644 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1645 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1646 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1647 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1648 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1649 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1650 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1651 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1652 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1653 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1654 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1655 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1656 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1657 -nthetyp-1:nthetyp+1,iparm))
1658 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1659 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1660 -nthetyp-1:nthetyp+1,iparm))
1661 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1662 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1663 -nthetyp-1:nthetyp+1,iparm))
1664 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1665 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1666 -nthetyp-1:nthetyp+1,iparm))
1667 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1668 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1669 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1670 allocate(bsc_all(maxlob,ntyp,iparm))
1671 !(maxlob,ntyp,max_parm)
1672 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1673 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1674 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1675 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1676 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1677 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1678 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1679 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1680 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1681 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1682 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1683 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1684 -ntortyp:ntortyp,2,iparm))
1685 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1686 -ntortyp:ntortyp,2,iparm))
1687 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1688 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1689 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1690 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1691 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1692 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1693 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1694 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1695 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1696 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1697 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1698 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1699 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1700 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1701 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1702 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1703 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1704 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1705 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1706 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1707 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1708 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1709 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1710 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1711 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1712 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1713 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1714 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1715 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1716 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1717 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1718 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1719 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1720 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1721 -ntortyp:ntortyp,2,iparm))
1722 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1723 -ntortyp:ntortyp,2,iparm))
1724 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1725 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1726 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1727 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1728 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1729 nsingle_all(iparm),&
1730 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1731 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1733 end subroutine alloc_enecalc_arrays
1734 !--------------------------------------------------------------------------
1735 !--------------------------------------------------------------------------