[gmx-users] Free energy of turning on restraints

[David Mobley]

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
>