From 8ed07f54db291e1e8e26fee1d3d6b88c4c955759 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Fri, 21 Jun 2024 11:20:38 -0700 Subject: [PATCH] Update the default aug_etb settings --- examples/df/01-auxbasis.py | 5 ++ pyscf/df/__init__.py | 3 +- pyscf/df/addons.py | 87 +++++++++----------------- pyscf/df/autoaux.py | 50 ++++++++------- pyscf/df/grad/casscf.py | 41 ------------ pyscf/df/test/test_addons.py | 18 ++++-- pyscf/pbc/df/__init__.py | 2 +- pyscf/pbc/df/test/test_gdf_builder.py | 2 +- pyscf/pbc/df/test/test_mdf.py | 4 +- pyscf/pbc/df/test/test_mdf_builder.py | 2 +- pyscf/pbc/df/test/test_rsdf_builder.py | 2 +- 11 files changed, 85 insertions(+), 131 deletions(-) diff --git a/examples/df/01-auxbasis.py b/examples/df/01-auxbasis.py index b28665ac61..4fa410366e 100644 --- a/examples/df/01-auxbasis.py +++ b/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() diff --git a/pyscf/df/__init__.py b/pyscf/df/__init__.py index 98cf9aee66..dacf73dd76 100644 --- a/pyscf/df/__init__.py +++ b/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 diff --git a/pyscf/df/addons.py b/pyscf/df/addons.py index bbbb2fc2be..91c5a48a3c 100644 --- a/pyscf/df/addons.py +++ b/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) diff --git a/pyscf/df/autoaux.py b/pyscf/df/autoaux.py index e43d55aefc..018b58190e 100644 --- a/pyscf/df/autoaux.py +++ b/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 diff --git a/pyscf/df/grad/casscf.py b/pyscf/df/grad/casscf.py index 243a100f83..a27b713602 100644 --- a/pyscf/df/grad/casscf.py +++ b/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) diff --git a/pyscf/df/test/test_addons.py b/pyscf/df/test/test_addons.py index 35007a56b9..ef0ce2ef42 100644 --- a/pyscf/df/test/test_addons.py +++ b/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): diff --git a/pyscf/pbc/df/__init__.py b/pyscf/pbc/df/__init__.py index 39c86b03fb..9c0aeeff4b 100644 --- a/pyscf/pbc/df/__init__.py +++ b/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 diff --git a/pyscf/pbc/df/test/test_gdf_builder.py b/pyscf/pbc/df/test/test_gdf_builder.py index f73633784f..8e8e836ae8 100644 --- a/pyscf/pbc/df/test/test_gdf_builder.py +++ b/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 diff --git a/pyscf/pbc/df/test/test_mdf.py b/pyscf/pbc/df/test/test_mdf.py index 7a332a3d8f..97a72a394f 100644 --- a/pyscf/pbc/df/test/test_mdf.py +++ b/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) diff --git a/pyscf/pbc/df/test/test_mdf_builder.py b/pyscf/pbc/df/test/test_mdf_builder.py index 45300403d4..ca519c587b 100644 --- a/pyscf/pbc/df/test/test_mdf_builder.py +++ b/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 diff --git a/pyscf/pbc/df/test/test_rsdf_builder.py b/pyscf/pbc/df/test/test_rsdf_builder.py index 8eaa43b61e..feaa1b4a71 100644 --- a/pyscf/pbc/df/test/test_rsdf_builder.py +++ b/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