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

[Dorota Jamróz]

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