cDFT - potential issue with translational invariance

[Gabriele Cesare Sosso]

Hi all.

I'm playing around with cDFT, with the idea of quantifying the extent of charge transfer within a certain dimer. 

In this tutorial, the coordinates of the (water) dimer are centred (via &CENTER_COORDINATES ON) when computing the energy of the system via either conventional DFT as well as cDFT. However, the coordinates are not centred (i.e. &CENTER_COORDINATES OFF) when computing the energy (and, most importantly, the electron and spin densities)  for the two fragments (i.e. the two water molecules) separately. 

To me (I would be very happy to be corrected!), this makes sense only if the center of mass of the dimer sits exactly in the middle of the simulation box. If that's not the case, at the stage of doing the actual cDFT calculation, when reading the electron/spin densities from the two .cube files (one for each fragment) these will not be centred correctly.

I have tested this, using the above mentioned water tutorial. Inputs/outputs attached. Here's what I got:

Tutorial (unmodified)

Target value of constraint  :                             7.994128187114

Current value of constraint :                             8.016256807586

Charge transfer (ct_energy) = 2.058885794043000, mHa

Shifted (atomic positions translated by 1 angsts in the x direction)
Target value of constraint  :                             9.681231144409
Current value of constraint :                             8.016256807586
ct_energy = 1552.849968090016000, mHa

Shifted_no_center (as Shifted, but without centring anything)
  Target value of constraint  :                             7.994129688685
  Current value of constraint :                             8.016241161786
  ct_energy = 2.056212410175000, mHa

Com (as Tutorial, but the center of mass of the dimer sits exactly in the middle of the box)
Target value of constraint  :                             7.999862917720
Current value of constraint :                             8.016256821935
ct_energy = 1.116896285772000, mHa

Com_no_center (as Com, but without centring anything)
Target value of constraint  :                             7.994129127081
Current value of constraint :                             8.016254935773
ct_energy = 2.058503499761000, mHa

It looks like:

• Centering the coordinates of the dimer (DFT and cDFTcalculations) but not those of the isolated fragments is dangerous, as the result is not translationally invariant (see Tutorial vs Shifted). It is clear that the target value of the constraint (which is computed from the superposition of the electron/spin density from the .cube files) is off. 
• Not centring any coordinate at any stage seems to be giving fairly consistent (but not exactly consistent) results.
• When adopting the Tutorial approach and moving the center of mass of the dimer in the middle of the box, we still get a result which appears to be incorrect. Apologies for probably being rather dense, here, but I don't immediately understand why, particularly because the target value of the constraint seems to be correct (and yet the charge transfer energy is not).

Any advice would be very much appreciated!

Best wishes,

Gabriele