2 c Read the PMFs from wham
5 include 'DIMENSIONS.PMF'
6 include 'COMMON.IOUNITS'
10 include 'COMMON.LANGEVIN.lang0.5diag'
12 include 'COMMON.LANGEVIN.lang0'
15 include 'COMMON.LANGEVIN'
17 include 'COMMON.QRESTR'
20 integer i,iumb,iiset,j,t,nbin
21 double precision beta_h,qtemp,htemp
23 write (iout,*)"PMFread: nT",nrep," delta_q",delta_q," nW",nset
24 print *,(remd_t(i),i=1,nrep)
27 read (inp,*) iiset,beta_h,nbin
28 write (iout,*) iiset,beta_h,nbin
29 if (iiset.ne.iumb) then
30 write(iout,*) "Error: inconsistency in US windows",
34 if (beta_h.ne.remd_t(j)) then
36 & "Error replica temperatures do not match PMF temperatures"
37 write (iout,*) beta_h,remd_T(j)
41 read (inp,*) qtemp,htemp
42 t = int(qtemp/delta_q+1.0d-4)
43 write (iout,*) qtemp,t,htemp
44 if (i.eq.1) tmin(j,iumb)=t
45 if (i.eq.nbin) tmax(j,iumb)=t
46 PMFtab(t,j,iumb)=dlog(htemp)*Rb*remd_t(j)
51 write (iout,*)"Input PMFs, restraint",iumb," q0",qinfrag(1,iumb)
52 write (iout,'(5x,20f10.1)') (remd_T(j),j=1,nrep)
53 do i=0,int(1.0/delta_q)
54 write (iout,'(f5.2,$)') i*delta_q
56 if (i.lt.tmin(j,iumb).or.i.gt.tmax(j,iumb)) then
57 write (iout,'(" ------",$)')
59 write (iout,'(f10.3$)') PMFtab(i,j,iumb)
67 c------------------------------------------------------------------------
68 subroutine PMF_energy(q,ePMF,ePMF_q)
69 c Calculate the energy and derivative in q due to the biasing PMF
70 c Caution! Only ONE q is handled, no multi-D q-restraints available!
73 include 'DIMENSIONS.PMF'
76 include 'COMMON.SETUP'
77 integer IERROR,ERRCODE
79 include 'COMMON.IOUNITS'
81 include 'COMMON.QRESTR'
84 integer i,iqmin,iqmax,irep
85 double precision q,qmin,qmax,ePMF,ePMF_q
86 c Determine the location of the q
88 do while (dabs(t_bath-remd_t(irep)).gt.1.0d-4)
91 if (irep.gt.nrep) then
92 write (iout,*) "CHUJ NASTAPIL in PMF_energy!!!"
94 & "Bath temperature does not match that of any replica"
95 write (iout,*) "t_bath",t_bath
96 write (iout,*) "remd_T",(remd_T(:nrep))
97 write (iout,*) "Terminating..."
100 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
104 iqmin=tmin(irep,iset)
105 iqmax=tmax(irep,iset)
109 write (iout,*) "PMF_energy q",q," qmin",qmin," qmax",qmax,
110 & " irep",irep," iset",iset," iqmin",iqmin," iqmax",iqmax
113 ePMF_q=(PMFtab(iqmin+1,irep,iset)-PMFtab(iqmin,irep,iset))
115 ePMF=PMFtab(iqmin,irep,iset)+ePMF_q*(q-qmin)
117 write (iout,*) "q<=qmin ePMF",ePMF," ePMF_q",ePMF_q
119 else if (q.ge.qmax) then
120 ePMF_q=(PMFtab(iqmax,irep,iset)-PMFtab(iqmax-1,irep,iset))/
122 ePMF=PMFtab(iqmax,irep,iset)+ePMF_q*(q-qmax)
124 write (iout,*) "q>=qmax ePMF",ePMF," ePMF_q",ePMF_q
129 if (q.ge.qmin .and. q.le.qmax) then
130 ePMF_q=(PMFtab(i,irep,iset)-PMFtab(i-1,irep,iset))/delta_q
131 ePMF=PMFtab(i-1,irep,iset)+ePMF_q*(q-qmin)
133 write (iout,*) "qmin<q<qmax bin",i," ePMF",ePMF," ePMF_q",ePMF_q