source/wham/src-HCD-5D/PMFprocess.F

   1       subroutine PMFread(*)
   2 c Read the PMFs from wham
   3       implicit none
   4       include 'DIMENSIONS'
   5       include 'DIMENSIONS.FREE'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.FREE'
   8       include 'COMMON.PMF'
   9       integer i,iumb,iiset,j,t,nbin,iparm,nRmax
  10       double precision beta_h_temp,qtemp,htemp
  11       read(inp,*,err=10,end=10) delta_q
  12       write(iout,*) "delta_q",delta_q
  13
  14       do iparm = 1, nparmset
  15
  16       write(iout,*) "PMFread: iparm",iparm," nT",nT_h(iparm),
  17      &  " nR",nR(:nT_h(iparm),iparm)
  18 c      print *,(beta_h(i),i=1,nT(iparm))
  19       do j=1,nT_h(iparm)
  20         do iumb=1,nR(j,iparm)
  21           read (inp,*,end=10,err=10) iiset,beta_h_temp,nbin
  22           write (iout,*) iiset,beta_h_temp,nbin
  23           if (iiset.ne.iumb) then
  24             write(iout,*) "Error: inconsistency in US windows",
  25      &         iiset,iumb
  26             return1
  27           endif
  28           beta_h_temp=1.0d0/(0.001987*beta_h_temp)
  29           if (dabs(beta_h_temp-beta_h(j,iparm)).gt.1.0d-6) then
  30             write (iout,*)
  31      &     "Error replica temperatures do not match PMF temperatures"
  32             write (iout,*) 1.0d0/(0.001987*beta_h_temp),
  33      &      1.0d0/(0.001987*beta_h(j,iparm))
  34             stop
  35           endif
  36           do i=1,nbin
  37             read (inp,*,end=10,err=10) qtemp,htemp
  38             t = int(qtemp/delta_q+1.0d-4)
  39             write (iout,*) qtemp,t,htemp
  40             if (i.eq.1) tmin(j,iumb,iparm)=t
  41             if (i.eq.nbin) tmax(j,iumb,iparm)=t
  42             PMFtab(t,j,iumb,iparm)=dlog(htemp)/beta_h_temp
  43           enddo ! i
  44         enddo ! iumb
  45       enddo ! j
  46
  47       nRmax=nR(1,iparm)
  48       do i=2,nT_h(iparm)
  49         if (nR(i,iparm).gt.nRmax) nRmax=nR(i,iparm)
  50       enddo
  51       do iumb=1,nR(j,iparm)
  52         write (iout,*)"Input PMFs, restraint",iumb,
  53      &   " q0",q0(1,iumb,:nT_h(iparm),iparm)
  54         write (iout,'(5x,20f10.1)') (1.0d0/(0.001987*beta_h(j,iparm)),
  55      &   j=1,nT_h(iparm))
  56         do i=0,int(1.0/delta_q)
  57           write (iout,'(f5.2,$)') i*delta_q
  58           do j=1,nT_h(iparm)
  59             if (i.lt.tmin(j,iumb,iparm).or.i.gt.tmax(j,iumb,iparm)) then
  60               write (iout,'("    ------",$)')
  61             else
  62               write (iout,'(f10.3$)') PMFtab(i,j,iumb,iparm)
  63             endif
  64           enddo
  65           write (iout,*)
  66         enddo ! i
  67       enddo ! iumb
  68
  69       enddo ! iparm
  70       return
  71    10 return1
  72       end
  73 c------------------------------------------------------------------------
  74       subroutine PMF_energy(q,irep,iset,iparm,ePMF,ePMF_q)
  75 c Calculate the energy and derivative in q due to the biasing PMF
  76 c Caution! Only ONE q is handled, no multi-D q-restraints available!
  77       implicit none
  78       include 'DIMENSIONS'
  79       include 'DIMENSIONS.FREE'
  80       include 'COMMON.IOUNITS'
  81       include 'COMMON.FREE'
  82       include 'COMMON.PMF'
  83       integer i,iqmin,iqmax,irep,iset,iparm
  84       double precision q,qmin,qmax,ePMF,ePMF_q
  85 c Determine the location of the q
  86       iqmin=tmin(irep,iset,iparm)
  87       iqmax=tmax(irep,iset,iparm)
  88       qmin=iqmin*delta_q
  89       qmax=iqmax*delta_q
  90 #ifdef DEBUG
  91       write (iout,*) "PMF_energy q",q," qmin",qmin," qmax",qmax,
  92      & " irep",irep," iset",iset
  93 #endif
  94       if (q.le.qmin) then
  95         ePMF_q=(PMFtab(iqmin+1,irep,iset,iparm)-
  96      &     PMFtab(iqmin,irep,iset,iparm))/delta_q
  97         ePMF=PMFtab(iqmin,irep,iset,iparm)+ePMF_q*(q-qmin)
  98 #ifdef DEBUG
  99        write (iout,*) "q<=qmin ePMF",ePMF," ePMF_q",ePMF_q
 100 #endif
 101       else if (q.ge.qmax) then
 102         ePMF_q=(PMFtab(iqmax,irep,iset,iparm)-
 103      &      PMFtab(iqmax-1,irep,iset,iparm))/delta_q
 104         ePMF=PMFtab(iqmax,irep,iset,iparm)+ePMF_q*(q-qmax)
 105 #ifdef DEBUG
 106        write (iout,*) "q>=qmax ePMF",ePMF," ePMF_q",ePMF_q
 107 #endif
 108       else
 109         do i=iqmin+1,iqmax
 110           qmax=i*delta_q
 111           if (q.ge.qmin .and. q.le.qmax) then
 112             ePMF_q=(PMFtab(i,irep,iset,iparm)-
 113      &       PMFtab(i-1,irep,iset,iparm))/delta_q
 114             ePMF=PMFtab(i-1,irep,iset,iparm)+ePMF_q*(q-qmin)
 115 #ifdef DEBUG
 116        write (iout,*) "qmin<q<qmax bin",i," ePMF",ePMF," ePMF_q",ePMF_q
 117 #endif
 118             exit
 119           endif
 120           qmin=qmax
 121         enddo
 122       endif
 123       return
 124       end