From f4b109e4e79a7fed5a12c14c685ad67fe8a99172 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Fri, 23 Aug 2024 14:27:35 -0700 Subject: [PATCH] Make autoaux match bse implementation --- pyscf/df/__init__.py | 2 +- pyscf/df/addons.py | 2 +- pyscf/df/autoaux.py | 120 +++++++++++++++++++++++++---------- pyscf/df/test/test_addons.py | 40 ++++++++++-- pyscf/gto/basis/bse.py | 83 +++++++++++++++++++----- 5 files changed, 190 insertions(+), 57 deletions(-) diff --git a/pyscf/df/__init__.py b/pyscf/df/__init__.py index dacf73dd76..9bf20700f3 100644 --- a/pyscf/df/__init__.py +++ b/pyscf/df/__init__.py @@ -34,7 +34,7 @@ from . import incore from . import outcore from . import addons -from .addons import (load, aug_etb, auto_aux, +from .addons import (load, aug_etb, autoaux, autoabs, DEFAULT_AUXBASIS, make_auxbasis, make_auxmol) from .df import DF, GDF, DF4C, GDF4C diff --git a/pyscf/df/addons.py b/pyscf/df/addons.py index 91c5a48a3c..3a22528009 100644 --- a/pyscf/df/addons.py +++ b/pyscf/df/addons.py @@ -24,7 +24,7 @@ from pyscf import ao2mo from pyscf.data import elements from pyscf.lib.exceptions import BasisNotFoundError -from pyscf.df.autoaux import auto_aux +from pyscf.df.autoaux import autoaux, autoabs from pyscf import __config__ DFBASIS = getattr(__config__, 'df_addons_aug_etb_beta', 'weigend') diff --git a/pyscf/df/autoaux.py b/pyscf/df/autoaux.py index 018b58190e..1fb7215bbe 100644 --- a/pyscf/df/autoaux.py +++ b/pyscf/df/autoaux.py @@ -25,7 +25,7 @@ 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]) +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]) BETA_SMALL = 1.8 @@ -45,18 +45,18 @@ def _primitive_emin_emax(basis): emax_by_l[l] = max(es.max(), emax_by_l[l]) emin_by_l[l] = min(es.min(), emin_by_l[l]) - 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]) + cs = e_c[:,1:] + cs = np.einsum('pi,p->pi', cs, gto.gto_norm(l, es)) # normalize GTOs + cs = gto.mole._nomalize_contracted_ao(l, es, cs) # prefactors in r_ints + 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) + # Eq (9) in the paper. r_exp**2 in the denominator may be a bug in bse + #e_eff = 2 * k**2 / (np.pi * r_exp) + e_eff = 2 * k**2 / (np.pi * r_exp**2) + e_eff_by_l[l] = max(e_eff.max(), 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): @@ -70,53 +70,107 @@ def _auto_aux_element(Z, basis, ecp_core=0): # TODO: handle ECP if Z <= 2: l_val = 0 - elif Z <= 18: + elif Z <= 20: l_val = 1 - elif Z <= 54: + elif Z <= 56: l_val = 2 else: l_val = 3 l_inc = 1 + if Z > 18: + l_inc = 2 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)] - a_max_by_l = [a_max_prim[liljsum==ll].max() for ll in range(l_max_aux+1)] - a_aux_by_l = [a_max_aux [liljsum==ll].max() for ll in range(l_max_aux+1)] + liljsub = abs(np.arange(l_max1)[:,None] - np.arange(l_max1)) + a_min_by_l = [a_min_prim[(liljsub<=ll) & (ll<=liljsum)].min() for ll in range(l_max_aux+1)] + a_max_by_l = [a_max_prim[(liljsub<=ll) & (ll<=liljsum)].max() for ll in range(l_max_aux+1)] + a_aux_by_l = [a_max_aux [(liljsub<=ll) & (ll<=liljsum)].max() for ll in range(l_max_aux+1)] - a_max_adjust = [max(F_LAUX[l] * a_aux_by_l[l], a_max_by_l[l]) + a_max_adjust = [min(F_LAUX[l] * a_aux_by_l[l], a_max_by_l[l]) for l in range(l_val*2+1)] 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 + emax = np.array(a_max_adjust) - ns_small = np.log(emax[:l_val*2+1] / emin[:l_val*2+1]) / np.log(BETA_SMALL) + ns_small = np.log(a_max_adjust[:l_val*2+1] / emin[:l_val*2+1]) / np.log(BETA_SMALL) etb = [] - for l, n in enumerate(np.ceil(ns_small).astype(int)): + # ns_small+1 to ensure the largest exponent in etb > emax + for l, n in enumerate(np.ceil(ns_small).astype(int) + 1): if n > 0: etb.append((l, n, emin[l], BETA_SMALL)) 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:l_max_aux+1])) - for l, n in enumerate(np.ceil(ns_big).astype(int)): + for l, n in enumerate(np.ceil(ns_big).astype(int) + 1): if n > 0: l = l + l_val*2+1 beta = BETA_BIG[l] etb.append((l, n, emin[l], beta)) return etb -def auto_aux(mol): - uniq_atoms = {a[0] for a in mol._atom} - newbasis = {} - for symb in uniq_atoms: +def autoaux(mol): + ''' + Create an auxiliary basis set for the given orbital basis set using + the Auto-Aux algorithm. + + See also: G. L. Stoychev, A. A. Auer, and F. Neese + Automatic Generation of Auxiliary Basis Sets + J. Chem. Theory Comput. 13, 554 (2017) + http://doi.org/10.1021/acs.jctc.6b01041 + ''' + from pyscf.gto.basis import bse + + def expand(symb): Z = gto.charge(symb) - basis = mol._basis[symb] - etb = _auto_aux_element(Z, basis) + etb = _auto_aux_element(Z, mol._basis[symb]) if etb: - newbasis[symb] = gto.expand_etbs(etb) - else: - raise RuntimeError(f'Failed to generate even-tempered auxbasis for {symb}') + return gto.expand_etbs(etb) + raise RuntimeError(f'Failed to generate even-tempered auxbasis for {symb}') + + uniq_atoms = {a[0] for a in mol._atom} + if isinstance(mol.basis, str): + try: + elements = [gto.charge(symb) for symb in uniq_atoms] + newbasis = bse.autoaux(mol.basis, elements) + except KeyError: + newbasis = {symb: expand(symb) for symb in uniq_atoms} + else: + newbasis = {} + for symb in uniq_atoms: + if symb in mol.basis and isinstance(mol.basis[symb], str): + try: + auxbs = bse.autoaux(mol.basis[symb], gto.charge(symb)) + newbasis[symb] = next(iter(auxbs.values())) + except KeyError: + newbasis[symb] = expand(symb) + else: + newbasis[symb] = expand(symb) + return newbasis + +def autoabs(mol): + ''' + Create a Coulomb fitting basis set for the given orbital basis set. + See also: + R. Yang, A. P. Rendell, and M. J. Frisch + Automatically generated Coulomb fitting basis sets: Design and accuracy for systems containing H to Kr + J. Chem. Phys. 127, 074102 (2007) + http://doi.org/10.1063/1.2752807 + ''' + from pyscf.gto.basis import bse + uniq_atoms = {a[0] for a in mol._atom} + if isinstance(mol.basis, str): + elements = [gto.charge(symb) for symb in uniq_atoms] + newbasis = bse.autoabs(mol.basis, elements) + else: + newbasis = {} + for symb in uniq_atoms: + Z = gto.charge(symb) + if symb in mol.basis and isinstance(mol.basis[symb], str): + auxbs = bse.autoabs(mol.basis[symb], gto.charge(symb)) + newbasis[symb] = next(iter(auxbs.values())) + else: + raise NotImplementedError return newbasis diff --git a/pyscf/df/test/test_addons.py b/pyscf/df/test/test_addons.py index ef0ce2ef42..32046f1b4c 100644 --- a/pyscf/df/test/test_addons.py +++ b/pyscf/df/test/test_addons.py @@ -16,13 +16,14 @@ # import unittest +import itertools import tempfile -import numpy +import numpy as np from pyscf import lib from pyscf import gto from pyscf import scf from pyscf import df - +from pyscf.gto.basis import bse class KnownValues(unittest.TestCase): def test_aug_etb(self): @@ -37,8 +38,8 @@ def test_aug_etb(self): self.assertEqual(len(auxbasis['O']), 42) self.assertEqual(len(auxbasis['H']), 13) - auxbasis = df.addons.auto_aux(mol) - self.assertEqual(len(auxbasis['O']), 54) + auxbasis = df.addons.autoaux(mol) + self.assertEqual(len(auxbasis['O']), 37) self.assertEqual(len(auxbasis['H']), 13) mol = gto.M( @@ -52,9 +53,9 @@ def test_aug_etb(self): 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) + auxbasis = df.addons.autoaux(mol) + self.assertEqual(len(auxbasis['O']), 47) + self.assertEqual(len(auxbasis['H']), 20) def test_make_auxbasis(self): mol = gto.M( @@ -88,6 +89,31 @@ def test_default_auxbasis(self): self.assertTrue(auxbasis['O'] == 'cc-pvdz-jkfit') self.assertTrue(isinstance(auxbasis['He'], list)) + @unittest.skipIf(bse.basis_set_exchange is None, "BSE library not installed.") + def test_auto_aux(self): + from pyscf.df.autoaux import autoaux, _auto_aux_element + ref = flatten(bse.autoaux('STO-3G', 'K')['K']) + mol = gto.M(atom='K', basis='STO-3G', spin=1) + etb = _auto_aux_element(mol.atom_charge(0), mol._basis['K']) + dat = flatten(gto.expand_etbs(etb)) + dat = np.array([float(f'{x:.6e}') for x in dat]) + self.assertAlmostEqual(abs(np.array(ref) - dat).max(), 0, 6) + + for key in ['H', 'Li', 'C', 'N', 'O', 'F', 'Si']: + ref = flatten(bse.autoaux('cc-pVTZ', key)[key]) + basis = bse.get_basis('cc-pVTZ', key)[key] + etb = _auto_aux_element(gto.charge(key), basis) + dat = flatten(gto.expand_etbs(etb)) + dat = np.array([float(f'{x:.6e}') for x in dat]) + self.assertAlmostEqual(abs(np.array(ref) - dat).max(), 0, 6) + +def flatten(lst): + if not isinstance(lst, list): + return [lst] + elif len(lst) == 0: + return lst + return list(itertools.chain(*[flatten(l) for l in lst])) + if __name__ == "__main__": print("Full Tests for df.addons") diff --git a/pyscf/gto/basis/bse.py b/pyscf/gto/basis/bse.py index 061d838a88..0e64ecbf83 100644 --- a/pyscf/gto/basis/bse.py +++ b/pyscf/gto/basis/bse.py @@ -11,13 +11,16 @@ basis_set_exchange = None -def _orbital_basis(basis): +def _orbital_basis(basis, make_general=False): '''Extracts the orbital basis from the BSE format in PySCF format''' r = {} - basis = manip.make_general(basis, False, True) - basis = sort.sort_basis(basis, False) + if make_general: + basis = manip.make_general(basis, False, use_copy=True) + basis = sort.sort_basis(basis, use_copy=False) + else: + basis = sort.sort_basis(basis, use_copy=True) # Elements for which we have electron basis electron_elements = [k for k, v in basis['elements'].items() if 'electron_shells' in v] @@ -40,21 +43,29 @@ def _orbital_basis(basis): ncontr = len(coefficients) nprim = len(exponents) am = shell['angular_momentum'] - assert len(am) == 1 - - shell_data = [am[0]] - for iprim in range(nprim): - row = [float(coefficients[ic][iprim]) for ic in range(ncontr)] - row.insert(0, float(exponents[iprim])) - shell_data.append(row) - atom_shells.append(shell_data) + if len(am) == len(coefficients): # SP basis + for l, c in zip(am, coefficients): + shell_data = [l] + for iprim in range(nprim): + row = [float(exponents[iprim]), float(c[iprim])] + shell_data.append(row) + atom_shells.append(shell_data) + else: + assert len(am) == 1 + shell_data = [am[0]] + for iprim in range(nprim): + row = [float(coefficients[ic][iprim]) for ic in range(ncontr)] + row.insert(0, float(exponents[iprim])) + shell_data.append(row) + atom_shells.append(shell_data) r[sym] = atom_shells # Collect the literature references - for ref in data['references']: - for key in ref['reference_keys']: - if key not in reference_list: - reference_list.append(key) + if 'references' in data: + for ref in data['references']: + for key in ref['reference_keys']: + if key not in reference_list: + reference_list.append(key) return r, reference_list @@ -110,6 +121,48 @@ def _print_basis_information(basis): print('Last revised on {}'.format(revision_date)) print('Revision description: {}'.format(revision_description)) +def get_basis(name, elements): + ''' + Obtain a basis set + + Args: + name : str + Name of the basis set, case insensitive. + elements : str, int or list + + Returns: + A dict of basis set in PySCF internal basis format. + ''' + basis = basis_set_exchange.api.get_basis(name, elements) + return _orbital_basis(basis)[0] + +def autoaux(name, elements): + ''' + Create an auxiliary basis set for the given orbital basis set using + the Auto-Aux algorithm. + + See also: G. L. Stoychev, A. A. Auer, and F. Neese + Automatic Generation of Auxiliary Basis Sets + J. Chem. Theory Comput. 13, 554 (2017) + http://doi.org/10.1021/acs.jctc.6b01041 + ''' + obs = basis_set_exchange.api.get_basis(name, elements) + auxbs = manip.autoaux_basis(obs) + return _orbital_basis(auxbs)[0] + +def autoabs(name, elements): + ''' + Create a Coulomb fitting basis set for the given orbital basis set. + See also: + R. Yang, A. P. Rendell, and M. J. Frisch + Automatically generated Coulomb fitting basis sets: Design and accuracy for systems containing H to Kr + J. Chem. Phys. 127, 074102 (2007) + http://doi.org/10.1063/1.2752807 + ''' + obs = basis_set_exchange.api.get_basis(name, elements) + auxbs = manip.autoabs_basis(obs) + return _orbital_basis(auxbs)[0] + if __name__ == '__main__': from basis_set_exchange import api, references