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

[Tamas Karpati]

Dear Justin,

I'm in the experimenting with [exclusions] for [thole_polarization].
Still not there but approaching... I'd like to ask more help as the
Manual is not detailed enough for me to comprehend the whole.

When using the core/shell model I expect the following extra
effects and functions beyond the normal simulations.
[Denotations: A = core atom, D = Drude (or shell) particle, AD_i =
the i_th A-D oscillator (dipole).]


  1. With a q_D charge on D, the "normal" q_A becomes q_C=q_A-q_D
to preserve the q_A=q_C+q_D atomic charge in a "c/s" simulation
(since D is attached to A). I assume that GROMPP expects q_C for
atoms, q_D for shell particles in the GRO files. Is this correct?


  2. I assume that GROMACS evaluates _Coulomb_potentials_ only between
atoms and only using q_A values, while q_D never ever occurs in such
expressions. (Of course, technically q_A is obtained as q_C+q_D.)

Is this correct or -on top of each (Ai,Aj) pair- all (Di,Dj), (Di,Aj)
and (Ai,Dj) contributions are added to the above Coulomb sum?
(See also the merge-comment below 3B.)


  3. I expect that the _polarization_potential_ is calculated by
either one of the following 2 possibilities (A or B) after the D positions
have been optimized in each EM or MD step (I guess forces are
evaluated similarly for that optimization).

A) The AD_i - AD_j pair potentials are summed using q_D and -q_D values
for the i_th and j_th dipoles. For a given pair one has Coulombic terms
for the (Ai,Aj), (Ai,Dj), (Aj,Di) and the (Di,Dj) charged point pairs.
I assume that the (Ai,Di) and (Aj,Dj) terms are only used in the shell
particle optimizations.

[Nothing is added twice, eg. in (Ai,Aj), as the "normal" (non c/s) part
used the q_A charges for the Coulomb interaction, while here we have
+/-q_D charges.]

B) The mu_i and mu_j dipole moments of the AD_i - AD_j pairs are summed
(of course, intrinsically calculated from the i_th and j_th q_D charges
and A-D distances).

Is A or B correct? Which one?
(I can imagine this polarization part merged with the Coulomb part
in 2. for efficiency but I'd like to make the conceptual distinction.)


  4. If the i_th and j_th dipoles are too close to one another then
the "polarization becomes catastrophic" which is indicated by too high
Coulombic forces, mu_AD dipole moments or A-D distances (not infinite
polarizabilities as alpha is an input variable).

Here is the entry point for Thole screening: for dipoles in one another's
proximity scale down the Coulomb forces between just the two Drude (shell)
particles. The scale factor is 1-(1+r12/2)*exp(-r12), where
r12=a*d/(alpha_1*alpha_2)^1/6 with
    a: magic number, calibrated for each material,
    d: distance of D_i and D_j and
    alpha_x: polarizabilities [nm^3].

Thus shell-shell (Di,Dj) catastrophe have been avoided but what is with
the remaining (Ai,Aj), (Ai,Dj), (Aj,Di) terms?

(I have tried using the [excludes] directive to remove these but my
simulation still couldn't converge. I also excluded the Di-Dj, Ai-Di
and Aj-Dj pairs with not more success. I'm not even convinced that
the inner-loop optimization for the shell particles completes so
I've set niter=10000.)


Also confusing is that -looking at the MDRUN output and the Fmax values
that increase from the initial 1e4 to even 1e16 during the EM process-
I notice both normal atoms and Drude particles having the highest force
acting on them in one or the other EM step.
Besides, Epot is gradually decreasing so I'm getting more
happy with the results. However, full convergence is not achieved.


My goal is to learn what to [exclude] and [thole_polarization] in my
TOP file after locating pairs of AD particle doublets that are too
close to each other. I imagine that for each AD_i-AD_j pairs I need
to do the same.
Could you please shed more light on what "same" means here?

I thank you for your continuous support.

With regards,
  toma

PS: invented a work-around for GROMPP to accept the
[exclusions] directive in the TOP file when no bonds are
present (the crystal is held together by Buckingham+Coulomb):
one needs to just write an empty [bonds] directive (without
a single bond specified). Otherwise it claims the [exclusions]
is at a wrong position in TOP.