From: Giacomo Fiorin (giacomo.fiorin_at_gmail.com)
Date: Thu Feb 23 2017 - 03:09:17 CST
Hi Carlo,
On Mon, Feb 20, 2017 at 2:20 PM, Carlo Guardiani <
carlo.guardiani_at_dsf.unica.it> wrote:
> Dear Giacomo,
> thanks for your help. Now I succeeded in reproducing
> the biasing potential and its gradients starting from
> the hills listed by the keepHills option. I simply hadn't
> realized that the widths reported by keepHills do not
> correspond to the standard deviations but twice the standard
> deviation.
>
Good, in newer versions you should also see that the standard deviation is
reported in the log file, labeled as half-width or "sigma".
>
> As I'm at it, I would like to ask another question about
> the hillwidth option. According to the manual, using the
> default value of sqrt(2*pi)/2 the volume of a gaussian
> should correspond to W*Prod(wi). However, there is
> something I don't understand in this calculation. If we
> consider 1D-gaussians the integral of W*exp[-(x-xc)^2/2*sigma^2]
> should be equal to W*sigma*sqrt(2*pi). If sigma=[sqrt(2*pi)/2]*wi/2
> then substituting in the expression of the integral, we get
> W*wi*pi/2 instead of W*wi. What did I get wrong ?
>
Indeed you are correct and that statement in the documentation is wrong:
thanks for pointing it out!
As you noticed, the value of "sigma" used will be reported during the
simulation in the units of the collective variables, and that can be used
to report the simulation details.
The current code design uses a single keyword "width" that all methods
use. However, its meaning is clearer in the case of ABF or histogram (the
spacing of the grid) than in metadynamics, where the grid spacing shouldn't
change the physical results and instead "sigma" is the meaningful
parameter. This is probably one case where uniformity of syntax between
methods generated a bit of confusion.
In an upcoming version of the code, the parameter sigma for metadynamics
will be given explicitly and will not have a default value. The current
syntax will still be honored to continue previous simulations.
Giacomo
>
> Thanks again for your help and kind regards,
>
> Carlo
>
>
>
>
> > 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
> >
>
>
>
-- 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 : Mon Dec 31 2018 - 23:20:07 CST