pyscf.mcpdft package#

Submodules#

pyscf.mcpdft.chkfile module#

pyscf.mcpdft.chkfile.dump_lpdft(mc, chkfile=None, key='pdft', e_tot=None, e_states=None, e_mcscf=None, ci=None, mcscf_key='mcscf')[source]#: Save L-PDFT calculation results in chkfile

pyscf.mcpdft.chkfile.dump_mcpdft(mc, chkfile=None, key='pdft', e_tot=None, e_ot=None, e_states=None, e_mcscf=None, mcscf_key='mcscf')[source]#: Save MC-PDFT calculation results in chkfile

pyscf.mcpdft.chkfile.load_pdft(chkfile)[source]#

pyscf.mcpdft.cmspdft module#

pyscf.mcpdft.cmspdft.coulomb_tensor(mc, mo_coeff=None, ci=None, h2eff=None, eris=None)[source]#

Compute w_IJKL = (tu|vx) D^IJ_tu D^KL_vx

Args:

mc : mcscf method instance

Kwargs:

mo_coeffndarray of shape (nao,nmo): Contains molecular orbital coefficients
cilist of ndarrays of shape (ndeta,ndetb): Contains CI vectors
h2effndarray of shape [ncas,]*4: Contains active-space ERIs
erismc_ao2mo.ERI object: Contains active-space ERIs. Ignored if h2eff is passed; if h2eff is not passed then it is constructed from eris.ppaa

Returns:

w : ndarray of shape [nroots,]*4

pyscf.mcpdft.cmspdft.e_coul(mc, mo_coeff=None, ci=None, h2eff=None, eris=None)[source]#

Compute the sum of active-space Coulomb energies (the diabatizer function for CMS-PDFT) and its first and second derivatives

Args:

mc : mcscf method instance

Kwargs:

mo_coeffndarray of shape (nao,nmo): Contains molecular orbital coefficients
cilist of ndarrays of shape (ndeta,ndetb): Contains CI vectors
h2effndarray of shape [ncas,]*4: Contains active-space ERIs
erismc_ao2mo.ERI object: Contains active-space ERIs. Ignored if h2eff is passed; if h2eff is not passed then it is constructed from eris.ppaa

Returns:

Qaafloat: sum of Coulomb energies
dQaandarray of shape npair = nroots*(nroots-1)/2: first derivatives of J wrt interstate rotation
d2Qaandarray of shape (npair,npair): Hessian of J wrt interstate rotation
Qaa_updatecallable: Takes a unitary matrix of shape (nroots, nroots) and returns Qaa, dQaa, and d2Qaa as above using the stored Coulomb tensor intermediate from this function.

pyscf.mcpdft.cmspdft.e_coul_o0(mc, ci)[source]#

pyscf.mcpdft.lpdft module#

pyscf.mcpdft.lpdft.get_lpdft_hconst(mc, E_ot, casdm1s_0, casdm2_0, hyb=1.0, ncas=None, ncore=None, veff1=None, veff2=None, mo_coeff=None)[source]#

Compute h_const for the L-PDFT Hamiltonian

Args:

mc : instance of class _PDFT

E_otfloat: On-top energy
casdm1s_0ndarray of shape (2, ncas, ncas): Spin-separated 1-RDM in the active space generated from expansion density.
casdm2_0ndarray of shape (ncas, ncas, ncas, ncas): Spin-summed 2-RDM in the active space generated from expansion density.

Kwargs:

hybfloat: Hybridization constant (lambda term)
ncasfloat: Number of active space MOs
ncore: float: Number of core MOs
veff1ndarray of shape (nao, nao): 1-body effective potential in the AO basis computed using the zeroth-order densities.
veff2pyscf.mcscf.mc_ao2mo._ERIS instance: Relevant 2-body effective potential in the MO basis.

Returns:

Constant term h_const for the expansion term.

pyscf.mcpdft.lpdft.get_transformed_h2eff_for_cas(mc, ncore=None, ncas=None)[source]#

Compute the CAS two-particle linear PDFT Hamiltonian

Args:

ncoreint: Number of core MOs
ncasint: Number of active space MOs

Returns:

ndarray of shape (ncas,ncas,ncas,ncas) which contain v_vwxy

pyscf.mcpdft.lpdft.kernel(mc, mo_coeff=None, ci0=None, ot=None, dump_chk=True, **kwargs)[source]#

pyscf.mcpdft.lpdft.linear_multi_state(mc, weights=(0.5, 0.5), **kwargs)[source]#

Build linearized multi-state MC-PDFT method object

Args:: mc : instance of class _PDFT
Kwargs:: weights : sequence of floats
Returns:: si : instance of class _LPDFT

pyscf.mcpdft.lpdft.linear_multi_state_mix(mc, fcisolvers, weights=(0.5, 0.5), **kwargs)[source]#

Build SA Mix linearized multi-state MC-PDFT method object

Args:

mc : instance of class _PDFT

fcisolvers : fcisolvers to construct StateAverageMixSolver with

Kwargs:

weights : sequence of floats

Returns:

si : instance of class _LPDFT

pyscf.mcpdft.lpdft.make_lpdft_ham_(mc, mo_coeff=None, ci=None, ot=None)[source]#

Compute the L-PDFT Hamiltonian

Args:

mo_coeffndarray of shape (nao, nmo): A full set of molecular orbital coefficients. Taken from self if not provided.
cilist of ndarrays of length nroots: CI vectors should be from a converged CASSCF/CASCI calculation

ot : an instance of on-top functional class - see otfnal.py

Returns:

lpdft_hamndarray of shape (nroots, nroots) or (nirreps, nroots, nroots): Linear approximation to the MC-PDFT energy expressed as a hamiltonian in the basis provided by the CI vectors. If StateAverageMix, then returns the block diagonal of the lpdft hamiltonian for each irrep.

pyscf.mcpdft.lpdft.transformed_h1e_for_cas(mc, E_ot, casdm1s_0, casdm2_0, hyb=1.0, mo_coeff=None, ncas=None, ncore=None)[source]#

Compute the CAS one-particle L-PDFT Hamiltonian

Args:

mc : instance of a _PDFT object

E_otfloat: On-top energy
casdm1s_0ndarray of shape (2,ncas,ncas): Spin-separated 1-RDM in the active space generated from expansion density
casdm2_0ndarray of shape (ncas,ncas,ncas,ncas): Spin-summed 2-RDM in the active space generated from expansion density
hybfloat: Hybridization constant (lambda term)
mo_coeffndarray of shape (nao,nmo): A full set of molecular orbital coefficients. Taken from self if not provided.
ncasint: Number of active space molecular orbitals
ncoreint: Number of core molecular orbitals

Returns:

A tuple, the first is the effective one-electron linear PDFT Hamiltonian defined in CAS space, the second is the modified core energy.

pyscf.mcpdft.lpdft.weighted_average_densities(mc, ci=None, weights=None)[source]#

Compute the weighted average 1- and 2-electron CAS densities. 1-electron CAS is returned as spin-separated.

Args:

mc : instance of class _PDFT

cilist of ndarrays of length nroots: CI vectors should be from a converged CASSCF/CASCI calculation
weightsndarray of length nroots: Weight for each state. If none, uses weights from SA-CASSCF calculation

Returns:

A tuple, the first is casdm1s and the second is casdm2 where they are weighted averages where the weights are given.

pyscf.mcpdft.mcpdft module#

pyscf.mcpdft.mcpdft.energy_dft(mc, mo_coeff=None, ci=None, ot=None, state=0, casdm1s=None, casdm2=None, max_memory=None, hermi=1)[source]#: Wrap to ot.energy_ot for subclassing.

pyscf.mcpdft.mcpdft.energy_elec(mc, *args, **kwargs)[source]#

pyscf.mcpdft.mcpdft.energy_mcwfn(mc, mo_coeff=None, ci=None, ot=None, state=0, casdm1s=None, casdm2=None, verbose=None)[source]#

Compute the parts of the MC-PDFT energy arising from the wave function

Args:

mcan instance of CASSCF or CASCI class: Note: this function does not currently run the CASSCF or CASCI calculation itself prior to calculating the MC-PDFT energy. Call mc.kernel () before passing to this function!

Kwargs:

mo_coeffndarray of shape (nao, nmo): contains molecular orbital coefficients
cilist or ndarray: contains ci vectors

ot : an instance of on-top functional class - see otfnal.py state : int

If mc describes a state-averaged calculation, select the state (0-indexed).

casdm1sndarray or compatible of shape (2,ncas,ncas): Contains spin-separated active-space 1RDM
casdm2ndarray or compatible of shape [ncas,]*4: Contains spin-summed active-space 2RDM

Returns:

e_mcwfnfloat: Energy from the multiconfigurational wave function: nuclear repulsion + 1e + coulomb

pyscf.mcpdft.mcpdft.energy_tot(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None)[source]#

Calculate MC-PDFT total energy

Args:

mcan instance of CASSCF or CASCI class: Note: this function does not currently run the CASSCF or CASCI calculation itself prior to calculating the MC-PDFT energy. Call mc.kernel () before passing to this function!

Kwargs:

mo_coeffndarray of shape (nao, nmo): Molecular orbital coefficients
cindarray or list of length (nroots): CI vector or vectors.

ot : an instance of on-top functional class - see otfnal.py state : int

If mc describes a state-averaged calculation, select the state (0-indexed).

verboseint: Verbosity of logger output; defaults to mc.verbose

Returns:

e_totfloat: Total MC-PDFT energy including nuclear repulsion energy
E_otfloat: On-top (cf. exchange-correlation) energy

pyscf.mcpdft.mcpdft.get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, grids_level=None, grids_attr=None, split_x_c=None, verbose=None)[source]#

Compute a decomposition of the MC-PDFT energy into nuclear potential (E0), one-electron (E1), Coulomb (E2c), exchange (EOTx), correlation (EOTc) terms, and additionally the nonclassical part (E2nc) of the MC-SCF energy:

E(MC-SCF) = E0 + E1 + E2c + Enc E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc

For hybrid functionals,

E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc + Enc

Where the Enc and EOTx/c terms are premultiplied by the hybrid factor. If mc.fcisolver.nroots > 1, lists are returned for everything except the nuclear potential energy.

Args:

mc : an instance of CASSCF or CASCI class

Kwargs:

mo_coeffndarray: Contains MO coefficients
cindarray or list of length nroots: Contains CI vectors

ot : an instance of (translated) on-top density fnal class otxc : string

identity of translated functional; overrides ot

grids_levelinteger: level preset for DFT quadrature grids
grids_attrdictionary: general attributes for DFT quadrature grids
split_x_clogical: whether to split the exchange and correlation parts of the ot functional into two separate contributions

Returns:

e_nucfloat: E0 = sum_A>B ZA*ZB/rAB
e_1efloat or list of length nroots: E1 = <T+sum_A ZA/rA>
e_coulfloat or list of length nroots: E2c = 1/2 int rho(1)rho(2)/r12 d1d2
e_otxcfloat or list of length nroots: EOTxc = translated functional energy if split_x_c == True, this is instead EOTx = exchange part of translated functional energy
e_otcfloat or list of length nroots: only returned if split_x_c == True EOTc = correlation part of translated functional
e_ncwfnfloat or list of length nroots: E2ncc = <H> - E0 - E1 - E2c If hybrid functional, this term is weighted appropriately. For pure functionals, it is the full NC component

pyscf.mcpdft.mcpdft.get_mcpdft_child_class(mc, ot, **kwargs)[source]#

pyscf.mcpdft.mcpdft.kernel(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None)#

Calculate MC-PDFT total energy

Args:

mcan instance of CASSCF or CASCI class: Note: this function does not currently run the CASSCF or CASCI calculation itself prior to calculating the MC-PDFT energy. Call mc.kernel () before passing to this function!

Kwargs:

mo_coeffndarray of shape (nao, nmo): Molecular orbital coefficients
cindarray or list of length (nroots): CI vector or vectors.

ot : an instance of on-top functional class - see otfnal.py state : int

If mc describes a state-averaged calculation, select the state (0-indexed).

verboseint: Verbosity of logger output; defaults to mc.verbose

Returns:

e_totfloat: Total MC-PDFT energy including nuclear repulsion energy
E_otfloat: On-top (cf. exchange-correlation) energy

pyscf.mcpdft.mspdft module#

pyscf.mcpdft.mspdft.get_diabfns(obj)[source]#

Interpret the name of the MS-PDFT method as a pair of functions which optimize the intermediate states and calculate the power series in the corresponding objective function to second order.

Args:

objstring: Specify particular MS-PDFT method. Currently, only “CMS” is supported. Not case-sensitive.

Returns:

diabatizercallable: Takes model-space CI vectors in a trial intermediate-state basis and returns the value and first and second derivatives of the objective function specified by obj
diabatizecallable: Takes model-space CI vectors and returns CI vectors in the optimized intermediate-state basis

pyscf.mcpdft.mspdft.make_heff_mcscf(mc, mo_coeff=None, ci=None)[source]#

Build Hamiltonian matrix in basis of ci vector

Args:

mc : an instance of MCPDFT class

Kwargs:

mo_coeffndarray of shape (nao, nmo): MO coefficients
cindarray or list of len (nroots): CI vectors describing the model space, presumed to be in the optimized intermediate-state basis

Returns:

heff_mcscfndarray of shape (nroots, nroots): Effective MC-SCF hamiltonian matrix in the basis of the provided CI vectors

pyscf.mcpdft.mspdft.multi_state(mc, weights=(0.5, 0.5), diabatization='CMS', **kwargs)[source]#

Build multi-state MC-PDFT method object

Args:: mc : instance of class _PDFT
Kwargs:: weights : sequence of floats diabatization : objective-function type

Currently supports only ‘cms’
Returns:: si : instance of class _MSPDFT

pyscf.mcpdft.mspdft.si_newton(mc, ci=None, objfn=None, max_cyc=None, conv_tol=None, sing_tol=None, nudge_tol=None)[source]#

Optimize the intermediate states describing the model space of an MS-PDFT calculation by maximizing the provided objective function using a gradient-ascent algorithm

Args:

mc : an instance of MSPDFT class

Kwargs:

cindarray or list of len (nroots): CI vectors spanning the model space
objfncallable: Takes CI vectors as a kwarg and returns the value, gradient, and Hessian of a chosen objective function wrt rotation between pairs of CI vectors
max_cycinteger: Maximum number of cycles of the gradient-ascent algorithm
conv_tolfloat: Maximum value of both gradient and step vectors at convergence
sing_tolfloat: Tolerance for determining when normal coordinate belongs to the null space (df = d2f = 0) or when the Hessian is singular (df != 0, d2f = 0).
nudge_tolfloat: Minimum step size along a normal coordinate when the surface is locally concave.

Returns:

convlogical: True if the optimization is converged
cilist of len (nroots): Optimized CI vectors describing intermediate states

pyscf.mcpdft.otfnal module#

class pyscf.mcpdft.otfnal.colle_salvetti_corr(mol, **kwargs)[source]#

Bases: otfnal

get_E_ot(rho, Pi, weights)[source]#: Colle & Salvetti, Theor. Chim. Acta 37, 329 (1975) see also Lee, Yang, Parr, Phys. Rev. B 37, 785 (1988) [Eq. (3)]

pyscf.mcpdft.otfnal.energy_ot(ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1)[source]#

Compute the on-top energy - the last term in

E_MCPDFT = h_pq l_pq + 1/2 v_pqrs l_pq l_rs + E_ot[rho,Pi]

Args:

ot : an instance of otfnal class casdm1s : ndarray of shape (2, ncas, ncas)

Contains spin-separated one-body density matrices in an active-orbital basis

casdm2ndarray of shape (ncas, ncas, ncas, ncas): Contains spin-summed two-body density matrix in an active- orbital basis
mo_coeffndarray of shape (nao, nmo): Contains molecular orbital coefficients for active-space orbitals. Columns ncore through ncore+ncas give the basis in which casdm1s and casdm2 are expressed.
ncoreinteger: Number of doubly occupied inactive “core” orbitals not explicitly included in casdm1s and casdm2

Kwargs:

max_memoryint or float: maximum cache size in MB default is 2000
hermiint: 1 if 1rdms are assumed hermitian, 0 otherwise

Returnsfloat

The MC-PDFT on-top (nonclassical) energy

pyscf.mcpdft.otfnal.ft_eval_xc(ni, xc_code, rho, spin=0, relativity=0, deriv=1, verbose=None)[source]#

For ‘fully translated’ functionals, otxc string = ‘ft’+xc string 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 part sets the correlation functional.

If “,” not appeared in string, the entire string is treated as the name of a compound functional (containing both the exchange and the correlation functional) which was declared in the functional aliases list. The full list of functional aliases can be obtained by calling the function pyscf.dft.xcfun.XC_ALIAS.keys() .

If the string was not found in the aliased functional list, it is treated as X functional.

To input only X functional (without C functional), leave the second part blank. E.g. description=’slater,’ means a functional with LDA contribution only.

To neglect the contribution of X functional (just apply C functional), leave blank in the first part, e.g. description=’,vwn’ means a functional with VWN only.

If compound XC functional is specified, no matter whether 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 names need to be separated by operators “+” or “-”. Blank spaces are ignored. NOTE the parser only reads operators “+” “-” “*”. / is not supported.

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). “HF” can be put in the correlation functional part (after comma). Putting “HF” in the correlation part is the same to putting “HF” in the exchange part.

String “RSH” means range-separated operator. Its format is RSH(omega, alpha, beta). 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 of GGA functional, in which the LDA contribution is included.

Args:

xc_codestr
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.

rhondarray
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:

spinint
spin polarized if spin > 0

relativityint
No effects.

verboseint 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, v2tau2, 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] | v2tau2[:,3] = (u_u, u_d, d_d) | v2rholapl[:,4] | v2rhotau[:,4] = (u_u, u_d, d_u, d_d) | v2lapltau[:,4] | v2sigmalapl[:,6] | v2sigmatau[:,6] = (uu_u, uu_d, ud_u, ud_d, dd_u, dd_d)

kxc for restricted case: (v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,

v3rho2lapl, v3rho2tau, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, 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) | v3rho2lapl[:,6] | v3rho2tau[:,6] = (u_u_u, u_u_d, u_d_u, u_d_d, d_d_u, d_d_d) | v3rhosigmalapl[:,12] | v3rhosigmatau[:,12] = (u_uu_u, u_uu_d, u_ud_u, u_ud_d, u_dd_u, u_dd_d,

d_uu_u, d_uu_d, d_ud_u, d_ud_d, d_dd_u, d_dd_d)

v3rholapl2[:,6]

v3rholapltau[:,8]

v3rhotau2[:,6] = (u_u_u, u_u_d, u_d_d, d_u_u, d_u_d, d_d_d)

v3sigma2lapl[:,12]

v3sigma2tau[:,12] = (uu_uu_u, uu_uu_d, uu_ud_u, uu_ud_d, uu_dd_u, uu_dd_d, ud_ud_u, ud_ud_d, ud_dd_u, ud_dd_d, dd_dd_u, dd_dd_d)

v3sigmalapl2[:,9]

v3sigmalapltau[:,12]

v3sigmatau2[:,9] = (uu_u_u, uu_u_d, uu_d_d, ud_u_u, ud_u_d, ud_d_d, dd_u_u, dd_u_d, dd_d_d)

v3lapl3[:,4]

v3lapl2tau[:,6]

v3lapltau2[:,6]

v3tau3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d)

see also libxc_itrf.c

pyscf.mcpdft.otfnal.ft_hybrid_coeff(ni, xc_code, spin=0)[source]#: For ‘fully translated’ functionals, otxc string = ‘ft’+xc string None

pyscf.mcpdft.otfnal.ft_nlc_coeff(ni, xc_code)[source]#: For ‘fully translated’ functionals, otxc string = ‘ft’+xc string None

pyscf.mcpdft.otfnal.ft_rsh_and_hybrid_coeff(ni, xc_code, spin=0)[source]#

For ‘fully translated’ functionals, otxc string = ‘ft’+xc string 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 | (1-erf(-omega r_{12}))/r_{12} | iq > LR_HFX = < pi | erf(-omega r_{12})/r_{12} | iq > alpha = c_LR beta = c_SR - c_LR

pyscf.mcpdft.otfnal.ft_rsh_coeff(ni, xc_code)[source]#: For ‘fully translated’ functionals, otxc string = ‘ft’+xc string None

pyscf.mcpdft.otfnal.ft_xc_type(ni, xc_code)[source]#: For ‘fully translated’ functionals, otxc string = ‘ft’+xc string None

class pyscf.mcpdft.otfnal.ftransfnal(ks, **kwargs)[source]#

Bases: transfnal

Parent class of on-top pair-density functional. The main callable is ``eval_ot,’’ which is comparable to pyscf.dft.libxc ``eval_xc.’’ A true ``kernel’’ method, which would take arbitrary 1- and 2-RDMs and return the total PDFT energy, awaits design decisions on how far I’m willing/able to generalize the otpd functions. For instance, in MP2 or CCSD, the 2-RDM spans the whole orbital space and it may not be possible to hold it in memory. At present, it’s all designed around MC-SCF, which is why the ``kernel’’ function that actually calculates the energy is in mcpdft.py instead of here.

Attributes:

mol : object of class pyscf.gto.mole grids : object of class pyscf.dft.gen_grid.Grids eval_ot : function with calling signature shown below _numint : object of class pyscf.dft.NumInt

member functions “hybrid_coeff”, “nlc_coeff, “rsh_coeff”, and “_xc_type” (at least) must be overloaded; see below

otxcstring: name of on-top pair-density exchange-correlation functional

``translated functional’’ of Li Manni et al., JCTC 10, 3669 (2014). The extra attributes are all callables; see their docstrings for more information.

Args:

ksobject of dft.RKS: ks.xc is the Kohn-Sham functional being ``translated’’

Extra attributes for ``fully-translated’’ extension of Carlson et al., JCTC 11, 4077 (2015):

R0float
connecting point to polynomial smoothing function; R0 <= 1.0. Default is 0.9.

R1float
endpoint of polynomial smoothing function, zeta(R1) = zeta’(R1) = zeta’’(R1) = 0.0; R1 >= 1.0. Default is 1.15.

Afloat
Quintic coefficient of polynomial smoothing function. Default = -475.60656009 is chosen to make zeta continuous through its second derivative at given the default R0 and R1.

Bfloat
Quartic coefficient of polynomial smoothing function. Default = -379.47331922 is chosen to make zeta continuous through its second derivative given the default R0 and R1.

Cfloat
Cubic coefficient of polynomial smoothing function. Default = -85.38149682 chosen to make zeta continuous through its second derivative given the default R0 and R1.

property Pi_deriv#

int([x]) -> integer int(x, base=10) -> integer

Convert a number or string to an integer, or return 0 if no arguments are given. If x is a number, return x.__int__(). For floating point numbers, this truncates towards zero.

If x is not a number or if base is given, then x must be a string, bytes, or bytearray instance representing an integer literal in the given base. The literal can be preceded by ‘+’ or ‘-’ and be surrounded by whitespace. The base defaults to 10. Valid bases are 0 and 2-36. Base 0 means to interpret the base from the string as an integer literal. >>> int(‘0b100’, base=0) 4

d_jT_op(x, rho, Pi, **kwargs)[source]#

Evaluate the x.(nabla j) contribution to the second density derivatives of the on-top energy in terms of the untranslated density and pair density

Args:

xndarray of shape (2,*,ngrids): Usually, a functional derivative of the on-top xc energy wrt translated densities
rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Returns: ndarray of shape (*,ngrids)

second derivative of the translation dotted with x 3 rows for tLDA and 5 rows for tGGA

get_rho_translated(Pi, rho)[source]#

Compute the “fully-translated” alpha and beta densities and their derivatives. This is the same as “translated” except

rho’_t^a += zeta’ * rho / 2 rho’_t^b -= zeta’ * rho / 2

And the functional form of “zeta” is changed (see “get_zeta”)

Args:

Pindarray of shape (*, ngrids): containing on-top pair density [and derivatives]
rhondarray of shape (2, *, ngrids): containing spin density [and derivatives]

Returns:

rho_ftndarray of shape (2,*,ngrids): Fully-translated spin density (and derivatives)

get_zeta(R, fn_deriv=1)[source]#

Compute the intermediate zeta used to compute the translated spin densities and its functional derivatives

From the “full” translation [Carlson et al., JCTC 11, 4077 (2015)]: zeta = (1-R)^(1/2) ; R < R0

= A*(R-R1)^5 + B*(R-R1)^4 + C*(R-R1)^3 ; R0 <= R < R1 = 0 ; otherwise

Args:

Rndarray of shape (*,ngrids): Ratio (4Pi/rho^2) and possibly its spatial derivatives Only the first row is used in this function

Kwargs:

fn_derivinteger: order of functional derivative (d^n z / dR^n) to return along with the value of zeta

Returns:

zeta : ndarray of shape (fn_deriv+1, ngrids)

jT_op(x, rho, Pi, **kwargs)[source]#

Evaluate jTx = (x.j)T where j is the Jacobian of the translated densities in terms of the untranslated density and pair density

Args:

xndarray of shape (2,*,ngrids): Usually, a functional derivative of the on-top xc energy wrt translated densities
rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Returns: ndarray of shape (*,ngrids)

Usually, a functional derivative of the on-top pair density exchange-correlation energy wrt to total density and its derivatives. The potential must be spin-symmetric in pair-density functional theory.

transl_prefix = 'ft'#

pyscf.mcpdft.otfnal.get_transfnal(mol, otxc)[source]#

pyscf.mcpdft.otfnal.make_hybrid_fnal(xc_code, hyb, hyb_type=1)[source]#

Convenience function to write “hybrid” xc functional in terms of only one parameter

Args:

xc_codestring: As used in pyscf.dft.libxc. An exception is raised if it is already a hybrid or contains a kinetic-energy functional component.
hybfloat: Parameter(s) defining the “hybridization” which is handled in various ways according to hyb_type

Kwargs:

hyb_typeint or string

The type of hybrid functional. Current options are: - 0 or ‘translation’: Hybrid fnal is

‘hyb*HF + (1-hyb)*x_code, hyb*HF + c_code’. Based on the idea that ‘exact exchange’ of the translated functional corresponds to exchange plus correlation energy of the underlying wave function. Requires len (hyb) == 1.

1 or ‘average’: Hybrid fnal is
‘hyb*HF + (1-hyb)*x_code, hyb*HF + (1-hyb)*c_code’. Based on the idea that hyb = 1 recovers the wave function energy itself. Requires len (hyb) == 1.
2 or ‘diagram’: Hybrid fnal is
‘hyb*HF + (1-hyb)*x_code, c_code’. Based on the idea that the exchange energy of the wave function somehow can be meaningfully separated from the correlation energy. Requires len (hyb) == 1.
3 or ‘lambda’: as in arXiv:1911.11162v1. Based on
existing ‘double-hybrid’ functionals. Requires len (hyb) == 1.
4 or ‘scaling’: Hybrid fnal is
‘a*HF + (1-a)*x_code, a*HF + (1-a**b)*c_code’ where a = hyb[0] and b = 1 + hyb[1]. Based on the scaling inequalities proven by Levy and Perdew in PRA 32, 2010 (1985): E_c[rho_a] < a*E_c[rho] if a < 1 and E_c[rho_a] > a*E_c[rho] if a > 1; BUT E_c[rho_a] ~/~ a^2 E_c[rho], implying that E_c[rho_a] ~ a^b E_c[rho] with b > 1 unknown. Requires len (hyb) == 2.

pyscf.mcpdft.otfnal.make_scaled_fnal(xc_code, hyb_x=0, hyb_c=0, fnal_x=None, fnal_c=None)[source]#

Convenience function to write the xc_code corresponding to a functional of the type

Exc = hyb_x*E_x[Psi] + fnal_x*E_x[rho] + hyb_c*E_c[Psi]

fnal_c*E_c[rho]

where E[Psi] is an energy from a wave function, and E[rho] is a density functional from libxc. The decomposition of E[Psi] into exchange (E_x) and correlation (E_c) components is arbitrary.

Args:

xc_codestring: As used in pyscf.dft.libxc. An exception is raised if it is already a hybrid or contains a kinetic-energy functional component.

Kwargs:

hyb_xfloat: fraction of wave function exchange to be included
hyb_cfloat: fraction of wave function correlation to be included
fnal_xfloat: fraction of density functional exchange to be included. Defaults to 1 - hyb_x.
fnal_cfloat: fraction of density functional correlation to be included. Defaults to 1 - hyb_c.

returns:

xc_codestring: If xc_code has exchange part x_code and correlation part c_code, the return value is ‘fnal_x * x_code + hyb_x * HF, fnal_c * c_code + hyb_c * HF’ You STILL HAVE TO PREPEND ‘t’ OR ‘ft’!!!

class pyscf.mcpdft.otfnal.otfnal(mol)[source]#

Bases: object

Parent class of on-top pair-density functional. The main callable is ``eval_ot,’’ which is comparable to pyscf.dft.libxc ``eval_xc.’’ A true ``kernel’’ method, which would take arbitrary 1- and 2-RDMs and return the total PDFT energy, awaits design decisions on how far I’m willing/able to generalize the otpd functions. For instance, in MP2 or CCSD, the 2-RDM spans the whole orbital space and it may not be possible to hold it in memory. At present, it’s all designed around MC-SCF, which is why the ``kernel’’ function that actually calculates the energy is in mcpdft.py instead of here.

Attributes:

mol : object of class pyscf.gto.mole grids : object of class pyscf.dft.gen_grid.Grids eval_ot : function with calling signature shown below _numint : object of class pyscf.dft.NumInt

member functions “hybrid_coeff”, “nlc_coeff, “rsh_coeff”, and “_xc_type” (at least) must be overloaded; see below

otxcstring: name of on-top pair-density exchange-correlation functional

Pi_deriv = 0#

property dens_deriv#

energy_ot(casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1)#

Compute the on-top energy - the last term in

E_MCPDFT = h_pq l_pq + 1/2 v_pqrs l_pq l_rs + E_ot[rho,Pi]

Args:

ot : an instance of otfnal class casdm1s : ndarray of shape (2, ncas, ncas)

Contains spin-separated one-body density matrices in an active-orbital basis

casdm2ndarray of shape (ncas, ncas, ncas, ncas): Contains spin-summed two-body density matrix in an active- orbital basis
mo_coeffndarray of shape (nao, nmo): Contains molecular orbital coefficients for active-space orbitals. Columns ncore through ncore+ncas give the basis in which casdm1s and casdm2 are expressed.
ncoreinteger: Number of doubly occupied inactive “core” orbitals not explicitly included in casdm1s and casdm2

Kwargs:

max_memoryint or float: maximum cache size in MB default is 2000
hermiint: 1 if 1rdms are assumed hermitian, 0 otherwise

Returnsfloat

The MC-PDFT on-top (nonclassical) energy

eval_ot(rho, Pi, dderiv=0, **kwargs)[source]#

Evaluate the on-dop energy and its functional derivatives

on a grid

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Kwargs:

dderivinteger: Order of derivatives to return

Returns:

eotndarray of shape (ngrids): integrand of the on-top exchange-correlation energy
vot(array_like (rho), array_like (Pi)) or None: first functional derivative of Eot wrt (density, pair- density) and their derivatives
fotndarray of shape (*,ngrids) or None: second functional derivative of Eot wrt density, pair- density, and derivatives; first dimension is lower- triangular matrix elements corresponding to the basis (rho, Pi, |drho|^2, drho'.dPi, |dPi|) stopping at Pi (3 elements) for t-LDA and |drho|^2 (6 elements) for t-GGA.

get_eff_1body(ao, weight, kern, non0tab=None, shls_slice=None, ao_loc=None, hermi=0)#

Contract the kern with d vrho/ dDpq.

Args:

aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to density (vrho)/ If not provided, it is calculated.

Kwargs:

non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether veff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective term corresponding to kernel times the AO’s, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

get_eff_2body(ao, weight, kern, aosym='s4', eff_ao=None)#

get_eff_2body_kl(ao_k, ao_l, weight, kern, symm=False)#

get_feff_1body(rho, Pi, crho, cPi, ao, weight, kern=None, non0tab=None, shls_slice=None, ao_loc=None, hermi=0, **kwargs)#

Get the terms [Delta F]_{pq}

Args:

rhondarray of shape (2,*,ngrids): Spin-density [and derivatives]
Pindarray with shape (*,ngrids): On-top pair density [and derivatives]
crhondarray of shape (2,*,ngrids): Spin-density [and derivatives] to contract the hessian with
cPindarray with shape (*,ngrids): On-top pair density [and derivatives] to contract Hessian with
aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

kernndarray of shape (*,ngrids): the hessian-vector product. If not provided, it is calculated.
non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether feff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective gradient response corresponding to this on-top pair density exchange-correlation functional, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

get_feff_2body(rho, Pi, crho, cPi, ao, weight, aosym='s4', kern=None, fao=None, **kwargs)#

Get the terms [Delta F]_{pqrs}

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
crhondarray of shape (2,*,ngrids): Spin-density [and derivatives] to contract the hessian with
cPindarray with shape (*,ngrids): On-top pair density [and derivatives] to contract Hessian with
aondarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals in which space to calculate the 2-body veff If a list of length 4, the corresponding set of eri-like elements are returned
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

aosymint or str: Index permutation symmetry of the desired integrals. Valid options are 1 (or ‘1’ or ‘s1’), 4 (or ‘4’ or ‘s4’), ‘2ij’ (or ‘s2ij’), and ‘2kl’ (or ‘s2kl’). These have the same meaning as in PySCF’s ao2mo module. Currently all symmetry exploitation is extremely slow and unparallelizable for some reason so trying to use this is not recommended until I come up with a C routine.
kernndarray of shape (*,ngrids): the hessian-vector product. If not provided, it is calculated.
faondarray of shape (*,ngrids,nao,nao) or: (*,ngrids,nao*(nao+1)//2). An intermediate in which the kernel and the k,l orbital indices have been contracted. Overrides kl_symm

Returnseri-like ndarray

The two-body effective gradient response corresponding to this on-top pair density exchange-correlation functional or elements thereof, in the provided basis.

get_veff_1body(rho, Pi, ao, weight, kern=None, non0tab=None, shls_slice=None, ao_loc=None, hermi=0, **kwargs)#

get the derivatives dEot / dDpq Can also be abused to get semidiagonal dEot / dPppqq if you pass the right kern and squared aos/mos

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to density (vrho)/ If not provided, it is calculated.
non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether veff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective potential corresponding to this on-top pair density exchange-correlation functional, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

get_veff_2body(rho, Pi, ao, weight, aosym='s4', kern=None, vao=None, **kwargs)#

get the derivatives dEot / dPijkl

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
aondarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals in which space to calculate the 2-body veff If a list of length 4, the corresponding set of eri-like elements are returned
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

aosymint or str: Index permutation symmetry of the desired integrals. Valid options are 1 (or ‘1’ or ‘s1’), 4 (or ‘4’ or ‘s4’), ‘2ij’ (or ‘s2ij’), and ‘2kl’ (or ‘s2kl’). These have the same meaning as in PySCF’s ao2mo module. Currently all symmetry exploitation is extremely slow and unparallelizable for some reason so trying to use this is not recommended until I come up with a C routine.
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to pair density (vot). If not provided, it is calculated.
vaondarray of shape (*,ngrids,nao,nao) or: (*,ngrids,nao*(nao+1)//2). An intermediate in which the kernel and the k,l orbital indices have been contracted. Overrides kl_symm

Returnseri-like ndarray

The two-body effective potential corresponding to this on-top pair density exchange-correlation functional or elements thereof, in the provided basis.

get_veff_2body_kl(rho, Pi, ao_k, ao_l, weight, symm=False, kern=None, **kwargs)#

get the two-index intermediate Mkl of dEot/dPijkl

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
ao_kndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index k
ao_lndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index l
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

symmlogical: Index permutation symmetry of the desired integral wrt k,l
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to pair density (vot). If not provided, it is calculated.

Returnsndarray of shape (*,ngrids,nao,nao)

or (*,ngrids,nao*(nao+1)//2). An intermediate for calculating the two-body effective potential corresponding to this on-top pair density exchange-correlation functional in the provided basis.

reset(mol=None)[source]#: Discard cached grid data and optionally update the mol

property xctype#

pyscf.mcpdft.otfnal.register_otfnal(xc_code, preset)[source]#

This function registers the new on-top functional if it hasn’t been registered previously. Args:

xc_code: str
The name of the on-top functional to be registered.

preset: dict
The dictionary containing the information about the on-top functional to be registered. xc_base: str

The name of the underylying KS-functional in the libxc library.

ext_params: dict, with LibXC exchange and correlation functional integer ID as key, and
an array-like object containing the functional parameters as value.

hyb: tuple
The hybrid functional parameters.

facs: tuple
The mixing factors.

kwargs: dict
The additional keyword arguments.

pyscf.mcpdft.otfnal.t_eval_xc(ni, xc_code, rho, spin=0, relativity=0, deriv=1, verbose=None)[source]#

For ‘translated’ functionals, otxc string = ‘t’+xc string 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 part sets the correlation functional.

If “,” not appeared in string, the entire string is treated as the name of a compound functional (containing both the exchange and the correlation functional) which was declared in the functional aliases list. The full list of functional aliases can be obtained by calling the function pyscf.dft.xcfun.XC_ALIAS.keys() .

If the string was not found in the aliased functional list, it is treated as X functional.

To input only X functional (without C functional), leave the second part blank. E.g. description=’slater,’ means a functional with LDA contribution only.

To neglect the contribution of X functional (just apply C functional), leave blank in the first part, e.g. description=’,vwn’ means a functional with VWN only.

If compound XC functional is specified, no matter whether 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 names need to be separated by operators “+” or “-”. Blank spaces are ignored. NOTE the parser only reads operators “+” “-” “*”. / is not supported.

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). “HF” can be put in the correlation functional part (after comma). Putting “HF” in the correlation part is the same to putting “HF” in the exchange part.

String “RSH” means range-separated operator. Its format is RSH(omega, alpha, beta). 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 of GGA functional, in which the LDA contribution is included.

Args:

xc_codestr
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.

rhondarray
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:

spinint
spin polarized if spin > 0

relativityint
No effects.

verboseint 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, v2tau2, 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] | v2tau2[:,3] = (u_u, u_d, d_d) | v2rholapl[:,4] | v2rhotau[:,4] = (u_u, u_d, d_u, d_d) | v2lapltau[:,4] | v2sigmalapl[:,6] | v2sigmatau[:,6] = (uu_u, uu_d, ud_u, ud_d, dd_u, dd_d)

kxc for restricted case: (v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,

v3rho2lapl, v3rho2tau, v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, v3rhotau2, v3sigma2lapl, v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, v3lapl2tau, v3lapltau2, 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) | v3rho2lapl[:,6] | v3rho2tau[:,6] = (u_u_u, u_u_d, u_d_u, u_d_d, d_d_u, d_d_d) | v3rhosigmalapl[:,12] | v3rhosigmatau[:,12] = (u_uu_u, u_uu_d, u_ud_u, u_ud_d, u_dd_u, u_dd_d,

d_uu_u, d_uu_d, d_ud_u, d_ud_d, d_dd_u, d_dd_d)

v3rholapl2[:,6]

v3rholapltau[:,8]

v3rhotau2[:,6] = (u_u_u, u_u_d, u_d_d, d_u_u, d_u_d, d_d_d)

v3sigma2lapl[:,12]

v3sigma2tau[:,12] = (uu_uu_u, uu_uu_d, uu_ud_u, uu_ud_d, uu_dd_u, uu_dd_d, ud_ud_u, ud_ud_d, ud_dd_u, ud_dd_d, dd_dd_u, dd_dd_d)

v3sigmalapl2[:,9]

v3sigmalapltau[:,12]

v3sigmatau2[:,9] = (uu_u_u, uu_u_d, uu_d_d, ud_u_u, ud_u_d, ud_d_d, dd_u_u, dd_u_d, dd_d_d)

v3lapl3[:,4]

v3lapl2tau[:,6]

v3lapltau2[:,6]

v3tau3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d)

see also libxc_itrf.c

pyscf.mcpdft.otfnal.t_hybrid_coeff(ni, xc_code, spin=0)[source]#: For ‘translated’ functionals, otxc string = ‘t’+xc string None

pyscf.mcpdft.otfnal.t_nlc_coeff(ni, xc_code)[source]#: For ‘translated’ functionals, otxc string = ‘t’+xc string None

pyscf.mcpdft.otfnal.t_rsh_and_hybrid_coeff(ni, xc_code, spin=0)[source]#

For ‘translated’ functionals, otxc string = ‘t’+xc string 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 | (1-erf(-omega r_{12}))/r_{12} | iq > LR_HFX = < pi | erf(-omega r_{12})/r_{12} | iq > alpha = c_LR beta = c_SR - c_LR

pyscf.mcpdft.otfnal.t_rsh_coeff(ni, xc_code)[source]#: For ‘translated’ functionals, otxc string = ‘t’+xc string None

pyscf.mcpdft.otfnal.t_xc_type(ni, xc_code)[source]#: For ‘translated’ functionals, otxc string = ‘t’+xc string None

class pyscf.mcpdft.otfnal.transfnal(ks, **kwargs)[source]#

Bases: otfnal

Parent class of on-top pair-density functional. The main callable is ``eval_ot,’’ which is comparable to pyscf.dft.libxc ``eval_xc.’’ A true ``kernel’’ method, which would take arbitrary 1- and 2-RDMs and return the total PDFT energy, awaits design decisions on how far I’m willing/able to generalize the otpd functions. For instance, in MP2 or CCSD, the 2-RDM spans the whole orbital space and it may not be possible to hold it in memory. At present, it’s all designed around MC-SCF, which is why the ``kernel’’ function that actually calculates the energy is in mcpdft.py instead of here.

Attributes:

mol : object of class pyscf.gto.mole grids : object of class pyscf.dft.gen_grid.Grids eval_ot : function with calling signature shown below _numint : object of class pyscf.dft.NumInt

member functions “hybrid_coeff”, “nlc_coeff, “rsh_coeff”, and “_xc_type” (at least) must be overloaded; see below

otxcstring: name of on-top pair-density exchange-correlation functional

``translated functional’’ of Li Manni et al., JCTC 10, 3669 (2014). The extra attributes are all callables; see their docstrings for more information.

Args:

ksobject of dft.RKS: ks.xc is the Kohn-Sham functional being ``translated’’

d_jT_op(x, rho, Pi)[source]#

Evaluate the x.(nabla j) contribution to the second density derivatives of the on-top energy in terms of the untranslated density and pair density

Args:

xndarray of shape (2,*,ngrids): Usually, a functional derivative of the on-top xc energy wrt translated densities
rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Returns: ndarray of shape (*,ngrids)

second derivative of the translation dotted with x 3 rows for tLDA and 5 rows for tGGA

eval_ot(rho, Pi, dderiv=1, weights=None, _unpack_vot=True)[source]#

Evaluate the on-dop energy and its functional derivatives

on a grid

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Kwargs:

dderivinteger: Order of derivatives to return

Returns:

eotndarray of shape (ngrids): integrand of the on-top exchange-correlation energy
vot(array_like (rho), array_like (Pi)) or None: first functional derivative of Eot wrt (density, pair- density) and their derivatives
fotndarray of shape (*,ngrids) or None: second functional derivative of Eot wrt density, pair- density, and derivatives; first dimension is lower- triangular matrix elements corresponding to the basis (rho, Pi, |drho|^2, drho'.dPi, |dPi|) stopping at Pi (3 elements) for t-LDA and |drho|^2 (6 elements) for t-GGA.

get_ratio(Pi, rho_avg)[source]#

R = Pi / [rho/2]^2 = Pi / rho_avg^2 An intermediate quantity when computing the translated spin densities

Note this function returns 1 for values and 0 for derivatives for every point where the charge density is close to zero (i.e., convention: 0/0 = 1)

Args:

Pindarray of shape (*,ngrids): Contains on-top pair density on a grid
rho_avgndarray of shape (*,ngrids): Contains the average of the spin-up and spin-down charge densities on a grid, (rho[0]+rho[1])/2

Returns:

Rndarray of shape (*,ngrids): on-top ratio

get_rho_translated(Pi, rho, _fn_deriv=0)[source]#

Compute the “translated” alpha and beta densities: For the unrestricted case, rho = [rho^a, rho^b] Here:

rho^a will have dim of 1,4 or 6 depends on the functional. For MGGA, rho^a = [rho_u,grad_xu, grad_yu, grad_zu, laplacian_u, tau_u] Similar for rho_b.

The translation is done as follows:

rho_t^a = (rho/2) * (1 + zeta) rho_t^b = (rho/2) * (1 - zeta) rho’_t^a = (rho’/2) * (1 + zeta) rho’_t^b = (rho’/2) * (1 - zeta) tau_t^a = (tau/2) * (1 + zeta) tau_t^b = (tau/2) * (1 - zeta)

See “get_zeta” for the meaning of “zeta”

Args:

Pindarray of shape (*, ngrids): containing on-top pair density [and derivatives]
rhondarray of shape (2, *, ngrids): containing spin density [and derivatives]

Kwargs:

_fn_derivinteger: Order of functional derivatives of zeta to compute. In “translated” functionals, no functional derivatives of zeta are used. This kwarg is used for convenience when calling from children classes. It changes the return signature and should not normally be touched by users.

Returns:

rho_tndarray of shape (2,*,ngrids): Translated spin density (and derivatives) in case of LDA or GGAs Translated spin density, derivatives, and kinetic energy density in case of MGGA

get_zeta(R, fn_deriv=0, _Rmax=1)[source]#

Compute the intermediate zeta used to compute the translated spin densities and its functional derivatives

From the original translation [Li Manni et al., JCTC 10, 3669 (2014)]:

zeta = (1-ratio)^(1/2) ; ratio < 1

= 0 ; otherwise

Args:

Rndarray of shape (*,ngrids): Ratio (4Pi/rho^2) and possibly its spatial derivatives Only the first row is used in this function

Kwargs:

fn_derivinteger: order of functional derivative (d^n z / dR^n) to return along with the value of zeta
_Rmaxfloat: maximum value of R for which to compute zeta or its derivatives; columns of zeta with R[0]>_Rmax are zero. This is a hook for the ``fully-translated’’ child class and should not be touched normally.

Returns:

zeta : ndarray of shape (fn_deriv+1, ngrids)

jT_op(x, rho, Pi)[source]#

Evaluate jTx = (x.j)T where j is the Jacobian of the translated densities in terms of the untranslated density and pair density

Args:

xndarray of shape (2,*,ngrids): Usually, a functional derivative of the on-top xc energy wrt translated densities
rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Returns: ndarray of shape (*,ngrids)

Usually, a functional derivative of the on-top pair density exchange-correlation energy wrt to total density and its derivatives. The potential must be spin-symmetric in pair-density functional theory. 2 rows for tLDA, 3 rows for tGGA, and 4 rows for meta-GGA

split_x_c()[source]#

Get one translated functional for just the exchange and one for just the correlation part of the energy.

Returns:

xfnalobject of transfnal: this functional, but only the exchange part
cfnalobject of transfnal: this functional, but only the correlation part

transl_prefix = 't'#

pyscf.mcpdft.otfnal.unregister_otfnal(xc_code)[source]#

This function unregisters the on-top functional if it has been registered previously. Args:

xc_code: str
The name of the on-top functional to be unregistered.

pyscf.mcpdft.otpd module#

pyscf.mcpdft.otpd.density_orbital_derivative(ot, ncore, ncas, casdm1s, cascm2, rho, mo, deriv=0, non0tab=None)[source]#

Compute the half-transformed density and 3/4-transformed pair- density matrix:

D_i(r) = sum_j D_ij phi_j(r) d_i(r) = sum_jkl d_ijkl phi_j(r) phi_k(r) phi_l(r)

so that, for instance, the derivative with respect to orbital rotations of the density and pair density are

d/dk_ij rho(r) = phi_i(r) D_j(r) - phi_j(r) D_k(r) d/dk_ij Pi(r) = i(r) d_j(r) - j(r) d_i(r)

and the derivatives with respect to nuclear displacement are

drho/dRA|(r not in A) = -sum_(mu in A) phi’_mu(r) D_mu(r) drho/dRA|(r in A) = +sum_(mu not in A) phi’_mu(r) D_mu(r) dPi/dRA|(r not in A) = -sum_(mu in A) phi’_mu(r) d_mu(r) dPi/dRA|(r in A) = +sum_(mu not in A) phi’_mu(r) d_mu(r)

There is a mismatch between the ndarray shape and the data layout in the arg mo and the two return arrays. For performance reasons, the grid index should be contiguous in memory for every array with a a grid dimension. However, by convention, in PySCF ndarrays with grid and orbital indices place the latter index in the last position, the former in the second-to-last position, and any other before that in row-major order. That is, for a four-index array of this type, transpose (2,3,1,0) gives a contiguous column-major array, and transpose (0,1,3,2) results in a contiguous row-major array.

Args:

otobject of otfnal: The on-top density functional containing grid information
ncoreinteger: Number of doubly-occupied core orbitals
ncasinteger: Number of active orbitals
casdm1sndarray of shape (2,ncas,ncas): Spin-separated one-body reduced density matrix
cascm2ndarray of shape [ncas,]*4: Spin-summed cumulant of two-body reduced density matrix
rhondarray of shape (2,*,ngrids): Spin-separated density (and derivatives) on a grid
mondarray of shape (*,ngrids,nmo): Molecular orbitals (and derivatives) evaluated on a grid. Data stride must be 0>2>1; i.e., mo.transpose (0,2,1) must be a contiguous row-major array.

Kwargs:

derivinteger: Order of derivatives of half-transformed density to compute; i.e., 0 for LDA and 1 for GGA
non0tab2D boolean array: Mask array to determine whether densities are considered to be numerically zero

Returns:

drhondarray of shape (2,*,ngrids,nmo): Half-transformed density matrix on a grid and its derivatives. Data stride is 0>1>3>2; i.e., drho.transpose (0,1,3,2) is a contiguous row-major array.
dPindarray of shape (*,ngrids,nmo): 3/4-transformed pair density matrix on a grid and its derivatives. Data stride is 0>2>1; i.e., dPi.transpose (0,2,1) is a contiguous row-major array.

pyscf.mcpdft.otpd.get_ontop_pair_density(ot, rho, ao, cascm2, mo_cas, deriv=0, non0tab=None)[source]#

Compute the on-top pair density and its derivatives on a grid:

Pi(r) = i(r)*j(r)*k(r)*l(r)*d_ijkl / 2

= rho[0](r)*rho[1](r) + i(r)*j(r)*k(r)*l(r)*l_ijkl / 2

Args:

ot : on-top pair density functional object rho : ndarray of shape (2,*,ngrids)

Contains spin-separated density [and derivatives]. The dm1s underlying these densities must correspond to the dm1s/dm1 in the expression for cascm2 below.

aondarray of shape (*, ngrids, nao)

contains values of aos [and derivatives]

cascm2ndarray of shape [ncas,]*4

contains spin-summed two-body cumulant density matrix in an active-orbital basis given by mo_cas:

cm2[u,v,x,y] = dm2[u,v,x,y] - dm1[u,v]*dm1[x,y]

dm1s[0][u,y]*dm1s[0][x,v]

dm1s[1][u,y]*dm1s[1][x,v]

where dm1 = dm1s[0] + dm1s[1]. The cumulant (cm2) has no nonzero elements for any index outside the active space, unlike the density matrix (dm2), which formally has elements involving uncorrelated, doubly-occupied ``core’’ orbitals which are not usually computed explicitly:

dm2[i,i,u,v] = dm2[u,v,i,i] = 2*dm1[u,v] dm2[u,i,i,v] = dm2[i,v,u,i] = -dm1[u,v]

mo_casndarray of shape (nao, ncas)

molecular-orbital coefficients for active-space orbitals

Kwargs:

derivderivative order through which to calculate.: deriv > 1 not implemented

non0tab : as in pyscf.dft.gen_grid and pyscf.dft.numint

Returnsndarray of shape (*,ngrids)

The on-top pair density and its derivatives if requested deriv = 0 : value (1d array) deriv = 1 : value, d/dx, d/dy, d/dz deriv = 2 : value, d/dx, d/dy, d/dz, d^2/d|r1-r2|^2_(r1=r2)

pyscf.mcpdft.pdft_eff module#

pyscf.mcpdft.pdft_eff.get_eff_1body(otfnal, ao, weight, kern, non0tab=None, shls_slice=None, ao_loc=None, hermi=0)[source]#

Contract the kern with d vrho/ dDpq.

Args:

aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to density (vrho)/ If not provided, it is calculated.

Kwargs:

non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether veff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective term corresponding to kernel times the AO’s, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

pyscf.mcpdft.pdft_eff.get_eff_2body(otfnal, ao, weight, kern, aosym='s4', eff_ao=None)[source]#

pyscf.mcpdft.pdft_eff.get_eff_2body_kl(ot, ao_k, ao_l, weight, kern, symm=False)[source]#

pyscf.mcpdft.pdft_feff module#

pyscf.mcpdft.pdft_feff.get_feff_1body(otfnal, rho, Pi, crho, cPi, ao, weight, kern=None, non0tab=None, shls_slice=None, ao_loc=None, hermi=0, **kwargs)[source]#

Get the terms [Delta F]_{pq}

Args:

rhondarray of shape (2,*,ngrids): Spin-density [and derivatives]
Pindarray with shape (*,ngrids): On-top pair density [and derivatives]
crhondarray of shape (2,*,ngrids): Spin-density [and derivatives] to contract the hessian with
cPindarray with shape (*,ngrids): On-top pair density [and derivatives] to contract Hessian with
aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

kernndarray of shape (*,ngrids): the hessian-vector product. If not provided, it is calculated.
non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether feff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective gradient response corresponding to this on-top pair density exchange-correlation functional, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

pyscf.mcpdft.pdft_feff.get_feff_2body(otfnal, rho, Pi, crho, cPi, ao, weight, aosym='s4', kern=None, fao=None, **kwargs)[source]#

Get the terms [Delta F]_{pqrs}

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
crhondarray of shape (2,*,ngrids): Spin-density [and derivatives] to contract the hessian with
cPindarray with shape (*,ngrids): On-top pair density [and derivatives] to contract Hessian with
aondarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals in which space to calculate the 2-body veff If a list of length 4, the corresponding set of eri-like elements are returned
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

aosymint or str: Index permutation symmetry of the desired integrals. Valid options are 1 (or ‘1’ or ‘s1’), 4 (or ‘4’ or ‘s4’), ‘2ij’ (or ‘s2ij’), and ‘2kl’ (or ‘s2kl’). These have the same meaning as in PySCF’s ao2mo module. Currently all symmetry exploitation is extremely slow and unparallelizable for some reason so trying to use this is not recommended until I come up with a C routine.
kernndarray of shape (*,ngrids): the hessian-vector product. If not provided, it is calculated.
faondarray of shape (*,ngrids,nao,nao) or: (*,ngrids,nao*(nao+1)//2). An intermediate in which the kernel and the k,l orbital indices have been contracted. Overrides kl_symm

Returnseri-like ndarray

The two-body effective gradient response corresponding to this on-top pair density exchange-correlation functional or elements thereof, in the provided basis.

pyscf.mcpdft.pdft_feff.get_feff_2body_kl(otfnal, rho, Pi, crho, cPi, ao_k, ao_l, weight, symm=False, kern=None, **kwargs)[source]#

get the two-index intermediate Mkl of [Delta cdot F]_{pqrs}

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
crhondarray of shape (2,*,ngrids): Spin-density [and derivatives] to contract the hessian with
cPindarray with shape (*,ngrids): On-top pair density [and derivatives] to contract Hessian with
ao_kndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index k
ao_lndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index l
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

symmlogical: Index permutation symmetry of the desired integral wrt k,l
kernndarray of shape (*,ngrids): the hessian-vector product. If not provided, it is calculated.

Returnsndarray of shape (*,ngrids,nao,nao)

or (*,ngrids,nao*(nao+1)//2). An intermediate for calculating the two-body effective gradient response corresponding to this on-top pair density exchange-correlation functional in the provided basis.

pyscf.mcpdft.pdft_feff.kernel(ot, dm1s, cascm2, c_dm1s, c_cascm2, mo_coeff, ncore, ncas, max_memory=2000, hermi=1, paaa_only=False, aaaa_only=False, jk_pc=False, delta=False)[source]#

Get the 1- and 2-body effective gradient responses from MC-PDFT. The $rho cdot mathbf{F}$ terms, or Hessian vector products.

Args:

ot : an instance of otfnal class dm1s : ndarray of shape (2, nao, nao)

Contains the spin-separated one-body density matrices to evaluate the kernel at

cascm2ndarray of shape (ncas, ncas, ncas, ncas): Spin-summed two-body cumulant density matrix in the active space to evaluate the kernel at
c_dm1sndarray of shape (2, nao, nao): Contains the spin-separated one-body density matrices to contract the kernel with.
c_cascm2ndarray of shape (ncas, ncas, ncas, ncas): Spin-summed two-body cumulant density matrix in the active space to contract the kernel with.
mo_coeffndarray of shape (nao, nmo): containing molecular orbital coefficients
ncoreinteger: number of inactive orbitals
ncasinteger: number of active orbitals

Kwargs:

max_memoryint or float: maximum cache size in MB default is 2000
hermiint: 1 if 1rdms are assumed hermitian, 0 otherwise
paaa_onlylogical: If true, only compute the paaa range of papa and ppaa (all other elements set to zero)
aaaa_onlylogical: If true, only compute the aaaa range of papa and ppaa (all other elements set to zero; overrides paaa_only)
jk_pclogical: If true, compute the ppii=pipi elements of veff2 (otherwise, these are set to zero)
deltalogical: If true, then contract with the delta density. The delta density is the c_dm1s-dm1s and similarly for the 2rdm element (though care is taken since 2rdm elements are expressed in cumulant form).

Returns:

feff1ndarray of shape (nao, nao): 1-body effective gradient response
feff2object of class pdft_eff._ERIS: 2-body effective gradient response

pyscf.mcpdft.pdft_feff.lazy_kernel(ot, dm1s, cascm2, c_dm1s, c_cascm2, mo_cas, hermi=1, max_memory=2000, delta=False)[source]#: 1- and 2-body gradient response (hessian-vector products) from MC-PDFT. This is the lazy way and doesn’t care about memory.

pyscf.mcpdft.pdft_veff module#

pyscf.mcpdft.pdft_veff.get_veff_1body(otfnal, rho, Pi, ao, weight, kern=None, non0tab=None, shls_slice=None, ao_loc=None, hermi=0, **kwargs)[source]#

get the derivatives dEot / dDpq Can also be abused to get semidiagonal dEot / dPppqq if you pass the right kern and squared aos/mos

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
aondarray or 2 ndarrays of shape (*,ngrids,nao): contains values and derivatives of nao. 2 different ndarrays can have different nao but not different ngrids
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to density (vrho)/ If not provided, it is calculated.
non0tabndarray of shape (nblk, nbas): Identifies blocks of grid points which are nonzero on each AO shell so as to exploit sparsity. If you want the “ao” array to be in the MO basis, just leave this as None. If hermi == 0, it only applies to the bra index ao array, even if the ket index ao array is the same (so probably always pass hermi = 1 in that case)
shls_slicesequence of integers of len 2: Identifies starting and stopping indices of AO shells
ao_locndarray of length nbas: Offset to first AO of each shell
hermiinteger or logical: Toggle whether veff is supposed to be a Hermitian matrix You can still pass two different ao arrays for the bra and the ket indices, for instance if one of them is supposed to be a higher derivative. They just have to have the same nao in that case.

Returnsndarray of shape (nao[0],nao[1])

The 1-body effective potential corresponding to this on-top pair density exchange-correlation functional, in the atomic-orbital basis. In PDFT this functional is always spin-symmetric.

pyscf.mcpdft.pdft_veff.get_veff_2body(otfnal, rho, Pi, ao, weight, aosym='s4', kern=None, vao=None, **kwargs)[source]#

get the derivatives dEot / dPijkl

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
aondarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals in which space to calculate the 2-body veff If a list of length 4, the corresponding set of eri-like elements are returned
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

aosymint or str: Index permutation symmetry of the desired integrals. Valid options are 1 (or ‘1’ or ‘s1’), 4 (or ‘4’ or ‘s4’), ‘2ij’ (or ‘s2ij’), and ‘2kl’ (or ‘s2kl’). These have the same meaning as in PySCF’s ao2mo module. Currently all symmetry exploitation is extremely slow and unparallelizable for some reason so trying to use this is not recommended until I come up with a C routine.
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to pair density (vot). If not provided, it is calculated.
vaondarray of shape (*,ngrids,nao,nao) or: (*,ngrids,nao*(nao+1)//2). An intermediate in which the kernel and the k,l orbital indices have been contracted. Overrides kl_symm

Returnseri-like ndarray

The two-body effective potential corresponding to this on-top pair density exchange-correlation functional or elements thereof, in the provided basis.

pyscf.mcpdft.pdft_veff.get_veff_2body_kl(otfnal, rho, Pi, ao_k, ao_l, weight, symm=False, kern=None, **kwargs)[source]#

get the two-index intermediate Mkl of dEot/dPijkl

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]
ao_kndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index k
ao_lndarray of shape (*,ngrids,nao): OR list of ndarrays of shape (,ngrids,) values and derivatives of atomic or molecular orbitals corresponding to index l
weightndarray of shape (ngrids): containing numerical integration weights

Kwargs:

symmlogical: Index permutation symmetry of the desired integral wrt k,l
kernndarray of shape (*,ngrids): the derivative of the on-top potential with respect to pair density (vot). If not provided, it is calculated.

Returnsndarray of shape (*,ngrids,nao,nao)

or (*,ngrids,nao*(nao+1)//2). An intermediate for calculating the two-body effective potential corresponding to this on-top pair density exchange-correlation functional in the provided basis.

pyscf.mcpdft.pdft_veff.kernel(ot, dm1s, cascm2, mo_coeff, ncore, ncas, max_memory=2000, hermi=1, paaa_only=False, aaaa_only=False, jk_pc=False)[source]#

Get the 1- and 2-body effective potential from MC-PDFT.

Args:

ot : an instance of otfnal class dm1s : ndarray of shape (2, nao, nao)

containing spin-separated one-body density matrices

cascm2ndarray of shape (ncas, ncas, ncas, ncas): containing spin-summed two-body cumulant density matrix in an active space
mo_coeffndarray of shape (nao, nmo): containing molecular orbital coefficients
ncoreinteger: number of inactive orbitals
ncasinteger: number of active orbitals

Kwargs:

max_memoryint or float: maximum cache size in MB default is 2000
hermiint: 1 if 1rdms are assumed hermitian, 0 otherwise
paaa_onlylogical: If true, only compute the paaa range of papa and ppaa (all other elements set to zero)
aaaa_onlylogical: If true, only compute the aaaa range of papa and ppaa (all other elements set to zero; overrides paaa_only)
jk_pclogical: If true, compute the ppii=pipi elements of veff2 (otherwise, these are set to zero)

Returns:

veff1ndarray of shape (nao, nao): 1-body effective potential
veff2object of class pdft_eff._ERIS: 2-body effective potential and related quantities

pyscf.mcpdft.pdft_veff.lazy_kernel(ot, dm1s, cascm2, mo_cas, max_memory=2000, hermi=1, veff2_mo=None)[source]#

Get the 1- and 2-body effective potential from MC-PDFT. Eventually I’ll be able to specify mo slices for the 2-body part

Args:

ot : an instance of otfnal class dm1s : ndarray of shape (2, nao, nao)

containing spin-separated one-body density matrices

cascm2ndarray of shape (ncas, ncas, ncas, ncas): containing spin-summed two-body cumulant density matrix in an active space
mo_casndarray of shape (nao, ncas): containing molecular orbital coefficients for active-space orbitals

Kwargs:

max_memoryint or float: maximum cache size in MB default is 2000
hermiint: 1 if 1rdms are assumed hermitian, 0 otherwise

Returnsfloat

The MC-PDFT on-top exchange-correlation energy

pyscf.mcpdft.tfnal_derivs module#

pyscf.mcpdft.tfnal_derivs.contract_fot(otfnal, fot, rho0, Pi0, rho1, Pi1, unpack=True, vot_packed=None)[source]#

Evaluate the product of a packed lower-triangular matrix with perturbed density, pair density, and derivatives.

Args:

fotndarray of shape (*,ngrids): Lower-triangular matrix elements corresponding to the basis (rho, Pi, |drho|^2, drho’.dPi, |dPi|^2) stopping at Pi (3 elements) for *tLDA and |drho|^2 (6 elements) for tGGA.
rho0ndarray of shape (2,*,ngrids): containing density [and derivatives] the density at which fot was evaluated
Pi0ndarray with shape (*,ngrids): containing on-top pair density [and derivatives] the density at which fot was evaluated
rho1ndarray of shape (2,*,ngrids): containing density [and derivatives] the density contracted with fot
Pi1ndarray with shape (*,ngrids): containing on-top pair density [and derivatives] the density contracted with fot

Kwargs:

unpacklogical

If True, returns vot1 in unpacked shape:

(ndarray of shape (*,ngrids),: ndarray of shape (*,ngrids))

corresponding to (density, pair density) and their derivatives. This requires vot_packed for *tGGA functionals Otherwise, returns vot1 in packed shape:

(rho, Pi, |rho’|^2, rho’.Pi’, |Pi’|^2)

stopping at Pi (3 elements) for *tLDA and |rho’|^2 (6 elements) for tGGA.

vot_packedndarray of shape (*,ngrids)

Returns:

vot1(ndarray of shape (*,ngrids),: ndarray of shape (*,ngrids))

product of fot wrt (density, pair density) and their derivatives

pyscf.mcpdft.tfnal_derivs.contract_vot(vot, rho, Pi)[source]#

Evalute the product of unpacked vot with perturbed density, pair density, and derivatives.

Args:

vot(ndarray of shape (,ngrids), ndarray of shape (, ngrids)): format is ([a, ngrids], [b, ngrids]) : (vrho, vPi) ftGGA: a=4, b=4 tGGA: a=4, b=1 *tLDA: a=1, b=1
rhondarray of shape (*,ngrids): containing density [and derivatives] the density contracted with vot
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives] the density contracted with vot

Returns:

cvotndarray of shape (ngrids): product of vot wrt (density, pair density) and their derivatives

pyscf.mcpdft.tfnal_derivs.eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True)[source]#

get the integrand of the on-top xc energy and its functional derivatives wrt rho and Pi

Args:

rhondarray of shape (2,*,ngrids): containing spin-density [and derivatives]
Pindarray with shape (*,ngrids): containing on-top pair density [and derivatives]

Kwargs:

dderivinteger: Order of derivatives to return
weightsndarray of shape (ngrids): used ONLY for debugging the total number of ``translated’’ electrons in the calculation of rho_t Not multiplied into anything!
_unpack_votlogical: If True, derivatives with respect to density gradients are reported as de/drho’ and de/dPi’; otherwise, they are reported as de/d|rho’|^2, de/d(rho’.Pi’), and de/d|Pi’|^2

Returns:

eotndarray of shape (ngrids): integrand of the on-top exchange-correlation energy
votndarrays of shape (*,ngrids) or None: first functional derivative of Eot wrt (density, pair density) and their derivatives. If _unpack_vot = True, shape and format is ([a, ngrids], [b, ngrids]) : (vrho, vPi); otherwise, [c, ngrids] : [rho,Pi,|rho’|^2,tau,rho’.Pi’,|Pi’|^2] ftGGA: a=4, b=4, c=5 (drop tau) tmGGA: a=5, b=1, c=4 (drop Pi’) tGGA: a=4, b=1, c=3 (drop Pi’, tau) *tLDA: a=1, b=1, c=2 (drop rho’, tau)
fotndarray of shape (*,ngrids) or None: second functional derivative of Eot wrt density, pair density, and derivatives; first dimension is lower- triangular matrix elements corresponding to the basis (rho, Pi, |rho’|^2, rho’.Pi’, |Pi’|^2) stopping at Pi (3 elements) for *tLDA and |rho’|^2 (6 elements) for tGGA.

pyscf.mcpdft.tfnal_derivs.unpack_vot(packed, rho, Pi)[source]#

pyscf.mcpdft.xmspdft module#

pyscf.mcpdft.xmspdft.diagonalize_safock(mc, mo_coeff=None, ci=None)[source]#: Diagonalizes the SA-Fock matrix. Returns the eigenvalues and eigenvectors.

pyscf.mcpdft.xmspdft.fock_h1e_for_cas(mc, sa_casdm1, mo_coeff=None, ncas=None, ncore=None)[source]#

Compute the CAS SA Fock Matrix 1-electron integrals.

Args:

mc : instance of class _PDFT

sa_casdm1ndarray of shape (ncas,ncas): 1-RDM in the active space generated from the state-average density.
mo_coeffndarray of shape (nao,nmo): A full set of molecular orbital coefficients. Taken from self if not provided.
ncasint: Number of active space molecular orbitals
ncoreint: Number of core molecular orbitals

Returns:

A tuple, the first is the one-electron integrals defined in the CAS space and the second is the constant, core part.

pyscf.mcpdft.xmspdft.make_fock_mcscf(mc, mo_coeff=None, ci=None, weights=None)[source]#

Compute the SA-Fock Hamiltonian/Matrix

Args:

mc : instance of class _PDFT

mo_coeffndarray of shape (nao,nmo): A full set of molecular orbital coefficients. Taken from self if not provided.
cindarray of shape (nroots): CI vectors should be from a converged CASSCF/CASCI calculation
weightsndarray of length nroots: Weight for each state. If none, uses weights from SA-CASSCF calculation

Returns:

safock_hamndarray of shape (nroots,nroots): State-average Fock matrix expressed in the basis provided by the CI vectors.

pyscf.mcpdft.xmspdft.safock_energy(mc, **kwargs)[source]#

Diabatizer Function

The “objective” function we are optimizing when solving for the SA-Fock eigenstates is that the SA-Fock energy (average) is minimized with the constraint to the final states being orthonormal. Its the whole saddle point thing similar to in SA-MC-PDFT gradients. For now I just omit this whole issue since the first derivative is by default 0 and don’t worry about the Hessian (or second derivatives) since I don’t care. It should fail if we try and do gradients but I can put a hard RaiseImplementation error

Returns:

SA-Fock Energyfloat: weighted sum of SA-Fock energies
dSA-Fock Energyndarray of shape npair = nroots*(nroots - 1)/2: first derivative of the SA-Fock energy wrt interstate rotation This is zero by default since we diagonalize a matrix
d2SA-Fock Energyndarray of shape (npair,npair): Should be the Lagrange multiplier terms. Currently returning None since we cannot do gradients yet.

pyscf.mcpdft.xmspdft.solve_safock(mc, mo_coeff=None, ci=None)[source]#

Diabatize Function. Finds the SA-Fock Eigenstates within the model space spanned by the CI vectors.

Args:

mc : instance of class _PDFT

mo_coeffndarray of shape (nao,nmo): A full set of molecular orbital coefficients. Taken from self if not provided.
cindarray of shape (nroots): CI vectors should be from a converged CASSCF/CASCI calculation

Returns:

convbool: Always true. If there is a convergence issue, scipy will raise an exception.
cilist of ndarrays of length = nroots: CI vectors of the optimized diabatic states

Module contents#

Multi-configuration pair-density functional theory#

Simple usage:

>>> from pyscf import gto, scf, mcpdft
>>> mol = gto.M(atom='N 0 0 0; N 0 0 1', basis='def2-tzvp')
>>> mf = scf.RHF(mol).run ()
>>> mc = mcpdft.CASSCF (mf, 'tPBE', 6, 6)
>>> mc.run()

pyscf.mcpdft.CASCI(mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, **kwargs)#

pyscf.mcpdft.CASCIPDFT(mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, **kwargs)[source]#

pyscf.mcpdft.CASSCF(mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, frozen=None, **kwargs)#

pyscf.mcpdft.CASSCFPDFT(mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, frozen=None, **kwargs)[source]#

class pyscf.mcpdft.MultiStateMCPDFTSolver[source]#: Bases: object