"""Identify the basis of maximally spin-like orbitals from given set of 1-RDMs and 2-RDMs."""
# Copyright © 2021-2025 HQS Quantum Simulations GmbH. All Rights Reserved.
import numpy as np
from typing import List, Tuple, Union, cast
from opt_einsum import contract
from scipy.optimize import bisect, minimize, OptimizeResult, approx_fprime
from scipy.linalg import expm, pinv
from hqs_spin_mapper.spin_mapper_protocols import (
Supports_Spin_Optimization,
Supports_SW_Transformation,
InsufficientProtocolError,
)
from hqs_spin_mapper.transformables import SpinFinder
__all__ = [
"spin_orbital_optimization",
"rotate_rdm1",
"rotate_rdm2",
"basis_optimization",
"global_spin_orbital_optimization"
]
MATRIX_CONTRACTION = "ij, iq, jp -> qp"
TENSOR_CONTRACTION = "ijkl, iq, jp, kr, ls -> qprs"
class SpinNotFoundException(Exception):
r"""Custom exception that is raised if the search for spin-like orbitals failed."""
# TODO: Work out whether the set of spin_indices should be emptied before the parity optimization
[docs]
def spin_orbital_optimization(
transformable_system: Union[SpinFinder, Supports_SW_Transformation],
space_reduction: bool = True,
fast_mode: bool = False,
target: str = "minimum",
_delta_sufficiently_spin_like: float = 1e-3,
) -> None:
r"""Identify spin-like orbitals using an optimization of the local parities.
Performs a sequence of pairwise rotation between orbital indices `i` and `j` trying to extremize
the resulting local parity of the orbital `i` after rotation.
Args:
transformable_system (Supports_Spin_Optimization): Container object for the system data.
space_reduction (bool): Flag for performing rotations in the reduced search space only.
fast_mode (bool): Flag for performing approximate rotations of the RDMs to increase speed.
(This will result in slight inaccuracies of the RDMs.)
target (str): Determine the extremum to be determined (default=minimum)
_delta_sufficiently_spin_like (float): Orbitals i with this value for :math:`\delta` in
:math:`P_i = (-1 + \delta)` are exempt from optimization.
Raises:
InsufficientProtocolError: The transformable_system does not satisfy
Supports_SW_Transformation.
ValueError: Invalid target was specified.
"""
if not isinstance(transformable_system, Supports_Spin_Optimization):
raise InsufficientProtocolError(transformable_system, Supports_Spin_Optimization)
# Determine number of orbitals of the system
dim = transformable_system.system_size
immutable_indices = transformable_system.immutable_indices
# Full transformation matrix for the complete 1-RDMs and 2-RDMs (and later the Hamiltonians).
transformation = np.eye(dim)
# Collect tuple into one rdm1
if isinstance(transformable_system.rdm1, tuple):
rdm1 = transformable_system.rdm1[0] + transformable_system.rdm1[1]
else:
rdm1 = cast(np.ndarray, transformable_system.rdm1)
rdm2 = transformable_system.rdm2
# Pre-screen indices for spin-like orbitals and set aside sufficiently spin-(un)like orbitals
for i in range(dim):
parity = 1.0 - 2.0 * rdm1[i, i] + 4.0 * rdm2[i, i, i, i]
if abs(parity + 1.0) < _delta_sufficiently_spin_like:
transformable_system.spin_indices.add(i)
immutable_indices.add(i)
print(f"Initial parity P({i}): {parity} - orbital {i} will not be modified")
if abs(parity - 1.0) < _delta_sufficiently_spin_like:
immutable_indices.add(i)
# List of sites not containing immutable indices
sites = [n for n in range(dim) if n not in immutable_indices]
# Check whether orbital rotations should only be performed in a subspace of the single
# particle Hilbert space.
if space_reduction:
# Project the RDMs down to the reduced search space.
# Indices of the tensors in the subspace
ij_r: Tuple = np.ix_(sites, sites)
ijkl_r: Tuple = np.ix_(sites, sites, sites, sites)
rdm1_redux = rdm1[ij_r]
rdm2_redux = rdm2[ijkl_r]
# Use temporary RDMs for the trial rotations
rdm1_tmp = rdm1.copy()
rdm2_tmp = rdm2.copy()
for _ in range(transformable_system.optimization_steps):
for n_i, i in enumerate(sites):
if fast_mode and space_reduction:
# Fast mode rotations introduce small errors. We therefore perform a full rotation
# using the correct transformation matrix from the original rdms to continue with
# correct rdms for the next loop.
rdm1_redux = cast(
np.ndarray,
contract(MATRIX_CONTRACTION, rdm1_tmp, transformation, transformation),
)[ij_r]
rdm2_redux = cast(
np.ndarray,
contract(
TENSOR_CONTRACTION,
rdm2_tmp,
transformation,
transformation,
transformation,
transformation,
),
)[ijkl_r]
# The index needs to be at least one higher than "i".
n_j = n_i + 1
for _, j in enumerate(sites[(n_i + 1) :]):
# Orbital optimization in search space only.
if space_reduction:
# Determine the rotation operator between orbitals "i" and "j" that minimizes
# the parity of the resulting orbital "i".
_, U_phi = basis_optimization(n_i, n_j, rdm1_redux, rdm2_redux, target)
U_cp = np.eye(dim)
U_cp[ij_r] = U_phi
if fast_mode:
# Performance-optimized, but imprecise rotations of the RDMs
# The rotation functions are only adding onto entries in the rdms. To avoid
# double counting the diagonal entries of the identity matrix need to be
# subtracted.
U_phi[n_i, n_i] -= 1.0
U_phi[n_j, n_j] -= 1.0
rdm1_redux = rotate_rdm1(rdm1_redux, U_phi, n_i, n_j)
rdm2_redux = rotate_rdm2(rdm2_redux, U_phi, n_i, n_j)
# Re-add the identity contribution to the transformation matrix to retrieve
# the full rotation matrix
U_phi[n_i, n_i] += 1.0
U_phi[n_j, n_j] += 1.0
else:
rdm1 = cast(np.ndarray, contract(MATRIX_CONTRACTION, rdm1, U_cp, U_cp))
rdm2 = cast(
np.ndarray,
contract(TENSOR_CONTRACTION, rdm2, U_cp, U_cp, U_cp, U_cp),
)
rdm1_redux = rdm1[ij_r]
rdm2_redux = rdm2[ijkl_r]
# Update the full transformation matrix by the current two-site rotation matrix.
transformation = transformation @ U_cp
n_j += 1
# Rotation of the complete RDMs
else:
# Determine the rotation operator between orbitals "i" and "j" that minimizes
# the parity of the resulting orbital "i".
_, U_phi = basis_optimization(i, j, rdm1, rdm2, target)
U_cp = U_phi
if fast_mode:
# The rotation functions are only adding onto entries in the rdms. To avoid
# double counting the diagonal entries of the identity matrix need to be
# subtracted.
U_phi[i, i] -= 1.0
U_phi[j, j] -= 1.0
rdm1 = rotate_rdm1(rdm1, U_phi, i, j)
rdm2 = rotate_rdm2(rdm2, U_phi, i, j)
# Re-add the identity contribution to the transformation matrix to retrieve
# the full rotation matrix
U_phi[i, i] += 1.0
U_phi[j, j] += 1.0
else:
# Performance-optimized rotations of the RDMs but in the full Hilbert space.
rdm1 = cast(np.ndarray, contract(MATRIX_CONTRACTION, rdm1, U_cp, U_cp))
rdm2 = cast(
np.ndarray,
contract(TENSOR_CONTRACTION, rdm2, U_cp, U_cp, U_cp, U_cp),
)
# Update the full transformation matrix by the current two-site rotation matrix.
transformation = transformation @ U_cp
# Add sites that were previously declared as static to range spins.
# If these are sites that are positively not-spinlike the set of static sites has to
# be removed from range spins afterwards. (This might be reworked in the future, depending on
# the predominant use case of the static spins feature.)
# Track the minimum parity achieved by optimization procedure.
(parity, min_parity) = (1.0, 1.0)
for n_i, i in enumerate(sites):
# In the small subspace the rdm entries are adressed by the order index of subspace
# sites.
# TODO: Generalize the rdm to spin-dependent rdms
if space_reduction:
parity = 1.0 - 2.0 * rdm1_redux[n_i, n_i] + 4.0 * rdm2_redux[n_i, n_i, n_i, n_i]
# In the full space entries are adressed by their full index.
else:
parity = 1.0 - 2.0 * rdm1[i, i] + 4.0 * rdm2[i, i, i, i]
# Update the minimum parity achieved.
if parity < min_parity:
min_parity = parity
# Accept as spinlike, if the parity is sufficiently close to P=-1,
# i.e. P + 1 < \epsilon
if (target == "minimum" or target == "extremum") and abs(
parity + 1.0
) < transformable_system.tolerance:
print(f"Parity P({i}): {parity}")
# TODO: Discuss whether this should replace existing values
transformable_system.spin_indices.add(i)
# If fast_mode was used overwrite the entries of the full RDMs in the seach space
# with the entries from the RDMs in the search space
if space_reduction:
for n_j, j in enumerate(sites):
rdm1[i, j] = rdm1_redux[n_i, n_j]
# Update every possible configuration of indices "i" and "j" for the 2RDM
rdm2[i, i, i, i] = rdm2_redux[n_i, n_i, n_i, n_i]
rdm2[i, i, i, j] = rdm2_redux[n_i, n_i, n_i, n_j]
rdm2[i, i, j, i] = rdm2_redux[n_i, n_i, n_j, n_i]
rdm2[i, j, i, i] = rdm2_redux[n_i, n_j, n_i, n_i]
rdm2[j, i, i, i] = rdm2_redux[n_j, n_i, n_i, n_i]
rdm2[i, i, j, j] = rdm2_redux[n_i, n_i, n_j, n_j]
rdm2[i, j, i, j] = rdm2_redux[n_i, n_j, n_i, n_j]
rdm2[j, i, i, j] = rdm2_redux[n_j, n_i, n_i, n_j]
rdm2[i, j, j, i] = rdm2_redux[n_i, n_j, n_j, n_i]
rdm2[j, i, j, i] = rdm2_redux[n_j, n_i, n_j, n_i]
rdm2[j, j, i, i] = rdm2_redux[n_j, n_j, n_i, n_i]
rdm2[j, j, j, i] = rdm2_redux[n_j, n_j, n_j, n_i]
rdm2[j, j, i, j] = rdm2_redux[n_j, n_j, n_i, n_j]
rdm2[j, i, j, j] = rdm2_redux[n_j, n_i, n_j, n_j]
rdm2[i, j, j, j] = rdm2_redux[n_i, n_j, n_j, n_j]
rdm2[j, j, j, j] = rdm2_redux[n_j, n_j, n_j, n_j]
# If no spin-like sites were found (aside from static sites) raise a RuntimeError,
# because further calculation becomes pointless.
if target == "minimum":
if not transformable_system.spin_indices:
print(f"Smallest identified local parity: {min_parity}")
# Enter the found spin-like sites in the configuration data.
if isinstance(transformable_system.rdm1, tuple):
rdm1_uu = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[0], transformation, transformation
)
rdm1_dd = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[1], transformation, transformation
)
rdm1_ud = np.zeros_like(rdm1_uu)
if len(transformable_system.rdm1) == 3:
rdm1_ud = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[2], transformation, transformation
)
transformable_system.rdm1 = (rdm1_uu, rdm1_dd, rdm1_ud)
else:
transformable_system.rdm1 = rdm1
transformable_system.rdm2 = rdm2
transformable_system.rotation_matrix = transformable_system.rotation_matrix @ transformation
[docs]
def rotate_rdm1(rdm1: np.ndarray, U: np.ndarray, i: int, j: int) -> np.ndarray:
r"""Rotate the 1-RDM for pairwise rotations.
Args:
rdm1 (np.ndarray): Matrix description of the 1-RDM.
U (np.ndarray): Rotation matrix :math:`U_{iq}`.
i (int): First orbital index of the pair.
j (int): Second orbital index of the pair.
Returns:
np.ndarray: Rotated 1-RDM
"""
# Rotate the first index of the matrix. Project the first index of the 1RDM and the
# transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 1RDM.
rdm1 += cast(np.ndarray, contract("ij, iq -> qj", rdm1[[i, j], :], U[[i, j], :]))
# Rotate the second index of the matrix. Project the second index of the 1RDM and the
# first index of transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 1RDM.
rdm1 += cast(np.ndarray, contract("ij, jp -> ip", rdm1[:, [i, j]], U[[i, j], :]))
return rdm1
[docs]
def rotate_rdm2(rdm2: np.ndarray, U: np.ndarray, i: int, j: int) -> np.ndarray:
r"""Rotate the 2-RDM for pairwise rotations.
Args:
rdm2 (np.ndarray): Tensor description of the 2-RDM.
U (np.ndarray): Rotation matrix :math:`U_{iq}`.
i (int): First orbital index of the pair.
j (int): Second orbital index of the pair.
Returns:
np.ndarray: Rotated 2-RDM
"""
# Rotate the first index of the tensor. Project the first index of the 2RDM and the
# first index of transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 2RDM.
rdm2 += cast(np.ndarray, contract("ijkl, iq -> qjkl", rdm2[[i, j], :, :, :], U[[i, j], :]))
# Rotate the second index of the tensor. Project the second index of the 2RDM and the
# first index of transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 2RDM.
rdm2 += cast(np.ndarray, contract("ijkl, jp -> ipkl", rdm2[:, [i, j], :, :], U[[i, j], :]))
# Rotate the third index of the tensor. Project the third index of the 2RDM and the
# first index of transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 2RDM.
rdm2 += cast(np.ndarray, contract("ijkl, kr -> ijrl", rdm2[:, :, [i, j], :], U[[i, j], :]))
# Rotate the fourth index of the tensor. Project the fourth index of the 2RDM and the
# first index of transformation matrix onto the pair of indices "i" and "j"/
# Does not modify the other indices of the 2RDM.
rdm2 += cast(np.ndarray, contract("ijkl, ls -> ijks", rdm2[:, :, :, [i, j]], U[[i, j], :]))
return rdm2
[docs]
def basis_optimization(
i: int, j: int, rdm1: np.ndarray, rdm2: np.ndarray, target: str = "minimum"
) -> Tuple[float, np.ndarray]:
r"""Optimize the basis for extremal local parities via a single pairwise rotation of orbitals.
Args:
i (int): First orbital index
j (int): Second orbital index
rdm1 (np.ndarray): Matrix of 1-RDM
rdm2 (np.ndarray): 4-Tensor of the 2-RDM
target (str): Determine the extremum to be determined, i.e ``minimum``, ``maximum`` or
``extremum``.
Returns:
float: Optimal rotation angle :math:`\phi_{ij}`
np.ndarray: Rotation matrix :math:`U_{iq}`.
Raises:
ValueError: Invalid optimization target
"""
# Evaluate the analytic expression of the first derivative of the parity of orbital "i"
# with respect to the rotation angle "phi".
if target not in set({"minimum", "maximum", "extremum"}):
raise ValueError("Invalid optimization target!")
dP_dphi = lambda phi: (
(
+2 * np.sin(2 * phi) * rdm1[i, i]
- 2 * np.sin(2 * phi) * rdm1[j, j]
- 2 * np.cos(2 * phi) * (rdm1[i, j] + rdm1[j, i])
- 16 * np.cos(phi) ** 3 * np.sin(phi) * rdm2[i, i, i, i]
+ 16 * np.sin(phi) ** 3 * np.cos(phi) * rdm2[j, j, j, j]
)
+ (
2
* np.sin(4 * phi)
* (
rdm2[i, i, j, j]
+ rdm2[i, j, j, i]
+ rdm2[j, j, i, i]
+ rdm2[j, i, i, j]
+ rdm2[i, j, i, j]
+ rdm2[j, i, j, i]
)
)
+ (
4
* (3 * np.cos(phi) ** 2 * np.sin(phi) ** 2 - np.sin(phi) ** 4)
* (rdm2[j, j, j, i] + rdm2[j, j, i, j] + rdm2[j, i, j, j] + rdm2[i, j, j, j])
)
+ (
4
* (np.cos(phi) ** 4 - 3 * np.sin(phi) ** 2 * np.cos(phi) ** 2)
* (rdm2[i, i, i, j] + rdm2[i, i, j, i] + rdm2[i, j, i, i] + rdm2[j, i, i, i])
)
)
# Evaluate the analytic expression of the second derivative of the parity of orbital "i"
# with respect to the rotation angle "phi".
d2P_dphi2 = lambda phi: (
+4 * np.cos(2 * phi) * rdm1[i, i]
- 4 * np.cos(2 * phi) * rdm1[j, j]
+ 4 * np.sin(2 * phi) * (rdm1[i, j] + rdm1[j, i])
+ 4
* (
((-4 * np.cos(phi) ** 4 + 12 * np.sin(phi) ** 2 * np.cos(phi) ** 2) * rdm2[i, i, i, i])
+ (
(-4 * np.sin(phi) ** 4 + 12 * np.cos(phi) ** 2 * np.sin(phi) ** 2)
* rdm2[j, j, j, j]
)
+ (
(
2 * np.cos(phi) ** 4
+ 2 * np.sin(phi) ** 4
- 12 * np.sin(phi) ** 2 * np.cos(phi) ** 2
)
* (
rdm2[i, i, j, j]
+ rdm2[i, j, j, i]
+ rdm2[j, j, i, i]
+ rdm2[j, i, i, j]
+ rdm2[i, j, i, j]
+ rdm2[j, i, j, i]
)
)
+ (
(6 * np.sin(phi) * np.cos(phi) ** 3 - 10 * np.sin(phi) ** 3 * np.cos(phi))
* (rdm2[j, j, j, i] + rdm2[j, j, i, j] + rdm2[j, i, j, j] + rdm2[i, j, j, j])
)
+ (
(6 * np.cos(phi) * np.sin(phi) ** 3 - 10 * np.sin(phi) * np.cos(phi) ** 3)
* (rdm2[i, i, i, j] + rdm2[i, i, j, i] + rdm2[i, j, i, i] + rdm2[j, i, i, i])
)
)
)
dim = rdm1.shape[0]
U = np.eye(dim)
# Collect phi and P(phi) values for all roots of the first derivative of P(phi)
x_value: List[float] = []
y_value: List[float] = []
# Extend the search space for roots of dP/dphi from small phi step-wise to -pi < phi < pi,
# since the derivative is 2 pi-periodic.
for n in range(4):
x_a = -0.2 * np.pi - 0.19999 * n * np.pi
x_b = +0.2 * np.pi + 0.19999 * n * np.pi
# Check whether the sign of dP/dphi is different for the two ends of the search interval.
# If that is not the case, there will likely be no new roots in the search interval.
# Typically P(phi) evolves slowly with phi and shows no discontinuities.
if np.sign(dP_dphi(x_a)) != np.sign(dP_dphi(x_b)):
# Search for roots of dP/dphi
phi, results = bisect(dP_dphi, x_a, x_b, full_output=True)
# Consider only roots with positive second derivative, i.e. only true minima of P(phi).
if results.converged and abs(d2P_dphi2(phi)) > 1e-10:
# Create potential transformation matrix.
# Subtract value of 1. for use in optimized rotation.
U[i, i] = np.cos(phi) - 1.0
U[j, j] = np.cos(phi) - 1.0
U[i, j] = -1.0 * np.sin(phi)
U[j, i] = +1.0 * np.sin(phi)
# Perform trial rotation to find actual parity after rotation.
rdm1_r = rotate_rdm1(rdm1.copy(), U, i, j)
rdm2_r = rotate_rdm2(rdm2.copy(), U, i, j)
# Store the root "phi" in x_value and its corresponding parity "P" in y_value.
x_value.append(phi)
# Relevant criterion
if target == "minimum" or target == "maximum":
y_value.append(1.0 - 2.0 * rdm1_r[i, i] + 4.0 * rdm2_r[i, i, i, i])
if target == "extremum":
y_value.append(abs(1.0 - 2.0 * rdm1_r[i, i] + 4.0 * rdm2_r[i, i, i, i]))
# Select the root corresponding to the smallest minimum.
if y_value:
if target == "minimum":
phi = x_value[cast(int, np.argmin(y_value))]
if target == "maximum" or target == "extremum":
phi = x_value[cast(int, np.argmax(y_value))]
U[i, i] = np.cos(phi)
U[j, j] = np.cos(phi)
U[i, j] = -1.0 * np.sin(phi)
U[j, i] = +1.0 * np.sin(phi)
# If no roots have been found return an identity matrix as transformation.
else:
U = np.eye(dim)
phi = 0.0
return (phi, U)
def general_generator_matrix(angles: np.ndarray, index: int = 0) -> np.ndarray:
r"""Returns the rotation generator matrix for a given orbital index and vector of angles.
The function returns the rotation matrix
.. math::
S = i\sum_{i} \phi_i \sigma^{y}_{a i}
where :math:`a` denotes the chosen orbital index and :math:`\sigma^y_{a i}` a Pauli-Y like
matrix with nonzero entries for the index pairs :math:`a` and :math:`i`. The :math:`phi_i`
denote the respective rotation angles between indices :math:`a` and :math:`i`.
If :math:`i >= a` then :math:`i := i+1`.
Args:
angles (mp.ndarray): The rotation angles between the chosen orbital index and the
remaining orbital indices.
index (int): The orbital index to be modified by the rotation.
Returns:
np.ndarray: The rotation generator matrix corresponding to the rotation angles.
"""
generator_matrix = np.zeros((len(angles) + 1, len(angles) + 1), dtype=complex)
for n, phi in enumerate(angles):
if n >= index:
n = n + 1
# Index convention: U = U_(i->q)
generator_matrix[n, index] += phi * (+1j)
generator_matrix[index, n] += phi * (-1j)
return 1j * generator_matrix
def general_rotation_matrix(angles: np.ndarray, index: int = 0) -> np.ndarray:
r"""Returns the general rotation matrix for a given orbital index and vector of angles.
The function returns the rotation matrix
.. math::
S = \exp\left(i\sum_{i} \phi_i \sigma^{y}_{a i}\right)
where :math:`a` denotes the chosen orbital index and :math:`\sigma^y_{a i}` a Pauli-Y like
matrix with nonzero entries for the index pairs :math:`a` and :math:`i`, where are the remaining
indices. The :math:`phi_i` denote the respective rotation angles between indices :math:`a` and
`i`. If :math:`i >= a` then :math:`i := i+1`.
Args:
angles (np.ndarray): The rotation angles between the chosen orbital index and the
remaining orbital indices.
index (int): The orbital index to be modified by the rotation.
Returns:
np.ndarray: The rotation matrix corresponding to the rotation angles.
"""
generator_matrix = general_generator_matrix(angles, index)
return np.real(expm(generator_matrix))
def orbital_parity(
angles: np.ndarray, rdm1: np.ndarray, rdm2: np.ndarray, index: int = 0
) -> float:
r"""Returns the local parity of the chosen orbital after rotation with the specified angles.
The local parity is calculated for the orbital of the specified `index` after rotation with the
set of chosen rotation `angles`
Args:
angles (np.ndarray): The rotation angles between the chosen orbital index and the
remaining orbital indices.
rdm1 (np.ndarray): The one-particle reduced density matrix (1-RDM) as 2-tensor
rdm2 (np.ndarray): The two-particle reduced density matrix (2-RDM) as 4-tensor
index (int): The index of the orbital to be modified by the rotation.
Returns:
np.ndaray: The local parity for the specified orbital after a general rotation with the
chosen angles.
"""
U = general_rotation_matrix(angles, index)[:, index]
return cast(
float,
(
1.0
- 2.0 * contract("i, j, ij", U, U, rdm1)
+ 4.0 * contract("i, j, k, l, ijkl", U, U, U, U, rdm2)
),
)
def orbital_parity_gradient(
angles: np.ndarray, rdm1: np.ndarray, rdm2: np.ndarray, index: int = 0
) -> np.ndarray:
r"""Returns the gradient of the local parity of the orbital with respect to the angles.
Calculates the gradient of the local parity of the specified orbital with respect to the
rotation angles :math:`\nabla_{\phi_i}` at the specified values for the rotation angles
Args:
angles (np.ndarray): The rotation angles between the chosen orbital index and the
remaining orbital indices.
rdm1 (np.ndarray): The one-particle reduced density matrix (1-RDM) as 2-tensor
rdm2 (np.ndarray): The two-particle reduced density matrix (1-RDM) as 4-tensor
index (int): The index of the orbital to be modified by the rotation.
Returns:
np.ndarray: The gradient of the local orbital parity of the specified orbital with respect
to the rotation angles.
"""
gradient = np.zeros_like(angles)
U = general_rotation_matrix(angles, index)
G = general_generator_matrix(angles, index)
PU = U[:, index]
for n, _ in enumerate(angles):
q = n
pauli_y = np.zeros((len(angles) + 1, len(angles) + 1), dtype=complex)
if n >= index:
n = n + 1
pauli_y[n, index] += -1
pauli_y[index, n] += +1
############################################################################################
# Nth order Zassenhaus correction:
# https://en.wikipedia.org/wiki/Derivative_of_the_exponential_map
############################################################################################
N = 10
V = (-1.0 / 2) * (G @ pauli_y - pauli_y @ G)
dU = (U @ pauli_y) + (U @ V)
for order in range(3, N):
V = (-1.0 / order) * (G @ V - V @ G)
dU += U @ V
dU = np.real(dU)[:, index]
gradient[q] = cast(
float,
(
-2.0 * (contract("i, j, ij", dU, PU, rdm1) + contract("i, j, ij", PU, dU, rdm1))
+ 4.0
* (
contract("i, j, k, l, ijkl", dU, PU, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, dU, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, dU, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, PU, dU, rdm2)
)
),
)
return gradient
def orbital_parity_hessian(
angles: np.ndarray, rdm1: np.ndarray, rdm2: np.ndarray, index: int = 0
) -> np.ndarray:
r"""Returns the Hessian of the local parity of the orbital with respect to the angles.
(Currently still faulty and yielding wrong results!!!)
Calculates the Hessian of the local parity of the specified orbital with respect to the
rotation angles :math:`\nabla_{\phi_i}\nabla_{\phi_j}` at the specified values for the rotation
angles
Args:
angles (np.ndarray): The rotation angles between the chosen orbital index and the
remaining orbital indices.
rdm1 (np.ndarray): The one-particle reduced density matrix (1-RDM) as 2-tensor
rdm2 (np.ndarray): The two-particle reduced density matrix (1-RDM) as 4-tensor
index (int): The index of the orbital to be modified by the rotation.
Returns:
np.ndarray: The Hessian of the local orbital parity of the specified orbital with respect to
the rotation angles.
"""
hessian = np.zeros((len(angles), len(angles)))
U = general_rotation_matrix(angles, index)
G = general_generator_matrix(angles, index)
PU = U[:, index]
for n, _ in enumerate(angles):
q = n
if n >= index:
n = n + 1
pauli_y_n = np.zeros((len(angles) + 1, len(angles) + 1), dtype=complex)
pauli_y_n[n, index] += -1
pauli_y_n[index, n] += +1
V = (-1.0 / 2) * (G @ pauli_y_n - pauli_y_n @ G)
Vtn = pauli_y_n + V
dUn = (U @ Vtn)[:, index]
for m, _ in enumerate(angles):
p = m
if m >= index:
m = m + 1
pauli_y_m = np.zeros((len(angles) + 1, len(angles) + 1), dtype=complex)
pauli_y_m[m, index] += -1
pauli_y_m[index, m] += +1
########################################################################################
# Derivative of the Zassenhaus expansian of the gradient to second order
########################################################################################
V = (-1.0 / 2) * (G @ pauli_y_m - pauli_y_m @ G)
Vtm = pauli_y_m + V
dUm = U @ Vtm
dVtmn = (
pauli_y_m @ G @ pauli_y_n
+ G @ pauli_y_m @ pauli_y_n
- 2.0 * pauli_y_m @ pauli_y_n @ G
- 2.0 * G @ pauli_y_n @ pauli_y_m
+ pauli_y_n @ pauli_y_m @ G
+ pauli_y_n @ G @ pauli_y_m
)
d2U = (dUm @ Vtn + U @ dVtmn)[:, index]
dUm = dUm[:, index]
hessian[p, q] = cast(
float,
-2.0
* (
contract("i, j, ij", d2U, PU, rdm1)
+ contract("i, j, ij", dUn, dUm, rdm1)
+ contract("i, j, ij", dUm, dUn, rdm1)
+ contract("i, j, ij", PU, d2U, rdm1)
)
+ 4.0
* (
contract("i, j, k, l, ijkl", d2U, PU, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", dUn, dUm, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", dUn, PU, dUm, PU, rdm2)
+ contract("i, j, k, l, ijkl", dUn, PU, PU, dUm, rdm2)
+ contract("i, j, k, l, ijkl", dUm, dUn, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, d2U, PU, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, dUn, dUm, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, dUn, PU, dUm, rdm2)
+ contract("i, j, k, l, ijkl", dUm, PU, dUn, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, dUm, dUn, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, d2U, PU, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, dUn, dUm, rdm2)
+ contract("i, j, k, l, ijkl", dUm, PU, PU, dUn, rdm2)
+ contract("i, j, k, l, ijkl", PU, dUm, PU, dUn, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, dUm, dUn, rdm2)
+ contract("i, j, k, l, ijkl", PU, PU, PU, d2U, rdm2)
),
)
return hessian
[docs]
def global_spin_orbital_optimization(
transformable_system: Union[SpinFinder, Supports_SW_Transformation],
_algorithm: str = "trust-ncg",
_delta_sufficiently_spin_like: float = 1e-3,
_maxiter: int = 150,
) -> None:
r"""Determine spin-like orbitals via an optimization of local parities.
Performs a sequence of pairwise rotation between orbital indices :math:`i` and :math:`j` trying
to minimize the resulting local parity of the orbital :math:`i`.
Args:
transformable_system (Supports_Spin_Optimization): Container object for the system data.
_algorithm (str): Name of the scipy optimization algorithm. Default is ``trust-ncg``
_delta_sufficiently_spin_like (float): Orbitals i with :math:`P_i = (-1 + \delta)` are not
optimized
_maxiter (int): Maximum number of attempted iterations for the algorithm
Raises:
InsufficientProtocolError: transformable_system does not satisfy Supports_SW_Transformation.
"""
if not isinstance(transformable_system, Supports_Spin_Optimization):
raise InsufficientProtocolError(transformable_system, Supports_Spin_Optimization)
# Determine number of orbitals of the system
dim = transformable_system.system_size
immutable_indices = transformable_system.immutable_indices
# Full transformation matrix for the complete 1-RDMs and 2-RDMs (and later the Hamiltonians).
transformation = np.eye(dim)
################################################################################################
# Collect tuple into one rdm1
################################################################################################
if isinstance(transformable_system.rdm1, tuple):
rdm1 = transformable_system.rdm1[0] + transformable_system.rdm1[1]
else:
rdm1 = cast(np.ndarray, transformable_system.rdm1)
rdm2 = transformable_system.rdm2
################################################################################################
# Pre-screen indices for spin-like orbitals and set aside sufficiently spin-(un)like orbitals
################################################################################################
for i in range(dim):
parity = 1.0 - 2.0 * rdm1[i, i] + 4.0 * rdm2[i, i, i, i]
if abs(parity + 1.0) < _delta_sufficiently_spin_like:
transformable_system.spin_indices.add(i)
immutable_indices.add(i)
print(
f"Initial orbital parity P({i}): {parity} - the orbital {i} will not be modified"
)
if abs(parity - 1.0) < _delta_sufficiently_spin_like:
immutable_indices.add(i)
# TODO: Reintroduce this with a verbose flag
# print(
# f"Initial orbital parity P({i}): {parity} - the orbital {i} will not be modified"
# )
additional_spin_like_orbital_possible = True
counter = 0
while additional_spin_like_orbital_possible:
# List of site indices not contained in the immutable indices set
counter += 1
optimization_space = [n for n in range(dim) if n not in immutable_indices]
# Project the RDMs down to the reduced optimization space.
# ij_p and ijkl_p are the indices of the tensors in the subspace
ij_p: Tuple = np.ix_(optimization_space, optimization_space)
ijkl_p: Tuple = np.ix_(
optimization_space, optimization_space, optimization_space, optimization_space
)
projected_rdm1 = rdm1[ij_p]
projected_rdm2 = rdm2[ijkl_p]
############################################################################################
# Scipy minimization routine
############################################################################################
# TODO A more general ansatz is required here
ansatz = (np.pi / 1000) * np.ones(len(optimization_space) - 1)
if len(optimization_space) > 1:
index = 0
if index == len(optimization_space) - 1:
additional_spin_like_orbital_possible = False
########################################################################################
par_phi = lambda phi: orbital_parity(phi, projected_rdm1, projected_rdm2, index)
grad_phi = lambda phi: orbital_parity_gradient(
phi, projected_rdm1, projected_rdm2, index
)
hess_phi = lambda phi: (
0.5 * approx_fprime(phi, grad_phi, epsilon=0.000005).T
+ 0.5 * approx_fprime(phi, grad_phi, epsilon=0.000005)
)
########################################################################################
result = minimize(
par_phi,
ansatz,
method=_algorithm,
jac=grad_phi,
hess=hess_phi,
options={"maxiter": _maxiter},
)
if result.success:
if (
abs(orbital_parity(result.x, projected_rdm1, projected_rdm2) + 1)
< transformable_system.tolerance
):
print(f"The spin-like orbital {optimization_space[index]} has been identified")
immutable_indices.add(optimization_space[index])
transformable_system.spin_indices.add(optimization_space[index])
rotation_matrix = np.eye(dim)
rotation_matrix[ij_p] = general_rotation_matrix(result.x)
transformation = transformation @ rotation_matrix
rdm1 = contract(MATRIX_CONTRACTION, rdm1, rotation_matrix, rotation_matrix)
rdm2 = contract(
TENSOR_CONTRACTION,
rdm2,
rotation_matrix,
rotation_matrix,
rotation_matrix,
rotation_matrix,
)
if len(optimization_space) == 1:
additional_spin_like_orbital_possible = False
else:
immutable_indices.add(optimization_space[index])
if len(optimization_space) == 1:
additional_spin_like_orbital_possible = False
else:
immutable_indices.add(optimization_space[index])
if len(optimization_space) == 1:
additional_spin_like_orbital_possible = False
else:
result = OptimizeResult(x=np.array([]), success=True)
if (
abs(orbital_parity(result.x, projected_rdm1, projected_rdm2) + 1)
< transformable_system.tolerance
):
immutable_indices.add(optimization_space[0])
transformable_system.spin_indices.add(optimization_space[0])
rotation_matrix = np.eye(dim)
rotation_matrix[ij_p] = general_rotation_matrix(result.x)
transformation = transformation @ rotation_matrix
rdm1 = contract(MATRIX_CONTRACTION, rdm1, rotation_matrix, rotation_matrix)
rdm2 = contract(
TENSOR_CONTRACTION,
rdm2,
rotation_matrix,
rotation_matrix,
rotation_matrix,
rotation_matrix,
)
additional_spin_like_orbital_possible = False
if isinstance(transformable_system.rdm1, tuple):
rdm1_uu = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[0], transformation, transformation
)
rdm1_dd = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[1], transformation, transformation
)
rdm1_ud = np.zeros_like(rdm1_uu)
if len(transformable_system.rdm1) == 3:
rdm1_ud = contract(
MATRIX_CONTRACTION, transformable_system.rdm1[2], transformation, transformation
)
transformable_system.rdm1 = (rdm1_uu, rdm1_dd, rdm1_ud)
else:
transformable_system.rdm1 = rdm1
transformable_system.rdm2 = rdm2
transformable_system.rotation_matrix = transformable_system.rotation_matrix @ transformation