[gmx-developers] Excluded charges and PME correction

[David van der Spoel]

On 5/3/10 2:54 PM, Berk Hess wrote:
> 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
>
Same with molecule displaced to smaller X yields negative force:
    box (3x3):
       box[    0]={ 2.00000e+00,  0.00000e+00,  0.00000e+00}
       box[    1]={ 0.00000e+00,  2.00000e+00,  0.00000e+00}
       box[    2]={ 0.00000e+00,  0.00000e+00,  2.00000e+00}
    x (2x3):
       x[    0]={ 5.00000e-01,  1.20000e+00,  1.50000e+00}
       x[    1]={ 5.00000e-01,  8.00000e-01,  1.50000e+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]={-8.38504e+01,  3.04622e+03,  8.38234e+01}
       f[    1]={-8.38604e+01, -3.04622e+03,  8.38232e+01}


    nkx                  = 23
    nky                  = 23
    nkz                  = 23
    pme_order            = 4
    ewald_rtol           = 1e-05
    ewald_geometry       = 0
    epsilon_surface      = 80
    rlist                = 0.9
    rlistlong            = 0.9
    rtpi                 = 0.05
    coulombtype          = PME
    rcoulomb_switch      = 0
    rcoulomb             = 0.9
    vdwtype              = Cut-off

OK, I now made a series with different X-coordinates for the molecule 
and I plotted the force for different grid spacing. Results are here:
http://folding.bmc.uu.se/images/forceX-all.png
http://folding.bmc.uu.se/images/forceX-all.xvg

It seems that the forces go down systematically with decreasing spacing 
(as they should of course). However in magnitude they are still big, 
even though the dipole is huge in this case (4D, 0.4 nm distance charge 
+/- 10). Then this probably means that all is fine anyway.

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