[gmx-developers] GROMACS 4.0.7 QM/MM exclusions

[Gerrit Groenhof]

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/