diff --git a/ccpy/drivers/diis.py b/ccpy/drivers/diis.py index 8768b6d..f4ae6b9 100644 --- a/ccpy/drivers/diis.py +++ b/ccpy/drivers/diis.py @@ -1,27 +1,30 @@ 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() @@ -29,39 +32,39 @@ def push(self, T, T_residuum, iteration): 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 diff --git a/ccpy/drivers/solvers.py b/ccpy/drivers/solvers.py index 3ec4b3f..bd15354 100644 --- a/ccpy/drivers/solvers.py +++ b/ccpy/drivers/solvers.py @@ -3,7 +3,7 @@ import numpy as np import h5py import signal -import os +import tempfile import sys from ccpy.utilities.printing import ( @@ -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): @@ -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?) @@ -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")