"""
Self-consistent field (SCF) implementation for Non-Equilibrium Green's Function calculations.
This module provides the base NEGF class for performing self-consistent DFT+NEGF
calculations using Gaussian quantum chemistry package. It implements the energy-independent
self-energy approach (Damle et al. [1]) for molecular transport calculations. The module
handles:
- Integration with Gaussian for DFT calculations
- SCF convergence with Pulay DIIS mixing [2]
- Contact self-energy calculations
- Voltage bias and electric field effects
- Spin-polarized calculations
References
----------
[1] Damle, P., Ghosh, A. W., & Datta, S. (2002). First-principles analysis of molecular
conduction using quantum chemistry software. Chemical Physics, 281(2-3), 171-187.
DOI: 10.1016/S0301-0104(02)00496-2
[2] Pulay, P. (1980). Convergence acceleration of iterative sequences. The case of SCF
iteration. Chemical Physics Letters, 73(2), 393-398.
DOI: 10.1016/0009-2614(80)80396-4
"""
# Python packages
import numpy as np
import numpy.linalg as LA
import jax
import jax.numpy as jnp
from scipy import io
import os
import time
import matplotlib.pyplot as plt
# Enable double precision for accurate comparisons with NumPy
jax.config.update("jax_enable_x64", True)
# Gaussian interface packages
from gauopen import QCBinAr as qcb
# Developed packages
from gauNEGF.matTools import *
from gauNEGF.density import *
from gauNEGF.spinTools import spinorTF
# Use JAX functions directly
from gauNEGF.config import (SCF_CONVERGENCE_TOL, SCF_DAMPING, SCF_MAX_CYCLES,
FERMI_CALCULATION_TOL, PULAY_MIXING_SIZE, ENERGY_MIN)
from gauNEGF.utils import fractional_matrix_power, inv, eig
# CONSTANTS:
har_to_eV = 27.211386 # eV/Hartree
V_to_au = 0.03675 # Volts to Hartree/elementary Charge
[docs]
class NEGF(object):
"""
Non-Equilibrium Green's Function calculator integrated with Gaussian DFT.
Manages the self-consistent field cycle between NEGF transport calculations and
Gaussian DFT. Uses constant (energy-independent) self-energies; for
energy-dependent calculations see the NEGFE subclass.
Parameters
----------
fn : str
Base filename for Gaussian input/output files (without extension)
basis : str, optional
Gaussian basis set name (default: 'chkbasis')
func : str, optional
DFT functional to use (default: 'hf')
spin : {'r', 'u', 'ro', 'g'}, optional
Spin configuration:
- 'r': restricted
- 'u': unrestricted
- 'ro': restricted open-shell
- 'g': generalized open-shell (non-collinear)
(default: 'r')
fullSCF : bool, optional
Whether to run full SCF or use Harris approximation (default: True)
route : str, optional
Additional Gaussian route commands (default: '')
section : str, optional
Gaussian input section specification (default: None)
nPulay : int, optional
Number of previous iterations to use in Pulay mixing
(default: PULAY_MIXING_SIZE from gauNEGF.config)
Attributes
----------
F : ndarray
Fock matrix in Hartree (multiply by 27.211386 for eV)
P : ndarray
Density matrix
S : ndarray
Overlap matrix
fermi : float
Fermi energy in eV
nelec : float
Number of electrons
"""
[docs]
def __init__(self, fn, basis="chkbasis", func="hf", spin="r", fullSCF=True, route=None, section=None, nPulay=PULAY_MIXING_SIZE):
"""
Initialize NEGF calculator and run initial DFT calculation.
Sets up the calculator with specified parameters and runs an initial
DFT calculation using Gaussian to obtain the starting Fock and
overlap matrices.
Parameters are documented in class docstring.
"""
# Set up variables
self.ifile = fn + ".gjf"
self.chkfile = fn + ".chk"
self.ofile = fn + ".log"
self.func = func
self.basis = basis
self.method= spin+func
self.otherRoute = route # Other commands that are needed in Gaussian
self.section = section
self.spin = spin
self.energyDep = False;
self.Total_E_Old=9999.0;
#Default Integration Limits
self.Eminf = ENERGY_MIN
self.TSW = None
self.fSearch = None
self.fermi = None
self.updFermi = False
# Spin-locking attributes (optional feature, disabled by default)
self.spinLockEnabled = False
# Start calculation: Load Initial Matrices from Gaussian
print('Calculation started at '+str(time.asctime()))
print('Filename:',fn)
self.start_time = time.time()
self.bar = qcb.BinAr(debug=False,lenint=8,inputfile=self.ifile)
self.runDFT(fullSCF)
self.nae = int(self.bar.ne/2 + (self.bar.multip-1)/2)
self.nbe = int(self.bar.ne/2 - (self.bar.multip-1)/2)
# Prepare self.F, Density, self.S, and TF matrices
self.P = getDen(self.bar, spin)
self.F, self.locs = getFock(self.bar, spin)
self.nsto = len(self.locs)
Omat = jnp.array(self.bar.matlist["OVERLAP"].expand())
if spin == "ro" or spin == "u":
self.S = jnp.block([[Omat, jnp.zeros(Omat.shape)],[jnp.zeros(Omat.shape),Omat]])
else:
self.S = Omat
self.X = fractional_matrix_power(self.S, -0.5)
# Set Emin/Emax from orbitals
orbs, _ = eig(self.X @ self.F @ self.X)
self.Emin = min(orbs.real)*har_to_eV - 5
self.Emax = max(orbs.real)*har_to_eV
self.convLevel = 9999
self.MaxDP = 9999
# Pulay Mixing Initialization
self.pList = np.array([self.P for i in range(nPulay)], dtype=complex)
self.DPList = np.ones((nPulay, self.nsto, self.nsto), dtype=complex)*1e4
self.pMat = np.ones((nPulay+1, nPulay+1), dtype=complex)*-1
self.pMat[-1, -1] = 0
self.pB = np.zeros(nPulay+1)
self.pB[-1] = -1
# DFT Info dump
print("ORBS:")
print(self.locs)
self.Total_E = self.bar.scalar("escf")
self.updateN()
print('Expecting', str(self.bar.ne), 'electrons')
print('Actual: ', str(self.nelec), 'electrons')
print('Charge is:', self.bar.icharg)
print('Multiplicity is:', self.bar.multip)
print("Initial SCF energy: ", self.Total_E)
print('###################################')
[docs]
def runDFT(self, fullSCF=True):
"""
Run DFT calculation using Gaussian.
Performs either a full SCF calculation or generates an initial Harris
guess using Gaussian. Updates the Fock matrix and orbital indices
after completion.
Parameters
----------
fullSCF : bool, optional
If True, runs full SCF to convergence.
If False, uses Harris guess only. (default: True)
Notes
-----
- Attempts to load from checkpoint file first
- Falls back to full calculation if checkpoint fails
- Updates self.F and self.locs with new Fock matrix
"""
if fullSCF:
try:
print('Checking '+self.chkfile+' for saved data...');
self.bar.update(model=self.method, basis=self.basis, toutput=self.ofile, dofock=True,chkname=self.chkfile, miscroute=self.otherRoute, add_section=self.section)
except:
print('Checkpoint not loaded, running full SCF...');
self.bar.update(model=self.method, basis=self.basis, toutput=self.ofile, dofock="scf",chkname=self.chkfile, miscroute=self.otherRoute, add_section=self.section)
print("Done!")
self.F, self.locs = getFock(self.bar, self.spin)
else:
print('Using default Harris DFT guess to initialize...')
self.bar.update(model=self.method, basis=self.basis, toutput=self.ofile, dofock="GUESS",chkname=self.chkfile, miscroute=self.otherRoute, add_section=self.section)
self.bar.update(model=self.method, basis=self.basis, toutput=self.ofile, dofock=True, miscroute=self.otherRoute, add_section=self.section)
print("Done!")
self.F, self.locs = getFock(self.bar, self.spin)
[docs]
def updateN(self):
"""
Update the total number of electrons from the density matrix.
Calculates the number of electrons by tracing the product of
the density and overlap matrices. For restricted calculations,
multiplies by 2 to account for spin degeneracy.
Returns
-------
float
Total number of electrons
"""
nOcc = np.real(np.trace(self.P @ self.S))
if self.spin == 'r':
self.nelec = 2*nOcc
else:
self.nelec = nOcc
return self.nelec
[docs]
def setFock(self, F_):
"""
Set the Fock matrix, converting from eV to Hartree units.
Parameters
----------
F_ : ndarray
Fock matrix in eV units
"""
self.F = np.array(F_)/har_to_eV
[docs]
def updatePulay(self, nPulay):
"""
Reset Pulay mixing buffers to use a new history length.
Allows changing the number of previous iterations used in Pulay DIIS
mixing without reinitializing the full NEGF object. Buffers are
re-seeded from the current density matrix, matching the initialization
in __init__. Any accumulated mixing history is discarded.
Parameters
----------
nPulay : int
New number of previous iterations to use in Pulay mixing.
"""
if nPulay < 1:
raise ValueError(f'nPulay must be >= 1, got {nPulay}')
self.pList = np.array([self.P for i in range(nPulay)], dtype=complex)
self.DPList = np.ones((nPulay, self.nsto, self.nsto), dtype=complex)*1e4
self.pMat = np.ones((nPulay+1, nPulay+1), dtype=complex)*-1
self.pMat[-1, -1] = 0
self.pB = np.zeros(nPulay+1)
self.pB[-1] = -1
print(f'Pulay mixing buffers reset, nPulay = {nPulay}')
[docs]
def setDen(self, P_, enableSpinLock=False, spinLockList=None):
"""
Set the density matrix and update dependent quantities.
Updates the density matrix, stores it in Gaussian format,
recalculates the number of electrons, and updates the Fock matrix.
Optionally extracts and stores spin information for spin-locking.
Parameters
----------
P_ : ndarray
New density matrix
enableSpinLock : bool, optional
If True, extract and lock spin orientations from the density matrix.
This is an advanced feature for specialized use cases. Default: False
spinLockList : list, optional
List of atoms to apply spin-locking to
"""
self.P = P_
storeDen(self.bar, self.P, self.spin)
self.updateN()
print(f'Density matrix loaded, nelec = {self.nelec:.2f} electrons')
# Enable/disable spin-locking
if enableSpinLock:
if self.spin in ['g', 'u', 'ro']:
self.spinLockEnabled = True
# Store the original density as the locked reference
# This ensures we always lock to the same spin structure, not a drifting one
self.P_locked = self.P.copy()
if hasattr(spinLockList, '__iter__') and hasattr(self, 'locs'):
inds = np.where(np.isin(abs(self.locs), spinLockList, invert=True))[0]
self.P_locked[np.ix_(inds, inds)] = 0
print('Spin-locking enabled: will use original density as reference for spin orientations')
elif self.spin == 'r':
print('Warning: Spin-locking requested for restricted calculation, disabling')
self.spinLockEnabled = False
else:
# Disable spin-locking if not explicitly enabled
self.spinLockEnabled = False
if hasattr(self, 'P_locked'):
delattr(self, 'P_locked')
self.PToFock()
[docs]
def getHOMOLUMO(self):
"""
Calculate HOMO and LUMO energies.
Diagonalizes the Fock matrix in orthogonalized basis to get
orbital energies, then identifies HOMO and LUMO based on
electron occupation and spin configuration.
Returns
-------
ndarray
Array of [HOMO, LUMO] energies in eV
"""
orbs, _ = eig(self.X @ self.F @ self.X)
orbs = np.sort(orbs)*har_to_eV
if self.spin=='r':
homo_lumo = orbs[self.nae-1:self.nae+1].real
else:
homo_lumo = orbs[self.nae+self.nbe-1:self.nae+self.nbe+1].real
return homo_lumo
[docs]
def setVoltage(self, qV, fermi=None, Emin=None, Eminf=None):
"""
Set voltage bias and Fermi energy, updating electric field.
Applies a voltage bias between contacts and updates the chemical
potentials and electric field. Can optionally set the Fermi energy
and integration limits.
Parameters
----------
qV : float
Voltage bias in eV
fermi : float, optional
Fermi energy in eV. If not provided, will be calculated or
use existing value (default: None)
Emin : float, optional
Minimum energy for integration in eV (default: None)
Eminf : float, optional
Lower bound energy in eV (default: None)
Notes
-----
- Requires contacts to be set first
- Updates chemical potentials as fermi ± qV/2
- Calculates and applies electric field between contacts
- If fermi not provided, uses (HOMO+LUMO)/2 for initial guess
"""
# Check to make sure contacts set
assert hasattr(self, 'rInd') and hasattr(self,'lInd'), "Contacts not set!"
# Set Fermi Energy
if fermi is None:
if self.fermi is None:
self.updFermi = True
# Set initial fermi energy as (HOMO + LUMO)/2
homo_lumo = self.getHOMOLUMO()
print(f'Setting initial Fermi energy between HOMO ({homo_lumo[0]:.2f} eV) and LUMO ({homo_lumo[1]:.2f} eV)')
fermi = sum(homo_lumo)/2
else:
fermi = self.fermi
else:
self.updFermi = False
# Set Integration limits
if Emin!=None:
self.Emin = Emin
if Eminf!=None:
self.Eminf = Eminf
self.fermi = fermi
self.qV = qV
self.mu1 = fermi + (qV/2)
self.mu2 = fermi - (qV/2)
# Calculate electric field to apply during SCF
lAtoms = [self.bar.c[int(i-1)*3:int(i)*3] for i in self.lContact]
rAtoms = [self.bar.c[int(i-1)*3:int(i)*3] for i in self.rContact]
lCenter = np.mean(lAtoms, axis=0)
rCenter = np.mean(rAtoms, axis=0)
vec = lCenter-rCenter
dist = LA.norm(vec)
if dist == 0:
print("WARNING: left and right contact atoms identical, E-field set to zero!")
field = [0,0,0]
else:
vecNorm = vec/dist
field = -1*vecNorm*qV*V_to_au/(dist*0.0001)
self.bar.scalar("X-EFIELD", round(field[0]))
self.bar.scalar("Y-EFIELD", round(field[1]))
self.bar.scalar("Z-EFIELD", round(field[2]))
if not self.updFermi:
print("E-field set to "+str(LA.norm(field))+" au")
# TODO: n>2 terminal device support
[docs]
def setSigma(self, lContact=None, rContact=None, sig=-0.1j, sig2=None):
"""
Set self-energies for left and right contacts.
Configures the contact self-energies, handling various input formats
and spin configurations. Self-energies can be scalar, vector, or matrix,
with automatic handling of spin degrees of freedom.
Parameters
----------
lContact : array-like
Atom numbers for left contact, all atoms if None (default: None)
rContact : array-like
Atom numbers for right contact, all atoms if None (default: None)
sig : scalar or array-like, optional
Self-energy for left contact. Can be:
- scalar: same value for all orbitals
- vector: one value per orbital
- matrix: full self-energy matrix
(default: -0.1j)
sig2 : scalar, array-like, or None, optional
Self-energy for right contact. If None, uses sig.
Same format options as sig. (default: None)
Returns
-------
tuple
Left and right contact orbital indices
Notes
-----
- Handles spin configurations ('r', 'u', 'ro', 'g')
- Automatically expands scalar/vector inputs to full matrices
- Verifies matrix dimensions match Fock matrix
- Updates self.sigma1, self.sigma2, self.sigma12
Raises
------
Exception
If matrix dimensions don't match or invalid input format
"""
lInd, rInd = self.setContacts(lContact, rContact)
#Is there a second sigma matrix? If not, copy the first one
if sig2 is None:
sig2 = sig + 0.0
# Sigma can be a value, list, or matrix
if np.ndim(np.array(sig)) == 0 and np.ndim(np.array(sig2)) == 0:
pass
elif np.ndim(sig) == 1 and np.ndim(sig2) == 1:
if len(sig) == len(lInd) and len(sig2) == len(rInd):
pass
elif len(sig) == len(lInd)/2 and len(sig2) == len(rInd)/2:
if self.spin=='g':
sig = np.kron(sig, [1, 1])
sig2 = np.kron(sig2, [1, 1])
elif self.spin=='ro' or self.spin=='u':
sig = np.kron([1, 1], sig)
sig2 = np.kron([1, 1], sig2)
else:
print("Debug shape (sig1, inds1, sig2, inds2): ", sig.shape, lInd.shape, sig2.shape, rInd.shape)
raise Exception('Sigma matrix dimension mismatch!')
elif np.ndim(sig) == 2 and np.ndim(sig2) == 2:
if len(sig) == len(lInd) and len(sig2) == len(rInd):
pass
elif len(sig) == len(rInd)/2 and len(sig2) == len(rInd)/2:
if self.spin=='g':
sig = np.kron(sig, np.eye(2))
sig2 = np.kron(sig2, np.eye(2))
elif self.spin=='ro' or self.spin=='u':
sig = np.kron(np.eye(2), sig)
sig2 = np.kron(np.eye(2), sig2)
else:
print("Debug shape (sig1, inds1, sig2, inds2): ", sig.shape, lInd.shape, sig2.shape, rInd.shape)
raise Exception('Sigma matrix dimension mismatch!')
else:
raise Exception('Sigma matrix dimension mismatch!')
# Store Variables
self.rInd = rInd
self.lInd = lInd
self.sigma1 = formSigma(lInd, sig, self.nsto, self.S)
self.sigma2 = formSigma(rInd, sig2, self.nsto, self.S)
if self.sigma1.shape != self.F.shape or self.sigma2.shape != self.F.shape:
raise Exception(f'Sigma size mismatch! F shape={self.F.shape},'+
f' sigma shapes={self.sigma1.shape}, {self.sigma2.shape}')
self.sigma12 = self.sigma1 + self.sigma2
print('Max imag sigma:', str(np.max(np.abs(np.imag(self.sigma12)))));
self.Gam1 = (self.sigma1 - self.sigma1.conj().T)*1j
self.Gam2 = (self.sigma2 - self.sigma2.conj().T)*1j
return lInd, rInd
[docs]
def getSigma(self, E=0): #E only used by NEGFE() object, function inherited
return (self.sigma1, self.sigma2)
# Calculate density matrix from stored Fock matrix
[docs]
def FockToP(self):
"""
Calculate density matrix from Fock matrix using energy-independent approach.
Diagonalizes the orthogonalized Fock+Sigma matrix, integrates analytically
for constant self-energies (Damle et al. 2002), and transforms back to the
non-orthogonal basis. Updates self.P and optionally self.fermi.
Returns
-------
tuple
(eigenvalues, occupations) sorted by energy
"""
# Prepare Variables for Analytical Integration
X = np.array(self.X)
self.F, self.locs = getFock(self.bar, self.spin)
Fbar = X @ (self.F*har_to_eV + self.sigma12) @ X
GamBar1 = X @ self.Gam1 @ X
GamBar2 = X @ self.Gam2 @ X
D,V = jnp.linalg.eig(np.array(Fbar))
Vc = jnp.linalg.inv(V.conj().T)
#Update Fermi Energy (if not constant)
if self.updFermi:
Nexp = self.bar.ne
conv = min(self.convLevel, FERMI_CALCULATION_TOL)
if self.spin=='r':
Nexp/=2
self.fermi = bisectFermi(V,Vc,D,GamBar1+GamBar2,Nexp,conv, self.Eminf)
self.setVoltage(self.qV)
print(f'Fermi Energy set to {self.fermi:.2f} eV')
#Integrate to get density matrix
if self.mu1 == self.mu2:
P = density(V, Vc, D, GamBar1+GamBar2, self.Eminf, self.fermi)
else:
P1 = density(V, Vc, D, GamBar1, self.Eminf, self.mu1)
P2 = density(V, Vc, D, GamBar2, self.Eminf, self.mu2)
P = P1 + P2
# Calculate Level Occupation, Lowdin TF, Return
pshift = V.conj().T @ P @ V
self.P = X @ P @ X
occList = np.diag(np.real(pshift))
EList = np.array(np.real(D)).flatten()
inds = np.argsort(EList)
return EList[inds], occList[inds]
[docs]
def PMix(self, damping, Pulay=False):
"""
Mix old and new density matrices using damping or Pulay DIIS method.
Applies simple linear damping by default. When Pulay=True, applies the
DIIS update using the accumulated history buffer (length set in __init__).
Parameters
----------
damping : float
Mixing parameter between 0 and 1
Pulay : bool, optional
Whether to use Pulay mixing (default: False)
Returns
-------
tuple
(RMSDP, MaxDP) - RMS and maximum density matrix differences
"""
# Apply spin-locking if enabled (before mixing)
Pback = getDen(self.bar, self.spin)
if self.spinLockEnabled:
# Validate compatibility
if self.spin not in ['g', 'u', 'ro']:
print('Warning: Spin-locking not compatible with spin type, disabling')
self.spinLockEnabled = False
elif not hasattr(self, 'P_locked'):
print('Warning: P_locked not found, disabling spin-locking')
self.spinLockEnabled = False
else:
# Apply spin-locking projection using ORIGINAL locked density as reference
# NOT Pback, which is the previous iteration's mixed density
n_orbitals = self.P.shape[0] // 2
U = np.eye(self.P.shape[0], dtype=complex)
for i in range(n_orbitals):
# Extract 2x2 blocks based on spin type
if self.spin == 'g':
# Generalized: interleaved structure
block_current = self.P[2*i:2*i+2, 2*i:2*i+2]
block_locked = self.P_locked[2*i:2*i+2, 2*i:2*i+2]
elif self.spin in ['u', 'ro']:
# Unrestricted: construct from diagonal elements
block_current = np.array([[self.P[i, i], 0],
[0, self.P[n_orbitals+i, n_orbitals+i]]], dtype=complex)
block_locked = np.array([[self.P_locked[i, i], 0],
[0, self.P_locked[n_orbitals+i, n_orbitals+i]]], dtype=complex)
trace_current = np.trace(block_current)
trace_locked = np.trace(block_locked)
if abs(trace_current) > 1e-10 and abs(trace_locked) > 1e-10:
U_i = spinorTF(block_current, block_locked)
else:
U_i = np.eye(2) #identity
newSpinor = U_i@block_current@U_i.conj().T
if self.spin == 'g':
U[2*i:2*i+2, 2*i:2*i+2] = U_i
self.P = self.P.at[2*i:2*i+2, 2*i:2*i+2].set(newSpinor)
elif self.spin in ['u', 'ro']:
idx = np.ix_([i, n_orbitals+i], [i, n_orbitals+i])
U[idx] = U_i
self.P = self.P.at[idx].set(newSpinor)
# Store Old Density Info
Dense_old = np.diag(Pback)
Dense_diff = abs(np.diag(self.P) - Dense_old)
self.pList[1:, :, :] = self.pList[:-1, :, :]
self.pList[0, :, :] = Pback + damping*(self.P - Pback)
self.DPList[1:, :, :] = self.DPList[:-1, :, :]
self.DPList[0, :, :] = self.P - Pback
# Pulay Mixing
for i, v1 in enumerate(self.DPList):
for j, v2 in enumerate(self.DPList):
self.pMat[i,j] = np.sum(v1*v2)
# Apply Damping, store to Gaussian matrix
if Pulay:
coeff = LA.solve(self.pMat, self.pB)[:-1]
print("Applying Pulay Coeff: ", coeff)
self.P = sum([self.pList[i, :, :]*coeff[i] for i in range(len(coeff))])
self.pList[0, :, :] = self.P
else:
print("Applying Damping value=", damping)
self.P = self.pList[0, :, :]
storeDen(self.bar, self.P, self.spin)
# Update counters, print data
self.updateN()
print(f'Total number of electrons (NEGF): {self.nelec:.2f}')
self.MaxDP = max(Dense_diff)
RMSDP = np.sqrt(np.mean(Dense_diff**2))
print(f'MaxDP: {self.MaxDP:.2E} | RMSDP: {RMSDP:.2E}')
return RMSDP, self.MaxDP
# Use Gaussian to calculate the density matrix
[docs]
def PToFock(self):
"""
Calculate new Fock matrix from current density matrix using Gaussian.
Returns
-------
float
Energy difference from previous iteration
"""
# Run Gaussian, update SCF Energy
dE = -1
try:
self.bar.update(model=self.method, basis=self.basis, toutput=self.ofile, dofock="DENSITY", miscroute=self.otherRoute, add_section=self.section)
dE = 0
except Exception as e:
print("WARNING: DFT METHOD HAD AN ERROR, CYCLE INVALID:")
print(e)
print("CONTINUING TO NEXT CYCLE...")
self.Total_E_Old = self.Total_E.copy()
self.Total_E = self.bar.scalar("escf")
print("SCF energy: ", self.Total_E)
# Convergence variables: dE, RMSDP and MaxDP
dE = self.Total_E-self.Total_E_Old if dE !=-1 else -1
print(f'Energy difference is: {dE:.3E}')
return dE
# Main SCF loop, runs Fock <-> Density cycle until convergence reached
# Convergence criteria: dE, RMSDP, and MaxDP < conv, or maxcycles reached
[docs]
def SCF(self, conv=SCF_CONVERGENCE_TOL, damping=SCF_DAMPING, maxcycles=SCF_MAX_CYCLES, checkpoint=True, pulay=True):
"""
Run self-consistent field calculation until convergence.
Alternates FockToP / PMix / PToFock until all three criteria (dE, RMSDP,
MaxDP) fall below conv, or maxcycles is reached. Pulay DIIS is applied every
nPulay iterations when pulay=True.
Parameters
----------
conv : float, optional
Convergence criterion for energy and density (default: SCF_CONVERGENCE_TOL)
damping : float, optional
Mixing parameter between 0 and 1 (default: SCF_DAMPING)
maxcycles : int, optional
Maximum number of SCF cycles (default: SCF_MAX_CYCLES)
checkpoint : bool, optional
Save density matrix at each iteration and load if job interrupted (default: True)
pulay : bool, optional
Whether to use Pulay DIIS mixing (default: True)
Returns
-------
tuple
(count, PP, TotalE) - cycle number, number of electrons, and DFT energy at each cycle
"""
# Check to make sure contacts and voltage set
assert hasattr(self, 'mu1') and hasattr(self, 'mu2'), "Voltage not set!"
assert hasattr(self, 'rInd') and hasattr(self,'lInd'), "Contacts not set!"
# Find saved data from midrun
checkpoint_file = self.ifile[:-4]+"_P.mat"
final_file = self.ifile[:-4]+"_Final.mat"
if os.path.exists(checkpoint_file) and checkpoint:
try:
print(f"Found checkpoint file {checkpoint_file}, loading...")
self.setDen(io.loadmat(checkpoint_file)['den'])
self.fermi=io.loadmat(checkpoint_file)['fermi'][0][0]
self.setVoltage(self.qV)
except Exception as e:
print(f"Warning: checkpoint loaded - Error: {e}")
#Main SCF Loop
Loop = True
Niter = 0
minConv = 9999
PP=[]
count=[]
TotalE=[]
print('Entering NEGF-SCF loop at: '+str(time.asctime()))
print('###################################')
while Loop:
print()
print('Iteration '+str(Niter)+':')
# Run Pulay Kick every nPulay iterations, if turned on
isPulay = pulay*((Niter+1)%(len(self.pList)+1)==0)
# Fock --> P --> Fock
EList, occList = self.FockToP()
RMSDP, MaxDP = self.PMix(damping, isPulay)
dE = self.PToFock()
# Write monitor variables
TotalE.append(self.Total_E)
if dE == -1:
#Niter -= 1
self.convLevel = 9999
else:
count.append(Niter)
PP.append(self.nelec)
self.convLevel = max(RMSDP, MaxDP, abs(dE)/damping)
# Check 3 convergence criteria
if self.convLevel<conv:
print('##########################################')
print('Convergence achieved after '+str(Niter)+' iterations!')
Loop = False
elif Niter >= maxcycles:
print('##########################################')
print('WARNING: Convergence criterion not met, maxcycles reached!')
Loop = False
# Save progress
if self.convLevel < minConv and checkpoint:
print('Saving density checkpoint...')
io.savemat(checkpoint_file, {'den':self.P, 'conv':self.convLevel, 'fermi':self.fermi})
minConv = self.convLevel + 0.0
Niter += 1
if self.convLevel < conv and checkpoint:
os.system(f'mv {checkpoint_file} {final_file}')
print("--- %s seconds ---" % (time.time() - self.start_time))
print('')
print('SCF Loop exited at', time.asctime())
homo_lumo = self.getHOMOLUMO()
print(f'Predicted HOMO: {homo_lumo[0]:.2f} eV , Predicted LUMO {homo_lumo[1]:.2f} eV, Fermi: {self.fermi:0.2f}')
print('=========================')
print('ENERGY LEVEL OCCUPATION:')
print('=========================')
for pair in zip(occList, EList):
print(f"Energy = {pair[1]:9.3f} eV | Occ = {pair[0]:5.3f}")
print('=========================')
return count, PP, TotalE
[docs]
def writeChk(self):
"""
Write current state to a Gaussian binary checkpoint (.chk) file.
Uses gauopen's bar.writefile(), which routes through the BAF/unfchk
path internally -- no .fchk intermediate is needed.
"""
print(f'Writing {self.chkfile} via gauopen.writefile (BAF -> unfchk)...')
self.bar.writefile(self.chkfile)
print(self.chkfile + ' written!')
[docs]
def saveMAT(self, matfile="out.mat"):
"""
Save calculation results to MATLAB format file.
Parameters
----------
matfile : str, optional
Output filename (default: "out.mat")
Returns
-------
ndarray
Fock matrix in orthogonalized basis
"""
(sigma1, sigma2) = self.getSigma(self.fermi)
# Save data in MATLAB .mat file
matdict = {"F":self.F*har_to_eV, "sig1": sigma1, "sig2": sigma2,
"S": self.S, "fermi": self.fermi, "qV": self.qV,
"spin": self.spin, "den": self.P, "conv": self.convLevel}
io.savemat(matfile, matdict)
return self.X @ self.F @ self.X