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

[Justin Lemkul]

On 7/29/14, 9:33 AM, Tamas Karpati wrote:
> Dear Justin,
>
> I've tried applying the [thole_polarization] section in
> the format you've suggested. I've noticed the followings.
>
>
> 1. First I made a mistake and forgot to remove the [polarization]
> section from the TOP file. Even this way MDRUN worked and indicated
> the usual ca. 10 kJ/mol contribution from polarization.
> Besides a new contribution of ca. 600 kJ/mol appeared from Thole pol.
>
> I switched the former off to have just one kind of pol effect.
> Is this correct?
>

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

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.

>
> 2. I assumed that in the [thole_polarization] section of the TOP
> file I needed to select which dipole pairs are to be included on
> basis of my "chemical intuition" -if I had any. I tried a few radii
> of "sphere-of-inclusion" but saw little differences.
>
> Is it a correct selection for Thole screening? (I have a crystal
> with separate ions held together by Buckingham forces and Thole
> polarization -there are no bonds at all.)
>
> Do I have a way to switch on/off intercell dipole interactions?
>

I don't know what the right answer is here.  What is the source of the force 
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.

>
> 3. Playing with the "a" factor of the Thole scheme resulted in more
> pronounced changes in the contribution from polarization but I still

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.

> see too strong Coulombic interactions. Omitting shell particles and
> thus polarization do scale Coulombics down by a factor of 100 or 1000.
> I guess that using [polarization] all dipoles interact with all the
> others (although no idea how dipoles interact through walls of the cell).
> I imagine that in the [thole_polarization] case only those dipole pairs
> will interact that I have mindlessly selected.
>
> In contrast with my expectation I noticed an order of magnitude
> greater Thole pol. values than Polarization values in the simple
> pol. case. (For your information, I have ca. 2000 atoms + 1000 shell
> particles in a box of ca. 3x3x3 nm. In case of the Thole scheme
> I have ca. 1e5 pcs. of dipole pairs with a 1 nm radius for inclusion
> while ca. 15000 for a radius of 0.5 nm.)
>

Sounds like it could be incorrect forces, but without seeing the Thole section 
of your topology, I can't say for sure.

0.5 - 1.0 nm sounds like far too large of a radius for Thole interactions.

-Justin

> The total potential is almost a pure Coulombic one and is too large.
> Do you have any hint on what could go wrong?
>
>
> Thanks for your attention.
>
> With best regards,
>    toma
>
>
> On Fri, Jul 25, 2014 at 9:31 PM, Tamas Karpati <tkarpati at gmail.com> wrote:
>> Dear Justin,
>>
>> Thank you for the valuable information. I'm starting experimenting with it.
>>
>> Have a nice weekend,
>>    toma
>>
>>
>> On Fri, Jul 25, 2014 at 7:53 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>>
>>> On 7/25/14, 12:40 PM, Tamas Karpati wrote:
>>>>
>>>> Dear Justin,
>>>>
>>>> Thanks for your educational answers.
>>>>
>>>>> Coulombic interactions fail at short distances; you probably need to
>>>>> apply
>>>>
>>>>
>>>> I was afraid of that... somehow removing shells from the cores in the
>>>> initial structure have let it "functionally" work (with not yet
>>>> reasonable results).
>>>>
>>>>> Thole screening to avoid polarization catastrophe.  Ions are particularly
>>>>> problematic in this regard.
>>>>
>>>>
>>>> I've seen this mentioned in the Manual, but hadn't ever hit GROMACS
>>>> specific details on how to apply polarization in the input files.
>>>> The only source I could locate is within the GROMACS package,
>>>> under the name "sw.itp". It does exclusively implement the
>>>> so called "water polarization" model -at least I think so.
>>>>
>>>
>>> The "water polarization" function is a water-specific anisotropy function.
>>> Don't try to use it for anything else; the interpretation of the atom
>>> numbers for local axis construction are very specific.
>>>
>>>
>>>> Can you please direct me to a source from which I could learn
>>>> how to polarize GROMACS? I was'not lucky on the Internet and,
>>>> indeed, even at the GROMACS site or Manual (neighter application
>>>> examples nor file format descriptions).
>>>>
>>>
>>> The Thole screening function is (in the released version) not used by
>>> anything, so it's not documented.  In its present incarnation, you need a
>>> [thole_polarization] directive that lists atom-shell/Drude pairs as follows:
>>>
>>> atom_i shell_i atom_j shell_j 2 a alpha_i alpha_j
>>>
>>> The "2" is a required function type.  My implementation of the CHARMM Drude
>>> FF is nearly done, and there are changes to the way the Thole directive is
>>> laid out in the future, but at the moment (up through version 5.0), this is
>>> the way it works.  The code is in src/gromacs/gmxlib/bondfree.c.
>>>
>>>
>>> -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
>>>
>>> ==================================================
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
>>> mail to gmx-users-request at gromacs.org.

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

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

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