DFT: LCAO¶
Introduction¶
QuantumATK can model the electronic properties of closed and open quantum systems within the framework of density functional theory (DFT) using numerical LCAO basis sets (linear combination of atomic orbitals). For periodic systems, QuantumATK can also use the planewave basis set for DFT calculations, as discussed in DFT: Plane Wave.
The key parameter in the selfconsistent calculation of the Kohn–Sham equations is the density matrix, which defines the electron density. For open systems, the density matrix is calculated using nonequilibrium Green’s functions (NEGFs), see NEGF: Device Calculators, while for closed or periodic systems it is calculated by diagonalization of the Kohn–Sham Hamiltonian, as described in the present chapter. The electron density then sets up an effective potential, which is given by the Hartree, exchangecorrelation, and external potentials. Knowing the effective potential allows us to obtain the Kohn–Sham Hamiltonian.
The next section describes the mathematical formalism behind the DFTLCAO model.
Background information¶
The LCAOCalculator provides a description of electronic structure using DFT and normconserving pseudopotentials. The method is based on an expansion of the singleparticle wave functions in a basis of numerical atomic orbitals with compact support. In this section, we describe the mathematical formalism of the DFTLCAO calculator, which closely follows the work of Soler et al. [SAG+02].
Kohn–Sham Hamiltonian¶
Within density functional theory, the manybody electronic structure of the system is described in terms of the oneelectron Kohn–Sham Hamiltonian:
In this equation, the first term is the kinetic energy of the electron, while the second term (the effective potential) is the potential energy of the electron moving in the mean field created by the other electrons as well as in any external potential field, e.g. the electrostatic potential of ions or any other external field. The electrons are described in terms of the total electron density, \(n=n(\mathbf{r})\).
The electron density is discussed in detail in the section Electron density, and the effective potential in the section Effective potential.
Solving the Kohn–Sham equations by means of a basis set expansion¶
We calculate the oneelectron eigenfunctions of the Kohn–Sham Hamiltonian, \(\psi_{\alpha}\), by solving the oneelectron Schrödinger equation,
This differential equation is called the Kohn–Sham equation within DFT. To solve it, we expand the eigenfunctions \(\psi_{\alpha}({\mathbf{r}})\) in a set of basis functions, \(\phi_i\):
This allows us to represent the differential equation as a matrix equation for determining the expansion coefficients, \(c_{\alpha i}\):
where the Hamiltonian matrix, \(H_{ij} = \langle \phi_i  \hat{H}_{\mathrm{1el}}  \phi_j \rangle\), and overlap matrix \(S_{ij} = \langle \phi_i  \phi_j \rangle\) are given by the multiple integrals with respect to the electron coordinates.
Electron density¶
The electron density of the manyelectron system is given by the occupied eigenstates of the Kohn–Sham Hamiltonian:
where \(f_\alpha\) is the occupation of the level denoted by \(\alpha\). For finite temperature calculation the occupations are determined by the FermiDirac distribution \(f_\alpha = \frac{1}{1 + e^{(\epsilon_\alpha  \epsilon_F)/kT}}\) but other smooth distributions may be introduced in order to help speed convergence (see Occupation Methods).
The electron density can then be expressed in terms of the density matrix:
where the density matrix is given by the basis set expansion coefficients \(c_{\alpha i}\):
Note
For open systems (DeviceConfiguration), the density matrix is calculated using nonequilibrium Green’s functions, see NEGF: Device Calculators.
Electron difference density¶
It is often convenient to compare the electron density of the manybody system to a superposition of individual atombased electron densities, \(n^{\mathrm{atom}}(\mathbf{r}\mathbf{R}_\mu)\), where \(\mathbf{R}_\mu\) is the position of atom \(\mu\) in the manybody system:
\(\Delta n\) is called the electron difference density, and it is calculated using the ElectronDifferenceDensity analysis object.
Effective potential¶
The effective potential, \(V^\mathrm{eff}[n]\), has three contributions:
The first two terms are due to electron–electron interactions, which depend on the electron density. The first term, \(V^{H}[n]\), is the Hartree potential due to the meanfield electrostatic interaction between the electrons, while the second term, \(V^\mathrm{xc}[n]\), is the exchangecorrelation potential, which arises from the quantum mechanical nature of the electrons.
The potential \(V^\mathrm{ext}\) represents any other electrostatic fields in the system. It can be separated into two contributions; the electrostatic potential of ions (given by normconserving pseudopotentials) and external electrostatic fields (given by one or more external sources).
External potential¶
The external potential is given by the pseudopotentials and an external electrostatic field. Such a field may arise from the inclusion of metallic gates:
The pseudopotential has two contributions:
where the first and second terms correspond to the local and nonlocal part of the pseudopotential, respectively.
The term \(V^\mathrm{gate}\) is the electrostatic potential due to external
gates, calculated with zero electron density. This term will be returned by the
Analysis object ExternalPotential
and is calculated using the
MultigridSolver
, DirectSolver
and
ParallelConjugateGradientSolver
Poisson solvers. The Analysis object
ElectrostaticDifferencePotential
also includes this term, among others.
Hartree potential¶
Using the electron density, we can calculate the classical electrostatic potential, the socalled Hartree potential. The actual calculation of the Hartree potential is described in detail in the section The Hartree Potential.
Exchangecorrelation potential¶
In the DFT method, the quantum mechanical part of the electron–electron interaction is approximated by the exchangecorrelation term, and a large number of different approximate exchangecorrelation density functionals exist. ATK supports many of these, see the section Exchangecorrelation energy.
The exchangecorrelation potential is defined as the functional derivative of the exchangecorrelation energy with respect to the electron density, which corresponds to a meanfield quantum mechanical interaction potential between the electrons:
Total energy and forces¶
The DFT total energy of a manyelectron system is a functional of the electron density, \(n\):
where \(T[n]\) is the kinetic energy of a noninteracting electron gas with density \(n\), \(E^\mathrm{xc}[n]\) the exchangecorrelation energy, \(E^{H}[n]\) the Hartree potential energy, and \(E^\mathrm{ext}[n]\) the interaction energy of the electrons in the electrostatic field created by ions and other external sources.
The electron kinetic energy may be defined as
The total energy is calculated using TotalEnergy.
Firstprinciples forces are calculated by differentiating the total energy with respect to the ionic coordinates of atom \(i\) at position \(\mathbf{R}_i\):
Pseudopotentials¶
QuantumATK uses normconserving pseudopotentials and PAW potentials, and is shipped with a database for the entire periodic table. See the full list here: Pseudopotentials. This database is reviewed and updated on a regular basis to provide increasingly accurate and generalpurpose pseudopotentials and PAW potentials.
QuantumATK uses the unified pseudopotential format (upf
) defined by the
Quantum ESPRESSO consortium. From their website,
one may download tools to convert several different pseudopotential formats into a
upf
file.
LCAO basis set¶
The eigenfunctions of the Kohn–Sham Hamiltonian can be expanded in a Linear Combination of Atomic Orbitals (LCAO’s):
where \(Y_{lm}\) are spherical harmonics, and \(R_{nl}\) are radial functions with compact support, being exactly zero outside a confinement radius.
Note
The basis set functions have a finite range, but the interaction range of the Hamiltonian is larger than that of the basis set. Theoretically, the Hamiltonian interaction range may exceed twice the basis set range, but is usually smaller in practice.
The basis orbitals have a number of parameters that determine the shape of the orbitals. It is possible to assemble the basis orbitals into your own basis set through the use of the BasisSet keyword.
QuantumATK comes with a number of prebuilt basis sets for each chemical element and for each type of pseudopotential. See the full list here: LCAO basis sets
Exchangecorrelation energy¶
The QuantumATK package use the libxc library for providing a large suite of exchangecorrelation functionals. For each functional it is possible to add a Hubbard correction, as described in the section XC+U meanfield Hubbard term.
The main families of exchangecorrelation functionals supported in QuantumATK are the Localdensity approximation (LDA), the Generalizedgradient approximation (GGA), and the MetaGGA. Each functional comes in four spin variants:
 nonpolarized,
 spinpolarized (collinear),
 noncollinear,
 noncollinear with spinorbit coupling.
Important
The choice of spin variant for the exchangecorrelation functional determines the spin type used in all of the DFTLCAO calculation.
Furthermore, a Hubbard correction can be added to all exchangecorrelation variants, regardless of spin type, see XC+U meanfield Hubbard term. Several standard GGA functionals can also be extended with the DFTD2 and DFTD3 dispersion corrections by Grimme and coworkers.
Initial spin¶
Both collinear and noncollinear calculations may have local minima of the total energy that correspond to different spin configurations. It is therefore important to prepare the system in the correct initial spin configuration, such that the DFT selfconsistent calculation will end up in the lowest energy spin configuration. This is done through the use of the methods InitialSpin and RandomSpin. The initial spin direction on each atom is specified in physical spherical coordinates \((r,\theta,\phi)\), where \(\theta\) is the angle with the zaxis, and \(\phi\) the polar angle in the xy plane relative to the xaxis. The collinear case is \(\theta=0\) Radians or \(\theta=\pi\) Radians.
Localdensity approximation (LDA)¶
In the LDA approximation, the exchangecorrelation energy is taken as a functional of the local electron density,
where \(\varepsilon^\mathrm{LDA}(n(\mathbf{r}))\) is the exchangecorrelation energy density of a homogeneous electron gas with density \(n(\mathbf{r})\).
It is possible to derive an exact, analytical expression for the exchange energy of the homogeneous electron gas, the socalled Dirac–Bloch exchange energy, which is used for all the LDA functionals. The correlation energy of the electron gas cannot be calculated exactly, so a number of different approximations have been proposed over the years. Most of them give almost identical results; the most commonly used variant is the parametrization by Perdew and Zunger, which is also the default LDA functional in DFT: LCAO.
For a list of the parametrizations for the LDA correlation functional implemented in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchangecorrelation functional. See also Table 21, Table 24, and Table 25.
Generalizedgradient approximation (GGA)¶
The GGA functionals is a large family of semilocal approximations for the exchangecorrelation energy, where the functional depends on both the local value and the local gradient of the electron density,
For a list of the GGA exchangecorrelation functionals available in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchangecorrelation functional. See also Table 22, Table 24, and Table 25.
MetaGGA¶
The metaGGA (MGGA) family of functionals extend the GGA approximation by additionally depending on one or both of the following quantities: the Laplacian of the density
and the socalled kinetic energy density
For a list of the MGGA exchangecorrelation functionals available in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchangecorrelation functional. See also Table 23, Table 24, and Table 25.
Different formulations of the MGGA can produce functionals with different uses. The SCAN functional has been found to give improved energetics over the LDA and GGA in several cases; however, band gaps are similar.
Conversely, the TB09 functional can provide a much more accurate description of band gaps than ordinary LDA and GGA. In fact, the accuracy of semiconductor band gaps obtained with TB09 are often comparable to those from calculations using GW or hybrid functionals, which have a significantly higher computational cost. It is even possible to tune the TB09 method such that it works optimally on both sides of an interface (see ExchangeCorrelation). Alternatively, one can use the Hubbard U model within the LDA or GGA with the U parameters for the semiconductor elements fitted against TB09 results.
However, TB09 is not meant for energetics, and, hence, geometry optimizations are disabled with this functional. For this, the GGA or other MGGA functionals should be used.
XC+U meanfield Hubbard term¶
The local approximations to the exchangecorrelation energy have a number of shortcomings. Two of these are of particular interest:
 Selfinteraction: The electron is formally allowed to interact with itself. This can prevent electrons from localizing properly.
 Excited states: The LDA and GGA description of conductionband energy levels is often poor, so band gaps are often too low.
The mean field Hubbard correction by Dudarev et al. [DBS+98] and Cococcioni et al. [CdG05], often denoted XC+U, DFT+U, LDA+U, or GGA+U, is a semiempirical correction which attempts to improve on these deficiencies of the local exchangecorrelation functionals by adding an extra term to the exchangecorrelation functional:
In this equation, \(n_\mu\) is the projection onto an atomic shell and \(U_\mu\) is the Hubbard U for that shell. The \(E_{U}\) energy term is zero for a fully occupied or unoccupied shell, while positive for a fractionally occupied shell.
The energy is thereby lowered if states become fully occupied or empty. This may happen if the energy levels move away from the Fermi Level, increasing the band gap, or if the broadening of the states is decreased, the electrons are then more localized. Thus, the Hubbard U method improves on the deficiencies of the local exchangecorrelation functionals listed above.
The value of \(U\) is often used as an empirical parameter, which is varied in order to improve the comparison between DFT and experimental data. However, it has also been suggested that the value of U may be obtained by minimizing the total energy of the system [CdG05].
Some caution must be taken when using the XC+U correction, since for large U values the electron density may have several local minima where some of them are unphysical, and it is often necessary to select an anisotropic initial electron state to reach the correct local minimum.
XC+U implementations in QuantumATK¶
The XC+U implementation in QuantumATK comes in four main variants; a standard localorbital formulation, socalled Onsite representation, and the Dual representation introduced by Han et al. [HOY06]. The difference between the two implementations is in the definition of the local occupation matrix. The occupations may be summed over a single orbital or over an entire angular momentum shell, the latter corresponds to the OnsiteShell and DualShell representations.
Onsite representation¶
In the onsite representation, the local occupation matrix is given by
where \(D\) is the density matrix, \(\mu\) a basis orbital index, and \(\sigma\) a spin index.
OnsiteShell representation¶
In the onsite shell representation, the occupation is obtained by summing together all the basis functions in each angular momentum shell, \(l\):
Dual representation¶
In the dual representation, the local occupation matrix is given by
where \(S\) is the overlap matrix.
DualShell representation¶
Similar to the onsite shell representation, the dual shell occupation is obtained by summing up all dual occupations in each angular momentum shell, \(l\):
Which model to use?¶
The onsite representation corresponds to projecting into a local orbital, and since the local orbitals are nonorthogonal, there are no sum rules for the occupation matrix and it is not guaranteed that a fully occupied shell has occupation 1.
In the dual representation, the occupation matrix has the form of a Mulliken population, and this form has the advantage that sum rules apply, i.e. the trace of the occupation matrix sums up to the total number of electrons.
The dual representation corresponds to projecting into orthogonalized orbitals. Such orbitals may have nonzero weights also on neighboring centers, and our experience shows that this kind of nonlocality can result in nonphysical results in some cases. The default in QuantumATK is therefore to use the onsite representation.
The onsite shell and dual shell representations often give a better description with multiplezeta basis sets compared to the onsite and dual representation, where the occupation is obtained by projecting into single basis functions.
Input format for XC+U¶
The Hubbard U approximation in the onsite representation may be selected through
keywords of the types LDAU.XCTYPE
, LSDAU.XCTYPE
, GGAU.XCTYP
, and SGGAU.XCTYPE
.
The keyword LSDAU.PZ
gives the default LSDA functional with the default Hubbard U method.
The specification
exchange_correlation = LSDAU.PZ
is therefore identical to
exchange_correlation = ExchangeCorrelation(
exchange=DiracBloch,
correlation=PerdewZunger,
hubbard_term=Onsite,
number_of_spins=2,
)
Use this form with the option hubbard_term=Dual
to select the dual representation.
The value of the U term is specified through the basis set:
basis_set = [LDABasis.Nickel_SingleZeta(hubbard_u=[4.6, 0.0]*eV,
filling_method=Anisotropic)]
which will give a singlezeta basis set for nickel, consisting of a 3d and a 4s
orbital, where the 3d orbital has U=4.6 eV. The 3d orbitals are filled
anisotropically, i.e, there is 1 electron in in orbitals with m=2,1,0,1,
leaving m=2 empty. The alternative is filling_method=SphericalSymmetric
,
which puts 4/5 electron in each orbital.
DFT1/2 method¶
The DFT1/2 method (often also denoted LDA1/2 or GGA1/2) is a semiempirical approach to correct the selfinteraction error in local and semilocal exchangecorrelation functionals for extended systems. As such, it can be viewed as an alternative to the XC+U meanfield Hubbard term with broadly the same aims: to improve the description of conductionband energy levels and band gaps.
This method, introduced by Ferreira et al. [FMT08], is based on the much older Slater halfoccupation scheme for molecules [SJ72] (also known as the transition state method). Slater’s original method consists in carrying out a selfconsistent calculation with half an electron removed from the system, and taking the eigenvalue of the halffilled state as an estimate for the ionization energy. This approach was later formalized in the theorem by Janak [Jan78]:
where \(E\) is the total energy of the system, \(f_\alpha\) is the occupation of state \(\alpha\) (between 0 and 1), and \(\varepsilon_\alpha\) is the eigenvalue of the state. The success of the halfoccupation scheme is based on the fact that \(\varepsilon_\alpha \left ( f_\alpha \right )\) is known to be almost precisely linear for many cases, which makes the relationship with the ionization energy exact.
The DFT1/2 method makes use of the theoretical insights from the halfoccupation scheme and to tackle the fundamental problem of the selfinteraction error for the case of extended systems with a band gap. The method attempts to correct for this by defining an atomic selfenergy potential which cancels the electronhole selfinteraction energy. This potential is calculated for atomic sites in the system, and is defined as the difference between the potential of the neutral atom and that of a charged ion resulting from the removal of a fraction of its charge, between 0 and 1 electrons. The total selfenergy potential is the sum of these atomic potentials.
The addition of the DFT1/2 selfenergy potential to the DFT Hamiltonian has been found to greatly improve band gaps for a wide range of semiconducting and insulating systems [FMT11].
It is important to note that the method is not entirely free of empirical parameters. Much like the value of \(U\) for the Hubbard U method, there are two values which much be fixed for every species in the system: the fractional charge \(f_\alpha\) removed from the neutral atom, and a cutoff radius \(r^\mathrm{cut}\) beyond which to trim the atomic selfenergy potential. The latter is needed to avoid excessive overlap of the potentials from different sites in the crystal.
The original authors of the method suggest to fix \(r^\mathrm{cut}\) variationally, by choosing the value which maximizes the band gap [FMT08] [FMT11].
The choice of \(f_\alpha\) is less welldefined: although a value of 0.5 is standard (hence the name of the method), this is known not to be the optimal choice for various materials (e.g., silicon). It may therefore be appropriate to treat it as an empirical parameter, to be varied by comparison with experiment.
It is also important to note that not all species in the system necessarily require the DFT1/2 correction; it is generally advisable only to add this to the anionic species, and leave the cationic species as normal [FMT08] [FMT11].
Input format for DFT1/2¶
The DFT1/2 method may be selected by specifying one of the predefined
exchangecorrelation functionals: LDAHalf
, LSDAHalf
, NCLDAHalf
,
SOLDAHalf
, GGAHalf
, SGGAHalf
, NCGGAHalf
, SOGGAHalf
.
The keyword LDAHalf.PZ
gives the default LDA functional with the DFT1/2
correction. The specification
exchange_correlation = LDAHalf.PZ
is therefore identical to
exchange_correlation = ExchangeCorrelation(
exchange=DiracBloch,
correlation=PerdewZunger,
dft_half_enabled=True
)
In general, the keyword dft_half_enabled=True
can be specified for any
userdefined exchangecorrelation functional.
The value of the DFT1/2 parameters are specified through the basis set:
dft_half_parameters = DFTHalfParameters(
element=Arsenic,
fractional_charge=[0.3, 0.0],
cutoff_radius=4.0*Bohr)
basis_set = [
LDABasis.Arsenic_DoubleZetaPolarized(
dft_half_parameters=dft_half_parameters),
LDABasis.Gallium_DoubleZetaPolarized(
dft_half_parameters=Disabled)
]
The above example shows a typical case for gallium arsenide, in which the anion
(As) is given a DFT1/2 selfenergy potential calculated with
\(f_\alpha = 0.3\) and \(r^\mathrm{cut} = 4~a_0\), while the cation (Ga)
has its DFT1/2 correction turned off with the keyword Disabled
.
If dft_half_parameters
is not specified for a species, the default
Automatic
is used. This will set the DFT1/2 default parameters
from the table of optimized values given below.
Warning
When the DFT1/2 method is enabled, the selfenergy potential is included in the calculation of the total energy, forces and stress. However, the method is only intended as a scheme for calculating band structures, not for structural relaxations. We recommend that structural properties be calculated before applying the DFT1/2 correction.
DFT1/2 default parameters¶
The default DFT1/2 parameters used by the Automatic
keyword have been
optimized against a wide range of materials, and should improve upon the
standard DFT band gap in most cases. The complete list is given in the table
below. As explained above, we only include elements which tend to have anionic
character in compounds of interest.
Important
All elements not present in the table default to Disabled
, i.e., no
DFT1/2 correction.
Important
For metaGGA exchangecorrelation functionals, all elements default to
Disabled
.
Element  LDA functionals  GGA functionals  

fractional_charge 
cutoff_radius 
fractional_charge 
cutoff_radius 

C  0.3  2.5 Bohr  0.4  2.5 Bohr 
N  0.4  3.0 Bohr  0.4  3.0 Bohr 
O  0.5  2.5 Bohr  0.4  3.0 Bohr 
F  0.7  3.0 Bohr  0.6  2.5 Bohr 
Si  0.2  4.0 Bohr  0.2  4.0 Bohr 
P  0.3  4.0 Bohr  0.3  3.5 Bohr 
S  0.6  3.5 Bohr  0.4  3.5 Bohr 
Cl  0.7  3.5 Bohr  0.7  3.5 Bohr 
Zn  0.5  2.0 Bohr  0.5  2.0 Bohr 
Ge  0.5  4.0 Bohr  0.6  4.0 Bohr 
As  0.3  4.0 Bohr  0.3  4.0 Bohr 
Se  0.5  3.5 Bohr  0.5  3.5 Bohr 
Br  0.7  4.0 Bohr  0.7  4.0 Bohr 
Sn  0.7  5.0 Bohr  0.9  4.5 Bohr 
Sb  0.3  4.5 Bohr  0.3  4.5 Bohr 
Te  0.4  4.0 Bohr  0.4  4.0 Bohr 
I  0.6  4.5 Bohr  0.6  4.0 Bohr 
Spinorbit coupling (SOC)¶
Background¶
The spinorbit coupling can lead to energy level splits of bandstructures and molecular energy spectra. As a rule of thumb, the magnitude of SOC scales as \(Z^{4}/n^{3}l^{2}\), where \(Z\), \(n\), and \(l\) are the atomic charge, principle, and angular quantum number of the considered energy level, respectively. The SOC is fundamentally a relativistic effect, and can be described by a secondorder expansion of the relativistic Hamiltonian in terms of the fine structure constant \(\alpha\) [FernandezSOSF06].
In DFT: LCAO, the SOC is taken into account via the NormConservingPseudoPotential. For example, the SG15 pseudopotentials, shipped with QuantumATK as of the 2016 release, come in two versions for each element; a standard scalarrelativistic pseudopotential (denoted “SG15”), and one that is generated by mapping the solution to the Dirac equation, which naturally includes the SOC, to a scalarrelativistic pseudopotential (denoted “SG15SO”):
with a local contribution \(V_{L}\) and nonlocal contributions from total angular momenta \(j=l+1/2\) and \(j=l1/2\).
Each nonlocal term is expanded in spinorbit projector functions \(P_{\alpha\beta}^{l \pm 1/2, \xi}\),
where \(\nu_{l\pm 1/2,\xi}\) are normalization constants. The indices \(\alpha,\beta\) denote the possible spin orientations (up, down). Each nonlocal pseudopotential term requires four projector functions per expansion order \(\xi\).
Tip
Use SG15SO or OpenMX pseudopotentials for DFT: LCAO calculations with spinorbit.
Usage¶
The SOC terms are enabled by choosing a pseudopotential containing spinorbit terms (SG15SO or OpenMX), and appropriate settings for the ExchangeCorrelation keyword. For example, an QuantumATK Python script for silicon with SOC could contain the following specification for the calculator:
#
# Basis Set
#
basis_set = [
BasisGGASG15SO.Silicon_Medium,
]
#
# ExchangeCorrelation
#
exchange_correlation = SOGGA.PBE
#
# Calculator
#
k_point_sampling = MonkhorstPackGrid(na=9,nb=9,nc=9)
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=k_point_sampling,
density_mesh_cutoff=100.0*Hartree,
)
calculator = LCAOCalculator(
basis_set=basis_set,
exchange_correlation=exchange_correlation,
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
Note the setting exchange_correlation=SOGGA.PBE
, which switches on
the SOC terms in the pseudopotential. Omitting this keyword or setting
it to a nonSOC method, e.g. NCGGA.PBE
, would still result in a calculation
with the SG15SO pseudopotential, but the SOC terms would be disregarded.
[CdG05]  (1, 2) M. Cococcioni and S. de Gironcoli. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B, 71:035105, Jan 2005. doi:10.1103/PhysRevB.71.035105. 
[DBS+98]  S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton. Electronenergyloss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B, 57:1505–1509, Jan 1998. doi:10.1103/PhysRevB.57.1505. 
[FernandezSOSF06]  L. FernándezSeivane, M. A. Oliveira, S. Sanvito, and J. Ferrer. Onsite approximation for spin–orbit coupling in linear combination of atomic orbitals density functional methods. J. Phys.: Condensed Matter, 18(34):7999, 2006. URL: http://stacks.iop.org/09538984/18/i=34/a=012. 
[FMT08]  (1, 2, 3) Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Approximation to density functional theory for the calculation of band gaps of semiconductors. Phys. Rev. B, 78:125116, Sep 2008. doi:10.1103/PhysRevB.78.125116. 
[FMT11]  (1, 2, 3) Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Slater halfoccupation technique revisited: the LDA1/2 and GGA1/2 approaches for atomic ionization energies and band gaps in semiconductors. AIP Adv., 1(3):032119, 2011. doi:10.1063/1.3624562. 
[HOY06]  M. J. Han, T. Ozaki, and J. Yu. O(N) LDA+U electronic structure calculation method based on the nonorthogonal pseudoatomic orbital basis. Phys. Rev. B, 73:045110, Jan 2006. doi:10.1103/PhysRevB.73.045110. 
[Jan78]  J. F. Janak. Proof that $\partial E / \partial n_i = \varepsilon \,$ in densityfunctional theory. Phys. Rev. B, 18:7165–7168, Dec 1978. doi:10.1103/PhysRevB.18.7165. 
[SJ72]  J. C. Slater and K. H. Johnson. Selfconsistentfield $X \alpha \,$ cluster method for polyatomic molecules and solids. Phys. Rev. B, 5:844–853, Feb 1972. doi:10.1103/PhysRevB.5.844. 
[SAG+02]  J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. SánchezPortal. The SIESTA method for ab initio orderN materials simulation. J. Phys.: Condensed Matter, 14(11):2745, 2002. URL: http://stacks.iop.org/09538984/14/i=11/a=302. 