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
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
http://lists.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 gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
