Merge branch 'master' into tdhf-diagonalization

.github/workflows/ci_linux/build_pyscf.sh

@@ -3,7 +3,7 @@
 set -e
 
 cd ./pyscf/lib
-curl -L "https://github.com/pyscf/pyscf-build-deps/blob/master/pyscf-2.4a-deps.tar.gz?raw=true" | tar xzf -
+curl -L "https://github.com/pyscf/pyscf-build-deps/blob/master/pyscf-2.8a-deps.tar.gz?raw=true" | tar xzf -
 mkdir build; cd build
 cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=ON -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
 make -j4

CHANGELOG

@@ -366,6 +366,7 @@ PySCF 1.7.6 (2021-03-28)
   - Auxiliary second-order Green's function perturbation theory (AGF2)
   - Smearing for molecules
   - Visscher small component correction approximation for DHF
+  - DFT+U
 * Improved
   - Threading safety in Mole temporary context
   - Basis parser to support arithmetic expressions in basis data

examples/dft/16-dft_d3.py

@@ -6,18 +6,21 @@
 '''
 D3 and D4 Dispersion.
 
-This is a simplified dispersion interface to
-d3 (https://github.com/dftd3/simple-dftd3) and
-d4 (https://github.com/dftd4/dftd4) libraries.
-This interface can automatically configure the necessary settings including
-dispersion, xc, and nlc attributes of PySCF mean-field objects. However,
-advanced features of d3 and d4 program are not available.
+This example shows the functionality of the pyscf-dispersion extension, which
+implements a simplified interface to d3 (https://github.com/dftd3/simple-dftd3)
+and d4 (https://github.com/dftd4/dftd4) libraries. This interface can
+automatically configure the necessary settings including dispersion, xc, and nlc
+attributes of PySCF mean-field objects. However, advanced features of d3 and d4
+program are not available. To execute the code in this example, you would need
+to install the pyscf-dispersion extension:
+
+pip install pyscf-dispersion
 
 If you need to access more features d3 and d4 libraries, such as overwriting the
 dispersion parameters, you can use the wrapper provided by the simple-dftd3 and
-dftd4 libraries. When using these libraries, please disable the .disp attribute
-of the underlying mean-field object, and properly set the .xc and .nlc attributes
-following this example.
+dftd4 libraries. When accessing dispersion through the APIs of these libraries,
+please disable the .disp attribute of the mean-field object, and properly set
+the .xc and .nlc attributes as shown in this example.
 '''
 
 mol = pyscf.M(
@@ -134,7 +137,7 @@
 mf.disp = 'd4'
 mf.kernel() # Crash
 
-# Please note, not every xc functional can be used for D3/D4 corrections.
+# Please note, NOT every xc functional can be used for D3/D4 corrections.
 # For the valid xc and D3/D4 combinations, please refer to
 # https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml
 # https://github.com/dftd4/dftd4/blob/main/assets/parameters.toml

examples/geomopt/16-apply_soscf.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+
+'''
+Automatically apply SOSCF for unconverged SCF calculations in geometry optimization.
+'''
+
+import pyscf
+from pyscf.geomopt import geometric_solver, as_pyscf_method
+
+# Force SCF to stop early at an unconverged state
+pyscf.scf.hf.RHF.max_cycle = 5
+
+mol = pyscf.M(
+    atom='''
+    O 0 0 0
+    H 0 .7 .8
+    H 0 -.7 .8
+    ''')
+mf = mol.RHF()
+try:
+    mol_eq = geometric_solver.optimize(mf)
+except RuntimeError:
+    print('geometry optimization should stop for unconverged calculation')
+
+# Apply SOSCF when SCF is not converged
+mf_scanner = mf.as_scanner()
+def apply_soscf_as_needed(mol):
+    mf_scanner(mol)
+    if not mf_scanner.converged:
+        mf_soscf = mf_scanner.newton().run()
+        for key in ['converged', 'mo_energy', 'mo_coeff', 'mo_occ', 'e_tot', 'scf_summary']:
+            setattr(mf_scanner, key, getattr(mf_soscf, key))
+    grad = mf_scanner.Gradients().kernel()
+    return mf_scanner.e_tot, grad
+
+mol_eq = geometric_solver.optimize(as_pyscf_method(mol, apply_soscf_as_needed))

examples/pbc/22-dft+u.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+
+'''
+DFT+U with kpoint sampling
+'''
+
+from pyscf.pbc import gto, dft
+cell = gto.Cell()
+cell.unit = 'A'
+cell.atom = 'C 0.,  0.,  0.; C 0.8917,  0.8917,  0.8917'
+cell.a = '''0.      1.7834  1.7834
+            1.7834  0.      1.7834
+            1.7834  1.7834  0.    '''
+
+cell.basis = 'gth-dzvp'
+cell.pseudo = 'gth-pade'
+cell.verbose = 4
+cell.build()
+
+kmesh = [2, 2, 2]
+kpts = cell.make_kpts(kmesh)
+# Add U term to the 2p orbital of the second Carbon atom
+U_idx = ['1 C 2p']
+U_val = [5.0]
+mf = dft.KRKSpU(cell, kpts, U_idx=U_idx, U_val=U_val, minao_ref='gth-szv')
+print(mf.U_idx)
+print(mf.U_val)
+mf.run()
+
+# When space group symmetry in k-point samples is enabled, the symmetry adapted
+# DFT+U method will be invoked automatically.
+kpts = cell.make_kpts(
+    kmesh, wrap_around=True, space_group_symmetry=True, time_reversal_symmetry=True)
+# Add U term to 2s and 2p orbitals
+U_idx = ['2p', '2s']
+U_val = [5.0, 2.0]
+mf = dft.KUKSpU(cell, kpts, U_idx=U_idx, U_val=U_val, minao_ref='gth-szv')
+print(mf.U_idx)
+print(mf.U_val)
+mf.run()

pyscf/cc/ccsd.py

@@ -76,10 +76,16 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
         normt = numpy.linalg.norm(tmpvec)
         tmpvec = None
         if mycc.iterative_damping < 1.0:
-            alpha = mycc.iterative_damping
-            t1new = (1-alpha) * t1 + alpha * t1new
-            t2new *= alpha
-            t2new += (1-alpha) * t2
+            alpha = numpy.asarray(mycc.iterative_damping)
+            if isinstance(t1, tuple): # e.g. UCCSD
+                t1new = tuple((1-alpha) * numpy.asarray(t1_part) + alpha * numpy.asarray(t1new_part)
+                    for t1_part, t1new_part in zip(t1, t1new))
+                t2new = tuple((1-alpha) * numpy.asarray(t2_part) + alpha * numpy.asarray(t2new_part)
+                    for t2_part, t2new_part in zip(t2, t2new))
+            else:
+                t1new = (1-alpha) * numpy.asarray(t1) + alpha * numpy.asarray(t1new)
+                t2new *= alpha
+                t2new += (1-alpha) * numpy.asarray(t2)
         t1, t2 = t1new, t2new
         t1new = t2new = None
         t1, t2 = mycc.run_diis(t1, t2, istep, normt, eccsd-eold, adiis)

pyscf/cc/test/test_uccsd.py

@@ -727,6 +727,20 @@ def test_ao2mo(self):
         self.assertAlmostEqual(abs(eri_incore.OVvo - eri_outcore.OVvo).max(), 0, 12)
         self.assertAlmostEqual(abs(eri_incore.OVvv - eri_outcore.OVvv).max(), 0, 12)
 
+    def test_damping(self):
+        mol = gto.M(
+            atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
+            basis = 'ccpvdz',
+            spin=1,
+            charge=-1,
+            symmetry = True,
+            verbose = 0
+        )
+        mf = scf.UHF(mol).run()
+        mycc = cc.UCCSD(mf)
+        mycc.iterative_damping = 0.5
+        mycc.run()
+        self.assertAlmostEqual(mycc.e_tot, -100.07710261186985, 7)
 
 if __name__ == "__main__":
     print("Full Tests for UCCSD")

pyscf/df/addons.py

@@ -151,8 +151,8 @@ def aug_etb_for_dfbasis(mol, dfbasis=DFBASIS, beta=ETB_BETA,
             if etb:
                 newbasis[symb] = gto.expand_etbs(etb)
                 for l, n, emin, beta in etb:
-                    logger.info(mol, 'l = %d, exps = %s * %g^n for n = 0..%d',
-                                l, emin, beta, n-1)
+                    logger.info(mol, 'ETB for %s: l = %d, exps = %s * %g^n , n = 0..%d',
+                                symb, l, emin, beta, n-1)
             else:
                 raise RuntimeError(f'Failed to generate even-tempered auxbasis for {symb}')
 

pyscf/df/autoaux.py

@@ -24,6 +24,7 @@
 from math import factorial
 import numpy as np
 from pyscf import gto
+from pyscf.lib import logger
 
 F_LAUX   = np.array([20 , 7.0, 4.0, 4.0, 3.5, 2.5, 2.0, 2.0])
 BETA_BIG = np.array([1.8, 2.0, 2.2, 2.2, 2.2, 2.3, 3.0, 3.0])
@@ -136,6 +137,9 @@ def expand(symb):
         Z = gto.charge(symb)
         etb = _auto_aux_element(Z, mol._basis[symb])
         if etb:
+            for l, n, emin, beta in etb:
+                logger.info(mol, 'ETB for %s: l = %d, exps = %s * %g^n , n = 0..%d',
+                            symb, l, emin, beta, n-1)
             return gto.expand_etbs(etb)
         raise RuntimeError(f'Failed to generate even-tempered auxbasis for {symb}')
 

pyscf/df/df_jk.py

@@ -236,7 +236,7 @@ def to_gpu(self):
         return lib.to_gpu(self, obj)
 
 
-def get_jk(dfobj, dm, hermi=1, with_j=True, with_k=True, direct_scf_tol=1e-13):
+def get_jk(dfobj, dm, hermi=0, with_j=True, with_k=True, direct_scf_tol=1e-13):
     assert (with_j or with_k)
     if (not with_k and not dfobj.mol.incore_anyway and
         # 3-center integral tensor is not initialized
@@ -266,8 +266,7 @@ def get_jk(dfobj, dm, hermi=1, with_j=True, with_k=True, direct_scf_tol=1e-13):
     if not with_k:
         for eri1 in dfobj.loop():
             # uses numpy.matmul
-            vj += (dmtril @ eri1.T) @ eri1
-
+            vj += dmtril.dot(eri1.T).dot(eri1)
 
     elif getattr(dm, 'mo_coeff', None) is not None:
         #TODO: test whether dm.mo_coeff matching dm
@@ -297,7 +296,8 @@ def get_jk(dfobj, dm, hermi=1, with_j=True, with_k=True, direct_scf_tol=1e-13):
             assert (nao_pair == nao*(nao+1)//2)
             if with_j:
                 # uses numpy.matmul
-                vj += (dmtril @ eri1.T) @ eri1
+                vj += dmtril.dot(eri1.T).dot(eri1)
+
             for k in range(nset):
                 nocc = orbo[k].shape[1]
                 if nocc > 0:
@@ -324,8 +324,7 @@ def get_jk(dfobj, dm, hermi=1, with_j=True, with_k=True, direct_scf_tol=1e-13):
             naux, nao_pair = eri1.shape
             if with_j:
                 # uses numpy.matmul
-                vj += (dmtril @ eri1.T) @ eri1
-
+                vj += dmtril.dot(eri1.T).dot(eri1)
 
             for k in range(nset):
                 buf1 = buf[0,:naux]
@@ -344,7 +343,7 @@ def get_jk(dfobj, dm, hermi=1, with_j=True, with_k=True, direct_scf_tol=1e-13):
     logger.timer(dfobj, 'df vj and vk', *t0)
     return vj, vk
 
-def get_j(dfobj, dm, hermi=1, direct_scf_tol=1e-13):
+def get_j(dfobj, dm, hermi=0, direct_scf_tol=1e-13):
     from pyscf.scf import _vhf
     from pyscf.scf import jk
     from pyscf.df import addons
@@ -437,7 +436,7 @@ def get_j(dfobj, dm, hermi=1, direct_scf_tol=1e-13):
     return numpy.asarray(vj).reshape(dm_shape)
 
 
-def r_get_jk(dfobj, dms, hermi=1, with_j=True, with_k=True):
+def r_get_jk(dfobj, dms, hermi=0, with_j=True, with_k=True):
     '''Relativistic density fitting JK'''
     t0 = (logger.process_clock(), logger.perf_counter())
     mol = dfobj.mol

pyscf/dft/numint.py

@@ -176,9 +176,9 @@ def eval_rho(mol, ao, dm, non0tab=None, xctype='LDA', hermi=0,
         if hermi:
             rho[1:4] *= 2  # *2 for + einsum('pi,ij,pj->p', ao[i], dm, ao[0])
         else:
+            c1 = _dot_ao_dm(mol, ao[0], dm.conj().T, non0tab, shls_slice, ao_loc)
             for i in range(1, 4):
-                c1 = _dot_ao_dm(mol, ao[i], dm, non0tab, shls_slice, ao_loc)
-                rho[i] += _contract_rho(c1, ao[0])
+                rho[i] += _contract_rho(c1, ao[i])
     else: # meta-GGA
         if with_lapl:
             # rho[4] = \nabla^2 rho, rho[5] = 1/2 |nabla f|^2

pyscf/dft/radi.py

@@ -47,14 +47,21 @@
 def murray(n, *args, **kwargs):
     raise RuntimeError('Not implemented')
 
-# Gauss-Chebyshev of the first kind,  and the transformed interval [0,\infty)
 def becke(n, charge, *args, **kwargs):
     '''Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033'''
     if charge == 1:
         rm = BRAGG_RADII[charge]
     else:
         rm = BRAGG_RADII[charge] * .5
-    t, w = numpy.polynomial.chebyshev.chebgauss(n)
+
+    # Points and weights for Gauss-Chebyshev quadrature points of the second kind
+    # The weights are adjusted to integrate a function on the interval [-1, 1] with uniform weighting
+    # instead of weighted by sqrt(1 - t^2) = sin(i pi / (n+1)).
+    i = numpy.arange(n) + 1
+    t = numpy.cos(i * numpy.pi / (n + 1))
+    w = numpy.pi / (n + 1) * numpy.sin(i * numpy.pi / (n + 1))
+
+    # Change of variables to map the domain to [0, inf)
     r = (1+t)/(1-t) * rm
     w *= 2/(1-t)**2 * rm
     return r, w

pyscf/dft/rks.py

@@ -322,6 +322,9 @@ class KohnShamDFT:
 
     _keys = {'xc', 'nlc', 'grids', 'disp', 'nlcgrids', 'small_rho_cutoff'}
 
+    # Use rho to filter grids
+    small_rho_cutoff = getattr(__config__, 'dft_rks_RKS_small_rho_cutoff', 1e-7)
+
     def __init__(self, xc='LDA,VWN'):
         # By default, self.nlc = '' and self.disp = None
         self.xc = xc
@@ -333,9 +336,6 @@ def __init__(self, xc='LDA,VWN'):
         self.nlcgrids = gen_grid.Grids(self.mol)
         self.nlcgrids.level = getattr(
             __config__, 'dft_rks_RKS_nlcgrids_level', self.nlcgrids.level)
-        # Use rho to filter grids
-        self.small_rho_cutoff = getattr(
-            __config__, 'dft_rks_RKS_small_rho_cutoff', 1e-7)
 ##################################################
 # don't modify the following attributes, they are not input options
         self._numint = numint.NumInt()

pyscf/dft/test/test_grids.py

@@ -91,7 +91,7 @@ def test_radi(self):
 
         grid.radi_method = radi.becke
         grid.build(with_non0tab=False)
-        self.assertAlmostEqual(lib.fp(grid.weights), 818061.875131255, 7)
+        self.assertAlmostEqual(lib.fp(grid.weights), 780.7183109298, 9)
 
     def test_prune(self):
         grid = gen_grid.Grids(h2o)

pyscf/dft/test/test_numint.py

@@ -177,10 +177,20 @@ def test_dot_ao_ao_high_cost(self):
         ao = dft.numint.eval_ao(mol, mf.grids.coords, deriv=1)
         nao = ao.shape[1]
         ao_loc = mol.ao_loc_nr()
+        cutoff = mf.grids.cutoff * 1e2
+        cutoff2 = numint.CUTOFF * 1e2
+        nbins = numint.NBINS * 2 - int(numint.NBINS * numpy.log(cutoff)
+                                       / numpy.log(mf.grids.cutoff))
+        pair_mask = mol.get_overlap_cond() < -numpy.log(cutoff2)
+        wv = numpy.ones(ao.shape[1])
         res0 = lib.dot(ao[0].T, ao[1])
         res1 = dft.numint._dot_ao_ao(mol, ao[0], ao[1], non0tab,
                                      shls_slice=(0,mol.nbas), ao_loc=ao_loc)
+        res2 = dft.numint._dot_ao_ao_sparse(ao[0], ao[1], wv, nbins,
+                                            non0tab, pair_mask, ao_loc,
+                                            hermi=0)
         self.assertAlmostEqual(abs(res0 - res1).max(), 0, 9)
+        self.assertAlmostEqual(abs(res0 - res2).max(), 0, 9)
 
     def test_eval_rho(self):
         numpy.random.seed(10)

pyscf/gto/basis/parse_cp2k.py

@@ -53,12 +53,14 @@ def parse(string, symb=None, optimize=False):
     '''
     if symb is not None:
         raw_data = list(filter(None, re.split(BASIS_SET_DELIMITER, string)))
-        string = _search_basis_block(raw_data, symb)
-        if not string:
+        line_data = _search_basis_block(raw_data, symb)
+        if not line_data:
             raise BasisNotFoundError(f'Basis not found for {symb}')
+    else:
+        line_data = string.splitlines()
 
     bastxt = []
-    for dat in string.splitlines():
+    for dat in line_data:
         x = dat.split('#')[0].strip()
         if (x and not x.startswith('END') and not x.startswith('BASIS')):
             bastxt.append(x)
@@ -122,8 +124,7 @@ def _parse(blines, optimize=False):
 def search_seg(basisfile, symb):
     with open(basisfile, 'r') as fin:
         fdata = re.split(BASIS_SET_DELIMITER, fin.read())
-    raw_basis = _search_basis_block(fdata[1:], symb)
-    if not raw_basis:
+    line_data = _search_basis_block(fdata[1:], symb)
+    if not line_data:
         raise BasisNotFoundError(f'Basis for {symb} not found in {basisfile}')
-    return [x.strip() for x in raw_basis.splitlines()
-            if x.strip() and 'END' not in x]
+    return [x for x in line_data if x and 'END' not in x]