[gmx-users] hints for core/shell optimization?

[Justin Lemkul]

On 7/30/14, 5:14 AM, Tamas Karpati wrote:
> Dear Justin,
>
>> No, they're not mutually exclusive.  The Thole function screens nearby
>> dipoles so they don't experience artificial forces.  The
>> "thole_polarization" name is a bit unfortunate; it should really be
>> "thole_screening" instead.
>
> I was probably confused by the two contribs in the MDRUN log file.
> Thank you for the clarification.
>
>> I suspect you're seeing large forces due to one (or more effects).  If your
>> model doesn't have any bonds, then you don't have any excluded interactions,
>> so instead of getting a screened Coulombic interaction via Thole screening,
>> you're getting both added together.  You should make use of [exclusions] for
>> neighboring pairs, likely, but I guess the bigger question is how your force
>> field was parametrized and what its expected behavior should be.
>
> Do you mean that for the Thole screening to work I need bonds
> in order to have 1-2- and 1-3 interactions? Or do I need the exclusions
> that pull bonds in? Should I declare bonds of an average (eg. 0.2 nm)
> length with k=0 force constants then add [exclusions] for farther atoms?
> Please confirm what I need to exclude (the neighboring, 1-2, pairs rather
> than the 1-3 or farther pairs?)
>
> FYI: I took the model and parameters (fitted to experimental
> geometries and diffusion enthalpies) from a (non-GROMACS) paper
> on some crystal. The model seems chaotic to me but is perfect
> to learn MM -coming from the QC world.
>

I'm hesitant to start making suggestions about modifications without knowing the 
exact details of the model.  There are often aspects left unexplained in such 
papers, and I suspect there are special aspects of whatever software the authors 
used that control these physical properties.  If they didn't use any screening 
function explicitly, then their code must be already handling it somehow.

But in short, no, bonds are not compulsory for using Thole, but exclusions 
definitely are.  I would exclude at least the nearest-neighbor interactions.

Can you post a snippet of your [polarization] and [thole_polarization] directives?

>> The other effect is one that I noticed this weekend.  I think the forces
>> coming out of the Thole function are incorrect; at least they are in the
>> case of our Drude polarizable force field.  I changed the force calculation
>> and everything works as expected now, in my development version of the code.
>
> Nice to hear that.  I'm looking forward to seeing it at your site for
> download :)
>
>> field parameters?  In our Drude model, 1-2 and 1-3 interactions (based on
>> atom connectivity, not Drude connectivity) are excluded and electrostatics
>> calculated via Thole.  You don't have any bonds, so the choice here is not
>> exactly clear and I don't want to pre-judge without knowing what the force
>> field is or how it should be expected to behave.
>
> As I said my FF is built from literature and I consider it chaotic.
> The crystal is held together by nonbonding interactions and polarization
> forces. In addition, I do also have some angle forces (without any
> bonds). The core/shell model means just some extra shell particles
> on selected atom types; c/s -s are "bonded" by [polarization] items.
> Recently -with your aid- I started screening them by the Thole method.
> (Of course, I did otherwise but I will correct my method, thanks :)
>
> No dummies are applied. Although I though about it, I haven't found
> in the Manual a way to put D on a single A (rather than on 2 or more atoms).
> In such a way I could have reduced charges on the atoms so that
> q_D = -q_S and [polarization] would have D-S instead of A-S "bonds".
> Only I couldn't figure out how to fix D on A sites.
>

Just a note for clarity here, since it might be a bit confusing.  A Drude and a 
shell are the same thing; "dummies" (virtual sites) haven't been mentioned yet 
and don't sound like they're relevant.  A-S is an atom-shell(Drude) bond and D-S 
is a dummy(virtual site)-shell bond.  The use of S and D to mean different 
things when they could be the same in other nomenclature systems is somewhat 
confusing.

> Is "Drude connectivity" something I also use (no bonds)?
> In a sense I then have both Drude and Buckingham connectivities.
>

What I meant was this (hopefully this shows up correctly)

D1   D2   D3
|    |    |
A1---A2---A3

The exclusions and Thole screening are based on atom connectivity only.  Drudes 
D1 and D3 are technically separated by 4 bonds, but atoms A2 and A3 are only 
separated by 2, so the 1-3 interactions between A1, D1, A3, and D3 are all 
excluded and treated by the screening function.  That is, a bond to a 
Drude/shell does not count as a bond when it comes to determining whether or not 
interactions should be excluded.

> Since undocumented, I do not really know in which (atom-atom,
> atom-shell and shell-shell) cases do we have Coulombic,
> point charge-dipole and dipole-dipole interactions and whether
> we can switch them on/off. Finally, only I expect from the c/s model
> to make atom-atom interactions more sophisticated: they
> have an extra channel of communication -between their
> polarizable aspect.
>

You can switch any of those off with exclusions, which work on nonbonded 
interactions.  Since the Thole screening is considered a bonded interaction in 
the code (very convenient), it selectively turns back on the (screened) 
interactions between neighboring dipoles that are otherwise excluded from other 
nonbonded interactions.

> GROMACS will crash if I put shells exactly at the atomic positions
> thus I use some 0.01 nm offset (in random direction to avoid
> metastable symmetric structures). I do not expect catastrophic
> polarization with such small offsets (small compared to typical
> bond lengths of 0.2 nm) but something goes wrong as both shell
> and atomic optimizations (EM run) stop after a few steps without
> approaching the force criterium.
>

A separation of 0.2 nm is very large for polarizable systems; our model 
considers polarization catastrophe if the Drude moves > 0.02 nm from the atom. 
I don't know what's expected for your system, but that's quite large.

>> As it should; the value of a should be the sum of the two atomic Thole
>> factors, or 2.6, whichever your force field uses.
>
> Thole screening was not mentioned in the paper I started from.
>
> I started from a=2.6 which was suggested in the GROMACS Manual.
> Later I noticed that this value was estabilished for ethanol and such.
> I've seen values between ca. 0.5 and 2.6 so I played between 0.1
> and 10 to see that it has a strong effect on forces and the final structure
> and see that it is to be calibrated for a new material.
> I understand that a controls the drop of interaction with distance.
> Maybe I need to use unusual values.
>

Indeed, sticking to 2.6 is usually just a default and needs refinement.  We 
typically don't have values much higher than 6 or lower than 0.25 or so.

>> Sounds like it could be incorrect forces, but without seeing the Thole
>> section of your topology, I can't say for sure.
>
> The Thole section in my TOP file lists all c/s pairs that are closer
> than R (eg. 0.5 nm) to one another (also specifying func=2, a and alpha values).
>
>> 0.5 - 1.0 nm sounds like far too large of a radius for Thole interactions.
>
> Should I use something like a range [0.3;0.5]? (Typical bonds are
> shorter than 0.3 nm.) Or should I use [0;0.3] to just save the neighbors
> from the catastrophe?
>

A reasonable start.

I can also send you a modified version of the Thole screening code to see if 
that improves the situation; I fought with a similar issue for days and it came 
down to incorrect scaling of forces (at least, inconsistent with how CHARMM does 
it, so that's why my force field was failing ;)

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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