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

[Eric Smoll]

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 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
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
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....

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
(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.
>