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 [email protected] 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 [email protected]. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
