[gmx-users] Re: decoupling charge while maintaining intramolecular potentials

[chris.neale at utoronto.ca]

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
> clarify?

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:
>> http://www.gromacs.org/pipermail/gmx-developers/2006-January/001498.html
>> 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,
> though.
>
>> 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
>> states).
>
> 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  
increased.

Still, I believe that it is not yet possible to define coulombic  
interactions in the pairs section so this would not work for the  
coulombic portion.

Thanks you,
I sincerely appreciate your time in assisting me,
Chris.