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

   1       subroutine PMFread
   2 c Read the PMFs from wham
   3       implicit none
   4       include 'DIMENSIONS'
   5       include 'DIMENSIONS.PMF'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.MD'
   8 #ifdef LANG0
   9 #ifdef FIVEDIAG
  10       include 'COMMON.LANGEVIN.lang0.5diag'
  11 #else
  12       include 'COMMON.LANGEVIN.lang0'
  13 #endif
  14 #else
  15       include 'COMMON.LANGEVIN'
  16 #endif
  17       include 'COMMON.QRESTR'
  18       include 'COMMON.PMF'
  19       include 'COMMON.REMD'
  20       integer i,iumb,iiset,j,t,nbin
  21       double precision beta_h,qtemp,htemp
  22       read(inp,*) delta_q
  23       write (iout,*)"PMFread: nT",nrep," delta_q",delta_q," nW",nset
  24       print *,(remd_t(i),i=1,nrep)
  25       do j=1,nrep
  26         do iumb=1,nset
  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",
  31      &         iiset,iumb
  32             stop
  33           endif
  34           if (beta_h.ne.remd_t(j)) then
  35             write (iout,*)
  36      &     "Error replica temperatures do not match PMF temperatures"
  37             write (iout,*) beta_h,remd_T(j)
  38             stop
  39           endif
  40           do i=1,nbin
  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)
  47           enddo ! i
  48         enddo ! iumb
  49       enddo ! j
  50       do iumb=1,nset
  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
  55           do j=1,nrep
  56             if (i.lt.tmin(j,iumb).or.i.gt.tmax(j,iumb)) then
  57               write (iout,'("    ------",$)')
  58             else
  59               write (iout,'(f10.3$)') PMFtab(i,j,iumb)
  60             endif
  61           enddo
  62           write (iout,*)
  63         enddo ! i
  64       enddo ! iumb
  65       return
  66       end
  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!
  71       implicit none
  72       include 'DIMENSIONS'
  73       include 'DIMENSIONS.PMF'
  74 #ifdef MPI
  75       include 'mpif.h'
  76       include 'COMMON.SETUP'
  77       integer IERROR,ERRCODE
  78 #endif
  79       include 'COMMON.IOUNITS'
  80       include 'COMMON.MD'
  81       include 'COMMON.QRESTR'
  82       include 'COMMON.REMD'
  83       include 'COMMON.PMF'
  84       integer i,iqmin,iqmax,irep
  85       double precision q,qmin,qmax,ePMF,ePMF_q
  86 c Determine the location of the q
  87       irep=1
  88       do while (dabs(t_bath-remd_t(irep)).gt.1.0d-4)
  89         irep=irep+1
  90       enddo
  91       if (irep.gt.nrep) then
  92         write (iout,*) "CHUJ NASTAPIL in PMF_energy!!!"
  93         write (iout,*)
  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..."
  98         call flush(iout)
  99 #ifdef MPI
 100         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 101 #endif
 102         stop
 103       endif
 104       iqmin=tmin(irep,iset)
 105       iqmax=tmax(irep,iset)
 106       qmin=iqmin*delta_q
 107       qmax=iqmax*delta_q
 108 #ifdef DEBUG
 109       write (iout,*) "PMF_energy q",q," qmin",qmin," qmax",qmax,
 110      & " irep",irep," iset",iset," iqmin",iqmin," iqmax",iqmax
 111 #endif
 112       if (q.le.qmin) then
 113         ePMF_q=(PMFtab(iqmin+1,irep,iset)-PMFtab(iqmin,irep,iset))
 114      &     /delta_q
 115         ePMF=PMFtab(iqmin,irep,iset)+ePMF_q*(q-qmin)
 116 #ifdef DEBUG
 117        write (iout,*) "q<=qmin ePMF",ePMF," ePMF_q",ePMF_q
 118 #endif
 119       else if (q.ge.qmax) then
 120         ePMF_q=(PMFtab(iqmax,irep,iset)-PMFtab(iqmax-1,irep,iset))/
 121      &      delta_q
 122         ePMF=PMFtab(iqmax,irep,iset)+ePMF_q*(q-qmax)
 123 #ifdef DEBUG
 124        write (iout,*) "q>=qmax ePMF",ePMF," ePMF_q",ePMF_q
 125 #endif
 126       else
 127         do i=iqmin+1,iqmax
 128           qmax=i*delta_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)
 132 #ifdef DEBUG
 133        write (iout,*) "qmin<q<qmax bin",i," ePMF",ePMF," ePMF_q",ePMF_q
 134 #endif
 135             exit
 136           endif
 137           qmin=qmax
 138         enddo
 139       endif
 140       return
 141       end