source/unres/src-HCD-5D/thermo_length.f

   1       subroutine thermo_length
   2 c Integrating |PMF| to claculate thermodynamic length
   3       implicit none
   4       include 'DIMENSIONS.PMF'
   5       include 'COMMON.IOUNITS'
   6       include 'COMMON.MD'
   7       include 'COMMON.PMF'
   8       double precision aL(0:maxHdim,maxT_h)
   9       double precision aLdelta,qdelta,qtemp,
  10       integer nT,tmin,tmax,i,j,t,tlow,tup,nW
  11       qdelta=(qmax-qmin)/(nW-1)
  12       do j=1,nT
  13         do t=tmin,tmax
  14           aL(t,j)=0.0d0
  15         enddo
  16         dH=hfin(tmin+1,j)-hfin(tmin,j)
  17         print *,"tmin",tmin,hfin(tmin+1,j),hfin(tmin,j)," dH",dH
  18         do t=tmin-1,0,-1
  19           hfin(t,j)=hfin(t+1,j)-dH
  20         enddo
  21         aL(0,j)=0.5d0*(dabs(hfin(0,j))+dabs(hfin(1,j)))
  22         do t=1,tmax-1
  23           aL(t,j)=aL(t-1,j)+0.5d0*(dabs(hfin(t,j))+dabs(hfin(t+1,j)))
  24         enddo
  25         aL(tmax,j)=aL(tmax-1,j)
  26         do t=0,tmax
  27          aL(t,j)=(aL(t,j)-0.5d0*(dabs(hfin(0,j))+dabs(hfin(1,j))))*delta
  28         enddo
  29         write (*,*) "T=",1.0d0/(1.987D-3*beta_h(j))," PMF and L"
  30         do t=0,tmax
  31           write (*,'(f5.2,2f10.5)') t*delta,hfin(t,j),aL(t,j)
  32         enddo
  33 c Determining window centers
  34         tlow=(qmin+1.0d-4)/delta
  35         tup=(qmax+1.0d-4)/delta
  36         aLspan=aL(tup,j)-aL(tlow,j)
  37         aLdelta=aLspan/(nW-1)
  38         print *,"tlow",tlow," tup",tup," aLspan",aLspan,
  39      &   " aLdelta",aLdelta
  40         q0(1,j)=qmin
  41         q0(nW,j)=qmax
  42         do i=2,nW
  43           aLw=aLdelta*(i-1)
  44           print *,"i",i," aLw",aLw
  45           do t=tlow,tup-1
  46             if (aLw.ge.aL(t,j).and.aLw.le.aL(t+1,j)) then
  47               print *,"i",i," t",t," aLw",aLw,aL(t,j),aL(t+1,j)
  48               q0(i,j)=t*delta+delta*(aLw-aL(t,j))/(aL(t+1,j)-aL(t,j))
  49               exit
  50             endif
  51           enddo
  52         enddo
  53         write (*,*) "q0"
  54         write (*,'(20f10.5)') (q0(i,j),i=1,nW)
  55         return
  56         end
  57 c Force constants
  58         do t=1,tmax-1
  59           Wprim=(hfin(t+1,j)-hfin(t-1,j))/(2*delta)
  60           Wbis=(hfin(t+1,j)+hfin(t-1,j)-2*hfin(t,j))/delta**2
  61           keff=0.5d0*Wprim**2+dabs(Wbis)
  62      &    +dsqrt((0.5d0*Wprim**2+dabs(Wbis))**2-Wbis**2)
  63           write (*,'(f5.2,3f10.1)') t*delta,Wprim,Wbis,keff
  64         enddo
  65       enddo ! j
  66       end
  67