Range separated functionals (RSF)¶
Introduction¶
Range separated functionals (RSF) are a subgroup of hybrid functionals. While conventional (global) hybrid functionals like PBE0 or B3LYP use fixed fractions of HartreeFock (HFT, E_{XX}) and DFT (E_{X}) exchange for exchange, f.e. 1/4 E_{XX} and 3/4 E_{X} in the case of PBE0, RSFs mix the two contributions by the spatial distance between two points, \(r_{12}\), using a soft function \(\omega_\mathrm{RSF}(\gamma, r_{12})\).
To achieve this, the coulomb interaction kernel, \(\frac{1}{r_{12}} = \frac{1}{r_1  r_2}\), which appears in the exchange integral from HFT is split into two parts:
\(\frac{1}{r_{12}} = \underbrace{\frac{1  [\alpha + \beta ( 1  \omega_\mathrm{RSF} (\gamma, r_{12}))]}{r_{12}}}_{\text{SR, DFT}} + \underbrace{\frac{\alpha + \beta ( 1  \omega_\mathrm{RSF} (\gamma, r_{12}))}{r_{12}}}_{\text{LR, HFT}}\),
the shortrange (SR) part is handled by the exchange from a (semi)local LDA
or GGA functional such as PBE, while the longrange part (LR) is handled by
the exchange from HFT. \(\alpha\) and \(\beta\) are functional dependent mixing
parameters. \(\alpha \ne 0\) and \(\beta = 0\) resembles conventional global
hybrids. RSFs with \(\alpha = 0\) and \(\beta \ne 0\) are usually denoted by
LC
and the name of the semilocal functional, f.e. LCPBE.
RSFs with \(\alpha \ne 0\) and \(\beta \ne 0\) are usually denoted by
CAM
and the name of the semilocal functional, f.e. CAMBLYP.
For the separating function \(\omega_\mathrm{RSF}\), two functions are in common use: either the complementary error function, \(\omega_\mathrm{RSF} = \mathrm{erfc}(\gamma r_{12})\), or the Slaterfunction, \(\omega_\mathrm{RSF} = e^{(\gamma r_{12})}\). While the use of the complementary error function is computationally fortunate for codes utilizing Gaussian type basis sets, the Slaterfunction give superior results in the calculation of Rydbergstate and charge transfer excitations. To distinguish between these both functions, functionals using the Slaterfunction append the letter “Y” to the RSF marker, f.e. LCYPBE or CAMYB3LYP, while functionals using the complementary error function keep the marker as it is, f.e. LCPBE or CAMB3LYP.
Besides \(r_{12}\), the separation functions use a second parameter, the screening factor \(\gamma\). The optional value for \(\gamma\) is under discussion. A density dependence is stated. For most RSF standard values for \(\gamma\) are defined, although it is possible to tune \(\gamma\) to optimal values for calculations investigating ionization potentials, charge transfer excitations and the binding curves of biradical cations.
Implementation¶
 The implementation of RSFs in gpaw consists of two parts:
 once the implementation of the semilocal functional part. This is done in libxc.
 once the implementation of the HartreeFock exchange. This is done
in
hybrid.py
.
As range separating function the Slaterfunction, \(\omega_\mathrm{RSF} = e^{(\gamma r_{12})}\), is used. Besides the possibility to set \(\gamma\) to an arbitrary value, the following functionals were implemented:
Functional  \(\alpha\)  \(\beta\)  \(\gamma\) (\(a_0^{1}\))  Reference 

CAMYBLYP  0.2  0.8  0.44  [AT08] 
CAMYB3LYP  0.19  0.46  0.34  [SZ12] 
LCYBLYP  0.0  1.0  0.75  [SZ12] 
LCYPBE  0.0  1.0  0.75  [SZ12] 
As the implementation of RSFs in gpaw is based on the finite difference exact exchange code (hybrid.py), the implementation inherits its positive and negative properties, in summary:
 selfconsistent calculations using RSFs
 calculations can only be done for the \(\Gamma\) point
 only nonperiodic boundary conditions can be used
 only RMMDIIS can be used as eigensolver
Important: As one of the major benefits of the RSF is to retain the \(\frac{1}{r}\) asymptote of the exchange potential, one has to use large boxes if neutral or anionic systems are considered. Large boxes start at 6Å vacuum around each atom. For anionic systems “large” should be extended.
Further information about the implementation and RSFs can be found in [WW18] and in detail in [Wu16].
Simple usage¶
In general calculations using RSF can simply be done choosing the appropriate functional as in the following snippet:
"""First example for using RSF."""
from ase import Atoms
from gpaw.eigensolvers import RMMDIIS
from gpaw.occupations import FermiDirac
h = 0.30
co = Atoms('CO', positions=[(0, 0, 0), (0, 0, 1.15)])
co.center(5)
# c = {'energy': 0.005, 'eigenstates': 1e4} # Usable values
c = {'energy': 0.1, 'eigenstates': 3, 'density': 3} # Values for test
calc = GPAW(txt='CO.txt', xc='LCYPBE', convergence=c,
eigensolver=RMMDIIS(), h=h,
occupations=FermiDirac(width=0.0), spinpol=False)
co.set_calculator(calc)
co.get_potential_energy()
Three main points can be seen already in this small snippet. Even if choosing
the RSF is quite simple by choosing xc=LCYPBE
, one has to choose RMMDIIS
as eigensolver, eigensolver=RMMDIIS()
, and has to decrease the
convergence criteria a little.
Improving results¶
However, there are a few drawbacks, at first in an SCF calculation the
contributions from the core electrons are also needed, which have to be
calculated during the generation of the PAW datasets. Second: for the
calculation of the exchange on the Cartesian grid, the (screened) Poisson
equation has to be solved numerically. For a charged system, as f.e. the
exchange of a state with itself, one has to neutralize the charge by
subtracting a Gaussian representing the “overcharge”, solve the
(screened) Poissonequation for the neutral system and add the solution
for the Gaussian to the solution for the neutral system. However, if the
charge to remove is “offcenter”, the center of the neutralizing charge
should match the center of the “overcharge”
preventing an artificial dipole. The latter is done by using a Poisson solver
which uses the charge center for removal:
poissonsolver=PoissonSolver(use_charge_center=True)
.
The next listing shows these two steps:
"""Calculations using RSF with dedicated datasets and Poissonsolver."""
from ase import Atoms
from gpaw import GPAW, setup_paths
from gpaw.poisson import PoissonSolver
from gpaw.eigensolvers import RMMDIIS
from gpaw.occupations import FermiDirac
from gpaw.test import gen
if setup_paths[0] != '.':
setup_paths.insert(0, '.')
for atom in ['C', 'O']:
gen(atom, xcname='PBE', scalarrel=True, exx=True,
yukawa_gamma=0.75)
h = 0.30
co = Atoms('CO', positions=[(0, 0, 0), (0, 0, 1.15)])
co.center(5)
# c = {'energy': 0.005, 'eigenstates': 1e4} # Usable values
c = {'energy': 0.1, 'eigenstates': 3, 'density': 3} # Values for test
calc = GPAW(txt='CO.txt', xc='LCYPBE', convergence=c,
eigensolver=RMMDIIS(), h=h,
poissonsolver=PoissonSolver(use_charge_center=True),
occupations=FermiDirac(width=0.0), spinpol=False)
co.set_calculator(calc)
co.get_potential_energy()
The generation of PAWdatasets can also be done by
gpawsetup f PBE x gamma=0.75 C O
Tuning \(\gamma\)¶
As stated in the introduction, the optimal value for \(\gamma\) is under
discussion. One way to find the optimal value for \(\gamma\) for ionization
potentials is to tune \(\gamma\) in a way, that the negative eigenvalue of the
HOMO matches the calculated IP. To use different values of \(\gamma\), one has
to pass the desired value of \(\gamma\) to the variable omega
.
"""Calculation utilizing RSF with optimized gamma."""
from ase import Atoms
from gpaw import GPAW, setup_paths
from gpaw.poisson import PoissonSolver
from gpaw.eigensolvers import RMMDIIS
from gpaw.occupations import FermiDirac
from gpaw.test import equal, gen
# IP for CO using LCYPBE with gamma=0.81 after
# dx.doi.org/10.1021/acs.jctc.8b00238
IP = 14.31
if setup_paths[0] != '.':
setup_paths.insert(0, '.')
for atom in ['C', 'O']:
gen(atom, xcname='PBE', scalarrel=True, exx=True,
yukawa_gamma=0.81)
h = 0.30
co = Atoms('CO', positions=[(0, 0, 0), (0, 0, 1.15)])
co.minimal_box(5)
# c = {'energy': 0.005, 'eigenstates': 1e4} # Usable values
c = {'energy': 0.1, 'eigenstates': 3, 'density': 3} # Values for test
calc = GPAW(txt='CO.txt', xc='LCYPBE:omega=0.81', convergence=c,
eigensolver=RMMDIIS(), h=h,
poissonsolver=PoissonSolver(use_charge_center=True),
occupations=FermiDirac(width=0.0), spinpol=False)
co.set_calculator(calc)
co.get_potential_energy()
(eps_homo, eps_lumo) = calc.get_homo_lumo()
equal(eps_homo, IP, 0.15)
linear response TDDFT¶
One of the major benefits of RSF is their ability to describe longrange charge transfer by
linear response timedependent DFT (lrTDDFT). If one uses RSF with lrTDDFT one has at least
to activate the use of the Fock operator (FO) on the unoccupied states. Also the charge
centered compensation of the over charge should be activated, see [Wu16] for details.
The use of the FO on the unoccupied states is activated by the keyword unocc=True
as in
the following code:
"""Check TDDFT ionizations with Yukawa potential."""
from ase.structure import molecule
from ase.units import Hartree
from gpaw import GPAW
from gpaw.mpi import world
from gpaw.cluster import Cluster
from gpaw.occupations import FermiDirac
from gpaw.test import equal
from gpaw.eigensolvers import RMMDIIS
from gpaw.lrtddft import LrTDDFT
h2o = Cluster(molecule('H2O'))
h2o.set_initial_magnetic_moments([2, 1, 1])
h2o.minimal_box(3.0, h=0.3)
h2o_plus = Cluster(molecule('H2O'))
h2o_plus.set_initial_magnetic_moments([2, 0.5, 0.5])
h2o_plus.minimal_box(3.0, h=0.3)
def get_paw():
"""Return calculator object."""
c = {'energy': 0.001, 'eigenstates': 0.001, 'density': 0.001}
return GPAW(convergence=c, eigensolver=RMMDIIS(),
xc='LCYPBE:omega=0.83:unocc=True',
parallel={'domain': world.size}, h=0.3,
occupations=FermiDirac(width=0.0, fixmagmom=True))
calc = get_paw()
calc.set(txt='H2O_LCY_PBE_083.log')
calc_plus = get_paw()
calc_plus.set(txt='H2O_plus_LCY_PBE_083.log', charge=1)
h2o.set_calculator(calc)
e_h2o = h2o.get_potential_energy()
h2o_plus.set_calculator(calc_plus)
e_h2o_plus = h2o_plus.get_potential_energy()
e_ion = e_h2o_plus  e_h2o
print(e_ion, 12.62)
equal(e_ion, 12.62, 0.1)
lr = LrTDDFT(calc_plus, txt='LCY_TDDFT_H2O.log', jend=4)
equal(lr.xc.omega, 0.83)
lr.write('LCY_TDDFT_H2O.ex.gz')
# reading is problematic with EXX on more than one core
if world.rank == 0:
lr2 = LrTDDFT('LCY_TDDFT_H2O.ex.gz')
lr2.diagonalize()
equal(lr2.xc.omega, 0.83)
for i, ip_i in enumerate([14.74, 18.51]):
ion_i = lr2[i].get_energy() * Hartree + e_ion
print(ion_i, ip_i)
equal(ion_i, ip_i, 0.6)
[AT08] 

[SZ12]  (1, 2, 3)

[Wu16]  (1, 2)

[WW18] 
