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

[Justin Lemkul]

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

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