Metadynamics state file

From: Carlo Guardiani (carlo.guardiani_at_dsf.unica.it)
Date: Fri Feb 17 2017 - 06:53:31 CST

Dear NAMD experts,

I a novice in running metadynamics simulations
with NAMD and I would like to ask some clarifications
on the content of the state file. My well-tempered
metadynamics simulation is biasing two collective
variables.

1) As I understand, the first list of numbers
   represents the biasing potential. I checked
   that the number of values in this list
   corresponds to the number of grid ponts in the
   pmf files. Could you please confirm whether these
   numbers represent the values of the biasing
   potential (sum of gaussians) in correspondence
   of the grid points ?

2) If this is the case, in which order are these
   grid point values listed ? Is it the same order
   employed in the pmf file, i.e.

   [CV1(1) CV2(1)] [CV1(1) CV2(2)] [CV1(1) CV2(3)]
          . . .
          . . .
          . . .

   [CV1(1) CV2(n2-2)] [CV1(1) CV2(n2-1)] [CV1(1) CV2(n2)]

          . . .
          . . .
          . . .

   [CV1(n1) CV2(1)] [CV1(n1) CV2(2)] [CV1(n1) CV2(3)]
          . . .
          . . .
          . . .

   [CV1(n1) CV2(n2-2)] [CV1(n1) CV2(n2-1)] [CV1(n1) CV2(n2)]

3) What is the exact meaning of the second list of numbers ?
   Do they represent the values of the gradients with respect
   to the two CV in corresponcence of the grid points ? In which
   order are they listed ? Is it the following (calling the two
   CV x1 and x2) ?

   dV[x1(1) x2(1)]/dx1 dV[x1(1) x2(1)]/dx2 dV[x1(1) x2(2)]/dx1
   dV[x1(1) x2(2)]/dx2 dV[x1(1) x2(3)]/dx1 dV[x1(1) x2(3)]/dx2
          . . .
          . . .
          . . .

   dV[x1(1) x2(n2-2)]/dx1 dV[x1(1) x2(n2-2)]/dx2 dV[x1(1) x2(n2-1)]/dx1
   dV[x1(1) x2(n2-1)]/dx2 dV[x1(1) x2(n2)]/dx1 dV[x1(1) x2(n2)]/dx2

          . . .
          . . .
          . . .

   dV[x1(n1) x2(1)]/dx1 dV[x1(n1) x2(1)]/dx2 dV[x1(n1) x2(2)]/dx1
   dV[x1(n1) x2(2)]/dx2 dV[x1(n1) x2(3)]/dx1 dV[x1(n1) x2(3)]/dx2
          . . .
          . . .
          . . .

   dV[x1(n1) x2(n2-2)]/dx1 dV[x1(n1) x2(n2-2)]/dx2 dV[x1(n1) x2(n2-1)]/dx1
   dV[x1(n1) x2(n2-1)]/dx2 dV[x1(n1) x2(n2)]/dx1 dV[x1(n1) x2(n2)]/dx2

4) How does NAMD compute the gradient ? Does it use analytical or numerical
   derivatives ? And in this latter case, how does it deal with the points at
   the boundary of the grid ?

5) Since I am using the keepHills option, I tried to recover the biasing
   potential by summing the guassians whose parameters are listed at the end
   of the state file. The PMF that I retrieved is qualitatively similar to
   the one computed by NAMD and shown on the pmf file but my free energy
values
   are significantly higher (it's not just the matter of an additive
constant,
   the distance between the deepest minimum and the highest maximum is
larger).
   The core of the fortran90 program that I wrote is the following:

      fact=(Temp + DeltaT)/DeltaT
      open(10,file=inpfile,status='old',form='formatted')
11 read (10,*,end=12) istep, weight, x1c, x2c, sig1, sig2

      do i=1,n1
       do j=1, n2
          x1=x1min+(i-1)*dx1
          x2=x2min+(j-1)*dx2
          exp1=(x1-x1c)**2
          exp1=exp1/(2.d0*sig1**2)
          exp2=(x2-x2c)**2
          exp2=exp2/(2.d0*sig2**2)
          exptot=exp1+exp2
          hills(i,j)=hills(i,j)+ weight*dexp(-exptot)
       end do
      end do

      goto 11
12 close(10)

      open(20,file=outfile1,status='unknown',form='formatted')
      do i=1,n1
       do j=1, n2
          x1=x1min+(i-1)*dx1
          x2=x2min+(j-1)*dx2
          pmf(i,j)=-hills(i,j)*fact
          write(20,*) x1, x2, hills(i,j)
       end do
       write(20,*) " "
      end do
      close(20)

      Could you please tell me what I did wrong ? Does NAMD employ
      any special trick to compute the biasing potential ?

      Thank you very much for your help and kind regards,

      Carlo Guardiani

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:06 CST