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

[Tamas Karpati]

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!

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

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

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

With regards,
  toma