From: Giacomo Fiorin (giacomo.fiorin_at_gmail.com)
Date: Fri Feb 17 2017 - 12:22:05 CST
Hi Carlo, you can refer to the Colvars reference doc here:
http://colvars.github.io/colvars-refman-namd/colvars-refman-namd.html
All derivatives (or forces, if you change the sign) of the Gaussian
functions are computed analytically at each grid point, and those from
different Gaussians are summed together. No numerical derivatives are used
in any case. When the total force of the Gaussian bias is needed at a
point that is not exactly the grid point, the closest grid point is taken.
Not knowing your configuration nor what do you mean by "significantly
higher" I would say that you may have a large grid spacing compared to the
Gaussian width. The default option is set so that the integral of one
Gaussian equals the volume of one grid bin, which is not a fine
interpolation. It works for many small Gaussians:
http://dx.doi.org/10.1080/00268976.2013.813594
but will not be accurate for few, large Gaussians.
I suggest you try changing the width and let the Colvars code recompute the
sum of all Gaussians on the grid, and compare them to your own
interpolation.
Giacomo
On Fri, Feb 17, 2017 at 7:53 AM, Carlo Guardiani <
carlo.guardiani_at_dsf.unica.it> wrote:
> 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
>
>
>
>
>
>
>
-- Giacomo Fiorin Associate Professor of Research, Temple University, Philadelphia, PA Contractor, National Institutes of Health, Bethesda, MD http://goo.gl/Q3TBQU https://github.com/giacomofiorin
This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:04 CST