[gmx-users] Shell (Drude) model for polarization in GROMACS
Justin Lemkul
jalemkul at vt.edu
Sun Jul 1 21:28:14 CEST 2018
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
==================================================
More information about the gromacs.org_gmx-users
mailing list