Frozen-Core Spin-Component-Scaled RMP2 Analytic Nuclear GradientChemistry
Task prompt
- Gemini 3.1 Pro Preview12.5%of 8 runs
- GPT-5.6 (sol)75%of 8 runs
Frozen-Core Spin-Component-Scaled RMP2 Analytic Nuclear Gradient
Objective
Implement compute_scs_mp2_gradient, an end-to-end driver that returns the analytic nuclear gradient of the frozen-core spin-component-scaled closed-shell MP2 energy, together with a diagnostics dict.
Context
The input Fock matrix fock_mo is in an MO basis but is not assumed diagonal — it may be mixed by an arbitrary orthogonal rotation within the active-occupied block and within the virtual block. The frozen-core block (first n_frozen orbitals) is never mixed with the active-occupied block.
The SCS energy is
and the gradient must be the exact analytic derivative of that scalar with respect to nuclear coordinates. Amplitudes must never touch frozen-core indices, and the Z-vector is solved only over the active-occupied × virtual subspace. The generalized Lagrangian is defined on the full MO space and in general has a nonzero frozen ↔ active-occupied coupling block.
The total gradient must be invariant under an orthogonal rotation of the input MO basis within the active-occupied block or within the virtual block (internal semi-canonicalization is cancelled by tensor-transformation invariance of the contractions).
Conventions
mo_eriandderivs_teiare in chemist notation(pq|rs).derivs_oei['S' | 'T' | 'V']each have shape(Natoms, 3, N, N).derivs_teihas shape(Natoms, 3, N, N, N, N).- The overlap contribution is
sum_{pq} I_{pq} S^{(x)}_{pq}, the one-electron contribution issum_{pq} P_{pq} (T+V)^{(x)}_{pq}, and the two-electron contribution issum_{pqrs} Γ_{pqrs} (pq|rs)^{(x)}.
Diagnostics
Return a dict with
'E_os','E_ss'— unscaled opposite-spin / same-spin MP2 correlation energies,'E_corr_scs' = c_os * E_os + c_ss * E_ss,'z_residual'— Euclidean norm ofG·Z − Xin the semi-canonical basis.
def compute_scs_mp2_gradient(grad_nuc: np.ndarray,
derivs_oei: dict,
derivs_tei: np.ndarray,
mo_eri: np.ndarray,
fock_mo: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
"""
# Frozen-Core Spin-Component-Scaled RMP2 Analytic Nuclear Gradient ## Objective Implement `compute_scs_mp2_gradient`, an end-to-end driver that returns the analytic nuclear gradient of the **frozen-core spin-component-s
Returns:
The computed result described by the prompt.
"""
return grad, infoBackground
Why each piece exists
Frozen core. Core orbitals are chemically inert; freezing them removes cheap-but-numerically-noisy correlation from very deep orbitals and slashes cost. Analytically, the correlation wave-function operator is restricted so that no excitation touches a frozen index; consequently, correlation corrections to the one- and two-particle density matrices vanish on any block involving a frozen index, and the CPHF response is solved only over active-occupied × virtual rotations. However, the generalized Lagrangian still has a nontrivial frozen ↔ active-occupied block, because the Brillouin condition between frozen and active occupied orbitals is broken once a frozen-core/SCS correlation energy is added.
Spin-component scaling. Grimme's SCS-MP2 rescales the same-spin and opposite-spin pair-correlation contributions by separate empirical factors. For the gradient to remain the exact analytic derivative of the SCS energy, the rescaling must be applied at the amplitude/density level (not only at the energy level), so that the unrelaxed densities, the Z-vector right-hand side, and the relaxed TPDM all carry consistent SCS weights.
Semi-canonicalization. MP2 amplitudes are defined with orbital-energy denominators ε_i + ε_j − ε_a − ε_b. When the input Fock matrix is not diagonal within the active-occupied and virtual sub-blocks, we diagonalize those blocks (without mixing them with the frozen-core block), transform all MO-basis tensors into that semi-canonical basis, do the MP2 machinery there, and rely on the fact that every physical contraction in the gradient is a full tensor trace — it is invariant under any orthogonal rotation applied consistently to all its tensor legs.
Z-vector. The unrelaxed MP2 density is not variational with respect to orbital rotations, so a CPHF-like linear system must be solved to obtain the orbital relaxation contribution. Under frozen core, this linear system lives only on the active-occupied × virtual subspace.
Assembly. The final gradient is the standard sum of four contractions: nuclear repulsion, overlap with the generalized Lagrangian, one-electron derivatives with the relaxed OPDM, and two-electron derivatives with the relaxed TPDM.
Parameters and returns
# Frozen-Core Spin-Component-Scaled RMP2 Analytic Nuclear Gradient ## Objective Implement compute_scs_mp2_gradient, an end-to-end driver that returns the analytic nuclear gradient of the **frozen-core spin-component-s
Returns:
The computed result described by the prompt.
Required dependencies
import numpy as np # numpy==2.4.4
import copyPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
import numpy as np
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def test_case_1():
"""Shapes, diagnostics keys, SCS energy relation, and n_frozen=0, (1,1) base case."""
grad_nuc, doei, dtei, eri, F, n_occ, n_fr = _make_random_system(
N=6, n_occ=3, n_frozen=0, seed=1, natoms=2)
grad, info = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, n_fr, c_os=1.0, c_ss=1.0)
assert grad.shape == (2, 3)
for k in ('E_os', 'E_ss', 'E_corr_scs', 'z_residual'):
assert k in info
assert isinstance(info[k], float)
assert np.isclose(info['E_corr_scs'], info['E_os'] + info['E_ss'], atol=1e-10)
assert info['z_residual'] < 1e-6
# SCS weighting relation
grad2, info2 = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, n_fr, c_os=1.2, c_ss=0.3333)
assert np.isclose(info2['E_corr_scs'],
1.2 * info2['E_os'] + 0.3333 * info2['E_ss'], atol=1e-10)
# Same unscaled components
assert np.isclose(info['E_os'], info2['E_os'], atol=1e-10)
assert np.isclose(info['E_ss'], info2['E_ss'], atol=1e-10)
import numpy as npTest case 2 of 3
import numpy as np
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def test_case_2():
"""Invariance under an orthogonal rotation of the active-occupied and virtual blocks of the input MO basis."""
grad_nuc, doei, dtei, eri, F, n_occ, n_fr = _make_random_system(
N=6, n_occ=3, n_frozen=1, seed=2, natoms=2)
grad_a, info_a = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, n_fr, c_os=1.2, c_ss=0.33)
# Build an orthogonal matrix that mixes active-occupied (n_occ-n_fr = 2) and virtual (3) blocks
rng = np.random.default_rng(42)
N = F.shape[0]
U = np.eye(N)
# Active-occupied rotation
Ao = rng.standard_normal((n_occ - n_fr, n_occ - n_fr))
Qo, _ = np.linalg.qr(Ao)
U[n_fr:n_occ, n_fr:n_occ] = Qo
# Virtual rotation
Av = rng.standard_normal((N - n_occ, N - n_occ))
Qv, _ = np.linalg.qr(Av)
U[n_occ:, n_occ:] = Qv
F_rot = U.T @ F @ U
eri_rot = np.einsum('pqrs,pP,qQ,rR,sS->PQRS', eri, U, U, U, U, optimize=True)
dS_rot = np.einsum('Axpq,pP,qQ->AxPQ', doei['S'], U, U, optimize=True)
dT_rot = np.einsum('Axpq,pP,qQ->AxPQ', doei['T'], U, U, optimize=True)
dV_rot = np.einsum('Axpq,pP,qQ->AxPQ', doei['V'], U, U, optimize=True)
doei_rot = {'S': dS_rot, 'T': dT_rot, 'V': dV_rot}
dtei_rot = np.einsum('Axpqrs,pP,qQ,rR,sS->AxPQRS',
dtei, U, U, U, U, optimize=True)
grad_b, info_b = compute_scs_mp2_gradient(grad_nuc, doei_rot, dtei_rot,
eri_rot, F_rot,
n_occ, n_fr, c_os=1.2, c_ss=0.33)
assert np.allclose(grad_a, grad_b, atol=1e-8), \
f"Gradient not invariant under active/virtual rotation; diff={np.max(np.abs(grad_a-grad_b))}"
assert np.isclose(info_a['E_corr_scs'], info_b['E_corr_scs'], atol=1e-10)
import numpy as npTest case 3 of 3
import numpy as np
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def _make_random_system(N=6, n_occ=3, n_frozen=0, seed=0, natoms=2):
rng = np.random.default_rng(seed)
# random symmetric Fock with realistic occupied/virtual ordering
eps = np.sort(rng.uniform(-2.0, -0.5, size=n_occ))
eps_v = np.sort(rng.uniform(0.3, 2.0, size=N - n_occ))
diag = np.concatenate([eps, eps_v])
F = np.diag(diag)
# ERI: build as chemist-symmetric random tensor
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.125 * (A + A.transpose(1, 0, 2, 3) + A.transpose(0, 1, 3, 2)
+ A.transpose(1, 0, 3, 2) + A.transpose(2, 3, 0, 1)
+ A.transpose(3, 2, 0, 1) + A.transpose(2, 3, 1, 0)
+ A.transpose(3, 2, 1, 0))
# Derivatives
dS = rng.standard_normal((natoms, 3, N, N))
dS = 0.5 * (dS + dS.transpose(0, 1, 3, 2))
dT = rng.standard_normal((natoms, 3, N, N))
dT = 0.5 * (dT + dT.transpose(0, 1, 3, 2))
dV = rng.standard_normal((natoms, 3, N, N))
dV = 0.5 * (dV + dV.transpose(0, 1, 3, 2))
derivs_oei = {'S': dS, 'T': dT, 'V': dV}
B = rng.standard_normal((natoms, 3, N, N, N, N)) * 0.05
derivs_tei = 0.25 * (B + B.transpose(0, 1, 3, 2, 4, 5)
+ B.transpose(0, 1, 2, 3, 5, 4)
+ B.transpose(0, 1, 4, 5, 2, 3))
grad_nuc = rng.standard_normal((natoms, 3))
return grad_nuc, derivs_oei, derivs_tei, eri, F, n_occ, n_frozen
def test_case_3():
"""Frozen-core sanity: freezing a deep core should change the gradient but keep diagnostics consistent."""
grad_nuc, doei, dtei, eri, F, n_occ, _ = _make_random_system(
N=6, n_occ=3, n_frozen=0, seed=7, natoms=2)
# Full (no frozen)
g0, i0 = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, 0, c_os=1.0, c_ss=1.0)
# 1 frozen core
g1, i1 = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, 1, c_os=1.0, c_ss=1.0)
# Correlation energy with frozen core must be smaller in magnitude than full
assert abs(i1['E_corr_scs']) <= abs(i0['E_corr_scs']) + 1e-10
# Z-vector residuals well converged
assert i0['z_residual'] < 1e-6
assert i1['z_residual'] < 1e-6
# And gradients differ (frozen core really changes the physics)
assert not np.allclose(g0, g1, atol=1e-4)
# SCS-weighted energy matches formula regardless of n_frozen
_, i2 = compute_scs_mp2_gradient(grad_nuc, doei, dtei, eri, F,
n_occ, 1, c_os=1.3, c_ss=0.4)
assert np.isclose(i2['E_corr_scs'], 1.3 * i2['E_os'] + 0.4 * i2['E_ss'],
atol=1e-10)Solution
import numpy as np
import copy
def semicanonicalize_active_blocks(fock_mo: np.ndarray,
mo_eri: np.ndarray,
derivs_oei: dict,
derivs_tei: np.ndarray,
n_occ: int,
n_frozen: int) -> tuple:
N = fock_mo.shape[0]
n_act = n_occ - n_frozen
n_vir = N - n_occ
U = np.eye(N)
eps = np.zeros(N)
if n_frozen > 0:
eps[:n_frozen] = np.diag(fock_mo)[:n_frozen]
if n_act > 0:
F_oo = fock_mo[n_frozen:n_occ, n_frozen:n_occ]
F_oo = 0.5 * (F_oo + F_oo.T)
w_o, V_o = np.linalg.eigh(F_oo)
U[n_frozen:n_occ, n_frozen:n_occ] = V_o
eps[n_frozen:n_occ] = w_o
if n_vir > 0:
F_vv = fock_mo[n_occ:, n_occ:]
F_vv = 0.5 * (F_vv + F_vv.T)
w_v, V_v = np.linalg.eigh(F_vv)
U[n_occ:, n_occ:] = V_v
eps[n_occ:] = w_v
fock_sc = U.T @ fock_mo @ U
eri_sc = np.einsum('pqrs,pP,qQ,rR,sS->PQRS', mo_eri, U, U, U, U, optimize=True)
dS = np.einsum('Axpq,pP,qQ->AxPQ', derivs_oei['S'], U, U, optimize=True)
dT = np.einsum('Axpq,pP,qQ->AxPQ', derivs_oei['T'], U, U, optimize=True)
dV = np.einsum('Axpq,pP,qQ->AxPQ', derivs_oei['V'], U, U, optimize=True)
dERI = np.einsum('Axpqrs,pP,qQ,rR,sS->AxPQRS', derivs_tei, U, U, U, U, optimize=True)
derivs_oei_sc = {'S': dS, 'T': dT, 'V': dV}
return eps, fock_sc, eri_sc, derivs_oei_sc, dERI, U
def compute_scs_frozen_core_unrelaxed_densities(mo_eri: np.ndarray,
epsilon: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
N = mo_eri.shape[0]
n_act = n_occ - n_frozen
n_vir = N - n_occ
o_slice = slice(n_frozen, n_occ)
v_slice = slice(n_occ, N)
eri_oovv_phys = mo_eri[o_slice, v_slice, o_slice, v_slice].transpose(0, 2, 1, 3)
e_o = epsilon[n_frozen:n_occ]
e_v = epsilon[n_occ:]
D = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
if n_act > 0 and n_vir > 0:
t2 = eri_oovv_phys / D
else:
t2 = np.zeros((n_act, n_act, n_vir, n_vir))
if n_act > 0 and n_vir > 0:
E_os = float(np.einsum('ijab,ijab->', eri_oovv_phys, t2, optimize=True))
E_ss = float(np.einsum('ijab,ijab->', eri_oovv_phys,
t2 - t2.swapaxes(2, 3), optimize=True))
else:
E_os = 0.0
E_ss = 0.0
E_corr_scs = c_os * E_os + c_ss * E_ss
t2_swap = t2.swapaxes(2, 3)
t2_tilde_scs = c_os * t2 + c_ss * (t2 - t2_swap)
opdm = np.zeros((N, N))
for i in range(n_occ):
opdm[i, i] = 2.0
if n_act > 0 and n_vir > 0:
Pij = -np.einsum('ikab,jkab->ij', t2, t2_tilde_scs, optimize=True)
Pij += -np.einsum('jkab,ikab->ij', t2, t2_tilde_scs, optimize=True)
Pab = np.einsum('ijac,ijbc->ab', t2, t2_tilde_scs, optimize=True)
Pab += np.einsum('ijbc,ijac->ab', t2, t2_tilde_scs, optimize=True)
opdm[o_slice, o_slice] += Pij
opdm[v_slice, v_slice] += Pab
opdm = 0.5 * (opdm + opdm.T)
opdm_ref = np.zeros((N, N))
for i in range(n_occ):
opdm_ref[i, i] = 2.0
tpdm = -0.25 * (2.0 * np.einsum('pr,qs->pqrs', opdm_ref, opdm_ref, optimize=True)
- np.einsum('ps,qr->pqrs', opdm_ref, opdm_ref, optimize=True))
if n_act > 0 and n_vir > 0:
tpdm[o_slice, o_slice, v_slice, v_slice] += t2_tilde_scs
tpdm[v_slice, v_slice, o_slice, o_slice] += t2_tilde_scs.transpose(2, 3, 0, 1)
return opdm, tpdm, E_os, E_ss, E_corr_scs
def solve_frozen_core_zvector_and_lagrangian(opdm: np.ndarray,
tpdm: np.ndarray,
mo_eri: np.ndarray,
fock_mo: np.ndarray,
epsilon: np.ndarray,
n_occ: int,
n_frozen: int) -> tuple:
N = opdm.shape[0]
nf = n_frozen
no = n_occ
nao = no - nf
nv = N - no
o_act = slice(nf, no)
v = slice(no, N)
f = slice(0, nf)
np_eri = mo_eri.swapaxes(1, 2)
Ip = np.zeros((N, N))
Ip += np.einsum('pr,rq->pq', fock_mo, opdm, optimize=True)
Ip += np.einsum('qr,rp->pq', opdm, fock_mo, optimize=True)
Ip += np.einsum('qrst,prst->pq', tpdm, np_eri, optimize=True)
Ip += np.einsum('rqst,rpst->pq', tpdm, np_eri, optimize=True)
Ip += np.einsum('rsqt,rspt->pq', tpdm, np_eri, optimize=True)
Ip += np.einsum('rstq,rstp->pq', tpdm, np_eri, optimize=True)
Ip[:, :n_occ] += 4.0 * np.einsum('rs,rpsq->pq', opdm,
np_eri[:, :, :, :n_occ], optimize=True)
Ip[:, :n_occ] -= 1.0 * np.einsum('rs,rpqs->pq', opdm,
np_eri[:, :, :n_occ, :], optimize=True)
Ip[:, :n_occ] -= 1.0 * np.einsum('rs,rqps->pq', opdm,
np_eri[:, :n_occ, :, :], optimize=True)
Ip *= -0.5
o_slice = slice(n_frozen, n_occ)
v_slice = slice(n_occ, N)
X = Ip[o_slice, v_slice].T - Ip[v_slice, o_slice]
e_o = epsilon[n_frozen:n_occ]
e_v = epsilon[n_occ:]
if nao > 0 and nv > 0:
G = 4.0 * np_eri[o_slice, o_slice, v_slice, v_slice]
G -= np_eri[o_slice, o_slice, v_slice, v_slice].swapaxes(2, 3)
G -= np_eri[o_slice, v_slice, o_slice, v_slice].swapaxes(1, 2)
G = G.swapaxes(1, 2)
eps_diag = e_v.reshape(-1, 1) - e_o
I_o = np.eye(nao)
I_v = np.eye(nv)
G += np.einsum('ai,ij,ab->iajb', eps_diag, I_o, I_v, optimize=True)
dim = nao * nv
G_flat = G.T.reshape(dim, dim)
X_flat = X.reshape(dim, 1)
Z_flat = np.linalg.solve(G_flat, X_flat)
Z = Z_flat.reshape(nv, nao)
z_residual = float(np.linalg.norm(G_flat @ Z_flat - X_flat))
else:
Z = np.zeros((nv, nao))
z_residual = 0.0
relaxed_opdm = opdm.copy()
if nao > 0 and nv > 0:
relaxed_opdm[o_slice, v_slice] = -Z.T
relaxed_opdm[v_slice, o_slice] = -Z
I_lag = copy.deepcopy(Ip)
I_lag[v_slice, o_slice] = Ip[o_slice, v_slice].T
if nao > 0 and nv > 0:
ERI_vooo = np_eri[v_slice, o_slice, o_slice, o_slice]
I_lag[o_slice, o_slice] += 2.0 * np.einsum('ak,aikj->ij', Z, ERI_vooo, optimize=True)
I_lag[o_slice, o_slice] -= 1.0 * np.einsum('ak,aijk->ij', Z, ERI_vooo, optimize=True)
I_lag[o_slice, o_slice] += 2.0 * np.einsum('ak,ajki->ij', Z, ERI_vooo, optimize=True)
I_lag[o_slice, o_slice] -= 1.0 * np.einsum('ak,ajik->ij', Z, ERI_vooo, optimize=True)
f_act_diag = np.diag(fock_mo)[n_frozen:n_occ]
I_lag[v_slice, o_slice] += Z * f_act_diag
I_lag[o_slice, v_slice] += (Z * f_act_diag).T
if nf > 0:
E1 = np_eri[:nf, v_slice, o_slice, o_slice]
I_lag[:nf, o_slice] += 2.0 * np.einsum('ak,fakj->fj', Z, E1, optimize=True)
I_lag[:nf, o_slice] -= 1.0 * np.einsum('ak,fajk->fj', Z, E1, optimize=True)
E2 = np_eri[:nf, o_slice, o_slice, v_slice]
I_lag[:nf, o_slice] += 2.0 * np.einsum('ak,fjka->fj', Z, E2, optimize=True)
I_lag[:nf, o_slice] -= 1.0 * np.einsum('ak,fjak->fj', Z, np_eri[:nf, o_slice, v_slice, o_slice], optimize=True)
I_lag[o_slice, :nf] = I_lag[:nf, o_slice].T
return relaxed_opdm, I_lag, Z, z_residual
def compute_scs_mp2_gradient(grad_nuc: np.ndarray,
derivs_oei: dict,
derivs_tei: np.ndarray,
mo_eri: np.ndarray,
fock_mo: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
if not (0 <= n_frozen <= n_occ):
raise ValueError("n_frozen must be between 0 and n_occ")
N = fock_mo.shape[0]
# Call Subproblem 1: Semi-canonicalize
epsilon, fock_sc, eri_sc, derivs_oei_sc, dERI, U = semicanonicalize_active_blocks(
fock_mo, mo_eri, derivs_oei, derivs_tei, n_occ, n_frozen
)
# Call Subproblem 2: Build unrelaxed densities
opdm, tpdm, E_os, E_ss, E_corr_scs = compute_scs_frozen_core_unrelaxed_densities(
eri_sc, epsilon, n_occ, n_frozen, c_os, c_ss
)
# Call Subproblem 3: Solve Z-vector and Lagrangian
relaxed_opdm, I_lag, Z, z_residual = solve_frozen_core_zvector_and_lagrangian(
opdm, tpdm, eri_sc, fock_sc, epsilon, n_occ, n_frozen
)
# Relax TPDM
relaxed_tpdm = tpdm.copy()
eye_o = np.eye(n_occ)
relaxed_tpdm[:, :n_occ, :, :n_occ] += 2.0 * np.einsum('pr,qs->pqrs',
relaxed_opdm, eye_o,
optimize=True)
relaxed_tpdm[:, :n_occ, :n_occ, :] -= 1.0 * np.einsum('ps,qr->pqrs',
relaxed_opdm, eye_o,
optimize=True)
# Assemble gradient
grad = np.array(grad_nuc, dtype=float, copy=True)
grad += np.einsum('pq,Axpq->Ax', I_lag, derivs_oei_sc['S'], optimize=True)
grad += np.einsum('pq,Axpq->Ax', relaxed_opdm, derivs_oei_sc['T'], optimize=True)
grad += np.einsum('pq,Axpq->Ax', relaxed_opdm, derivs_oei_sc['V'], optimize=True)
grad += np.einsum('pqrs,Axprqs->Ax', relaxed_tpdm, dERI, optimize=True)
info = {
'E_os': float(E_os),
'E_ss': float(E_ss),
'E_corr_scs': float(E_corr_scs),
'z_residual': float(z_residual),
}
return grad, infoPrompt
Objective
Implement semicanonicalize_active_blocks that produces a block-diagonal unitary rotation which diagonalizes the active-occupied and virtual sub-blocks of a (possibly non-diagonal) Fock matrix, while leaving the frozen-core block untouched.
This semi-canonicalization is a prerequisite for defining MP2 amplitude energy denominators when the input MO Fock matrix is not already diagonal, as can happen after freezing core orbitals, under localization, or in non-canonical reference treatments.
Function signature
def semicanonicalize_active_blocks(fock_mo: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
"""
**Objective** Implement `semicanonicalize_active_blocks` that produces a block-diagonal unitary rotation which diagonalizes the active-occupied and virtual sub-blocks of a (possibly non-diagonal) Fock matrix, while leavi
Returns:
The computed result described by the prompt.
"""
return U, epsilon_scRequirements
Umust be orthogonal (to tolerance1e-10) and exactly the identity on indices[0, n_frozen).Umust not mix frozen, active-occupied, and virtual subspaces.U.T @ fock_mo @ Umust be diagonal within the active-occupied block (rows/cols[n_frozen, n_occ)) and within the virtual block (rows/cols[n_occ, N)).epsilon_scpacks the frozen-core diagonal, then sorted active-occupied eigenvalues, then sorted virtual eigenvalues.
Background
In post-Hartree-Fock methods such as MP2, the standard working equations assume that orbital energies are well-defined numbers attached to individual orbitals. This is automatic when the Fock matrix is diagonal in the MO basis, which is the canonical Hartree-Fock choice.
However, several common situations break that assumption while leaving the reference determinant itself unchanged. Freezing core orbitals partitions the occupied space into a frozen-core block and an active-occupied block; if the frozen and active occupied orbitals were not separately canonical to begin with, the active-occupied sub-block of the Fock matrix need not be diagonal. Orbital localization, natural-orbital bases, and ROHF-like treatments can leave similar off-diagonal structure. In all of these cases, MP2 denominators like ε_i + ε_j − ε_a − ε_b are ill-defined until each relevant orbital subspace is rotated to its own canonical form.
The standard remedy is semi-canonicalization: within each subspace that participates in the correlation treatment independently (here, the active-occupied block and the virtual block), perform an internal orthogonal rotation that diagonalizes the Fock matrix restricted to that subspace. Crucially, no rotation is allowed between subspaces — for example, mixing a frozen-core orbital with an active-occupied orbital would alter the reference determinant and invalidate the frozen-core approximation.
The resulting overall rotation is an orthogonal matrix that is block-diagonal with respect to the {frozen, active-occupied, virtual} partition. It leaves the Hartree-Fock energy and all physical observables invariant, but gives a basis in which orbital-energy denominators are well-defined so that MP2 amplitudes, Z-vector equations, and related machinery can be written in their usual form.
Function signature
def semicanonicalize_active_blocks(fock_mo: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
return U, epsilon_scPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
import numpy as np
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def test_case_1():
# Basic shapes, orthogonality, block-diagonal structure, and diagonalization.
N, n_occ, n_frozen = 6, 3, 1
F = _random_symmetric(N, seed=7)
# Make occupied diagonal dominant to keep realistic
np.fill_diagonal(F, np.array([-20.0, -5.0, -3.0, 1.0, 2.0, 4.0]))
F = 0.5 * (F + F.T)
U, eps_sc = semicanonicalize_active_blocks(F, n_occ, n_frozen)
assert U.shape == (N, N), "U must have shape (N, N)"
assert eps_sc.shape == (N,), "epsilon_sc must have shape (N,)"
# Orthogonality
assert np.allclose(U.T @ U, np.eye(N), atol=1e-10), "U must be orthogonal"
# Identity on frozen-core rows/cols
assert np.allclose(U[:n_frozen, :], np.eye(N)[:n_frozen, :], atol=1e-12), \
"U must be identity on frozen-core rows"
assert np.allclose(U[:, :n_frozen], np.eye(N)[:, :n_frozen], atol=1e-12), \
"U must be identity on frozen-core cols"
# Block-diagonal across {frozen, active-occ, virtual}
assert np.allclose(U[:n_frozen, n_frozen:], 0.0, atol=1e-12)
assert np.allclose(U[n_frozen:n_occ, :n_frozen], 0.0, atol=1e-12)
assert np.allclose(U[n_frozen:n_occ, n_occ:], 0.0, atol=1e-12)
assert np.allclose(U[n_occ:, :n_occ], 0.0, atol=1e-12)
# Active-occ and virtual sub-blocks of U^T F U must be diagonal
Frot = U.T @ F @ U
act = Frot[n_frozen:n_occ, n_frozen:n_occ]
vir = Frot[n_occ:, n_occ:]
assert np.allclose(act - np.diag(np.diag(act)), 0.0, atol=1e-10), \
"Active-occupied block must be diagonal after rotation"
assert np.allclose(vir - np.diag(np.diag(vir)), 0.0, atol=1e-10), \
"Virtual block must be diagonal after rotation"
# epsilon_sc consistency
assert np.allclose(eps_sc[:n_frozen], np.diag(F)[:n_frozen], atol=1e-12)
assert np.allclose(eps_sc[n_frozen:n_occ], np.diag(act), atol=1e-10)
assert np.allclose(eps_sc[n_occ:], np.diag(vir), atol=1e-10)
# Sorted ascending within each block
assert np.all(np.diff(eps_sc[n_frozen:n_occ]) >= -1e-12)
assert np.all(np.diff(eps_sc[n_occ:]) >= -1e-12)
import numpy as npTest case 2 of 3
import numpy as np
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def test_case_2():
# Already-diagonal Fock: U should be identity, eps_sc should equal diag(F).
N, n_occ, n_frozen = 5, 2, 1
diag_vals = np.array([-10.0, -1.5, 0.7, 2.3, 5.1])
F = np.diag(diag_vals)
U, eps_sc = semicanonicalize_active_blocks(F, n_occ, n_frozen)
assert np.allclose(U, np.eye(N), atol=1e-10), \
"U must be identity when F is already diagonal"
assert np.allclose(eps_sc, diag_vals, atol=1e-12), \
"epsilon_sc must equal diag(F) when F is diagonal"
import numpy as npTest case 3 of 3
import numpy as np
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def _random_symmetric(n, seed):
rng = np.random.default_rng(seed)
A = rng.standard_normal((n, n))
return 0.5 * (A + A.T)
def test_case_3():
# No frozen core and no virtuals edge-cases; plus invariance of eigenvalues.
N, n_occ, n_frozen = 4, 2, 0
F = _random_symmetric(N, seed=13)
U, eps_sc = semicanonicalize_active_blocks(F, n_occ, n_frozen)
# With n_frozen=0, still block-diagonal across occ/vir only
assert np.allclose(U[:n_occ, n_occ:], 0.0, atol=1e-12), \
"No mixing between occupied and virtual blocks"
assert np.allclose(U[n_occ:, :n_occ], 0.0, atol=1e-12)
# Orthogonality
assert np.allclose(U.T @ U, np.eye(N), atol=1e-10)
# Eigenvalues of active-occupied block in F must match eps_sc[:n_occ]
w_occ = np.linalg.eigvalsh(F[:n_occ, :n_occ])
w_vir = np.linalg.eigvalsh(F[n_occ:, n_occ:])
assert np.allclose(np.sort(eps_sc[:n_occ]), np.sort(w_occ), atol=1e-10), \
"Active-occupied semi-canonical energies must be eigenvalues of the occ block"
assert np.allclose(np.sort(eps_sc[n_occ:]), np.sort(w_vir), atol=1e-10), \
"Virtual semi-canonical energies must be eigenvalues of the virtual block"
# Rotated F must be diagonal on both blocks
Frot = U.T @ F @ U
occ_blk = Frot[:n_occ, :n_occ]
vir_blk = Frot[n_occ:, n_occ:]
assert np.allclose(occ_blk - np.diag(np.diag(occ_blk)), 0.0, atol=1e-10)
assert np.allclose(vir_blk - np.diag(np.diag(vir_blk)), 0.0, atol=1e-10)Solution
import numpy as np
def semicanonicalize_active_blocks(fock_mo: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
N = fock_mo.shape[0]
if not (0 <= n_frozen <= n_occ <= N):
raise ValueError("Require 0 <= n_frozen <= n_occ <= N")
U = np.eye(N)
epsilon_sc = np.zeros(N)
# Frozen-core block: identity, diagonal values of fock_mo
if n_frozen > 0:
epsilon_sc[:n_frozen] = np.diag(fock_mo)[:n_frozen]
# Active-occupied block
n_act_occ = n_occ - n_frozen
if n_act_occ > 0:
F_act = fock_mo[n_frozen:n_occ, n_frozen:n_occ]
F_act = 0.5 * (F_act + F_act.T)
evals_act, evecs_act = np.linalg.eigh(F_act)
# eigh returns ascending eigenvalues already
U[n_frozen:n_occ, n_frozen:n_occ] = evecs_act
epsilon_sc[n_frozen:n_occ] = evals_act
# Virtual block
n_vir = N - n_occ
if n_vir > 0:
F_vir = fock_mo[n_occ:, n_occ:]
F_vir = 0.5 * (F_vir + F_vir.T)
evals_vir, evecs_vir = np.linalg.eigh(F_vir)
U[n_occ:, n_occ:] = evecs_vir
epsilon_sc[n_occ:] = evals_vir
return U, epsilon_scObjective
Implement compute_scs_frozen_core_unrelaxed_densities to build the unrelaxed one- and two-particle density matrices (OPDM, TPDM) and the spin-component-resolved correlation energies for a frozen-core, spin-component-scaled RMP2 calculation.
Setup
- The MO basis is assumed semi-canonical: the Fock matrix is diagonal with entries given by
epsilon, andmo_eriis in chemist notation(pq|rs)in that basis. - Orbital partition (indexing
[0, N)): - frozen-core:[0, n_frozen)- active-occupied:[n_frozen, n_occ)- virtual:[n_occ, N) - Same-spin and opposite-spin correlation pieces are weighted by
c_ssandc_osrespectively when forming the amplitudes whose analytic gradient the driver will eventually compute. E_osandE_ssreturned by this function are the unscaled spin components of the ordinary RMP2 correlation energy (so thatc_os*E_os + c_ss*E_ssis the SCS correlation energy).
Requirements
- No amplitude involving any frozen-core orbital may contribute: all frozen-core rows/columns of
opdmandtpdmmust equal their plain HF reference values. - When
n_frozen == 0and(c_os, c_ss) == (1.0, 1.0), the output must reduce to the standard closed-shell unrelaxed RMP2 OPDM/TPDM, andE_os + E_ssmust equal the standard RMP2 correlation energy. - The returned
opdmmust be symmetric within1e-10. - Do not call
semicanonicalize_active_blocksorsolve_frozen_core_zvector_and_lagrangian.
def compute_scs_frozen_core_unrelaxed_densities(mo_eri: np.ndarray,
epsilon: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
"""
### Objective Implement `compute_scs_frozen_core_unrelaxed_densities` to build the unrelaxed one- and two-particle density matrices (OPDM, TPDM) and the spin-component-resolved correlation energies for a frozen-core, sp
Returns:
The computed result described by the prompt.
"""
return opdm, tpdm, E_os, E_ssBackground
In post-Hartree-Fock gradient theory, the energy can be written as a single trace of density matrices against integrals once the response of the orbitals is folded in. The quantities that carry the correlation physics are the unrelaxed one-particle and two-particle density matrices (OPDM and TPDM). The OPDM has a reference Hartree-Fock part (doubly occupied orbitals) plus correlation corrections on the occupied-occupied and virtual-virtual blocks; the TPDM has a reference antisymmetrized product of occupied densities plus correlation pieces built from the MP2 amplitudes. Frozen-core approximations exclude chemically inert core electrons from the correlation treatment: no excitation is allowed out of a frozen orbital. This keeps the core fully described at the Hartree-Fock level and removes core-correlation contributions from both the energy and the density matrices. Spin-component-scaled MP2 (SCS-MP2), due to Grimme, separately rescales the opposite-spin and same-spin pair-correlation energies by empirical factors (commonly 6/5 and 1/3). Writing the analytic gradient of the SCS energy requires propagating those same scale factors into the MP2 amplitudes that enter the density matrices, so that the densities describe the derivative of the SCS energy rather than of plain MP2.
Function signature
def compute_scs_frozen_core_unrelaxed_densities(mo_eri: np.ndarray,
epsilon: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
return opdm, tpdm, E_os, E_ssPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
import numpy as np
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def test_case_1():
# Reduction to plain RMP2 when n_frozen=0 and c_os=c_ss=1
N, n_occ = 5, 2
mo_eri, eps = _random_system(N, n_occ, seed=1)
opdm, tpdm, E_os, E_ss = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0)
# Compute plain RMP2 correlation energy independently
o = slice(0, n_occ); v = slice(n_occ, N)
iajb = mo_eri[o, v, o, v]
e_o = eps[:n_occ]; e_v = eps[n_occ:]
D = (e_o[:, None, None, None] + e_o[None, None, :, None]
- e_v[None, :, None, None] - e_v[None, None, None, :])
t_iajb = iajb / D
t2 = t_iajb.transpose(0, 2, 1, 3)
g = iajb.transpose(0, 2, 1, 3)
E_ref = np.einsum('ijab,ijab->', g, 2 * t2 - t2.transpose(0, 1, 3, 2))
assert np.isclose(E_os + E_ss, E_ref, atol=1e-10), \
f"E_os+E_ss={E_os+E_ss} should equal plain RMP2 E={E_ref}"
# OPDM symmetry
assert np.allclose(opdm, opdm.T, atol=1e-10), "OPDM must be symmetric"
# HF reference trace: doubly occupied -> 2*n_occ + correlation trace
# The HF piece gives trace 2*n_occ; correlation trace should be ~0 to leading order
# but need not be exactly zero. Just check diagonal frozen/occupied baseline > 0.
assert np.all(np.diag(opdm)[:n_occ] > 1.0), "Occupied OPDM diagonal should be near 2"
import numpy as npTest case 2 of 3
import numpy as np
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def test_case_2():
# Frozen-core: no amplitudes touch frozen indices -> frozen rows/cols of opdm/tpdm
# must match plain HF reference.
N, n_occ = 5, 3
n_frozen = 1
mo_eri, eps = _random_system(N, n_occ, seed=2)
opdm, tpdm, E_os, E_ss = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=n_frozen, c_os=1.2, c_ss=0.3)
# Build HF reference densities
D_ref = np.zeros((N, N))
for p in range(n_occ):
D_ref[p, p] = 2.0
tpdm_ref = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Frozen rows/cols of OPDM must equal D_ref
fc = slice(0, n_frozen)
assert np.allclose(opdm[fc, :], D_ref[fc, :], atol=1e-10), \
"Frozen-core rows of OPDM must match HF reference"
assert np.allclose(opdm[:, fc], D_ref[:, fc], atol=1e-10), \
"Frozen-core cols of OPDM must match HF reference"
# Any TPDM element with at least one frozen index must match reference
mask = np.zeros((N, N, N, N), dtype=bool)
mask[fc, :, :, :] = True
mask[:, fc, :, :] = True
mask[:, :, fc, :] = True
mask[:, :, :, fc] = True
assert np.allclose(tpdm[mask], tpdm_ref[mask], atol=1e-10), \
"TPDM frozen-index entries must equal HF reference"
import numpy as npTest case 3 of 3
import numpy as np
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def _random_system(N=5, n_occ=2, seed=0):
rng = np.random.default_rng(seed)
# Random symmetric ERI with chemist 8-fold-ish symmetry enforced partially
A = rng.standard_normal((N, N, N, N)) * 0.1
# symmetrize (pq|rs) = (qp|rs) = (pq|sr) = (rs|pq)
A = 0.5 * (A + A.transpose(1, 0, 2, 3))
A = 0.5 * (A + A.transpose(0, 1, 3, 2))
A = 0.5 * (A + A.transpose(2, 3, 0, 1))
eps = np.array([-1.0 - 0.3 * i for i in range(n_occ)] +
[1.0 + 0.3 * i for i in range(N - n_occ)])
return A, eps
def test_case_3():
# SCS weighting: E_os, E_ss returned are the UNSCALED components and do not
# depend on c_os, c_ss. Also OPDM scales linearly in (c_os, c_ss) in its
# correlation part.
N, n_occ = 5, 2
mo_eri, eps = _random_system(N, n_occ, seed=3)
opdm_a, tpdm_a, Eos_a, Ess_a = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0)
opdm_b, tpdm_b, Eos_b, Ess_b = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=1.2, c_ss=1.0 / 3.0)
# Spin-component energies are unscaled => should be identical
assert np.isclose(Eos_a, Eos_b, atol=1e-12), "E_os must not depend on c_os/c_ss"
assert np.isclose(Ess_a, Ess_b, atol=1e-12), "E_ss must not depend on c_os/c_ss"
# Reference HF OPDM
D_ref = np.zeros((N, N))
for p in range(n_occ):
D_ref[p, p] = 2.0
# Correlation part of OPDM = opdm - D_ref; check linearity in (c_os, c_ss).
corr_a = opdm_a - D_ref
corr_b = opdm_b - D_ref
# Rerun with c_os only and c_ss only to decompose
_, _, _, _ = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=0.0)
opdm_os, _, _, _ = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=0.0)
opdm_ss, _, _, _ = compute_scs_frozen_core_unrelaxed_densities(
mo_eri, eps, n_occ, n_frozen=0, c_os=0.0, c_ss=1.0)
corr_os = opdm_os - D_ref
corr_ss = opdm_ss - D_ref
# corr_a == corr_os + corr_ss
assert np.allclose(corr_a, corr_os + corr_ss, atol=1e-10), \
"Correlation OPDM must be linear in (c_os, c_ss)"
# corr_b == 1.2*corr_os + (1/3)*corr_ss
assert np.allclose(corr_b, 1.2 * corr_os + (1.0 / 3.0) * corr_ss, atol=1e-10), \
"Correlation OPDM must scale linearly with SCS weights"Solution
import numpy as np
def compute_scs_frozen_core_unrelaxed_densities(mo_eri: np.ndarray,
epsilon: np.ndarray,
n_occ: int,
n_frozen: int,
c_os: float,
c_ss: float) -> tuple:
N = mo_eri.shape[0]
n_act = n_occ - n_frozen
n_vir = N - n_occ
# active-occupied slice and virtual slice
o = slice(n_frozen, n_occ)
v = slice(n_occ, N)
e_o = epsilon[n_frozen:n_occ]
e_v = epsilon[n_occ:N]
# Chemist (pq|rs) -> physicist <pr|qs> = (pq|rs) with swap of middle axes
# For MP2 we want <ij|ab> = (ia|jb)
# mo_eri[i,a,j,b] with i,j in active-occ and a,b in vir
iajb = mo_eri[o, v, o, v] # shape (nact, nvir, nact, nvir)
# Denominator D_{ijab} = e_i + e_j - e_a - e_b
D = (e_o[:, None, None, None] + e_o[None, None, :, None]
- e_v[None, :, None, None] - e_v[None, None, None, :])
# shape (nact, nvir, nact, nvir) matching iajb indexing (i,a,j,b)
if n_act > 0 and n_vir > 0:
t_iajb = iajb / D # t[i,a,j,b]
else:
t_iajb = np.zeros_like(iajb)
# Transform to t2[i,j,a,b]
# t_iajb is indexed (i,a,j,b); we want t2[i,j,a,b]
t2 = np.transpose(t_iajb, (0, 2, 1, 3)) # (i,j,a,b)
# Integrals (ia|jb) indexed (i,j,a,b):
g_ijab = np.transpose(iajb, (0, 2, 1, 3)) # (i,j,a,b)
# Opposite-spin and same-spin energies (closed shell RMP2):
# E_corr = sum_{ijab} (ia|jb) * [ 2 t[i,j,a,b] - t[i,j,b,a] ]
# Split: E_os = sum (ia|jb) * t[i,j,a,b] (counts the two OS spin combinations)
# E_ss = sum (ia|jb) * ( t[i,j,a,b] - t[i,j,b,a] )
# so that E_os + E_ss = sum (ia|jb) * (2 t[ijab] - t[ijba]) = plain RMP2.
if n_act > 0 and n_vir > 0:
E_os = float(np.einsum('ijab,ijab->', g_ijab, t2, optimize=True))
E_ss = float(np.einsum('ijab,ijab->', g_ijab, t2 - t2.transpose(0, 1, 3, 2), optimize=True))
else:
E_os = 0.0
E_ss = 0.0
# SCS-weighted amplitudes used for the density matrices.
# Define antisymmetrized combination T[i,j,a,b] = 2 t[i,j,a,b] - t[i,j,b,a].
# The OS part of this is 2*(1/2)*(t + t_swap) -> t + t_swap contribution split:
# Standard closed-shell decomposition:
# t_os[i,j,a,b] = t[i,j,a,b]
# t_ss[i,j,a,b] = t[i,j,a,b] - t[i,j,b,a]
# Then 2 t - t_swap = t_os + t_ss (since t_os - t_swap + t = 2t - t_swap) -- check:
# t_os + t_ss = t + (t - t_swap) = 2 t - t_swap. OK.
# SCS-weighted spin-summed amplitude:
if n_act > 0 and n_vir > 0:
T_scs = c_os * t2 + c_ss * (t2 - t2.transpose(0, 1, 3, 2))
else:
T_scs = np.zeros_like(t2)
# ------- OPDM -------
opdm = np.zeros((N, N))
# HF reference: doubly occupied orbitals (both frozen and active)
for p in range(n_occ):
opdm[p, p] = 2.0
if n_act > 0 and n_vir > 0:
# Standard closed-shell MP2 unrelaxed OPDM (active part only):
# P_ij = -2 sum_{kab} T_bar[i,k,a,b] * t[j,k,a,b]
# P_ab = 2 sum_{ijc} T_bar[i,j,a,c] * t[i,j,b,c]
# where T_bar = 2 t - t_swap. With SCS, replace T_bar by T_scs.
# Reference: standard RMP2 relaxed-density formulas.
P_ij = -2.0 * np.einsum('ikab,jkab->ij', T_scs, t2, optimize=True)
# symmetrize occ-occ
P_ij = 0.5 * (P_ij + P_ij.T)
P_ab = 2.0 * np.einsum('ijac,ijbc->ab', T_scs, t2, optimize=True)
P_ab = 0.5 * (P_ab + P_ab.T)
opdm[o, o] += P_ij
opdm[v, v] += P_ab
# enforce exact symmetry
opdm = 0.5 * (opdm + opdm.T)
# ------- TPDM (chemist notation) -------
tpdm = np.zeros((N, N, N, N))
# Reference TPDM for closed-shell HF in chemist notation (pq|rs):
# Gamma_ref[p,q,r,s] = 2 * delta_{pq in occ} delta_{rs in occ} * (... )
# Standard form: Gamma[p,q,r,s] = 2 D[p,q] D[r,s] - D[p,s] D[r,q]
# where D is the HF OPDM in chemist-style index pairs.
# Here D_ref[p,q] = 2 delta_{pq} for p in occupied (including frozen).
D_ref = np.zeros((N, N))
for p in range(n_occ):
D_ref[p, p] = 2.0
tpdm += 2.0 * np.einsum('pq,rs->pqrs', D_ref, D_ref, optimize=True) * 0.5
tpdm -= np.einsum('ps,rq->pqrs', D_ref, D_ref, optimize=True) * 0.5
# The factor 1/2 convention: use the common Gamma[p,q,r,s] = D[p,q] D[r,s] - 0.5 D[p,s] D[r,q]
# We'll reset and use that cleaner convention:
tpdm[:] = 0.0
tpdm += np.einsum('pq,rs->pqrs', D_ref, D_ref, optimize=True)
tpdm -= 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref, optimize=True)
# Correlation TPDM contribution (active/virtual only).
# In chemist notation (pq|rs), the MP2 unrelaxed TPDM correlation piece has
# Gamma_corr[i,a,j,b] += T_scs[i,j,a,b] (and symmetrizations to fill
# (a,i,j,b), (i,a,b,j), (a,i,b,j), plus the (r,s,p,q) symmetry).
if n_act > 0 and n_vir > 0:
# Build full-size correlation block indexed in chemist notation.
G = np.zeros((N, N, N, N))
# Place T_scs into (i,a,j,b):
G[o, v, o, v] += T_scs.transpose(0, 2, 1, 3) # from (i,j,a,b) to (i,a,j,b)
# Symmetrize under (p,q,r,s) <-> (r,s,p,q)
G = 0.5 * (G + G.transpose(2, 3, 0, 1))
tpdm += G
return opdm, tpdm, E_os, E_ssFrozen-Core CPHF Z-Vector Solve and Generalized Lagrangian for SCS-RMP2
You are given the unrelaxed one- and two-particle density matrices (OPDM and TPDM) from a frozen-core, spin-component-scaled RMP2 calculation, the chemist-notation MO two-electron integrals (pq|rs), a diagonal Fock matrix in the semi-canonical basis (its diagonal is also supplied as epsilon), and the frozen/active-occupied/virtual partitioning. Implement the coupled-perturbed HF (CPHF) Z-vector solve restricted to the active-occupied × virtual block, form the orbital-response (relaxed) OPDM, and build the generalized Lagrangian on the full MO space — including the frozen ↔ active-occupied coupling block that is generally nonzero under frozen-core SCS-RMP2.
def solve_frozen_core_zvector_and_lagrangian(opdm: np.ndarray, tpdm: np.ndarray, mo_eri: np.ndarray, fock_mo: np.ndarray, epsilon: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
"""
### Frozen-Core CPHF Z-Vector Solve and Generalized Lagrangian for SCS-RMP2 You are given the unrelaxed one- and two-particle density matrices (OPDM and TPDM) from a frozen-core, spin-component-scaled RMP2 calculation,
Returns:
The computed result described by the prompt.
"""
return relaxed_opdm, lagrangian, Z, z_residualBackground
Why a Z-vector?
The MP2 energy itself depends only on orbital energies and two-electron integrals in the MO basis, but its gradient with respect to nuclear coordinates requires knowing how the MOs respond to the perturbation. Rather than computing the orbital response (the U-matrix) for every perturbation, Handy and Schaefer's Z-vector trick reformulates the problem so that one solves a single linear system, independent of the number of perturbations. The unknown is a set of coefficients coupling occupied and virtual orbitals.
Frozen core
In frozen-core MP2 the deepest core orbitals are held fixed: no amplitudes involve them, and their density corrections vanish. Consequently, the CPHF Z-vector equation need only be solved on the active-occupied × virtual subspace. However, once core orbitals are frozen while correlation acts on the active occupied space, the Brillouin condition (which makes the occupied-occupied block of the generalized Fock matrix automatically symmetric at the HF level) no longer protects the frozen ↔ active-occupied block of the Lagrangian, and that coupling block must be built explicitly so that the final gradient contracts correctly against overlap derivatives.
Generalized Lagrangian
The generalized (or energy-weighted) Lagrangian is the quantity that multiplies overlap derivatives in analytic gradient theory. For correlated methods it is constructed from the relaxed OPDM, the TPDM, and the MO integrals. It is defined on the full MO space even when orbital response is restricted to an active subspace.
Function signature
def solve_frozen_core_zvector_and_lagrangian(opdm: np.ndarray, tpdm: np.ndarray, mo_eri: np.ndarray, fock_mo: np.ndarray, epsilon: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
return relaxed_opdm, lagrangian, Z, z_residualPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
import numpy as np
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def test_case_1():
# Basic shape + residual check for n_frozen=0
N, n_occ = 6, 3
eps, fock, eri = _build_random_system(N=N, n_occ=n_occ, seed=1)
opdm, tpdm, _, _ = _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0)
relaxed_opdm, lagrangian, Z, z_res = solve_frozen_core_zvector_and_lagrangian(
opdm, tpdm, eri, fock, eps, n_occ, 0)
nv = N - n_occ
nao = n_occ - 0
assert relaxed_opdm.shape == (N, N), "relaxed_opdm shape wrong"
assert lagrangian.shape == (N, N), "lagrangian shape wrong"
assert Z.shape == (nv, nao), f"Z shape wrong: got {Z.shape}"
assert z_res < 1e-8, f"z_residual too large: {z_res}"
import numpy as npTest case 2 of 3
import numpy as np
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def test_case_2():
# Frozen-core: Z must not touch frozen indices; relaxed_opdm frozen/virtual block zero.
N, n_occ, nf = 7, 4, 2
eps, fock, eri = _build_random_system(N=N, n_occ=n_occ, seed=2)
# Build an unrelaxed density consistent with frozen core by zeroing amplitudes touching frozen
# We re-use the n_frozen=0 builder but then zero the frozen blocks of the correlation correction.
opdm, tpdm, _, _ = _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0)
# Zero out correlation in frozen rows/cols of opdm: keep only the reference 2.0 on diag.
D_ref = np.zeros((N, N))
for i in range(n_occ):
D_ref[i, i] = 2.0
# frozen correction = opdm[frozen,:] - D_ref[frozen,:]; set to zero
for i in range(nf):
opdm[i, :] = D_ref[i, :]
opdm[:, i] = D_ref[:, i]
# Zero tpdm correction wherever any index is frozen: tpdm -> tpdm_ref on those entries
tpdm_ref = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
for i in range(nf):
tpdm[i, :, :, :] = tpdm_ref[i, :, :, :]
tpdm[:, i, :, :] = tpdm_ref[:, i, :, :]
tpdm[:, :, i, :] = tpdm_ref[:, :, i, :]
tpdm[:, :, :, i] = tpdm_ref[:, :, :, i]
relaxed_opdm, lagrangian, Z, z_res = solve_frozen_core_zvector_and_lagrangian(
opdm, tpdm, eri, fock, eps, n_occ, nf)
nv = N - n_occ
nao = n_occ - nf
assert Z.shape == (nv, nao), f"Z shape wrong under frozen core: {Z.shape}"
# frozen <-> virtual block of relaxed_opdm must be zero
assert np.allclose(relaxed_opdm[:nf, n_occ:], 0.0), "frozen-virtual block nonzero"
assert np.allclose(relaxed_opdm[n_occ:, :nf], 0.0), "virtual-frozen block nonzero"
# residual small
assert z_res < 1e-7, f"z_residual too large: {z_res}"
# Lagrangian must be full-shape and not identically zero in frozen rows (coupling block)
assert lagrangian.shape == (N, N)
import numpy as npTest case 3 of 3
import numpy as np
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def _build_random_system(N=6, n_occ=3, n_frozen=0, seed=0):
rng = np.random.default_rng(seed)
# Semi-canonical Fock: diagonal with occupied below virtual
eps = np.concatenate([np.sort(rng.uniform(-2.0, -0.3, size=n_occ)),
np.sort(rng.uniform(0.2, 2.5, size=N - n_occ))])
fock = np.diag(eps)
# Symmetric chemist-notation ERIs (pq|rs) = (rs|pq), (pq|rs)=(qp|rs)=(pq|sr)
A = rng.standard_normal((N, N, N, N)) * 0.05
eri = 0.5 * (A + A.transpose(1, 0, 2, 3))
eri = 0.5 * (eri + eri.transpose(0, 1, 3, 2))
eri = 0.5 * (eri + eri.transpose(2, 3, 0, 1))
return eps, fock, eri
def _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0, c_os=1.0, c_ss=1.0):
N = len(eps)
no = n_occ
nv = N - no
nf = n_frozen
nao = no - nf
o = slice(nf, no)
v = slice(no, N)
# Physicist integrals <ij||ab> etc. We use chemist -> physicist: <ij|ab> = (ia|jb)
iajb = eri[o, v, :, :][:, :, o, v] # shape (i,a,j,b)
# t2[i,j,a,b] amplitude with denominator e_i+e_j-e_a-e_b
e_o = eps[nf:no]
e_v = eps[no:N]
denom = (e_o[:, None, None, None] + e_o[None, :, None, None]
- e_v[None, None, :, None] - e_v[None, None, None, :])
# <ij|ab> in physicist = (ia|jb)
g_ijab = np.einsum('iajb->ijab', iajb)
t2 = g_ijab / denom
# Same-spin and opposite-spin energies (closed-shell formulation)
# E_os = sum (ia|jb) * t2[i,j,a,b]
# E_ss = sum [(ia|jb) - (ib|ja)] * t2[i,j,a,b]
E_os = float(np.einsum('iajb,ijab->', iajb, t2))
ibja = eri[o, v, :, :][:, :, o, v].transpose(0, 3, 2, 1) # (i,b,j,a) -> rearr to (i,a,j,b) order
# Simpler: compute via t2 antisym
t2_as = t2 - t2.transpose(0, 1, 3, 2)
E_ss = float(np.einsum('iajb,ijab->', iajb, t2_as))
# Build OPDM
opdm = np.zeros((N, N))
for i in range(no):
opdm[i, i] = 2.0
# Closed-shell MP2 unrelaxed OPDM:
# D_ij corr = -2 * sum_{kab} t_ik^ab * (2 t_jk^ab - t_jk^ba)
# D_ab corr = +2 * sum_{ijc} t_ij^ac * (2 t_ij^bc - t_ij^cb)
# With SCS scaling: the standard SCS gradient uses scaled amplitudes; for c_os=c_ss=1 this reduces to RMP2.
# For simplicity (and this is what the reference uses), scale t2 in the density contractions:
# Build "t_tilde" = c_os * t2_os + c_ss * t2_ss decomposition. For c_os=c_ss=1 this is just t2.
# We implement the c_os=c_ss=1 case here (test helper).
assert c_os == 1.0 and c_ss == 1.0, "test helper only supports unscaled"
T = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (2t - t_swap)
D_oo_corr = -2.0 * np.einsum('ikab,jkab->ij', t2, T)
D_vv_corr = 2.0 * np.einsum('ijac,ijbc->ab', t2, T)
opdm[nf:no, nf:no] += D_oo_corr
opdm[no:N, no:N] += D_vv_corr
# Symmetrize
opdm = 0.5 * (opdm + opdm.T)
# Build TPDM in chemist notation (pq|rs)
# Reference: Gamma_ref[p,q,r,s] = D_ref_pq D_ref_rs - 0.5 D_ref_ps D_ref_rq
D_ref = np.zeros((N, N))
for i in range(no):
D_ref[i, i] = 2.0
tpdm = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
# Correlation TPDM for closed-shell MP2 in chemist notation:
# Gamma_corr[i,a,j,b] = t2[i,j,a,b] scaled, etc.
# A compact standard form: gamma2[i,a,j,b] = 2*t2[i,j,a,b] - t2[i,j,b,a]
# (with appropriate symmetrizations). We add a symmetric piece so that
# 2*sum (pr|st) Gamma_corr pattern reproduces the MP2 energy.
Gcorr_iajb = 2.0 * t2 - t2.transpose(0, 1, 3, 2) # (i,j,a,b)
# Place into chemist slots (i,a,j,b) and its symmetry partner (j,b,i,a) and (a,i,b,j) (b,j,a,i)
tmp = np.zeros((N, N, N, N))
ia_o = list(range(nf, no))
va = list(range(no, N))
# Fill (i,a,j,b)
for ii, i in enumerate(ia_o):
for jj, j in enumerate(ia_o):
for aa, a in enumerate(va):
for bb, b in enumerate(va):
val = Gcorr_iajb[ii, jj, aa, bb]
tmp[i, a, j, b] += val
tmp[j, b, i, a] += val
tmp[a, i, b, j] += val
tmp[b, j, a, i] += val
tmp *= 0.25 # average the four equivalent placements
tpdm = tpdm + tmp
return opdm, tpdm, E_os, E_ss
def test_case_3():
# Placement of -Z in relaxed OPDM on active-occ/virtual off-diagonal block.
N, n_occ, nf = 6, 3, 1
eps, fock, eri = _build_random_system(N=N, n_occ=n_occ, seed=3)
opdm, tpdm, _, _ = _build_unrelaxed_mp2_densities(eri, eps, n_occ, n_frozen=0)
# Zero frozen blocks
D_ref = np.zeros((N, N))
for i in range(n_occ):
D_ref[i, i] = 2.0
for i in range(nf):
opdm[i, :] = D_ref[i, :]
opdm[:, i] = D_ref[:, i]
tpdm_ref = (np.einsum('pq,rs->pqrs', D_ref, D_ref)
- 0.5 * np.einsum('ps,rq->pqrs', D_ref, D_ref))
for i in range(nf):
tpdm[i, :, :, :] = tpdm_ref[i, :, :, :]
tpdm[:, i, :, :] = tpdm_ref[:, i, :, :]
tpdm[:, :, i, :] = tpdm_ref[:, :, i, :]
tpdm[:, :, :, i] = tpdm_ref[:, :, :, i]
relaxed_opdm, lagrangian, Z, z_res = solve_frozen_core_zvector_and_lagrangian(
opdm, tpdm, eri, fock, eps, n_occ, nf)
# Check: relaxed_opdm[active_occ, virtual] should equal -Z.T (shape (nao, nv))
block_ov = relaxed_opdm[nf:n_occ, n_occ:N]
assert np.allclose(block_ov, -Z.T, atol=1e-10), "relaxed_opdm active-occ/virtual block mismatch"
block_vo = relaxed_opdm[n_occ:N, nf:n_occ]
assert np.allclose(block_vo, -Z, atol=1e-10), "relaxed_opdm virtual/active-occ block mismatch"
# Diagonal blocks unchanged (match opdm)
assert np.allclose(relaxed_opdm[nf:n_occ, nf:n_occ], opdm[nf:n_occ, nf:n_occ], atol=1e-12)
assert np.allclose(relaxed_opdm[n_occ:N, n_occ:N], opdm[n_occ:N, n_occ:N], atol=1e-12)Solution
import numpy as np
def solve_frozen_core_zvector_and_lagrangian(opdm: np.ndarray, tpdm: np.ndarray, mo_eri: np.ndarray, fock_mo: np.ndarray, epsilon: np.ndarray, n_occ: int, n_frozen: int) -> tuple:
N = opdm.shape[0]
nf = n_frozen
no = n_occ
nao = no - nf # active occupied count
nv = N - no
o_act = slice(nf, no)
v = slice(no, N)
f = slice(0, nf)
o_all = slice(0, no)
# ---- Build generalized Fock-like intermediate F_pq = sum_rs [ 2(pq|rs) - (pr|qs) ] * D_rs_ref ...
# Actually we need the Fock contribution from the unrelaxed OPDM correction and the TPDM.
# The standard formulation: define I'_pq (the "Lagrangian intermediate") as
# I'_pq = sum_r h_pr D_rq + 2 sum_{rst} (pr|st) Gamma_{qr,st}
# where h is the one-electron part. But we don't have h directly; we have fock_mo.
# We will reconstruct via: the full Lagrangian from correlated theory is
# L_pq = F_pq^{eff} with appropriate contractions of opdm and tpdm and integrals.
# We adopt the standard closed-shell MP2 formulation using a generic construction:
# I_pq = sum_r D_pr F_rq + 2 sum_{rst} (ps|rt) Gamma_{qs,rt}
# but since fock includes all mean-field pieces, we use a cleaner split:
#
# Define the 2-index intermediate:
# Y_pq = 2 * sum_{rst} (pr|st) * tpdm_corr[q,r,s,t]
# plus orbital-energy pieces from opdm. For canonical HF reference, F is diagonal and
# F_pp = eps_p. We implement the generic form below.
# Split OPDM into reference + correction
opdm_ref = np.zeros_like(opdm)
for i in range(no):
opdm_ref[i, i] = 2.0
opdm_corr = opdm - opdm_ref
# Split TPDM into reference + correction.
# Closed-shell HF reference TPDM in chemist notation (pq|rs):
# Gamma_ref[p,q,r,s] = 2 D_ref_pq D_ref_rs - D_ref_ps D_ref_rq
# where D_ref = opdm_ref (so D_ref_ii = 2 for occupied).
# This gives Gamma_ref[i,i,j,j] = 2*2*2 - 2*2*delta_ij ... but let's just compute it.
tpdm_ref = (np.einsum('pq,rs->pqrs', opdm_ref, opdm_ref)
- 0.5 * np.einsum('ps,rq->pqrs', opdm_ref, opdm_ref))
tpdm_corr = tpdm - tpdm_ref
# Build the "I' intermediate" = generalized Fock-like matrix from correlated densities.
# In chemist (pq|rs) convention, the MP2 Lagrangian intermediate is:
# I_pq = sum_r F_pr opdm_corr[r,q] (the orbital-energy / OEI piece via fock)
# + 2 * sum_{rst} (pr|st) * tpdm_corr[q,r,s,t]
# We use fock_mo for the one-electron-like part contracted with opdm_corr.
I_int = np.einsum('pr,rq->pq', fock_mo, opdm_corr)
I_int += 2.0 * np.einsum('prst,qrst->pq', mo_eri, tpdm_corr)
# ---- Build CPHF RHS X on active-occ x vir: X_{a,i} = I_int[a,i] - I_int[i,a]
# (i active-occ, a virtual). Size (nv, nao).
X = I_int[v, o_act] - I_int[o_act, v].T
# ---- Build electronic Hessian G acting on active-occ x vir rotations.
# G_{ai,bj} Z_{bj} = (eps_a - eps_i) Z_{ai} + A_{ai,bj} Z_{bj}
# where A is the standard CPHF coupling:
# A_{ai,bj} = 4(ai|bj) - (ab|ij) - (aj|bi) [chemist notation (pq|rs)]
# We solve G Z = X.
eps = epsilon
e_o = eps[nf:no]
e_v = eps[no:N]
# denom[a,i] = e_v[a] - e_o[i]
denom = e_v[:, None] - e_o[None, :]
# Coupling: reshape into matrix of size (nv*nao, nv*nao)
# A[a,i,b,j] = 4(ai|bj) - (ab|ij) - (aj|bi)
eri_vo_vo = mo_eri[v, o_act, :, :][:, :, v, o_act] # (a,i,b,j)
eri_vv_oo = mo_eri[v, v, :, :][:, :, o_act, o_act] # (a,b,i,j)
eri_vo_ov = mo_eri[v, o_act, :, :][:, :, o_act, v] # (a,j? careful). We need (aj|bi)
# (aj|bi): indices p=a(v), q=j(o_act), r=b(v), s=i(o_act) -> shape (nv,nao,nv,nao)
eri_ajbi = mo_eri[v, o_act, :, :][:, :, v, o_act] # (a,j,b,i)
A = 4.0 * eri_vo_vo \
- np.einsum('abij->aibj', eri_vv_oo) \
- np.einsum('ajbi->aibj', eri_ajbi)
# Full operator G: diag part + A
dim = nv * nao
G = A.reshape(dim, dim).copy()
diag_add = denom.reshape(dim)
G[np.arange(dim), np.arange(dim)] += diag_add
X_flat = X.reshape(dim)
if dim == 0:
Z = np.zeros((nv, nao))
z_residual = 0.0
else:
Z_flat = np.linalg.solve(G, X_flat)
resid = G @ Z_flat - X_flat
z_residual = float(np.linalg.norm(resid))
Z = Z_flat.reshape(nv, nao)
# ---- Build relaxed OPDM
relaxed_opdm = opdm.copy()
# Place -Z on active-occ/virtual off-diagonal blocks (symmetric placement).
# Z is indexed (a,i) with a=virtual, i=active-occ.
# Convention: relaxed_opdm[i,a] = -Z[a,i], relaxed_opdm[a,i] = -Z[a,i]
relaxed_opdm[o_act, v] = -Z.T
relaxed_opdm[v, o_act] = -Z
# Frozen <-> virtual block: zero (no response there under frozen core)
relaxed_opdm[f, v] = 0.0
relaxed_opdm[v, f] = 0.0
# ---- Build generalized Lagrangian on the full MO space.
# L_pq = sum_r F_pr D^relaxed_rq + 2 sum_{rst} (pr|st) Gamma_corr[q,r,s,t]
# where we include the response piece through relaxed_opdm.
# Actually the standard Lagrangian uses the relaxed OPDM in the one-electron piece
# and the unrelaxed correlation TPDM (plus reference) in the two-electron piece.
# We build it generically:
lagrangian = np.einsum('pr,rq->pq', fock_mo, relaxed_opdm)
lagrangian += 2.0 * np.einsum('prst,qrst->pq', mo_eri, tpdm_corr)
# The frozen <-> active-occupied coupling block is automatically nonzero above because
# tpdm_corr and opdm_corr (through relaxed_opdm) have frozen indices with zero corrections,
# but the off-diagonal coupling arises from the TPDM contraction: elements where p is frozen
# and q is active-occ pick up contributions from (pr|st) Gamma_corr[q,r,s,t] with q active-occ.
# We leave this block as computed (do NOT zero out frozen rows/cols).
return relaxed_opdm, lagrangian, Z, z_residual