[lammps-users] About the PMF in SMD

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I have a question about the PMF calculation in SMD. With “fix smd”, it can give us the accumulated PMF (the sum of pulling forces times displacement). If I use this command working on one rigid group in my NVT system with constant steering velocity, what’s the relation between the output PMF and the PMF of the system people usually talk about?

I read Park’s paper (J. Chem. Phys. 2004) mentioned in the manual, and found that it essentially need the exponential average estimator of <exp(-W/KBT)> or the cumulant estimator to get the estimate of PMF(W is the work done on the system). Also, they discussed the error of the PMF calculation based on different number of trajectories and different sized samplings. Then what about here in “fix smd”? Or there is any approximation or something else?

Your kind instruction would be most appreciated.

Best regards
Yajie

2

I'll let Axel answer this, as fix smd is his code.

Steve

3

I'll let Axel answer this, as fix smd is his code.

Steve

> Dear LAMMPS users,
>
> I have a question about the PMF calculation in SMD. With "fix smd", it can
> give us the accumulated PMF (the sum of pulling forces times displacement).
> If I use this command working on one rigid group in my NVT system with
> constant steering velocity, what's the relation between the output PMF and
> the PMF of the system people usually talk about?

here is what is accumulated into pmf:

pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;

so it is the integral over time of the force projected
on the pulling/pushing direction. this is done in every
step so a sophisticated numerical integration is not needed.

nevertheless, it is provided primarily for convenience
and debugging. you have access to all kinds of properties
computed during the SMD step and print them to a file and
post-process them as you wish.

> I read Park's paper (J. Chem. Phys. 2004) mentioned in the manual, and found
> that it essentially need the exponential average estimator of <exp(-W/KBT)>
> or the cumulant estimator to get the estimate of PMF(W is the work done on
> the system). Also, they discussed the error of the PMF calculation based on
> different number of trajectories and different sized samplings. Then what
> about here in "fix smd"? Or there is any approximation or something else?

fix smd simply implements the pushing/pulling. the sampling
issues are independent from how you compute the PMF.

of course the pulling/pushing speed matters as well as
the spring constant and how many (uncorrelated!!!) steering
calculations you do.

thus, for any trivial reaction coordinate, umbrella integration
via fix spring might be a competitive alternative, as it is
easy to adapt the sampling for each bin.

i've been working on cleaning up and simplifying the code
last fall, but then could not complete the work since. with
a bit of luck, i might get to it this summer.

in the mean time, just think of it as being implemented pretty
much straight from the paper and don't get confused by r0.
it has to be set to 0.0 for practically anything that i have
come across and will be eliminated in the next version of
fix smd, since it causes more confusion than anything else.

cheers,
   axel.

4

Hi, Axel,

Can I ask you something to clarify your email:
pmf += (fxxn + fyyn + fz*zn) * v_smd * update->dt;

Based on your explanation here, the xn is the normal unit vector in x direction, right? Not the pushing distance in the x direction?

Thank you so much,
-Fei

5

Hi, Axel,

Can I ask you something to clarify your email:
pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;

Based on your explanation here, the xn is the normal unit vector in x
direction, right? Not the pushing distance in the x direction?

the answer depends on whether you operate in tether or couple mode.

for tether, xn, yn, and zn are the components of the normalized vector
of the pushing or pulling direction. they are set to 0.0 this dimension
is ignored in SMD. those components are computed once at the beginning
and kept constant thoughout.

for couple, the pushing/pulling direction is recomputed in every step
and so are the vector components.

as a consequence, the term (fx*xn + fy*yn + fz*zn) is the spring force
projected on the current pushing/pulling direction. this is true in both
cases and the reason for the choice.

hope that helps,
    axel.

p.s.: sorry for the late reply. just noticed that i was posting from
the wrong e-mail account. please send mailing list related stuff to
[email protected] and use my temple e-mail only for personal mail.

6

Hi, Axel,

Thank you for your explanation.

But I still don’t quite understand your comment

“for couple, the pushing/pulling direction is recomputed in every step
and so are the vector components.”

If I set

fix stretch cterm smd cvel 20.0 0.0001 couple nterm 0 0 -1 0.0

which should menat the couple is done in the z direction only. Is the pushing/pulling force still recomputed in every step for its direction?

Thank you so much.

Best,
Fei

7

Hi, Axel,

Thank you for your explanation.

But I still don't quite understand your comment

"for couple, the pushing/pulling direction is recomputed in every step
and so are the vector components."

If I set

fix stretch cterm smd cvel 20.0 0.0001 couple nterm 0 0 -1 0.0

which should menat the couple is done in the z direction only. Is the
pushing/pulling force still recomputed in every step for its
direction?

this is a special case and quite untypical for
the use of smd in couple mode. with 0.0 0.0 -1.0
you pull the two groups apart, but in a way that
they are aligned in z-direction. also, the force
is - of course - always recomputed in every step,
what xn,yn,zn refers to is the orientation of
the pushing vector to which the fix group(s) are
attached.

the normalized vector it is indeed recomputed all
the time, but because of the fixed alignment of
the pushing vector, the result won't change.

a push in only z-direction with x and y being free
would be implemented through couple with a pushing
vector or NULL NULL -1.0

hope that clears up matters,
    axel.