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

Tamas Karpati tkarpati at gmail.com
Wed Jul 30 11:15:21 CEST 2014

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.

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

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

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.

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.

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

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

I rewrite my code to co-generate [polarization] with [thole_polarization]
and experiment a little with the what-to-include and the value of a.

Thanks for your co-operation.

With best regards,

More information about the gromacs.org_gmx-users mailing list