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

[Justin Lemkul]

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) interactions 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!
>>
>> * 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

==================================================