[gmx-users] A method to scale Coulombic 1-4 interactions seperately
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Wed Sep 13 05:08:45 CEST 2006
There is a recent publication about combining the Berger lipids with OPLS-AA and
getting the correct 1-4 scaling factor. New dihedrals were developed because it
is possible to set special 1-4 parameters for LJ, but not coulomb. However,
there is another workaround to this that should be generally applicable to 1-4
coulombic scaling.
The lipids.itp file has already defined [ pairtypes ] so the LJ 1-4 is scaled
correctly (FudgeLJ will not be applied to the values in [ pairtypes ]). The
problem now is that OPLS-AA requires a FudgeQQ of 0.5 and the Berger lipids
desire a FudgeQQ of one.
The solution comes in two parts:
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.
2. Every entry in the [ pairs ] section of pope.itp should appear twice. Now the
LJ and coulombic 1-4 interactions are both exactly what they should be.
This solution should work for an arbitrary 1-4 scaling problem, just make sure
that you get the ratios correct. A detailed method is outlined below and that is
followed by the method that I used to test this procedure. I have also run a
0.5ns simulation and the system "looks" fine (when I say looks fine I mean by
visual inspection or by g_density). 0.5ns is not nearly long enough to properly
evaluate area per lipid, but for what it's worth I get 0.53nm^2 per lipid with
regular combination rules and 0.52nm^2 per lipid with the
half-epsilon-double-pair method outlined here.
Perhaps the Gromacs crew could tell us whether a double [ pairs ] entry might
run into problems in any special situations?
#####################################################
Method to change your files:
1. I assume that OPLS-AA + Berger lipid users have done something similar to this:
http://www.gromacs.org/pipermail/gmx-users/2006-May/021416.html
2. copy ffoplsaa.itp to ffoplsaa_mod.itp and edit the new file to contain
entries like this:
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 0.5 0.5
#include "/dir/ffoplsaanb_mod.itp"
#include "/dir/ffoplsaabon.itp"
3. copy ffoplsaanb.itp to ffoplsaanb_mod.itp and edit the new file to contain
entries like this:
;Here are the unmodified atomtypes from lipid.itp
LO LO 1 15.9994 0.000 A 2.96000e-01 8.87864e-01
;carbonyl O, OPLS
LOM LOM 1 15.9994 0.000 A 2.96000e-01 8.87864e-01
;carbonyl O, OPLS
;...
;Here are the original parameters from lipid.itp that have been moved directly
to this file. Note that these parameters are commented out and could be removed.
;[ pairtypes ]
; i j funct sigma epsilon
; LO LO 1 2.96E-01 1.10E-01
; LO LOM 1 2.96E-01 1.10E-01
;....
;Here are the parameters from above where epsilon has been divided by 2.0
[ pairtypes ]
; i j funct sigma epsilon
LO LO 1 2.96E-01 5.50E-02
LO LOM 1 2.96E-01 5.50E-02
;....
4. Copy your lipid topology file (e.g. pope.itp) and, in the new version, copy
the [ pairs ] section to right underneath itself (without the [ pairs ] header)
If it looked like this first (pope.itp):
[ pairs ]
4 7 1
5 8 1
6 9 1
It should look like this after (pope_mod.itp):
[ pairs ]
;Here is the first copy
4 7 1
5 8 1
6 9 1
; Here is the second copy. Only use this with the modified LJ-14 epsilon values
4 7 1
5 8 1
6 9 1
5. Your run .top file should look like this:
; Include forcefield parameters
#include "/dir/ffoplsaa_mod.itp"
; Include topologies
#include "/dir/pope_mod.itp"
[ system ]
...
#####################################################
TESTING:
Use a system with a protein inserted into a membrane and solvated. OPLS-AA rules.
1. Do a run with nsteps=0 for the original setup and for the new setup
2. gmxdump the .trr file
3. diff the two gmxdumps
The only differences are in the forces section and they only apply to atoms in
the [ pairs ] section that have been intentionally modified. For example, in
pope.itp I get different forces for all atoms except:
H1,H2,H3, C17, C20-C31, C36, C39-CA2
There is no difference for some of these because they are not in the [ pairs ]
list at all:
H1,H2,H3,C20-C23,C25-C31,C39-CA2
There is no difference for some others because one of the parteners in the pair
has a charge of zero:
C17 in [ pair ] 13 17 (LH1 LP2; LP2 has zero charge)
C36 in [ pair ] 32 36 (LC2 LP2; LP2 has zero charge)
C24 in [ pair ] 24 27 (LP2 LH1; LP2 has zero charge)
C25 in [ pair ] 22 25 (LH1 LP2; LP2 has zero charge)
Remember that the LJepsilon was divided by 2 so with double inclusion the Lj
portion is not changed at all. This is a decent internal test also because if we
had accidentally modified the effective LJ1-4 parameters [(regular/2.0 and
include twice) vs. (regular)] then we should have seen force differences for
these last four atoms (C17,C36,C24,C25).
#####################################################
Here is the .mdp from the test outlined above:
title = seriousMD
cpp = /home/cneale/exe/cpp
integrator = md
nsteps = 0
tinit = 0
dt = 0.002
comm_mode = linear
nstcomm = 1
comm_grps = System
nstxout = 1
nstvout = 1
nstfout = 1
nstlog = 1
nstlist = 10
nstenergy = 1
nstxtcout = 1
ns_type = grid
pbc = xyz
coulombtype = PME
rcoulomb = 0.9
fourierspacing = 0.12
pme_order = 4
vdwtype = cut-off
rvdw_switch = 0
rvdw = 1.4
rlist = 0.9
DispCorr = no
Pcoupl = Berensden
pcoupltype = semiisotropic
compressibility = 4.5e-5 4.5e-5
ref_p = 1. 1.
tau_p = 4.0 4.0
tcoupl = Berendsen
tc_grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300. 300.
annealing = no
gen_vel = yes
gen_temp = 300.
gen_seed = 9896
constraints = all-bonds
constraint_algorithm= lincs
lincs-iter = 1
lincs-order = 4
More information about the gromacs.org_gmx-users
mailing list