Shape and Energy Consistent Pseudopotentials for Correlated Electron systems

TL;DR

This paper introduces energy consistent correlated-electron pseudopotentials (eCEPPs) that improve the accuracy of molecular geometries and dissociation energies in correlated-electron calculations, outperforming previous pseudopotentials.

Contribution

The paper develops a novel method combining shape and energy consistency for generating pseudopotentials tailored for correlated-electron systems, with demonstrated accuracy improvements.

Findings

eCEPPs significantly improve molecular geometry predictions.

Dissociation energy errors are an order of magnitude smaller than previous pseudopotentials.

Validated with coupled cluster calculations on various elements.

Abstract

A method is developed for generating pseudopotentials for use in correlated-electron calculations. The paradigms of shape and energy consistency are combined and defined in terms of correlated-electron wave-functions. The resulting energy consistent correlated electron pseudopotentials (eCEPPs) are constructed for H, Li--F, Sc--Fe, and Cu. Their accuracy is quantified by comparing the relaxed molecular geometries and dissociation energies they provide with all electron results, with all quantities evaluated using coupled cluster singles doubles and triples calculations. Errors inherent in the pseudopotentials are also compared with those arising from a number of approximations commonly used with pseudopotentials. The eCEPPs provide a significant improvement in optimised geometries and dissociation energies for small molecules, with errors for the latter being an order-of-magnitude smaller than for Hartree-Fock-based pseudopotentials available in the literature. Gaussian basis sets are optimised for use with these pseudopotentials.

Tables

Table 1. Table 1: Summary of pseudopotential errors for TNDF pseudopotentials, CEPPs, and eCEPPs. Values are for the mean, absolute mean, and maximum error taken over the test set, denoted Δ ​ x i ¯ ¯ Δ subscript 𝑥 𝑖 \overline{\Delta x_{i}} , | Δ ​ x i | ¯ ¯ Δ subscript 𝑥 𝑖 \overline{|\Delta x_{i}|} , and Max [ Δ ​ x i ] delimited-[] Δ subscript 𝑥 𝑖 [\Delta x_{i}] for geometry parameters, and Δ ​ D e ¯ ¯ Δ subscript 𝐷 𝑒 \overline{\Delta D_{e}} , | Δ ​ D e | ¯ ¯ Δ subscript 𝐷 𝑒 \overline{|\Delta D_{e}|} , and Max [ Δ ​ D e ] delimited-[] Δ subscript 𝐷 𝑒 [\Delta D_{e}] for dissociation energies. The test set of CEPPs does not include LiB, LiC, LiF, CuC, CuO, or CuF as these molecules were not included in the previous work. The test set for the TNDF pseudopotentials excludes the same molecules, as well as H 2 O 2 , CrF, and FeF due to lack of convergence.

	($\times {10}^{-3}$ Å)			(meV)
Pseudopotential	$\overline{\Delta x_i}$	$\overline{\left|\Delta x_i\right|}$	Max$\left[\Delta x_i\right]$	$\overline{\Delta D_e}$	$\overline{\left|\Delta D_e\right|}$	Max$\left[\Delta D_e\right]$
TNDF	$1.2$	$4.3$	$31.1$	$13.65$	$117.05$	$431.70$
CEPP	$-3.3$	$4.1$	$14.1$	$-8.89$	$30.72$	$99.11$
eCEPP	$-0.6$	$1.2$	$-7.5$	$1.67$	$11.04$	$42.02$

Table 2. Table 2: Bond lengths ( r e subscript 𝑟 𝑒 r_{e} ) and dissociation energies ( D e subscript 𝐷 𝑒 D_{e} ) for monoxides of Sc–Fe and Cu. Results are shown evaluated using CCSD(T) with all electrons present, and with eCEPPs. Experimental values from the literature are shown for all oxides, and previously published composite results for TiO and CrO are evaluated by correcting CCSD(T) energies by including spin-orbit coupling and higher-order correlation effects (with geometries optimised using CCSD(T) only).

Method	ScO	TiO	VO	CrO	MnO	FeO	CuO
	$r_e$ (Å)
eCEPP	$1.669\left(6\right)$	$1.620\left(3\right)$	$1.584\left(5\right)$	$1.606\left(3\right)$	$1.633\left(5\right)$	$1.616\left(2\right)$	$1.730\left(2\right)$
AE	$1.669\left(7\right)$	$1.619\left(3\right)$	$1.582\left(6\right)$	$1.605\left(5\right)$	$1.633\left(5\right)$	$1.617\left(2\right)$	$1.736\left(3\right)$
Composite	$-$	$1.619$111 From [Pan_2017, ]. Bond lengths are optimised using CCSD(T), and dissociation energies are CCSD(T) corrected to include spin-orbit coupling and higher-order correlation.	$-$	$1.6185$222 From [Chan_2012, ]. As [Pan_2017, ].	$-$	$-$	$-$
Expt.333 Sum of experimental atomisation energies and zero-point vibrational energies. Data as cited by [Miliordos_SctoFe, ] for Sc–Fe, and [Harrison_2000, ] for Cu.	$1.6656$	$1.6203$	$1.589$	$1.6213$	$1.6477$	$1.6194$	$1.724$
	$D_e$ (eV)
eCEPP	$7.019\left(8\right)$	$7.097\left(10\right)$	$6.765\left(10\right)$	$4.540\left(5\right)$	$3.712\left(9\right)$	$3.804\left(3\right)$	$2.813\left(2\right)$
AE	$7.024\left(11\right)$	$7.080\left(11\right)$	$6.749\left(10\right)$	$4.521\left(4\right)$	$3.696\left(12\right)$	$3.791\left(5\right)$	$2.776\left(3\right)$
Composite	$-$	$6.916$111 From [Pan_2017, ]. Bond lengths are optimised using CCSD(T), and dissociation energies are CCSD(T) corrected to include spin-orbit coupling and higher-order correlation.	$-$	$4.660$222 From [Chan_2012, ]. As [Pan_2017, ].	$-$	$-$	$-$
Expt.333 Sum of experimental atomisation energies and zero-point vibrational energies. Data as cited by [Miliordos_SctoFe, ] for Sc–Fe, and [Harrison_2000, ] for Cu.	$6.981\left(9\right)$	$6.93\left(7\right)$	$6.55\left(9\right)$	$4.63\left(7\right)$	$3.88\left(8\right)$	$4.235$	$2.928$
				$4.80\left(9\right)$	$3.8\left(2\right)$

Equations

$${\hat{V}}^{\text{pp}}$$

$${\tilde{V}}_l^{\text{pp-cpp}}$$

$$V^{\text{pp}}_{l}=\left\{\begin{array}[]{rl}\tilde{V}^{\text{pp-cpp}}_{l}-Z_{v}/r+\tilde{V}^{\text{pp-cpp}}_{l_{max}}+V^{\text{cpp}}&l\neq l_{max}\\ -Z_{v}/r+\tilde{V}^{\text{pp-cpp}}_{l_{max}}+V^{\text{cpp}}&l=l_{max}.\end{array}\right.$$

$$V^{\text{cpp}}=-\frac{1}{2}{f}^2(r/r_0)\frac{\alpha }{r^4}.$$

$$A_{ql}$$

$$a_{ql}$$

$$P_{1}(\{\lambda^{A}_{l}\},\{\lambda^{a}_{ql}\},r_{0})=\sum_{l=0}^{l_{max}}\left(\left[\phi_{l}(r^{c}_{l})-\rho^{\frac{1}{2}}_{l}(r^{c}_{l})\right]^{2}+1-\left[\int_{r^{c}_{l}}^{\infty}\phi_{l}\rho^{\frac{1}{2}}_{l}dr\right]^{2}\Big{/}\left[\int_{r^{c}_{l}}^{\infty}\phi_{l}^{2}dr\int_{r^{c}_{l}}^{\infty}\rho_{l}dr\right]\right)$$

$$P_2(\{{\lambda }_l^A\},\{{\lambda }_{ql}^a\},r_0)=\frac{1}{{\beta }^2}\times \sum _{j=1}^N{\left[\left(E_j^{\text{AE}}-E_0^{\text{AE}}\right)-\left(E_j^{\text{pp}}-E_0^{\text{pp}}\right)\right]}^2$$

$$\hat{V}=V_{l_{max}}^{pp}+\sum _l\sum _{ij}|\delta V_l{\varphi }_{li}⟩B_{l,ij}⟨{\varphi }_{lj}\delta V_l|,$$

$$\Delta$$

$${\Delta }_{\text{pp}}$$

$${\Delta }_{\text{cpp}}$$

$${\Delta }_{\text{DFT}}$$

$${\Delta }_{\text{KB}}$$

Full text

I Introduction

Electronic structure methods in solid state theory and quantum chemistry provide a powerful hierarchy of ab initio tools for the accurate description of interacting atomic systems. There are fundamental limits to the accuracy at each level of theory. Pseudopotentials, or effective core potentials, provide an approximation that replaces the interaction of largely inactive core electrons and valence electrons with a potential, which reduces the number of particles required.Dolg_2012 Such an approximation is often necessary due to limitations on computational resources, and it also allows the elimination of analytic and numerical problems associated with singularities in electron-nuclear interaction potentials.

Methods for generating pseudopotentials within density functional theory (DFT) are well established. Errors associated with pseudopotentials are often systematically controlled, and reliably estimated to be small and secondary to the errors due to approximate exchange correlation functionals.Milman_2000 However, the accuracy of DFT is not sufficient for many systems. Methods that involve explicitly correlated wave functions, such as coupled clusterccsd (CC) or quantum Monte CarloCeperley-Alder_1980 , Foulkes_QMC_review , Casino_reference , Lester_review_2012 (QMC) theory, are often required.

The cost of coupled-cluster calculations with single, double, and perturbative triple excitations [CCSD(T)] scales with the seventh power of the number of electrons, or as $O(n^{3}N^{4})$,ccsd with $n$ the number of electrons and $N$ the number of basis functions. The use of pseudopotentials reduces both $n$ and $N$ by removing core electrons and requiring fewer basis functions than the nucleus-electron interaction.

The computational cost of fermion QMC calculations increases more slowly, as the third or fourth power of the number of particles. However, the scaling of the error with atomic number, $Z$, is $Z^{5}-Z^{6.5}$.Ceperley_1986 , Ma_2005 The use of pseudopotentials reduces the effective value of $Z$, making QMC calculations feasible for systems with heavy atoms.

At these higher levels of theory estimates of, and methods to control, pseudopotential errors are not usually available. Furthermore, available pseudopotentials often reproduce atomic properties at the DFT or HF level of theory. This introduces an uncontrolled error due to the inconsistent application of theory to the definition and application of the pseudopotentials.

In this study we generate pseudopotentials from explicitly correlated atomic wave-functions, with no recourse to reducing these many-body wave-functions to orbitals arising from a self-consistent potential. The influence of core electrons on valence electrons is recreated by consistently combining three distinct representations of the core-valence interaction, namely core scattering, core polarisation, and differences between ground state energies.

Scattering properties of core electrons are defined in terms of density matrices of atoms or ions with one valence electron. Long-range effects from dynamic core relaxation (arising from electron-electron correlation) can be included in the pseudopotential via a core polarisation potential (CPP).Mitroy_2010

Pseudopotentials defined to reproduce the one-body core-valence interaction and core polarisation are not unique. Consequently, the accuracy may be further improved by searching over all pseudopotentials that reproduce such properties, and seeking that which most accurately reproduces ionisation and excitation energies, and electron affinities for an isolated atom.

A combined reproduction of core scattering, core polarisation, and atomic excitation energies allows the generation of a new pseudopotential from correlated electron calculations, referred to as an energy consistent correlated electron pseudopotential (eCEPP).

Our procedure for constructing eCEPPs is an extension of the correlated electron pseudopotential (CEPP) method described in two previous papers, Trail_2013_pseudopotentials , Trail_2015_pseudopotentials which are summarised here. Only two of the criteria given above are used to generate CEPPs, the preservation of core scattering and the representation of core polarizability. An all-electron (AE) multi-determinant wave-function for an ion/atom with one valence electron provides an AE charge density. From this a single-electron charge density is constructed by removing an ab initio core electron charge density. Within a ‘core radius’ of the nuclei the single-electron charge density is replaced by a standard form that preserves the continuity of the value and first four derivatives at the core radius. Direct inversion of a one-electron Schrödinger equation then provides the pseudopotential. The analytic basis for this approach is an extension of norm-conservationHamann_pseudopotentials_1979 to the many-body Hamiltonian, which conserves the scattering properties of core electrons to first-order by reproducing the one-body density matrix outside of the core region.Acioli_1994

Two subtle issues arise in the process of constructing CEPPs. The multi-determinant wave-function is composed of orbitals that are not unique, and without associated energy eigenvalues. A separation of the density matrix into a core and valence components is achieved using natural orbitals,Davidson_1972 and defining the core in terms of the natural orbitals with the largest occupation numbers. While core natural orbitals have no direct effect on the pseudo wave-function in the core region they validate the core radius used and remove a small but significant influence of core electrons at the core radius.

The long-range, one-body part of the CPP naturally emerges in the process of generating a CEPP. This is represented using well-known parameterised CPP formsShirley_1993 , Muller_1984 such that the CEPP is the sum of an effective potential and the CPP, and the one-body part of the total potential is ab initio. The ion-ion and many-body part of the CPP remains semi-empirical, but with parameters that are consistent with the one-body component.

The CEPPs generated for H, Li–F, and Sc–Fe, together with CCSD(T), were shown to reproduce AE results for small molecules. Perhaps the strongest indicator that the CEPP provides a good starting point is that it is a pseudopotential generated from a highly ionic state which is accurate for a neutral molecule. This is in contrast to the well known poor transferability of KS-DFT pseudopotentials between different oxidation states, which is primarily caused by a pseudopotential definition that takes the self-consistent potential to be fixed and uses de-screening to remove the interaction between valence electrons. Neither of these approximations is used for the CEPPs. Iron provides a good example, with a CEPP generated from a Fe15+ ion resulting in properties for molecules containing iron which are as accurate as for molecules with few valence electrons such as Li2.

Before defining eCEPPs we note two alternative strategies for improving the accuracy of pseudopotentials that have not been investigated. The definition of the CEPPs could be extended to reproduce two-body core scattering properties. This would require the construction of a pseudo-Hamiltonian that reproduces both first and second-order density matrices outside the core region of an isolated ion/atom with two valence electrons. Such a pseudo-Hamiltonian would contain an extra two-body potential that describes the electron-electron interaction mediated by core electrons. This two-body interaction potential is approximately included in the CPP part of the CEPPs and is small, so it is unlikely that a more accurate description would significantly improve the pseudopotential accuracy.

Another, related, strategy would be to reproduce one-body core scattering properties to higher than first-order. An extension of norm-conserving KS-DFT pseudopotentials to reproduce higher order scattering propertiesShirley_1989 provides a modest reduction in transfer error, suggesting that this approach will not reliably provide a more transferable pseudopotential.

These observations suggest that an eCEPP that conserves the CEPP properties while enforcing further ‘energy consistency’ between a number of states may provide more accurate pseudopotentials. Such a combination of extrapolation and interpolation of atomic properties combines the strengths of shape-consistent and energy consistent ab initio pseudopotentials. Electron correlation can be included throughout, without the need for an independent electron approximation.

The article is organised as follows. Section II describes the basis and implementation of the combined energy and shape consistent strategy outlined above. The ability of eCEPPs to reproduce AE results in correlated electron calculations is analysed in Sec. III. The performance of the pseudopotentials for titanium are addressed in Sec. III.1, and errors in optimised geometries and dissociation energies for a test set of molecules composed of first row and $3d$ transition metal atoms are assessed in Sec. III.2. Section III.3 provides a comparison of errors arising from the eCEPPs with those from other approximations commonly used in ab initio calculations. Gaussian basis sets optimised for use with these potentials are described and validated in Sec. III.4.

II Shape and energy consistent correlated pseudopotentials

II.1 Theoretical basis

The pseudopotentials are formulated as the sum of a semi-local operator and a many-body potential,

[TABLE]

where each $V^{\text{pp-cpp}}_{l}$ and the $V^{\text{pp-cpp}}_{l_{max}}$ act on the spherical harmonic projection to channels $l<l_{max}$ and $l\geq l_{max}$, respectively. The potential $V^{\text{cpp}}$ is the many-body CPP.

The $V^{\text{pp-cpp}}_{l}$ terms are represented using the parameterised form

[TABLE]

For the local channel ($l=l_{max}$), $n_{ql}=-1$ for $q=1$, $n_{ql}=+1$ for $q=6$, and $n_{ql}=0$ otherwise. For other channels $n_{ql}=0$. The pseudo-atomic number is $Z_{v}=Z-Z_{c}$, where $Z_{c}$ is the number of core electrons represented by the pseudopotential. The local channel is $l_{max}=max\{l_{occ}\}$, where $\{l_{occ}\}$ is the set of all $l$ values occurring in the electronic configurations used to generate the pseudopotential.

In terms of the parameterised components of the pseudopotential, $\tilde{V}^{\text{pp-cpp}}_{l}$, the complete pseudopotential is

[TABLE]

The $V^{\text{cpp}}$ potential is composed of the sum of a one-body core-electron potential, a two-body electron-electron interaction, and two further terms involving the positions of both electrons and nuclei which are expressed in the form provided by Shirley and Martin.Shirley_1993 For a single atom with one valence electron only the local component of the CPP remains,

[TABLE]

The short range truncation function $f(x)=(1-e^{-x^{2}})^{2}$ removes the non-physical singularity at $r=0$, $r_{0}$ is a cutoff radii for this function, and $\alpha$ is the dipole polarizability of the core.

The previously available CEPPs are parameterised in the same form. An eCEPP candidate is expressed in terms of scaled CEPP parameters,

[TABLE]

and $r_{0}$ is varied without scaling, so that the candidate potentials are equal to the CEPP for scaling parameters $(\{\lambda^{A}_{l}\},\{\lambda^{a}_{ql}\})$ all equal to unity and $r_{0}=r_{0}^{CEPP}$.

Overcompleteness, and the accompanying emergence of undesirable oscillations in the potential, is prevented by the uniform scaling of $\{A_{ql}^{CEPP}\}$ parameters. The CPP cutoff radii, $r_{0}$, is taken to be a free parameter, while the core polarizability, $\alpha$ in Eq. (6), is taken from previously published highly accurate values available for all the atoms considered.Shirley_1993 , Patil_1985

The pseudopotential is required to be finite and have a first- and second derivative equal to zero at $r=0$. Taking these constraints into account results in $8+7l_{max}$ free parameters (for first row and $3d$ transition metal atoms $l_{max}=$2 and 3, resulting in 22 and 29 parameters, respectively).

In order to seek an eCEPP we start with all scaling parameters equal to unity. Scaling parameters are then varied to search for new values that preserve the CEPP charge density outside of the core region for single-valence electron states, and reproduce AE energy differences between atomic states of different symmetry and charge. The relaxation of the one-electron pseudo charge density within the core region provides the freedom for a candidate potential to reproduce the AE energy differences while preserving the core-valence interaction that defines the CEPP.

We define a penalty function that is zero for an eCEPP, and expressed as a sum of two terms, $P=P_{1}+P_{2}$. The first of these, $P_{1}$, measures the degree to which a candidate potential reproduces the one-electron density outside of the core region. This is given by

[TABLE]

where the sum is taken over the single-valence-electron ${}^{2}L$ states corresponding to $l=L$ for $l\leq l_{max}$, and the core radius for each channel, $r^{c}_{l}$, defines the sphere within which the one-electron pseudo charge density is allowed to vary.

In Eq. 9, $\rho_{l}$ is the AE electron density with the core electrons removed, and $\phi_{l}$ is the one-body wave-function corresponding to the candidate pseudopotential. For the initial parameter values, $\rho_{l}=\phi^{2}_{l}$ and $P_{1}=0$.

Equation (9) is constructed as the simplest form for which $P_{1}$ is non-negative, the lowest order variation of $P_{1}$ about the minimum is second order, and $\rho_{l}=\phi^{2}_{l}$ outside of the core region results in $P_{1}=0$.

A minimum with $P_{1}=0$ also occurs for parameter values that are not equal to those for the CEPP. It is these parameters that we seek, as they correspond to the many pseudopotentials for which the one-electron charge densities are equal to the AE valence charge density outside of the core region only.

The second component of the penalty function, $P_{2}$, measures the degree to which the candidate potential reproduces atomic excitation energies. A set of $N+1$ states is considered, which is generally composed of the neutral atom ground state, the anion ground state, and several ground and excited states for neutral atoms and ions. We index the set of states using $j$, with $j=0$ the neutral ground state, and define the second penalty function as

[TABLE]

where $E_{j}^{\text{AE}}$ and $E_{j}^{\text{pp}}$ are the total energies of the $j^{th}$ atomic state evaluated with all electrons present and the candidate pseudopotential, respectively. The prefactor $1/\beta^{2}$ sets the importance of the second penalty function relative to the first, with $\beta$ setting the accuracy with which the eCEPP is required to reproduce AE total energy differences. (No prefactor is used in $P_{1}$ as it is expected to be bounded by 0 and ${\sim}1$ for all parameter values.)

Minimising the full penalty function, $P(\{\lambda^{A}_{l}\},\{\lambda^{a}_{ql}\},r_{0})$, with respect to $(\{\lambda^{A}_{l}\},\{\lambda^{a}_{ql}\},r_{0})$ provides the eCEPP that preserves core scattering properties as well as reproducing the AE energy spectrum of an isolated atom.

II.2 Implementation

The evaluation of the two components of the penalty function, $P=P_{1}+P_{2}$, are distinct from each other.

The function $P_{1}$ includes correlation effects if the target charge density, $\rho_{l}$, is provided by a correlated-electron calculation. We take the target density outside of the core region to be that provided by the CEPP as, by definition, this is equal to the AE charge density provided by Multi-Configuration Hartree Fock (MCHF)Fischer_1997 , atsp2k calculations with the core contribution removed. Consequently, the target and candidate charge densities outside of the core are both obtained by direct solution of the one-body Schrödinger equation. The target charge density is provided using the CEPP, while the candidate charge density is provided by the candidate pseudopotential. Each term in $P_{1}$ is evaluated using summation and numerical integration.

The evaluation of $P_{2}$ is performed by summation of the squared differences between CCSD(T) energies, with energies evaluated using Gaussian basis sets and extrapolation to the complete basis set (CBS) limit.

All CCSD(T) energies are evaluated using the MOLPROmolpro code for both AE and candidate eCEPP atoms to provide the target and candidate energies, respectively. The active space is defined to include all electron excitations, including core electrons when they are present. Uncontracted (relativistic) aug-cc-pV$n$Z(-DK)basis , basisDK basis sets are used to evaluate (AE) energies, where $n$ is the number of correlating functions present in the basis. Differences between total energies, such as electron affinities and ionisation energies, are provided in the CBS limit using extrapolation Trail_2013_pseudopotentials , Trail_2015_pseudopotentials , Feller_2010 and CPP corrections when required. Extrapolating energies from $n=Q,5$ basis sets and correcting to include CPP energies from $n=T$ basis sets is referred to as ‘$n=(TQ5)$ extrapolation’.Supplemental

II.3 eCEPPs for H, Li–F, Sc–Fe, and Cu

Energy-consistent correlated electron pseudopotentials are constructed for the atoms H, Li–F, Sc–Fe and Cu. This is the set of atoms for which CEPPs are available, with Cu added due to its many properties and applications. There is no core for H, for Li–F a [He] core is represented by the pseudopotentials, and for Sc–Fe and Cu the core is [Ne]. The initial CEPP parameter values are taken from previous publications, except for Cu for which a CEPP is generated using the same method. Core radii for each pseudopotential channel are the same as for the CEPPs.Trail_2013_pseudopotentials , Trail_2015_pseudopotentials The total of core natural orbital occupation numbers arising when generating CEPPs are close to the total number of core electrons, deviating by $7$–$1\times 10^{-3}$ and $6$–$2\times 10^{-3}$ for the first row and $3d$ transition metal atoms, respectively. Across each row this small deviation decreases monotonically with $Z$, and is correlated with the CEPP (and eCEPP) core radii.

For first row atoms $P_{1}$ is defined using the ${}^{2}S$, ${}^{2}P$, ${}^{2}D$ single-valence ground-states, while for the $3d$ transition metal atoms the ${}^{2}F$ state is also included.

A finite number of states are used to define $P_{2}$. Selection rules provide $N+1=$5–8 states for the first row atoms, and $N+1=$8–9 states for the $3d$ transition metal atoms.Supplemental

Our target accuracy for the eCEPPs is better than chemical accuracy, so we take the parameter $\beta$ in Eq. (10) to be 43 meV. The combination of parameters, numerical integration for $P_{1}$, and CCSD(T) data extrapolated to the CBS limit for $P_{2}$ provides the penalty function for a given set of free parameters, $P(\{\lambda^{A}_{l}\},\{\lambda^{a}_{ql}\},r_{0})$.

Optimisation of the penalty function is carried out using the Broyden$-$Fletcher$-$Goldfarb$-$Shanno (BFGS) method,Kelley_1999 with the transformed parameters $(\{\lambda^{A}_{l}\},\{\ln\lambda^{a}_{ql}\},r_{0})$ used to improve stability. Numerical derivatives are defined with a finite difference of $\Delta=10^{-5}$, the smallest value possible given the precision of the CCSD(T) calculations.

Energies and derivatives for the first row atoms are evaluated using the $n=(TQ5)$ extrapolation. However, for the $3d$ transition metal atoms a slightly different approach was required. Convergence with basis set is slow and extrapolation using $n=(TQ5)$ is required to achieve the target accuracy. However, although this choice of basis sets provides CCSD(T) energies converged to the required precision, a significant numerical error in the finite difference prevents the optimisation method from succeeding reliably.

Reliable convergence to the required accuracy is achieved by evaluating the penalty function using $n=(TQ5)$ extrapolation, but evaluating the gradient of the penalty function using energies provided by the $n=Q$ basis only (for both AE and candidate eCEPP total energies). This stabilises and accelerates the optimisation process.

The optimisation was terminated for $P<N/400$. This criterion provides a mean absolute agreement between eCEPP and AE energy differences of 3 meV, with a maximum disagreement of 11 meV occurring for ionisation to the ${}^{6}S$ state of Cr*+*. This optimisation error is small when compared with the estimated CCSD(T) basis set error for all atoms.

Optimisation provides an eCEPP for which the penalty function is small, but non-zero. This corresponds to the charge density of the single valence electron ion relaxing outside of the core region and not exactly reproducing the target MCHF AE density. The degree of relaxation is small and adequately controlled by the penalty function chosen, with the overlap term differing from unity by $<2\times 10^{-4}$ for all atoms and states (the maximum occurs for the ${}^{2}D$ state of Sc). The mean-absolute error for $\phi(r^{c}_{l})$ is $0.6\%$, with all errors less than $3.5\%$.

Pseudopotentials for H are unusual in that the CEPP is equivalent to a norm-conserving HF pseudopotential since neither correlation or exchange occur and no CPP arises. Furthermore, only a few bound atomic states are available and the CEPP accurately reproduces the energies of these states. Consequently, for H the CEPP and eCEPP are identical and equal to a norm-conserving HF pseudopotential.

III Results

III.1 Titanium atom and Titanium Oxide molecule

Two basic requirements of an accurate pseudopotential are to reproduce AE energy differences between atomic states and to reproduce charge densities outside of the core region.

Figure 1 shows the deviation of pseudo-atomic energies from AE energies for several Ti pseudopotentials (including the eCEPP). Energy differences shown are ionisation energies and the electron affinity, evaluated as differences between CCSD(T) energies extrapolated to the CBS limit.

The eCEPP energies agree with AE results to better than both CBS error and chemical accuracy. This is expected in light of the generation procedure. Invoking the eCEPP with no CPP included introduces a small error due to the low polarizability of the [Ne] core.

The figure also provides a useful comparison of the accuracy of the eCEPP with the Trail-Needs Dirac-FockTrail_2005_pseudopotentials , TNDF_website (TNDF) norm-conserving pseudopotential, the Burkatzki-Filippi-DolgBurkatzki_2007 , BFD_website (BFD) energy consistent HF pseudopotential, and the CEPP. For all four pseudopotentials the CBS error consistently increases with ionisation, as expected given that the aug-cc-pV$n$Z(-DK) basis sets are optimised for neutral atoms. Although the CBS error is largely removed by de-contracting basis sets, the width parameters of the Gaussian basis become increasingly sub-optimal as charge densities contract with increasing total charge.

All three non-eCEPPs show an increase in pseudopotential error with increasing charge. This trend is caused by the cancellation of errors in energy differences becoming less perfect with increasing energy difference, and is particularly apparent for the large ionisation energies where electrons are removed from the $4s$ sub-shell, (IP3 and IP4).

Comparing errors arising for all pseudopotential types suggests an interpretation of the errors removed from the eCEPPs. The large errors for the TNDF pseudopotentials suggests that a [Ne] core is necessary for $3d$ transition metal atoms.Pacios_1988 The large errors for the BFD pseudopotential suggest that energy-consistency at the CCSD(T) level of theory is required, and HF is not sufficient. The opposite sign, and similar magnitude, of the errors for the CEPP and BFD pseudopotentials suggests that a combination of CEPP and energy consistency is likely to be successful. The accuracy provided by the eCEPP confirms that it successfully removes these errors.

The eCEPP is not constrained to reproduce the charge density for any system except the single-valence ion, so we compare the AE and pseudo charge densities for a small, neutral molecule. Figure 2 shows such charge densities for the TiO molecule, with an experimental bond length.Harrison_2000 Densities are evaluated using configuration interaction singles and doubles (CISD) calculations, and plotted along the primary symmetry axis. The densities provided by the AE and eCEPP calculations agree well outside of the core region for each atom, demonstrating that good transferability is not limited to atomic ionisation energies.

These results confirm that the eCEPP for titanium accurately reproduces the atomic properties it is designed to recreate, and that the pseudopotentials transfer well to the TiO molecule, accurately reproducing AE charge densities outside of the core regions.

III.2 Small molecules geometries and dissociation energies

To assess the accuracy of the eCEPPs in reproducing molecular dissociation energies and equilibrium geometries we consider a moderately large set of small molecules. A test set is constructed starting with the neutral members of the G1 setG1 which contain only the atoms H and Li–F. We then add LiB, LiC, LiF, H2, BH, Be2, B2, C2, and NO2, to give 38 first row molecules (in 39 states). To this set we add the diatomic molecules composed of Sc–Fe and Cu together with H, C, N, O, and F atoms, selecting those for which the ground state term and dominant configuration is not in doubt, as described by Harrison,Harrison_2000 and for which single-reference CCSD(T) is stable. Titanium is of particular interest,Titanium_tech , Trail_2017 so we also include TiH4, TiO2, and three excited states of TiN. The final test set contains 65 molecules in 69 states.Supplemental

Geometry optimisation is performed for each molecule by minimisation of the AE and pseudopotential CCSD(T) energies with respect to bond lengths and bond–angles that characterise the geometry of each molecule, $x_{i}$. All $x_{i}$ are corrected to include the influences of the CPP, and the basis set error estimated.Supplemental

The differences between molecular geometry spatial parameters arising from the eCEPP and AE Hamiltonians are shown in Fig. 3. The figure also shows the geometry parameters previously reported for the CEPP and TNDF pseudopotentials when available.

All eCEPP geometry parameters fall within chemical accuracy (of 0.01Å) of the relativistic AE values, unlike the results for both the TNDF and CEPP pseudopotentials.

The eCEPP geometry parameters for the mean absolute deviation (MAD) from AE results is ${\sim}1/4\times$ that for either the CEPPs or TNDF pseudopotentials. The overall bias is small for the eCEPPs, with an average deviation of $6\times 10^{-4}$Å.

For all molecules considered the estimated basis set error is negligible, except for H2O2 and H4N2 for which a particularly shallow Born-Oppenheimer surface magnifies this error for certain angular degrees of freedom. There is a weak correlation between errors for the TNDF pseudopotentials and CEPPs, but no correlation remains in the eCEPP errors.

For molecules containing Li, Be, and B, the bond lengths for the eCEPPs are distinctly underestimated, although by less than chemical accuracy. Given the spatially large and highly polarizable [He] cores of these atoms, it seems likely that the dominant source of this error is the incomplete description of core relaxation provided by the CPP.

The negligible error for H2 (of $1(3)\times 10^{-4}$Å) confirms that most of the error for other molecules is due to the difficulty of representing core-valence interaction with a pseudopotential rather than limiting the reproduction of scattering properties to linear order (all three pseudopotentials are equivalent for H).

The improvement in optimised geometries provided by the transition metal eCEPPs is greater than that for the first row atoms. This is particularly apparent for TiO2 for which the notable error in bond-angle present for the other two pseudopotentials is completely removed.

Of all the transition metal molecules, the three containing Cu show the largest underestimation of bond lengths, which is probably due to the eventual emergence of an expected increase in pseudopotential error with the increasing number of valence electrons.

The differences between molecular dissociation energies evaluated using CCSD(T) for the eCEPP and AE Hamiltonians are shown in Fig. 4. The figure also shows the dissociation energies previously reported for the CEPP and TNDF pseudopotentials when available.

All eCEPP dissociation energies fall within chemical accuracy (of 43 meV) of the relativistic AE values, unlike the results for the TNDF and CEPP pseudopotentials. This improvement in accuracy is greater than for optimised geometries, with the MAD for eCEPP ${\sim}1/10\times$ that for TNDF pseudopotentials. For all molecules considered, including H2O2 and H4N2, the estimated basis set error is negligible. Unlike the optimised geometries, errors in dissociation energies for molecules containing Li, Be, or B are not unusually large and the CPP contribution to energies is accurate. The overall bias is small for the eCEPPs when compared with other pseudopotentials, with an average deviation of 2 meV.

There is only a weak correlation between the errors for the TNDF and CEPP potentials, and no correlation is apparent with the eCEPP errors. No consistent trend is apparent for the eCEPPs except for a slow increase in error with the number of valence electrons. As for optimised geometries, the absence of core-valence electron interaction for H2 results in a negligible error (of 0.043(4) meV) for all three pseudopotentials.

Mean, absolute mean, and maximum deviations from AE results for the geometry parameters and dissociation energies evaluated with the three pseudopotentials types are shown in Table 1.

In order to assess the energetics absent from CCSD(T) the available monoxide molecules are considered in more detail. Table 2 show the equilibrium bond-lengths and dissociation energies evaluated with all electrons present, and with core electrons represented by the eCEPPs. We compare this with the experimental data selected by Miliordos et al. for ScO–FeOMiliordos_SctoFe and Harrison for CuO.Harrison_2000 Results are also quoted for TiOPan_2017 and CrOChan_2012 , calculated using a composite method that corrects accurate CCSD(T) energies to include spin-orbit coupling and higher-order correlation.

The disagreement between AE and eCEPP results is consistently negligible compared with the significant disagreement with experimental results, for both geometries and energies.

Bond lengths agree with experiment to within 0.01 Å for all oxides except CrO and MnO, with the largest underestimate of -0.016 Å occurring for CrO. Chan et al. provide a CCSD(T) CrO bond length closer to experiment, but this agreement is due to the small cc-pwCVTZ basis set and disappears with a more complete basis. This suggests that when an error in bond length occurs for these oxides it is primarily due to missing energetics rather than pseudopotential error.

Only the ScO dissociation energy (for AE or eCEPP) agrees with experiment to within 43 meV. The maximum error, of -0.44 eV, occurs for FeO, with CCSD(T) overestimating $D_{e}$ for ScO–VO and underestimating for CrO–CuO. The composite dissociation energies for TiO and CrO are in much better agreement with experiment, with most of the improvement provided by spin-orbit corrections. Higher-order correlation corrections are of secondary importance for these two molecules.

So far we have not considered the accuracy with which the pseudopotentials reproduce molecular potential energy surfaces beyond the minima, that is energy differences considerably smaller than dissociation energies. Hydrogen peroxide is of particular interest due to the relatively large error in equilibrium spatial variables apparent in Fig. 3. This is primarily due to the shallow variation of the total energy with rotation of O-H2 groups around the O-O axis, an internal rotation that plays an important role in conformational analysis.Song_2005

Geometries for H2O2 are optimised at the CCSD(T) level of theory using the uncontracted aug-cc-pVTZ(-DK) basis sets (for AE calculations), with the dihedral angle between O-H2 groups, $\phi$, held constant. Relaxed geometries are evaluated for several values of $\phi$, followed by an extrapolation of the total energy to the CBS limit using $n=(TQ5)$. Energy profiles evaluated with AEs present and using eCEPP pseudopotentials are shown in Fig. 5. Agreement is good, with a maximum pseudopotential and basis set error of 2.5 and 2.7 meV, respectively, and both occurring for $\phi=0$.

Except for TiO and CrO, no attempt has been made to estimate the error due to correlation neglected at the CCSD(T) level of theory. While this error has been quantified as small when compared with basis set errors for first-row molecules,Feller_2006 it is expected to be significant for some molecules containing $3d$ transition metal atoms.Jiang_2012 However, this ‘missing correlation’ is usually found to be between valence electrons in molecules, rather than isolated atoms, suggesting that the accuracy of the eCEPPs will be preserved for calculations including correlation beyond that present in CCSD(T).

III.3 Comparison of errors for a theory hierarchy

An important application of the eCEPPs is expected to be to QMC methods for bulk systems, particularly Diffusion Monte Carlo (DMC). DMC requires a fermionic nodal surface to be supplied, and for bulk systems this is normally approximated as the nodal surface of a determinant constructed from orbitals resulting from a KS-DFT calculation with a plane-wave basis.

In light of this we begin with an AE CCSD(T), then quantify and compare the error introduced by each approximation involved in an eCEPP KS-DFT calculation with a plane-wave basis set. The purpose of this analysis is to compare the relative size of the various components of the error in the final calculation.

If the error due to the eCEPP alone is small when compared with other errors then it may be expected that the contribution of the eCEPP to the error in the approximate nodal surface is also small. While this measure of the accuracy of the nodal surface is not rigorous (an exact exchange-correlation functional would not provide an exact nodal surface, and in principle, an inaccurate pseudopotential or exchange-correlation functional could provide a more accurate nodal surface), it is physically reasonable that less accurate KS-DFT molecular properties will correspond to a less accurate nodal surface.

A separation of errors naturally arises by considering the application of several levels of approximation. The most accurate description considered is AE CCSD(T) with a Gaussian basis. The next level of approximation is to replace atomic cores with the full eCEPP. We then increase the level of approximation progressively by removing the CPP component of the eCEPP, followed by replacing CCSD(T) with KS-DFT, and finally by replacing the full eCEPP with a Kleinman-BylanderKleinman_1982 representation (of the eCEPP) together with a plane-wave basis. Both KS-DFT calculations are performed using the PBE exchange-correlation functional.Perdew_1996

The KB pseudopotential representation of the eCEPP is given by

[TABLE]

with $\delta V_{l}=V^{pp}_{l}-V^{pp}_{l_{max}}$. The $\phi_{li}$ functions used for each projector are supplied for each channel of each pseudopotential, and the matrices for each channel are given by $B_{l,ij}=\left<\phi_{li}|\delta V_{l}|\phi_{lj}\right>^{-1}$. This separable form is less accurate than the semi-local spherical harmonic projector representation used to define the eCEPPs. However, the KB representation provides a considerable reduction in computational cost, and the additional error introduced is often small.

Three projection orbitals are used for first row atoms, and six projection orbitals for the $3d$ transition metal atoms. For all atoms the local channel is taken to be that with the largest $l$ value, $l_{max}$. In each case orbitals are taken from the neutral atom, and supplemented by orbitals from the lowest-energy excited states that provide projectors for all pseudopotential channels.Supplemental For example, for Ti the ground state provides projectors from the $3s$, $4s$, $3p$ and $3d$ orbitals, to which we add $4p$ and $4f$ projectors from the [Ar]$3d^{2}4s4p$ and [Ar]$3d^{2}4f$ states, respectively.

Note that projectors for $l=l_{max}$ are evaluated but unused. These are provided to allow flexibility in the choice of the local channel of the KB representation, should it be required. All projectors are generated using the PBE functional and the atomic KS-DFT code contained in the Quantum Espresso package.Giannozzi_espresso_2009

Errors are assessed for the subset of diatomic molecules present in the full test set, excluding the excited TiN states. Bond lengths are relaxed and dissociation energies evaluated for each of these molecules, and at each level of theory.

All electron and eCEPP results with the CPP are identical to those included in Sec. III.2. Dissociation energies with the CPP removed are evaluated with the same data but excluding CPP corrections.

Results for KS-DFT without the KB representation are generated using the uncontracted aug-pp-pV$5$Z Gaussian basis sets. Extrapolation is not required as exponential convergence results from correlation being contained in the exchange-correlation functional, rather than the wave function itself. Estimated basis set error is $<$8 meV, from differences between dissociation energies evaluated using the $n=(Q5)$ basis sets.

For KS-DFT with a plane-wave basis set the isolated molecule is modelled by a single molecule in a $(15$Å$)^{3}$ periodic unit cell, and using the CASTEPcastep plane-wave pseudopotential code (with a small modification to allow for more than one projector per pseudopotential channel). Calculations for all atoms and molecules are performed using the plane-wave basis set defined by an energy cutoff of 150 Ha, resulting in a basis set error of $<$9 meV. The maximum error occurs for CuF due to the number of valence electrons present and the depth of the Cu and F potentials. The same accuracy could be achieved with considerably smaller basis sets for most of the molecules in our test set.

The total error arising from all approximations taken together is

[TABLE]

Components of the total error due to pseudopotentials, the absent CPP, density functional theory (with PBE), and the KB representation are given by

[TABLE]

respectively. The superscripts for each dissociation energy in Eq. (III.3) describe the level of approximation. Superscripts ‘AE’, ‘pp’, and ‘pp-cpp’ denote calculations including all electrons, describing core electrons using an eCEPP, and describing core electrons using an eCEPP with core polarisation removed. Superscripts ‘CCSD(T)’, ‘DFT’, and ‘DFT+KB’ denote the use of CCSD(T), KS-DFT, and KS-DFT with a KB representation of the eCEPPs.

The error breakdown is shown in Fig. 6. The overall trend is for KS-DFT errors to be dominant when compared with the three errors associated with the pseudopotential. The negligible error that arises for the CrF, CH, and OH molecules is unlikely to be more than fortuitous. The largest KS-DFT errors arise for molecules containing transition metal atoms, which is unsurprising given the strong correlation and highly inhomogeneous charge densities for these molecules.

Errors due to the absence of a CPP are almost negligible compared to the KS-DFT error. Molecules containing Li, with the most polarizable core, show the largest contribution from the CPP, although the net polarizability stands out as particularly large for CN, N2, CO and C2.

The KB representation results in a relatively small error. No overall pattern is apparent except that the KB representation is significantly less accurate for molecules containing transition metal atoms, and the error is greatest in TiC. This is as expected, since the non-local and local channels differ the most for these atoms, with the eCEPPs composed a deep $d$ channels and relatively shallow $s$, $p$, and $f$ channels.

Physically unrealistic ‘ghost’ eigenstates can arise when the KB representation is used.Gonze_1991 Agreement between KS-DFT energy differences evaluated with and without the KB representation (see Fig. 6) provide a useful check that such ghost states do not occur for eCEPPs with the local channel chosen as $l=l_{max}$, and this is further confirmed by good agreement between total energies (not shown). This absence of ghost states for the $3d$ transition metals is very unlike the TNDF pseudopotentials, Drummond_2016 and can be ascribed to the presence of a physically realistic eCEPP channel for $l\geq 3$.

The breakdown of errors in bond lengths (not shown) shows similar behaviour, except that the error due to the absence of a CPP is comparable to the KS-DFT error for molecules containing Li, Be, and B, and KB error is dominant for ScH and ScF. Unlike the dissociation energies, bond lengths for CrF, CH, and OH are not notably accurate.

III.4 Optimised basis sets

Contracted correlation consistent Gaussian basis sets are well established to be computationally efficient for correlated electron calculations, providing controlled convergence behaviour that allows accurate extrapolation to the CBS limit.

The eCEPP pseudo atom wave-functions are very different from AE wave-functions in the core region, hence such a basis set must be reconstructed for each eCEPP. We name basis sets optimised for the eCEPPs ‘aug-cc-pV$n$Z-eCEPP’. As for the AE case, these are constructed by optimising a Gaussian basis within HF, systematically introducing correlation consistent components of the basis at higher angular momentum quantum numbers, and augmenting the basis with further Gaussians optimised to accurately represent the more diffuse anionic state.

The basis sets are optimised using the procedure developed for first row atoms by Xu et al.,Xu_2013 and the $3d$ transition metal approach of Balabanov et al.,basisDK adapted to cores described by pseudopotentials.Supplemental

In the presence of pseudopotentials, valence states become smoother and core states are absent, hence fewer Gaussians and contractions are required than for the AE case. This results in a significant improvement in computational efficiency, and the average CCSD(T) run-time is $\sim 1/5\times$ that for uncontracted standard basis sets. This speedup should increase with molecular size.

Figure 7 shows the errors in eCEPP CCSD(T) dissociation energies evaluated using both the uncontracted standard aug-cc-pV$n$Z basis set, and the new contracted aug-cc-pV$n$Z-eCEPP basis set. The error is quantified as the deviation from the baseline of AE CCSD(T) energies evaluated using the uncontracted standard aug-cc-pV$n$Z-DK basis set.

For molecules containing first row atoms the new contracted basis sets reproduce the uncontracted aug-cc-pV$n$Z basis set errors accurately. Molecules containing $3d$ transition metal atoms show a less consistent agreement. However, the additional error is small and marginally greater than chemical accuracy only for the three molecules containing copper. For most cases the error due to the contracted aug-cc-pV$n$Z-eCEPP basis is similar or less than that due to the pseudopotential itself.

In Fig. 7 we also show the error introduced to AE CCSD(T) by using contracted standard aug-cc-pV$n$Z-DK basis sets. This is quantified as the difference between dissociation energies evaluated using the contracted and uncontracted aug-cc-pV$n$Z-DK basis sets, with the latter being the more complete basis, which provides a useful scale for assessing the accuracy of the aug-cc-pV$n$Z-eCEPP basis sets.

The basis set error introduced by contraction is significantly greater for the AE calculations than for the eCEPP calculations (except for H2 for which both are negligible), especially for the molecules containing $3d$ transition metal atoms. We conclude that the construction of aug-cc-pV$n$Z-eCEPP basis sets and application with the eCEPPs is successful, with an impressively low error for the test set considered.

Finally, we note that the relatively large error in the relativistic AE dissociation energy calculated with the contracted aug-cc-pV$n$Z-DK basis sets is not unexpected, and is not a deficiency of these basis sets per se. These standard contractions are constructed to represent correlations between electrons outside of a [He] core for Li–F, and outside of an [Ar] core for Sc–Cu, so do not describe correlation between these cores and valence electrons well, despite the CCSD(T) active space including excitation of all electrons.

This explanation is consistent with the very small error achieved for H2, with the larger errors for the $3d$ transition metal molecules, and the accuracy of the aug-cc-pwCV$n$Z basis setsbasisDK for which contractions include core-valence correlation effects.

IV Conclusions

The paradigms of shape and energy consistency in the definition of electronic pseudopotentials are successfully and consistently combined to generate accurate effective potentials. The resulting eCEPPs reproduce core electron scattering properties to first order, reproduce atomic excitation and ionisation energies, and partially reproduce core relaxation effects via a dipole polarisation potential. This approach is an extension of the CEPPs (that reproduce only core electron scattering and core polarisation). Electron correlation is included throughout, no recourse is made to independent electron orbitals, and the generation process is ab initio.

Comparing AE and eCEPPs results for relaxed geometries and dissociation energies, and for a test set of 65 molecules (in 69 states) at the CCSD(T) level of theory, demonstrates a reduction in error by ${\sim}1/4\times$ for geometry parameters, and by ${\sim}1/10\times$ in dissociation energies when compared with Hartree-Fock pseudopotentials.

When compared with the previously published CEPPs, the eCEPPs introduce an interpolation of the effect of the core on several atomic states in addition to the extrapolation of core scattering. This extension is significant, reducing errors in geometry parameters and dissociation energies by a factor of ${\sim}1/4\times$ and ${\sim}2/5\times$.

The eCEPPs will primarily be useful for correlated wave-function methods, particularly for a CCSD(T) description of molecular systems, and for the application of QMC methods to extended systems. For both cases pseudopotentials are often unavoidable and the eCEPPs provided here are validated by small errors for a reasonably large test set.

The improved accuracy is particularly important when $3d$ transition metal atoms are present as the eCEPP error is sufficiently small to resolve the effect of spin-orbit coupling and aspects of electron-electron correlation often absent in ab initio many-body methods.

Further work naturally suggests itself. It would be useful and desirable to extend the eCEPP generation process to more atomic types. Such an extension is however limited by the size and stability of MCHF and CCSD(T) calculations, and the absence of correlation consistent basis sets for heavier atoms. However, completing the H–Kr set of atoms can be expected to be straightforward.

It is perhaps more interesting to consider the accuracy of the eCEPPs, and how much further improvement is possible. In principle, the accuracy could be improved by extending the reproduction of core scattering properties to higher order and many valence electrons. As discussed in the introduction, implementing this is non-trivial as it will result in a many-body effective potential.

Another possible strategy would be to extend the reproduction of atomic ionisation and excitation energies to valence charge densities outside of the core region. Given the success of the eCEPP in reproducing these densities for TiO, such an extension appears to be unnecessary.

The most direct approach to improving the eCEPP accuracy lies in the choice of excited/ionic states for which energy consistency is enforced. Here we have considered the atomic states with the lowest energy. A more detailed consideration of the molecular bonds each atom is prone to form should, in principle, allow a selection of atomic states which are more appropriate for each atom. This could also be extended to molecular geometries and dissociation energies. The primary drawback of this approach (and why it was not used here) is that the pseudopotential generation process would become more biased.

A simple, and possibly effective, approach to reducing pseudopotential error would be a more careful choice of core radii. These radii parameterise a balance between errors due to defining the eCEPP with a finite core charge density outside of the core region, and errors due to replacing the valence charge density with a pseudo density within the core region. Optimised core radii, different from those used in this paper, may well decrease errors significantly. This is likely to be most influential for the lightest atoms with the most diffuse cores.

The most important issue that has not been addressed is the accuracy of the eCEPPs when used in QMC calculations. CCSD(T) results are accurate, and DMC is known to provide a similar or better accuracy than CCSD(T). However, DMC involves a further ‘pseudopotential localisation’ approximation that expresses the semi-local one-body pseudopotential(s) as a local many-body potential via a trial wave-function. The resulting error is second order in the accuracy of the trial wave-function.Foulkes_QMC_review , Casino_reference A measure of this error would be desirable, and the molecules considered here are useful candidates for such tests.

Supplementary Material

See supplementary material for details of CCSD(T) calculations, atomic and molecular states, KB projectors, eCEPP basis set optimisation, and the eCEPPs themselves.

Acknowledgements.

R.J.N. and J.R.T. acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the U.K. [EP/J017639/1].