Topic Options
#6265 - 04/18/05 10:08 PM About the PMF program
label Offline
Forum Member

Registered: 12/30/03
Posts: 80
Dear all ,
I obtained the program from prof. T. Allen. But I didn't know how to use it .I posted the program here .
I wish someone tell me how to use it .
THX

PMF.f90
______________________________________________________
program wham ! One dimensional PMF. T.Allen 01/03
! Designed for using the Charmm v30 IUMMFP output files of an MMFP GEO.
! See "FORMAT" below to check or change file formatting.
! Create a file "wham.lst" listing file names and parameters for every window.
! The file should have a line with the filename followed by a line with
! the window number to which the data belongs, z0 and F_constant of each window.
! A .le.0 window number will lead to automatic number generation.
! If the file "wham.lst" is absent, I will try to create it for you
! by assuming files pmf_IW_IR.dat, where IW is the window number and IR is
! the run number, exist.
! All other input parameters are defined in the file "wham.inp".
! See "INPUT" below for description of inputs.
implicit none

integer, parameter :: maxwin=202 ! Windows, runs limits
integer, parameter :: maxbin=5001 ! # bins limit
integer, parameter :: maxiter=10000 ! Max # iterations
integer, parameter :: nstmin=100 ! Ignore nearly empty data files

integer :: nwin ! # of windows
integer :: nsteps(maxwin) ! Total # steps in windows
integer :: nstepp ! Nstep count for Periodic
integer :: nskip ! Omit for equilibration
integer :: testfrq ! Freq. to test convergence,
integer :: nrhob(maxwin,-maxbin:maxbin) ! Biased histograms
integer :: ntotal(maxwin) ! Total biased histogram
integer :: nerr,ierror ! Convergence integers
integer :: nbin,ibin,kbin,npbin ! Bin integers
integer :: ifile,len3w,len3r,iw,ir ! Temporary integers
integer :: iter,i,j,k,t,iwin,jwin ! Temp integers
integer :: jincl,nrun,iwtot ! Temp integers
integer :: kwin ! # of windows temp variable

real*8, parameter :: temp=330.d0 ! K (could make input)
real*8, parameter :: kT=1.381d-23*temp ! Joules
real*8, parameter :: kjmol=1.660578d-21/kT
real*8, parameter :: kcal=4.1868d0*kjmol ! Multiply kcal/mol to get kT
real*8, parameter :: fmix=1.d0 ! Mixing ratio

real*8 :: f(maxwin),f0(maxwin) ! Free eenergy constants
real*8 :: z ! reaction coordinate values
real*8 :: z0(maxwin),fcon(maxwin) ! Window function parameters
real*8 :: zmean(maxwin),z2(maxwin) ! Coordinate statistics
real*8 :: tolw ! Convergence tolerance
real*8 :: test0(-maxbin:maxbin) ! Convergence test
real*8 :: error,errormax ! Convergence errors
real*8 :: zmin,zmax ! Desired reaction coordinate range
real*8 :: dz ! Bin size
real*8 :: rhob(-maxbin:maxbin) ! Biased density
real*8 :: rhob_pr(-maxbin:maxbin) ! Normalised biased density
real*8 :: zbin(-maxbin:maxbin) ! Plotting grid
real*8 :: rho(-maxbin:maxbin) ! Unbiased rho
real*8 :: pmf(-maxbin:maxbin) ! PMF
real*8 :: w(maxwin,-maxbin:maxbin) ! Window functions
real*8 :: ztest,rtest,wminval,sum ! Temporary reals
real*8 :: umin,umax,bmin,bmax ! Temp reals
real*8 :: dz2,zmed,zlen ! Temp reals
real*8 :: zval,testval,tempf,tempfcon ! Temp reals
real*8 :: pi,degrad ! Pi, deg-rad conversion

logical :: winclude(maxwin) ! Include window flags
logical :: lrest ! Restart flag
logical :: readf,oex,epmf,inpx ! File exist flags
logical :: pmfmon ! Monitor the PMF?
logical :: converged ! Has the run converged?
logical :: zerorho,zeromrho ! Zero density tests
logical :: periodic ! Apply periodic boundaries
logical :: nowrap ! Use for angles Not periodic
! I.e. Keep continueous functions, but don't impose wrap around edges.
logical :: angle ! An angular PMF
logical :: symbias ! Symmetrize the biased density?
logical :: debugbias ! Internal debug flag
logical :: lfiles ! Shall time-series be read?

character*80 :: ftest ! IO variables
character*3 :: rtemp,wtemp ! IO variables


debugbias = .true. ! for debug

pi=dacos(-1.d0)
degrad= pi/180.d0

print *, 'WHAM UNBIASING OF MMFP UMBRELLA SAMPLING PMF - T.ALLEN 2002'
print *
if(maxwin>999) stop 'maxwin too high.'

! Defaults for inputs:
dz=0.01d0
zmin=-25.d0
zmax=25.d0
tolw= 0.001d0 ! kcal/mol
testfrq = 100
pmfmon=.false.
nskip=1 ! 1 is for the title line, use more for equilibration
periodic=.false.
angle=.false.
symbias=.false.
nowrap=.false.
! INPUT:
inquire(file='wham.inp',exist=inpx)
if(inpx) then
open(unit=97,file='wham.inp',status='old')
read(97,*)
read(97,*) dz ! Ignored on restart
read(97,*) zmin,zmax ! zmin,zmax: Ignored on restart
read(97,*) tolw ! Tolerance on PMF
read(97,*) testfrq ! Freq. to test convergence
read(97,*) pmfmon ! Output large PMF convergence plot file?
read(97,*) nskip ! Equilibration
read(97,*) periodic ! Periodic boundaries?
read(97,*) angle ! Use deg/kcal/mol/rad^2,
! ! NOT Angstrom,kcal/mol/A^2
read(97,*) symbias ! Symmetrize biased density?
read(97,*,end=1) nowrap ! no periodic with Angles
1 close(unit=97)
endif

write(*,'(a)') 'Input parameters:'
if(inpx) write(*,'(a)') '(Taken from wham.inp)'
write(*,'(a,f16.10)') 'Bin size =',dz
write(*,'(a,2(1x,f12.6))') &
'Plot range =',zmin,zmax
write(*,'(a,f16.10)') 'PMF convergence tolerance =',tolw
write(*,'(a,i8)') 'Test PMF convergence freq. =',testfrq
if(pmfmon) write(*,'(a)') 'Monitoring PMF in pmf.test!'
nskip=nskip+1 ! to allow for title line
write(*,'(a,i8)') 'Nskip equilibration =',nskip
print *, 'periodic = ',periodic
if(periodic) write(*,*) 'Applying periodic boundaries -'
print *, 'angle = ',angle
if(angle) then
write(*,*) 'This is an Anglular PMF!'
endif
if(periodic) then
write(*,'(a,f12.6,a,f12.6)') &
'Periodic range ',zmin,' to ',zmax
zlen=zmax-zmin
endif
if(nowrap) write(*,*) 'Angles with no periodicity.'
if(angle) then
dz=dz*degrad
zmin=zmin*degrad
zmax=zmax*degrad
zlen=zlen*degrad
endif
print *, 'symbias = ',symbias
if(symbias) write(*,*) 'Symmetrizing biased density!'

write(*,*)

! Is this a restart?
! If BOTh fe and fcon files exist, a restart will occur:
lrest=.true.
inquire(file='fcon.dat',exist=readf)
if(.not. readf) lrest=.false.
inquire(file='fe.dat',exist=readf)
if(.not. readf) lrest=.false.
! If bias.dat does NOT exist, the timeseries will be read again
lfiles=.false.
inquire(file='bias.dat',exist=readf)
if(.not. readf) lfiles=.true.
! Note, a restart can occur without pmf.dat.
! However, if NOT a restart, then we MUST generate bias.dat from scratch:
if(.not. lrest) lfiles=.true.

! Read restart files if requested:
if(lrest) then
write(*,'(a)')'This is a Restart run!'
write(*,*)
print *, 'Reading Force constants from file...'
open(unit=88,file='fcon.dat',status='old')
read(88,61,end=21,err=21) nwin
print *, 'Number of windows = ',nwin
if(nwin.gt.maxwin) Stop 'Nwin is too large!'
iwin=0
20 iwin=iwin+1
read(88,62,end=21,err=21) &
jwin,z0(iwin),tempfcon,nsteps(iwin),jincl
! Note that z0 already corrected for a periodic system
61 format(1x,i6)
62 format(1x,i6,2(2x,e16.8),2x,i10,i1)
if(jincl.eq.1) then
winclude(iwin)=.true.
write(*,26) iwin,nsteps(iwin),z0(iwin),tempfcon
else
winclude(iwin)=.false.
write(*,23) iwin,nsteps(iwin),z0(iwin),tempfcon
endif
fcon(iwin)=kcal*tempfcon ! now in kT/A^2
if(angle) then
z0(iwin)=z0(iwin)*degrad ! now in radians
endif
goto 20
21 close(unit=88)
if(iwin-1.ne.nwin) stop 'Wrong number of entries in fcon.dat.'

27 format(a20)
do i=1,nwin
f(i)=0.d0
enddo
print *, 'Reading FE constants from file...'
open(unit=88,file='fe.dat',status='old')
iwin=0
30 iwin=iwin+1
read(88,66,end=31,err=31) jwin,tempf
66 format(1x,i6,2x,e16.8)
f(iwin)=kcal*tempf
goto 30
31 close(unit=88)
if(iwin-1.ne.nwin) stop 'Wrong number of entries in fe.dat.'

! If available, Use final PMF as reference for convergence on first test:
do ibin=-maxbin,maxbin
test0(ibin)=0.d0
enddo
npbin=0
inquire(file='pmf.dat',exist=oex)
if(oex) then
print *, 'Reading last PMF for first convergence test...'
open(unit=9,file='pmf.dat',status='old')
read(9,108,end=45,err=45) npbin
if(npbin.gt.maxbin) Stop 'npbin is too large!'
do ibin=-npbin,npbin
read(9,120,end=45,err=45) zval,testval
test0(ibin)=testval*kcal ! now in kT
enddo
45 close(unit=9)
endif

if(.not.lfiles) then
print *, 'Reading biased density And Grid from file...'
open(unit=89,file='bias.dat',status='old')
read(89,64,end=2,err=2) nbin,dz ! dz overrides any input
if(angle) then
dz=dz*degrad
endif
if(nbin.gt.maxbin) Stop 'Nbin is too large!'
if(npbin.gt.0 .and. npbin.ne.nbin) Stop 'npbin .ne. nbin!'
do i=-nbin,nbin
read(89,68,end=2,err=2) j,zbin(i),rhob(i)
if(angle) then
zbin(i)=zbin(i)*degrad ! now in radians
endif
enddo
print *, 'Total bins read = ', 2*nbin+1
goto 3
2 stop 'Error reading restart biased density.'
3 close(unit=89)
endif

endif ! on restart test

! If no bias.dat, we must generate histograms:
if(lfiles) then
print *,'Generating biased density from time-series...'
print *

! Create bins for PMF and find biased probabilities n_i rho_i _b:
print *,'Using bin size = ',dz
print *
nbin=int((zmax-zmin-dz)/2.d0/dz) +1
zmed=0.5d0*(zmax+zmin)
bmin=zmed-dfloat(nbin)*dz
bmax=zmed+dfloat(nbin)*dz
if(angle) then
print *, 'Position of centres of 1st and last bins = ', &
bmin/degrad,bmax/degrad
else
print *, 'Position of centres of 1st and last bins = ', &
bmin,bmax
endif
! print *, 'DEBUG: bmin-zmin = ',bmin-zmin
! print *, 'DEBUG: bmax-zmax = ',bmax-zmax
if(dabs(bmin-zmin).gt.1.d-5) stop 'Grid definition error.'
if(dabs(bmax-zmax).gt.1.d-5) stop 'Grid definition error.'
print *, 'Total number of bins = ', 2*nbin+1
print *
if(nbin.gt.maxbin) stop 'Too many bins!'
do ibin=-nbin,nbin
zbin(ibin)=zmin + dfloat(ibin+nbin)*dz
enddo
if(periodic) then
! Put the z0 in the right range:
do i=1,nwin
if(z0(i).gt.zmax) z0(i)= z0(i)-zlen
if(z0(i).lt.zmin) z0(i)= z0(i)+zlen
enddo
endif

! Check/create file list:
inquire(file='wham.lst',exist=inpx)

if(.not.inpx) then
! There is no file list, so we try to create it:
open(unit=12,file='wham.lst',status='new')
do iw=1,maxwin
call nam3(iw,wtemp,len3w)
do ir=1,99 ! Can change this to isolate runs
call nam3(ir,rtemp,len3r)
ftest='pmf_'//trim(wtemp(1:len3w))//'_'// &
trim(rtemp(1:len3r))//'.dat'
inquire(file=ftest,exist=oex)
if(oex) then
open(unit=8,file=ftest,status='old')
read(8,8,end=13,err=13) z0(iw),fcon(iw)
8 format(10x,f14.8,5x, f14.8)
write(12,11) ftest
write(12,12) iw,z0(iw),fcon(iw)
13 close(unit=8)
11 format(a80)
12 format(i6,1x,f12.8,1x,f12.8)
endif
enddo
enddo
close(unit=12)
endif

! Searching and reading files:
nwin=0 ! is the number of the last included window
umin=9.d9
umax=-9.d9
iwtot=0
do iw=1,maxwin
nsteps(iw)=0
winclude(iw)=.false.
zmean(iw)=0.d0
z2(iw)=0.d0
do ibin=-nbin,nbin
nrhob(iw,ibin)=0
enddo
enddo
open(unit=12,file='wham.lst',status='old')
16 read(12,11,end=19,err=19) ftest
inquire(file=trim(ftest),exist=oex)
if(.not.oex) goto 16
open(unit=8,file=ftest,status='old')
read(12,*,end=19,err=19) iwin,zval,tempfcon
if(iwin.le.0) iwin=nwin+1 ! Give it a new window if unknown
if(iwin.gt.nwin) nwin=iwin
z0(iwin)=zval
if(angle) then
z0(iwin)=z0(iwin)*degrad ! Now in radians
endif
fcon(iwin)=tempfcon
do k=1,nskip
read(8,*,end=200)
enddo
! FORMAT:
10 read(8,192,end=200,err=200) z
192 format(11x,f14.8)
nsteps(iwin)=nsteps(iwin)+1
if(angle) z=z*degrad ! Now in radians
if(periodic) then
! Put the data in the right range:
if(z.gt.zmax) z=z-zlen
if(z.lt.zmin) z=z+zlen
endif
if(z.lt.umin) umin=z
if(z.gt.umax) umax=z
zmean(iwin)=zmean(iwin)+z
z2(iwin)=z2(iwin)+z**2
kbin=idnint((z-zmin)/dz-nbin)
if(iabs(kbin).gt.nbin) stop 'Data exists out of range!'
nrhob(iwin,kbin)=nrhob(iwin,kbin)+1
goto 10
200 close(unit=8)
print *, trim(ftest)
goto 16 ! next file test
19 close(unit=12)

print *
print *, 'Summary of Windows: '
do iw=1,nwin
zmean(iw)=zmean(iw)/dfloat(nsteps(iw))
z2(iw)=dsqrt(z2(iw)/dfloat(nsteps(iw))-zmean(iw)**2)
if(nsteps(iw).ge.nstmin) then
winclude(iw)=.true.
else
iwtot=iwtot+1
print *, 'Window does not contain enough data: ',iw,' excluded.'
endif
if(angle) then
if(winclude(iw)) then
write(*,25) iw,nsteps(iw),z0(iw)/degrad,fcon(iw), &
zmean(iw)/degrad,z2(iw)/degrad
else
write(*,23) iw,nsteps(iw),z0(iw)/degrad,fcon(iw)
endif
else
if(winclude(iw)) then
write(*,25) iw,nsteps(iw),z0(iw),fcon(iw), &
zmean(iw),z2(iw)
else
write(*,23) iw,nsteps(iw),z0(iw),fcon(iw)
endif
endif
fcon(iw)=fcon(iw)*kcal ! now in kT/A^2
enddo

23 format(1x,'Win=',i3,', steps=',i6,', z0=',f8.4, &
', K=',f10.6,' EXCLUDED')
25 format(1x,'Win=',i3,', steps=',i6,', z0=',f8.4, &
', K=',f10.6,' z=',f8.4,'+-',f8.4)
26 format(1x,'Win=',i3,', steps=',i6,', z0=',f8.4, &
', K=',f10.6)
print *
print *, 'Number of last window = ',nwin
print *
print *, 'Observed range: ',umin,umax
if(umin.lt.zmin) stop 'Trajectory < zmin exists!'
if(umax.gt.zmax) stop 'Trajectory > zmax exists!'

print *
print *, 'Number windows excluded due to lack of data = ',iwtot
print *, 'Total number windows included in calculation = ',nwin-iwtot
!
if(symbias) then
print *
print *, 'Symmetrizing biased density by duplicating windows...'
if(2*nwin.gt.maxwin) stop 'maxwin too small for symbias.'
kwin=nwin
do i=1,kwin
nwin=nwin+1
z0(nwin)=-z0(i)
fcon(nwin)=fcon(i)
nsteps(nwin)=nsteps(i)
winclude(nwin)=winclude(i)
do ibin=-nbin,nbin
nrhob(nwin,ibin)=nrhob(i,-ibin)
enddo
enddo
print *, 'Total number of windows is now = ',nwin
print *
endif

! Write restart file force constants:
print *, 'Writing window information to file...'
open(unit=88,file='fcon.dat',status='unknown')
write(88,61) nwin
do i=1,nwin
jincl=0
if(winclude(i)) jincl=1
if(angle) then
write(88,62) i,z0(i)/degrad,fcon(i)/kcal,nsteps(i),jincl
else
write(88,62) i,z0(i),fcon(i)/kcal,nsteps(i),jincl
endif
enddo
close(unit=88)
print *
print *, 'Finding total biased density...'
open(unit=89,file='bias.dat',status='unknown')
if(angle) then
write(89,64) nbin, dz/degrad
else
write(89,64) nbin, dz
endif
do ibin=-nbin,nbin
rhob(ibin)=0.d0
do i=1,nwin
if(winclude(i)) rhob(ibin)=rhob(ibin)+dfloat(nrhob(i,ibin))
enddo
if(angle) then
write(89,68) ibin,zbin(ibin)/degrad,rhob(ibin)
else
write(89,68) ibin,zbin(ibin),rhob(ibin)
endif
enddo
close(unit=89)
64 format(1x,i8,2x,f12.8)
68 format(1x,i6,2(2x,e16.8))
endif ! on restart

! Actual biased density for printing
sum=0.d0
do j=1,nwin
if(winclude(j)) sum=sum+nsteps(j)
enddo
do ibin=-nbin,nbin
rhob_pr(ibin)=rhob(ibin)/sum
enddo
print *
print *, 'Finding window functions...'
do ibin=-nbin,nbin
do i=1,nwin
if(winclude(i)) then
if(periodic .and. .not. nowrap) then
zval=zbin(ibin)-z0(i)
if(zval.gt.0.5d0*zlen) zval=zval-zlen
if(zval.lt.-0.5d0*zlen) zval=zval+zlen
! Assume all data within one period of centre
w(i,ibin)=0.5d0*fcon(i)*zval**2
else
w(i,ibin)=0.5d0*fcon(i)*(zbin(ibin)-z0(i))**2
endif
endif
enddo
enddo

if(debugbias) then
do i=1,nwin
call nam3(i,wtemp,len3w)
open(unit=80,file='nrhob.'//trim(wtemp(1:len3w)),status='unknown')
do ibin=-nbin,nbin
write(80,33) zbin(ibin),nrhob(i,ibin)
33 format(1x,f12.6,2x,i8)
enddo
close(unit=80)
enddo
endif

! Iterate equations:
print *
print *, 'Iterating equations (kill at anytime and restart)...'
print *
print *,'Using PMF tolerance (kcal/mol) = ',tolw
print *,'Test convergence each = ',testfrq,' steps.'
print *
tolw = tolw*kcal ! Now in kT
iter=0
converged=.false.
100 iter=iter+1
! Set up guess at unbiased density:
do ibin=-nbin,nbin
sum=0.d0
do j=1,nwin
if(winclude(j)) then
! To prevent underflow:
if(w(j,ibin)-f(j).le.100.d0) then
sum=sum+nsteps(j)*dexp(-(w(j,ibin)-f(j)))
endif
endif
enddo
! To prevent divi by zero:
if(sum.lt.1.d-100) then
rho(ibin)=0.d0 ! because rhob(ibin) will be absolutely zero
if(rhob(ibin).gt.1.d-8) then
print *, 'ibin,rhob = ',ibin,rhob(ibin)
stop 'Why isnt rhob zero here?'
endif
else
rho(ibin)=rhob(ibin)/sum
endif
enddo ! on ibin loop

if(mod(iter,testfrq).eq.0) then
! Test for convergence on W(z):
errormax=0.d0
nerr=0
wminval=9.d9
converged=.true.
if(pmfmon) then
! (re-)Open file that will have plot of PMF every testfrq steps.
inquire(file='pmf.test',exist=epmf)
if(epmf) then
open(unit=29,file='pmf.test',position='append',status='old')
else
open(unit=29,file='pmf.test',status='new')
endif
endif
! Open plotting PMF for restart and final result.
open(unit=9,file='pmf.dat',status='unknown')
write(9,110) nbin
110 format('# Z ',3x,' PMF(kcal/mol): nbin= ',i8)
108 format(34x,i8)
do ibin=-nbin,nbin
if(rho(ibin).ge.1.d-20) then
pmf(ibin)=-dlog(rho(ibin)/dz)
else
pmf(ibin)=100.d0 * kcal ! So limit will be 100 kcal on print
endif
if(pmf(ibin).lt. wminval) wminval=pmf(ibin)
enddo
do ibin=-nbin,nbin
pmf(ibin)=pmf(ibin)-wminval
if(angle) then
if(pmfmon) write(29,120) zbin(ibin)/degrad,pmf(ibin)/kcal
write(9,120) zbin(ibin)/degrad,pmf(ibin)/kcal ! now in kcal/mol
else
if(pmfmon) write(29,120) zbin(ibin),pmf(ibin)/kcal ! in kcal/mol
write(9,120) zbin(ibin),pmf(ibin)/kcal ! now in kcal/mol
endif
error=pmf(ibin)-test0(ibin)
if(dabs(error)>tolw) nerr=nerr+1
if(dabs(error)>dabs(errormax)) then
errormax=error
ierror=ibin
endif
test0(ibin)=pmf(ibin)
enddo
if(pmfmon) write(29,*)
if(pmfmon) close(unit=29)
close(unit=9)
write(*,92) iter,ierror,pmf(ierror)/kcal,errormax/kcal,nerr
92 format(1x,i8,' Bin ',i6,' Worst: ',f10.6,' +- ',f10.6,' #bad= ',i6)
if(nerr.gt.0) converged=.false.

! Write current FE constants to file in case of kill (overwrite):
open(unit=88,file='fe.dat',status='unknown')
do i=1,nwin
jincl=0
if(winclude(i)) jincl=1
write(88,66) i,f(i)/kcal
enddo
close(unit=88)
endif ! testfrq test

! If not converged we must iterate the equations:
if(.not. converged) then

! Integrate to find new guess at F:
do i=1,nwin
f0(i)=f(i)
sum=0.5d0*(rho(-nbin)*dexp(-w(i,-nbin))+rho(nbin)*dexp(-w(i,nbin)))
do ibin=-nbin+1,nbin-1
if(winclude(i).and. w(i,ibin).le.100.d0) then
sum=sum+rho(ibin)*dexp(-w(i,ibin))
endif
enddo
f(i)=-fmix*dlog(sum)+(1.d0-fmix)*f0(i)
enddo

if(iter goto 100
else
stop 'Solution not converged in maxiter iteractions.'
endif

endif ! not converged

! Converged:
print *
print *, 'Solution converged in ',iter,' iterations.'
print *

! Write rho:
if(angle) then
! Convert bins back to degree for the remainder of program:
do ibin=-nbin,nbin
zbin(ibin)=zbin(ibin)/degrad
enddo
endif
open(unit=10,file='rho.dat',status='unknown')
open(unit=11,file='rho_biased.dat',status='unknown')
write(10,112)
write(11,112)
do ibin=-nbin,nbin
rho(ibin)=rho(ibin)/dz
rhob_pr(ibin)=rhob_pr(ibin)/dz
write(10,120) zbin(ibin),rho(ibin)
write(11,120) zbin(ibin),rhob_pr(ibin)
enddo
112 format('# Z ',3x,'Rho(units^-1)')
120 format(1x,f10.4,3x,f10.4)
close(unit=9)
close(unit=10)
close(unit=11)
! Find symmetrised PMF to write:
open(unit=9,file='pmf.sym',status='unknown')
open(unit=10,file='rho.sym',status='unknown')
open(unit=11,file='rho_biased.sym',status='unknown')
write(9,110) nbin
write(10,112)
write(11,112)
do ibin=0,nbin
! The zerorho vars ensure symmetrizing is avoided with missing data
zerorho=.true.
zeromrho=.true.
if(rho(ibin).gt.1.d-20) zerorho=.false.
if(rho(-ibin).gt.1.d-20) zeromrho=.false.
if(zerorho) then
if(.not. zeromrho) then
write(9,120) zbin(ibin),pmf(-ibin)/kcal
write(10,120) zbin(ibin),rho(-ibin)
endif
else
if(zeromrho) then
write(9,120) zbin(ibin),pmf(ibin)/kcal
write(10,120) zbin(ibin),rho(ibin)
else
write(9,120) zbin(ibin),(pmf(ibin)+pmf(-ibin))*0.5d0/kcal
write(10,120) zbin(ibin),(rho(ibin)+rho(-ibin))*0.5d0
endif
endif
write(11,120) zbin(ibin),(rhob_pr(ibin)+rhob_pr(-ibin))*0.5d0
enddo
close(unit=9)
close(unit=10)
close(unit=11)
stop 'Calculation complete.'
end


subroutine nam3(num,fch,len3)
implicit none
integer, intent(in) :: num
character*1 :: c(3)
character*3, intent(out) :: fch
integer :: n1,k(3),i
integer, intent(out) :: len3
c(1)=' '
c(2)=' '
c(3)=' '
if(num.gt.999) stop 'Num > 999 in nam3!'
k(1)=num/100
n1=num-k(1)*100
k(2)=n1/10
k(3)=n1-k(2)*10
do i=1,3
if(k(1).eq.0) then
k(1)=k(2)
k(2)=k(3)
k(3)=-1
endif
enddo
len3=0
if(k(1).ge.0) then
c(1)=char(48+k(1))
len3=1
endif
if(k(2).ge.0) then
c(2)=char(48+k(2))
len3=2
endif
if(k(3).ge.0) then
c(3)=char(48+k(3))
len3=3
endif
fch=c(1)//c(2)//c(3)
if(fch.eq.' ') then
fch='0 '
len3=1
endif
return
end subroutine nam3

Top
#6266 - 05/29/05 10:01 PM Re: About the PMF program [Re: label]
charmm31 Offline
Forum Member

Registered: 03/17/05
Posts: 31
Loc: Singapore
Hi,

What this program is used for? Does PMF means Potential of Mean Force?

Thanks,

Top
#6267 - 05/30/05 12:07 PM Re: About the PMF program [Re: label]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8407
Loc: 39 03 48 N, 77 06 54 W
With all the formatting removed, posting a Fortran program is somewhat problematic; you should consider using the HTML markup tags, or else post such items as an Attachment. This isn't necessarily the correct forum for this, either.
_________________________
Rick Venable
computational chemist


Top
#6268 - 03/30/06 12:35 PM Re: About the PMF program [Re: rmv]
noskov Offline
Forum Member

Registered: 06/18/04
Posts: 162
Loc: Calgary, AB
This staff is doing WHAM analysis in presence of the biasing force applied through the MMFP. Works fine. PMF still should be obtained through MD simulation for many windows.


Edited by noskov (03/30/06 12:38 PM)
_________________________
Sergei Noskov snoskov@ucalgary.ca Dept.BioSciences University of Calgary

Top

Moderator:  chmgr, John Legato, petrella