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.
_______________________________________________
gmx-users mailing list [email protected]
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
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
