From c67deaa39900d053935ab4d3bf944fb3e2470ae0 Mon Sep 17 00:00:00 2001
From: Qiming Sun <osirpt.sun@gmail.com>
Date: Sat, 7 Sep 2024 09:49:08 -0700
Subject: [PATCH] Fix bug in GHF-X2C dipole function (#2396)

* Fix bug in GHF-X2C dipole function (fix #2393)

* typo
---
 pyscf/x2c/test/test_x2c.py |  1 +
 pyscf/x2c/x2c.py           | 19 ++++++++++++-------
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/pyscf/x2c/test/test_x2c.py b/pyscf/x2c/test/test_x2c.py
index f1d963f88a..ddca2b3206 100644
--- a/pyscf/x2c/test/test_x2c.py
+++ b/pyscf/x2c/test/test_x2c.py
@@ -215,6 +215,7 @@ def test_ghf(self):
         mf.kernel(dm0=dm)
         self.assertTrue(mf.converged)
         self.assertAlmostEqual(mf.e_tot, ref.e_tot, 9)
+        self.assertAlmostEqual(abs(mf.dip_moment() - ref.dip_moment()).max(), 0, 9)
 
     def test_undo_x2c(self):
         mf = mol.RHF().x2c().density_fit()
diff --git a/pyscf/x2c/x2c.py b/pyscf/x2c/x2c.py
index 4bd83298d3..d1b0e34ed6 100644
--- a/pyscf/x2c/x2c.py
+++ b/pyscf/x2c/x2c.py
@@ -774,17 +774,22 @@ def dip_moment(self, mol=None, dm=None, unit='Debye', verbose=logger.NOTE,
         charges = mol.atom_charges()
         coords  = mol.atom_coords()
         nucl_dip = numpy.einsum('i,ix->x', charges, coords)
-        with mol.with_common_orig(nucl_dip):
-            r = mol.intor_symmetric('int1e_r')
-            ao_dip = numpy.array([_block_diag(x) for x in r])
+        with mol.with_common_orig((0,0,0)):
             if picture_change:
                 xmol = self.with_x2c.get_xmol()[0]
                 nao = xmol.nao
-                prp = xmol.intor_symmetric('int1e_sprsp').reshape(3,4,nao,nao)[:,0]
-                prp = numpy.array([_block_diag(x) for x in prp])
-                ao_dip = self.with_x2c.picture_change((ao_dip, prp))
+                r = xmol.intor_symmetric('int1e_r')
+                r = numpy.array([_block_diag(x) for x in r])
+                c1 = 0.5/lib.param.LIGHT_SPEED
+                prp = xmol.intor_symmetric('int1e_sprsp').reshape(3,4,nao,nao)
+                prp = numpy.array([_sigma_dot(x*c1**2) for x in prp])
+                ao_dip = self.with_x2c.picture_change((r, prp))
+            else:
+                r = mol.intor_symmetric('int1e_r')
+                ao_dip = numpy.array([_block_diag(x) for x in r])
 
-        mol_dip = -numpy.einsum('xij,ji->x', ao_dip, dm).real
+        el_dip = numpy.einsum('xij,ji->x', ao_dip, dm).real
+        mol_dip = nucl_dip - el_dip
 
         if unit.upper() == 'DEBYE':
             mol_dip *= nist.AU2DEBYE