Improve docstring and other fixes (pyscf#2377)

* Add shls_slice kwarg for pyscf.df.incore (fix pyscf#1729)

* Fix example 31-xc_value_on_grid.py (fix pyscf#1725)

* Fix FCI wfn_symmetry example (fix pyscf#1726); simplify the fci direct_spin0_symm implementation

* Update GHF-X2C docstrings (fix pyscf#1707)

* Generating XC_CODES from libxc (fix pyscf#1692, pyscf#48)

* Lint issues

* Update copyright line (issue pyscf#1421)

* Add example for BCCD (fix pyscf#1591)

NOTICE

@@ -1,4 +1,4 @@
-Copyright 2014-2021 The PySCF Developers.
+Copyright 2014-2024 The PySCF Developers.
 PySCF (www.pyscf.org) was developed by:
 
 Qiming Sun

examples/cc/02-bccd.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+
+'''
+Brueckner coupled-cluster doubles (BCCD) and BCCD(T) calculations for RHF, UHF and GHF.
+See also the relevant discussions in https://github.com/pyscf/pyscf/issues/1591 .
+'''
+
+from pyscf import gto, scf, cc
+
+mol = gto.M(
+    atom = 'H 0 0 0; F 0 0 1.1',
+    basis = 'ccpvdz',
+    verbose = 0,
+    spin = 0,
+)
+myhf = mol.HF().run()
+
+mycc = cc.BCCD(myhf).run()
+e_r = mycc.e_tot
+e_ccsd_t = mycc.ccsd_t()
+PSI4_reference = -0.002625521337000
+print(f'RHF-BCCD total energy {mycc.e_tot}.')
+print(f'BCCD(T) correlation energy {e_ccsd_t}. Difference to Psi4 {e_ccsd_t - PSI4_reference}')
+print(f'Max. value in BCCD T1 amplitudes {abs(mycc.t1).max()}')
+
+mygcc = cc.BCCD(myhf.to_ghf()).run()
+e_g = mygcc.e_tot
+print(f'GHF-BCCD total energy {mygcc.e_tot}.')
+
+# Run BCCD with frozen orbitals
+myhf = mol.UHF().run()
+myucc = cc.BCCD(myhf)
+myucc.frozen = [0]
+myucc.kernel()
+print(f'UHF-BCCD total energy {myucc.e_tot}.')

examples/dft/31-xc_value_on_grid.py

@@ -77,7 +77,6 @@
 aoL_value = mol.eval_gto('GTOval_spinor', coords)
 # Small components
 aoS_value = 1/(2*lib.param.LIGHT_SPEED) * mol.eval_gto('GTOval_sp_spinor', coords)
-# mL, mS are the spin-magentic moment at each point
 rho_m_L = r_numint.eval_rho(mol, aoL_value, dmLL)
 rho_m_S = r_numint.eval_rho(mol, aoS_value, dmSS)
 rho = rho_m_L[0] + rho_m_S[0]

examples/fci/13-wfn_symmetry.py

@@ -4,7 +4,11 @@
 #
 
 '''
-Assign FCI wavefunction symmetry
+Assign FCI wavefunction symmetry.
+
+This example demonstrates only the solver that utlizes the D2h (and its
+subgroup) symmetry. For molecules with cylindrical symmetry, a special solver is
+available, as shown in 03-cylindrical_symmetry.py
 '''
 
 import numpy
@@ -42,7 +46,7 @@
 ncore = mo_core.shape[1]
 nelec = mol.nelectron - ncore * 2
 fs = fci.direct_spin0_symm.FCI(mol)
-fs.wfnsym = 'A2'
+fs.wfnsym = 'A1'
 e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym, verbose=5)
 print('Energy of %s state %.12f' % (fs.wfnsym, e))
 
@@ -61,7 +65,7 @@
 h2 = ao2mo.kernel(mol, mo_cas)
 orbsym = scf.hf_symm.get_orbsym(mol, m.mo_coeff)[cas_idx]
 fs = fci.direct_spin0_symm.FCI(mol)
-fs.wfnsym = 'A2'
+fs.wfnsym = 'A1'
 e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym)
 print('Energy of %s state %.12f' % (fs.wfnsym, e))
 
@@ -79,7 +83,7 @@
 h2 = mc.get_h2eff(mo)
 orbsym = scf.hf_symm.get_orbsym(mol, m.mo_coeff)[cas_idx]
 fs = fci.direct_spin0_symm.FCI(mol)
-fs.wfnsym = 'A2'
+fs.wfnsym = 'A1'
 e, c = fs.kernel(h1, h2, norb, nelec, ecore=ecore, orbsym=orbsym)
 print('Energy of %s state %.12f' % (fs.wfnsym, e))
 

examples/scf/21-x2c.py

@@ -39,3 +39,9 @@
 mf.with_x2c = False
 energy = mf.kernel()
 print('E = %.12f, ref = %.12f' % (energy, scf.UKS(mol).kernel()))
+
+# Note, applying the x2c method for GHF/GKS classes results in a calculation
+# with the two-component relativistic Hamiltonian (including the spin-orbit
+# coupling effects). Although solving in spherical GTO bases, the results in
+# theory are equivalent to the X2C methods computed in the spinor bases.
+# See relevant examples in examples/x2c/03-x2c_ghf.py

pyscf/cc/__init__.py

@@ -216,3 +216,21 @@ def _finalize(self):
         return self
     mycc._finalize = _finalize.__get__(mycc, mycc.__class__)
     return mycc
+
+def BCCD(mf, frozen=None, u=None, conv_tol_normu=1e-5, max_cycle=20, diis=True,
+         canonicalization=True):
+    from pyscf.cc.bccd import bccd_kernel_
+    from pyscf.lib import StreamObject
+    mycc = CCSD(mf, frozen=frozen)
+
+    class BCCD(mycc.__class__):
+        def kernel(self):
+            obj = self.view(mycc.__class__)
+            obj.conv_tol = 1e-3
+            obj.kernel()
+            bccd_kernel_(obj, u, conv_tol_normu, max_cycle, diis,
+                         canonicalization, self.verbose)
+            self.__dict__.update(obj.__dict__)
+            return self.e_tot
+
+    return mycc.view(BCCD)

pyscf/cc/bccd.py

@@ -313,49 +313,3 @@ def A2u(A):
         mycc.t2 = t2
 
     return mycc
-
-if __name__ == "__main__":
-    import pyscf
-    from pyscf import cc
-
-    np.set_printoptions(3, linewidth=1000, suppress=True)
-    mol = pyscf.M(
-        atom = 'H 0 0 0; F 0 0 1.1',
-        basis = 'ccpvdz',
-        verbose = 4,
-        spin = 0,
-    )
-
-    myhf = mol.HF()
-    myhf.kernel()
-    E_ref = myhf.e_tot
-    rdm1_mf = myhf.make_rdm1()
-
-    mycc = cc.CCSD(myhf, frozen=None)
-    #mycc.frozen = [0]
-    mycc.kernel()
-    mycc.conv_tol = 1e-3
-
-    mycc = bccd_kernel_(mycc, diis=True, verbose=4)
-    e_r = mycc.e_tot
-    e_ccsd_t = mycc.ccsd_t()
-    # PSI4 reference
-    assert abs(e_ccsd_t - -0.002625521337000) < 1e-5
-
-    print (la.norm(mycc.t1))
-    assert la.norm(mycc.t1) < 1e-5
-
-    myhf = mol.UHF()
-    myhf.kernel()
-    myucc = cc.CCSD(myhf, frozen=None)
-    myucc.frozen = [0]
-    myucc.kernel()
-
-    mybcc = bccd_kernel_(myucc)
-    e_u = mybcc.e_tot
-
-    mygcc = cc.addons.convert_to_gccsd(mycc)
-    mybcc = bccd_kernel_(mygcc)
-    e_g = mybcc.e_tot
-    print (abs(e_g - e_r))
-    assert abs(e_g - e_r) < 1e-6

pyscf/df/incore.py

@@ -38,22 +38,34 @@
 
 
 def aux_e2(mol, auxmol_or_auxbasis, intor='int3c2e', aosym='s1', comp=None, out=None,
-           cintopt=None):
+           cintopt=None, shls_slice=None):
     '''3-center AO integrals (ij|L), where L is the auxiliary basis.
 
     Kwargs:
-        cintopt : Libcint-3.14 and newer version support to compute int3c2e
-            without the opt for the 3rd index.  It can be precomputed to
-            reduce the overhead of cintopt initialization repeatedly.
+        cintopt :
+            Precomputing certain pair-shell data. It can be created by
 
             cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, 'int3c2e')
+
+        shls_slice : 6-element tuple
+            Label the start-stop shells for each index in the integral tensor.
+            For the (ij|aux) = intor('int3c2e'), the tuple should be given as
+            (ish_start, ish_end, jsh_start, jsh_end, aux_start, aux_end)
     '''
     if isinstance(auxmol_or_auxbasis, gto.MoleBase):
         auxmol = auxmol_or_auxbasis
     else:
         auxbasis = auxmol_or_auxbasis
         auxmol = addons.make_auxmol(mol, auxbasis)
-    shls_slice = (0, mol.nbas, 0, mol.nbas, mol.nbas, mol.nbas+auxmol.nbas)
+    if shls_slice is None:
+        shls_slice = (0, mol.nbas, 0, mol.nbas,
+                      mol.nbas, mol.nbas+auxmol.nbas)
+    else:
+        assert len(shls_slice) == 6
+        assert shls_slice[5] < auxmol.nbas
+        shls_slice = list(shls_slice)
+        shls_slice[4] += mol.nbas
+        shls_slice[5] += mol.nbas
 
     # Extract the call of the two lines below
     #  pmol = gto.mole.conc_mol(mol, auxmol)