ALPB and ddCOSMO not rotationally invariant

[alanchienSDGR]

Describe the bug
xtb version 6.7.1 (d42779f) with cpcm-x has small, but significant, energy differences (0.0001 Hartree in the attached example) with respect to rotation of a simple molecule. The energies match up to the 7th decimal when not using cpcm-x.

To Reproduce
Run xtb {input.sdf} --cpcmx water on both the attached ref.sdf and rotated.sdf input files (I can't upload .sdf files apparently so I converted them to txt)
ref_sdf.txt
rotated_sdf.txt

rotated.sdf is just a rotation of the molecule in ref.sdf around the x-axis by 30 degrees.

ref.sdf output with cpcmx
| TOTAL ENERGY -17.598183210842 Eh |
rotated.sdf output with cpcmx
| TOTAL ENERGY -17.598280276499 Eh |

ref.sdf output without cpcmx
| TOTAL ENERGY -17.591235057712 Eh |
rotated.sdf output without cpcmx
| TOTAL ENERGY -17.591235064338 Eh |

Expected behaviour
I expected the total energies to match up to about the 6th or 7th decimal place, similar to the accuracy of the dispersion/nuclear repulsion energies and when running without cpcmx.

[lukaswittmann]

Thank you for reporting.

Upon review, we noticed the files you provided were not only rotated but also exhibited an RMSD of 1E-6. For comparison, we conducted a test using the following structures with an RMSD of 1E-11:
struc.xyz.txt
rotated.xyz.txt

Test Results:

(GFN2, solvent: water, energies in Eh)

	xTB gas	xTB ddcosmo	xTB cpcmx	xTB alpb	tblite ddcpcm	tblite alpb
ref	-28.30391089	-28.45055007	--	-28.42249445	-28.36792386	-28.44749431
rot	-28.30391087	-28.45043739	-28.47818823	-28.42243534	-28.36814497	-28.44749430
Difference	1.65E-08	1.13E-04	--	5.91E-05	-2.21E-04	5.09E-09

Observations:

The issues already arise using only the underlying ddCOSMO/ddCPCM which is likely carried over into CPCMX.
We already plan to revise our implementation by implementing ddX and thus replacing our native ddCOSMO/ddCPCM implementations.

Unfortunately, ALPB also does not show rotational invariance in xTB either.
The tblite version is fine only because it does not yet have the corresponding ALPB CDS term implemented.
The CDS term includes an integration over a Lebedev grid surface, which seems to be the issue. ddCOSMO/ddCPCM also utilizes the corresponding routines/grids.

We will further look into this issue.

CC: @thfroitzheim

[alanchienSDGR]

Hi @lukaswittmann - thanks for confirming the bug! Looking forward to the fix.

In the meantime, are you aware any recommended settings that might mitigate this error? Maybe an option to increase the grid size?

[haneug]

You can try to increase the grid size: https://xtb-docs.readthedocs.io/en/latest/gbsa.html#available-grids

[lukaswittmann]

Increasing the grid does improve the result: Using ddCPCM in tblite with the largest grid (ld5810, equiv. to Gridlevel=extreme in xTB), the error is around 5E-6, which is more than sufficient.

In principle, the issue is known for apparent surface charge PCMs.
Rotational invariance is not strictly preserved since each octahedral Lebedev grid is generated in the laboratory frame before being translated to the target atomic center. While rotational invariance is achievable with an infinitely dense grid. Using Lebedev grids generally results in errors that are minimal enough for practical purposes, i.e. far below <0.1 kcal/mol, normally even much lower (Fig. 4. in the SwiG paper).
However, >1E-4 Eh difference is too high, we would like something more like ~1E-5 Eh; we will look into grid settings and test some ideas.

Regarding the ALPB issue, we have to look more into the details, but the ORCA SMD CDS part is rotationally invariant, which means this should also be applicable to the ALPB CDS part also.

[corinwagen]

This problem can be solved if you just specify a "standard" molecular orientation, see Johnson/Gill/Pople 1994