Modules

Wu and Yang

Summary This is a program for preforming Kohn-Sham inversion with algorithm introduced by Wu and Yang [WuYang2003].

Written for Python 3.7.4 Lots of other text here, like details about how this uses the original Wu Yang idea [WuYang2003], and how the other ref is used for Regularization [HBY2007] . Any useful info about the algorithm should also be given.

Module author: Seungsoo Nam <skaclitz@yonsei.ac.kr> <http://tccl.yonsei.ac.kr/mediawiki/index.php/Main_Page> ORCID: 000-0001-9948-6140 Ryan J. McCarty <rmccarty@uci.edu> <http://ryanjmccarty.com> ORCID: 0000-0002-2417-8917

References
WuYang2003(1,2)

Qin Wu and Weitao Yang. A direct optimization method for calculating density functionals and exchange–correlation potentials from electron densities (2003) <https://doi.org/10.1063/1.1535422> Journal of Chemical Physics, 118(6).

HBY2007

Tim Heaton-Burgess, Felipe A. Bulat, and Weitao Yang. Optimized Effective Potentials in Finite Basis Sets (2007) <https://doi.org/10.1103/PhysRevLett.98.256401> Physical review letters, 98(25).

GNGK2020(1,2,3)

Rachel Garrick, Amir Natan, Tim Gould, and Leeor Kronik. Exact Generalized Kohn-Sham Theory for Hybrid Functionals (2020) <https://doi.org/10.1103/PhysRevX.10.021040> Physical Review X, 10(2).

Funding

This research was made possible by funding from the National Research Foundation of Korea (NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737).

wy.project_b(mw1, mw2)[source]

Summary: Transfer b vector from one WY calculation (mw1) to other (mw2)

Parameters

  • mw1 – WY object to project b

  • mw2 – WY object to be projected

wy.numint_3c2b(mol, pbas, level=5)[source]

Summary: Three-center overlap integral with different atomic- and potential- basis sets with numerical integration

Parameters

  • mol – an instance of Mole

  • pbas (list) – potential basis set for WY, given as a list of functions

  • level (int 0~9) – PySCF preset mesh grid used for numerical integration. Default is 5

Returns

Sijt three-center overlap integral with shape ((n1, n1, n2)) where n1 is the length of atomic orbital basis set defined in mol, and n2 is the length of potential basis set

Return type

(ndarray)

wy.time_profile(mw)[source]

Summary: Print to terminal a collection of time data in seconds

Terminal prints:

  • initialize : time for initialization

  • total_opt : total elasped time

  • solve_eig : time to construct fock matrix and diagonalize it

  • Ws_eval : time to evaluate objective function Ws

  • Gd_eval : time to evaluate gradient of objective function

  • Hs_eval : time to evaluate hessian of objective function

wy.run(mw, xbas=None)[source]

Summary: Minimize an objective function -Ws

Parameters

mw – RWY or UWY object

wy.wyscf(mw, ddmtol=1e-06)[source]

Summary: Run WY until its self-consistent. This is required for hybrid guiding potentials

Note

This function performs very simplied generalization of WY for non-local guiding potential, mimicking [GNGK2020]. There is no theoretical proof that the result of this function can be trusted. The optimization of WY objective functional is only theoretically proven for local multiplicative potentials, not non-local potential.

Parameters

  • mw – RWY or UWY object

  • ddmtol (float, optional) – Convergence criteria of density matrix.

wy.info(mw)[source]

Summary: Print summary of minimization

Parameters

mw – RWY or UWY object

wy.basic(mw, mol, pbas, Sijt, tbas, Smnt)[source]

Summary: Common basic initialization function for RWY and UWY objects

Parameters

  • mw – RWY or UWY object

  • mol – an instance of Mole

  • pbas (dict or str) – potential basis set for WY

  • Sijt (ndarray) – three-center overlap integral with shape ((no, no, np)) where no is the length of atomic orbital basis set defined in mol, while np is the length of potential basis set Only effective for model systems and molecular systems with user-defined pbs. Otherwise, this is ignored and analytical integration is performed.

  • tbas (dict or str) – basis set of target density matrix

  • Smnt (ndarray) – three-center overlap integral with shape ((nt, nt, np)) where nt is the length of atomic orbital basis set that defines target density matrix, while np is the length of potential basis set Only effective for model systems and molecular systems with user-defined pbs. Otherwise, this is ignored and analytical integration is performed.

class wy.RWY(mol, dm_tar, pbas=None, Sijt=None, tbas=None, Smnt=None, dm_aux=None)[source]

Bases: object

Summary: Perform WY calculation in restricted scheme

Variables

  • mol (object) – an instance of Mole

  • dm_tar (ndarray) – Density matrix of target density in atomic orbital basis representation

  • pbas (dict or str or list) – Potential basis set for WY. If not given, same with atomic orbital basis. If this is given as a list of function, those functions are recognized as potential basis set.

  • Sijt (ndarray) – Three-center overlap integral. Ignored (calculated analitically) when pbas is given as dict or str (i.e., pyscf supported formats)

  • dm_aux (ndarray) – Auxilary density matrix to construct density-dependent part of guiding potential. Default is dm_tar

  • guide (str) –

    Guiding potential. Can be set as

    None : no guiding potential except external potential
    ’faxc’ : Exchange-correlation part of Fermi-Amaldi potential
    xc : ks.xc attribute in pyscf DFT

    or any combination of them, for example: b3lyp-0.2*hf+0.2*faxc

  • reg (float) – strength of regularization. Default is 0

  • tol (float) – tolarance of optimization (maximum gradient element). Default is 1e-6

  • method (str) – optimization algorithm used in SciPy. Can be (case insensitive) ‘cg’, ‘bfgs’, ‘newton-cg’, ‘dogleg’, ‘trust-ncg’, ‘trust-krylov’, ‘trust-exact’. Defult is ‘trust-exact’

Ivar

  • converged (bool) – optimization converged or not

  • mo_energy (ndarray) – molecular orbital energies

  • mo_occ (ndarray) – molecular orbital occupation numbers

  • mo_coeff (ndarray) – molecular orbital coefficients

  • dm (ndarray) – density matrix in atomic orbital basis representation

  • b (ndarray) – optimized linear combination coefficients

run(xbas=None)

Summary: Minimize an objective function -Ws

Parameters

mw – RWY or UWY object

time_profile()

Summary: Print to terminal a collection of time data in seconds

Terminal prints:

  • initialize : time for initialization

  • total_opt : total elasped time

  • solve_eig : time to construct fock matrix and diagonalize it

  • Ws_eval : time to evaluate objective function Ws

  • Gd_eval : time to evaluate gradient of objective function

  • Hs_eval : time to evaluate hessian of objective function

info()

Summary: Print summary of minimization

Parameters

mw – RWY or UWY object

wyscf(ddmtol=1e-06)

Summary: Run WY until its self-consistent. This is required for hybrid guiding potentials

Note

This function performs very simplied generalization of WY for non-local guiding potential, mimicking [GNGK2020]. There is no theoretical proof that the result of this function can be trusted. The optimization of WY objective functional is only theoretically proven for local multiplicative potentials, not non-local potential.

Parameters

  • mw – RWY or UWY object

  • ddmtol (float, optional) – Convergence criteria of density matrix.

initialize(xbas=None)[source]

Summary: Construct a fixed part of fock matrix F0

\[F_{0} = T + V_{ext} + V_{h} + V_{g}\]
solve(b)[source]

Summary: Under a given b vector, construct a fock matrix and diagonalize it

F = F0 + V(b)

FC = SCE

resulting mo_coeff(C), mo_energy(E), and density matrix are stored as instance attributes

eval_Ws(b)[source]

Summary: Calculate objective function Ws under given b vector

When mw.reg > 0, regularization is added

eval_Gd(b)[source]

Summary: -d(Ws)/b_t

eval_Hs(b)[source]

Summary: -d^2(Ws)/(b_t)(b_u)

Dvb(b=None)[source]

Summary: Calcuate the norm of the correction potential derivative

\[\| \nabla v_C ( \textbf{r} ) \|^2 = \int v_C ( \textbf{r} ) \nabla^2 v_C ( \textbf{r} ) d \textbf{r}\]
class wy.UWY(mol, dm_tar, pbas=None, Sijt=None, tbas=None, Smnt=None, dm_aux=None)[source]

Bases: object

Summary: Perform WY calculation in unrestricted scheme

Variables

  • mol (object) – an instance of Mole

  • dm_tar (ndarray) – Density matrix of target density in atomic orbital basis representation

  • pbas (dict or str or list) – Potential basis set for WY. If not given, same with atomic orbital basis. If this is given as a list of function, those functions are recognized as potential basis set.

  • Sijt (ndarray) – Three-center overlap integral. Ignored (calculated analitically) when pbas is given as dict or str (i.e., pyscf supported formats)

  • dm_aux (ndarray) – Auxilary density matrix to construct density-dependent part of guiding potential. Default is dm_tar

  • guide (str) –

    Guiding potential. Can be set as

    None : no guiding potential except external potential
    ’faxc’ : Exchange-correlation part of Fermi-Amaldi potential
    xc : ks.xc attribute in pyscf DFT

    or any combination of them, for example: b3lyp-0.2*hf+0.2*faxc

  • reg (float) – strength of regularization. Default is 0

  • tol (float) – tolarance of optimization (maximum gradient element). Default is 1e-6

  • method (str) – optimization algorithm used in SciPy. Can be (case insensitive) ‘cg’, ‘bfgs’, ‘newton-cg’, ‘dogleg’, ‘trust-ncg’, ‘trust-krylov’, ‘trust-exact’. Defult is ‘trust-exact’

Ivar

  • converged (bool) – optimization converged or not

  • mo_energy (ndarray) – molecular orbital energies

  • mo_occ (ndarray) – molecular orbital occupation numbers

  • mo_coeff (ndarray) – molecular orbital coefficients

  • dm (ndarray) – density matrix in atomic orbital basis representation

  • b (ndarray) – optimized linear combination coefficients.

run(xbas=None)

Summary: Minimize an objective function -Ws

Parameters

mw – RWY or UWY object

time_profile()

Summary: Print to terminal a collection of time data in seconds

Terminal prints:

  • initialize : time for initialization

  • total_opt : total elasped time

  • solve_eig : time to construct fock matrix and diagonalize it

  • Ws_eval : time to evaluate objective function Ws

  • Gd_eval : time to evaluate gradient of objective function

  • Hs_eval : time to evaluate hessian of objective function

info()

Summary: Print summary of minimization

Parameters

mw – RWY or UWY object

wyscf(ddmtol=1e-06)

Summary: Run WY until its self-consistent. This is required for hybrid guiding potentials

Note

This function performs very simplied generalization of WY for non-local guiding potential, mimicking [GNGK2020]. There is no theoretical proof that the result of this function can be trusted. The optimization of WY objective functional is only theoretically proven for local multiplicative potentials, not non-local potential.

Parameters

  • mw – RWY or UWY object

  • ddmtol (float, optional) – Convergence criteria of density matrix.

initialize(xbas=None)[source]

Summary: Construct a fixed part of fock matrix F0

\[F_{0} = T + V_{ext} + V_{h} + V_{g}\]
solve(b)[source]

Summary: Under a given b vector, construct a fock matrix and diagonalize it

F = F0 + V(b)

FC = SCE

resulting mo_coeff(C), mo_energy(E), and density matrix are stored as instance attributes

eval_Ws(b)[source]

Summary: Calculate objective function Ws under given b vector

When mw.reg > 0, regularization is added

eval_Gd(b)[source]

Summary: -d(Ws)/b_t

eval_Hs(b)[source]

Summary: -d^2(Ws)/(b_t)(b_u)

Dvb(b=None)[source]

Summary: Calcuate the norm of correction potential derivative

\[\| \nabla v_C ( \textbf{r} ) \| = \frac{1}{2} \sum_{\sigma} \int v_{C,\sigma} ( \textbf{r} ) \nabla^2 v_{C,\sigma} ( \textbf{r} ) d \textbf{r}\]

Zhao-Morrison-Parr

Summary This script preforms a Zhao-Morrison-Parr [ZMP1994] Kohn Sham inversion.

Written for Python 3.7.4

References
ZMP1994(1,2)

Qingsheng Zhao, Robert C Morrison, and Robert G Parr. From electron densities to Kohn-Sham kinetic energies, orbital energies, exchange-correlation potentials, and exchange-correlation energies. (1994) <https://doi.org/10.1103/PhysRevA.50.2138> Physical Review A, 50(3) 2138.

THG1997

David J Tozer, Nicholas C Handy, and William H Green. Exchange-correlation functionals from ab initio electron densities (1997) <https://doi.org/10.1016/S0009-2614(97)00586 and-1> Chemical Physics Letters, 273(3-4) 183-194

Pulay1980

P Pulay. Convergence acceleration of iterative sequences. the case of SCF iteration (1980) <https://doi.org/10.1016/0009-2614(80)80396-4> Chemical Physics Letters, 73 (2): 393–398

Module author: Seungsoo Nam <skaclitz@yonsei.ac.kr> <http://tccl.yonsei.ac.kr/mediawiki/index.php/Main_Page> ORCID: 000-0001-9948-6140

Funding

This research was made possible by funding from the National Research Foundation of Korea (NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737).

class zmp.DIIS(S, diis_space)[source]

Bases: object

Summary: Class for DIIS extrapolation used in ZMP

See [Pulay1980] for some extra context.

extrapolate(iteration, fock, dm)[source]

Summary: New fock matrix by linear combination of previous fock matrices

Parameters

  • iteration (integer) – present SCF iteration

  • fock (ndarray) – fock matrix

  • dm (ndarray) – density matrix obtained from previous step

Returns

tuple containing:

(ndarray): newfock extrapolated fock matrix

(float): diis_error DIIS error used in convergence criteria

Return type

(tuple)

zmp.basic(mz, mol)[source]

Summary: Common basic initialization function for RZMP and UZMP objects

Parameters

  • mz – RZMP or UZMP object

  • mol (object) – an instance of Mole

class zmp.RZMP(mol, dm_tar, dm_aux=None)[source]

Bases: object

Summary: Perform ZMP calculation in restricted scheme, see [ZMP1994] for detail.

Variables

  • mol (object) – an instance of Mole

  • dm_tar (ndarray) – Density matrix of target density in atomic orbital basis representation

  • dm_aux (ndarray) – Auxilary density matrix to construct a fixed part of fock matrix. Default is dm_tar

  • guide (str) –

    Guiding potential. Can be set as

    None : no guiding potential except external potential
    ’faxc’ : Exchange-correlation part of Fermi-Amaldi potential
    xc : ks.xc attribute in pyscf DFT

    or any combination of them, for example: b3lyp-0.2*hf+0.2*faxc

  • diis_space (int) – DIIS space size. Default is 40

  • level_shift (float) – Level shift (in AU) for virtual space. Default is 0.2

  • max_cycle (int) – max number of zscf iterations. Defalut is 400

  • conv_tol_dm (float) – converge threshold for density matrix. Default is 1e-7

  • conv_tol_diis (float) – converge threshold for DIIS error. Default is 1e-5

Ivar

  • converged (bool) – zscf converged or not

  • mo_energy (ndarray) – molecular orbital energies (Note that energies when level_shift = 0)

  • mo_occ (ndarray) – molecular orbital occupation numbers

  • mo_coeff (ndarray) – molecular orbital coefficients

  • dm (ndarray) – density matrix in atomic orbital basis representation

  • l (float) – the last given lambda

initialize()[source]

Summary: Construct a fixed part of fock matrix F0

\[F_{0} = T + V_{ext} + V_{h} + V_{g}\]
zscf(l)[source]

Summary: Run self-consistent ZMP equation under given lambda (l). Prints lambda, HOMO-LUMO, dN and C to terminal.

Parameters

l (float) – Lagrange multiplier lambda

class zmp.UZMP(mol, dm_tar, dm_aux=None)[source]

Bases: object

Summary: Perform ZMP calculation in unrestricted scheme, see [THG1997].

Variables

  • mol (object) – an instance of Mole

  • dm_tar (ndarray) – Density matrix of target density in atomic orbital basis representation

  • dm_aux (ndarray) – Auxilary density matrix to construct a fixed part of fock matrix. Default is dm_tar

  • guide (str) –

    Guiding potential. Can be set as

    None : no guiding potential except external potential
    ’faxc’ : Exchange-correlation part of Fermi-Amaldi potential
    xc : ks.xc attribute in pyscf DFT

    or any combination of them, for example: b3lyp-0.2*hf+0.2*faxc

  • diis_space (int) – DIIS space size. Default is 40

  • level_shift (float) – Level shift (in AU) for virtual space. Default is 0.2

  • max_cycle (int) – max number of zscf iterations. Defalut is 400

  • conv_tol_dm (float) – converge threshold for density matrix. Default is 1e-7

  • conv_tol_diis (float) – converge threshold for DIIS error. Default is 1e-5

Ivar

  • converged (bool) – zscf converged or not

  • mo_energy (ndarray) – molecular orbital energies (Note that energies when level_shift = 0)

  • mo_occ (ndarray) – molecular orbital occupation numbers

  • mo_coeff (ndarray) – molecular orbital coefficients

  • dm (ndarray) – density matrix in atomic orbital basis representation

  • l (float) – the last given lambda

initialize()[source]

Summary: Construct a fixed part of fock matrix F0

\[F_{0} = T + V_{ext} + V_{h} + V_{g}\]
zscf(l)[source]

Summary: Run self-consistent ZMP equation under given lambda (l)

Args: l (float): Lagrange multiplier lambda

Utility Functions

Summary Utility functions for KSPies.

References
Mirko2014

Mirko Franchini, Pierre Herman Theodoor Philipsen, Erik van Lenthe, and Lucas Visscher. Accurate Coulomb Potentials for Periodic and Molecular Systems through Density Fitting. (2014) <https://doi.org/10.1021/ct500172n> Journal of Chemical Theory and Computation, 10(5), 1994-2004.

Module author: Seungsoo Nam <skaclitz@yonsei.ac.kr> <http://tccl.yonsei.ac.kr/mediawiki/index.php/Main_Page> ORCID: 000-0001-9948-6140 Hansol Park <hansol7954@yonsei.ac.kr> <http://tccl.yonsei.ac.kr/mediawiki/index.php/Main_Page> ORCID: 000-0001-8706-5015

functions

mo2ao

wfnreader

eval_vh

eval_vxc

Todo

  • Make IF statment for kspies_fort

Internal Log

2020-06-06 SN made edits

2020-07-26 Updates on eval_vh (Hansol ver.)

2020-08-21 SN corrected typos, minor changes in attribute names, etc.

2020-11-11 Added wfnreader

util.mo2ao(mo, p1, p2=None)[source]

Summary: Convert mo-basis density matrices to basis-set representation density matrices

Parameters

  • mo (ndarray) – molecular orbital coefficients

  • p1 (ndarray) – mo-basis one-particle density matrices

  • p2 (ndarray) – mo-basis two-particle density matrices, optional

Returns

tuple containing:

(ndarray): dm1 ao-basis one-particle density matrices

(ndarray): dm2 ao-basis two-particle density matrices, returned only when p2 is given

Return type

(tuple)

util.readwfn(filename, mol, makerdm=False)[source]

Summary: Convert wfn file to PySCF format

Parameters

  • filename (string) – name of wfn file to convert

  • mol (object) – an instance of Mole

Returns

tuple containing:

(ndarray): mo_coeff molecular orbital coefficient from wfn

(ndarray): mo_occ molecular orbital occupation number from wfn

(ndarray): mo_energy molecular orbital energy from wfn

Return type

(tuple)

util.parse_guide(description)[source]

Summary: Guiding potential parser for ZMP and WY

Parameters

description (str) – guiding potential description for inversion

Returns

tuple containing:

(float): fac_faxc factor for Fermi-Amaldi potential (faxc)

(string): dft_xc description of dft part of xc

Return type

(tuple)

util.eval_vh(mol, coords, dm, Lvl=3, ang_lv=2)[source]

Summary: Calculate real-space Hartree potential from given density matrix. See [Mirko2014] for some extra context.

Parameters

  • mol (object) – an instance of Mole

  • coords (ndarray) – grids space used for calculating Hartree potential

  • dm (ndarray) – one-particle reduced density matrix in basis set representation

  • Lvl (integer) – Interpolation grids space level (input : 0 ~ 9, default 3)

  • ang_lv (integer) – setting a limit of angular momentum of spherical harmonics (input : 0 ~ 4, default 2)

Returns

vh Pointwise Hartree potential

Return type

(ndarray)

util.eval_vxc(mol, dm, xc_code, coords, delta=1e-07)[source]

Summary: Calculate real-space exchange-correlation potential for GGA from given density matrix

Parameters

  • mol (object) – an instance of Mole

  • coords (ndarray) – grids space used for calculating XC potential

  • dm (ndarray) – one-particle reduced density matrix in basis set representation

  • xc_code (str) – XC functional description.

  • delta (float) – amount of finite difference to calculate numerical differentiation, default is 1e-7 a.u.

Returns

tuple containing:

(ndarray) Pointwise XC potential, vxc(r) for RKS dm, vxc_alpha(r) for UKS dm

(ndarray) Pointwise XC potential, vxc(r) for RKS dm, vxc_beta(r) for UKS dm

Return type

(tuple)