Fix bug in GHF-X2C dipole function (pyscf#2396)

* Fix bug in GHF-X2C dipole function (fix pyscf#2393) * typo
sunqm · Sep 7, 2024 · c67deaa · c67deaa
1 parent 5417cc2
commit c67deaa
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 7 deletions.
diff --git 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
@@ -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