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

Eric Smoll ericsmoll at gmail.com
Wed Jul 4 00:52:50 CEST 2018 

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?

Best,
Eric

On Sun, Jul 1, 2018 at 6:52 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 7/1/18 6:37 PM, Eric Smoll wrote:
>
>> Justin,
>>
>> Thank you for the response!  I deeply appreciate the effort you put into
>> this community.
>>
>> I am still very confused *but* I learn more with every email.  I am now
>> certain that I need to add drude coordinates in the gro file.  This was
>> not
>> specified anywhere in the manual.  I have also learned that one should not
>> completley remove coulomb interactions between nearby intramolecular
>> core-shell pairs. Instead, I should leave them in with thole screening.
>> You have also clarified that this is only necessary at short range but how
>> short is not particularly clear.  For an OPLSAA style forcefield, I assume
>> it is reasonable to start by thole screening 1-2 and 1-3 coulomb
>> interactions, then scaling 1-4 non-bonded interactions by 0.5 without
>> thole
>> screening, then 1-5 non-bonded interactions and beyond can be a full
>> strength.
>>
>
> You can develop a model with whatever convention you like. The advice I'm
> giving you is what we do in our Drude force field, which works quite well.
>
> You mention that you don't really know how the polarization directive works
>> since you don't use it.  Can I implement drude polarizability with thole
>> screening without the polarization directive using your method of just
>> defining the core-shell pairs and manually attaching them with bonds?  I
>>
>
> Yes.
>
> like your method since it is less opaque.  Do I still need to label the
>> shell particles in each pair with an 'S' in the topology (also with zero
>> mass and LJ params)?  I suspect the shell label is still necessary to
>>
>
> Yes.
>
>> instruct gromacs to iteratively optimize shell positions at every
>> timestep.
>> Your method also causes complications since there is no harmonic bond type
>> that does not create a connectivity link used to auto-generate exclusions.
>> Which brings me to my next cluster of questions....
>>
>
> In this case, you *do* want to use my Drude branch, because I coded
> pdb2gmx to do precisely that. Please see the .rtp entries in the tarball
> from http://mackerell.umaryland.edu/charmm_drude_ff.shtml
>
> If you specify in the header of the .rtp file that the force field is
> polarizable, and specify alpha and Thole for each atom, pdb2gmx will do the
> rest, with all proper pairs and exclusions. Perhaps before diving into
> this, do something really easy, like ethane or something, for which you
> could compute the energy even by hand so you know what the result should be.
>
> Is there a clean strategy that will prevent me from making a mistake on the
>> intramolecular nonbonded interactions when I am using your method (no
>> polarization directive)?  I am struggling to understand what this strategy
>> is.  What needs to be excluded and how?  There are many pieces here
>>
>
> In our method, the Drudes inherit the same pair and exclusion list as
> their parent atom. The Drude-atom bond does not count for the purposes of
> determining 1-2, 1-3, and 1-4 interactions.
>
> -Justin
>
>
> (defaults directive settings, exclusion directive settings, the necessity
>> to scale 1-4 all LJ and coulomb interactions, the pairs directive, the
>> moleculetype directive, and the thole _polarization directive).
>>
>> Sorry for the additional questions.  I suspect this conversation will drag
>> on for a long time.
>>
>> Best,
>> Eric
>>
>> On Sun, Jul 1, 2018 at 1:27 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>> On 7/1/18 3:13 PM, Eric Smoll wrote:
>>>
>>> Thank you,
>>>>
>>>> I have compiled the "Drude" branch and am preparing my topology.
>>>> Unfortunately, the best reference for implementing shell molecular
>>>> dynamics
>>>> in Gromacs is an extremely long and confusing 2014 gmx-users thread
>>>> between
>>>> you and an extremely confused Gromacs user (
>>>> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users
>>>> /2014-July/091066.html).
>>>> I  have more questions.
>>>>
>>>> I am interested in adding polarization to an OPLS-AA style model. Let's
>>>> build a model system to add clarity to our communication.  Consider a
>>>> linear, 5 physical-atom (not gromacs-atom) artificial molecule.  To
>>>> describe this artificial molecule in a shell molecular dynamics
>>>> simulation
>>>> we need to define 10 gromacs-particles.  5 of these gromacs-particles
>>>> will
>>>> be of the gromacs-atom type (marked with an A in the topology) and 5 of
>>>> these gromacs-particle will be of the gromacs-shell type (marked with
>>>> an S
>>>> in the topology).  Each gromacs-shell will be "attached" to a
>>>> corresponding
>>>> gromacs-atom type as discussed in (1) below.
>>>>
>>>> (1) As you clarified earlier, I create a polarization directive
>>>> specifying
>>>>
>>>> I will begin my remarks by saying I do not do any simulations that use a
>>> [polarization] directive. In our Drude model, everything is just an atom
>>> defined under [atoms] with explicit bonds listed in [bonds]. So I may not
>>> have exact answers to everything you want to know. Further, you don't
>>> necessarily need my Drude branch to deal with anything you've mentioned
>>> here.
>>>
>>> gromacs-atom/gromacs-shell pairs.  For each entry in the polarization
>>>
>>>> directive, I specify the function number as 1 and the polarizability in
>>>> units of nm^3.  This "attaches" the gromacs-shell to the gromacs-atom
>>>> with
>>>> a harmonic force.  This interaction is not considered a "bonded
>>>> interaction" meaning that the resulting connectivity
>>>> gromacs-shell/gromacs-atom
>>>> connectivity is not considered when excluding interactions at the
>>>> distance
>>>> specified in the molecule type directive.  Only the network defined by
>>>> the
>>>> gromacs-atom/gromacs-atom connectivity is considered when excluding
>>>> interactions.  Assuming the exclusion distance is three bonds and I have
>>>> scaled 1-4 interactions specified (defaults directive settings and pairs
>>>> directive entries), what happens with the gromacs-particle interactions
>>>> as
>>>> a function of distance?
>>>> 1 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>> excluded
>>>> 2 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>> excluded
>>>>
>>>> I will note here that indeed nrexcl *should* work this way (e.g. the
>>> Drude/shell inherits its parent atom's exclusions) but it's actually
>>> important in polarizable models to *include* all of these neighboring
>>> interactions, and this is what Thole screening is for. It's really
>>> important for obtaining a physically realistic molecular polarizability
>>> and
>>> is one of the significant advantages over additive models.
>>>
>>> 3 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>
>>>> included but scaled by 0.5. I am correct in assuming that
>>>> gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and
>>>> gromacs-shell/gromacs-shell coulomb interactions are all scaled down by
>>>> a
>>>> factor of 0.5?
>>>>
>>>> Presumably, but do check with an easily definable single-point energy.
>>> Our Drude model does work this way, but I am not 100% sure what a
>>> declaration via [polarization] does. I think that's how it works, too.
>>>
>>> 4 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>
>>>> included (gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and
>>>> gromacs-shell/gromacs-shell)
>>>>
>>>> (2) Do I need to create .gro coordinates for the 10 gromacs-particles or
>>>> just the 5 gromacs-atom?  If I need to declare the coordinates for all
>>>> 10
>>>> gromacs-particles, can I specify each gromacs-shell particle with the
>>>> same
>>>> coordinates as the gromacs-atom particle that it is attached to?
>>>>
>>>> Yes, you need the coordinates to be present. Initial assignment to the
>>> coordinates of the parent atom is fine. The dipoles will adjust as a
>>> function of the electric field in the system.
>>>
>>> (3) Now, according to the manual, thole polarization is intended to
>>>
>>>> screen/damp/diminish coulomb interactions in shell molecular dynamics
>>>> simulations.  The wording suggests that thole polarization only impact
>>>> gromacs-shell/gromacs-shell interactions (not gromacs-atom/gromacs-atom
>>>> or
>>>> gromacs-shell/gromacs-shell interactions).  Am I correct?  Each entry in
>>>>
>>>> No. Thole screening is applied between neighboring dipoles, not
>>> neighboring charges. Explicit charge-charge interactions are computed for
>>> all possible atom-Drude, atom-atom, and Drude-Drude pairs in the listed
>>> interaction and are then screened via the specified constant.
>>>
>>> the associated directive requires listing 4 numbers. The first two
>>> numbers
>>>
>>>> specify a gromacs-atom and its attached gromacs-shell.  The last two
>>>> numbers specify a different gromacs-atom and its attached gromacs-shell.
>>>> Given my assumption in (1) above, I would only need to include these for
>>>> all 1-4, 1-5, etc. intramolecular gromacs-atom/gromacs-shell pairs.  I
>>>> assume that 1-4 interactions are still diminished by 0.5.  Am I correct?
>>>>
>>>> See above. You should likely not do anything special to 1-4 and 1-5
>>> interactions. Closer (by bonds) interact
>>> <https://maps.google.com/?q=interactions.+Closer+(by+bonds)+interact&entry=gmail&source=g>ions
>>> should be screened. But you
>>> can develop whatever convention you like, as long as you can justify it.
>>> But the Thole effect really is a short-range effect.
>>>
>>> -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
>>> <https://maps.google.com/?q=iling_Lists/GMX-Users_List+before+posting&entry=gmail&source=g>
>>> !
>>>
>>> * 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.
>

More information about the gromacs.org_gmx-users mailing list