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

[Justin Lemkul]

On 7/30/14, 9:37 AM, Tamas Karpati wrote:
>> But in short, no, bonds are not compulsory for using Thole, but exclusions
>> definitely are.  I would exclude at least the nearest-neighbor interactions.
>
> Does this mean that by [thole_polarization] I define bonded interactions
> which can be dropped by [exclusions]? Would it not be easier to not
> define Thole at all for the nearest-neighbors in the first place?
>

No, exclusions work on nonbonded interactions.  Thole screening is, in theory, a 
nonbonded interaction, but it is handled as a bonded interaction in the code, 
which allows it to act on otherwise excluded atom pairs.

>> Can you post a snippet of your [polarization] and [thole_polarization]
>> directives?
>
> Sure. Shell particles directly follow their corresponding core (part. type "A").
> Atoms 1 and 3 are metallic, 11 and 13 are nonmetallic elements. For the
> rest of the atom types (ions) there are no polarizations prescribed.
>
> [ polarization ]
> ;   ai    aj funct   alpha
>       1     2     1   6.7
>       3     4     1   6.7
>      11    12     1   1.5
>      13    14     1   1.5
> ...
>
> [ thole_polarization ]
> ; atom_i shell_i atom_j shell_j func a alpha_i alpha_j
> 1 2  3  4 2 2.6 6.7 6.7
> 1 2 11 12 2 2.6 6.7 1.5
> ...
>

I think your units are wrong here.  The alpha values should be in nm^3.  Those 
values of alpha are about 1000x larger than I think they should be.

> I assume it is irrelevant whether I use S or D particle type, thus I kept S.
> In my TOP file there are [moleculetype], [atoms], [angles], [polarization],
> [thole_polarization], [system] and [molecules] sections, nothing else.
>

Keep S.  A "dummy" (D) particle is a virtual site; very different!

> I quasi-manually put polarization on 2 types and put Thole on all pairs
> of A-S pairs that are within 0.35 nm.
> Should I then purge via [exclusions] all that are within "one bond distance"?
> Or should I just not put Thole there?
>

Anything that is treated by Thole should be listed as an exclusion.

> The Manual sais "The exclusions for non-bonded interactions are
> generated by grompp for neighboring atoms up to a certain number
> of bonds away..." which is embarassing for two reasons:
>    1. [thole_polarization] is considered as bonded interaction

I don't see how that is inconsistent with the manual.

>    2. I create GRO and TOP files then feed them (with MDP) to GROMPP
>        to generate the TPR file. This is to control (=understand) every
>        possible details of what's goind on.
> Do I miss something happening in the TOP --[GROMPP]--> TPR conversion?
>
>> Just a note for clarity here, since it might be a bit confusing.  A Drude
>> and a shell are the same thing; "dummies" (virtual sites) haven't been
>> mentioned yet and don't sound like they're relevant.  A-S is an
>> atom-shell(Drude) bond and D-S is a dummy(virtual site)-shell bond.  The use
>> of S and D to mean different things when they could be the same in other
>> nomenclature systems is somewhat confusing.
>
> Thanks for clarifying that.
>
>> What I meant was this (hopefully this shows up correctly)
>>
>> D1   D2   D3
>> |    |    |
>> A1---A2---A3
>>
>> The exclusions and Thole screening are based on atom connectivity only.
>> Drudes D1 and D3 are technically separated by 4 bonds, but atoms A2 and A3
>> are only separated by 2, so the 1-3 interactions between A1, D1, A3, and D3
>> are all excluded and treated by the screening function.  That is, a bond to
>> a Drude/shell does not count as a bond when it comes to determining whether
>> or not interactions should be excluded.
>
> Thanks for the figure. (Do you mean A1 and A3 are separated by 2 bonds?)
>

Yes.

> As I understand, Thole scr. isn't based on but rather it creates a kind of
> atom connectivity which [exclusions] can partially remove. Eg. A1-D3
> and D1-A3 get deleted for being 1-3 interactions. This is done by GROMPP
> automagically applying [exclusions] to my TOP file which doesn't have any.

Have you verified this by using gmxdump?  Exclusions are generated from bonds. 
You said you don't have any bonds, so there shouldn't be any exclusions aside 
from those you are adding.

> This, however, remained hidden from me (only appeared in the binary TPR file).

That's what gmxdump is for.

> Then MDRUN's EM job applied [thole_polarization] rather than [polarization]
> on the excluded connections.
>
> On the other hand, the Manual explains in "4.4.3. Thole polarization"
> that V_Thole is defined between 2 shell particles. In this case, I don't
> understand why to include atom indices in the [thole_polarization] section.
> I could simply define [polarization] for all my A-S pairs and
> [thole_polarization] for all the S-S pairs that are too close -and done.
>

Thole screening is applied to atom-shell(Drude) pairs to account for unphysical 
dipole-dipole interactions.  The manual is in need of updating.

> Even if I misunderstand something, for my simple scheme
>
> XYZ --[my code]-> GRO + TOP --[GROMPP]--> TPR --[MDRUN]--> ...
>
> the automatisms in GROMACS seem to be overcomplicated.
> Thank you for helping me grabbing concepts and usage together:)
>
> Do you think there is a way for my code to directly create everything
> necessary in the TOP file and avoid GROMPP doing something in the dark?
>

I still don't understand what's "in the dark" here.

>> You can switch any of those off with exclusions, which work on nonbonded
>> interactions.  Since the Thole screening is considered a bonded interaction
>> in the code (very convenient), it selectively turns back on the (screened)
>> interactions between neighboring dipoles that are otherwise excluded from
>> other nonbonded interactions.
>
> So first [exclusions] then [thole_polarization] is applyied?
> Does this mean that [exclusions] turn Buckinghams off then
> [thole_polarization] takes place defining the screening function?
> I do not know if killing down Buckingham is of any good...
>

That's why I have said a few times that I am hesitant to be making all of these 
recommendations.  I don't know what paper you're referring to, what model, or 
what software.  Ad hoc changes are probably a bad idea, but in the polarizable 
models I have used, nearby interactions are extremely sensitive to dipole-dipole 
interactions.

-Justin

>> A separation of 0.2 nm is very large for polarizable systems; our model
>> considers polarization catastrophe if the Drude moves > 0.02 nm from the
>> atom. I don't know what's expected for your system, but that's quite large.
>
> Thank you. I took it down to 0.005.
>
>> Indeed, sticking to 2.6 is usually just a default and needs refinement.  We
>> typically don't have values much higher than 6 or lower than 0.25 or so.
>
> It is nice to know in advance, thanks a lot.
>
>> I can also send you a modified version of the Thole screening code to see if
>> that improves the situation; I fought with a similar issue for days and it
>> came down to incorrect scaling of forces (at least, inconsistent with how
>> CHARMM does it, so that's why my force field was failing ;)
>
> We all wish the bests for bug hunting :)
>
> I appreciate your aid pretty much.
>
> Best regards,
>    toma
>

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

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

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