[gmx-users] need help with free energy of solvation calculation

[David Mobley]

Jonathan,

I responded to this four days ago but my e-mail bounced because it was
too long, so let me try again. Here's what I wrote before:

Jonathan,

I can't speak too much to the details of the system you're working on,
but there are a couple things I would try checking. I haven't read the
paper, though...


>The reported value in Table 3 of that paper is 393 kJ per mol.  I tried to
>calculate the value with Gromacs 3.3rc3 using thermodynamic integration with
>soft cores (since the hydrogens that change to methyls have no sigma or
>epsilon) and 20 lambda values between 0 and 1 (50 ps of equilibration
>followed by 150 ps of production at each lambda value).  I have reasonable
>confidence that I have the force field implemented correctly because I
>calculated exactly the same values for the various energy terms as one of
>the paper's authors did for a configuration of an isolated cellobiose
>molecule provided by the author.  However, the free energy difference that I
>calculated was the wrong sign and an order of magnitude too large, so I must
>have done something very wrong.  One thing that has me perplexed is that if
>I do a gmxdump of the .tpr file, it looks like the A and B states are
>exactly the same for things like atom type, charge, and mass...suggesting
>that I'm not specifying the B state successfully.  But if the A and B states
>are exactly the same, then why am I getting anything at all for dgdl?  I'm
>relatively new to GROMACS and have never done a coupling-parameter-type
>calculation before.

One thing to be aware of is that the minute you switch on soft-core,
the potential energy becomes a function of lambda. Hence you can
actually do no perturbation and still get nonzero dgdl at most of your
lambda values. It's only the integral that's zero. Odd, I know, but
that's how it works. I am not sure why your A and B states would be
the same, though; your toplogy file (mostly) looks OK (see below). If
you can send the full topology and coordinate files and an mdp file as
attachments I can take a more detailed look.

Also, I would do this in two steps: First change the charges, then
change the vdW. Then you can use linear (not soft core) scaling for
the electrostatics part and only change the vdW with soft core. You
can do the electrostatics part with relatively few lambda values this
way, as it's often pretty smooth. Maybe this isn't necessary here
since you shouldn't have any charge sites with zero vdW radii, though
(which is usually the reason you want to do this in two steps) so
maybe this isn't a problem.

Also (and here's where my ignorance of the paper comes in), are you
actually intending to change some of the atoms from hydrogens to
carbons, and appear new sites? Your topology doesn't look like it's
appearing any new sites right now, I don't think (don't you need dummy
atoms initially for that? I don't have any exp with those sorts of
calculations). And if you are changing atom types which involve a
change in mass you'll probably have a dEkin/dlambda term which is
nonzero also. Since the free energy difference really requires
integrating dH/dlambda you probably need to take that into account too
(or maybe you're already doing so?).

I can probably help you troubleshoot a little more if you'll send
topologies, etc. Otherwise, I would try starting with something which
has been set up in GROMACS before and make sure you can get the right
answer. I can provide a topology for something simple like
disappearing methane in water (and a right answer) if that would be
helpful.

David