lumeq.opt package

Submodules

lumeq.opt.grassmann module

class Grassmann(mf=None, S=None, Q=None, Y=None, P=None, update_method='steepest_descent')

Bases: object

barzilai_borwein_step(G0, G1, T0)

init_guess()

kernel(C=None, tol=8)

polak_ribiere_step(G0, G1, T0)

class Involution_Grassmann(mf=None, S=None, Q=None, Y=None, P=None, update_method='steepest_descent')

Bases: Grassmann

init_guess()

project_tangent_space(V, G)

update(C, T)

class Projection_Grassmann(mf=None, S=None, Q=None, Y=None, P=None, update_method='steepest_descent')

Bases: Grassmann

init_guess()

project_tangent_space(P, G)

class Quotient_Grassmann(mf=None, S=None, Q=None, Y=None, P=None, update_method='steepest_descent')

Bases: Grassmann

conjugate_gradient_step(Y, G0, G1, T0)

init_guess()

project_tangent_space(Y, G)

tangent_parallel_transport(Y, G, T)

update(Y, T)

cs_decompose(T, full=True, scale=1.0)

geodesic_exp(T, full=True, scale=1.0)

geodesic_svd(B, scale=1.0)

geodesic_svd_compact(Y, T)

lumeq.opt.monte_carlo module

diffusion_walker(potential, n_walkers=1000, dt=0.01, n_steps=500, seed=None)

Run a simple diffusion Monte Carlo walker simulation.

The algorithm alternates drift-diffusion proposals, branching weights, and walker resampling with a reference-energy update.

Parameters:

• potential (callable) – Potential-energy function.

• n_walkers (int) – Number of walkers.

• dt (float) – Time step.

• n_steps (int) – Number of diffusion steps.

• seed (int, optional) – Random seed for reproducibility.

Returns:

Final walker positions and estimated energies.

Return type:

tuple

Notes

The diffusion Monte Carlo propagation follows - d psi / dt = - (T + V - E_ref) psi.

The update is split into three stages:

1. drift-diffusion move from the kinetic term: x' = x + chi * sqrt(dt),

2. branching weight from the potential term: w = exp(-(V(x') - E_ref) * dt),

3. resampling walkers according to the branching weights.

States with local energy closer to E_ref are therefore favored.

importance_sampling(log_prob, proposal_log_prob)

Importance sampling generates samples from a target distribution p(x) using a proposal distribution q(x) that is easier. <f> = int f(x) p(x) dx = int f(x) (p(x)/q(x)) q(x) dx = <f(x) w(x)> where w(x) = p(x)/q(x) are the importance weights.

Parameters:

• log_prob (callable) – Logarithmic target probability density function.

• proposal_log_prob (callable) – Logarithmic proposal density function.

Returns:

Importance weights.

Return type:

numpy.ndarray

log_prob(action)

Abstract function to compute the logarithmic probability density p(x) = e^{-S(x)} such that log p(x) = -S(x)

metropolis(log_prob, x, step_size=1.0, n_steps=5000, seed=None)

Metropolis sampler generates new points based on Gaussian proposals (random walk) provided a target distribution with the acceptance ratio: A = min(1, p(x’)/p(x))

Parameters:

• log_prob (callable) – Logarithmic probability density function.

• x – Initial position.

• step_size (float) – Step size for the proposal distribution.

• n_steps (int) – Number of steps to sample.

• seed (int, optional) – Random seed for reproducibility.

Returns:

Sampled positions.

Return type:

numpy.ndarray

metropolis_pi(log_prob, x, step_size=0.5, n_steps=5000, seed=None, **kwargs)

Metropolis walker for path integral method generates new points based on Gaussian proposals (random walk) provided a target distribution with the acceptance ratio: A = min(1, p(x’)/p(x))

Parameters:

• log_prob (callable) – Logarithmic probability density function.

• x – Initial bead positions.

• step_size (float) – Step size for bead updates.

• n_steps (int) – Number of steps to sample.

• seed (int, optional) – Random seed for reproducibility.

• **kwargs – Additional keyword arguments passed to log_prob.

Returns:

Sampled bead positions.

Return type:

numpy.ndarray

run_dmc(model, n_walkers=3000, dt=0.01, n_steps=7000, n_burn=2000)

run a Diffusion Monte Carlo simulation

run_pimc(model, beta=1.0, M=100, step_size=1.0, n_steps=10000, n_burn=1000)

Run a Path Integral Monte Carlo simulation

Parameters:

• model – Quantum model with log_prob and local_energy methods.

• beta (float) – Inverse temperature.

• M (int) – Number of beads on the integration path.

• step_size (float) – Step size for Metropolis sampling.

• n_steps (int) – Number of Metropolis steps.

• n_burn (int) – Number of initial samples to discard.

Returns:

Sampled paths and estimated energies.

Return type:

tuple

run_vmc(model, step_size=1.0, n_steps=10000, n_burn=1000)

Run a Variational Monte Carlo simulation

Parameters:

• model – Quantum model with log_prob and local_energy methods.

• step_size (float) – Step size for Metropolis sampling.

• n_steps (int) – Number of Metropolis steps.

• n_burn (int) – Number of initial samples to discard.

Returns:

Sampled points and estimated energies.

Return type:

tuple

class toy_qm(alpha=0.5)

Bases: object

A trivial quantum mechanics example using a Gaussian wavefunction

action(x, beta=1.0, M=100)

Action S(x) = int m/2 * dot(x)^2 + V(x)

= sum m/(2*dt) * (x(t+dt) - x(t))^2 + V(x) dt

Parameters:

• x (numpy.ndarray) – Discretized points on the integration path.

• beta (float) – Inverse temperature.

• M (int) – Number of beads (Trotter slices).

kinetic_energy(x)

-1/2 * (d^2 psi / dx^2) / psi

local_energy(x)

log_prob(x)

log_psi(x)

potential_energy(x)

psi(x)

lumeq.opt.mrsf_roks module

MO_BASE = 1

Multi-reference spin-flip (MRSF) DFT uses higher spin reference state and calculates the real ground-state and excited-states with spin flipping. references from Cheol Ho Choi JCP 2019. 10.1063/1.5086895 JCTC 2021. 10.1021/acs.jctc.0c01074 JPCA 2024. 10.1021/acs.jpca.4c04521

MRSF_CIS

alias of MRSF_TDA

class MRSF_TDA(mf, frozen=None)

Bases: TDA

analyze(verbose=None)

gen_vind(mf=None)

Generate function to compute Ax

init_guess(mf, nstates=None, wfnsym=None, return_symmetry=False)

kernel(x0=None, nstates=None)

TDA diagonalization solver

singlet = True

to_gpu(out=None)

Convert a method to its corresponding GPU variant, and recursively converts all attributes of a method to cupy objects or gpu4pyscf objects.

analyze(tdobj, verbose=None)

gen_tda_hop(mf, fock_ao=None, singlet=True, wfnsym=None)

A x

Kwargs:

wfnsymint or str

Point group symmetry irrep symbol or ID for excited CIS wavefunction.

gen_tda_operation(mf, fock_ao=None, singlet=True, wfnsym=None)

A x

Kwargs:

wfnsymint or str

Point group symmetry irrep symbol or ID for excited CIS wavefunction.

mrsf_dimension_transform(nocc, nvir, nbas, singlet=True)

Build the MRSF amplitude-dimension transform matrices.

The returned matrices reorder and combine spin-flip excitation blocks for the singlet or triplet MRSF working equations.

Parameters:

• nocc – Numbers of alpha and beta occupied orbitals.

• nvir – Numbers of alpha and beta virtual orbitals.

• nbas (int) – Number of AO basis functions.

• singlet (bool) – Whether to calculate the singlet (True) or triplet (False).

Returns:

Transform matrix for singlet and triplet excited states.

The rows collect the relevant excitation blocks for the bra side.

Ut (numpy.ndarray): Reordered transform matrix for the ket side.

This matrix is the reordered partner used on the ket side.

Return type:

U (numpy.ndarray)

Notes

Reference: JCP 2019, 10.1063/1.5086895.

The spin-pairing structure of the A matrix follows Fig. 1:

[U S  ] [alpha->beta,  alpha->beta ] [U S]
[UCO  ] [alpha->alpha, alpha->alpha] [UOV]
[UOV  ] [beta ->beta,  beta ->beta ] [UCO]
[UCOCO] [alpha->beta,  beta ->alpha] [UCOCO]
[UOVOV] [beta ->alpha, alpha->beta ] [UOVOV]

spin_square(tdobj)

Calculate <S^2> expectation value for each excited state. S^2 = Sz^2 + 0.5 * (S+ S- + S- S+)

lumeq.opt.mrsf_uks module

MRSF_CIS

alias of MRSF_TDA

class MRSF_TDA(mf, frozen=None)

Bases: TDA

analyze(verbose=None)

gen_vind(mf=None)

Generate function to compute Ax

init_guess(mf, nstates=None, wfnsym=None, return_symmetry=False)

kernel(x0=None, nstates=None)

TDA diagonalization solver

positive_eig_threshold = -0.3

singlet = True

to_gpu(out=None)

Convert a method to its corresponding GPU variant, and recursively converts all attributes of a method to cupy objects or gpu4pyscf objects.

analyze(tdobj, verbose=None)

convert_roks_to_uks(mf)

Convert ROKS object to UKS object for TDDFT calculation.

gen_tda_hop(mf, fock_ao=None, singlet=True, wfnsym=None)

A x

Kwargs:

wfnsymint or str

Point group symmetry irrep symbol or ID for excited CIS wavefunction.

gen_tda_operation(mf, fock_ao=None, singlet=True, wfnsym=None)

A x

Kwargs:

wfnsymint or str

Point group symmetry irrep symbol or ID for excited CIS wavefunction.

lumeq.opt.optimization module

cg_fletcher_reeves(transport, preconditioner, x1, g0, g1)

cg_polak_ribiere(transport, preconditioner, x1, g0, g1)

choose_direction(method='fletcher_reeves')

conjugate_gradient(func, gradf, retraction, transport, preconditioner, x0, ls_method='armijo', cg_method='fletcher_reeves', nmax=50, thresh=1e-08, *args, **kwargs)

gradient_descent(func, gradf, retraction, x0, ls_method='armijo', nmax=50, thresh=1e-08, *args, **kwargs)

line_search(method='armijo')

method options: float, newton, backtracking, armijo, steepest

ls_armijo(func, gradf, retraction, x0, tau=0.5, r=0.0001, alpha=1.0)

alpha is the step size we need to determine reduce it from a large value

ls_backtracking(func, gradf, retraction, x0, tau=0.5, r=0.0001, alpha=1.0)

alpha is the step size we need to determine reduce it from a large value

ls_newton_raphson(func, gradf, retraction, x0, **kwargs)

ls_steepest(step=1.0)

newton_2nd(func, gradf, hessf, tangent_solver, geodesic, x0, dt=1.0, nmax=50, thresh=1e-08, *args, **kwargs)

trust_region(func, gradf, hessf, tangent_solver, geodesic, x0, dt=1.0, nmax=50, thresh=1e-08, *args, **kwargs)

lumeq.opt.riemannian module

Grassmann(**kwargs)

class Grassmann_Quotient(ndim=None, x0=None, A=None, B=None, retraction='qr', **kwargs)

Bases: Stiefel

Grassmann manifold represented as a quotient of the Stiefel manifold.

This implementation works with the orbital coefficient matrix p rather than the projector built from it.

Notes

The quotient-space relations used here are

Gr(k, n) = St(k, n) / O(k)
         = {p in F^(n x k) | p^dagger B p = I_k, P B P = P}

P = p p^dagger
T_p Gr(k, n) = {V in F^(n x k) | p^dagger B V = 0}

check_sanity()

check_tangent(x, v)

verify the obtained vector v is truly a tangent

property dimension

property dist

calculate the distance between two tangent vectors at one point

exp(x, v, dt=1.0)

exp_p (V) = p V cos(s) V^dagger + U sin(s) V^dagger qr is needed for numerically stablity

get_tangent(x, v)

horizontal_lift(x, v, lift=True)

Delta = (0 & - B^T \ B & 0) Delta_p^horizontal = Delta * p = p_perp * B in Stiefel representation

inverse_retraction_polar(x1, x2, dt=1.0)

log(x1, x2, dt=1.0)

adopted from alg. 1 from 10.1007/s10444-023-10090-8 recovers Stiefel from exp operation and avoids matrix inverse

newton_2nd(gradf, hessf, x0, dt=1.0, nmax=50, thresh=1e-08)

projection(x, v)

the resulted tangent vector is perpendicular to BX <Delta, BX> = 0

retraction_polar(x, v, dt=1.0)

Retr_p V = U V^dagger where UsV^dagger = svd(p+V)

retraction_qr(x, v, dt=1.0)

Retr_p V = QD where QR = qr(p+V) and D=diag(sgn(diag(R)+.5))

tangent_solver(x, grad, hess, method='direct')

transport(x0, vg, v0, dt=1.0)

parallel transport v0 at x0 along geodesic defined by vg to tangent space of x1 refer to 10.1137/S0895479895290954

weingarten(x, v, grad, normal)

connection between tangent vectors

class OrthogonalGroup(ndim=None, x0=None, A=None, B=None, retraction='qr', **kwargs)

Bases: Riemannian

the points are n*n square matrix, n-dimensions and n-planes

check_sanity()

property eigenvalue

exp(x, v, dt=1.0)

geodesic mapping tangent vector v of x to other points on manifold exp_p: T_p M -> M

func(x)

assume x is normalized

func_grad(x)

assume x is normalized

log(x1, x2, dt=1.0)

geodesic mapping points on manifold to the tangent space of the first point log_p: M -> T_p M

projection(x, v)

project a general matrix v to the tangent space of a point x on the manifold

retraction_cayley(x, v, dt=1.0)

retraction_norm(x, v, dt=1.0)

put (x+v) on the sphere

retraction_polar(x, v, dt=1.0)

retraction_qr(x, v, dt=1.0)

class Riemannian(ndim=None, x0=None, A=None, B=None, retraction='qr', **kwargs)

Bases: object

check_tangent(x, v)

verify the obtained vector v is truly a tangent

conjugate_gradient(method=None, cg_method='fletcher_reeves', nmax=50, thresh=1e-08, *args, **kwargs)

property dimension

property dist

calculate the distance between two tangent vectors at one point

exp(x, v, dt=1.0)

geodesic mapping tangent vector v of x to other points on manifold exp_p: T_p M -> M

geodesic(x, v, dt=1.0)

geodesic mapping tangent vector v of x to other points on manifold exp_p: T_p M -> M

gradient_descent(method=None, nmax=50, thresh=1e-08, *args, **kwargs)

inverse_retraction()

approximation to logrithrim geodesic mapping project points to get the connecting tangent vector at first point

log(x1, x2, dt=1.0)

geodesic mapping points on manifold to the tangent space of the first point log_p: M -> T_p M

property norm

calculate the length of a tangent vector v of point x use euclidean norm by default

projection(x, v)

project a general matrix v to the tangent space of a point x on the manifold

retraction(x, v)

approximation to exponential geodesic mapping project new vector x+v from point x to the manifold

riemannian_gradient(x, grad)

project euclidean gradient to riemannian gradient of a function

riemannian_hessian(x, v, grad, hess, normal=None)

project euclidean hessian to riemannian hessian of a function

to_tangent_space(x, v)

project a general matrix v to the tangent space of a point x on the manifold

transport(x0, x1, v0)

parallel or vector transport moves tangent vector v0 at point x0 to point x1 T_{x0->x1} (v0): v0 in T_{x0} M -> v1 in T_{x1} M = Proj_M (x1, v0) here x0 is a dummy variable project v0 to the tangent space of x1

weingarten(x, v, *args, **kwargs)

connection between tangent vectors

class Stiefel(ndim=None, x0=None, A=None, B=None, retraction='qr', **kwargs)

Bases: OrthogonalGroup

St(k,n) = {p in F^{n*k} | p^dagger B p = I_k} the point p is a n*k matrix, n-dimensions and k-planes quotient space from OrthogonalGroup where k<n whose tangent vector V is on the tangent space of point p T_p St(k,n) = {V in F^{n*k} | p^dagger B V + V^dagger B p = 0_k}

check_sanity()

check_tangent(x, v)

verify the obtained vector v is truly a tangent

property dimension

property dist

calculate the distance between two tangent vectors at one point

exp(x, v, dt=1.0)

exp_p V = (p \ V) exp((p^dagger V & - V^dagger V \ I_n & p^dagger V)) (exp(-p^dagger V) \ 0_n)

inverse_retraction_polar(x1, x2, dt=1.0)

inverse_retraction_qr(x1, x2, dt=1.0)

log(x1, x2, dt=1.0)

geodesic mapping points on manifold to the tangent space of the first point log_p: M -> T_p M

projection(x, v)

Proj(p,V) = V - p Sym(p^dagger B V)

retraction_polar(x, v, dt=1.0)

Retr_p V = U V^dagger where UsV^dagger = svd(p+V)

retraction_qr(x, v, dt=1.0)

Retr_p V = QD where QR = qr(p+V) and D=diag(sgn(diag(R)+.5))

weingarten(x, v, grad, normal)

connection between tangent vectors

get_random_matrix(ndim, seed=None, sym=True)

solve_sylvester(A, C, B=None, solver='iter')

solve X from AX + XB = C aka. Lyapunov, Stein equation return X

Module contents