1616# Author: Qiming Sun <[email protected] > 1717#
1818
19+
1920import sys
2021import numpy
2122from pyscf .lib import logger
2223from pyscf import gto
2324from pyscf import ao2mo
2425from pyscf .data import elements
2526from pyscf .lib .exceptions import BasisNotFoundError
27+ from pyscf .df .autoaux import autoaux , autoabs
2628from pyscf import __config__
2729
2830DFBASIS = getattr (__config__ , 'df_addons_aug_etb_beta' , 'weigend' )
2931ETB_BETA = getattr (__config__ , 'df_addons_aug_dfbasis' , 2.0 )
3032FIRST_ETB_ELEMENT = getattr (__config__ , 'df_addons_aug_start_at' , 36 ) # 'Rb'
3133
34+ # TODO: Switch to other default scheme for auxiliary basis generation.
35+ # The auxiliary basis set generated by version 2.6 (and earlier) lacks compact
36+ # functions. It may cause higher errors in ERI integrals.
37+ USE_VERSION_26_AUXBASIS = True
38+
3239# Obtained from http://www.psicode.org/psi4manual/master/basissets_byfamily.html
3340DEFAULT_AUXBASIS = {
3441 # AO basis JK-fit MP2-fit
@@ -72,6 +79,58 @@ class load(ao2mo.load):
7279 def __init__ (self , eri , dataname = 'j3c' ):
7380 ao2mo .load .__init__ (self , eri , dataname )
7481
82+ def _aug_etb_element (nuc_charge , basis , beta ):
83+ l_max = max (b [0 ] for b in basis )
84+ emin_by_l = [1e99 ] * (l_max + 1 )
85+ emax_by_l = [0 ] * (l_max + 1 )
86+ for b in basis :
87+ l = b [0 ]
88+ if isinstance (b [1 ], int ):
89+ e_c = numpy .array (b [2 :])
90+ else :
91+ e_c = numpy .array (b [1 :])
92+ es = e_c [:,0 ]
93+ cs = e_c [:,1 :]
94+ es = es [abs (cs ).max (axis = 1 ) > 1e-3 ]
95+ emax_by_l [l ] = max (es .max (), emax_by_l [l ])
96+ emin_by_l [l ] = min (es .min (), emin_by_l [l ])
97+
98+ conf = elements .CONFIGURATION [nuc_charge ]
99+ # 1: H - Be, 2: B - Ca, 3: Sc - La, 4: Ce -
100+ max_shells = 4 - conf .count (0 )
101+
102+ if USE_VERSION_26_AUXBASIS :
103+ # This is the method that version 2.6 (and earlier) generates auxiliary
104+ # basis. It estimates the exponents ranges by geometric average.
105+ # This method is not recommended because it tends to generate diffuse
106+ # functions. Important compact functions might be improperly excluded.
107+ l_max = min (l_max , max_shells )
108+ l_max_aux = l_max * 2
109+ l_max1 = l_max + 1
110+ emin_by_l = numpy .array (emin_by_l [:l_max1 ])
111+ emax_by_l = numpy .array (emax_by_l [:l_max1 ])
112+ emax = (emax_by_l [:,None ] * emax_by_l ) ** .5 * 2
113+ emin = (emin_by_l [:,None ] * emin_by_l ) ** .5 * 2
114+ else :
115+ # Using normal average, more auxiliary functions, especially compact
116+ # functions, will be generated.
117+ l_max_aux = min (l_max , max_shells ) * 2
118+ l_max1 = l_max + 1
119+ emin_by_l = numpy .array (emin_by_l )
120+ emax_by_l = numpy .array (emax_by_l )
121+ emax = emax_by_l [:,None ] + emax_by_l
122+ emin = emin_by_l [:,None ] + emin_by_l
123+
124+ liljsum = numpy .arange (l_max1 )[:,None ] + numpy .arange (l_max1 )
125+ emax_by_l = numpy .array ([emax [liljsum == ll ].max () for ll in range (l_max_aux + 1 )])
126+ emin_by_l = numpy .array ([emin [liljsum == ll ].min () for ll in range (l_max_aux + 1 )])
127+
128+ ns = numpy .log ((emax_by_l + emin_by_l )/ emin_by_l ) / numpy .log (beta )
129+ etb = []
130+ for l , n in enumerate (numpy .ceil (ns ).astype (int )):
131+ if n > 0 :
132+ etb .append ((l , n , emin_by_l [l ], beta ))
133+ return etb
75134
76135def aug_etb_for_dfbasis (mol , dfbasis = DFBASIS , beta = ETB_BETA ,
77136 start_at = FIRST_ETB_ELEMENT ):
@@ -86,50 +145,14 @@ def aug_etb_for_dfbasis(mol, dfbasis=DFBASIS, beta=ETB_BETA,
86145 nuc_charge = gto .charge (symb )
87146 if nuc_charge < nuc_start :
88147 newbasis [symb ] = dfbasis
89- #?elif symb in mol._ecp:
90148 else :
91- conf = elements .CONFIGURATION [nuc_charge ]
92- max_shells = 4 - conf .count (0 )
93- emin_by_l = [1e99 ] * 8
94- emax_by_l = [0 ] * 8
95- l_max = 0
96- for b in mol ._basis [symb ]:
97- l = b [0 ]
98- l_max = max (l_max , l )
99- if l >= max_shells + 1 :
100- continue
101-
102- if isinstance (b [1 ], int ):
103- e_c = numpy .array (b [2 :])
104- else :
105- e_c = numpy .array (b [1 :])
106- es = e_c [:,0 ]
107- cs = e_c [:,1 :]
108- es = es [abs (cs ).max (axis = 1 ) > 1e-3 ]
109- emax_by_l [l ] = max (es .max (), emax_by_l [l ])
110- emin_by_l [l ] = min (es .min (), emin_by_l [l ])
111-
112- l_max1 = l_max + 1
113- emin_by_l = numpy .array (emin_by_l [:l_max1 ])
114- emax_by_l = numpy .array (emax_by_l [:l_max1 ])
115-
116- # Estimate the exponents ranges by geometric average
117- emax = numpy .sqrt (numpy .einsum ('i,j->ij' , emax_by_l , emax_by_l ))
118- emin = numpy .sqrt (numpy .einsum ('i,j->ij' , emin_by_l , emin_by_l ))
119- liljsum = numpy .arange (l_max1 )[:,None ] + numpy .arange (l_max1 )
120- emax_by_l = [emax [liljsum == ll ].max () for ll in range (l_max1 * 2 - 1 )]
121- emin_by_l = [emin [liljsum == ll ].min () for ll in range (l_max1 * 2 - 1 )]
122- # Tune emin and emax
123- emin_by_l = numpy .array (emin_by_l ) * 2 # *2 for alpha+alpha on same center
124- emax_by_l = numpy .array (emax_by_l ) * 2 #/ (numpy.arange(l_max1*2-1)*.5+1)
125-
126- ns = numpy .log ((emax_by_l + emin_by_l )/ emin_by_l ) / numpy .log (beta )
127- etb = []
128- for l , n in enumerate (numpy .ceil (ns ).astype (int )):
129- if n > 0 :
130- etb .append ((l , n , emin_by_l [l ], beta ))
149+ basis = mol ._basis [symb ]
150+ etb = _aug_etb_element (nuc_charge , basis , beta )
131151 if etb :
132152 newbasis [symb ] = gto .expand_etbs (etb )
153+ for l , n , emin , beta in etb :
154+ logger .info (mol , 'l = %d, exps = %s * %g^n for n = 0..%d' ,
155+ l , emin , beta , n - 1 )
133156 else :
134157 raise RuntimeError (f'Failed to generate even-tempered auxbasis for { symb } )
135158
@@ -182,12 +205,13 @@ def make_auxbasis(mol, mp2fit=False):
182205 auxbasis .update (auxdefault )
183206 aux_etb = set (auxbasis ) - set (auxdefault )
184207 if aux_etb :
185- logger .info (mol , 'Even tempered Gaussians are generated as '
208+ logger .warn (mol , 'Even tempered Gaussians are generated as '
186209 'DF auxbasis for %s' , ' ' .join (aux_etb ))
187210 for k in aux_etb :
188211 logger .debug (mol , ' ETB auxbasis for %s %s' , k , auxbasis [k ])
189212 return auxbasis
190213
214+ # TODO: add auxbasis keyword etb and auto
191215def make_auxmol (mol , auxbasis = None ):
192216 '''Generate a fake Mole object which uses the density fitting auxbasis as
193217 the basis sets. If auxbasis is not specified, the optimized auxiliary fitting
@@ -198,26 +222,30 @@ def make_auxmol(mol, auxbasis=None):
198222 even-tempered Gaussian basis set will be generated.
199223
200224 See also the paper JCTC, 13, 554 about generating auxiliary fitting basis.
225+
226+ Kwargs:
227+ auxbasis : str, list, tuple
228+ Similar to the input of orbital basis in Mole object.
201229 '''
202230 pmol = mol .copy (deep = False )
203231
204232 if auxbasis is None :
205233 auxbasis = make_auxbasis (mol )
206- elif isinstance (auxbasis , str ) and '+etb' in auxbasis :
207- dfbasis = auxbasis [:- 4 ]
208- auxbasis = aug_etb_for_dfbasis (mol , dfbasis )
209234 pmol .basis = auxbasis
210235
211236 if isinstance (auxbasis , (str , list , tuple )):
212237 uniq_atoms = {a [0 ] for a in mol ._atom }
213238 _basis = {a : auxbasis for a in uniq_atoms }
214- elif 'default' in auxbasis :
215- uniq_atoms = {a [0 ] for a in mol ._atom }
216- _basis = {a : auxbasis ['default' ] for a in uniq_atoms }
217- _basis .update (auxbasis )
218- del (_basis ['default' ])
219239 else :
220- _basis = auxbasis
240+ assert isinstance (auxbasis , dict )
241+ if 'default' in auxbasis :
242+ uniq_atoms = {a [0 ] for a in mol ._atom }
243+ _basis = {a : auxbasis ['default' ] for a in uniq_atoms }
244+ _basis .update (auxbasis )
245+ del (_basis ['default' ])
246+ else :
247+ _basis = auxbasis
248+
221249 pmol ._basis = pmol .format_basis (_basis )
222250
223251 # Note: To pass parameters like gauge origin, rsh-omega to auxmol,
0 commit comments