[gmx-users] Shell (Drude) model for polarization in GROMACS

[Eric Smoll]

Hello Justin,

Thank you for the guidance.

Although more testing is needed, molecules where shells are attached to
every atom appear to be working properly.

The issue appears to be with models where shells are only attached to heavy
atoms and bonds between hydrogen atoms and heavy atoms are constrained.
Initial energy minimization tests (emtol = 100) show that bonds between
heavy atoms and hydrogen in organic molecules slowly stretch or compress
far more than they should while the network of non-hydrogen atoms maintains
a sensible geometry.  All atom-drude displacements on the heavy atoms
converge quickly and are stable.

I am still hunting for whatever topology problem is causing this.  All
exclusions out to 1-4 interactions between atoms and drudes should be
properly included (combined action of the moleculetype directive and
additional exclusions directives).  Thole screening only applies to
atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H
atoms (no attached shells).  It may be difficult to provide any useful
guidance without details but troubleshooting suggestions are welcome if you
have any.

Are there issues associated with adding shells to specific atoms in a
molecule?  Are simulations that place shells on all atoms (hydrogen and
heavy) more stable for some reason?  I am going to build a new model where
shells are attached to all atoms to see if bonds to hydrogen atoms still
slowly compress/stretch during a tight-emtol energy minimization.

Best,
Eric

On Mon, Jul 9, 2018 at 6:28 AM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 7/9/18 1:31 AM, Eric Smoll wrote:
>
>> Thanks Justin,
>>
>> That is very helpful.  I can run the water example provided in your paper
>> with the GROMACS "drude" branch.
>>
>> I have built my own tools to construct a core-shell MD topology file for
>> my
>> system.  Grompp produces a tpr without complaint but mdrun crashes after
>> the first step with a segmentation fault and no other useful errors.  I am
>> trying to troubleshoot.
>>
>> Is it required that each drude/shell immediately follow the atom it is
>> attached to in the atoms directive and the input coordinate file?  For
>> convenience, I have attached them to the end of each molecule?
>>
>
> There is no such requirement (you'll note that the water Drude is at the
> end of the atom listing in SWM4-NDP, to avoid changes in the SETTLE code).
>
> In my topology, I manually account for the fact that drude/shell particles
>> are attached with bonds which count when evaluating exclusions.
>> Specifically, I assume a "3" in the molecules directive will take care of
>> all exclusions up to 3 bonds away.  To account for the extra bond lengths
>> introduced by drude/shell particles, I add manually add certain exclusions
>> for particles that should be 3 bonds away according to the atom-atom
>> connectivity.  Will this work or is there some "gotcha" I don't
>> understand.  For example, do you invalidate the moleculetype directive by
>> providing an exclusions directive?
>>
>
> Sounds right, but without an actual example, I can't tell you anything
> really useful.
>
> -Justin
>
>
>
>> Best,
>> Eric
>>
>> On Tue, Jul 3, 2018 at 5:39 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>> On 7/3/18 6:52 PM, Eric Smoll wrote:
>>>
>>> Justin,
>>>>
>>>> Thanks again for the response and the guidance.  Your 2015 JCC paper on
>>>> the
>>>> subject was very helpful.
>>>>
>>>> Your code uses a different thole_polarization format.  Currently,
>>>> Gromacs
>>>> uses the following format
>>>> [ thole_polarization ]
>>>> ;   i   j   k   l   function-number  thole-a-parameter   alpha_ij
>>>>   alpha_kl
>>>>
>>>> Your code seems to use
>>>> [ thole_polarization ]
>>>> ;   i   j   k   l   function-number   alpha_ij   alpha_kl
>>>> thole-a-parameter-1   thole-a-parameter-2
>>>>
>>>> If I would use 2.6 in the current Gromacs thole_polarization format,
>>>> what
>>>> would I use in your thole_polarization format?
>>>>
>>>> 1.3 for both atoms to get 1.3 + 1.3 = 2.6.  Thole's magic number quickly
>>> breaks down for complicated molecules, so our convention is that it is a
>>> per-atom value that is summed to give a new value of "a" per pair, e.g.
>>> a_i
>>> + a_j = a_total.
>>>
>>> -Justin
>>>
>>> --
>>> ==================================================
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Assistant Professor
>>> Virginia Tech Department of Biochemistry
>>>
>>> 303 Engel Hall
>>> 340 West Campus Dr.
>>> Blacksburg, VA 24061
>>>
>>> jalemkul at vt.edu | (540) 231-3129
>>> http://www.thelemkullab.com
>>>
>>> ==================================================
>>>
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at http://www.gromacs.org/Support
>>> /Mailing_Lists/GMX-Users_List before posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-request at gromacs.org.
>>>
>>>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Assistant Professor
> Virginia Tech Department of Biochemistry
>
> 303 Engel Hall
> 340 West Campus Dr.
> Blacksburg, VA 24061
>
> jalemkul at vt.edu | (540) 231-3129
> http://www.thelemkullab.com
>
> ==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support
> /Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>