Effective Core Potentials¶
Overview¶
Effective core potentials (ECPs) are a useful means of replacing the core electrons in a calculation with an effective potential, thereby eliminating the need for the core basis functions, which usually require a large set of Gaussians to describe them. In addition to replacing the core, they may be used to represent relativistic effects, which are largely confined to the core. In this context, both the scalar (spin-free) relativistic effects and spin-orbit (spin-dependent) relativistic effects may be included in effective potentials. NWChem has the facility to use both, and these are described in the next two sections.
A brief recapitulation of the development of RECPs is given here, following L.F. Pacios and P.A. Christiansen, J. Chem. Phys. 82, 2664 (1985). The process can be viewed as starting from an atomic Dirac-Hartree-Fock calculation, done in jj coupling, and producing relativistic effective potentials (REPs) for each \(l\) and \(j\) value, \(U^{\rm REP}_{lj}\) which for example contains the Coulomb potential of the core electrons balanced by the part of the nuclear attraction which cancels the core electron charge. The residue is expressed in a semi-local form,
where \(L\) is one larger than the maximum angular momentum in the atom. The scalar potential is obtained by averaging the REPs for each \(j\) for a given \(l\) to give an averaged relativistic effective potential, or AREP
These are summed into the full potential
The spin-orbit potential is obtained from the difference between the REPs for the two \(j\) values for a given \(l\), and may be represented in terms of an effective spin-orbit operator,
where
The relavistic potential \(U^{REP}\) is the sum of \(H^{SO}\) and \(U^{AREP}\).
The spin-orbit integrals generated by NWChem are the integrals over the sum, including the factor of \(2/(2l+1)\) as an effective spin-orbit operator without further factors introduced.
The effective potentials, both scalar and spin-orbit, are fitted to Gaussians with the form
where \(A_{lk}\) is the contraction coefficient, \(n_{lk}\) is the exponent of the \(r\) term (r-exponent), and \(B_{lk}\) is the Gaussian exponent. The \(n_{lk}\) exponent is shifted by 2, in accordance with most of the ECP literature and implementations, i.e., \(n_{lk} = 0\) implies \(r^{-2}\). The current implementation allows \(n_{lk}\) values of only 0, 1, or 2.
Scalar ECPs¶
The optional directive ECP allows the user to describe an effective core
potential (ECP) in terms of contracted Gaussian functions as given
above. Potentials using these functions must be specified explicitly by
user input in the ECP directive. This directive has essentially the same
form and properties as the standard BASIS directive, except for
essential differences required for ECPs. Because of this, the ECP is
treated internally as a basis set. The form of the input for the ECP
directive is as
follows:
ECP [<string name default "ecp basis">] \
[print || noprint default print]
<string tag> library [<string tag_in_lib>] \
<string standard_set> [file <filename>] \
[except<string tag list>]
<string tag> [nelec] <integer number_of_electrons_replaced>
...
<string tag> <string shell_type>
<real r-exponent> <real Gaussian-exponent> <real list_of_coefficients>
...
END
ECPs are automatically segmented, even if general contractions are
input. The projection operators defined in an ECP are spherical by
default, so there is no need to include the CARTESIAN or SPHERICAL
keyword as there is for a standard basis set. ECPs are associated with
centers in geometries through tags or names of centers. These tags must
match in the same manner as for basis sets the tags in a GEOMETRY and
ECP directives, and are limited to sixteen (16) characters. Each center
with the same tag will have the same ECP. By default, the input module
prints each ECP that it encounters. The NOPRINT option can be used to
disable printing. There can be only one active ECP, even though several
may exist in the input deck. The ECP modules load ecp basis inputs
along with any ao basis inputs present. ECPs may be used in both
energy and gradient calculations.
ECPs are named in the same fashion as geometries or regular basis sets, with the default name being “ecp basis”. It should be clear from the above discussion on geometries and database entries how indirection is supported. All directives that are in common with the standard Gaussian basis set input have the same function and syntax.
As for regular basis sets, ECPs may be obtained from the standard library. For a complete list of basis sets and associated ECPs in the NWChem library see the available basis sets or the Basis Set Exchange for naming conventions and their specifications.
The keyword nelec allows the user to specify the number of core
electrons replaced by the ECP. Additional input lines define the
specific coefficients and exponents. The variable <shell_type> is used
to specify the components of the ECP. The keyword ul entered for
<shell_type> denotes the local part of the ECP. This is equivalent to
the highest angular momentum functions specified in the literature for
most ECPs. The standard entries (s, p, d, etc.) for shell_type specify
the angular momentum projector onto the local function. The shell type
label of s indicates the ul-s projector input, p indicates the ul-p,
etc.
For example, the Christiansen, Ross and Ermler ARECPs are available in
the standard basis set library named crenbl_ecp. To perform a
calculation on uranyl UO22+ with all-electron oxygen
(aug-cc-pvdz basis), and uranium with an ARECP and using the
corresponding basis the following input can be used
geometry
U 0 0 0
O 0 0 1.65
O 0 0 -1.65
end
basis
U library crenbl_ecp
O library aug-cc-pvdz
end
ecp
U library crenbl_ecp
end
The following is an example of explicit input of an ECP for H2CO.
It defines an ECP for the carbon and oxygen atoms in the molecule.
ecp
C nelec 2 # ecp replaces 2 electrons on C
C ul # d
1 80.0000000 -1.60000000
1 30.0000000 -0.40000000
2 0.5498205 -0.03990210
C s # s - d
0 0.7374760 0.63810832
0 135.2354832 11.00916230
2 8.5605569 20.13797020
C p # p - d
2 10.6863587 -3.24684280
2 23.4979897 0.78505765
O nelec 2 # ecp replaces 2 electrons on O
O ul # d
1 80.0000000 -1.60000000
1 30.0000000 -0.40000000
2 1.0953760 -0.06623814
O s # s - d
0 0.9212952 0.39552179
0 28.6481971 2.51654843
2 9.3033500 17.04478500
O p # p - s
2 52.3427019 27.97790770
2 30.7220233 -16.49630500
end
Various ECPs without a local function are available, including those of
the Stuttgart group. For those, no ul part needs to be defined. To
define the absence of the local potential, simply specify one
contraction with a zero coefficient:
<string tag> ul
2 1.00000 0.00000
Spin-orbit ECPs¶
The Spin-orbit ECPs can be used with the Density Functional Approach, but one has to run the calculations without symmetry. Note: when a Hartree-Fock method is specified the spin-orbit input will be ignored.
Spin-orbit ECPs are fitted in precisely the same functional form as the scalar RECPs and have the same properties, with the exception that there is no local potential ul, no s potential and no effective charge has to be defined. Spin-orbit potentials are specified in the same way as ECPs except that the directive SO is used instead of ECP. Note that there currently are no spin-orbit ECPs defined in the standard NWChem library. The SO directive is as follows:
SO [<string name default "so basis">] \
[print || noprint default print]
<string tag> library [<string tag_in_lib>] \
<string standard_set> [file <filename>]
[except `<string tag list>]
...
<string tag> <string shell_type>
<real r-exponent> <real Gaussian-exponent> <real list_of_coefficients>
...
END
Note: in the literature the coefficients of the spin-orbit potentials are NOT always defined in the same manner. The NWChem code assumes that the spin-orbit potential defined in the input is of the form: \(\(\Delta U^{NWChem}_{l} = \,\! \frac{2}{2l+1} \Delta U_{l}\)\)
For example, in the literature (most of) the Stuttgart potentials are
defined as \(\Delta U_{l}\) and, hence, have to be multiplied by \(2/(2l+1)\)
(Note: On the Stuttgart/Köln web pages
https://www.tc.uni-koeln.de/PP/clickpse.en.html,
spin-orbit potentials have already been corrected by the appropriate scaling factor and can be used as is). On the other hand, the CRENBL
potentials in the published papers are defined as
\(\,\!\frac{l}{2l+1} \Delta U_{l}\)
have been corrected with the \(2/l\) factor, so make sure the
appropriate scaling is applied).
For example, to use the Stuttgart/Köln ECP and SO-ECP for Hg (ECP60MDF) in NWChem.
The following URL will display bot the the ECP and SO parts.
http://www.tc.uni-koeln.de/cgi-bin/pp.pl?language=en,format=molpro,element=Hg,job=getecp,ecp=ECP60MDF
The highlighted section (last four lines) below is the SO part.
The un-highlighted part (first five lines) is the ECP.
! Q=20., MEFIT, MCDHF+Breit, Ref 37.
ECP,Hg,60,5,4;
1; 2,1.000000,0.000000;
2; 2,12.413071,275.774797; 2,6.897913,49.267898;
4; 2,11.310320,80.506984; 2,10.210773,161.034824; 2,5.939804,9.083416; 2,5.019755,18.367773;
4; 2,8.407895,51.137256; 2,8.214086,76.707459; 2,4.012612,6.561821; 2,3.795398,9.818070;
2; 2,3.273106,9.429001; 2,3.208321,12.494856;
2; 2,4.485296,-6.338414; 2,4.513200,-8.099863;
4; 2,11.310320,-161.013967;2,10.210773,161.034824;2,5.939804,-18.166832;2,5.019755,18.367773;
4; 2,8.407895,-51.137256; 2,8.214086,51.138306; 2,4.012612,-6.561821; 2,3.795398,6.545380;
2; 2,3.273106,-6.286001; 2,3.208321,6.247428;
2; 2,4.485296,3.169207; 2,4.513200,-3.239945;
! References:
! [37] D. Figgen, G. Rauhut, M. Dolg, H. Stoll, Chem. Phys. 311, 227 (2005).
The corresponding NWChem input is
ecp
Hg nelec 60
Hg ul
2 1.0000000 0.0000000
Hg S
2 12.4130710 275.7747970
2 6.8979130 49.2678980
Hg P
2 11.3103200 80.5069840
2 10.2107730 161.0348240
2 5.9398040 9.0834160
2 5.0197550 18.3677730
Hg D
2 8.4078950 51.1372560
2 8.2140860 76.7074590
2 4.0126120 6.5618210
2 3.7953980 9.8180700
Hg F
2 3.2731060 9.4290010
2 3.2083210 12.4948560
Hg G
2 4.4852960 -6.3384140
2 4.5132000 -8.0998630
end
so
Hg P
2 11.310320 161.013967
2 10.210773 161.034824
2 5.939804 -18.166832
2 5.019755 18.367773
Hg D
2 8.407895 -51.137256
2 8.214086 51.138306
2 4.012612 -6.561821
2 3.795398 6.545380
Hg F
2 3.273106 -6.286001
2 3.208321 6.247428
Hg G
2 4.485296 3.169207
2 4.513200 -3.239945
end