[gmx-developers] Excluded charges and PME correction

[Florian Dommert]

On 03.05.2010, at 13:20, 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.
> 

A further question is the accuracy of the SPME settings. The forces in X and Z direction are in the order of 10^1 in Y-direction, while in the Y direction you the forces is in the order of magnitude of 10^3. Increasing the grid spacing does not automatically provide more accurate results as a big error in real space can occur. The point is that you have to adapt the splitting parameter, which is not directly accessible in an mdp file.

We are already working on this issue and a paper should be publish in the very near future. However I have already implemented an a priori error estimate together with Carsten in to g_tune_pme, but currently the development of the code is somehow stucked. The last time when we tested it there were some problems with OpenMPI and the g_tune_pme code. Perhaps this works already. A further point is that the parallelization of the code is not very good so far, because I have somehow problems on machines with many cpus. It seems that long blocking times occur, when using the MPI_Barrier command. As soon as I have solved this problem I will make the code accessible for everybody.

For example if I use the standard SPME parameters in the water tutorial I have an error of about 5 kJ / (mol nm), which completely stems from the error in reciprocal space. When I keep all parameters apart from the splitting parameter constant I can tune the error to about 3 kJ / mol. As this systems is not completely homogeneous the error will be slightly smaller, but it will be still too large I think.

As I assume many people refer to the standard parameters of SPME,

order 4, rtol 1e-5, rcut 0.9, spacing 0.12

the force calculation can sometimes really be an issue.

What do you think ?

/Flo 


> 
> -- 
> 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
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-developers-request at gromacs.org.

--
Florian Dommert
Dipl.-Phys.

Institute for Computational Physics

University Stuttgart

Pfaffenwaldring 27
70569 Stuttgart

Phone: +49(0)711/685-6-3613
Fax:   +49-(0)711/685-6-3658 

EMail: dommert at icp.uni-stuttgart.de
Home: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 194 bytes
Desc: This is a digitally signed message part
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20100503/5af641ce/attachment.sig>