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

[chris.neale at utoronto.ca]

Hi Dorota, please keep all such correspondence on the gromacs mailing  
list. One can then send somebody a personal email indicating that you  
would appreciate it if they took a look. I have copied this answer to  
the list. Please respond there.

Quoting Dorota Jamróz <jamroz at chemia.uj.edu.pl>:

> Hallo Chris,
>
> While searching for information about a possibility to combine the Berger
> lipid FF with the OPLSA-AA I came across your posts at the gmx-users forum
> and also your text on the subject you made available on the net.
> I am a GROMACS beginner and I haven't yet mastered all the nuances of the
> program, which may be the reason of my problem, but I must admit I cannot
> understand the philosophy of your method. I would be very grateful if you
> could answer my questions and clear my doubts.
>
> To begin with - I presume, the main idea of the procedure is to modify the
> non-bonded parameters from the Berger FF so that the potential energy for
> the lipid, as calculated with the OPLS-AA formulas, gives the proper value
> (i.e same as it would be if calculated with Berger FF). Am I right on this
> point?

Yes, as calculated with the OPLS-AA *1-4 rules*

> Next - I understand the main problem is the scaling of the 1-4 potential,
> which is 1 for Berger (i.e no scaling) and 0.5 for OPLS-AA (i.e. only half
> of the 1-4 potential is included into the total potential).
> If it is so then I would say, in order to get the proper 1-4 LJ potential
> you should either MULTIPLY the Berger epsilons by 2 or leave them as they
> are but insert a second copy of [ pairs ] into your topology file. Your
> solution was to DIVIDE the Berger epsilons by 2 AND include two [ pairs ]
> section. With the FudgeLJ = 0.5, as defined by OPLS-AA, the calculated LJ
> 1-4 potential will be then half of the proper value (taking half of the
> 1-4 LJ potential calculated with epsilons divided by 2 gives you 0.25 of
> the proper value. You sum these interaction twice, so finally you get the
> V(OPLS)=0.5V(Berger) )

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 don't take my word for it. Do a zero-step mdrun and look at your  
energy breakdown via g_energy for both lipid.itp/ffgmx and the hedp  
method. This is one of the things that I did to ensure that not only  
the idea was sound but there were not any unexpected bugs that would  
come to the surface with a double pairs list.

> 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.

> In what way will the Coulombic
> potential be changed? It does not depend on the epsilon value and, as I
> understand, there is no way to define another subset of charges for the
> 1-4 Coulomb interaction. The only way GROMACS allows you to modify this
> potential for specified pairs of atoms is to include it twice in the total
> potential, if a given pair is found in the [ pairs ] section.
> Summarizing, it seems to me that the 1-4 potential calculated with the
> parameters changed according to your method will give the right Coulombic
> potential but only half of the proper LJ potential.

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.

>
> As I wrote, it's pretty possible I am a victim to some misconception about
> how GROMACS works. In any case I would very grateful for you comment.
>
> Cheers,
> Dorota
>
>
> --
> Dorota Jamróz, Ph.D
> Faculty of Chemistry
> Jagiellonian University
> ul.Ingardena 3, 30-060 Cracow, POLAND
> Phone: (+48) 012 6632263
> e-mail: jamroz at chemia.uj.edu.pl
>
>
>
>
>