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,21)
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
534 if (.not.check_conf) then
538 call int_from_cart1(.false.)
540 if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. &
541 (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
543 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
544 " for conformation",ii
545 if (iprint.gt.1) then
546 write (iout,*) "The Cartesian geometry is:"
547 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
548 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
549 write (iout,*) "The internal geometry is:"
550 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
551 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
552 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
553 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
554 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
555 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
557 if (iprint.gt.0) write (iout,*) &
558 "This conformation WILL NOT be added to the database."
565 if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. &
566 (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
568 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
569 " for conformation",ii
570 if (iprint.gt.1) then
571 write (iout,*) "The Cartesian geometry is:"
572 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
573 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
574 write (iout,*) "The internal geometry is:"
575 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
576 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
577 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
578 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
579 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
580 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
582 if (iprint.gt.0) write (iout,*) &
583 "This conformation WILL NOT be added to the database."
589 if (theta(j).le.0.0d0) then
591 write (iout,*) "Zero theta angle(s) in conformation",ii
592 if (iprint.gt.1) then
593 write (iout,*) "The Cartesian geometry is:"
594 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
595 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
596 write (iout,*) "The internal geometry is:"
597 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
598 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
599 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
600 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
601 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
602 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
604 if (iprint.gt.0) write (iout,*) &
605 "This conformation WILL NOT be added to the database."
609 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
612 ! write (iout,*) "conf_check passed",ii
614 end function conf_check
615 !-----------------------------------------------------------------------------
617 !-----------------------------------------------------------------------------
618 subroutine store_parm(iparm)
620 ! Store parameters of set IPARM
621 ! valence angles and the side chains and energy parameters.
624 ! include 'DIMENSIONS'
625 ! include 'DIMENSIONS.ZSCOPT'
626 ! include 'DIMENSIONS.FREE'
627 ! include 'COMMON.IOUNITS'
628 ! include 'COMMON.CHAIN'
629 ! include 'COMMON.INTERACT'
630 ! include 'COMMON.GEO'
631 ! include 'COMMON.LOCAL'
632 ! include 'COMMON.TORSION'
633 ! include 'COMMON.FFIELD'
634 ! include 'COMMON.NAMES'
635 ! include 'COMMON.SBRIDGE'
636 ! include 'COMMON.SCROT'
637 ! include 'COMMON.SCCOR'
638 ! include 'COMMON.ALLPARM'
639 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
641 call alloc_enecalc_arrays(iparm)
642 !el allocate(ww_all(n_ene,iparm))
646 ww_all(3,iparm)=welec
647 ww_all(4,iparm)=wcorr
648 ww_all(5,iparm)=wcorr5
649 ww_all(6,iparm)=wcorr6
650 ww_all(7,iparm)=wel_loc
651 ww_all(8,iparm)=wturn3
652 ww_all(9,iparm)=wturn4
653 ww_all(10,iparm)=wturn6
654 ww_all(11,iparm)=wang
655 ww_all(12,iparm)=wscloc
656 ww_all(13,iparm)=wtor
657 ww_all(14,iparm)=wtor_d
658 ww_all(15,iparm)=wstrain
659 ww_all(16,iparm)=wvdwpp
660 ww_all(17,iparm)=wbond
661 ww_all(19,iparm)=wsccor
662 ! Store bond parameters
663 vbldp0_all(iparm)=vbldp0
666 nbondterm_all(i,iparm)=nbondterm(i)
668 vbldsc0_all(j,i,iparm)=vbldsc0(j,i)
669 aksc_all(j,i,iparm)=aksc(j,i)
670 abond0_all(j,i,iparm)=abond0(j,i)
673 ! Store bond angle parameters
676 a0thet_all(i,iparm)=a0thet(i)
680 athet_all(j,i,ichir1,ichir2,iparm)=athet(j,i,ichir1,ichir2)
681 bthet_all(j,i,ichir1,ichir2,iparm)=bthet(j,i,ichir1,ichir2)
686 polthet_all(j,i,iparm)=polthet(j,i)
689 gthet_all(j,i,iparm)=gthet(j,i)
691 theta0_all(i,iparm)=theta0(i)
692 sig0_all(i,iparm)=sig0(i)
693 sigc0_all(i,iparm)=sigc0(i)
696 nthetyp_all(iparm)=nthetyp
697 ntheterm_all(iparm)=ntheterm
698 ntheterm2_all(iparm)=ntheterm2
699 ntheterm3_all(iparm)=ntheterm3
700 nsingle_all(iparm)=nsingle
701 ndouble_all(iparm)=ndouble
702 nntheterm_all(iparm)=nntheterm
704 ithetyp_all(i,iparm)=ithetyp(i)
707 do i=-maxthetyp1,maxthetyp1
708 do j=-maxthetyp1,maxthetyp1
709 do k=-maxthetyp1,maxthetyp1
710 aa0thet_all(i,j,k,iblock,iparm)=aa0thet(i,j,k,iblock)
712 aathet_all(l,i,j,k,iblock,iparm)=aathet(l,i,j,k,iblock)
716 bbthet_all(m,l,i,j,k,iblock,iparm)= &
717 bbthet(m,l,i,j,k,iblock)
718 ccthet_all(m,l,i,j,k,iblock,iparm)= &
719 ccthet(m,l,i,j,k,iblock)
720 ddthet_all(m,l,i,j,k,iblock,iparm)= &
721 ddthet(m,l,i,j,k,iblock)
722 eethet_all(m,l,i,j,k,iblock,iparm)= &
723 eethet(m,l,i,j,k,iblock)
729 if (iblock.eq.1) then
730 ffthet_all1(mm,m,l,i,j,k,iparm)=&
731 ffthet(mm,m,l,i,j,k,iblock)
732 ggthet_all1(mm,m,l,i,j,k,iparm)=&
733 ggthet(mm,m,l,i,j,k,iblock)
735 ffthet_all2(mm,m,l,i,j,k,iparm)=&
736 ffthet(mm,m,l,i,j,k,iblock)
737 ggthet_all2(mm,m,l,i,j,k,iparm)=&
738 ggthet(mm,m,l,i,j,k,iblock)
749 ! Store the sidechain rotamer parameters
752 !! write (iout,*) i,"storeparm1"
754 nlob_all(iii,iparm)=nlob(iii)
756 bsc_all(j,iii,iparm)=bsc(j,iii)
758 censc_all(k,j,i,iparm)=censc(k,j,i)
762 gaussc_all(l,k,j,i,iparm)=gaussc(l,k,j,i)
770 sc_parmin_all(j,i,iparm)=sc_parmin(j,i)
774 ! Store the torsional parameters
776 do i=-ntortyp+1,ntortyp-1
777 do j=-ntortyp+1,ntortyp-1
778 v0_all(i,j,iblock,iparm)=v0(i,j,iblock)
779 nterm_all(i,j,iblock,iparm)=nterm(i,j,iblock)
780 nlor_all(i,j,iblock,iparm)=nlor(i,j,iblock)
781 do k=1,nterm(i,j,iblock)
782 v1_all(k,i,j,iblock,iparm)=v1(k,i,j,iblock)
783 v2_all(k,i,j,iblock,iparm)=v2(k,i,j,iblock)
785 do k=1,nlor(i,j,iblock)
786 vlor1_all(k,i,j,iparm)=vlor1(k,i,j)
787 vlor2_all(k,i,j,iparm)=vlor2(k,i,j)
788 vlor3_all(k,i,j,iparm)=vlor3(k,i,j)
793 ! Store the double torsional parameters
795 do i=-ntortyp+1,ntortyp-1
796 do j=-ntortyp+1,ntortyp-1
797 do k=-ntortyp+1,ntortyp-1
798 ntermd1_all(i,j,k,iblock,iparm)=ntermd_1(i,j,k,iblock)
799 ntermd2_all(i,j,k,iblock,iparm)=ntermd_2(i,j,k,iblock)
800 do l=1,ntermd_1(i,j,k,iblock)
801 v1c_all(1,l,i,j,k,iblock,iparm)=v1c(1,l,i,j,k,iblock)
802 v1c_all(2,l,i,j,k,iblock,iparm)=v1c(2,l,i,j,k,iblock)
803 v2c_all(1,l,i,j,k,iblock,iparm)=v2c(1,l,i,j,k,iblock)
804 v2c_all(2,l,i,j,k,iblock,iparm)=v2c(2,l,i,j,k,iblock)
806 do l=1,ntermd_2(i,j,k,iblock)
807 do m=1,ntermd_2(i,j,k,iblock)
808 v2s_all(l,m,i,j,k,iblock,iparm)=v2s(l,m,i,j,k,iblock)
815 ! Store parameters of the cumulants
816 do i=-nloctyp,nloctyp
818 b1_all(j,i,iparm)=b1(j,i)
819 b1tilde_all(j,i,iparm)=b1tilde(j,i)
820 b2_all(j,i,iparm)=b2(j,i)
824 cc_all(k,j,i,iparm)=cc(k,j,i)
825 ctilde_all(k,j,i,iparm)=ctilde(k,j,i)
826 dd_all(k,j,i,iparm)=dd(k,j,i)
827 dtilde_all(k,j,i,iparm)=dtilde(k,j,i)
828 ee_all(k,j,i,iparm)=ee(k,j,i)
832 ! Store the parameters of electrostatic interactions
835 app_all(j,i,iparm)=app(j,i)
836 bpp_all(j,i,iparm)=bpp(j,i)
837 ael6_all(j,i,iparm)=ael6(j,i)
838 ael3_all(j,i,iparm)=ael3(j,i)
841 ! Store sidechain parameters
844 aa_all(j,i,iparm)=aa(j,i)
845 bb_all(j,i,iparm)=bb(j,i)
846 r0_all(j,i,iparm)=r0(j,i)
847 sigma_all(j,i,iparm)=sigma(j,i)
848 chi_all(j,i,iparm)=chi(j,i)
849 augm_all(j,i,iparm)=augm(j,i)
850 eps_all(j,i,iparm)=eps(j,i)
854 chip_all(i,iparm)=chip(i)
855 alp_all(i,iparm)=alp(i)
857 ! Store the SCp parameters
860 aad_all(i,j,iparm)=aad(i,j)
861 bad_all(i,j,iparm)=bad(i,j)
864 ! Store disulfide-bond parameters
873 ! Store SC-backbone correlation parameters
874 do i=-nsccortyp,nsccortyp
875 do j=-nsccortyp,nsccortyp
877 nterm_sccor_all(j,i,iparm)=nterm_sccor(j,i)
881 do k=1,nterm_sccor(j,i)
882 v1sccor_all(k,l,j,i,iparm)=v1sccor(k,l,j,i)
883 v2sccor_all(k,l,j,i,iparm)=v2sccor(k,l,j,i)
888 write(iout,*)"end of store_parm"
890 end subroutine store_parm
891 !--------------------------------------------------------------------------
892 subroutine restore_parm(iparm)
894 ! Store parameters of set IPARM
895 ! valence angles and the side chains and energy parameters.
898 ! include 'DIMENSIONS'
899 ! include 'DIMENSIONS.ZSCOPT'
900 ! include 'DIMENSIONS.FREE'
901 ! include 'COMMON.IOUNITS'
902 ! include 'COMMON.CHAIN'
903 ! include 'COMMON.INTERACT'
904 ! include 'COMMON.GEO'
905 ! include 'COMMON.LOCAL'
906 ! include 'COMMON.TORSION'
907 ! include 'COMMON.FFIELD'
908 ! include 'COMMON.NAMES'
909 ! include 'COMMON.SBRIDGE'
910 ! include 'COMMON.SCROT'
911 ! include 'COMMON.SCCOR'
912 ! include 'COMMON.ALLPARM'
913 integer :: i,j,k,l,m,mm,iparm,ichir1,ichir2,iblock,iii
918 welec=ww_all(3,iparm)
919 wcorr=ww_all(4,iparm)
920 wcorr5=ww_all(5,iparm)
921 wcorr6=ww_all(6,iparm)
922 wel_loc=ww_all(7,iparm)
923 wturn3=ww_all(8,iparm)
924 wturn4=ww_all(9,iparm)
925 wturn6=ww_all(10,iparm)
926 wang=ww_all(11,iparm)
927 wscloc=ww_all(12,iparm)
928 wtor=ww_all(13,iparm)
929 wtor_d=ww_all(14,iparm)
930 wstrain=ww_all(15,iparm)
931 wvdwpp=ww_all(16,iparm)
932 wbond=ww_all(17,iparm)
933 wsccor=ww_all(19,iparm)
934 ! Restore bond parameters
935 vbldp0=vbldp0_all(iparm)
938 nbondterm(i)=nbondterm_all(i,iparm)
940 vbldsc0(j,i)=vbldsc0_all(j,i,iparm)
941 aksc(j,i)=aksc_all(j,i,iparm)
942 abond0(j,i)=abond0_all(j,i,iparm)
945 ! Restore bond angle parameters
948 a0thet(i)=a0thet_all(i,iparm)
952 athet(j,i,ichir1,ichir2)=athet_all(j,i,ichir1,ichir2,iparm)
953 bthet(j,i,ichir1,ichir2)=bthet_all(j,i,ichir1,ichir2,iparm)
958 polthet(j,i)=polthet_all(j,i,iparm)
961 gthet(j,i)=gthet_all(j,i,iparm)
963 theta0(i)=theta0_all(i,iparm)
964 sig0(i)=sig0_all(i,iparm)
965 sigc0(i)=sigc0_all(i,iparm)
968 nthetyp=nthetyp_all(iparm)
969 ntheterm=ntheterm_all(iparm)
970 ntheterm2=ntheterm2_all(iparm)
971 ntheterm3=ntheterm3_all(iparm)
972 nsingle=nsingle_all(iparm)
973 ndouble=ndouble_all(iparm)
974 nntheterm=nntheterm_all(iparm)
976 ithetyp(i)=ithetyp_all(i,iparm)
979 do i=-maxthetyp1,maxthetyp1
980 do j=-maxthetyp1,maxthetyp1
981 do k=-maxthetyp1,maxthetyp1
982 aa0thet(i,j,k,iblock)=aa0thet_all(i,j,k,iblock,iparm)
984 aathet(l,i,j,k,iblock)=aathet_all(l,i,j,k,iblock,iparm)
988 bbthet(m,l,i,j,k,iblock)= &
989 bbthet_all(m,l,i,j,k,iblock,iparm)
990 ccthet(m,l,i,j,k,iblock)= &
991 ccthet_all(m,l,i,j,k,iblock,iparm)
992 ddthet(m,l,i,j,k,iblock)= &
993 ddthet_all(m,l,i,j,k,iblock,iparm)
994 eethet(m,l,i,j,k,iblock)= &
995 eethet_all(m,l,i,j,k,iblock,iparm)
1001 if (iblock.eq.1) then
1002 ffthet(mm,m,l,i,j,k,iblock)= &
1003 ffthet_all1(mm,m,l,i,j,k,iparm)
1004 ggthet(mm,m,l,i,j,k,iblock)= &
1005 ggthet_all1(mm,m,l,i,j,k,iparm)
1007 ffthet(mm,m,l,i,j,k,iblock)= &
1008 ffthet_all2(mm,m,l,i,j,k,iparm)
1009 ggthet(mm,m,l,i,j,k,iblock)= &
1010 ggthet_all2(mm,m,l,i,j,k,iparm)
1020 ! Restore the sidechain rotamer parameters
1025 nlob(iii)=nlob_all(iii,iparm)
1027 bsc(j,iii)=bsc_all(j,iii,iparm)
1029 censc(k,j,i)=censc_all(k,j,i,iparm)
1033 gaussc(l,k,j,i)=gaussc_all(l,k,j,i,iparm)
1041 sc_parmin(j,i)=sc_parmin_all(j,i,iparm)
1045 ! Restore the torsional parameters
1047 do i=-ntortyp+1,ntortyp-1
1048 do j=-ntortyp+1,ntortyp-1
1049 v0(i,j,iblock)=v0_all(i,j,iblock,iparm)
1050 nterm(i,j,iblock)=nterm_all(i,j,iblock,iparm)
1051 nlor(i,j,iblock)=nlor_all(i,j,iblock,iparm)
1052 do k=1,nterm(i,j,iblock)
1053 v1(k,i,j,iblock)=v1_all(k,i,j,iblock,iparm)
1054 v2(k,i,j,iblock)=v2_all(k,i,j,iblock,iparm)
1056 do k=1,nlor(i,j,iblock)
1057 vlor1(k,i,j)=vlor1_all(k,i,j,iparm)
1058 vlor2(k,i,j)=vlor2_all(k,i,j,iparm)
1059 vlor3(k,i,j)=vlor3_all(k,i,j,iparm)
1064 ! Restore the double torsional parameters
1066 do i=-ntortyp+1,ntortyp-1
1067 do j=-ntortyp+1,ntortyp-1
1068 do k=-ntortyp+1,ntortyp-1
1069 ntermd_1(i,j,k,iblock)=ntermd1_all(i,j,k,iblock,iparm)
1070 ntermd_2(i,j,k,iblock)=ntermd2_all(i,j,k,iblock,iparm)
1071 do l=1,ntermd_1(i,j,k,iblock)
1072 v1c(1,l,i,j,k,iblock)=v1c_all(1,l,i,j,k,iblock,iparm)
1073 v1c(2,l,i,j,k,iblock)=v1c_all(2,l,i,j,k,iblock,iparm)
1074 v2c(1,l,i,j,k,iblock)=v2c_all(1,l,i,j,k,iblock,iparm)
1075 v2c(2,l,i,j,k,iblock)=v2c_all(2,l,i,j,k,iblock,iparm)
1077 do l=1,ntermd_2(i,j,k,iblock)
1078 do m=1,ntermd_2(i,j,k,iblock)
1079 v2s(l,m,i,j,k,iblock)=v2s_all(l,m,i,j,k,iblock,iparm)
1086 ! Restore parameters of the cumulants
1087 do i=-nloctyp,nloctyp
1089 b1(j,i)=b1_all(j,i,iparm)
1090 b1tilde(j,i)=b1tilde_all(j,i,iparm)
1091 b2(j,i)=b2_all(j,i,iparm)
1095 cc(k,j,i)=cc_all(k,j,i,iparm)
1096 ctilde(k,j,i)=ctilde_all(k,j,i,iparm)
1097 dd(k,j,i)=dd_all(k,j,i,iparm)
1098 dtilde(k,j,i)=dtilde_all(k,j,i,iparm)
1099 ee(k,j,i)=ee_all(k,j,i,iparm)
1103 ! Restore the parameters of electrostatic interactions
1106 app(j,i)=app_all(j,i,iparm)
1107 bpp(j,i)=bpp_all(j,i,iparm)
1108 ael6(j,i)=ael6_all(j,i,iparm)
1109 ael3(j,i)=ael3_all(j,i,iparm)
1112 ! Restore sidechain parameters
1115 aa(j,i)=aa_all(j,i,iparm)
1116 bb(j,i)=bb_all(j,i,iparm)
1117 r0(j,i)=r0_all(j,i,iparm)
1118 sigma(j,i)=sigma_all(j,i,iparm)
1119 chi(j,i)=chi_all(j,i,iparm)
1120 augm(j,i)=augm_all(j,i,iparm)
1121 eps(j,i)=eps_all(j,i,iparm)
1125 chip(i)=chip_all(i,iparm)
1126 alp(i)=alp_all(i,iparm)
1128 ! Restore the SCp parameters
1131 aad(i,j)=aad_all(i,j,iparm)
1132 bad(i,j)=bad_all(i,j,iparm)
1135 ! Restore disulfide-bond parameters
1137 d0cm=d0cm_all(iparm)
1138 akcm=akcm_all(iparm)
1139 akth=akth_all(iparm)
1140 akct=akct_all(iparm)
1141 v1ss=v1ss_all(iparm)
1142 v2ss=v2ss_all(iparm)
1143 v3ss=v3ss_all(iparm)
1144 ! Restore SC-backbone correlation parameters
1145 do i=-nsccortyp,nsccortyp
1146 do j=-nsccortyp,nsccortyp
1148 nterm_sccor(j,i)=nterm_sccor_all(j,i,iparm)
1150 do k=1,nterm_sccor(j,i)
1151 v1sccor(k,l,j,i)=v1sccor_all(k,l,j,i,iparm)
1152 v2sccor(k,l,j,i)=v2sccor_all(k,l,j,i,iparm)
1158 end subroutine restore_parm
1159 !--------------------------------------------------------------------------
1161 !--------------------------------------------------------------------------
1162 subroutine make_ensembles(islice,*)
1163 ! construct the conformational ensembles at REMD temperatures
1164 use geometry_data, only:c
1165 use io_base, only:ilen
1166 use io_wham, only:pdboutW
1168 ! include "DIMENSIONS"
1169 ! include "DIMENSIONS.ZSCOPT"
1170 ! include "DIMENSIONS.FREE"
1173 ! include "COMMON.MPI"
1174 integer :: ierror,errcode,status(MPI_STATUS_SIZE)
1176 ! include "COMMON.IOUNITS"
1177 ! include "COMMON.CONTROL"
1178 ! include "COMMON.FREE"
1179 ! include "COMMON.ENERGIES"
1180 ! include "COMMON.FFIELD"
1181 ! include "COMMON.INTERACT"
1182 ! include "COMMON.SBRIDGE"
1183 ! include "COMMON.CHAIN"
1184 ! include "COMMON.PROTFILES"
1185 ! include "COMMON.PROT"
1186 real(kind=4) :: csingle(3,nres*2)
1187 real(kind=8),dimension(6) :: fT,fTprim,fTbis
1188 real(kind=8) :: quot,quotl1,quotl,kfacl,&
1189 eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
1190 real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
1191 escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
1192 eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
1193 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
1194 real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
1195 character(len=80) :: bxname
1196 character(len=2) :: licz1,licz2
1197 character(len=3) :: licz3,licz4
1198 character(len=5) :: ctemper
1201 real(kind=4) :: Fdimless(MaxStr),Fdimless_(MaxStr)
1202 real(kind=8) :: enepot(MaxStr)
1203 integer :: iperm(MaxStr)
1205 integer,dimension(0:nprocs) :: scount_
1208 if (me.eq.Master) then
1210 write (licz2,'(bz,i2.2)') islice
1211 if (nslice.eq.1) then
1212 if (.not.separate_parset) then
1213 bxname = prefix(:ilen(prefix))//".bx"
1215 write (licz3,'(bz,i3.3)') myparm
1216 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
1219 if (.not.separate_parset) then
1220 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1222 write (licz3,'(bz,i3.3)') myparm
1223 bxname = prefix(:ilen(prefix))//"par_"//licz3// &
1224 "_slice_"//licz2//".bx"
1227 open (ientout,file=bxname,status="unknown",&
1228 form="unformatted",access="direct",recl=lenrec1)
1233 if (iparm.ne.iparmprint) exit
1234 call restore_parm(iparm)
1237 write (iout,*) "iparm",iparm," ib",ib
1239 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1240 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1247 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1249 !el old rescale weights
1251 ! if (rescale_mode.eq.1) then
1252 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1254 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1255 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
1256 #elif defined(FUNCT)
1267 ! fT(l)=kfacl/(kfacl-1.0d0+quotl)
1269 ! else if (rescale_mode.eq.2) then
1270 ! quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
1271 !#if defined(FUNCTH)
1272 ! tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
1273 ! ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
1274 !#elif defined(FUNCT)
1282 ! fT(l)=1.12692801104297249644d0/ &
1283 ! dlog(dexp(quotl)+dexp(-quotl))
1285 ! write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
1286 ! else if (rescale_mode.eq.0) then
1292 ! "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
1296 ! el end old rescale weihgts
1297 call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
1304 evdw=enetb(1,i,iparm)
1305 ! evdw_t=enetb(21,i,iparm)
1306 evdw_t=enetb(20,i,iparm)
1308 ! evdw2_14=enetb(17,i,iparm)
1309 evdw2_14=enetb(18,i,iparm)
1310 evdw2=enetb(2,i,iparm)+evdw2_14
1312 evdw2=enetb(2,i,iparm)
1316 ees=enetb(3,i,iparm)
1317 evdw1=enetb(16,i,iparm)
1319 ees=enetb(3,i,iparm)
1322 ecorr=enetb(4,i,iparm)
1323 ecorr5=enetb(5,i,iparm)
1324 ecorr6=enetb(6,i,iparm)
1325 eel_loc=enetb(7,i,iparm)
1326 eello_turn3=enetb(8,i,iparm)
1327 eello_turn4=enetb(9,i,iparm)
1328 eello_turn6=enetb(10,i,iparm)
1329 ebe=enetb(11,i,iparm)
1330 escloc=enetb(12,i,iparm)
1331 etors=enetb(13,i,iparm)
1332 etors_d=enetb(14,i,iparm)
1333 ehpb=enetb(15,i,iparm)
1335 estr=enetb(17,i,iparm)
1336 ! estr=enetb(18,i,iparm)
1337 ! esccor=enetb(19,i,iparm)
1338 esccor=enetb(21,i,iparm)
1339 ! edihcnstr=enetb(20,i,iparm)
1340 edihcnstr=enetb(19,i,iparm)
1342 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1344 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1345 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1346 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1347 ! +ft(2)*wturn3*eello_turn3 &
1348 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
1349 ! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1352 ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1353 ! +ft(1)*welec*(ees+evdw1) &
1354 ! +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1355 ! +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1356 ! +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1357 ! +ft(2)*wturn3*eello_turn3 &
1358 ! +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1359 ! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1364 etot=wsc*evdw+wscp*evdw2+welec*ees &
1366 +wang*ebe+wtor*etors+wscloc*escloc &
1367 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1368 +wcorr6*ecorr6+wturn4*eello_turn4 &
1369 +wturn3*eello_turn3 &
1370 +wturn6*eello_turn6+wel_loc*eel_loc &
1371 +edihcnstr+wtor_d*etors_d+wsccor*esccor &
1374 etot=wsc*evdw+wscp*evdw2 &
1375 +welec*(ees+evdw1) &
1376 +wang*ebe+wtor*etors+wscloc*escloc &
1377 +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
1378 +wcorr6*ecorr6+wturn4*eello_turn4 &
1379 +wturn3*eello_turn3 &
1380 +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
1381 +wtor_d*etors_d+wsccor*esccor &
1387 beta_h(ib,iparm)*etot-entfac(i)
1390 write (iout,*) i,indstart(me)+i-1,ib,&
1391 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
1392 -entfac(i),Fdimless(i)
1395 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
1401 Fdimless_(i)=Fdimless(i)
1403 call MPI_Gatherv(Fdimless_(1),scount(me),&
1404 MPI_REAL,Fdimless(1),&
1405 scount_(0),idispl(0),MPI_REAL,Master,&
1408 call MPI_Gatherv(potE(1,iparm),scount_(me),&
1409 MPI_DOUBLE_PRECISION,potE(1,iparm),&
1410 scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1412 call MPI_Gatherv(entfac(1),scount(me),&
1413 MPI_DOUBLE_PRECISION,entfac(1),&
1414 scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
1417 if (me.eq.Master) then
1419 write (iout,*) "The FDIMLESS array before sorting"
1421 write (iout,*) i,fdimless(i)
1428 call mysort1(ntot(islice),Fdimless,iperm)
1430 write (iout,*) "The FDIMLESS array after sorting"
1432 write (iout,*) i,iperm(i),fdimless(i)
1437 qfree=qfree+exp(-fdimless(i)+fdimless(1))
1439 ! write (iout,*) "qfree",qfree
1442 do i=1,min0(ntot(islice),ensembles)
1443 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
1445 write (iout,*) i,ib,beta_h(ib,iparm),&
1446 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
1447 potE(iperm(i),iparm),&
1448 -entfac(iperm(i)),fdimless(i),sumprob
1450 if (sumprob.gt.0.99d0) goto 122
1456 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,&
1458 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,&
1463 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
1466 if (iproc.ge.nprocs) then
1467 write (iout,*) "Fatal error: processor out of range",iproc
1472 close (ientout,status="delete")
1476 ik=ii-indstart(iproc)+1
1477 if (iproc.ne.Master) then
1478 if (me.eq.iproc) then
1480 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
1481 " energy",potE(ik,iparm)
1483 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
1484 Master,i,WHAM_COMM,IERROR)
1485 else if (me.eq.Master) then
1486 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,&
1487 WHAM_COMM,STATUS,IERROR)
1489 else if (me.eq.Master) then
1490 enepot(i)=potE(ik,iparm)
1495 enepot(i)=potE(iperm(i),iparm)
1499 if (me.eq.Master) then
1501 write(licz3,'(bz,i3.3)') iparm
1502 write(licz2,'(bz,i2.2)') islice
1503 if (temper.lt.100.0d0) then
1504 write(ctemper,'(f3.0)') temper
1505 else if (temper.lt.1000.0) then
1506 write (ctemper,'(f4.0)') temper
1508 write (ctemper,'(f5.0)') temper
1510 if (nparmset.eq.1) then
1511 if (separate_parset) then
1512 write(licz4,'(bz,i3.3)') myparm
1513 pdbname=prefix(:ilen(prefix))//"_par"//licz4
1515 pdbname=prefix(:ilen(prefix))
1518 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
1520 if (nslice.eq.1) then
1521 pdbname=pdbname(:ilen(pdbname))//"_T_"// &
1522 ctemper(:ilen(ctemper))//"pdb"
1524 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// &
1525 ctemper(:ilen(ctemper))//"pdb"
1527 open(ipdb,file=pdbname)
1529 read (ientout,rec=iperm(i)) &
1530 ((csingle(l,k),l=1,3),k=1,nres),&
1531 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1532 nss,(ihpb(k),jhpb(k),k=1,nss),&
1533 eini,efree,rmsdev,iscor
1540 call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev)
1550 close(ientout,status="delete")
1553 scount(i)=scount_(i)
1556 end subroutine make_ensembles
1557 !--------------------------------------------------------------------------
1558 subroutine mysort1(n, x, ipermut)
1560 integer :: i,j,imax,ipm,n
1561 real(kind=4) :: x(n)
1562 integer :: ipermut(n)
1563 real(kind=4) :: xtemp
1568 if (x(j).lt.xtemp) then
1576 ipermut(imax)=ipermut(i)
1580 end subroutine mysort1
1581 !--------------------------------------------------------------------------
1582 subroutine alloc_enecalc_arrays(iparm)
1585 use geometry_data, only:maxlob
1587 !---------------------------
1588 ! COMMON.ENERGIES form wham_data
1590 allocate(potE(MaxStr_Proc,iparm)) !(MaxStr_Proc,Max_Parm)
1591 allocate(entfac(MaxStr_Proc)) !(MaxStr_Proc)
1592 allocate(q(nQ+2,MaxStr_Proc)) !(MaxQ+2,MaxStr_Proc)
1593 allocate(enetb(max_ene,MaxStr_Proc,iparm)) !(max_ene,MaxStr_Proc,Max_Parm)
1595 ! allocate ENECALC arrays
1596 !---------------------------
1599 allocate(ww_all(max_eneW,iparm)) !(max_ene,max_parm) ! max_eneW
1600 allocate(vbldp0_all(iparm),akp_all(nParmSet)) !(max_parm)
1601 allocate(vbldsc0_all(maxbondterm,ntyp,iparm),&
1602 aksc_all(maxbondterm,ntyp,iparm),&
1603 abond0_all(maxbondterm,ntyp,iparm)) !(maxbondterm,ntyp,max_parm)
1604 allocate(a0thet_all(-ntyp:ntyp,iparm)) !(-ntyp:ntyp,max_parm)
1605 allocate(athet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm),&
1606 bthet_all(2,-ntyp:ntyp,-1:1,-1:1,iparm)) !(2,-ntyp:ntyp,-1:1,-1:1,max_parm)
1607 allocate(polthet_all(0:3,-ntyp:ntyp,iparm)) !(0:3,-ntyp:ntyp,max_parm)
1608 allocate(gthet_all(3,-ntyp:ntyp,iparm)) !(3,-ntyp:ntyp,max_parm)
1609 allocate(theta0_all(-ntyp:ntyp,iparm),&
1610 sig0_all(-ntyp:ntyp,iparm),sigc0_all(-ntyp:ntyp,nParmSet)) !(-ntyp:ntyp,max_parm)
1611 allocate(aa0thet_all(-maxthetyp1:maxthetyp1,&
1612 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1613 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1614 allocate(aathet_all(maxtheterm,-maxthetyp1:maxthetyp1,&
1615 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1616 !(maxtheterm,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1617 allocate(bbthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1618 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1619 allocate(ccthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1620 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1621 allocate(ddthet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1622 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1623 allocate(eethet_all(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1624 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,iparm))
1625 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,
1626 ! & -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2,max_parm)
1627 allocate(ffthet_all1(maxdouble,maxdouble,maxtheterm3,&
1628 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
1629 -maxthetyp1:maxthetyp1,iparm))
1630 allocate(ggthet_all1(maxdouble,maxdouble,maxtheterm3,&
1631 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
1632 -maxthetyp1:maxthetyp1,iparm))
1633 allocate(ffthet_all2(maxdouble,maxdouble,maxtheterm3,&
1634 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
1635 -maxthetyp1:maxthetyp1,iparm))
1636 allocate(ggthet_all2(maxdouble,maxdouble,maxtheterm3,&
1637 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,&
1638 -maxthetyp1:maxthetyp1,iparm))
1639 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1640 !-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,max_parm)
1641 allocate(dsc_all(ntyp1,iparm),dsc0_all(ntyp1,nParmSet)) !(ntyp1,max_parm)
1642 allocate(bsc_all(maxlob,ntyp,iparm))
1643 !(maxlob,ntyp,max_parm)
1644 allocate(censc_all(3,maxlob,-ntyp:ntyp,iparm)) !(3,maxlob,-ntyp:ntyp,max_parm)
1645 allocate(gaussc_all(3,3,maxlob,-ntyp:ntyp,iparm)) !(3,3,maxlob,-ntyp:ntyp,max_parm)
1646 allocate(sc_parmin_all(65,ntyp,iparm)) !(65,ntyp,max_parm)
1647 allocate(v0_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
1648 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1649 allocate(v1_all(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,iparm))
1650 allocate(v2_all(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,iparm))
1651 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1652 allocate(vlor1_all(maxlor,maxtor,maxtor,iparm))
1653 allocate(vlor2_all(maxlor,maxtor,maxtor,iparm))
1654 allocate(vlor3_all(maxlor,maxtor,maxtor,iparm)) !(maxlor,maxtor,maxtor,max_parm)
1655 allocate(v1c_all(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,&
1656 -maxtor:maxtor,2,iparm))
1657 allocate(v1s_all(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,&
1658 -maxtor:maxtor,2,iparm))
1659 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1660 allocate(v2c_all(maxtermd_2,maxtermd_2,-maxtor:maxtor,&
1661 -maxtor:maxtor,-maxtor:maxtor,2,iparm))
1662 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1663 allocate(v2s_all(maxtermd_2,maxtermd_2,-maxtor:maxtor,&
1664 -maxtor:maxtor,-maxtor:maxtor,2,iparm))
1665 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1666 allocate(b1_all(2,-maxtor:maxtor,iparm))
1667 allocate(b2_all(2,-maxtor:maxtor,iparm)) !(2,-maxtor:maxtor,max_parm)
1668 allocate(cc_all(2,2,-maxtor:maxtor,iparm))
1669 allocate(dd_all(2,2,-maxtor:maxtor,iparm))
1670 allocate(ee_all(2,2,-maxtor:maxtor,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1671 allocate(ctilde_all(2,2,-maxtor:maxtor,iparm))
1672 allocate(dtilde_all(2,2,-maxtor:maxtor,iparm)) !(2,2,-maxtor:maxtor,max_parm)
1673 allocate(b1tilde_all(2,-maxtor:maxtor,iparm)) !(2,-maxtor:maxtor,max_parm)
1674 allocate(app_all(2,2,iparm),bpp_all(2,2,nParmSet),&
1675 ael6_all(2,2,iparm),ael3_all(2,2,nParmSet)) !(2,2,max_parm)
1676 allocate(aad_all(ntyp,2,iparm),bad_all(ntyp,2,nParmSet)) !(ntyp,2,max_parm)
1677 allocate(aa_all(ntyp,ntyp,iparm),bb_all(ntyp,ntyp,nParmSet),&
1678 augm_all(ntyp,ntyp,iparm),eps_all(ntyp,ntyp,nParmSet),&
1679 sigma_all(ntyp,ntyp,iparm),r0_all(ntyp,ntyp,nParmSet),&
1680 chi_all(ntyp,ntyp,iparm)) !(ntyp,ntyp,max_parm)
1681 allocate(chip_all(ntyp,iparm),alp_all(ntyp,nParmSet)) !(ntyp,max_parm)
1682 allocate(ebr_all(iparm),d0cm_all(nParmSet),akcm_all(nParmSet),&
1683 akth_all(iparm),akct_all(nParmSet),v1ss_all(nParmSet),&
1684 v2ss_all(iparm),v3ss_all(nParmSet)) !(max_parm)
1685 allocate(v1sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1686 allocate(v2sccor_all(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,iparm))
1687 !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp,max_parm)
1688 allocate(nlob_all(ntyp1,iparm)) !(ntyp1,max_parm)
1689 allocate(nlor_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
1690 allocate(nterm_all(-maxtor:maxtor,-maxtor:maxtor,2,iparm))
1691 !(-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1692 allocate(ntermd1_all(-maxtor:maxtor,-maxtor:maxtor,&
1693 -maxtor:maxtor,2,iparm))
1694 allocate(ntermd2_all(-maxtor:maxtor,-maxtor:maxtor,&
1695 -maxtor:maxtor,2,iparm))
1696 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2,max_parm)
1697 allocate(nbondterm_all(ntyp,iparm)) !(ntyp,max_parm)
1698 allocate(ithetyp_all(-ntyp1:ntyp1,iparm)) !(-ntyp1:ntyp1,max_parm)
1699 allocate(nthetyp_all(iparm),ntheterm_all(nParmSet),&
1700 ntheterm2_all(iparm),ntheterm3_all(nParmSet),&
1701 nsingle_all(iparm),&
1702 ndouble_all(iparm),nntheterm_all(nParmSet)) !(max_parm)
1703 allocate(nterm_sccor_all(-ntyp:ntyp,-ntyp:ntyp,iparm)) !(-ntyp:ntyp,-ntyp:ntyp,max_parm)
1705 end subroutine alloc_enecalc_arrays
1706 !--------------------------------------------------------------------------
1707 !--------------------------------------------------------------------------