[gmx-users] Shell (Drude) model for polarization in GROMACS
ericsmoll at gmail.com
Wed Jul 4 00:52:50 CEST 2018
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
If I would use 2.6 in the current Gromacs thole_polarization format, what
would I use in your thole_polarization format?
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:
>> 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
>> 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
>> screening, then 1-5 non-bonded interactions and beyond can be a full
> 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
> 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
>> 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.
> (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.
>> 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
>>>> in Gromacs is an extremely long and confusing 2014 gmx-users thread
>>>> you and an extremely confused Gromacs user (
>>>> 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
>>>> we need to define 10 gromacs-particles. 5 of these gromacs-particles
>>>> 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
>>>> gromacs-atom type as discussed in (1) below.
>>>> (1) As you clarified earlier, I create a polarization directive
>>>> 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
>>> 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
>>>> a harmonic force. This interaction is not considered a "bonded
>>>> interaction" meaning that the resulting connectivity
>>>> connectivity is not considered when excluding interactions at the
>>>> specified in the molecule type directive. Only the network defined by
>>>> 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
>>>> a function of distance?
>>>> 1 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>> 2 bond away: gromacs-particle/gromacs-particle non-bonded interactions
>>>> 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
>>> 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
>>>> 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
>>>> (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
>>>> gromacs-particles, can I specify each gromacs-shell particle with the
>>>> 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
>>>> 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
>>>> 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
>>> 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 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
>>> 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
> 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