From: Ignacio Fernández Galván <[EMAIL PROTECTED]>
Reply-To: Discussion list for GROMACS users <gmx-users@gromacs.org>
To: gmx-users@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.
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [EMAIL PROTECTED]
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
