[gmx-developers] GROMACS 4.0.7 QM/MM exclusions
Lee-Ping
leeping at MIT.EDU
Fri Nov 19 20:17:55 CET 2010
Hi there,
>From our discussions, I guess a function to copy the moltype doesn't
currently exist. I decided to implement one by myself (see attached);
the only difference is that my version allocates memory for the ilist
and assigns each value using a loop. This is designed to protect the MM
bonded interactions from being modified by the generate_qmexcl_moltype
function.
I'm including some output from gmxdump (trimmed to only show the
molecule types). As you can see in the attachments, before I added the
fix, the two bonds in the MM water molecule were "0 2" and "0 2"; after
the fix, the two bonds are now "0 1" and "0 2".
I tested this by running a QM/MM (H2O/H2O) calculation and a pure-MM
calculation (2H2O) without the QM/MM interaction potential (i.e. MM
atomic charges set to zero). Before the fix, the calculations gave
different forces for the second water molecule; after the fix, they now
give the same force.
- Lee-Ping
On Fri, 2010-11-19 at 08:53 -0500, Lee-Ping Wang wrote:
> Hi there,
>
> Berk:
>
> I do use the energy group exclusion already; otherwise the VdW interaction between atoms in the QM group are counted, such that a pure-QM calculation still has a nonzero LJ energy.
>
> The two problems of the segfault, and of the unintended modification of MM bonded interactions, exists even when the energy group exclusions are used (e.g. energygrps = System, energygrp_excl = System System).
>
> Gerrit:
>
> I was under the impression that I had already put the QM atoms and the MM atoms into the same topology. In my topology I include a single itp file at the top, "spce.itp", and my topology has two water molecules in the [ molecules ] section. At the start of generate_qmexcl, there is only a single molblock (with 2 molecules) and a single moltype. After it is called, there are two molblocks with 2 moltypes, but the original moltype of the MM water was altered.
>
> Thanks a lot,
>
> - Lee-Ping
>
> On Nov 19, 2010, at 5:08 AM, Berk Hess wrote:
>
> > On 11/19/2010 08:21 AM, Gerrit Groenhof wrote:
> >> Hi,
> >>
> >> This is a problem when the QM atoms are not in the same topology. I
> >> have not yet found a satisfactory solution, also because a very easy
> >> work-around is to put all QM atoms into the same topology file.
> >>
> >> But I agree that´s not the best solution and this needs to be solved.
> >>
> >> Gerrit
> >>
> > Isn't the energy group exclusion that we discussed a simple solution?
> > We could do something like that automated as well (if it works).
> >
> > Berk
> >>
> >> On 11/19/2010 01:01 AM, Lee-Ping wrote:
> >>> ver, there is still another problem. The entire moltype struct is
> >>> not actually copied at line 1000; line 1004 copies the exclusions over
> >>> and the comment claims that "this is the only thing that needs to be
> >>> modified for QMMM." But this is untrue; generate_qmexcl_moltype not
> >>> only modifies the exclusions, but also deletes entries from the bonded
> >>> interactions (ilist) within moltype!
> >>>
> >>> This will modify the bonded interactions for the MM molecules if they
> >>> were originally the same type as the QM molecules. I have verified
> >>> this, and it's a serious problem that will break QM/MM calculations
> >>> where QM and MM molecules start from the same type (e.g. if they are
> >>> both water).
> >>>
> >>> I think this is a bug that can be avoided by copying over the ilist
> >>> struct (or the entire moltype struct) when the QM/MM moltype is created,
> >>> but no such functions seem to be built-in to GROMACS 4.0.7. I could
> >>> write these functions myself if people would find it useful. Please
> >>
> >
> > --
> > gmx-developers mailing list
> > gmx-developers at gromacs.org
> > http://lists.gromacs.org/mailman/listinfo/gmx-developers
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-developers-request at gromacs.org.
>
-------------- next part --------------
moltype (0):
name="SOL"
atoms:
atom (3):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 1.59994e+01, q= 0.00000e+00, mB= 1.59994e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 8}
atom[ 1]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 2]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom (3):
atom[0]={name="O1"}
atom[1]={name="H2"}
atom[2]={name="H3"}
type (3):
type[0]={name="O",nameB="O"}
type[1]={name="H",nameB="H"}
type[2]={name="H",nameB="H"}
residue (1):
residue[0]={name="SOL"}
cgs:
nr=3
cgs[0]={0..0}
cgs[1]={1..1}
cgs[2]={2..2}
excls:
nr=3
nra=9
excls[0][0..2]={0, 1, 2}
excls[1][3..5]={0, 1, 2}
excls[2][6..8]={0, 1, 2}
Bond:
nr: 6
iatoms:
0 type=4 (BONDS) 0 1
1 type=4 (BONDS) 0 2
Angle:
nr: 4
iatoms:
0 type=5 (ANGLES) 1 0 2
moltype (1):
name="SOL"
atoms:
atom (3):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 1.59994e+01, q= 0.00000e+00, mB= 1.59994e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 8}
atom[ 1]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 2]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom (3):
atom[0]={name="O1"}
atom[1]={name="H2"}
atom[2]={name="H3"}
type (3):
type[0]={name="O",nameB="O"}
type[1]={name="H",nameB="H"}
type[2]={name="H",nameB="H"}
residue (1):
residue[0]={name="SOL"}
cgs:
nr=3
cgs[0]={0..0}
cgs[1]={1..1}
cgs[2]={2..2}
excls:
nr=3
nra=9
excls[0][0..2]={0, 1, 2}
excls[1][3..5]={0, 1, 2}
excls[2][6..8]={0, 1, 2}
Connect Bonds:
nr: 6
iatoms:
0 type=0 (LJ_SR) 0 1
1 type=0 (LJ_SR) 0 2
-------------- next part --------------
moltype (0):
name="SOL"
atoms:
atom (3):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 1.59994e+01, q= 0.00000e+00, mB= 1.59994e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 8}
atom[ 1]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 2]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom (3):
atom[0]={name="O1"}
atom[1]={name="H2"}
atom[2]={name="H3"}
type (3):
type[0]={name="O",nameB="O"}
type[1]={name="H",nameB="H"}
type[2]={name="H",nameB="H"}
residue (1):
residue[0]={name="SOL"}
cgs:
nr=3
cgs[0]={0..0}
cgs[1]={1..1}
cgs[2]={2..2}
excls:
nr=3
nra=9
excls[0][0..2]={0, 1, 2}
excls[1][3..5]={0, 1, 2}
excls[2][6..8]={0, 1, 2}
Bond:
nr: 6
iatoms:
0 type=4 (BONDS) 0 2
1 type=4 (BONDS) 0 2
Angle:
nr: 4
iatoms:
0 type=5 (ANGLES) 1 0 2
moltype (1):
name="SOL"
atoms:
atom (3):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 1.59994e+01, q= 0.00000e+00, mB= 1.59994e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 8}
atom[ 1]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 2]={type= 1, typeB= 1, ptype= Atom, m= 1.00790e+00, q= 0.00000e+00, mB= 1.00790e+00, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom (3):
atom[0]={name="O1"}
atom[1]={name="H2"}
atom[2]={name="H3"}
type (3):
type[0]={name="O",nameB="O"}
type[1]={name="H",nameB="H"}
type[2]={name="H",nameB="H"}
residue (1):
residue[0]={name="SOL"}
cgs:
nr=3
cgs[0]={0..0}
cgs[1]={1..1}
cgs[2]={2..2}
excls:
nr=3
nra=9
excls[0][0..2]={0, 1, 2}
excls[1][3..5]={0, 1, 2}
excls[2][6..8]={0, 1, 2}
Connect Bonds:
nr: 6
iatoms:
0 type=0 (LJ_SR) 0 1
1 type=0 (LJ_SR) 0 2
-------------- next part --------------
A non-text attachment was scrubbed...
Name: topio_upload.c
Type: text/x-csrc
Size: 31073 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20101119/51bf087e/attachment.bin>
More information about the gromacs.org_gmx-developers
mailing list