[gmx-users] Re: Some questions about OPLS/Berger combination

[chris.neale at utoronto.ca]

The flaw in your reasoning is a misunderstanding of [ pairtypes ].  
Look to section 5.3.4, especially page 99, and 5.7.1, especially page  
112, of the gromacs 4 manual. The fudge is only applied to generated  
pairs. Since we define [ pairtypes ] for all lipid and lipid-water  
pairs, fudgeLJ is not applied to those (although fudgeQQ is). When I  
originally posted this method, Berk had a nice suggestion that it  
would be simpler not to modify epsilon at all and instead to add a  
second pairs list that contained all zeroes for the LJ portion and  
thus this second pairs section would only double the coulomb part to  
restore it to full strength. While this is a cleaner solution, the  
result should be the exact same and I had already fully tested my .itp  
files so there was no motivation for me to retest this new, cleaner,  
method.

If you still disagree, please use your system to do zero-step mdruns  
(i.e. nstep=0) with i) the original ffgmx/berger combination from the  
Tieleman website, ii) my method, and iii) with whatever method you  
believe to be superior. I posit that you will find the energies of my  
method and the ffgmx/berger combination to be identical excepting  
slight differences in rounding errors due to a difference in order of  
operations, and further that you will find your suggestion to yield  
different energies. This is in fact one of the many tests that I did  
at the outset.

The remaining issue is how will the lipids see the opls protein for  
generated nonbonded interactions. This obviously has nothing to do  
with 1-4 interactions, and is a problem inherent in mixing  
forcefields. I don't propose that my method addresses this problem to  
any degree. In fact, I state quite explicitly that I don't even  
attempt to address this problem in the pdf that I have posted on my  
website. The Berger lipids are, quite simply, not parameterized  
consistently with any protein forcefied. The problem that I solved is  
only that of getting the Berger lipids to maintain their intended  
intramolecular interactions.

Chris.

Thank you Chris for your answer.
It didn't though clarify the things to me, so I'll try to explain again
where my objection is.

> The method is sound. A key realization is that the [ pairs ] section
> defines both the LJ and the Q pairs, so we will have double of each.
> We actually want double Q, since it was cut to 0.5 and 0.5*2=1.0. What
> we don't actually want is double epsilon, since the epsilon in
> question is already a special pairtypes epsilon that has been properly
> modified for 1-4 interactions. We thus cut epsilon in half to
> counteract the fact that we are forced to add it again twice in the
> pairs section.

But why should the epsilons be doubled? For OPLS-AA both the FudgeLJ and
FudgeQQ parameters are set to 0.5, so both the Coulombic and LJ 1-4
interactions are cut by 50%. We want them both to be 100% so, to
counteract this reduction, one can put a second copy of [ pairs ] section
to make GROMACS calculate both potentials twice. Dividing the epsilons by
2 will excessively reduce the 1-4 LJ potential (see below)

>> Another point concerns the 1-4 Coulombic interaction; you state that
>> dividing the Berger epsilons by 2 will modify both the 1-4 interactions
>> types, the LJ and the Coulombic ones.
>
> I don't believe that I do say that anywhere.

Ok, I misinterpreted what you wrote in
http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html
"1. Divide the epsilon entries in the [ pairtypes ] of lipid.itp by 2.0.
Now the LJ and coulombic 1-4 interactions are both exactly half of what
they should be."
>

> While I am always open to the possibility that I have made a mistake,
> I actually put a lot of time into developing this and distributing it,
> so I'm not very motivated to go look into it again before I see a bit
> of data or perhaps a proper mathematical derivation of whatever
> problem you propose. I still believe that the method is properly
> derived and implemented.
>
> Chris.
>

All right, I'll try to explain what I mean with the exact formulas.
I understand the OPLS-AA LJ potential for 1-4 interaction between atoms A
and B is calculated as:
V(OPLS) = 0.5*4*epsilon*[(sigma/r)^12 - (sigma/r)^6] = 0.5*V(Berger)
where 0.5 is the scaling factor as defined by FudgeLJ and the epsilon and
sigma values are taken from the [ pairtypes ] section, if the pair A-B has
been specified there, otherwise they are calculated from the epsilons and
sigmas for the A and B atoms by the combination rule.
So, if you put into this formula the Berger epsilon for a given pair of
atoms, as taken from the [ pairtypes ], you get V(OPLS) which is a half of
what the potential should be. What you need is to double this value, which
may be achieved by putting a second copy of [ pairs ].
Now, if you divide the Berger epsilons by 2, the calculated OPLS potential
will be
V(OPLS) = 0.5*4*(epsilon/2)*[(sigma/r)^12 - (sigma/r)^6]=0.25*V(Berger)
Putting a second copy of [ pairs ] will double this value, so you end with
50% of the proper potential value instead of 100%.
So, where is the flaw in my reasoning?

All my best,
Dorota