**From:** Carlo Guardiani (*carlo.guardiani_at_dsf.unica.it*)

**Date:** Mon Feb 20 2017 - 13:20:05 CST

**Next message:**Obayed Ullah: "Regarding Error on HB2"**Previous message:**Oscar Bastidas: "Charm++ Fatal Error Question"**In reply to:**Giacomo Fiorin: "Re: Metadynamics state file"**Next in thread:**Giacomo Fiorin: "Re: Metadynamics state file"**Reply:**Giacomo Fiorin: "Re: Metadynamics state file"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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.

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 ?

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
*

*>
*

**Next message:**Obayed Ullah: "Regarding Error on HB2"**Previous message:**Oscar Bastidas: "Charm++ Fatal Error Question"**In reply to:**Giacomo Fiorin: "Re: Metadynamics state file"**Next in thread:**Giacomo Fiorin: "Re: Metadynamics state file"**Reply:**Giacomo Fiorin: "Re: Metadynamics state file"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

*
This archive was generated by hypermail 2.1.6
: Sun Dec 31 2017 - 23:21:05 CST
*