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

[Justin Lemkul]

On 7/31/14, 2:17 PM, Tamas Karpati wrote:
> Dear Justin,
>
>>>     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)
>
> If I have a Mg++ ion with q_D=+3.2 then for the shell particle I'd set
> +3.2 and for the Mg atom I should stay with +2? So far -in order to maintain
> charge neutrality of the system- I've been using -1.2 for Mg to have the
> correct +2e for the AD duett. Big difference!
>

Yes, the Drude and core atom will have charges of opposite sign to yield the 
total charge on the pair.  A positively charged auxiliary particle is a bit 
weird, but not wrong, if the model was parametrized that way.  Normally the 
auxiliary particle represents deformation of the electron cloud, so it is 
negative while the core atom is positive.

>> Anything with a charge contributes to the Coulomb potential.
>
> OK, I see now why to [exclude] lots of things.
>
>> 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.
>
> Thanks. I'm going to check the source code.
>

Small correction (sorry for the typo) - the file is src/gromacs/gmxlib/bondfree.c.

>> Thole screening is applied to all possible pairs.  Remember, it's
>
> With "possible" meaning close-enough dipoles?
>
>> dipole-dipole screening, so all of those terms are screened, hence the
>
> I think I took words too sharply... I begin to realize that dip-dip screening
> means point charge-point charge screening. Thus, eg. with AD_i and AD_j
> close enough we have 4 particles and 4*3/2=6 point-to-point interactions
> to screen. (I did it, not yet there.)
>
>> 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'm sure it worths looking at. Appreciated.
>
>> 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.
>
> I'm looking forward to seeing the update :)
>

It will be a while.  It works, but is not yet compatible with DD (coding 
ongoing) and needs to be tested quite a bit.  Hopefully soon...

>> 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
>
> Ooops... I expected finding something like that in a tutorial.
> But I see that you specify pairs twice: 1-2 and 2-1 etc.
> Maybe I just "eliminated" half of the polarization and screening
> potentials?
>

The way I have done the exclusions is redundant and in principle not totally 
necessary.  Atom-Drude pairs connected via bonds or within [polarization] 
directives (e.g. 1-2 and 3-4) should already be excluded from one another, so in 
reality:

[ exclusions ]
1 3 4
2 3 4

should cover it.

-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

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