[gmx-users] Free energy, frozen atoms and exclusions

[Berk Hess]

>From: Ignacio Fernández Galván <jellby at yahoo.com>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: gmx-users at gromacs.org
>Subject: [gmx-users] Free energy, frozen atoms and exclusions
>Date: Mon, 11 Sep 2006 15:18:01 +0100 (BST)
>
>Hi all,
>
>I'm trying to calculate solvation free energies and I've found
>something "strange".
>
>The thing is I'm interested in free energy of rigid molecules (I
>consider other contributions separately), typically a small solute
>molecule in water. I'm making a test with methanol. The setup is a
>frozen methanol molecule in a box of of water molecules. Since I'm
>considering rigid molecules, I create an .itp file with no bonded
>terms, just atoms (with charges and LJ). The methanol molecule is
>called MOH, and I have this in the .mdp:
>
>energy_grps          = MOH
>energygrp_excl       = MOH MOH
>freezegrps           = MOH
>freezedim            = y y y
>
>which should freeze the methanol and disregard all nonbonded terms
>between the methanol atoms (and there are no bonded terms).
>
>Well, after the usual equilibration and such, I get this output for the
>initial configuration with lambda=1:
>
>LJ (SR): 5.70341e+03
>Coulomb (SR): -3.65377e+04
>...
>dVpot/dlambda: 2.81161e+02
>
>It seems right, but after simulation for several lambda values I get a
>free energy of around 10 times what I would expect.
>
>Now I repeat exactly the same simulation, but adding all pairs of
>methanol atoms in the [exclusions] of the .itp file. For the same
>initial configuration as before, I get this:
>
>LJ (SR): 5.70341e+03
>Coulomb (SR): -3.65377e+04
>...
>dVpot/dlambda: 2.14385e+00
>
>(the only difference is in dVpot/dlambda). This (if I extrapolate)
>would give a much more sensible value for the free energy.
>
>As far as I could understand, having "energygrp_excl = MOH MOH" should
>be the same as putting all atoms in the [exclusions] section, but it
>obviously isn't, at least when it comes to calculating dVpot/dlambda.
>
>Is there something I'm missing? A bug/feature in the code? A
>limitation?

If you are using PME, the mesh contribution is not excluded when
using energygrp_excl, grompp should give you a warning about this.
Adding all the pairs as exclusions solves this problem.

>
>Thanks in advance
>
>PS. As a side question, I'd say that, if I have a frozen molecule it
>would be best *not* to remove the center of mass motion (not even of
>the rest of the system). Am I right?

Yes.
Note that there are still issues with freeze-groups and pressure scaling
as the rest of the system scales around your frozen atoms.
There is currently no solution for this,
although slow coupling would help.

But why do you want to freeze the methonal?
Things would be much easier with a free methanol.

Berk.