[gmx-developers] GROMACS 4.0.7 QM/MM exclusions

[Gerrit Groenhof]

  Hi Lee-Ping,

I am sorry the fixes are that dirty, but in the end you can do what you 
want. I have looked into this several time, but have not found a 
satisfactory solution yet.

As Berk suggested, the exclusions also works.

Concering the interface to Q-chem, we'd be very interested in that. Also 
we like to have the modification you did to the non-grid based nbsearch.

Best,

Gerrit



On 11/22/2010 06:53 PM, Lee-Ping Wang wrote:
> Hi Gerrit,
>
> I agree now that including all of the molecules in a single topology will fix the moltype copying problem with generating exclusions; using the original build using --with-qmmm_gaussian now gives the correct topology.  I agree with you that the fix would still be nice, since including several water molecules in a single topology is a bit unusual. :)
>
> Also, for small systems with no MM region, I agree that there are better solutions for running MD than QM/MM in GROMACS (e.g. BOMD or CPMD).  I intended it to be a minimal test system but perhaps it was too "minimal".
>
> You are right that my fix in topio.c does not exclude the intra-QM classical interactions.  I'm still not sure how to do that in the source code.  To get the energy calculations to work, I still need to set the energy group exclusions in my .mdp file - e.g. "energygrp_excl = qm qm".  The combination of these two things seems to give the correct behavior; in my tests the intra-QM classical interactions are all zero, even when the QM region is big enough to feel intra-QM VdW and Coulomb interactions.
>
> - Lee-Ping
>
> PS: As a side note, we would be happy to contribute the gmx/Q-Chem interface to the official GROMACS after everything works, if you are interested.
>
> -----Original Message-----
> From: gmx-developers-bounces at gromacs.org [mailto:gmx-developers-bounces at gromacs.org] On Behalf Of Gerrit Groenhof
> Sent: Monday, November 22, 2010 4:39 AM
> To: Discussion list for GROMACS development
> Subject: Re: [gmx-developers] GROMACS 4.0.7 QM/MM exclusions
>
>    Hi Lee-Ping,
>
>
>
>> 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.
> What I meant is that if you have two waters, you make an top or itp like:
>
> [ atoms ]
> O  blablab
> H
> H
> O
> H
> H
>
> [ bonds]
> 1 2 5
> 1 3 5
> 4 5 5
> 4 6 5
>
> So that the  complete QM region is in one topology.
>
> If you have small QM subsystems, with no MM present, switch to
> QMMMscheme=  ONIOM.
>
> Again, not a nice solution, but if you have only QM atoms, you are
> better off running Born-Oppenheimer dynamics in the QM package.
>
> Hope this helps,
>
> I have tried your topio.c bigfix. It does now indeed generate a tpr
> file, but neither the intra-QM region LJ interactions, nor the coulomb
> interactions are excluded, and therefore double counted.
>
> I agree this needs fixing, but for the moment, the trick above might get
> you further.
>
> Gerrit
>
>
>> 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/
>>>