[gmx-users] Re: decoupling charge while maintaining intramolecular potentials
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Wed Jan 28 06:07:43 CET 2009
Thank you Michael for the detailed reply, I have posted clarifications
and additional questions or remarks inline below.
> Hi, Chris-
> I unfortunately can't be too much help, because my free energy
> calculations are done through a modified version of Gromacs 3.1.4, and
> I am currently working with Berk and Erik to get the important
> modifications into the 4.0 branch.
>> I am a new user of the free energy code. I am somewhat confused
>> regarding the method that should be applied to decouple the long range
>> interactions of the solvent from the solute while still maintaining
>> intramolecular long-range interactions for the solute.
> By long range interactions, you mean the ones outside the cutoff,
> correct? Because there's another set of difficulties dealing with the
> ones that are shortange nonbonded interactions as well. Could you
When I said "long-range" I should have said "non-bonded" instead
throughout this entire document. Sorry for the incorrect terminology.
Therefore I am under the impression that one wants to decouple the
intermolecular solvent-solute nonbonded interactions without affecting
any intramolecular nonbonded interactions. I have seen this applied in
your papers and some others, many of them done via charmm.
>> I was able to find some information in this thread:
>> but it is still unclear to me.
> I'm not surprised, it was unclear for us at the time as well :)
>> I have reproduced the methane and tip3p energies of solvation based on
>> the tutorial that is on the gromacs wiki. In this case, I simply
>> assigned new charge values of 0 in the B state without making any
>> special considerations for intramolecular O-H values. However, tip3p
>> is rigid and perhaps maintaining the intramolecular q-q and LJ
>> components is not essential in this case.
> This is where I get confused -- are you talking about long range, or
> any intramolecule interactions. Certainly for methane and tip3p,
> there aren't any, because all atoms are 1,3 neighbors.
I am talking about any nonbonded intramolecular interactions. Further,
your point regarding 1-2 and 1-3 interactions is well taken. My
molecule of interest is a dodecylphosphocholine detergent and
therefore there will be substantial intramolecular nonbonded
interactions that I believe it would be best to maintain.
>> 1. Is it necessary to maintain the intramolecular long-range
>> interactions for the solute while decoupling LJ or charge? If not
>> absolutely required, does it affect the rate of convergence?
> If you change the intramolecular nonbonded interactions, then you
> would have to perform a second vacuum calculation in which you turn
> them back on. In terms of rates of convergence -- nobody really
> knows. If it's a intramolecular hydrogen bonding system, then turning
> off the intramolecular interactions might be faster.
What I am actually doing is annihilating one DPC detergent from a DPC
micelle and, separately, annihilating one DPC detergent monomer from
bulk water. The dG(monomer) minus dG(micelle) should then be related
to the relative probabilities of monomeric and aggregated states. I
have left out the monomer restraints and the volume correction here
for simpicity. It is my impression that this monomer-in-bulk
annihilation provides the second simulation that you mentioned above,
although I will still require this second simulation even if I
maintain my intramolecular non-bonded interactions since I am
interested in the relative probabilities of the monomeric and
aggregated state, both in solution, and not in solvation free energies.
The reason that I thought one would want to not change the
intramolecular interactions is that it seems to me that this will make
the d(pot)d(lambda) energies larger and then subtracting 200-190=10
seems like it will have more uncertainty than if the dvdl integrated
energies were all smaller. This is just based on the idea that a small
difference of two large numbers is often difficult to obtain precisely.
However, I see your point about relatively stable intramolecular
interactions. I will need to think further about that one.
>> 2. Is this already handled by the free-energy code?
> I can't speak for the 4.0 code. Berk was introducing nonbonded pair
> terms such that these pair terms would overrule the 'alchemical'
> transformation, resulting in unchanged intramolelcular nonbonded
> interactions. I actually don't know the current state of this change,
>> 3. If not, how might one go about doing this? My confusion with some
>> additional [ pairs ] entry is how gromacs would get the right
>> combination for lambda=0 and lambda=1 (not to mention intermediate
> I'm a bit confused. I would think that you would just want to set
> them to the lambda=0 state, so that intramolecular interactions were
> preserved. Am I thinking of something different than you are?
If you're just talking about how to handle the lambda=1 simulation by
setting pairs according to the nonbonded values that are achieved in
the lambda=0 state, then I believe that we are thinking about the same
thing. However, it gets messy for intermediate values of lambda (not
to mention impossible for the charging simulations) so I'll outline my
logic more elaborately below.
Assume that the 4.0.3 free-energy code uses lambda to scale all
nonbonded interactions including intreamolecular ones for the
annihilated molecule. Further assume that I create a new [ pairs ]
section that entirely accounts for all intramolecular LJ interactions
of the DPC to be annihilated (let's talk LJ annihilation only to make
this part simpler). If I include this new pairs section for lambda=1,
then I expect no solute-solvent nonbonded interactions but to maintain
regular intramolecular nonbonded interactions from the pairs section.
However, if I include this new pairs section for lambda=0.5, then I
expect to have 1.5x strength nonbonded intramolecular interactions for
the annihilated solute (once from pairs and one-half from regular
nonbonded). This would require different .itp files for each lambda
value with pairs set according to:
pair(lambda)=pair(Astate)*(1-lambda) + pair(Bstate)*(lambda)
if the conversion from A state to B state is linear and thus the new
pairs section would be different for each lambda and would compensate
for the loss of intramolecular nonbonded interactions as lambda is
Still, I believe that it is not yet possible to define coulombic
interactions in the pairs section so this would not work for the
I sincerely appreciate your time in assisting me,
More information about the gromacs.org_gmx-users