[gmx-developers] Excluded charges and PME correction

[Berk Hess]

David van der Spoel wrote:
> On 5/3/10 11:39 AM, Berk Hess wrote:
>> David van der Spoel wrote:
>>> On 5/3/10 11:15 AM, Berk Hess wrote:
>>>> If you have exclusions in your topology, the ewald_LRcorrection will
>>>> completely remove
>>>> the mesh electrostatic forces, up to the mesh precision.
>>>> But a single dipole in a box will feel the effect of the periodic
>>>> images.
>>>> Depending on the size of your box, this can be a large force.
>>>
>>> OK, but isn't it strange that this leads to forces perpendicular to
>>> the dipole as well (the box is cubic)? And shouldn't the dipole
>>> correction (epsilon_surface != 0) take care of this force?
>> It all depends on your box matrix and PME accuracy settings.
>> The dipole correction never fully takes care of this, it simply reduces
>> it by 1/epsilon.
>> The 3dc geometry will remove most of it, but not everything.
>>
>> I guess the perpendicular forces should be (close to) zero, unless your
>> box is triclinic.
> It seems that it is worse, forces are position-dependent. If my
> molecule is in the center of a cubic box of 2 nm I get the
> perpendicular force close to zero (step 0):
>    x (2x3):
>       x[    0]={ 1.00000e+00,  1.20000e+00,  1.00000e+00}
>       x[    1]={ 1.00000e+00,  8.00000e-01,  1.00000e+00}
>    v (2x3):
>       v[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>       v[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
>    f (2x3):
>       f[    0]={ 0.00000e+00,  3.17667e+03,  8.47711e-04}
>       f[    1]={-2.28882e-03, -3.17667e+03, -8.47711e-04}
>
> If I change the grid spacing from 20 (0.1 nm) to 21
>    f (2x3):
>       f[    0]={-4.50611e-03,  3.12076e+03,  6.21270e-03}
>       f[    1]={-5.40733e-03, -3.12078e+03,  2.74640e-03}
>
> to 22
>    f (2x3):
>       f[    0]={ 6.34670e-03,  3.37936e+03,  0.00000e+00}
>       f[    1]={ 2.51770e-03, -3.37933e+03,  0.00000e+00}
>
> to 23
>    f (2x3):
>       f[    0]={-3.29018e-03,  3.05586e+03,  5.31634e-03}
>       f[    1]={ 1.07188e-02, -3.05586e+03,  5.31634e-03}
>
> Now if I move the coordinates by 0.3 nm in the X-direction I get:
>    f (2x3):
>       f[    0]={-3.65754e+01,  3.04537e+03, -4.07583e-04}
>       f[    1]={-3.65607e+01, -3.04536e+03, -4.07583e-04}
>
> Or by 0.5 nm in X:
>    f (2x3):
>       f[    0]={ 8.38358e+01,  3.05105e+03,  1.18744e-03}
>       f[    1]={ 8.38223e+01, -3.05105e+03,  1.18744e-03}
>
> Or  by 05 nm in X and 0.5 nm in Z:
>    f (2x3):
>       f[    0]={ 8.38596e+01,  3.04622e+03,  8.38235e+01}
>       f[    1]={ 8.38329e+01, -3.04621e+03,  8.38232e+01}
>
>
> So the force in the Y direction depends on grid spacing and in X and Z
> on the position.
> As you see, the above force will lead to a net translation.
>
>
Can you post the box and all PME parameters?
These are all things I would expect, but the magnitude depends on your
settings.

Berk