2 !-----------------------------------------------------------------------------
6 use geometry_data, only:nres
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
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)
214 anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
215 q(nQ+1,iii+1)=rmsnat(iii+1)
217 q(nQ+2,iii+1)=gyrate(iii+1)
218 ! write(iout,*)"wczyt",anatemp,q(nQ+2,iii+1) !el
219 ! fT=T0*beta_h(ib,ipar)*1.987D-3
220 ! ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
221 ! EL start old rescale
222 ! if (rescale_modeW.eq.1) then
223 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
225 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
226 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
227 !#elif defined(FUNCT)
237 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
239 ! else if (rescale_modeW.eq.2) then
240 ! quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
242 ! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
243 ! ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
244 !#elif defined(FUNCT)
252 ! fT(l)=1.12692801104297249644d0/ &
253 ! dlog(dexp(quotl)+dexp(-quotl))
255 ! else if (rescale_modeW.eq.0) then
260 ! write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",&
266 ! write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
267 ! & " kfac",kfac,"quot",quot," fT",fT
269 write(iout,*)"weights"
270 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
271 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
280 call int_from_cart1(.false.)
283 ! call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
286 write (iout,*) "before restore w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
287 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
288 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
291 call restore_parm(iparm)
292 call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
294 write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
295 write(iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
296 wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
299 ! call etotal(energia(0),fT)
300 call etotal(energia(0))
301 !write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
303 !write(iout,*)k,"c=",(c(l,k),l=1,3)
304 !write(iout,*)k,"dc=",(dc(l,k),l=1,3)
305 !write(iout,*)k,"dc_norm=",(dc_norm(l,k),l=1,3)
308 !write(iout,*)k,"vbld=",vbld(k)
309 !write(iout,*)k,"vbld_inv=",vbld_inv(k)
312 !write(iout,*)"energia",(energia(j),j=0,n_ene)
313 !write(iout,*)"enerprint tuz po call etotal"
314 call enerprint(energia(0))
316 write (iout,*) "Conformation",i
317 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
318 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
319 ! call enerprint(energia(0),fT)
320 call enerprint(energia(0))
321 write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49)
322 write (iout,*) "ftors",ftors
323 !el call briefout(i,energia(0))
324 temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
325 write (iout,*) "temp", temp
326 call pdboutW(i,temp,energia(0),energia(0),0.0d0,0.0d0)
328 if (energia(0).ge.1.0d20) then
329 write (iout,*) "NaNs detected in some of the energy",&
330 " components for conformation",ii+1
331 write (iout,*) "The Cartesian geometry is:"
332 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
333 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
334 write (iout,*) "The internal geometry is:"
336 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
337 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
338 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
339 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
340 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
341 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
342 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
343 write (iout,*) "The components of the energy are:"
344 ! call enerprint(energia(0),fT)
345 call enerprint(energia(0))
347 "This conformation WILL NOT be added to the database."
352 if (ipar.eq.iparm) write (iout,*) i,iparm,&
353 1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
355 if (ipar.eq.iparm .and. einicheck.gt.0 .and. &
356 dabs(eini-energia(0)).gt.tole) then
357 if (errmsg_count.le.maxerrmsg_count) then
358 write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') &
359 "Warning: energy differs remarkably from ",&
360 " the value read in: ",energia(0),eini," point",&
361 iii+1,indstart(me1)+iii," T",&
362 1.0d0/(1.987D-3*beta_h(ib,ipar))
364 call pdboutW(indstart(me1)+iii,&
365 1.0d0/(1.987D-3*beta_h(ib,ipar)),&
366 energia(0),eini,0.0d0,0.0d0)
367 write (iout,*) "The Cartesian geometry is:"
368 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
369 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
370 write (iout,*) "The internal geometry is:"
372 ! call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
373 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
374 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
375 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
376 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
377 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
378 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
379 call enerprint(energia(0))
380 ! call enerprint(energia(0),fT)
381 errmsg_count=errmsg_count+1
382 if (errmsg_count.gt.maxerrmsg_count) &
383 write (iout,*) "Too many warning messages"
384 if (einicheck.gt.1) then
385 write (iout,*) "Calculation stopped."
388 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
395 potE(iii+1,iparm)=energia(0)
397 enetb(k,iii+1,iparm)=energia(k)
399 ! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
400 ! write (iout,*) "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
402 write (iout,'(2i5,f10.1,3e15.5)') i,iii,&
403 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
404 ! call enerprint(energia(0),fT)
407 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
408 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
409 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
410 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
411 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
412 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
413 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
414 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
415 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
416 write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
417 write (iout,'(f10.5,i10)') rmsdev,iscor
418 ! call enerprint(energia(0),fT)
419 call enerprint(energia(0))
420 call pdboutW(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
427 if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
428 write (ientout,rec=iii) &
429 ((csingle(l,k),l=1,3),k=1,nres),&
430 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
431 nss,(ihpb(k),jhpb(k),k=1,nss),&
432 potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
433 ! write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
435 if (separate_parset) then
436 snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
438 snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
440 ! write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
441 ! & " snk",snk_p(iR,ib,ipar)
443 snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
449 write (iout,*) "Me",me," scount",scount(me)
451 ! Master gathers updated numbers of conformations written by all procs.
452 call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
453 MPI_INTEGER, WHAM_COMM, IERROR)
457 indstart(i)=indend(i-1)+1
458 indend(i)=indstart(i)+scount_(i)-1
461 write (iout,*) "Revised conformation counts"
463 write (iout,'(a,i5,a,i7,a,i7,a,i7)') &
464 "Processor",i," indstart",indstart(i),&
465 " indend",indend(i)," count",scount_(i)
468 call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),&
469 MaxR*MaxT_h*nParmSet,&
470 MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
476 stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
480 write (iout,*) "Revised SNK"
483 write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') &
484 iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),&
485 (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
486 write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
489 write (iout,'("Total",i10)') stot(islice)
495 101 write (iout,*) "Error in scratchfile."
499 end subroutine enecalc
500 !------------------------------------------------------------------------------
501 logical function conf_check(ii,iprint)
503 use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,rad2deg,vbld
504 use geometry, only:int_from_cart1
505 ! include "DIMENSIONS"
506 ! include "DIMENSIONS.ZSCOPT"
507 ! include "DIMENSIONS.FREE"
511 ! include "COMMON.MPI"
513 ! include "COMMON.CHAIN"
514 ! include "COMMON.IOUNITS"
515 ! include "COMMON.PROTFILES"
516 ! include "COMMON.NAMES"
517 ! include "COMMON.VAR"
518 ! include "COMMON.SBRIDGE"
519 ! include "COMMON.GEO"
520 ! include "COMMON.FFIELD"
521 ! include "COMMON.ENEPS"
522 ! include "COMMON.LOCAL"
523 ! include "COMMON.WEIGHTS"
524 ! include "COMMON.INTERACT"
525 ! include "COMMON.FREE"
526 ! include "COMMON.ENERGIES"
527 ! include "COMMON.CONTROL"
528 ! include "COMMON.TORCNSTR"
531 integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
533 integer :: j,k,l,ii,itj,iprint,mnum,mnum1
534 if (.not.check_conf) then
538 call int_from_cart1(.false.)
543 if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
544 (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
546 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
547 " for conformation",ii
548 if (iprint.gt.1) then
549 write (iout,*) "The Cartesian geometry is:"
550 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
551 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
552 write (iout,*) "The internal geometry is:"
553 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
554 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
555 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
556 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
557 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
558 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
560 if (iprint.gt.0) write (iout,*) &
561 "This conformation WILL NOT be added to the database."
570 if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. &
571 (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
573 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
574 " for conformation",ii
575 if (iprint.gt.1) then
576 write (iout,*) "The Cartesian geometry is:"
577 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
578 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
579 write (iout,*) "The internal geometry is:"
580 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
581 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
582 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
583 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
584 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
585 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
587 if (iprint.gt.0) write (iout,*) &
588 "This conformation WILL NOT be added to the database."
594 if (theta(j).le.0.0d0) then
596 write (iout,*) "Zero theta angle(s) in conformation",ii
597 if (iprint.gt.1) then
598 write (iout,*) "The Cartesian geometry is:"
599 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
600 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
601 write (iout,*) "The internal geometry is:"
602 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
603 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
604 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
605 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
606 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
607 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
609 if (iprint.gt.0) write (iout,*) &
610 "This conformation WILL NOT be added to the database."
614 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
617 ! write (iout,*) "conf_check passed",ii
619 end function conf_check
620 !-----------------------------------------------------------------------------
622 !-----------------------------------------------------------------------------
623 subroutine store_parm(iparm)
625 ! Store parameters of set IPARM
626 ! valence angles and the side chains and energy parameters.
629 ! include 'DIMENSIONS'
630 ! include 'DIMENSIONS.ZSCOPT'
631 ! include 'DIMENSIONS.FREE'
632 ! include 'COMMON.IOUNITS'
633 ! include 'COMMON.CHAIN'
634 ! include 'COMMON.INTERACT'
635 ! include 'COMMON.GEO'
636 ! include 'COMMON.LOCAL'
637 ! include 'COMMON.TORSION'
638 ! include 'COMMON.FFIELD'
639 ! include 'COMMON.NAMES'
640 ! include 'COMMON.SBRIDGE'
641 ! include 'COMMON.SCROT'
642 ! include 'COMMON.SCCOR'
643 ! include 'COMMON.ALLPARM'
644 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
646 call alloc_enecalc_arrays(iparm)
647 !el allocate(ww_all(n_ene,iparm))
651 ww_all(3,iparm)=welec
652 ww_all(4,iparm)=wcorr
653 ww_all(5,iparm)=wcorr5
654 ww_all(6,iparm)=wcorr6
655 ww_all(7,iparm)=wel_loc
656 ww_all(8,iparm)=wturn3
657 ww_all(9,iparm)=wturn4
658 ww_all(10,iparm)=wturn6
659 ww_all(11,iparm)=wang
660 ww_all(12,iparm)=wscloc
661 ww_all(13,iparm)=wtor
662 ww_all(14,iparm)=wtor_d
663 ww_all(15,iparm)=wstrain
664 ww_all(16,iparm)=wvdwpp
665 ww_all(17,iparm)=wbond
666 ww_all(19,iparm)=wsccor
667 ! Store bond parameters
668 vbldp0_all(iparm)=vbldp0
671 nbondterm_all(i,iparm)=nbondterm(i)
673 vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
674 aksc_all(j,i,iparm)=aksc(j,i)
675 abond0_all(j,i,iparm)=abond0(j,i)
678 ! Store bond angle parameters
681 a0thet_all(i,iparm)=a0thet(i)
685 athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
686 bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
691 polthet_all(j,i,iparm)=polthet(j,i)
694 gthet_all(j,i,iparm)=gthet(j,i)
696 theta0_all(i,iparm)=theta0(i)
697 sig0_all(i,iparm)=sig0(i)
698 sigc0_all(i,iparm)=sigc0(i)
701 nthetyp_all(iparm)=nthetyp
702 ntheterm_all(iparm)=ntheterm
703 ntheterm2_all(iparm)=ntheterm2
704 ntheterm3_all(iparm)=ntheterm3
705 nsingle_all(iparm)=nsingle
706 ndouble_all(iparm)=ndouble
707 nntheterm_all(iparm)=nntheterm
709 ithetyp_all(i,iparm)=ithetyp(i)
712 do i=-nthetyp-1,nthetyp+1
713 do j=-nthetyp-1,nthetyp+1
714 do k=-nthetyp-1,nthetyp+1
715 aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
717 aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
721 bbthet_all(m,l,i,j,k,iblock,iparm)= &
722 bbthet(m,l,i,j,k,iblock)
723 ccthet_all(m,l,i,j,k,iblock,iparm)= &
724 ccthet(m,l,i,j,k,iblock)
725 ddthet_all(m,l,i,j,k,iblock,iparm)= &
726 ddthet(m,l,i,j,k,iblock)
727 eethet_all(m,l,i,j,k,iblock,iparm)= &
728 eethet(m,l,i,j,k,iblock)
734 if (iblock.eq.1) then
735 ffthet_all1(mm,m,l,i,j,k,iparm)=&
736 ffthet(mm,m,l,i,j,k,iblock)
737 ggthet_all1(mm,m,l,i,j,k,iparm)=&
738 ggthet(mm,m,l,i,j,k,iblock)
740 ffthet_all2(mm,m,l,i,j,k,iparm)=&
741 ffthet(mm,m,l,i,j,k,iblock)
742 ggthet_all2(mm,m,l,i,j,k,iparm)=&
743 ggthet(mm,m,l,i,j,k,iblock)
754 ! Store the sidechain rotamer parameters
757 !! write (iout,*) i,"storeparm1"
759 nlob_all(iii,iparm)=nlob(iii)
761 bsc_all(j,iii,iparm)=bsc(j,iii)
763 censc_all(k,j,i,iparm)=censc(k,j,i)
767 gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
775 sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
779 ! Store the torsional parameters
781 do i=-ntortyp+1,ntortyp-1
782 do j=-ntortyp+1,ntortyp-1
783 v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
784 nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
785 nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
786 do k=1,nterm(i,j,iblock)
787 v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
788 v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
790 do k=1,nlor(i,j,iblock)
791 vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
792 vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
793 vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
798 ! Store the double torsional parameters
800 do i=-ntortyp+1,ntortyp-1
801 do j=-ntortyp+1,ntortyp-1
802 do k=-ntortyp+1,ntortyp-1
803 ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
804 ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
805 do l=1,ntermd_1(i,j,k,iblock)
806 v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
807 v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
808 v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
809 v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
811 do l=1,ntermd_2(i,j,k,iblock)
812 do m=1,ntermd_2(i,j,k,iblock)
813 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
820 ! Store parameters of the cumulants
821 do i=-nloctyp,nloctyp
823 b1_all(j,i,iparm)=b1(j,i)
824 b1tilde_all(j,i,iparm)=b1tilde(j,i)
825 b2_all(j,i,iparm)=b2(j,i)
829 cc_all(k,j,i,iparm)=cc(k,j,i)
830 ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
831 dd_all(k,j,i,iparm)=dd(k,j,i)
832 dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
833 ee_all(k,j,i,iparm)=ee(k,j,i)
837 ! Store the parameters of electrostatic interactions
840 app_all(j,i,iparm)=app(j,i)
841 bpp_all(j,i,iparm)=bpp(j,i)
842 ael6_all(j,i,iparm)=ael6(j,i)
843 ael3_all(j,i,iparm)=ael3(j,i)
846 ! Store sidechain parameters
849 aa_all(j,i,iparm)=aa_aq(j,i)
850 bb_all(j,i,iparm)=bb_aq(j,i)
851 r0_all(j,i,iparm)=r0(j,i)
852 sigma_all(j,i,iparm)=sigma(j,i)
853 chi_all(j,i,iparm)=chi(j,i)
854 augm_all(j,i,iparm)=augm(j,i)
855 eps_all(j,i,iparm)=eps(j,i)
859 chip_all(i,iparm)=chip(i)
860 alp_all(i,iparm)=alp(i)
862 ! Store the SCp parameters
865 aad_all(i,j,iparm)=aad(i,j)
866 bad_all(i,j,iparm)=bad(i,j)
869 ! Store disulfide-bond parameters
878 ! Store SC-backbone correlation parameters
879 do i=-nsccortyp,nsccortyp
880 do j=-nsccortyp,nsccortyp
882 nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
886 do k=1,nterm_sccor(j,i)
887 v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
888 v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
893 write(iout,*)"end of store_parm"
895 end subroutine store_parm
896 !--------------------------------------------------------------------------
897 subroutine restore_parm(iparm)
899 ! Store parameters of set IPARM
900 ! valence angles and the side chains and energy parameters.
903 ! include 'DIMENSIONS'
904 ! include 'DIMENSIONS.ZSCOPT'
905 ! include 'DIMENSIONS.FREE'
906 ! include 'COMMON.IOUNITS'
907 ! include 'COMMON.CHAIN'
908 ! include 'COMMON.INTERACT'
909 ! include 'COMMON.GEO'
910 ! include 'COMMON.LOCAL'
911 ! include 'COMMON.TORSION'
912 ! include 'COMMON.FFIELD'
913 ! include 'COMMON.NAMES'
914 ! include 'COMMON.SBRIDGE'
915 ! include 'COMMON.SCROT'
916 ! include 'COMMON.SCCOR'
917 ! include 'COMMON.ALLPARM'
918 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
923 welec=ww_all(3,iparm)
924 wcorr=ww_all(4,iparm)
925 wcorr5=ww_all(5,iparm)
926 wcorr6=ww_all(6,iparm)
927 wel_loc=ww_all(7,iparm)
928 wturn3=ww_all(8,iparm)
929 wturn4=ww_all(9,iparm)
930 wturn6=ww_all(10,iparm)
931 wang=ww_all(11,iparm)
932 wscloc=ww_all(12,iparm)
933 wtor=ww_all(13,iparm)
934 wtor_d=ww_all(14,iparm)
935 wstrain=ww_all(15,iparm)
936 wvdwpp=ww_all(16,iparm)
937 wbond=ww_all(17,iparm)
938 wsccor=ww_all(19,iparm)
939 ! Restore bond parameters
940 vbldp0=vbldp0_all(iparm)
943 nbondterm(i)=nbondterm_all(i,iparm)
945 vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
946 aksc(j,i)=aksc_all(j,i,iparm)
947 abond0(j,i)=abond0_all(j,i,iparm)
950 ! Restore bond angle parameters
953 a0thet(i)=a0thet_all(i,iparm)
957 athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
958 bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
963 polthet(j,i)=polthet_all(j,i,iparm)
966 gthet(j,i)=gthet_all(j,i,iparm)
968 theta0(i)=theta0_all(i,iparm)
969 sig0(i)=sig0_all(i,iparm)
970 sigc0(i)=sigc0_all(i,iparm)
973 nthetyp=nthetyp_all(iparm)
974 ntheterm=ntheterm_all(iparm)
975 ntheterm2=ntheterm2_all(iparm)
976 ntheterm3=ntheterm3_all(iparm)
977 nsingle=nsingle_all(iparm)
978 ndouble=ndouble_all(iparm)
979 nntheterm=nntheterm_all(iparm)
981 ithetyp(i)=ithetyp_all(i,iparm)
984 do i=-nthetyp-1,nthetyp+1
985 do j=-nthetyp-1,nthetyp+1
986 do k=-nthetyp-1,nthetyp+1
987 aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
989 aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
993 bbthet(m,l,i,j,k,iblock)= &
994 bbthet_all(m,l,i,j,k,iblock,iparm)
995 ccthet(m,l,i,j,k,iblock)= &
996 ccthet_all(m,l,i,j,k,iblock,iparm)
997 ddthet(m,l,i,j,k,iblock)= &
998 ddthet_all(m,l,i,j,k,iblock,iparm)
999 eethet(m,l,i,j,k,iblock)= &
1000 eethet_all(m,l,i,j,k,iblock,iparm)
1006 if (iblock.eq.1) then
1007 ffthet(mm,m,l,i,j,k,iblock)= &
1008 ffthet_all1(mm,m,l,i,j,k,iparm)
1009 ggthet(mm,m,l,i,j,k,iblock)= &
1010 ggthet_all1(mm,m,l,i,j,k,iparm)
1012 ffthet(mm,m,l,i,j,k,iblock)= &
1013 ffthet_all2(mm,m,l,i,j,k,iparm)
1014 ggthet(mm,m,l,i,j,k,iblock)= &
1015 ggthet_all2(mm,m,l,i,j,k,iparm)
1025 ! Restore the sidechain rotamer parameters
1030 nlob(iii)=nlob_all(iii,iparm)
1032 bsc(j,iii)=bsc_all(j,iii,iparm)
1034 censc(k,j,i)=censc_all(k,j,i,iparm)
1038 gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1046 sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1050 ! Restore the torsional parameters
1052 do i=-ntortyp+1,ntortyp-1
1053 do j=-ntortyp+1,ntortyp-1
1054 v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1055 nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1056 nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1057 do k=1,nterm(i,j,iblock)
1058 v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1059 v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1061 do k=1,nlor(i,j,iblock)
1062 vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1063 vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1064 vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1069 ! Restore the double torsional parameters
1071 do i=-ntortyp+1,ntortyp-1
1072 do j=-ntortyp+1,ntortyp-1
1073 do k=-ntortyp+1,ntortyp-1
1074 ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1075 ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1076 do l=1,ntermd_1(i,j,k,iblock)
1077 v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1078 v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1079 v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1080 v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1082 do l=1,ntermd_2(i,j,k,iblock)
1083 do m=1,ntermd_2(i,j,k,iblock)
1084 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1091 ! Restore parameters of the cumulants
1092 do i=-nloctyp,nloctyp
1094 b1(j,i)=b1_all(j,i,iparm)
1095 b1tilde(j,i)=b1tilde_all(j,i,iparm)
1096 b2(j,i)=b2_all(j,i,iparm)
1100 cc(k,j,i)=cc_all(k,j,i,iparm)
1101 ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1102 dd(k,j,i)=dd_all(k,j,i,iparm)
1103 dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1104 ee(k,j,i)=ee_all(k,j,i,iparm)
1108 ! Restore the parameters of electrostatic interactions
1111 app(j,i)=app_all(j,i,iparm)
1112 bpp(j,i)=bpp_all(j,i,iparm)
1113 ael6(j,i)=ael6_all(j,i,iparm)
1114 ael3(j,i)=ael3_all(j,i,iparm)
1117 ! Restore sidechain parameters
1120 aa_aq(j,i)=aa_all(j,i,iparm)
1121 bb_aq(j,i)=bb_all(j,i,iparm)
1122 r0(j,i)=r0_all(j,i,iparm)
1123 sigma(j,i)=sigma_all(j,i,iparm)
1124 chi(j,i)=chi_all(j,i,iparm)
1125 augm(j,i)=augm_all(j,i,iparm)
1126 eps(j,i)=eps_all(j,i,iparm)
1130 chip(i)=chip_all(i,iparm)
1131 alp(i)=alp_all(i,iparm)
1133 ! Restore the SCp parameters
1136 aad(i,j)=aad_all(i,j,iparm)
1137 bad(i,j)=bad_all(i,j,iparm)
1140 ! Restore disulfide-bond parameters
1142 d0cm=d0cm_all(iparm)
1143 akcm=akcm_all(iparm)
1144 akth=akth_all(iparm)
1145 akct=akct_all(iparm)
1146 v1ss=v1ss_all(iparm)
1147 v2ss=v2ss_all(iparm)
1148 v3ss=v3ss_all(iparm)
1149 ! Restore SC-backbone correlation parameters
1150 do i=-nsccortyp,nsccortyp
1151 do j=-nsccortyp,nsccortyp
1153 nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1155 do k=1,nterm_sccor(j,i)
1156 v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1157 v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1163 end subroutine restore_parm
1164 !--------------------------------------------------------------------------
1166 !--------------------------------------------------------------------------
1167 subroutine make_ensembles(islice,*)
1168 ! construct the conformational ensembles at REMD temperatures
1169 use geometry_data, only:c
1170 use io_base, only:ilen
1171 use io_wham, only:pdboutW
1173 ! include "DIMENSIONS"
1174 ! include "DIMENSIONS.ZSCOPT"
1175 ! include "DIMENSIONS.FREE"
1178 ! include "COMMON.MPI"
1179 integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1181 ! include "COMMON.IOUNITS"
1182 ! include "COMMON.CONTROL"
1183 ! include "COMMON.FREE"
1184 ! include "COMMON.ENERGIES"
1185 ! include "COMMON.FFIELD"
1186 ! include "COMMON.INTERACT"
1187 ! include "COMMON.SBRIDGE"
1188 ! include "COMMON.CHAIN"
1189 ! include "COMMON.PROTFILES"
1190 ! include "COMMON.PROT"
1191 real(kind=4) :: csingle(3,nres*2)
1192 real(kind=8),dimension(6) :: fT,fTprim,fTbis
1193 real(kind=8) :: quot,quotl1,quotl,kfacl,&
1194 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1195 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1196 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1197 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
1198 ecation_prot, ecationcation
1199 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1200 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1201 character(len=80) :: bxname
1202 character(len=2) :: licz1,licz2
1203 character(len=3) :: licz3,licz4
1204 character(len=5) :: ctemper
1207 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1208 real(kind=8) :: enepot(MaxStr)
1209 integer :: iperm(MaxStr)
1211 integer,dimension(0:nprocs) :: scount_
1214 if (me.eq.Master) then
1216 write (licz2,'(bz,i2.2)') islice
1217 if (nslice.eq.1) then
1218 if (.not.separate_parset) then
1219 bxname = prefix(:ilen(prefix))//".bx"
1221 write (licz3,'(bz,i3.3)') myparm
1222 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1225 if (.not.separate_parset) then
1226 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1228 write (licz3,'(bz,i3.3)') myparm
1229 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1230 "_slice_"//licz2//".bx"
1233 open (ientout,file=bxname,status="unknown",&
1234 form="unformatted",access="direct",recl=lenrec1)
1239 if (iparm.ne.iparmprint) exit
1240 call restore_parm(iparm)
1243 write (iout,*) "iparm",iparm," ib",ib
1245 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1246 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1253 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1255 !el old rescale weights
1257 ! if (rescale_mode.eq.1) then
1258 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1260 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1261 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1262 #elif defined(FUNCT)
1273 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1275 ! else if (rescale_mode.eq.2) then
1276 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1277 !#if defined(FUNCTH)
1278 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1279 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1280 !#elif defined(FUNCT)
1288 ! fT(l)=1.12692801104297249644d0/ &
1289 ! dlog(dexp(quotl)+dexp(-quotl))
1291 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1292 ! else if (rescale_mode.eq.0) then
1298 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1302 ! el end old rescale weihgts
1303 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1310 evdw=enetb(1,i,iparm)
1311 ! evdw_t=enetb(21,i,iparm)
1312 evdw_t=enetb(20,i,iparm)
1314 ! evdw2_14=enetb(17,i,iparm)
1315 evdw2_14=enetb(18,i,iparm)
1316 evdw2=enetb(2,i,iparm)+evdw2_14
1318 evdw2=enetb(2,i,iparm)
1322 ees=enetb(3,i,iparm)
1323 evdw1=enetb(16,i,iparm)
1325 ees=enetb(3,i,iparm)
1328 ecorr=enetb(4,i,iparm)
1329 ecorr5=enetb(5,i,iparm)
1330 ecorr6=enetb(6,i,iparm)
1331 eel_loc=enetb(7,i,iparm)
1332 eello_turn3=enetb(8,i,iparm)
1333 eello_turn4=enetb(9,i,iparm)
1334 eello_turn6=enetb(10,i,iparm)
1335 ebe=enetb(11,i,iparm)
1336 escloc=enetb(12,i,iparm)
1337 etors=enetb(13,i,iparm)
1338 etors_d=enetb(14,i,iparm)
1339 ehpb=enetb(15,i,iparm)
1341 estr=enetb(17,i,iparm)
1342 ! estr=enetb(18,i,iparm)
1343 ! esccor=enetb(19,i,iparm)
1344 esccor=enetb(21,i,iparm)
1345 ! edihcnstr=enetb(20,i,iparm)
1346 edihcnstr=enetb(19,i,iparm)
1348 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1350 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1351 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1352 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1353 ! +ft(2)*wturn3*eello_turn3 &
1354 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1355 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1358 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1359 ! +ft(1)*welec*(ees+evdw1) &
1360 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1361 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1362 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1363 ! +ft(2)*wturn3*eello_turn3 &
1364 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1365 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1370 etot=wsc*evdw+wscp*evdw2+welec*ees &
1372 +wang*ebe+wtor*etors+wscloc*escloc &
1373 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1374 +wcorr6*ecorr6+wturn4*eello_turn4 &
1375 +wturn3*eello_turn3 &
1376 +wturn6*eello_turn6+wel_loc*eel_loc &
1377 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1378 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1380 etot=wsc*evdw+wscp*evdw2 &
1381 +welec*(ees+evdw1) &
1382 +wang*ebe+wtor*etors+wscloc*escloc &
1383 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1384 +wcorr6*ecorr6+wturn4*eello_turn4 &
1385 +wturn3*eello_turn3 &
1386 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1387 +wtor_d*etors_d+wsccor*esccor &
1388 +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
1393 beta_h(ib,iparm)*etot-entfac(i)
1396 write (iout,*) i,indstart(me)+i-1,ib,&
1397 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1398 -entfac(i),Fdimless(i)
1401 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1407 Fdimless_(i)=Fdimless(i)
1409 call MPI_Gatherv(Fdimless_(1),scount(me),&
1410 MPI_REAL,Fdimless(1),&
1411 scount_(0),idispl(0),MPI_REAL,Master,&
1414 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1415 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1416 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1418 call MPI_Gatherv(entfac(1),scount(me),&
1419 MPI_DOUBLE_PRECISION,entfac(1),&
1420 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1423 if (me.eq.Master) then
1425 write (iout,*) "The FDIMLESS array before sorting"
1427 write (iout,*) i,fdimless(i)
1434 call mysort1(ntot(islice),Fdimless,iperm)
1436 write (iout,*) "The FDIMLESS array after sorting"
1438 write (iout,*) i,iperm(i),fdimless(i)
1443 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1445 ! write (iout,*) "qfree",qfree
1448 do i=1,min0(ntot(islice),ensembles)
1449 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1451 write (iout,*) i,ib,beta_h(ib,iparm),&
1452 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1453 potE(iperm(i),iparm),&
1454 -entfac(iperm(i)),fdimless(i),sumprob
1456 if (sumprob.gt.0.99d0) goto 122
1462 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1464 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1469 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1472 if (iproc.ge.nprocs) then
1473 write (iout,*) "Fatal error: processor out of range",iproc
1478 close (ientout,status="delete")
1482 ik=ii-indstart(iproc)+1
1483 if (iproc.ne.Master) then
1484 if (me.eq.iproc) then
1486 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1487 " energy",potE(ik,iparm)
1489 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1490 Master,i,WHAM_COMM,IERROR)
1491 else if (me.eq.Master) then
1492 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1493 WHAM_COMM,STATUS,IERROR)
1495 else if (me.eq.Master) then
1496 enepot(i)=potE(ik,iparm)
1501 enepot(i)=potE(iperm(i),iparm)
1505 if (me.eq.Master) then
1507 write(licz3,'(bz,i3.3)') iparm
1508 write(licz2,'(bz,i2.2)') islice
1509 if (temper.lt.100.0d0) then
1510 write(ctemper,'(f3.0)') temper
1511 else if (temper.lt.1000.0) then
1512 write (ctemper,'(f4.0)') temper
1514 write (ctemper,'(f5.0)') temper
1516 if (nparmset.eq.1) then
1517 if (separate_parset) then
1518 write(licz4,'(bz,i3.3)') myparm
1519 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1521 pdbname=prefix(:ilen(prefix))
1524 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1526 if (nslice.eq.1) then
1527 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1528 ctemper(:ilen(ctemper))//"pdb"
1530 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1531 ctemper(:ilen(ctemper))//"pdb"
1533 open(ipdb,file=pdbname)
1535 read (ientout,rec=iperm(i)) &
1536 ((csingle(l,k),l=1,3),k=1,nres),&
1537 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1538 nss,(ihpb(k),jhpb(k),k=1,nss),&
1539 eini,efree,rmsdev,iscor
1546 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1556 close(ientout,status="delete")
1559 scount(i)=scount_(i)
1562 end subroutine make_ensembles
1563 !--------------------------------------------------------------------------
1564 subroutine mysort1(n, x, ipermut)
1566 integer :: i,j,imax,ipm,n
1567 real(kind=4) :: x(n)
1568 integer :: ipermut(n)
1569 real(kind=4) :: xtemp
1574 if (x(j).lt.xtemp) then
1582 ipermut(imax)=ipermut(i)
1586 end subroutine mysort1
1587 !--------------------------------------------------------------------------
1588 subroutine alloc_enecalc_arrays(iparm)
1591 use geometry_data, only:maxlob
1593 !---------------------------
1594 ! COMMON.ENERGIES form wham_data
1596 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1597 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1598 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1599 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1601 ! allocate ENECALC arrays
1602 !---------------------------
1605 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1606 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1607 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1608 aksc_all(maxbondterm,ntyp,iparm),&
1609 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1610 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1611 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1612 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1613 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1614 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1615 allocate(theta0_all(-ntyp:ntyp,iparm),&
1616 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1617 allocate(aa0thet_all(-nthetyp-1:nthetyp+1,&
1618 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1619 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1620 allocate(aathet_all(ntheterm,-nthetyp-1:nthetyp+1,&
1621 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1622 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1623 allocate(bbthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1624 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1625 allocate(ccthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1626 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1627 allocate(ddthet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1628 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1629 allocate(eethet_all(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1630 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2,iparm))
1631 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1632 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1633 allocate(ffthet_all1(ndouble,ndouble,ntheterm3,&
1634 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1635 -nthetyp-1:nthetyp+1,iparm))
1636 allocate(ggthet_all1(ndouble,ndouble,ntheterm3,&
1637 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1638 -nthetyp-1:nthetyp+1,iparm))
1639 allocate(ffthet_all2(ndouble,ndouble,ntheterm3,&
1640 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1641 -nthetyp-1:nthetyp+1,iparm))
1642 allocate(ggthet_all2(ndouble,ndouble,ntheterm3,&
1643 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,&
1644 -nthetyp-1:nthetyp+1,iparm))
1645 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1646 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1647 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1648 allocate(bsc_all(maxlob,ntyp,iparm))
1649 !(maxlob,ntyp,max_parm)
1650 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1651 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1652 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1653 allocate(v0_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1654 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1655 allocate(v1_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1656 allocate(v2_all(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1657 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1658 allocate(vlor1_all(maxlor,ntortyp,ntortyp,iparm))
1659 allocate(vlor2_all(maxlor,ntortyp,ntortyp,iparm))
1660 allocate(vlor3_all(maxlor,ntortyp,ntortyp,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1661 allocate(v1c_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1662 -ntortyp:ntortyp,2,iparm))
1663 allocate(v1s_all(2,maxtermd_1,-ntortyp:ntortyp,-ntortyp:ntortyp,&
1664 -ntortyp:ntortyp,2,iparm))
1665 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1666 allocate(v2c_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1667 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1668 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1669 allocate(v2s_all(maxtermd_2,maxtermd_2,-ntortyp:ntortyp,&
1670 -ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1671 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1672 allocate(b1_all(2,-nloctyp:nloctyp,iparm))
1673 allocate(b2_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1674 allocate(cc_all(2,2,-nloctyp:nloctyp,iparm))
1675 allocate(dd_all(2,2,-nloctyp:nloctyp,iparm))
1676 allocate(ee_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1677 allocate(ctilde_all(2,2,-nloctyp:nloctyp,iparm))
1678 allocate(dtilde_all(2,2,-nloctyp:nloctyp,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1679 allocate(b1tilde_all(2,-nloctyp:nloctyp,iparm)) !(2,-maxtor:maxtor,max_parm)
1680 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1681 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1682 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1683 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1684 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1685 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1686 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1687 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1688 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1689 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1690 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1691 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1692 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1693 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1694 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1695 allocate(nlor_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1696 allocate(nterm_all(-ntortyp:ntortyp,-ntortyp:ntortyp,2,iparm))
1697 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1698 allocate(ntermd1_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1699 -ntortyp:ntortyp,2,iparm))
1700 allocate(ntermd2_all(-ntortyp:ntortyp,-ntortyp:ntortyp,&
1701 -ntortyp:ntortyp,2,iparm))
1702 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1703 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1704 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1705 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1706 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1707 nsingle_all(iparm),&
1708 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1709 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1711 end subroutine alloc_enecalc_arrays
1712 !--------------------------------------------------------------------------
1713 !--------------------------------------------------------------------------