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
