GMX dipoles for charged molecules - center of mass?

User discussions
1

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.

GMX dipoles - center of mass
GMX dipoles - center of mass1496×1560 112 KB

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

2

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.

3

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

4

my 5 cents:
Regarding what you write above, fact is that in most papers in the literature the center of mass is used as a reference point when calculating the dipole moment of a protein (or any large molecule with a net charge). The center of geometry might be a better idea as you suggest in your example, but I have not seen this being applied anywhere in the literature. Also, there is at least one more alternative, the so-called center-of-diffusion, see:

not sure how much difference this makes, but the authors claim that this should be used.
Anyway, my impression is that the gmx dipoles tool could perhaps do with some revision - I just posted this: Gmx dipoles and index groups vs molecules
which suggests that there might be undocumented limitations of this tool.
michael