"""
Matrix manipulation utilities for quantum transport calculations.
This module provides helper functions for handling matrices in quantum transport
calculations, with a focus on:
- Self-energy matrix construction
- Density and Fock matrix manipulation
- Spin treatment (restricted, unrestricted, and generalized)
- Integration with Gaussian's matrix formats
The functions handle three types of spin treatments:
- 'r': Restricted (same orbitals for alpha and beta electrons)
- 'ro'/'u': Unrestricted (separate alpha and beta orbitals)
- 'g': Generalized (non-collinear spin treatment)
"""
# Python Packages
import numpy as np
# Gaussian interface packages
from gauopen import QCBinAr as qcb
from gauopen import QCOpMat as qco
# Matrix Headers
AlphaDen = "ALPHA DENSITY MATRIX"
BetaDen = "BETA DENSITY MATRIX"
AlphaSCFDen = "ALPHA SCF DENSITY MATRIX"
BetaSCFDen = "BETA SCF DENSITY MATRIX"
AlphaFock = "ALPHA FOCK MATRIX"
BetaFock = "BETA FOCK MATRIX"
AlphaEnergies = "ALPHA ORBITAL ENERGIES"
BetaEnergies = "BETA ORBITAL ENERGIES"
# CONSTANTS:
har_to_eV = 27.211386 # eV/Hartree
#### HELPER FUNCTIONS ####
[docs]
def getDen(bar, spin):
"""
Build density matrix from Gaussian checkpoint file.
Extracts the density matrix from a Gaussian checkpoint file based on
the specified spin treatment. Handles restricted, unrestricted, and
generalized cases.
Parameters
----------
bar : QCBinAr
Gaussian interface object
spin : str
Spin treatment to use:
- 'r': Restricted (same orbitals for alpha/beta)
- 'ro'/'u': Unrestricted (separate alpha/beta)
- 'g': Generalized (non-collinear)
Returns
-------
ndarray
Density matrix in appropriate format for spin treatment:
- Restricted: Single block
- Unrestricted: Two blocks (alpha/beta)
- Generalized: Single block with complex entries
Raises
------
ValueError
If spin treatment is not recognized
"""
if spin == "r" or spin == "g":
P = np.array(bar.matlist[AlphaSCFDen].expand())
elif spin == "ro" or spin == "u":
PA = np.array(bar.matlist[AlphaSCFDen].expand())
PB = np.array(bar.matlist[BetaSCFDen].expand())
P = np.block([[PA, np.zeros(PA.shape)], [np.zeros(PB.shape), PB]])
else:
raise ValueError("Spin treatment not recognized!")
return P
[docs]
def getFock(bar, spin):
"""
Build Fock matrix from Gaussian checkpoint file.
Extracts the Fock matrix and orbital indices from a Gaussian checkpoint
file based on the specified spin treatment. Handles restricted,
unrestricted, and generalized cases.
Parameters
----------
bar : QCBinAr
Gaussian interface object
spin : str
Spin treatment to use:
- 'r': Restricted (same orbitals for alpha/beta)
- 'ro'/'u': Unrestricted (separate alpha/beta)
- 'g': Generalized (non-collinear)
Returns
-------
tuple
(Fock, locs) where:
- Fock: ndarray, Fock matrix in appropriate format for spin treatment
- locs: ndarray, Orbital indices with positive for alpha/paired and
negative for beta orbitals
Raises
------
ValueError
If spin treatment is not recognized
"""
if spin == "r":
locs = bar.ibfatm
Fock = np.array(bar.matlist[AlphaFock].expand())
elif spin == "ro" or spin == "u":
locs = np.concatenate((bar.ibfatm, bar.ibfatm*-1))
AFock = np.array(bar.matlist[AlphaFock].expand())
BFock = np.array(bar.matlist[BetaFock].expand())
Fock = np.block([[AFock, np.zeros(AFock.shape)], [np.zeros(BFock.shape), BFock]])
elif spin == "g":
locs = [loc for pair in zip(bar.ibfatm, bar.ibfatm*-1) for loc in pair]
Fock = np.array(bar.matlist[AlphaFock].expand())
else:
raise ValueError("Spin treatment not recognized!")
locs = np.array(locs)
return Fock,locs
[docs]
def getEnergies(bar, spin):
"""
Get orbital energies from Gaussian checkpoint file.
Extracts orbital energies from a Gaussian checkpoint file based on
the specified spin treatment. Converts energies from Hartrees to eV.
Parameters
----------
bar : QCBinAr
Gaussian interface object
spin : str
Spin treatment to use:
- 'r': Restricted (same orbitals for alpha/beta)
- 'ro'/'u': Unrestricted (separate alpha/beta)
- 'g': Generalized (non-collinear)
Returns
-------
ndarray
All orbital energies in eV, sorted in ascending order. For
restricted/unrestricted spin, alpha and beta levels are interleaved
before the final sort, so each energy level appears twice; for
generalized spin, each level appears once.
Raises
------
ValueError
If spin treatment is not recognized
"""
if spin =="r":
Alevels = np.sort(bar.matlist[AlphaEnergies].expand())
levels = [level for pair in zip(Alevels, Alevels) for level in pair]
elif spin=="ro" or spin == "u":
Alevels = np.sort(bar.matlist[AlphaEnergies].expand())
Blevels = np.sort(bar.matlist[BetaEnergies].expand())
levels = [level for pair in zip(Alevels, Blevels) for level in pair]
elif spin=="g":
levels = np.sort(bar.matlist[AlphaEnergies].expand())
else:
raise ValueError("Spin treatment not recognized!")
return np.sort(levels)*har_to_eV
[docs]
def storeDen(bar, P, spin):
"""
Store density matrix in Gaussian checkpoint format.
Converts the density matrix to the appropriate format based on spin
treatment and stores it in the Gaussian checkpoint file. Handles
compression and proper typing of matrices.
Parameters
----------
bar : QCBinAr
Gaussian interface object
P : ndarray
Density matrix to store
spin : str
Spin treatment to use:
- 'r': Restricted (same orbitals for alpha/beta)
- 'ro'/'u': Unrestricted (separate alpha/beta)
- 'g': Generalized (non-collinear)
Notes
-----
For restricted calculations, the density is divided by 2 to account
for double occupation. For generalized calculations, complex matrices
are used.
Raises
------
ValueError
If spin treatment is not recognized
"""
nsto = len(bar.ibfatm)
if spin=="r":
P = np.real(np.array(P))
PaO = qco.OpMat(AlphaSCFDen,P/2,dimens=(nsto,nsto))
PaO.compress()
bar.addobj(PaO)
elif spin=="ro" or spin=="u":
P = np.real(np.array(P))
Pa = P[0:nsto, 0:nsto]
Pb = P[nsto:, nsto:]
PaO = qco.OpMat(AlphaSCFDen,Pa,dimens=(nsto,nsto))
PbO = qco.OpMat(BetaSCFDen,Pb,dimens=(nsto,nsto))
PaO.compress()
PbO.compress()
bar.addobj(PaO)
bar.addobj(PbO)
elif spin=="g":
P = np.complex128(np.array(P))
PaO = qco.OpMat(AlphaSCFDen,P,dimens=(nsto*2,nsto*2), typed='c')
PaO.compress()
bar.addobj(PaO)
else:
raise ValueError("Spin treatment not recognized!")