[gmx-developers] QMMM, exclusions and neighbor search

hess at sbc.su.se hess at sbc.su.se
Mon Mar 8 12:47:38 CET 2010 

Hi,

Gerrit, are you sure?
Did 3.3 really not add exclusions between qm atoms in different molecules?

Anyhow, with PME this still would not be correct.
The PME exclusion forces are determined as erf(beta*r)/r,
which is accurate at short distance, but quite wrong at a distance
of half the box length.
Only doing an extra mesh calculation would fix this.
I am currently considering if that would be worth implementing.

Berk

> Hi Pia,
>
> It has indeed always been like this. Therefore, the fix of incuding the
> QM atoms in the energy group exclusions will have the same effect as
> explicitly removing the interactions.
>
>
> Gerrit
>
>
> Pia Toelle wrote:
>> Hi Berk.
>>
>> It just came to my mind, that there are already some QM-QM interactions
>> excluded in the code. (My first remark in July 2009 was, that ONLY the
>> nonbonded interaction of QM atoms belonging to the SAME molecule
>> (Gromacs 4.x.x).)
>>
>> So what about the "original" Gromacs 4.0.5 code and the use of QM/MM and
>> PME? And the Gromacs 3.x.x code where QM-nonbonded interaction has also
>>  been excluded as far as I remember.
>>
>> Pia
>>
>>
>> Pia Toelle wrote:
>>
>>> Hi.
>>>
>>> So what's about my idea just to set the QM-QM interaction in
>>> "sys->ffparams" (in topio.c) to zero by including special QM-atomtypes.
>>> Do you see any problem there? For me it seems to work.
>>>
>>> Pia
>>>
>>>
>>> Berk Hess wrote:
>>>
>>>> Hi,
>>>>
>>>> With PME you can not easily exclude Coulomb interactions between atom
>>>> pairs that come from
>>>> the PME mesh part. You would need one or multiple (I don't recall now)
>>>> extra mesh calculations
>>>> to do that. It can be done, but requires some coding effort.
>>>>
>>>> Berk
>>>>
>>>> Pia Toelle wrote:
>>>>
>>>>> Hello.
>>>>>
>>>>>
>>>>> Actually, we are using PME (and our QM-calculation also uses pbc with
>>>>> Ewald summation). What problem do you refer to?
>>>>>
>>>>>
>>>>> Best wishes,
>>>>> Pia
>>>>>
>>>>>
>>>>> Berk Hess wrote:
>>>>>
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I suggested a simple solution for this problem to Gerrit, which does
>>>>>> not
>>>>>> require code changes,
>>>>>> except for a check in grompp.
>>>>>>
>>>>>> I think you can simply use energy group exclusions by putting all
>>>>>> the qm
>>>>>> atoms in a separate
>>>>>> energy group and using: energygrp_excl = qm qm
>>>>>> This should do what you want, except when using PME and to a certain
>>>>>> extent also reaction field.
>>>>>> But for those methods it is quite complicated to come up with a
>>>>>> proper
>>>>>> scheme.
>>>>>>
>>>>>> Berk
>>>>>>
>>>>>> Pia Toelle wrote:
>>>>>>
>>>>>>
>>>>>>> Hi Gerrit.
>>>>>>>
>>>>>>> I made the following changes in  src/kernel/topio.c:
>>>>>>>
>>>>>>>
>>>>>>> - "generate_qmexcl": "molblock" filling changed (I am not sure, if
>>>>>>> this
>>>>>>> was related to the QM-QM interaction problem)
>>>>>>>
>>>>>>>
>>>>>>> - added "generate_qm_atom_type": new QM-atomtype for all atomtypes
>>>>>>> added
>>>>>>>  ( sys->atomtypes.nr were doubled)
>>>>>>> The added QM-atomtypes have zeros in sys->ffparams for the QM-QM
>>>>>>> non-bonded interaction.
>>>>>>>
>>>>>>>
>>>>>>> I hope, that the changes in the attached file topio.c are
>>>>>>> understandable
>>>>>>> and that it will help you.
>>>>>>>
>>>>>>>
>>>>>>> Greetings,
>>>>>>> Pia
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Gerrit Groenhof wrote:
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> Hallo Pia,
>>>>>>>>
>>>>>>>> Sorry for not getting back on this issue earlier.
>>>>>>>>
>>>>>>>> Have you been able to fix this bug?
>>>>>>>>
>>>>>>>> BEst wishes,
>>>>>>>>
>>>>>>>> Gerrit
>>>>>>>>
>>>>>>>>
>>>>>>>> On 11 Aug 2009, at 13:10, Pia Toelle wrote:
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>> Okay.
>>>>>>>>>
>>>>>>>>> Thanks for your reply. I will try to do so.
>>>>>>>>> Please, tell me if their are any news concerning this problem.
>>>>>>>>>
>>>>>>>>> Cheers,
>>>>>>>>> Pia
>>>>>>>>>
>>>>>>>>> Berk Hess wrote:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> Hi,
>>>>>>>>>>
>>>>>>>>>> That would work.
>>>>>>>>>> But we should certainly fix the code.
>>>>>>>>>> I did not think of this situation when introducing the moltype
>>>>>>>>>> data
>>>>>>>>>> structure in version 4.0.
>>>>>>>>>>
>>>>>>>>>> Berk
>>>>>>>>>>
>>>>>>>>>> Pia Toelle wrote:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> Hello.
>>>>>>>>>>>
>>>>>>>>>>> To avoid a different handling of QM-atoms of different molecule
>>>>>>>>>>> types,
>>>>>>>>>>> would it be possible to put all QM atoms in one molblock and
>>>>>>>>>>> moltype, as
>>>>>>>>>>> the molblocks ans moltypes are filled and defined anyway in
>>>>>>>>>>> generate_qmexcl? Or would this influence other parts of the
>>>>>>>>>>> program?
>>>>>>>>>>>
>>>>>>>>>>> greetings
>>>>>>>>>>> Pia
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Gerrit Groenhof wrote:
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>> Hi,
>>>>>>>>>>>>
>>>>>>>>>>>> I had not realized your previous question had not been
>>>>>>>>>>>> answered
>>>>>>>>>>>> completely, sorry.
>>>>>>>>>>>>
>>>>>>>>>>>> In 3.3 it was a bit more intuitive, also for me.
>>>>>>>>>>>>
>>>>>>>>>>>> In the generate_qmexcl there is a loop over all molblocks.
>>>>>>>>>>>> Thus all qm
>>>>>>>>>>>> atoms are found, even if they are in different molecules.
>>>>>>>>>>>> However, you
>>>>>>>>>>>> might have a point that between molecules the qm atoms do not
>>>>>>>>>>>> seem to
>>>>>>>>>>>> get excluded. Such situation could happen, for example if part
>>>>>>>>>>>> of the
>>>>>>>>>>>> protein and a few waters are treated at the QM level, which is
>>>>>>>>>>>> typical
>>>>>>>>>>>> for enzymes. To avoid problems (if any), you can include such
>>>>>>>>>>>> waters in
>>>>>>>>>>>> the .itp file of the protein.  But it should work also with
>>>>>>>>>>>> the
>>>>>>>>>>>> different QM atom in different topologies and I will  look
>>>>>>>>>>>> into this .
>>>>>>>>>>>>
>>>>>>>>>>>> Best,
>>>>>>>>>>>>
>>>>>>>>>>>> Gerrit
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> On 10 Aug 2009, at 18:27, Pia Toelle wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>> Hello.
>>>>>>>>>>>>>
>>>>>>>>>>>>> I still have problems with "generate_qmexcl" in version
>>>>>>>>>>>>> 4.0.5.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Comparing to earlier version 3.3.3, one notices that the data
>>>>>>>>>>>>> structure
>>>>>>>>>>>>> moltype and molblock are included. In the earlier version, qm
>>>>>>>>>>>>> exclusions
>>>>>>>>>>>>> were generated for each QM-atom inside the system, now in
>>>>>>>>>>>>> version 4.0.5,
>>>>>>>>>>>>> they are generated per atom inside each molecule type.
>>>>>>>>>>>>>
>>>>>>>>>>>>> I still wounder whether the exclusions between QM-atoms of
>>>>>>>>>>>>> different
>>>>>>>>>>>>> molecule types are considered here and - if yes - how. Or -
>>>>>>>>>>>>> if not -
>>>>>>>>>>>>> where this is done.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Is the splitting of the QM part into molblock's correct?
>>>>>>>>>>>>>
>>>>>>>>>>>>> I would be glad for any comment or reply! Please tell me if
>>>>>>>>>>>>> my question
>>>>>>>>>>>>> is not precise.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Pia
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> Pia Toelle wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>> Hello Gerrit.
>>>>>>>>>>>>>> Thanks for the reply, it helped a lot, but I still have some
>>>>>>>>>>>>>> problems to
>>>>>>>>>>>>>> understand.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> As far as I understood in "generate_qmexcl" the QM molecules
>>>>>>>>>>>>>> are put
>>>>>>>>>>>>>> into separate molblock and each molecule has its own
>>>>>>>>>>>>>> moleculetype.
>>>>>>>>>>>>>> Then "generate_qmexcl_moltype" is called and the exclusions
>>>>>>>>>>>>>> are set.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> At first an array "qm_arr" is generated. As the loop
>>>>>>>>>>>>>> "for(i=0;i<molt->atoms.nr;i++)" only passes through the
>>>>>>>>>>>>>> atoms with the
>>>>>>>>>>>>>> same moleculetype.
>>>>>>>>>>>>>> "sys->moltype[--].excls" contains the exclusions between
>>>>>>>>>>>>>> Atoms of the
>>>>>>>>>>>>>> same topologie.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Are the nonbonded interactions listed in the ilist and
>>>>>>>>>>>>>> excluded by
>>>>>>>>>>>>>> deleting them in this list?
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Where is the ilist generated in first place?
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Please correct my estimates.
>>>>>>>>>>>>>> Cheers,
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Pia
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Gerrit Groenhof wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> The QM - QM interactions, both Coulomb and van der Waals
>>>>>>>>>>>>>>> are explicitly
>>>>>>>>>>>>>>> exluded.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> The code for this begins with the comment
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> /* creating the exclusion block for the QM atoms. Each QM
>>>>>>>>>>>>>>> atom has
>>>>>>>>>>>>>>>  * as excluded elements all the other QMatoms (and itself).
>>>>>>>>>>>>>>>  */
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> In generate_qmexcl_moltype (topio.c).
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Thus the exclusionlists are expanded with all QM- QM
>>>>>>>>>>>>>>> interactions.
>>>>>>>>>>>>>>> The coulomb interactions are handled a bit more dirty: In
>>>>>>>>>>>>>>> the
>>>>>>>>>>>>>>> initialisation (qmmm.c), these charges are put to zero.
>>>>>>>>>>>>>>> Thus the standard nbsearching will look at the exclusions
>>>>>>>>>>>>>>> and check
>>>>>>>>>>>>>>> if a
>>>>>>>>>>>>>>> particle is charged. This way it will never include the QM
>>>>>>>>>>>>>>> QM  atom
>>>>>>>>>>>>>>> pairs, or QM MM coulomb interactions.
>>>>>>>>>>>>>>> The QMMM nbsearching in the next step takes the QM atoms as
>>>>>>>>>>>>>>> i-particles
>>>>>>>>>>>>>>> and searches all MM neighbours.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Best wishes,
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Gerrit
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> On 15 Jul 2009, at 20:10, Pia Toelle wrote:
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Hello
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> (I am using Gromacs 4.0.5 (QM/MM))
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> I have a question about the exclusion of nonbonded
>>>>>>>>>>>>>>>> (specially VdW)
>>>>>>>>>>>>>>>> interaction in QM/MM simulation. The coulomb and VdW
>>>>>>>>>>>>>>>> interaction
>>>>>>>>>>>>>>>> between
>>>>>>>>>>>>>>>> two QM Atoms should not be calculated by Gromacs as the
>>>>>>>>>>>>>>>> QM-QM
>>>>>>>>>>>>>>>> interaction is done by the QM-program.
>>>>>>>>>>>>>>>> Where (and how) in the code is the VdW (and also coulomb)
>>>>>>>>>>>>>>>> interaction
>>>>>>>>>>>>>>>> between two QM atoms excluded?
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Please also comment on the following:
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> (topio.c called by grompp.c)
>>>>>>>>>>>>>>>> "generate_qmexcl" and "generate_qmexcl_moltype", here the
>>>>>>>>>>>>>>>> exclusions for
>>>>>>>>>>>>>>>> the INTRA-molecular interaction are set. Is this right?
>>>>>>>>>>>>>>>> What is about
>>>>>>>>>>>>>>>> the Interaction between different QM-Molecules?
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> (ns.c)
>>>>>>>>>>>>>>>> In "search_neighbours" the "nsgrid_core" is called to set
>>>>>>>>>>>>>>>> up the
>>>>>>>>>>>>>>>> neighbor lists (for all atoms QM and MM !?), then it is
>>>>>>>>>>>>>>>> called to
>>>>>>>>>>>>>>>> set up
>>>>>>>>>>>>>>>> the QM/MM neighbor list.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> (force.c)
>>>>>>>>>>>>>>>> In "do_forces" the "do_force_lowlevel" is called, there
>>>>>>>>>>>>>>>> the
>>>>>>>>>>>>>>>> QM-calculation and "do_nonbonded" takes place. The later
>>>>>>>>>>>>>>>> does not
>>>>>>>>>>>>>>>> contain any information about the atom being QM or not, as
>>>>>>>>>>>>>>>> far as I
>>>>>>>>>>>>>>>> can
>>>>>>>>>>>>>>>> see. Therefore, I believe, that the decision to calculate
>>>>>>>>>>>>>>>> VdW was made
>>>>>>>>>>>>>>>> by setting up the lists. Is this right???
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Thanks in advance,
>>>>>>>>>>>>>>>> Pia Tölle
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> BCCMS, university of Bremen, Germany
>>>>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>>>>> 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/
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>>>> 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
>>>>>>>>>>>>>>> thewww
>>>>>>>>>>>>>>> interface or send it to gmx-developers-request at gromacs.org.
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>>> 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.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>>> 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/
>>>>>>>>>>>>
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> 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
>>>>>>>>>>>> thewww
>>>>>>>>>>>> interface or send it to gmx-developers-request at gromacs.org.
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> 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.
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>>> 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.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> 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/
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>
>>>>>
>>>
>>
>>
>>
>
> --
> 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.
>

More information about the gromacs.org_gmx-developers mailing list