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