[gmx-users] posting an answer to my previous question on how to turn off only intermolecular vdw interactions

[David Mobley]

This post is for the record -- it's not a question, but an answer.

Previously I wrote to the list (this thread:
http://www.gromacs.org/pipermail/gmx-users/2005-November/018307.html,
and here: http://www.gromacs.org/pipermail/gmx-users/2005-November/018367.html)
having problems figuring out how to turn off only intermolecular vdW
interactions. I was told that it is possible to do this but it isn't
documented. Below will follow a brief explanation for anyone who needs
to know how to do this.

Particularly, for free energy perturbation/TI, I wanted to turn off a
small molecule's vdW interactions with a protein. However, I would
prefer not to turn off its vdW interactions entirely (partly because
in some cases this allows it to adopt really strange states and, I
think, also have funny interactions with the thermostat, especially in
the totally uncoupled state) and would rather leave its intramolecular
vdW on. That is, I wanted to turn off only its interactions with its
environment.

Anyway, it turns out that the way to do this involves several
modifications in your topology file:
1) As usual for turning off vDW interactions, define dummy atom types
with vdW sigma and epsilon of 0.
2) Specify [ nonbond_params ] between the atoms whose interaction you
want to maintain. This list is by pairs of atom types, so for every
pair of atom types whose interaction you want to maintain you need to
have an entry here specifying the vdW params for their interaction
(that is, you need to specify that the "dummy" atoms are supposed to
still continue interacting with one another). This list will get quite
long if you are using a large molecule with lots of different atom
types.
The way the parameters work, if you do NOT specify nonbonded
parameters for a pair here, the parameters are automatically generated
using whatever CombRule your force field uses for the vdw parameters.
To override this for your dummy atoms you need to explicitly state the
parameters here.
3) If your molecule uses a pairlist to add 1-4 interactions, you need
to explicitly state the interaction parameters in the pairlist,
otherwise they will often be automatically generated using the
combination rule and fudgeLJ. (Note that if you explicitly specify
parameters here you need to take fudgeLJ for your force field into
account when you compute them).

That's it. To summarize: You explicitly add back in nonbonded
interactions for any interactions you want to maintain, and you need
to make sure all of your pair interactions you want to maintain have
parameters specified in your topology file.

There is a fairly simple check you can do to make sure this is working
right: Compute the free energy of turning off vdW interactions for
your small molecule in vacuum. If you are correctly maintaining the
intramolecular interactions, the free energy of doing this should be
zero.

There is, however, a caveat: If you are using soft core, the free
energy of doing this will only be zero in total. That is, if you are
doing TI, you will have to integrate dV/dlambda to get zero (and if
you're using FEP/BAR you'll have to sum up your energy differences).
The reason for this is that the soft core method sort of changes the
effective distances between molecules, so soft core changes the
potential energy with lambda even when you aren't actually doing
anything.

I hope if anyone else needs to know how to do this, this will be sufficient.

David Mobley