12. dft — Density functional theory

12.1. Customizing XC functional

XC functional of DFT methods can be customized. The simplest way to customize XC functional is to assigned a string expression to mf.xc:

from pyscf import gto, dft
mol = gto.M(atom='H  0  0  0; F  0.9  0  0', basis='6-31g')
mf = dft.RKS(mol)
mf.xc = 'HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN'
mf.kernel()
mf.xc = 'HF*0.5 + .08*LDA + .42*B88, .81*LYP + .19*VWN'
mf.kernel()
mf.xc = 'HF*0.8 + .08*LDA + .12*B88, .81*LYP + .19*VWN'
mf.kernel()
mf.xc = 'HF'
mf.kernel()

The XC functional string is parsed against the rules, as described below.

  • The given functional description must be a one-line string.

  • The functional description is case-insensitive.

  • The functional description string has two parts, separated by ,. The first part describes the exchange functional, the second is the correlation functional. - If , was not appeared in string, the entire string is considered as

    X functional.

    • To neglect X functional (just apply C functional), leave blank in the first part, eg mf.xc=',vwn' for pure VWN functional
  • The functional name can be placed in arbitrary order. Two names needs to be separated by operators + or -. Blank spaces are ignored. NOTE the parser only reads operators + - *. / is not supported.

  • A functional name is associated with one factor. If the factor is not given, it is assumed equaling 1.

  • String 'HF' stands for exact exchange (HF K matrix). It is allowed to put 'HF' in C (correlation) functional part.

  • Be careful with the libxc convention on GGA functional, in which the LDA contribution is included.

Another way to customize XC functional is to redefine the eval_xc() method of numerical integral class:

mol = gto.M(atom='H 0 0 0; F 0.9 0 0', basis = '6-31g')
mf = dft.RKS(mol)
def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, verbose=None):
    # A fictitious XC functional to demonstrate the usage
    rho0, dx, dy, dz = rho
    gamma = (dx**2 + dy**2 + dz**2)
    exc = .01 * rho0**2 + .02 * (gamma+.001)**.5
    vrho = .01 * 2 * rho0
    vgamma = .02 * .5 * (gamma+.001)**(-.5)
    vlapl = None
    vtau = None
    vxc = (vrho, vgamma, vlapl, vtau)
    fxc = None  # 2nd order functional derivative
    kxc = None  # 3rd order functional derivative
    return exc, vxc, fxc, kxc
dft.libxc.define_xc_(mf._numint, eval_xc, xctype='GGA')
mf.kernel()

By calling dft.libxc.define_xc_() function, the customized eval_xc() function is patched to the numerical integration class mf._numint dynamically.

More examples of customizing DFT XC functional can be found in examples/dft/24-custom_xc_functional.py and examples/dft/24-define_xc_functional.py.

12.2. Program reference

Non-relativistic restricted Kohn-Sham

class pyscf.dft.rks.RKS(mol)[source]

Restricted Kohn-Sham SCF base class. non-relativistic RHF.

Attributes:
verbose
: int
Print level. Default value equals to Mole.verbose
max_memory
: float or int
Allowed memory in MB. Default equals to Mole.max_memory
chkfile
: str
checkpoint file to save MOs, orbital energies etc. Writing to chkfile can be disabled if this attribute is set to None or False.
conv_tol
: float
converge threshold. Default is 1e-9
conv_tol_grad
: float
gradients converge threshold. Default is sqrt(conv_tol)
max_cycle
: int
max number of iterations. Default is 50
init_guess
: str
initial guess method. It can be one of ‘minao’, ‘atom’, ‘hcore’, ‘1e’, ‘chkfile’. Default is ‘minao’
DIIS
: DIIS class
The class to generate diis object. It can be one of diis.SCF_DIIS, diis.ADIIS, diis.EDIIS.
diis
: boolean or object of DIIS class defined in scf.diis.
Default is the object associated to the attribute self.DIIS. Set it to None/False to turn off DIIS. Note if this attribute is inialized as a DIIS object, the SCF driver will use this object in the iteration. The DIIS informations (vector basis and error vector) will be held inside this object. When kernel function is called again, the old states (vector basis and error vector) will be reused.
diis_space
: int
DIIS space size. By default, 8 Fock matrices and errors vector are stored.
diis_start_cycle
: int
The step to start DIIS. Default is 1.
diis_file: ‘str’
File to store DIIS vectors and error vectors.
level_shift
: float or int
Level shift (in AU) for virtual space. Default is 0.
direct_scf
: bool
Direct SCF is used by default.
direct_scf_tol
: float
Direct SCF cutoff threshold. Default is 1e-13.
callback
: function(envs_dict) => None
callback function takes one dict as the argument which is generated by the builtin function locals(), so that the callback function can access all local variables in the current envrionment.
conv_check
: bool
An extra cycle to check convergence after SCF iterations.
check_convergence
: function(envs) => bool
A hook for overloading convergence criteria in SCF iterations.

Saved results:

converged
: bool
SCF converged or not
e_tot
: float
Total HF energy (electronic energy plus nuclear repulsion)
mo_energy :
Orbital energies
mo_occ
Orbital occupancy
mo_coeff
Orbital coefficients

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> mf = scf.hf.SCF(mol)
>>> mf.verbose = 0
>>> mf.level_shift = .4
>>> mf.scf()
-1.0811707843775884
Attributes for RKS:
xc
: str
‘X_name,C_name’ for the XC functional. Default is ‘lda,vwn’
nlc
: str
‘NLC_name’ for the NLC functional. Default is ‘’ (i.e., None)
omega
: float
Omega of the range-separated Coulomb operator e^{-omega r_{12}^2} / r_{12}
grids
: Grids object

grids.level (0 - 9) big number for large mesh grids. Default is 3

radii_adjust
radi.treutler_atomic_radii_adjust (default)
radi.becke_atomic_radii_adjust
None : to switch off atomic radii adjustment
grids.atomic_radii
radi.BRAGG_RADII (default)
radi.COVALENT_RADII
None : to switch off atomic radii adjustment
grids.radi_method scheme for radial grids
radi.treutler (default)
radi.delley
radi.mura_knowles
radi.gauss_chebyshev
grids.becke_scheme weight partition function
gen_grid.original_becke (default)
gen_grid.stratmann
grids.prune scheme to reduce number of grids
gen_grid.nwchem_prune (default)
gen_grid.sg1_prune
gen_grid.treutler_prune
None : to switch off grids pruning

grids.symmetry True/False to symmetrize mesh grids (TODO)

grids.atom_grid Set (radial, angular) grids for particular atoms. Eg, grids.atom_grid = {‘H’: (20,110)} will generate 20 radial grids and 110 angular grids for H atom.

small_rho_cutoff
: float
Drop grids if their contribution to total electrons smaller than this cutoff value. Default is 1e-7.

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', verbose=0)
>>> mf = dft.RKS(mol)
>>> mf.xc = 'b3lyp'
>>> mf.kernel()
-76.415443079840458
energy_elec(ks, dm=None, h1e=None, vhf=None)

Electronic part of RKS energy.

Args:

ks : an instance of DFT class

dm
: 2D ndarray
one-partical density matrix
h1e
: 2D ndarray
Core hamiltonian
Returns:
RKS electronic energy and the 2-electron part contribution
get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)

Coulomb + XC functional

Note

This function will change the ks object.

Args:
ks
: an instance of RKS
XC functional are controlled by ks.xc attribute. Attribute ks.grids might be initialized.
dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference Vxc potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
Returns:
matrix Veff = J + Vxc. Veff can be a list matrices, if the input dm is a list of density matrices.
pyscf.dft.rks.energy_elec(ks, dm=None, h1e=None, vhf=None)[source]

Electronic part of RKS energy.

Args:

ks : an instance of DFT class

dm
: 2D ndarray
one-partical density matrix
h1e
: 2D ndarray
Core hamiltonian
Returns:
RKS electronic energy and the 2-electron part contribution
pyscf.dft.rks.get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Coulomb + XC functional

Note

This function will change the ks object.

Args:
ks
: an instance of RKS
XC functional are controlled by ks.xc attribute. Attribute ks.grids might be initialized.
dm
: ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last
: ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last
: ndarray or a list of ndarrays or 0
The reference Vxc potential matrix.
hermi
: int

Whether J, K matrix is hermitian

0 : no hermitian or symmetric
1 : hermitian
2 : anti-hermitian
Returns:
matrix Veff = J + Vxc. Veff can be a list matrices, if the input dm is a list of density matrices.

Non-relativistic Unrestricted Kohn-Sham

class pyscf.dft.uks.UKS(mol)[source]

Unrestricted Kohn-Sham See pyscf/dft/rks.py RKS class for document of the attributes

get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)

Coulomb + XC functional for UKS. See pyscf/dft/rks.py get_veff() fore more details.

pyscf.dft.uks.get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1)[source]

Coulomb + XC functional for UKS. See pyscf/dft/rks.py get_veff() fore more details.

Generate DFT grids and weights, based on the code provided by Gerald Knizia <>

Reference for Lebedev-Laikov grid:
V. I. Lebedev, and D. N. Laikov “A quadrature formula for the sphere of the 131st algebraic order of accuracy”, Doklady Mathematics, 59, 477-481 (1999)
class pyscf.dft.gen_grid.Grids(mol)[source]

DFT mesh grids

Attributes for Grids:
level
: int (0 - 9)
big number for large mesh grids, default is 3
atomic_radii
: 1D array
radi.BRAGG_RADII (default)
radi.COVALENT_RADII
None : to switch off atomic radii adjustment
radii_adjust
: function(mol, atomic_radii) => (function(atom_id, atom_id, g) => array_like_g)
Function to adjust atomic radii, can be one of | radi.treutler_atomic_radii_adjust | radi.becke_atomic_radii_adjust | None : to switch off atomic radii adjustment
radi_method
: function(n) => (rad_grids, rad_weights)
scheme for radial grids, can be one of | radi.treutler (default) | radi.delley | radi.mura_knowles | radi.gauss_chebyshev
becke_scheme
: function(v) => array_like_v
weight partition function, can be one of | gen_grid.original_becke (default) | gen_grid.stratmann
prune
: function(nuc, rad_grids, n_ang) => list_n_ang_for_each_rad_grid
scheme to reduce number of grids, can be one of | gen_grid.nwchem_prune (default) | gen_grid.sg1_prune | gen_grid.treutler_prune | None : to switch off grid pruning
symmetry
: bool
whether to symmetrize mesh grids (TODO)
atom_grid
: dict
Set (radial, angular) grids for particular atoms. Eg, grids.atom_grid = {‘H’: (20,110)} will generate 20 radial grids and 110 angular grids for H atom.
level
: int
To control the number of radial and angular grids. The default level 3 corresponds to (50,302) for H, He; (75,302) for second row; (80~105,434) for rest.

Examples:

>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.level = 4
>>> grids.build()
gen_atomic_grids(mol, atom_grid=None, radi_method=None, level=None, prune=None, **kwargs)[source]

Generate number of radial grids and angular grids for the given molecule.

Returns:
A dict, with the atom symbol for the dict key. For each atom type, the dict value has two items: one is the meshgrid coordinates wrt the atom center; the second is the volume of that grid.
gen_partition(mol, atom_grids_tab, radii_adjust=None, atomic_radii=array([ 0., 0.66140414, 2.64561657, 2.74010288, 1.98421243, 1.60626721, 1.32280829, 1.22832198, 1.13383567, 0.94486306, 2.83458919, 3.40150702, 2.83458919, 2.36215766, 2.07869874, 1.88972612, 1.88972612, 1.88972612, 3.40150702, 4.15739747, 3.40150702, 3.0235618, 2.64561657, 2.55113027, 2.64561657, 2.64561657, 2.64561657, 2.55113027, 2.55113027, 2.55113027, 2.55113027, 2.45664396, 2.36215766, 2.17318504, 2.17318504, 2.17318504, 3.59047964, 4.44085639, 3.77945225, 3.40150702, 2.92907549, 2.74010288, 2.74010288, 2.55113027, 2.45664396, 2.55113027, 2.64561657, 3.0235618, 2.92907549, 2.92907549, 2.74010288, 2.74010288, 2.64561657, 2.64561657, 3.96842486, 4.91328792, 4.06291117, 3.68496594, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 2.92907549, 2.74010288, 2.55113027, 2.55113027, 2.45664396, 2.55113027, 2.55113027, 2.55113027, 2.83458919, 3.59047964, 3.40150702, 3.0235618, 3.59047964, 2.74010288, 3.96842486, 3.40150702, 4.06291117, 3.68496594, 3.40150702, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072]), becke_scheme=<function original_becke>)[source]

Generate the mesh grid coordinates and weights for DFT numerical integration. We can change radii_adjust, becke_scheme functions to generate different meshgrid.

Returns:
grid_coord and grid_weight arrays. grid_coord array has shape (N,3); weight 1D array has N elements.
make_mask(mol=None, coords=None, relativity=0, shls_slice=None, verbose=None)[source]

Mask to indicate whether a shell is zero on grid

Args:

mol : an instance of Mole

coords
: 2D array, shape (N,3)
The coordinates of grids.
Kwargs:
relativity
: bool
No effects.
shls_slice
: 2-element list
(shl_start, shl_end). If given, only part of AOs (shl_start <= shell_id < shl_end) are evaluated. By default, all shells defined in mol will be evaluated.
verbose
: int or object of Logger
No effects.
Returns:
2D mask array of shape (N,nbas), where N is the number of grids, nbas is the number of shells.
pyscf.dft.gen_grid.gen_atomic_grids(mol, atom_grid={}, radi_method=<function gauss_chebyshev>, level=3, prune=<function nwchem_prune>, **kwargs)[source]

Generate number of radial grids and angular grids for the given molecule.

Returns:
A dict, with the atom symbol for the dict key. For each atom type, the dict value has two items: one is the meshgrid coordinates wrt the atom center; the second is the volume of that grid.
pyscf.dft.gen_grid.gen_partition(mol, atom_grids_tab, radii_adjust=None, atomic_radii=array([ 0., 0.66140414, 2.64561657, 2.74010288, 1.98421243, 1.60626721, 1.32280829, 1.22832198, 1.13383567, 0.94486306, 2.83458919, 3.40150702, 2.83458919, 2.36215766, 2.07869874, 1.88972612, 1.88972612, 1.88972612, 3.40150702, 4.15739747, 3.40150702, 3.0235618, 2.64561657, 2.55113027, 2.64561657, 2.64561657, 2.64561657, 2.55113027, 2.55113027, 2.55113027, 2.55113027, 2.45664396, 2.36215766, 2.17318504, 2.17318504, 2.17318504, 3.59047964, 4.44085639, 3.77945225, 3.40150702, 2.92907549, 2.74010288, 2.74010288, 2.55113027, 2.45664396, 2.55113027, 2.64561657, 3.0235618, 2.92907549, 2.92907549, 2.74010288, 2.74010288, 2.64561657, 2.64561657, 3.96842486, 4.91328792, 4.06291117, 3.68496594, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 2.92907549, 2.74010288, 2.55113027, 2.55113027, 2.45664396, 2.55113027, 2.55113027, 2.55113027, 2.83458919, 3.59047964, 3.40150702, 3.0235618, 3.59047964, 2.74010288, 3.96842486, 3.40150702, 4.06291117, 3.68496594, 3.40150702, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072]), becke_scheme=<function original_becke>)[source]

Generate the mesh grid coordinates and weights for DFT numerical integration. We can change radii_adjust, becke_scheme functions to generate different meshgrid.

Returns:
grid_coord and grid_weight arrays. grid_coord array has shape (N,3); weight 1D array has N elements.
pyscf.dft.gen_grid.make_mask(mol, coords, relativity=0, shls_slice=None, verbose=None)[source]

Mask to indicate whether a shell is zero on grid

Args:

mol : an instance of Mole

coords
: 2D array, shape (N,3)
The coordinates of grids.
Kwargs:
relativity
: bool
No effects.
shls_slice
: 2-element list
(shl_start, shl_end). If given, only part of AOs (shl_start <= shell_id < shl_end) are evaluated. By default, all shells defined in mol will be evaluated.
verbose
: int or object of Logger
No effects.
Returns:
2D mask array of shape (N,nbas), where N is the number of grids, nbas is the number of shells.
pyscf.dft.gen_grid.nwchem_prune(nuc, rads, n_ang, radii=array([ 0., 0.66140414, 2.64561657, 2.74010288, 1.98421243, 1.60626721, 1.32280829, 1.22832198, 1.13383567, 0.94486306, 2.83458919, 3.40150702, 2.83458919, 2.36215766, 2.07869874, 1.88972612, 1.88972612, 1.88972612, 3.40150702, 4.15739747, 3.40150702, 3.0235618, 2.64561657, 2.55113027, 2.64561657, 2.64561657, 2.64561657, 2.55113027, 2.55113027, 2.55113027, 2.55113027, 2.45664396, 2.36215766, 2.17318504, 2.17318504, 2.17318504, 3.59047964, 4.44085639, 3.77945225, 3.40150702, 2.92907549, 2.74010288, 2.74010288, 2.55113027, 2.45664396, 2.55113027, 2.64561657, 3.0235618, 2.92907549, 2.92907549, 2.74010288, 2.74010288, 2.64561657, 2.64561657, 3.96842486, 4.91328792, 4.06291117, 3.68496594, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.49599333, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 2.92907549, 2.74010288, 2.55113027, 2.55113027, 2.45664396, 2.55113027, 2.55113027, 2.55113027, 2.83458919, 3.59047964, 3.40150702, 3.0235618, 3.59047964, 2.74010288, 3.96842486, 3.40150702, 4.06291117, 3.68496594, 3.40150702, 3.40150702, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072, 3.30702072]))[source]

NWChem

Args:
nuc
: int
Nuclear charge.
rads
: 1D array
Grid coordinates on radical axis.
n_ang
: int
Max number of grids over angular part.
Kwargs:
radii
: 1D array
radii (in Bohr) for atoms in periodic table
Returns:
A list has the same length as rads. The list element is the number of grids over angular part for each radial grid.
pyscf.dft.gen_grid.original_becke(g)[source]

Becke, JCP, 88, 2547 (1988)

pyscf.dft.gen_grid.sg1_prune(nuc, rads, n_ang, radii=array([ 0., 1., 0.5882, 3.0769, 2.0513, 1.5385, 1.2308, 1.0256, 0.8791, 0.7692, 0.6838, 4.0909, 3.1579, 2.5714, 2.1687, 1.875, 1.6514, 1.4754, 1.3333]))[source]

SG1, CPL, 209, 506

Args:
nuc
: int
Nuclear charge.
rads
: 1D array
Grid coordinates on radical axis.
n_ang
: int
Max number of grids over angular part.
Kwargs:
radii
: 1D array
radii (in Bohr) for atoms in periodic table
Returns:
A list has the same length as rads. The list element is the number of grids over angular part for each radial grid.
pyscf.dft.gen_grid.stratmann(g)[source]

Stratmann, Scuseria, Frisch. CPL, 257, 213 (1996)

pyscf.dft.gen_grid.treutler_prune(nuc, rads, n_ang, radii=None)[source]

Treutler-Ahlrichs

Args:
nuc
: int
Nuclear charge.
rads
: 1D array
Grid coordinates on radical axis.
n_ang
: int
Max number of grids over angular part.
Returns:
A list has the same length as rads. The list element is the number of grids over angular part for each radial grid.
pyscf.dft.numint.cache_xc_kernel(ni, mol, grids, xc_code, mo_coeff, mo_occ, spin=0, max_memory=2000)[source]

Compute the 0th order density, Vxc and fxc. They can be used in TDDFT, DFT hessian module etc.

pyscf.dft.numint.eval_ao(mol, coords, deriv=0, shls_slice=None, non0tab=None, out=None, verbose=None)[source]

Evaluate AO function value on the given grids.

Args:

mol : an instance of Mole

coords
: 2D array, shape (N,3)
The coordinates of the grids.
Kwargs:
deriv
: int
AO derivative order. It affects the shape of the return array. If deriv=0, the returned AO values are stored in a (N,nao) array. Otherwise the AO values are stored in an array of shape (M,N,nao). Here N is the number of grids, nao is the number of AO functions, M is the size associated to the derivative deriv.
relativity
: bool
No effects.
shls_slice
: 2-element list
(shl_start, shl_end). If given, only part of AOs (shl_start <= shell_id < shl_end) are evaluated. By default, all shells defined in mol will be evaluated.
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling make_mask()
out
: ndarray
If provided, results are written into this array.
verbose
: int or object of Logger
No effects.
Returns:
2D array of shape (N,nao) for AO values if deriv = 0. Or 3D array of shape (:,N,nao) for AO values and AO derivatives if deriv > 0. In the 3D array, the first (N,nao) elements are the AO values, followed by (3,N,nao) for x,y,z compoents; Then 2nd derivatives (6,N,nao) for xx, xy, xz, yy, yz, zz; Then 3rd derivatives (10,N,nao) for xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz; ...

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz')
>>> coords = numpy.random.random((100,3))  # 100 random points
>>> ao_value = eval_ao(mol, coords)
>>> print(ao_value.shape)
(100, 24)
>>> ao_value = eval_ao(mol, coords, deriv=1, shls_slice=(1,4))
>>> print(ao_value.shape)
(4, 100, 7)
>>> ao_value = eval_ao(mol, coords, deriv=2, shls_slice=(1,4))
>>> print(ao_value.shape)
(10, 100, 7)
pyscf.dft.numint.eval_mat(mol, ao, weight, rho, vxc, non0tab=None, xctype='LDA', spin=0, verbose=None)[source]

Calculate XC potential matrix.

Args:

mol : an instance of Mole

ao
: ([4/10,] ngrids, nao) ndarray
2D array of shape (N,nao) for LDA, 3D array of shape (4,N,nao) for GGA or (10,N,nao) for meta-GGA. N is the number of grids, nao is the number of AO functions. If xctype is GGA, ao[0] is AO value and ao[1:3] are the real space gradients. If xctype is meta-GGA, ao[4:10] are second derivatives of ao values.
weight
: 1D array
Integral weights on grids.
rho
: ([4/6,] ngrids) ndarray

Shape of ((,N)) for electron density (and derivatives) if spin = 0; Shape of ((,N),(,N)) for alpha/beta electron density (and derivatives) if spin > 0; where N is number of grids. rho (,N) are ordered as (den,grad_x,grad_y,grad_z,laplacian,tau) where grad_x = d/dx den, laplacian = nabla^2 den, tau = 1/2(nabla f)^2 In spin unrestricted case, rho is ((den_u,grad_xu,grad_yu,grad_zu,laplacian_u,tau_u)

(den_d,grad_xd,grad_yd,grad_zd,laplacian_d,tau_d))
vxc
: ([4,] ngrids) ndarray
XC potential value on each grid = (vrho, vsigma, vlapl, vtau) vsigma is GGA potential value on each grid. If the kwarg spin != 0, a list [vsigma_uu,vsigma_ud] is required.
Kwargs:
xctype
: str
LDA/GGA/mGGA. It affects the shape of ao and rho
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling make_mask()
spin
: int
If not 0, the returned matrix is the Vxc matrix of alpha-spin. It is computed with the spin non-degenerated UKS formula.
Returns:
XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.
pyscf.dft.numint.eval_rho(mol, ao, dm, non0tab=None, xctype='LDA', hermi=0, verbose=None)[source]

Calculate the electron density for LDA functional, and the density derivatives for GGA functional.

Args:

mol : an instance of Mole

ao
: 2D array of shape (N,nao) for LDA, 3D array of shape (4,N,nao) for GGA
or (5,N,nao) for meta-GGA. N is the number of grids, nao is the number of AO functions. If xctype is GGA, ao[0] is AO value and ao[1:3] are the AO gradients. If xctype is meta-GGA, ao[4:10] are second derivatives of ao values.
dm
: 2D array
Density matrix
Kwargs:
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling make_mask()
xctype
: str
LDA/GGA/mGGA. It affects the shape of the return density.
hermi
: bool
dm is hermitian or not
verbose
: int or object of Logger
No effects.
Returns:
1D array of size N to store electron density if xctype = LDA; 2D array of (4,N) to store density and “density derivatives” for x,y,z components if xctype = GGA; (6,N) array for meta-GGA, where last two rows are nabla^2 rho and tau = 1/2(nabla f)^2

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz')
>>> coords = numpy.random.random((100,3))  # 100 random points
>>> ao_value = eval_ao(mol, coords, deriv=0)
>>> dm = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> dm = dm + dm.T
>>> rho, dx_rho, dy_rho, dz_rho = eval_rho(mol, ao, dm, xctype='LDA')
pyscf.dft.numint.eval_rho2(mol, ao, mo_coeff, mo_occ, non0tab=None, xctype='LDA', verbose=None)[source]

Calculate the electron density for LDA functional, and the density derivatives for GGA functional. This function has the same functionality as eval_rho() except that the density are evaluated based on orbital coefficients and orbital occupancy. It is more efficient than eval_rho() in most scenario.

Args:

mol : an instance of Mole

ao
: 2D array of shape (N,nao) for LDA, 3D array of shape (4,N,nao) for GGA
or (5,N,nao) for meta-GGA. N is the number of grids, nao is the number of AO functions. If xctype is GGA, ao[0] is AO value and ao[1:3] are the AO gradients. If xctype is meta-GGA, ao[4:10] are second derivatives of ao values.
dm
: 2D array
Density matrix
Kwargs:
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling make_mask()
xctype
: str
LDA/GGA/mGGA. It affects the shape of the return density.
verbose
: int or object of Logger
No effects.
Returns:
1D array of size N to store electron density if xctype = LDA; 2D array of (4,N) to store density and “density derivatives” for x,y,z components if xctype = GGA; (6,N) array for meta-GGA, where last two rows are nabla^2 rho and tau = 1/2(nabla f)^2
pyscf.dft.numint.get_rho(ni, mol, dm, grids, max_memory=2000)[source]

Density in real space

pyscf.dft.numint.nr_fxc(mol, grids, xc_code, dm0, dms, spin=0, relativity=0, hermi=0, rho0=None, vxc=None, fxc=None, max_memory=2000, verbose=None)[source]

Contract XC kernel matrix with given density matrices

... math:

a_{pq} = f_{pq,rs} * x_{rs}
pyscf.dft.numint.nr_rks(ni, mol, grids, xc_code, dms, relativity=0, hermi=0, max_memory=2000, verbose=None)[source]

Calculate RKS XC functional and potential matrix on given meshgrids for a set of density matrices

Args:

ni : an instance of NumInt

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: 2D array a list of 2D arrays
Density matrix or multiple density matrices
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
Returns:
nelec, excsum, vmat. nelec is the number of electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.

Examples:

>>> from pyscf import gto, dft
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.coords = numpy.random.random((100,3))  # 100 random points
>>> grids.weights = numpy.random.random(100)
>>> nao = mol.nao_nr()
>>> dm = numpy.random.random((nao,nao))
>>> ni = dft.numint.NumInt()
>>> nelec, exc, vxc = ni.nr_rks(mol, grids, 'lda,vwn', dm)
pyscf.dft.numint.nr_rks_fxc(ni, mol, grids, xc_code, dm0, dms, relativity=0, hermi=0, rho0=None, vxc=None, fxc=None, max_memory=2000, verbose=None)[source]

Contract RKS XC (singlet hessian) kernel matrix with given density matrices

Args:

ni : an instance of NumInt

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: 2D array a list of 2D arrays
Density matrix or multiple density matrices
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
rho0
: float array
Zero-order density (and density derivative for GGA). Giving kwargs rho0, vxc and fxc to improve better performance.
vxc
: float array
First order XC derivatives
fxc
: float array
Second order XC derivatives
Returns:
nelec, excsum, vmat. nelec is the number of electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.

Examples:

pyscf.dft.numint.nr_rks_fxc_st(ni, mol, grids, xc_code, dm0, dms_alpha, relativity=0, singlet=True, rho0=None, vxc=None, fxc=None, max_memory=2000, verbose=None)[source]

Associated to singlet or triplet Hessian Note the difference to nr_rks_fxc, dms_alpha is the response density matrices of alpha spin, alpha+/-beta DM is applied due to singlet/triplet coupling

Ref. CPL, 256, 454

pyscf.dft.numint.nr_rks_vxc(ni, mol, grids, xc_code, dms, relativity=0, hermi=0, max_memory=2000, verbose=None)

Calculate RKS XC functional and potential matrix on given meshgrids for a set of density matrices

Args:

ni : an instance of NumInt

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: 2D array a list of 2D arrays
Density matrix or multiple density matrices
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
Returns:
nelec, excsum, vmat. nelec is the number of electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.

Examples:

>>> from pyscf import gto, dft
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.coords = numpy.random.random((100,3))  # 100 random points
>>> grids.weights = numpy.random.random(100)
>>> nao = mol.nao_nr()
>>> dm = numpy.random.random((nao,nao))
>>> ni = dft.numint.NumInt()
>>> nelec, exc, vxc = ni.nr_rks(mol, grids, 'lda,vwn', dm)
pyscf.dft.numint.nr_uks(ni, mol, grids, xc_code, dms, relativity=0, hermi=0, max_memory=2000, verbose=None)[source]

Calculate UKS XC functional and potential matrix on given meshgrids for a set of density matrices

Args:

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: a list of 2D arrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
Returns:
nelec, excsum, vmat. nelec is the number of (alpha,beta) electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix for (alpha,beta) spin.

Examples:

>>> from pyscf import gto, dft
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.coords = numpy.random.random((100,3))  # 100 random points
>>> grids.weights = numpy.random.random(100)
>>> nao = mol.nao_nr()
>>> dm = numpy.random.random((2,nao,nao))
>>> ni = dft.numint.NumInt()
>>> nelec, exc, vxc = ni.nr_uks(mol, grids, 'lda,vwn', dm)
pyscf.dft.numint.nr_uks_fxc(ni, mol, grids, xc_code, dm0, dms, relativity=0, hermi=0, rho0=None, vxc=None, fxc=None, max_memory=2000, verbose=None)[source]

Contract UKS XC kernel matrix with given density matrices

Args:

ni : an instance of NumInt

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: 2D array a list of 2D arrays
Density matrix or multiple density matrices
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
rho0
: float array
Zero-order density (and density derivative for GGA). Giving kwargs rho0, vxc and fxc to improve better performance.
vxc
: float array
First order XC derivatives
fxc
: float array
Second order XC derivatives
Returns:
nelec, excsum, vmat. nelec is the number of electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.

Examples:

pyscf.dft.numint.nr_uks_vxc(ni, mol, grids, xc_code, dms, relativity=0, hermi=0, max_memory=2000, verbose=None)

Calculate UKS XC functional and potential matrix on given meshgrids for a set of density matrices

Args:

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: a list of 2D arrays
A list of density matrices, stored as (alpha,alpha,...,beta,beta,...)
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
Returns:
nelec, excsum, vmat. nelec is the number of (alpha,beta) electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix for (alpha,beta) spin.

Examples:

>>> from pyscf import gto, dft
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.coords = numpy.random.random((100,3))  # 100 random points
>>> grids.weights = numpy.random.random(100)
>>> nao = mol.nao_nr()
>>> dm = numpy.random.random((2,nao,nao))
>>> ni = dft.numint.NumInt()
>>> nelec, exc, vxc = ni.nr_uks(mol, grids, 'lda,vwn', dm)
pyscf.dft.numint.nr_vxc(mol, grids, xc_code, dms, spin=0, relativity=0, hermi=0, max_memory=2000, verbose=None)[source]

Evaluate RKS/UKS XC functional and potential matrix on given meshgrids for a set of density matrices. See nr_rks() and nr_uks() for more details.

Args:

mol : an instance of Mole

grids
: an instance of Grids
grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
xc_code
: str
XC functional description. See parse_xc() of pyscf/dft/libxc.py for more details.
dms
: 2D array or a list of 2D arrays
Density matrix or multiple density matrices
Kwargs:
hermi
: int
Input density matrices symmetric or not
max_memory
: int or float
The maximum size of cache to use (in MB).
Returns:
nelec, excsum, vmat. nelec is the number of electrons generated by numerical integration. excsum is the XC functional value. vmat is the XC potential matrix in 2D array of shape (nao,nao) where nao is the number of AO functions.

Examples:

>>> from pyscf import gto, dft
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> grids = dft.gen_grid.Grids(mol)
>>> grids.coords = numpy.random.random((100,3))  # 100 random points
>>> grids.weights = numpy.random.random(100)
>>> nao = mol.nao_nr()
>>> dm = numpy.random.random((2,nao,nao))
>>> nelec, exc, vxc = dft.numint.nr_vxc(mol, grids, 'lda,vwn', dm, spin=1)

XC functional, the interface to libxc (http://www.tddft.org/programs/octopus/wiki/index.php/Libxc)

pyscf.dft.libxc.define_xc(ni, description, xctype='LDA', hyb=0, rsh=(0, 0, 0))[source]

Define XC functional. See also eval_xc() for the rules of input description.

Args:

ni : an instance of NumInt

description
: str
A string to describe the linear combination of different XC functionals. The X and C functional are separated by comma like ‘.8*LDA+.2*B86,VWN’. If “HF” was appeared in the string, it stands for the exact exchange.
Kwargs:
xctype
: str
‘LDA’ or ‘GGA’ or ‘MGGA’
hyb
: float
hybrid functional coefficient
rsh
: float
coefficients for range-separated hybrid functional

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz')
>>> mf = dft.RKS(mol)
>>> define_xc_(mf._numint, '.2*HF + .08*LDA + .72*B88, .81*LYP + .19*VWN')
>>> mf.kernel()
-76.3783361189611
>>> define_xc_(mf._numint, 'LDA*.08 + .72*B88 + .2*HF, .81*LYP + .19*VWN')
>>> mf.kernel()
-76.3783361189611
>>> def eval_xc(xc_code, rho, *args, **kwargs):
...     exc = 0.01 * rho**2
...     vrho = 0.01 * 2 * rho
...     vxc = (vrho, None, None, None)
...     fxc = None  # 2nd order functional derivative
...     kxc = None  # 3rd order functional derivative
...     return exc, vxc, fxc, kxc
>>> define_xc_(mf._numint, eval_xc, xctype='LDA')
>>> mf.kernel()
48.8525211046668
pyscf.dft.libxc.define_xc_(ni, description, xctype='LDA', hyb=0, rsh=(0, 0, 0))[source]

Define XC functional. See also eval_xc() for the rules of input description.

Args:

ni : an instance of NumInt

description
: str
A string to describe the linear combination of different XC functionals. The X and C functional are separated by comma like ‘.8*LDA+.2*B86,VWN’. If “HF” was appeared in the string, it stands for the exact exchange.
Kwargs:
xctype
: str
‘LDA’ or ‘GGA’ or ‘MGGA’
hyb
: float
hybrid functional coefficient
rsh
: float
coefficients for range-separated hybrid functional

Examples:

>>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz')
>>> mf = dft.RKS(mol)
>>> define_xc_(mf._numint, '.2*HF + .08*LDA + .72*B88, .81*LYP + .19*VWN')
>>> mf.kernel()
-76.3783361189611
>>> define_xc_(mf._numint, 'LDA*.08 + .72*B88 + .2*HF, .81*LYP + .19*VWN')
>>> mf.kernel()
-76.3783361189611
>>> def eval_xc(xc_code, rho, *args, **kwargs):
...     exc = 0.01 * rho**2
...     vrho = 0.01 * 2 * rho
...     vxc = (vrho, None, None, None)
...     fxc = None  # 2nd order functional derivative
...     kxc = None  # 3rd order functional derivative
...     return exc, vxc, fxc, kxc
>>> define_xc_(mf._numint, eval_xc, xctype='LDA')
>>> mf.kernel()
48.8525211046668
pyscf.dft.libxc.eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, verbose=None)[source]

Interface to call libxc library to evaluate XC functional, potential and functional derivatives.

  • The given functional xc_code must be a one-line string.
  • The functional xc_code is case-insensitive.
  • The functional xc_code string has two parts, separated by ”,”. The first part describes the exchange functional, the second is the correlation functional.
    • If ”,” not appeared in string, the entire string is considered as X functional.
    • To neglect X functional (just apply C functional), leave blank in the first part, eg description=’,vwn’ for pure VWN functional
  • The functional name can be placed in arbitrary order. Two name needs to be separated by operators “+” or “-”. Blank spaces are ignored. NOTE the parser only reads operators “+” “-” “*”. / is not in support.
  • A functional name is associated with one factor. If the factor is not given, it is assumed equaling 1.
  • String “HF” stands for exact exchange (HF K matrix). It is allowed to put in C functional part.
  • Be careful with the libxc convention on GGA functional, in which the LDA contribution is included.
Args:
xc_code
: str
A string to describe the linear combination of different XC functionals. The X and C functional are separated by comma like ‘.8*LDA+.2*B86,VWN’. If “HF” (exact exchange) is appeared in the string, the HF part will be skipped. If an empty string “” is given, the returns exc, vxc,... will be vectors of zeros.
rho
: ndarray

Shape of ((,N)) for electron density (and derivatives) if spin = 0; Shape of ((,N),(,N)) for alpha/beta electron density (and derivatives) if spin > 0; where N is number of grids. rho (,N) are ordered as (den,grad_x,grad_y,grad_z,laplacian,tau) where grad_x = d/dx den, laplacian = nabla^2 den, tau = 1/2(nabla f)^2 In spin unrestricted case, rho is ((den_u,grad_xu,grad_yu,grad_zu,laplacian_u,tau_u)

(den_d,grad_xd,grad_yd,grad_zd,laplacian_d,tau_d))
Kwargs:
spin
: int
spin polarized if spin > 0
relativity
: int
No effects.
verbose
: int or object of Logger
No effects.
Returns:

ex, vxc, fxc, kxc

where

  • vxc = (vrho, vsigma, vlapl, vtau) for restricted case
  • vxc for unrestricted case | vrho[:,2] = (u, d) | vsigma[:,3] = (uu, ud, dd) | vlapl[:,2] = (u, d) | vtau[:,2] = (u, d)
  • fxc for restricted case: (v2rho2, v2rhosigma, v2sigma2, v2lapl2, vtau2, v2rholapl, v2rhotau, v2lapltau, v2sigmalapl, v2sigmatau)
  • fxc for unrestricted case: | v2rho2[:,3] = (u_u, u_d, d_d) | v2rhosigma[:,6] = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd) | v2sigma2[:,6] = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd) | v2lapl2[:,3] | vtau2[:,3] | v2rholapl[:,4] | v2rhotau[:,4] | v2lapltau[:,4] | v2sigmalapl[:,6] | v2sigmatau[:,6]
  • kxc for restricted case: (v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)
  • kxc for unrestricted case: | v3rho3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d) | v3rho2sigma[:,9] = (u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd) | v3rhosigma2[:,12] = (u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd) | v3sigma3[:,10] = (uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd)

see also libxc_itrf.c

pyscf.dft.libxc.hybrid_coeff(xc_code, spin=0)[source]

Support recursively defining hybrid functional

pyscf.dft.libxc.nlc_coeff(xc_code)[source]

Get NLC coefficients

pyscf.dft.libxc.parse_xc(description)[source]

Rules to input functional description:

  • The given functional description must be a one-line string.
  • The functional description is case-insensitive.
  • The functional description string has two parts, separated by ”,”. The first part describes the exchange functional, the second is the correlation functional.
    • If ”,” was not in string, the entire string is considered as a compound XC functional (including both X and C functionals, such as b3lyp).
    • To input only X functional (without C functional), leave the second part blank. E.g. description=’slater,’ means pure LDA functional.
    • To neglect X functional (just apply C functional), leave the first part blank. E.g. description=’,vwn’ means pure VWN functional.
    • If compound XC functional is specified, no matter whehter it is in the X part (the string in front of comma) or the C part (the string behind comma), both X and C functionals of the compound XC functional will be used.
  • The functional name can be placed in arbitrary order. Two name needs to be separated by operators “+” or “-”. Blank spaces are ignored. NOTE the parser only reads operators “+” “-” “*”. / is not in support.
  • A functional name can have at most one factor. If the factor is not given, it is set to 1. Compound functional can be scaled as a unit. For example ‘0.5*b3lyp’ is equivalent to ‘HF*0.1 + .04*LDA + .36*B88, .405*LYP + .095*VWN’
  • String “HF” stands for exact exchange (HF K matrix). Putting “HF” in correlation functional part is the same to putting “HF” in exchange part.
  • String “RSH” means range-separated operator. Its format is RSH(alpha; beta; omega). Another way to input RSH is to use keywords SR_HF and LR_HF: “SR_HF(0.1) * alpha_plus_beta” and “LR_HF(0.1) * alpha” where the number in parenthesis is the value of omega.
  • Be careful with the libxc convention on GGA functional, in which the LDA contribution has been included.
Args:
xc_code
: str
A string to describe the linear combination of different XC functionals. The X and C functional are separated by comma like ‘.8*LDA+.2*B86,VWN’. If “HF” was appeared in the string, it stands for the exact exchange.
rho
: ndarray

Shape of ((,N)) for electron density (and derivatives) if spin = 0; Shape of ((,N),(,N)) for alpha/beta electron density (and derivatives) if spin > 0; where N is number of grids. rho (,N) are ordered as (den,grad_x,grad_y,grad_z,laplacian,tau) where grad_x = d/dx den, laplacian = nabla^2 den, tau = 1/2(nabla f)^2 In spin unrestricted case, rho is ((den_u,grad_xu,grad_yu,grad_zu,laplacian_u,tau_u)

(den_d,grad_xd,grad_yd,grad_zd,laplacian_d,tau_d))
Kwargs:
spin
: int
spin polarized if spin > 0
relativity
: int
No effects.
verbose
: int or object of Logger
No effects.
Returns:

ex, vxc, fxc, kxc

where

  • vxc = (vrho, vsigma, vlapl, vtau) for restricted case
  • vxc for unrestricted case | vrho[:,2] = (u, d) | vsigma[:,3] = (uu, ud, dd) | vlapl[:,2] = (u, d) | vtau[:,2] = (u, d)
  • fxc for restricted case: (v2rho2, v2rhosigma, v2sigma2, v2lapl2, vtau2, v2rholapl, v2rhotau, v2lapltau, v2sigmalapl, v2sigmatau)
  • fxc for unrestricted case: | v2rho2[:,3] = (u_u, u_d, d_d) | v2rhosigma[:,6] = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd) | v2sigma2[:,6] = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd) | v2lapl2[:,3] | vtau2[:,3] | v2rholapl[:,4] | v2rhotau[:,4] | v2lapltau[:,4] | v2sigmalapl[:,6] | v2sigmatau[:,6]
  • kxc for restricted case: v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, v3rho2tau, v3rhosigmatau, v3rhotau2, v3sigma2tau, v3sigmatau2, v3tau3
  • kxc for unrestricted case: | v3rho3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d) | v3rho2sigma[:,9] = (u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd) | v3rhosigma2[:,12] = (u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd) | v3sigma3[:,10] = (uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd) | v3rho2tau | v3rhosigmatau | v3rhotau2 | v3sigma2tau | v3sigmatau2 | v3tau3

see also libxc_itrf.c

pyscf.dft.libxc.parse_xc_name(xc_name='LDA, VWN')[source]

Convert the XC functional name to libxc library internal ID.

pyscf.dft.libxc.rsh_coeff(xc_code)[source]

Range-separated parameter and HF exchange components: omega, alpha, beta

Exc_RSH = c_SR * SR_HFX + c_LR * LR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec
= alpha * HFX + beta * SR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec = alpha * LR_HFX + hyb * SR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec

SR_HFX = < pi | e^{-omega r_{12}}/r_{12} | iq > LR_HFX = < pi | (1-e^{-omega r_{12}})/r_{12} | iq > alpha = c_LR beta = c_SR - c_LR = hyb - alpha