[gmx-developers] Excluded charges and PME correction

[David van der Spoel]

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.


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se