Skip to content

Commit

Permalink
Use tempfile names for out-of-core DIIS and Davidson.
Browse files Browse the repository at this point in the history
  • Loading branch information
kgururangan committed Dec 18, 2024
1 parent fa1ee45 commit 404d18f
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 43 deletions.
61 changes: 32 additions & 29 deletions ccpy/drivers/diis.py
Original file line number Diff line number Diff line change
@@ -1,67 +1,70 @@
import numpy as np
import h5py
import tempfile
from ccpy.utilities.utilities import remove_file


class DIIS:
def __init__(self, T, diis_size, out_of_core):

ftmp = tempfile.NamedTemporaryFile()
self.diis_size = diis_size
self.out_of_core = out_of_core
self.ndim = T.ndim
self.data_type = T.a.dtype
self.file_name = ftmp.name

if self.out_of_core:
remove_file("cc-diis-vectors.hdf5")
f = h5py.File("cc-diis-vectors.hdf5", "w")
self.T_list = f.create_dataset("t-vectors", (self.diis_size, self.ndim), dtype=np.float64)
self.T_residuum_list = f.create_dataset("resid-vectors", (self.diis_size, self.ndim), dtype=np.float64)
remove_file(self.file_name)
f = h5py.File(self.file_name, "w")
self.T_list = f.create_dataset("t-vectors", (self.diis_size, self.ndim), dtype=self.data_type)
self.T_residuum_list = f.create_dataset("resid-vectors", (self.diis_size, self.ndim), dtype=self.data_type)
else:
self.T_list = np.zeros((self.diis_size, self.ndim), dtype=np.float64)
self.T_residuum_list = np.zeros((self.diis_size, self.ndim), dtype=np.float64)
self.T_list = np.zeros((self.diis_size, self.ndim), dtype=self.data_type)
self.T_residuum_list = np.zeros((self.diis_size, self.ndim), dtype=self.data_type)

def cleanup(self):
if self.out_of_core:
remove_file("cc-diis-vectors.hdf5")
remove_file(self.file_name)

def push(self, T, T_residuum, iteration):
self.T_list[iteration % self.diis_size, :] = T.flatten()
self.T_residuum_list[iteration % self.diis_size, :] = T_residuum.flatten()

def extrapolate(self):
B_dim = self.diis_size + 1
B = -1.0 * np.ones((B_dim, B_dim))
B = -1.0 * np.ones((B_dim, B_dim), self.data_type)
for i in range(self.diis_size):
for j in range(i, self.diis_size):
B[i, j] = np.dot(self.T_residuum_list[i, :].T, self.T_residuum_list[j, :])
B[i, j] = np.dot(self.T_residuum_list[i, :].T.conj(), self.T_residuum_list[j, :])
B[j, i] = B[i, j]
B[-1, -1] = 0.0

rhs = np.zeros(B_dim)
rhs = np.zeros(B_dim, dtype=self.data_type)
rhs[-1] = -1.0

# TODO: replace with numpy.linalg.solve
# TODO: replace with scipy.linalg.lu
coeff = solve_gauss(B, rhs)
x_xtrap = np.zeros(self.ndim)
coeff = self.solve_gauss(B, rhs)
x_xtrap = np.zeros(self.ndim, dtype=self.data_type)
for i in range(self.diis_size):
x_xtrap += coeff[i] * self.T_list[i, :]

return x_xtrap

def solve_gauss(A, b):
"""DIIS helper function. Solves the linear system Ax=b using
Gaussian elimination"""
n = A.shape[0]
for i in range(n - 1):
for j in range(i + 1, n):
m = A[j, i] / A[i, i]
A[j, :] -= m * A[i, :]
b[j] -= m * b[i]
x = np.zeros(n)
k = n - 1
x[k] = b[k] / A[k, k]
while k >= 0:
x[k] = (b[k] - np.dot(A[k, k + 1 :], x[k + 1 :])) / A[k, k]
k = k - 1
def solve_gauss(self, A, b):
"""DIIS helper function. Solves the linear system Ax=b using
Gaussian elimination"""
n = A.shape[0]
for i in range(n - 1):
for j in range(i + 1, n):
m = A[j, i] / A[i, i]
A[j, :] -= m * A[i, :]
b[j] -= m * b[i]
x = np.zeros(n, dtype=self.data_type)
k = n - 1
x[k] = b[k] / A[k, k]
while k >= 0:
x[k] = (b[k] - np.dot(A[k, k + 1 :], x[k + 1 :])) / A[k, k]
k = k - 1

return x
return x
29 changes: 15 additions & 14 deletions ccpy/drivers/solvers.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import numpy as np
import h5py
import signal
import os
import tempfile
import sys

from ccpy.utilities.printing import (
Expand All @@ -15,16 +15,16 @@
from ccpy.drivers.diis import DIIS
from ccpy.utilities.utilities import remove_file

# Define a signal handler function to handle SIGINT
def signal_handler(sig, frame):
print("\nCtrl+C detected. Cleaning up...", end=" ")
remove_file("eomcc-vectors.hdf5")
remove_file("cc-diis-vectors.hdf5")
print("Cleanup complete.")
sys.exit(0)

# Register the signal handler for SIGINT
signal.signal(signal.SIGINT, signal_handler)
# # Define a signal handler function to handle SIGINT
# def signal_handler(sig, frame):
# print("\nCtrl+C detected. Cleaning up...", end=" ")
# remove_file("eomcc-vectors.hdf5")
# remove_file("cc-diis-vectors.hdf5")
# print("Cleanup complete.")
# sys.exit(0)
#
# # Register the signal handler for SIGINT
# signal.signal(signal.SIGINT, signal_handler)

# [TODO]: Add biorthogonal L and R single-root Davidson solver (non-Hermitian Hirao-Nakatsuji algorithm)
def eomcc_nonlinear_diis(HR, update_r, B0, R, dR, omega, T, H, X, fock, system, options):
Expand Down Expand Up @@ -101,9 +101,10 @@ def eomcc_davidson(HR, update_r, B0, R, dR, omega, T, H, system, options, t3_exc
print_eomcc_iteration_header()

# Create new HDF5 file by first checking if one exists and if so, remove it
remove_file("eomcc-vectors.hdf5")
# remove_file("eomcc-vectors.hdf5")
ftmp = tempfile.NamedTemporaryFile()
if options["davidson_out_of_core"]:
f = h5py.File("eomcc-vectors.hdf5", "w")
f = h5py.File(ftmp.name, "w")

# Maximum subspace size
nrest = 1 # number of previous vectors used to restart (>1 does not work, why?)
Expand Down Expand Up @@ -219,7 +220,7 @@ def eomcc_davidson(HR, update_r, B0, R, dR, omega, T, H, system, options, t3_exc
# store the actual root you've solved for
R.unflatten(r)
# remove HDF5 file
remove_file("eomcc-vectors.hdf5")
remove_file(ftmp.name)
# print the time taken for the root
minutes, seconds = divmod(time.perf_counter() - t_root_start, 60)
print(f" Completed in {minutes:.1f}m {seconds:.1f}s")
Expand Down

0 comments on commit 404d18f

Please sign in to comment.