GMX dipoles for charged molecules - center of mass?

[Jacek]

Dear all,

the documentation for GMX DIPOLES says that “For molecules with a net charge, the net charge is subtracted at center of mass of the molecule”. (gmx dipoles - GROMACS 2024.3 documentation)

If I’m not mistaken this subtraction of net charge at the COM is incorrect and it should be the center of geometry - see uploaded file with a quick derivation for a two atomic molecule.

I am now wondering, if it’s just the documentation or if the COM is really wrongly implemented in the code (or if I’m having a trivial mistake in my thinking). Thanks!

Best,
Jacek

[MagnusL]

Thanks! Would you be able to open an new issue at Issues · GROMACS / GROMACS · GitLab? It seems like this is something that should be fixed, or at least discussed whether the current implementation is correct or not.

[Jacek]

Hi Magnus,

thanks for looking at it, I just wrote a ticket.

Just for everybody who’s interested:
I tested this (below) and there indeed seems to be problem in gmx dipoles, that for charged molecules/systems a wrong dipole is determined.

I used a trajectory a two-atomic CO molecule (GROMACS 2021.5; same result for earlier versions; the newest version still seems to use COM - gromacs/src/gromacs/gmxana/gmx_dipoles.cpp at 3ffa4ddb141f847edbcd42c9fa88dc820ad70b0b · gromacs/gromacs · GitHub). Keeping the coordinates of the initial trajectory, I added or subtracted the net charge. In that way, having the same geometry, I had two cases with similar charge separation:

Case 1 - net charge = 0:
q_C = 0.633601
q_O = -0.633601

Case 2 - net charge = -1:
q_C = 0.133601
q_O = -1.133601

Assuming correct treatment of dipoles, GROMACS should give the same result since net charge is subtracted. However, results are different:

For case 1, GMX DIPOLES gives me:

Dipole moment (Debye)

Average = 3.6924 Std. Dev. = 0.0656 Error = 0.0001

The following averages for the complete trajectory have been calculated:

Total < M_x > = 0.365268 Debye
Total < M_y > = 0.08411 Debye
Total < M_z > = -0.0795452 Debye

Total < M_x^2 > = 4.49683 Debye^2
Total < M_y^2 > = 4.21801 Debye^2
Total < M_z^2 > = 4.92327 Debye^2

Total < |M|^2 > = 13.6381 Debye^2
Total |< M >|^2 = 0.146823 Debye^2

< |M|^2 > - |< M >|^2 = 13.4913 Debye^2

Finite system Kirkwood g factor G_k = 0.989547
Infinite system Kirkwood g factor g_k = 0.9745

Epsilon = 1.0478

In Case 2, I obtain:

Dipole moment (Debye)

Average = 7.4619 Std. Dev. = 0.1802 Error = 0.0002

The following averages for the complete trajectory have been calculated:

Total < M_x > = 0.775687 Debye
Total < M_y > = 0.105891 Debye
Total < M_z > = -0.0110759 Debye

Total < M_x^2 > = 18.5724 Debye^2
Total < M_y^2 > = 17.3047 Debye^2
Total < M_z^2 > = 19.8354 Debye^2

Total < |M|^2 > = 55.7125 Debye^2
Total |< M >|^2 = 0.613027 Debye^2

< |M|^2 > - |< M >|^2 = 55.0995 Debye^2

Finite system Kirkwood g factor G_k = 0.989573
Infinite system Kirkwood g factor g_k = 0.935698

Epsilon = 1.19521

Best,
Jacek