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

[Justin Lemkul]

On 7/31/14, 12:35 PM, Tamas Karpati wrote:
> 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?
>

No.  Every atom in the topology is assigned its charge.  grompp doesn't 
manipulate them.  You need to assign q_A and q_D, e.g. as is done in sw.itp. 
The contents of [polarization] are only used to calculate the force constant of 
the bond between the atom and Drude (see equation in sw.itp or code in 
src/gromacs/gmxlib/bondfree.c)

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

Anything with a charge contributes to the Coulomb potential.

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

Neither.  The polarization energy is the output of the polarize() function 
(again, see src/gromacs/bondfree.c) and is a simple bonded potential using the 
atom - Drude distance and the value of k (force constant) that is computed from 
q_D and alpha.

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

Thole screening is applied to all possible pairs.  Remember, it's dipole-dipole 
screening, so all of those terms are screened, hence the reason you need to 
specify 4 atoms in the [thole_polarization] directive.  See thole_pol() and 
do_1_thole() in bondfree.c.  For once, the CHARMM documentation is actually 
better than Gromacs (though I am in the process of fixing that):

http://www.charmm.org/documentation/c34b1/drude.html

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

This still suggests topological instability.  Or, as I said before, incorrect 
force calculation that is in do_1_thole() that I found to cause failures in 
otherwise sensible input files for our FF.

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

Exclude nearest neighbors, apply Thole screening to them.  Depending on the 
geometry of the system, perhaps second neighbors, as well.  For instance, if 
atom 1 is an atom, 2 is its Drude, 3 is an atom, and 4 is its Drude:

[ exclusions ]
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3

[ thole_polarization ]
; ai  aj  ak  al  func  etc.
    1   2   3   4    2   a  alpha_12  alpha_34

-Justin

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

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

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

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