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

[Tamas Karpati]

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

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

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?

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

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.
This, however, remained hidden from me (only appeared in the binary TPR file).
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.

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?

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

> 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