Re: Metadynamics state file

From: Carlo Guardiani (carlo.guardiani_at_dsf.unica.it)
Date: Mon Feb 20 2017 - 13:20:05 CST

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
>

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