[gmx-developers] GROMACS 4.0.7 QM/MM exclusions

Mon Nov 22 00:13:51 CET 2010

Hi Gerrit,

Thanks.

You mentioned that putting everything into a single .itp file will solve
the problem, but I think I'm already doing this.  My test system
consists of one QM water and one MM water, and their force fields are
both "spce.itp".  I can actually put everything into a single topol.top
file (attached), and I still encounter the segmentation fault.

To reproduce the bug, extract the tar, compile GMX 4.0.7 with
--with-qmmm_gaussian and run "grompp -n" in the directory.  The
index.ndx file is to distinguish the QM region.

On a different note, I think the neighbor-searching code (ns.c) might
not be working for extremely small QM/MM systems (e.g. for a single
diatomic molecule with no MM region).  In these cases, one of the atoms
gets shift vector 22, but the other atom gets shift vector 0; the QM
calculation will then crash because the atoms are several box vectors
apart.  I think this is happening at line 407 in ns.c, where a new
i-atom is assigned only if the current i-atom accumulates a nonzero
number of j-atoms.  For very small systems, an i-atom might have zero
j-atoms, so the i-atom gets overwritten; I don't think this is the
correct behavior for QM/MM.  Can somebody help look into this?

Thank you,

- Lee-Ping

On Sun, 2010-11-21 at 09:06 +0100, Gerrit Groenhof wrote:
> Hi Lee-Ping,
> 
> Thanks for you effort. Tomorrow, I will also test your code.
> 
> Coming back to my first reply: putting everything into a single .itp file does the trick as well. So if we have waters in the QM region, we add them to the topology file in which also the rest of the QM region is defined.
> 
> best,
> 
> Gerrit
> 
> On 19 Nov 2010, at 20:17, Lee-Ping wrote:
> 
> > 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.
> >> 
> > 
> > <moltype_after.txt><moltype_before.txt><topio_upload.c>-- 
> > 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.
> 
> --
> Gerrit Groenhof
> MPI biophysical chemistry
> Goettingen
> Germany
> http://wwwuser.gwdg.de/~ggroenh/
> 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: water-dimer.tar.gz
Type: application/x-compressed-tar
Size: 824 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20101121/d20b8087/attachment.bin>