[gmx-users] Free energy of turning on restraints

[Robert Johnson]

Thanks for your response David.

Are you assuming that the lambda dependence of the restraint follows
from equation 4.113? Or do you use the expression in equation 4.100
for harmonic potentials?  If the latter, I'm unsure how to apply this
equation for the dihedral restraints (equation 4.70) - again, because
of that MOD in there.

Thanks,
Bob


On Thu, Feb 28, 2008 at 10:21 AM, David Mobley <dmobley at gmail.com> wrote:
> Bob,
>
>
>  >  Unfortunately, I'm having problems including the correct restraints in
>  >  my topology file. First off, the form of the angle restraint in
>  >  Gromacs involves 4 (2 pairs) atoms. The restraint in the Boresch paper
>  >  assumes that the angle is between 3 atoms. Is there anyway to apply
>  >  the correct restraint in Gromacs?
>
>  Yes. Check out the angle calculation expression in the manual (if I
>  remember correctly) or the source code and then see what it reduces to
>  in the following case, for example:
>
>  [ angle_restraints ]
>  ; i j k l            type    theta0     fc             mult
>  1393 1391 2601 1391  1      94       fc_angle     1
>
>  Notice that I have the index 1391 repeated, which means I'm getting
>  the angle between the 1393-1391 and 2601-1391 vectors. IIRC, there is
>  something that I thought was funky about how gromacs calculates the
>  angle such that I would get back the angle I thought I wanted only if
>  I took the A-B C-B angle rather than the A-B B-C angle. But I believe
>  this is all in the manual. Let me know if you have any trouble finding
>  it.
>
>
>  >  Also, it doesn't seem like the [ angle_restraints ] or [
>  >  dihedral_restraints ] sections can be perturbed in a free energy
>  >  calculation. For example, when I include two sets of restraint
>  >  parameters (one for topology A and one for topology B) for, say,
>  >  dihedral restraints, I get the following error:
>  >  -------------------------------------------------------
>  >  Program grompp, VERSION 3.3.1
>  >  Source code file: toppush.c, line: 1180
>  >
>  >  Fatal error:
>  >  Incorrect number of parameters - found 9, expected 5 or 5.
>  >  -------------------------------------------------------
>
>  I believe they can, but I am not sure having never done this myself.
>  Actually the way I do it is without using the free energy code: I just
>  put a name like fc_angle in my topology (see above) and then use a
>  define statement in my mdp files to vary the force constant with
>  "lambda", i.e.:
>  define= -Dfc_angle=0.4184 -Dfc_dist=41.84 -Dfc_dihed=0.4184 in one mdp file.
>
>  Since you know the functional form of your restraining potential, it's
>  fairly straightforward to either (a) compute dV/dlambda using the
>  restraining potential you get back (for TI), or (b) compute the
>  potential energy at other lambda values by evaluating the angle at
>  each snapshot (i.e. for BAR).
>
>  However, this is of course less straightforward than using the free
>  energy code and I really only do it this way because when I started,
>  there was no way to do it with the free energy code. As I said, I
>  believe it can now be done with the free energy code but I am not
>  sure. Developers?
>
>  Thanks!
>  David
>
>
>
>  >  How do I perturb (i.e. turn on) these restraints over the course of a
>  >  simulation?
>  >
>  >  Thanks,
>  >  Bob Johnson
>  >  _______________________________________________
>  >  gmx-users mailing list    gmx-users at gromacs.org
>  >  http://www.gromacs.org/mailman/listinfo/gmx-users
>  >  Please search the archive at http://www.gromacs.org/search before posting!
>  >  Please don't post (un)subscribe requests to the list. Use the
>  >  www interface or send it to gmx-users-request at gromacs.org.
>  >  Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>  >
>  _______________________________________________
>  gmx-users mailing list    gmx-users at gromacs.org
>  http://www.gromacs.org/mailman/listinfo/gmx-users
>  Please search the archive at http://www.gromacs.org/search before posting!
>  Please don't post (un)subscribe requests to the list. Use the
>  www interface or send it to gmx-users-request at gromacs.org.
>  Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>