Update the default aug_etb settings

examples/df/01-auxbasis.py

@@ -65,3 +65,8 @@
 mf.with_df.auxbasis = df.aug_etb(mol, beta=1.7)
 mf.kernel()
 
+# Generating auxiliary basis using the AutoAux algorithms proposed by Stoychev
+# (JCTC, 13, 554)
+mf = scf.RHF(mol).density_fit()
+mf.with_df.auxbasis = df.auto_aux(mol)
+mf.kernel()

pyscf/df/__init__.py

@@ -34,7 +34,8 @@
 from . import incore
 from . import outcore
 from . import addons
-from .addons import load, aug_etb, DEFAULT_AUXBASIS, make_auxbasis, make_auxmol
+from .addons import (load, aug_etb, auto_aux,
+                     DEFAULT_AUXBASIS, make_auxbasis, make_auxmol)
 from .df import DF, GDF, DF4C, GDF4C
 
 from . import r_incore

pyscf/df/addons.py

@@ -24,15 +24,13 @@
 from pyscf import ao2mo
 from pyscf.data import elements
 from pyscf.lib.exceptions import BasisNotFoundError
+from pyscf.df.autoaux import auto_aux
 from pyscf import __config__
 
 DFBASIS = getattr(__config__, 'df_addons_aug_etb_beta', 'weigend')
 ETB_BETA = getattr(__config__, 'df_addons_aug_dfbasis', 2.0)
 FIRST_ETB_ELEMENT = getattr(__config__, 'df_addons_aug_start_at', 36)  # 'Rb'
 
-LVAL_SCHEME = 'PySCF'  # or 'Gaussian', 'ORCA'
-L_INC = 1
-
 # Obtained from http://www.psicode.org/psi4manual/master/basissets_byfamily.html
 DEFAULT_AUXBASIS = {
     # AO basis       JK-fit                     MP2-fit
@@ -77,36 +75,11 @@ def __init__(self, eri, dataname='j3c'):
         ao2mo.load.__init__(self, eri, dataname)
 
 def _aug_etb_element(nuc_charge, basis, beta):
-    # lval
-    if LVAL_SCHEME == 'Gaussian':
-        # Yang, JCP, 127, 074102
-        # 0: H - Be, 1: B - Ca, 2: Sc - Ba, 3: La -
-        if nuc_charge <= 2:
-            max_shells = 0
-        elif nuc_charge <= 18:
-            max_shells = 1
-        elif nuc_charge <= 54:
-            max_shells = 2
-        else:
-            max_shells = 3
-    elif LVAL_SCHEME == 'ORCA':
-        # Stoychev, JCTC, 13, 554
-        # 0: H - Be, 1: B - Ca, 2: Sc - Ba, 3: La -
-        conf = elements.NRSRHFS_CONFIGURATION[nuc_charge]
-        max_shells = max(3 - conf.count(0), 0)
-    else:
-        # PySCF original settings
-        # 1: H - Be, 2: B - Ca, 3: Sc - La, 4: Ce -
-        conf = elements.CONFIGURATION[nuc_charge]
-        max_shells = 4 - conf.count(0)
-
-    emin_by_l = [1e99] * 8
-    emax_by_l = [0] * 8
-    l_max = 0
+    l_max = max(b[0] for b in basis)
+    emin_by_l = [1e99] * (l_max+1)
+    emax_by_l = [0] * (l_max+1)
     for b in basis:
         l = b[0]
-        l_max = max(l_max, l)
-
         if isinstance(b[1], int):
             e_c = numpy.array(b[2:])
         else:
@@ -117,21 +90,35 @@ def _aug_etb_element(nuc_charge, basis, beta):
         emax_by_l[l] = max(es.max(), emax_by_l[l])
         emin_by_l[l] = min(es.min(), emin_by_l[l])
 
-    l_max1 = l_max + 1
-    emin_by_l = numpy.array(emin_by_l[:l_max1])
-    emax_by_l = numpy.array(emax_by_l[:l_max1])
+    conf = elements.CONFIGURATION[nuc_charge]
+    # 1: H - Be, 2: B - Ca, 3: Sc - La, 4: Ce -
+    max_shells = 4 - conf.count(0)
+
+    if 0:
+        # This is the method that PySCF-2.6 (and older versions) generates
+        # auxiliary basis. It estimates the exponents ranges by geometric average.
+        # This method is not recommended because it tends to generate diffused
+        # functions. Important compact functions might be improperly excluded.
+        l_max = min(l_max, max_shells)
+        l_max_aux = l_max * 2
+        l_max1 = l_max + 1
+        emin_by_l = numpy.array(emin_by_l[:l_max1])
+        emax_by_l = numpy.array(emax_by_l[:l_max1])
+        emax = (emax_by_l[:,None] * emax_by_l) ** .5 * 2
+        emin = (emin_by_l[:,None] * emin_by_l) ** .5 * 2
+    else:
+        # Using normal average, more auxiliary functions, especially compact
+        # functions, will be generated.
+        l_max_aux = min(l_max, max_shells) * 2
+        l_max1 = l_max + 1
+        emin_by_l = numpy.array(emin_by_l)
+        emax_by_l = numpy.array(emax_by_l)
+        emax = emax_by_l[:,None] + emax_by_l
+        emin = emin_by_l[:,None] + emin_by_l
 
-    # Estimate the exponents ranges by geometric average
-    emax = numpy.sqrt(numpy.einsum('i,j->ij', emax_by_l, emax_by_l))
-    emin = numpy.sqrt(numpy.einsum('i,j->ij', emin_by_l, emin_by_l))
     liljsum = numpy.arange(l_max1)[:,None] + numpy.arange(l_max1)
-    # Setting l_max_aux as JCTC, 13, 554
-    l_max_aux = max(max_shells*2, l_max + L_INC)
-    emax_by_l = [emax[liljsum==ll].max() for ll in range(l_max_aux+1)]
-    emin_by_l = [emin[liljsum==ll].min() for ll in range(l_max_aux+1)]
-    # Tune emin and emax
-    emin_by_l = numpy.array(emin_by_l) * 2
-    emax_by_l = numpy.array(emax_by_l) * 2
+    emax_by_l = numpy.array([emax[liljsum==ll].max() for ll in range(l_max_aux+1)])
+    emin_by_l = numpy.array([emin[liljsum==ll].min() for ll in range(l_max_aux+1)])
 
     ns = numpy.log((emax_by_l+emin_by_l)/emin_by_l) / numpy.log(beta)
     etb = []
@@ -236,18 +223,11 @@ def make_auxmol(mol, auxbasis=None):
 
     if auxbasis is None:
         auxbasis = make_auxbasis(mol)
-    elif auxbasis == 'etb':
-        auxbasis = aug_etb_for_dfbasis(mol, start_at=0)
-    elif '+etb' in auxbasis:
-        dfbasis = auxbasis[:-4]
-        auxbasis = aug_etb_for_dfbasis(mol, dfbasis)
     pmol.basis = auxbasis
 
     if isinstance(auxbasis, (str, list, tuple)):
         uniq_atoms = {a[0] for a in mol._atom}
         _basis = {a: auxbasis for a in uniq_atoms}
-        uniq_atoms = set([a[0] for a in mol._atom])
-        _basis = dict([(a, auxbasis) for a in uniq_atoms])
     else:
         assert isinstance(auxbasis, dict)
         if 'default' in auxbasis:
@@ -257,11 +237,6 @@ def make_auxmol(mol, auxbasis=None):
             del (_basis['default'])
         else:
             _basis = auxbasis
-        if 'etb' in _basis.values():
-            etb_abs = aug_etb_for_dfbasis(mol, start_at=0)
-            for elem, bas in _basis.items():
-                if bas == 'etb':
-                    _basis[elem] = etb_abs[elem]
 
     pmol._basis = pmol.format_basis(_basis)
 

pyscf/df/autoaux.py

@@ -22,19 +22,18 @@
 '''
 
 from math import factorial
-import np as np
+import numpy as np
 from pyscf import gto
 
 F_LAUX = np.array([20, 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])
 BETA_SMALL = 1.8
-L_INC = 1
 
 def _primitive_emin_emax(basis):
     l_max = max(b[0] for b in basis)
     emin_by_l = [1e99] * (l_max+1)
     emax_by_l = [0] * (l_max+1)
-    a_eff_by_l = [0] * (l_max+1)
+    e_eff_by_l = [0] * (l_max+1)
 
     for b in basis:
         l = b[0]
@@ -43,25 +42,32 @@ def _primitive_emin_emax(basis):
         else:
             e_c = np.array(b[1:])
         es = e_c[:,0]
-        cs = e_c[:,1:]
-        cs = np.einsum('pi,p->pi', cs, gto.gto_norm(l, es))
         emax_by_l[l] = max(es.max(), emax_by_l[l])
         emin_by_l[l] = min(es.min(), emin_by_l[l])
 
-        ee = es[:,None] + es
-        r_ints = gto.gaussian_int(l*2+3, ee)
-        r_exp = np.einsum('pi,pq,qi->i', cs, r_ints, cs)
-
-        k = 2**(2*l+1) * (factorial(l+1))**2 / factorial(2*l+2)
-        a_eff_by_l[l] = max(2 * k**2 / (np.pi * r_exp.min()),
-                            a_eff_by_l)
-    return np.array(emax_by_l), np.array(emin_by_l), np.array(a_eff_by_l)
-
-def _auto_aux_element(Z, basis):
+        if 0 and es.size == 1:
+            e_eff_by_l[l] = max(e_eff_by_l[l], emax_by_l[l])
+        else:
+            cs = e_c[:,1:]
+            cs = np.einsum('pi,p->pi', cs, gto.gto_norm(l, es))
+            ee = es[:,None] + es
+            r_ints = gto.gaussian_int(l*2+3, ee) # \int \chi^2 r dr
+            r_exp = np.einsum('pi,pq,qi->i', cs, r_ints, cs)
+
+            k = 2**(2*l+1) * factorial(l+1)**2 / factorial(2*l+2)
+            e_eff = 2 * k**2 / (np.pi * r_exp)
+            e_eff_by_l[l] = max(e_eff.min(), e_eff_by_l[l])
+    return np.array(emax_by_l), np.array(emin_by_l), np.array(e_eff_by_l)
+
+def _auto_aux_element(Z, basis, ecp_core=0):
     a_max_by_l, a_min_by_l, a_eff_by_l = _primitive_emin_emax(basis)
-    l_max1 = a_max_by_l.size
+    a_min_prim = a_min_by_l[:,None] + a_min_by_l
+    a_max_prim = a_max_by_l[:,None] + a_max_by_l
+    a_max_aux = a_eff_by_l[:,None] + a_eff_by_l
 
+    l_max1 = a_max_by_l.size
     l_max = l_max1 - 1
+    # TODO: handle ECP
     if Z <= 2:
         l_val = 0
     elif Z <= 18:
@@ -70,11 +76,8 @@ def _auto_aux_element(Z, basis):
         l_val = 2
     else:
         l_val = 3
-    l_max_aux = max(l_val*2, l_max+L_INC)
-
-    a_min_prim = a_min_by_l[:,None] + a_min_by_l
-    a_max_prim = a_max_by_l[:,None] + a_max_by_l
-    a_max_aux = a_eff_by_l[:,None] + a_eff_by_l
+    l_inc = 1
+    l_max_aux = min(max(l_val*2, l_max+l_inc), l_max*2)
 
     liljsum = np.arange(l_max1)[:,None] + np.arange(l_max1)
     a_min_by_l = [a_min_prim[liljsum==ll].min() for ll in range(l_max_aux+1)]
@@ -86,6 +89,7 @@ def _auto_aux_element(Z, basis):
     a_max_adjust = a_max_adjust + a_aux_by_l[l_val*2+1 : l_max_aux+1]
 
     emin = np.array(a_min_by_l)
+    # To ensure emax > emin even if there is only one GTO function
     emax = np.array(a_max_adjust) + 1e-3
 
     ns_small = np.log(emax[:l_val*2+1] / emin[:l_val*2+1]) / np.log(BETA_SMALL)
@@ -96,7 +100,7 @@ def _auto_aux_element(Z, basis):
 
     if l_max_aux > l_val*2:
         ns_big = (np.log(emax[l_val*2+1:] / emin[l_val*2+1:])
-                  / np.log(BETA_BIG[l_val*2+1:]))
+                  / np.log(BETA_BIG[l_val*2+1:l_max_aux+1]))
         for l, n in enumerate(np.ceil(ns_big).astype(int)):
             if n > 0:
                 l = l + l_val*2+1
@@ -112,7 +116,7 @@ def auto_aux(mol):
         basis = mol._basis[symb]
         etb = _auto_aux_element(Z, basis)
         if etb:
-            newbasis[symb] = etb
+            newbasis[symb] = gto.expand_etbs(etb)
         else:
             raise RuntimeError(f'Failed to generate even-tempered auxbasis for {symb}')
     return newbasis

pyscf/df/grad/casscf.py

@@ -227,44 +227,3 @@ def _finalize(self):
     to_gpu = lib.to_gpu
 
 Grad = Gradients
-
-#from pyscf import mcscf
-#mcscf.mc1step.CASSCF.Gradients = lib.class_as_method(Gradients)
-
-
-if __name__ == '__main__':
-    from pyscf import scf
-    from pyscf import mcscf
-    from pyscf import df
-    #from pyscf.grad import numeric
-
-    mol = gto.Mole()
-    mol.atom = 'N 0 0 0; N 0 0 1.2; H 1 1 0; H 1 1 1.2'
-    mol.basis = '631g'
-    mol.build()
-    aux = df.aug_etb (mol)
-    mf = scf.RHF(mol).density_fit (auxbasis=aux).run()
-    mc = mcscf.CASSCF(mf, 4, 4).run()
-    mc.conv_tol = 1e-10
-    de = Gradients (mc).kernel()
-    #de_num = numeric.Gradients (mc).kernel ()
-    #print(lib.finger(de) - 0.019602220578635747)
-    #print(lib.finger(de) - lib.finger (de_num))
-
-    mol = gto.Mole()
-    mol.verbose = 0
-    mol.atom = 'N 0 0 0; N 0 0 1.2'
-    mol.basis = 'sto3g'
-    mol.build()
-    mf = scf.RHF(mol).density_fit (auxbasis=aux).run()
-    mc = mcscf.CASSCF(mf, 4, 4)
-    mc.conv_tol = 1e-10
-    mc.kernel ()
-    de = Gradients (mc).kernel()
-
-    mcs = mc.as_scanner()
-    mol.set_geom_('N 0 0 0; N 0 0 1.201')
-    e1 = mcs(mol)
-    mol.set_geom_('N 0 0 0; N 0 0 1.199')
-    e2 = mcs(mol)
-    print(de[1,2], (e1-e2)/0.002*lib.param.BOHR)

pyscf/df/test/test_addons.py

@@ -33,7 +33,13 @@ def test_aug_etb(self):
                       1     0    0.757    0.587''',
             basis = 'cc-pvdz',
         )
-        df.addons.aug_etb(mol)
+        auxbasis = df.addons.aug_etb(mol)
+        self.assertEqual(len(auxbasis['O']), 42)
+        self.assertEqual(len(auxbasis['H']), 13)
+
+        auxbasis = df.addons.auto_aux(mol)
+        self.assertEqual(len(auxbasis['O']), 54)
+        self.assertEqual(len(auxbasis['H']), 13)
 
         mol = gto.M(
             verbose = 0,
@@ -43,8 +49,12 @@ def test_aug_etb(self):
             basis = ('cc-pvdz', [[4, (1., 1.)]])
         )
         auxbasis = df.addons.aug_etb(mol)
-        self.assertEqual(len(auxbasis['O']), 36)
-        self.assertEqual(len(auxbasis['H']), 12)
+        self.assertEqual(len(auxbasis['O']), 59)
+        self.assertEqual(len(auxbasis['H']), 16)
+
+        auxbasis = df.addons.auto_aux(mol)
+        self.assertEqual(len(auxbasis['O']), 58)
+        self.assertEqual(len(auxbasis['H']), 16)
 
     def test_make_auxbasis(self):
         mol = gto.M(
@@ -69,7 +79,7 @@ def test_make_auxbasis(self):
                      'O': ('631g', [[0, 0, (1., 1.)]])}
         )
         auxbasis = df.addons.make_auxbasis(mol)
-        self.assertEqual(len(auxbasis['O']), 32)
+        self.assertEqual(len(auxbasis['O']), 35)
         self.assertEqual(len(auxbasis['H']), 3)
 
     def test_default_auxbasis(self):

pyscf/pbc/df/__init__.py

@@ -23,7 +23,7 @@
 from .mdf import MDF
 from .aft import AFTDF
 from .fft import FFTDF
-from pyscf.df.addons import aug_etb
+from pyscf.df.addons import aug_etb, auto_aux
 from .incore import make_auxcell
 
 # For backward compatibility

pyscf/pbc/df/test/test_gdf_builder.py

@@ -19,7 +19,7 @@
 from pyscf import lib
 import pyscf.pbc
 from pyscf.pbc import gto as pgto
-from pyscf.pbc.df import df, aug_etb, FFTDF
+from pyscf.pbc.df import df, FFTDF
 from pyscf.pbc.df import gdf_builder
 from pyscf.pbc.df import ft_ao
 from pyscf.pbc.tools import pbc as pbctools

pyscf/pbc/df/test/test_mdf.py

@@ -145,8 +145,8 @@ def test_get_eri_gamma_1(self):
         odf.mesh = [11]*3
         eri0000 = odf.get_eri()
         self.assertTrue(eri0000.dtype == numpy.double)
-        self.assertAlmostEqual(eri0000.real.sum(), 27.27275005755319, 7)
-        self.assertAlmostEqual(lib.fp(eri0000), 1.0613327424886785, 8)
+        self.assertAlmostEqual(eri0000.real.sum(), 27.27280884567336, 7)
+        self.assertAlmostEqual(lib.fp(eri0000), 1.0613205275015447, 8)
 
         eri1111 = kmdf1.get_eri((kpts[0],kpts[0],kpts[0],kpts[0]))
         self.assertTrue(eri1111.dtype == numpy.double)

pyscf/pbc/df/test/test_mdf_builder.py

@@ -20,7 +20,7 @@
 from pyscf import ao2mo
 import pyscf.pbc
 from pyscf.pbc import gto as pgto
-from pyscf.pbc.df import df, aug_etb, FFTDF
+from pyscf.pbc.df import df, FFTDF
 from pyscf.pbc.df import mdf
 from pyscf.pbc.df import ft_ao
 from pyscf.pbc.tools import pbc as pbctools

pyscf/pbc/df/test/test_rsdf_builder.py

@@ -19,7 +19,7 @@
 from pyscf import lib
 import pyscf.pbc
 from pyscf.pbc import gto as pgto
-from pyscf.pbc.df import df, aug_etb, FFTDF, make_auxcell
+from pyscf.pbc.df import df, FFTDF, make_auxcell
 from pyscf.pbc.df import rsdf_builder
 from pyscf.pbc.df import ft_ao
 from pyscf.pbc.tools import pbc as pbctools