[gmx-developers] Excluded charges and PME correction

[Berk Hess]

David van der Spoel wrote:
> 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.
>
Yes, the length of your dipole is only a fifth of your box size.
Make the box ten times as large and the effect will be orders of
magnitude smaller.

Berk