Hi,
What I want to do is basically to get the free energy map of a linkage bond, i.e. F(\phi, \psi). According to my understanding from literature, I need to increment the \phi, and \psi and obtain a local energy minimum. Then I don’t know who to apply these procedures in LAMMPS. I checked “fix smd” which seems to be applying a force upon a bond while what I need to vary is the local torsion. I also don’t know how to output the local energy.
Any suggestion will be highly appreciated!
-Fei
I'm not clear on what you are asking to do, but you can change
the dihedral coefficients of a particular dihedral (assuming it has
its own type), by re-using the dihedral_coeff command. It could
be put in a loop: run, change, run, change, etc. There is no simple
way to get an output of the energy of a single dihedral interaction,
other than to post-process it.
Steve
Hi,
hi,
What I want to do is basically to get the free energy map of a linkage bond,
i.e. F(\phi, \psi). According to my understanding from literature, I need to
increment the \phi, and \psi and obtain a local energy minimum. Then I don't
there are multiple ways to compute free energy profiles.
what you describe source a lot like umbrella sampling.
there is also umbrella integration, steered MD, metadynamics
ABF and others that i have forgotten.
know who to apply these procedures in LAMMPS. I checked "fix smd" which
seems to be applying a force upon a bond while what I need to vary is the
local torsion. I also don't know how to output the local energy.
you don't need a "local" energy (what would that be?). if you
are after a free energy, then all contributions of all degrees of
freedom matter, since they contribute to the entropy part of your
free energy. otherwise you could just plot the potential energy
of your force field components...
indeed fix smd pushes or pulls a system into a specific direction.
our group has implemented a free energy module that has many different
free energy methods integrated, but i could not find the time, yet, to
interface it to LAMMPS. it is currently only integrated into NAMD.
there also is our competition PLUMED,
http://merlino.mi.infn.it/plumed
it is interfaced to LAMMPS, but i have never tried it.
Any suggestion will be highly appreciated!
doing free energy calculations successfully requires a decent
understanding of statistical mechanics. it is usually fairly easy to
set those free energy calculations up, but can be difficult to tell
when they are converged and if you have a correct result.
it is strongly recommended to first do a few practice runs
on simple systems with known results.
cheers,
axel.
Hi, Axel & Steve,
Thank you for the reply.
If I want to do the umbrella sampling, can I just use the way Steve suggested (vary the dihedral coefficient), then do the run + change loops?
Also, what does it mean by the difficulty to judge whether the system has converged? Shouldn’t I just wait till the energy goes to a equilibrium and then record the avg of the equilibrium trajectory?
-Fei
By thinking of it again, I don’t think I should vary the dihedral coeff. If so, the free energy will be calculated using a different coeff, and it won’t be the same as if I sample the dihedral space. How do I vary the dihedral angles then?
-Fei
Hi, Axel & Steve,
Thank you for the reply.
If I want to do the umbrella sampling, can I just use the way Steve
suggested (vary the dihedral coefficient), then do the run + change loops?
yes, almost.
Also, what does it mean by the difficulty to judge whether the system has
converged? Shouldn't I just wait till the energy goes to a equilibrium and
then record the avg of the equilibrium trajectory?
in umbrella sampling you put a harmonic restraint (or multiple) on your
system and then you record and average over the _restraint_ force until
it is converged. you do this for equilibrium angles. if you want a 2d-map,
you do this for all combination of angle pairs.
with the histogram of the averaged forces you can recover
the potential of mean force by (numerical) integration.
you have to keep in mind that this restraint goes on top of the
regular interaction potentials and has to be fairly "stiff", so that
your system stays very close to the desired reference point.
depending on the shape of your PMF, you may need more or
less bins for sampling forces and more or less trajectory frames
to get a well enough converged average force.
cheers,
axel.
To change the dihedral angle, you’d have to change the atom positions. You can do this using the set command:
http://lammps.sandia.gov/doc/set.html
Paul
If I change the dihedral angle by changing the 4 atom positions, I am also changing the relative positions of the these 4 atoms compared with other parts of the molecule. I doubt that is a right way to go.
That is the right way to go. You want to include all the effects of the unrestrained motion of every other atom in the system. This is how the entropy makes it's way into the free energy. I don't know that you even need to do umbrella sampling or other biasing technique. If the "reaction coordinate" you are looking at (the dihedral angle) is well sampled, you can get the free energy by binning the populations along a regular dynamical trajectory.
Matt
Quoting Fei Liang <[email protected]...>:
no, that is exactly the right thing to do.
in a typical umbrella sampling scenario, you would
change your restraint by a given delta, run until the
system is equilibrated, sample the restraint force
until converged and then move on.
you can turn this also into an embarrassingly
parallel problem, by first running an MD where
you slowly move the restraint around, and take
snapshot configurations at the appropriate
frames and then use those to start the initial
MD simulations for the various bins, which
are now independent and can run in parallel
first until they are well equilibrated and then
with force sampling until you have well enough
converged forces.
there should be quite a bit of literature on the
subject and also some more detailed explanations
on MD simulation text books.
cheers,
axel.