source/unres/src-5hdiag-tmp/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       include 'COMMON.QRESTR'
   9 #ifndef LANG0
  10       include 'COMMON.LANGEVIN'
  11 #else
  12       include 'COMMON.LANGEVIN.lang0'
  13 #endif
  14       include 'COMMON.PMF'
  15       include 'COMMON.REMD'
  16       integer i,iumb,iiset,j,t,nbin
  17       double precision beta_h,qtemp,htemp
  18       read(inp,*) delta_q
  19       write (iout,*)"PMFread: nT",nrep," delta_q",delta_q," nW",nset
  20       print *,(remd_t(i),i=1,nrep)
  21       do j=1,nrep
  22         do iumb=1,nset
  23           read (inp,*) iiset,beta_h,nbin
  24           write (iout,*) iiset,beta_h,nbin
  25           if (iiset.ne.iumb) then
  26             write(iout,*) "Error: inconsistency in US windows",
  27      &         iiset,iumb
  28             stop
  29           endif
  30           if (beta_h.ne.remd_t(j)) then
  31             write (iout,*)
  32      &     "Error replica temperatures do not match PMF temperatures"
  33             write (iout,*) beta_h,remd_T(j)
  34             stop
  35           endif
  36           do i=1,nbin
  37             read (inp,*) 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)=t
  41             if (i.eq.nbin) tmax(j,iumb)=t
  42             PMFtab(t,j,iumb)=dlog(htemp)*Rb*remd_t(j)
  43           enddo ! i
  44         enddo ! iumb
  45       enddo ! j
  46       do iumb=1,nset
  47         write (iout,*)"Input PMFs, restraint",iumb," q0",qinfrag(1,iumb)
  48         write (iout,'(5x,20f10.1)') (remd_T(j),j=1,nrep)
  49         do i=0,int(1.0/delta_q)
  50           write (iout,'(f5.2,$)') i*delta_q
  51           do j=1,nrep
  52             if (i.lt.tmin(j,iumb).or.i.gt.tmax(j,iumb)) then
  53               write (iout,'("    ------",$)')
  54             else
  55               write (iout,'(f10.3$)') PMFtab(i,j,iumb)
  56             endif
  57           enddo
  58           write (iout,*)
  59         enddo ! i
  60       enddo ! iumb
  61       return
  62       end
  63 c------------------------------------------------------------------------
  64       subroutine PMF_energy(q,ePMF,ePMF_q)
  65 c Calculate the energy and derivative in q due to the biasing PMF
  66 c Caution! Only ONE q is handled, no multi-D q-restraints available!
  67       implicit none
  68       include 'DIMENSIONS'
  69       include 'DIMENSIONS.PMF'
  70 #ifdef MPI
  71       include 'mpif.h'
  72       include 'COMMON.SETUP'
  73       integer IERROR,ERRCODE
  74 #endif
  75       include 'COMMON.IOUNITS'
  76       include 'COMMON.MD'
  77       include 'COMMON.REMD'
  78       include 'COMMON.QRESTR'
  79       include 'COMMON.PMF'
  80       integer i,iqmin,iqmax,irep
  81       double precision q,qmin,qmax,ePMF,ePMF_q
  82 c Determine the location of the q
  83       irep=1
  84       do while (dabs(t_bath-remd_t(irep)).gt.1.0d-4)
  85         irep=irep+1
  86       enddo
  87       if (irep.gt.nrep) then
  88         write (iout,*) "CHUJ NASTAPIL in PMF_energy!!!"
  89         write (iout,*)
  90      &    "Bath temperature does not match that of any replica"
  91         write (iout,*) "t_bath",t_bath
  92         write (iout,*) "remd_T",(remd_T(:nrep))
  93         write (iout,*) "Terminating..."
  94         call flush(iout)
  95 #ifdef MPI
  96         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
  97 #endif
  98         stop
  99       endif
 100       iqmin=tmin(irep,iset)
 101       iqmax=tmax(irep,iset)
 102       qmin=iqmin*delta_q
 103       qmax=iqmax*delta_q
 104 #ifdef DEBUG
 105       write (iout,*) "PMF_energy q",q," qmin",qmin," qmax",qmax,
 106      & " irep",irep," iset",iset," iqmin",iqmin," iqmax",iqmax
 107 #endif
 108       if (q.le.qmin) then
 109         ePMF_q=(PMFtab(iqmin+1,irep,iset)-PMFtab(iqmin,irep,iset))
 110      &     /delta_q
 111         ePMF=PMFtab(iqmin,irep,iset)+ePMF_q*(q-qmin)
 112 #ifdef DEBUG
 113        write (iout,*) "q<=qmin ePMF",ePMF," ePMF_q",ePMF_q
 114 #endif
 115       else if (q.ge.qmax) then
 116         ePMF_q=(PMFtab(iqmax,irep,iset)-PMFtab(iqmax-1,irep,iset))/
 117      &      delta_q
 118         ePMF=PMFtab(iqmax,irep,iset)+ePMF_q*(q-qmax)
 119 #ifdef DEBUG
 120        write (iout,*) "q>=qmax ePMF",ePMF," ePMF_q",ePMF_q
 121 #endif
 122       else
 123         do i=iqmin+1,iqmax
 124           qmax=i*delta_q
 125           if (q.ge.qmin .and. q.le.qmax) then
 126             ePMF_q=(PMFtab(i,irep,iset)-PMFtab(i-1,irep,iset))/delta_q
 127             ePMF=PMFtab(i-1,irep,iset)+ePMF_q*(q-qmin)
 128 #ifdef DEBUG
 129        write (iout,*) "qmin<q<qmax bin",i," ePMF",ePMF," ePMF_q",ePMF_q
 130 #endif
 131             exit
 132           endif
 133           qmin=qmax
 134         enddo
 135       endif
 136       return
 137       end