Numerical instabilities in the computation of pseudopotential matrix elements

Steep high angular momentum Gaussian basis functions in the vicinity of a nucleus whose inner electrons are replaced by an effective core potential may lead to numerical instabilities when calculating matrix elements of the core potential. Numerical roundoff errors may be amplified to an extent that spoils any result obtained in such a calculation. Effective core potential matrix elements for a model problem are computed with high numerical accuracy using the standard algorithm used in quantum chemical codes and compared to results of the MOLPRO program. Thus, it is demonstrated how the relative and absolute errors depend an basis function angular momenta, basis function exponents and the distance between the off-center basis function and the center carrying the effective core potential. Then, the problem is analyzed and closed expressions are derived for the expected numerical error in the limit of large basis function exponents. It is briefly discussed how other algorithms would behave in the critical case, and they are found to have problems as well. The numerical stability could be increased a little bit if the type 1 matrix elements were computed without making use of a partial wave expansion.