[gmx-users] Re: Problems with non-bonded interactions using OPLSAA
Berk Hess
gmx3 at hotmail.com
Tue Jul 14 12:55:06 CEST 2009
Hi,
Again, to avoid confusion:
there are two (or three) completely separate parts to this issue in Gromacs.
1) The generation of the actual pair interactions in the [ pairs ] section
in the topology (.top) file. pdb2gmx does this for you.
If you write a topology by hand, you also have to add the 1-4 interactions
in the [ pairs ] section by hand.
2) The generation of the interaction parameters for the interactions
in the [ pairs ] section. This is what gen-pairs switches and this is what
the comment in the manual is about. Pair parameters are taken from
the [ pairtypes ] section in the force field file or topology file.
If parameters are not found for a certain pair of atom types in [ pairtypes ],
they will be generated when gen-pairs is true, a fatal error will be
generated when gen-pairs is false.
(I will change "interactions" to "interaction parameters" in the quoted
line of the manual.)
3) pairs are bonded interactions and bonded interactions in no way affect
the calculation of non-bonded interactions in Gromacs. Thus to not calculate
interactions twice, you need to use nrexcl=3 when you put 1-4 interactions
in the [ pairs ] section.
Berk
Date: Tue, 14 Jul 2009 13:43:19 +0300
From: ntelll at gmail.com
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Re: Problems with non-bonded interactions using OPLSAA
Hi.
I found the same behavior with both 3 (3.3.3) and 4 (4.0.5) versions on
x86_64.
What I found is :
1. gen-pair correctly generates the 1-4 parameters (verified with
gmxdump -s topol.tpr)
But with no [pairs] defined, 1-4 interactions
are not calculated no matter what the scale factor is.
One has to define the [pairs] to calculate scaled 1-4.
2. It seems that in the case of nrexcl 2 and pairs defined, the 1-4
interactions are calculated twice
(once with scaling 1 in LJ(SR) and once scaled by whatever the
scaling factor is in LJ-14)
With nrexcl 3 and pairs defined only the scaled 1-4 are calculated.
Probably the manual is not clear at this point (p.99 of v4 manual),
where it states :
The GROMOS force fields list all these interactions explicitly, but
this section might be
empty for force fields like OPLS that calculate the 1-4 interactions by
scaling.
To verify this I did :
zero steps with a single molecule in a large box, varying nrexcl and
scaling factors with pairs defined.
What I found is :
NREXCL SCALE FACTOR LJ(SR) LJ-14
2 0.0 3.22186
0
2 0.5 3.22186
2.58086
2 1.0 3.22186
5.16172
3 0.0 -1.93987
0
3 0.5 -1.93987
2.58086
3 1.0 -1.93987
5.16172
Scaling in 1-4 works as it is expected.
The total LJ (LJ(SR) + LJ-14) in the case of nrexcl 3 and scale factor
1.0 is 3.22185 that is equal to
the LJ(SR) in the case of nrexcl 2 and scale factor 0. This means that
with nrexcl 2,
the 1-4 are calculated unscaled no matter the value of scaling factor
AND scaled by scaling factors,
i.e. scaling 1.5 instead of 0.5 as Mike notes,
Similar for Coulombic interactions.
What I use is nrexcl 3 and pairs defined.
DD
Berk Hess wrote:
Hi,
There is surely no problem with the gen-pairs option, as that would
cause
serious trouble with nearly all simulations performed with Gromacs.
gen-pairs only sets the generation of pair parameters, not of the actual
pair interactions in the topology.
I have no clue to what it going wrong in your system.
You can file a bugzilla at bugzilla.gromacs.org.
Berk
> Date: Mon, 13 Jul 2009 21:09:37 +0200
> From: mikewykes at gmail.com
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Re: Problems with non-bonded interactions
using OPLSAA
>
> Dear all
>
> I have since found the problem. For some reason gen-pairs = yes
option
> is not working, so my 1-4 interactions were being completely
excluded
> and not scaled by 0.5 as they should have been. Defining all 1-4
pairs
> in a [ pairs ] section in conjunction with nexcl =3 in the topology
> solved the problem.
>
> I was not sure whether to use nexcl =2 or nexcl =3 in the case of
> scaled 1-4 interactions, but by trial and error I found that one
> should use nexcl =3, otherwise you have the full 1-4 interaction +
the
> scaled one, so in total a scaling of 1.5 instead of 0.5.
>
> Any ideas why gen-pairs = yes option is not working would be much
appreciated.
>
> Thanks to those who responded to my first email,
>
> Mike
>
>
>
> On Thu, Jul 2, 2009 at 6:54 PM, Mike
Wykes<mikewykes at gmail.com> wrote:
> > Dear all
> >
> > I would like to simulate beta cyclodextrin in various organic
> > solvents with the OPLSAA FF for Carbohydrates (J Comput Chem
18:
> > 1955-1970, 1997) but am having problems with the short range
Coulomb
> > and LJ interactions.
> > This FF (and OPLS in general) does not assign LJ parameters to
> > hydrogen atoms in OH groups, relying on the repulsion between
oxygens
> > to keep the hydrogen (charge +0.435 ) of one OH group getting
too
> > close to the O (charge -0.7 ) of another.
> > However in MD simulations, the hydrogen of one OH group
collides with
> > the O of another, and shortly after the system explodes.
> >
> > Obviously, this could be a mistake of how I converted the
parameters
> > in the paper into gromacs parameters, so I have checked this
and found
> > no mistakes. Some of the parameters are taken from regular
OPLS,
> > allowing me to check my conversion by comparing to the
parameters in
> > gromacs/share/top/ffoplsaa*.itp files.
> >
> > Out of curiosity I implemented the same forcefield in the
tinker md
> > package and the O...H system was stable during MD, with no
O..H
> > collisions. Comparing the energies of exactly the same
geometry of
> > beta cyclodextrin with the same OPLS parameters shows
identical bonded
> > interactions, but differences in the non-bonded interactions:
> >
> > tinker(kcal/mol) gmx (kj/mol)
gmx(kcal/mol) difference
> > Total Potential energy 312 -1416.98 -338.67 650.67
> > Bond Stretching 23.32 97.58 23.32 0
> > Angle Bending 30.38 127.11 30.38 0
> > Torsional Angle 312.06 1305.64 312.06 0
> > Van der Waals 13.69 -119.03 -28.45 42.14
> > Charge-Charge -67.45 -2828.29 -675.98 608.52
> >
> > In both cases the molecule was in the gas phase, all
non-bonded
> > interactions being treated with a cutoff of 1.5 nm.
> >
> > Any suggestions as to what could be going wrong in my gromacs
> > calculations would be much appreciated.
> >
> > Please find my mdp, top and itp files below. I am using
version 4.0.5.
> >
> > Many thanks,
> >
> > Mike
> >
> >
> _______________________________________________
> gmx-users mailing list gmx-users at 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Express yourself instantly with MSN Messenger! MSN Messenger
_______________________________________________
gmx-users mailing list gmx-users at 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-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090714/c672be69/attachment.html>
More information about the gromacs.org_gmx-users
mailing list