[gmx-users] RE: Gibbs Energy Calculation and charges

Tue Nov 5 23:50:36 CET 2013

Thank you for the pointer Michael.

couple-intramol = no
What a diff of the output from gmxdump of the two tpr files shows is in both of these cases (normal and double charged), when:
	lambda is set to 1.0 (atoms within both molecules will have zero charge)
	lambda is set to 0.00 and 0.50, respectively (both will have the same charge)
There are the following differences:
	functype[] = LJC14_Q, qi and qj are set to the original charges, not the ones scaled by lambda
	functype[] = LJC_NB, qi and qj are set to the original charges, not the ones scaled by lambda
	atom[]     = q is set to the original charges, not the ones scaled by lambda

So this explains why I do not see the two topologies giving the same value at the same atomic charges, since the topologies being simulated are not the same.  The 1-4 charge interactions are still at the original charges, as are the 1-5 and beyond.

What is the reason that the charges are being left untouched here with changing lambda?  I can understand when we are dealing with LJ/van der Waals, since the 1-4 are important for the proper dihedrals, but what is the reason for charges being left untouched?  Having thought this through, answered it myself, it is because here we are interested in molecule - external environment interactions being turned off, we are moving the entire molecule from full interacting with its external environment to non-interacting.  The molecule itself should be left alone.

Turning to the side issue of turning couple-intramol on -

couple-intramol = yes
diff of the equivalent files, there are the following differences:
	functype[] = LJC14_Q, qi and qj are set to the original charges, not the ones scaled by lambda
	atom[]     = q is set to the original charges, not the ones scaled by lambda

Which confirms what Michael mentioned earlier about couple-intramol only affecting those 1-5 and beyond i.e. LJC_NB

Which then begs the question, why does the value of dH/dl change so dramatically when this option is turned on, as I observed at http://ozreef.org/stuff/gromacs/couple-intramol.png  The only thing being changed is the fact that LJC_NB is now being scaled with lambda.

Catch ya,

Dr. Dallas Warren
Drug Delivery, Disposition and Dynamics
Monash Institute of Pharmaceutical Sciences, Monash University
381 Royal Parade, Parkville VIC 3052
dallas.warren at monash.edu
+61 3 9903 9304
---------------------------------
When the only tool you own is a hammer, every problem begins to resemble a nail. 

> -----Original Message-----
> From: gmx-users-bounces at gromacs.org [mailto:gmx-users-
> bounces at gromacs.org] On Behalf Of Michael Shirts
> Sent: Thursday, 31 October 2013 1:52 PM
> To: Discussion list for GROMACS users
> Subject: Re: [gmx-users] RE: Gibbs Energy Calculation and charges
> 
> I likely won't have much time to look at it tonight, but you can see
> exactly what the option is doing to the topology.  run gmxdump on the
> tpr.  All of the stuff that couple-intramol does is in grompp, so the
> results will show up in the detailed listings of the interactions, and
> which ones have which values set for the A and B states.
> 
> On Wed, Oct 30, 2013 at 5:36 PM, Dallas Warren
> <Dallas.Warren at monash.edu> wrote:
> > Michael, thanks for taking the time to comment and have a look.
> >
> > The real issue I am having is a bit deeper into the topic than that,
> my last reply was just an observation on something else.  Will
> summarise what I have been doing etc.
> >
> > I have a molecule that are calculating the Gibbs energy of hydration
> and solvation (octanol).  In a second topology the only difference is
> that the atomic charges have been doubled.  Considering that charges
> are scaled linearly with lambda, the normal charge values of dH/dl from
> lambda 0 to 1 obtained should reproduce that of the double charged
> molecule from lambda 0.5 to 1.0.  Is that a correct interpretation?
> Since the only difference should be that charge of the atoms and over
> that range the charge will be identical.
> >
> > I was using couple-intramol = no and the following are the results
> from those simulations.
> >
> > For the OE atom within the molecule, I have plotted the following
> graphs of dH/dl versus charge of that atom for both of the topologies.
> >         octanol - http://ozreef.org/stuff/octanol.gif
> >         water - http://ozreef.org/stuff/water.gif
> >         mdp file - http://ozreef.org/stuff/gromacs/mdout.mdp
> >
> > The mismatch between the two topologies is the real issue that I am
> having.  I was hoping to get the two to overlap.
> >
> > My conclusion based on this is that there is actually something else
> being changed with the topology by GROMACS when the simulations are
> being run.  The comments in the manual allude to that, but not entirely
> sure what is going on.
> >
> >> From the manual:
> >>________________________
> >>couple-intramol:
> >>
> >>no
> >>         All intra-molecular non-bonded interactions for moleculetype
> >>couple-moltype are replaced by exclusions and explicit pair
> >>interactions. In this manner the decoupled state of the molecule
> >>corresponds to the proper vacuum state without periodicity effects.
> >>yes
> >>         The intra-molecular Van der Waals and Coulomb interactions
> are
> >>also turned on/off. This can be useful for partitioning free-energies
> >>of relatively large molecules, where the intra-molecular non-bonded
> >>interactions might lead to kinetically trapped vacuum conformations.
> >>The 1-4 pair interactions are not turned off.
> >
> > Chris Neale commented:
> >>Ah, I see. I guess that you are using couple-intramol = no (the
> default in v4.6.3 at least). That means
> >>that the intramolecular charge-charge interactions are always at
> full-strength (and therefore different).
> >>I would expect that normal at lambda=0 should be the same as double
> at lambda=0.5 only for
> >>couple-intramol = yes
> >>If you were using couple-intramol = yes already, then I am as
> confused as you are.
> >
> > >From Chris' comment, I then went back and repeated the calculation
> using couple-intramol = yes with the results of this in the below image
> (plus the previous when set to no for comparison)
> >
> > http://ozreef.org/stuff/gromacs/couple-intramol.png
> >
> > My rather garbled comment that you saw was my attempt to make sense
> of the fact that this setting made things worse (and try to add
> something to get this issue to have another pass on the list, hoping
> someone like yourself will comment), and the value of dH/dl for c-i=yes
> at lambda=1 matches with c-i=no at lambad=0.  This you can see with the
> previously linked to graphs.
> >
> > Hopefully that helps understand what I have been doing and the
> issues.
> >
> > Catch ya,
> >
> > Dr. Dallas Warren
> > Drug Delivery, Disposition and Dynamics
> > Monash Institute of Pharmaceutical Sciences, Monash University
> > 381 Royal Parade, Parkville VIC 3052
> > dallas.warren at monash.edu
> > +61 3 9903 9304
> > ---------------------------------
> > When the only tool you own is a hammer, every problem begins to
> resemble a nail.
> > --
> > gmx-users mailing list    gmx-users at gromacs.org
> > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> > * Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists