From a5ba0e56e4d3a4d68bca65e50bbf107190f113ce Mon Sep 17 00:00:00 2001 From: TianLin-228 Date: Fri, 16 Aug 2024 15:16:42 +0200 Subject: [PATCH 1/4] Davidson works when doubles folding in ADC(2) --- adcc/AdcMatrix.py | 52 ++++ adcc/solver/davidson.py | 449 ++++++++++++++++++++++++++-- adcc/solver/test_davidson_folded.py | 92 ++++++ adcc/test_workflow.py | 20 ++ adcc/workflow.py | 62 +++- 5 files changed, 634 insertions(+), 41 deletions(-) create mode 100644 adcc/solver/test_davidson_folded.py diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 016475b48..719073d06 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -716,3 +716,55 @@ def block_view(self, block): "projected ADC matrices.") # TODO The way to implement this is to ask the inner matrix to # a block_view and then wrap that in an AdcMatrixProjected. + + +class AdcMatrixFolded(AdcMatrix): + def __init__(self, matrix): + """ + Initialise a ADC(2) matrix when using doubles-folding. + + Parameters + ---------- + matrix : AdcMatrix + ADC(2) matrix + """ + super().__init__(matrix.method, matrix.ground_state, + block_orders=matrix.block_orders, + intermediates=matrix.intermediates) + + def update_omega(self, omega): + self.omega = omega + + def diagonal(self): + """ + Return the approximate diagonal of the ADC(2) matrix with doubles-folding. + """ + return AmplitudeVector(ph=super().diagonal().ph) + + def matvec(self, other): + """ + Compute the doubles-folded matrix-vector product of the singles vector with + an effective ADC matrix which depends on the eigenvalue ω. + """ + diag = super().diagonal().pphh + e = diag.ones_like() + u2 = self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag) + return AmplitudeVector(ph=self.block_apply("ph_ph", other.ph) + + self.block_apply("ph_pphh", u2)) + + def __matmul__(self, other): + if isinstance(other, AmplitudeVector): + return self.matvec(other) + if isinstance(other, list): + if all(isinstance(elem, AmplitudeVector) for elem in other): + return [self.matvec(v) for v in other] + return NotImplemented + + def unfold(self, u1): + """ + recompute the doubles component and return the complete vector. + """ + diag = super().diagonal().pphh + e = diag.ones_like() + u2 = self.block_apply("pphh_ph", u1.ph) / (e * self.omega - diag) + return AmplitudeVector(ph=u1.ph, pphh=u2) diff --git a/adcc/solver/davidson.py b/adcc/solver/davidson.py index 4819de69c..ad2a91823 100644 --- a/adcc/solver/davidson.py +++ b/adcc/solver/davidson.py @@ -27,13 +27,14 @@ import scipy.sparse.linalg as sla from adcc import evaluate, lincomb -from adcc.AdcMatrix import AdcMatrixlike +from adcc.AdcMatrix import AdcMatrixlike, AdcMatrixFolded from adcc.AmplitudeVector import AmplitudeVector from .common import select_eigenpairs from .preconditioner import JacobiPreconditioner from .SolverStateBase import EigenSolverStateBase from .explicit_symmetrisation import IndexSymmetrisation +from itertools import product class DavidsonState(EigenSolverStateBase): @@ -42,6 +43,23 @@ def __init__(self, matrix, guesses): self.residuals = None # Current residuals self.subspace_vectors = guesses.copy() # Current subspace vectors self.algorithm = "davidson" + self.DIIS_iter = 0 # Total number of DIIS iterations + self.folded = False # Folded or normal ADC matrix + + +class FoldedDavidsonState(DavidsonState): + def __init__(self, matrix, guesses): + super().__init__(matrix, guesses) # with folded ADC matrix + self.folded = True # Folded or normal ADC matrix + self.n_state = None # Current state + self.macro_iter = 0 # Number of macro iterations + self.history_rval = [] # Previous Ritz values + self.energy_diff = None # Difference between Ritz values + self.residual_norm = None # Current residual norm of one state + self.residual = None # Current residual of one state + self.eigenvector = None # Current eigenvector of one state + self.converged_macro = False # Convergence of macro iteration + self.converged_diis = False # Convergence of DIIS def default_print(state, identifier, file=sys.stdout): @@ -53,15 +71,27 @@ def default_print(state, identifier, file=sys.stdout): # TODO Use colour! if identifier == "start" and state.n_iter == 0: - print("Niter n_ss max_residual time Ritz values", - file=file) + if not state.folded: + print("Niter n_ss max_residual time Ritz values", file=file) + else: + print("Niter n_ss energy_difference time Ritz values", + file=file) elif identifier == "next_iter": time_iter = state.timer.current("iteration") - fmt = "{n_iter:3d} {ss_size:4d} {residual:12.5g} {tstr:5s}" - print(fmt.format(n_iter=state.n_iter, tstr=strtime_short(time_iter), - ss_size=len(state.subspace_vectors), - residual=np.max(state.residual_norms)), - "", state.eigenvalues[:7], file=file) + if not state.folded: + fmt = "{n_iter:3d} {ss_size:4d} {residual:12.5g} {tstr:5s}" + print(fmt.format(n_iter=state.n_iter, + tstr=strtime_short(time_iter), + ss_size=len(state.subspace_vectors), + residual=np.max(state.residual_norms)), + "", state.eigenvalues[:10], file=file) + else: + fmt = "{n_iter:3d} {ss_size:4d} {energy_diff:12.5g} {tstr:5s}" + print(fmt.format(n_iter=state.n_iter, + tstr=strtime_short(time_iter), + ss_size=len(state.subspace_vectors), + energy_diff=state.energy_diff), + "", state.eigenvalues[:10], file=file) if hasattr(state, "subspace_orthogonality"): print(33 * " " + "nonorth: {:5.3g}" "".format(state.subspace_orthogonality)) @@ -72,6 +102,43 @@ def default_print(state, identifier, file=sys.stdout): print(" Total solver time: ", strtime(soltime), file=file) elif identifier == "restart": print("=== Restart ===", file=file) + # For folded matrix + elif identifier == "folded_start": + print(f"============ State {state.n_state} ============", file=file) + print("folded_matrix.omega(initial):", state.matrix.omega, file=file) + elif identifier == "micro": + print(f"--macro iteration {state.macro_iter}, ", + f"Number of micro iterations: {state.n_iter}, " + f"Ritz_value: {state.eigenvalues[state.n_state]}, ", + f"residual_norm: {state.residual_norm}, ", + f"Converged: {state.converged}--", file=file) + elif identifier == "macro_stop": + print("== Summary of macro iterations ==", file=file) + print(" Number of macro-iterations:", state.macro_iter, + " n_applies:", state.n_applies, file=file) + print(" Ritz value:", state.eigenvalues[state.n_state], + " residual norm:", state.residual_norm, + " Converged or not:", state.converged_macro, file=file) + print(" time:", strtime_short(state.timer.current + ("folded iterations")), file=file) + elif identifier == "DIIS_steps": + print(f"--DIIS, Omega: {state.matrix.omega}, " + f"residual_norm: {state.residual_norm}--", file=file) + elif identifier == "DIIS_stop": + print("== Summary of DIIS ==", file=file) + print(" Number of DIIS iterations:", state.DIIS_iter, file=file) + print(" Omega:", state.eigenvalues[state.n_state], + " residual norm:", state.residual_norm, + " Converged or not:", state.converged_diis, file=file) + print(" time:", strtime_short(state.timer.current + ("folded iterations")), file=file) + elif identifier == "sum_folded": + print("========= Converged (folded matrix) =========", file=file) + print(" Number of matrix applies: ", state.n_applies, file=file) + print(" Total solver time: ", strtime(state.timer.total + ("folded iterations")), file=file) + print(" Number of Davidson iterations: ", state.n_iter, file=file) + print(" Number of DIIS: ", state.DIIS_iter, file=file) # TODO This function should be merged with eigsh @@ -127,6 +194,9 @@ def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep, if callback is None: def callback(state, identifier): pass + if state.folded: # initialization for each micro iteration + state.history_rval = [matrix.omega] + state.n_iter = 0 # The problem size n_problem = matrix.shape[1] @@ -196,26 +266,31 @@ def form_residual(rval, rvec): epair_mask = select_eigenpairs(rvals, n_ep, which) state.eigenvalues = rvals[epair_mask] state.residuals = [residuals[i] for i in epair_mask] - state.residual_norms = np.array([r @ r for r in state.residuals]) - # TODO This is misleading ... actually residual_norms contains - # the norms squared. That's also the used e.g. in adcman to - # check for convergence, so using the norm squared is fine, - # in theory ... it should just be consistent. I think it is - # better to go for the actual norm (no squared) inside the code - # - # If this adapted, also change the conv_tol to tol conversion - # inside the Lanczos procedure. + if not state.folded: + state.residual_norms = np.array([np.sqrt(r @ r) + for r in state.residuals]) + else: + state.energy_diff = np.abs(state.eigenvalues[state.n_state] + - state.history_rval[-1]) + state.history_rval.append(state.eigenvalues[state.n_state]) callback(state, "next_iter") state.timer.restart("iteration") if is_converged(state): # Build the eigenvectors we desire from the subspace vectors: - state.eigenvectors = [lincomb(v, SS, evaluate=True) - for i, v in enumerate(np.transpose(rvecs)) - if i in epair_mask] - + if not state.folded: + state.eigenvectors = [lincomb(v, SS, evaluate=True) + for i, v in enumerate(np.transpose(rvecs)) + if i in epair_mask] + callback(state, "is_converged") + else: + # update guesses vectors for next macro iteration + state.subspace_vectors = [lincomb(v, SS, evaluate=True) + for v in np.transpose(rvecs)] + assert len(state.subspace_vectors) == n_block + state.eigenvector = state.subspace_vectors[epair_mask + [state.n_state]] state.converged = True - callback(state, "is_converged") state.timer.stop("iteration") return state @@ -223,9 +298,17 @@ def form_residual(rval, rvec): warnings.warn(la.LinAlgWarning( f"Maximum number of iterations (== {max_iter}) " "reached in davidson procedure.")) - state.eigenvectors = [lincomb(v, SS, evaluate=True) - for i, v in enumerate(np.transpose(rvecs)) - if i in epair_mask] + if not state.folded: + state.eigenvectors = [lincomb(v, SS, evaluate=True) + for i, v in enumerate(np.transpose(rvecs)) + if i in epair_mask] + else: + # update guesses vectors for next macro iteration + state.subspace_vectors = [lincomb(v, SS, evaluate=True) + for v in np.transpose(rvecs)] + assert len(state.subspace_vectors) == n_block + state.eigenvector = state.subspace_vectors[epair_mask + [state.n_state]] state.timer.stop("iteration") state.converged = False return state @@ -298,9 +381,18 @@ def form_residual(rval, rvec): if n_ss_added == 0: state.timer.stop("iteration") state.converged = False - state.eigenvectors = [lincomb(v, SS, evaluate=True) - for i, v in enumerate(np.transpose(rvecs)) - if i in epair_mask] + if not state.folded: + state.eigenvectors = [lincomb(v, SS, evaluate=True) + for i, v in enumerate(np.transpose(rvecs)) + if i in epair_mask] + else: + # Compute all eigenvectors as guesses vectors + # for the next macro iteration. + state.subspace_vectors = [lincomb(v, SS, evaluate=True) + for v in np.transpose(rvecs)] + assert len(state.subspace_vectors) == n_block + state.eigenvector = state.subspace_vectors[epair_mask + [state.n_state]] warnings.warn(la.LinAlgWarning( "Davidson procedure could not generate any further vectors for " "the subspace. Iteration cannot be continued like this and will " @@ -330,7 +422,7 @@ def eigsh(matrix, guesses, n_ep=None, max_subspace=None, max_subspace : int or NoneType, optional Maximal subspace size conv_tol : float, optional - Convergence tolerance on the l2 norm squared of residuals to consider + Convergence tolerance on the l2 norm of residuals to consider them converged which : str, optional Which eigenvectors to converge to (e.g. LM, LA, SM, SA) @@ -402,6 +494,232 @@ def convergence_test(state): return state +def eigsh_folded(matrix, guesses, omegas=None, n_ep=None, max_subspace=None, + conv_tol=1e-9, which="SA", max_iter=70, + callback=None, preconditioner=None, + preconditioning_method="Davidson", + debug_checks=False, residual_min_norm=None, + explicit_symmetrisation=IndexSymmetrisation, + macro_conv_tol=1e-3, macro_max_iter=30, + num_diis_vecs=50, diis_max_iter=200): + """Davidson eigensolver for ADC problems with doubles-folding + + Parameters + ---------- + matrix + ADC(2) matrix instance + guesses : list + Guess vectors (fixes also the Davidson block size) + n_ep : int or NoneType, optional + Number of eigenpairs to be computed + max_subspace : int or NoneType, optional + Maximal subspace size + conv_tol : float, optional + Convergence tolerance on the l2 norm of residuals to consider + them converged during the final DIIS iterations. + which : str, optional + Which eigenvectors to converge to (e.g. LM, LA, SM, SA) + max_iter : int, optional + Maximal number of Davidson iterations + callback : callable, optional + Callback to run after each iteration + preconditioner + Preconditioner (type or instance) + preconditioning_method : str, optional + Precondititoning method. Valid values are "Davidson" + or "Sleijpen-van-der-Vorst" + explicit_symmetrisation + Explicit symmetrisation to apply to new subspace vectors before + adding them to the subspace. Allows to correct for loss of index + or spin symmetries (type or instance) + debug_checks : bool, optional + Enable some potentially costly debug checks + (Loss of orthogonality etc.) + residual_min_norm : float or NoneType, optional + Minimal norm a residual needs to have in order to be accepted as + a new subspace vector + (defaults to 2 * len(matrix) * machine_expsilon) + macro_conv_tol : float, optional (default=1e-3) + Convergence tolerance on the l2 norm of residuals to consider + them converged during macro iterations. + macro_max_iter : int, optional (default=30) + Maximal number of macro iterations + num_diis_vecs: int, optional (default=50) + Maximal number of DIIS vectors to keep + diis_max_iter : int, optional (default=200) + Maximal number of DIIS iterations + """ + if callback is None: + def callback(state, identifier): + pass + if not isinstance(matrix, AdcMatrixlike): + raise TypeError("matrix is not of type AdcMatrixlike") + for guess in guesses: + if not isinstance(guess, AmplitudeVector): + raise TypeError("One of the guesses is not of type AmplitudeVector") + + if explicit_symmetrisation is not None and \ + isinstance(explicit_symmetrisation, type): + explicit_symmetrisation = explicit_symmetrisation(matrix) + + if n_ep is None: + n_ep = len(guesses) + elif n_ep > len(guesses): + raise ValueError("n_ep cannot exceed the number of guess vectors.") + if not max_subspace: + # TODO Arnoldi uses this: + # max_subspace = max(2 * n_ep + 1, 20) + max_subspace = max(6 * n_ep, 20, 5 * len(guesses)) + + def convergence_test(state): + state.residuals_converged = state.residual_norms < conv_tol + state.converged = np.all(state.residuals_converged) + return state.converged + + def convergence_micro(state): # really rough + state.converged = np.abs(state.history_rval[0] + - state.history_rval[1]) > state.energy_diff + return state.converged + + def residualNorm_folded(state, diis_omegaUpdate=False): + state.eigenvector /= np.sqrt(state.eigenvector @ state.eigenvector) + Av = folded_matrix @ state.eigenvector + state.n_applies += 1 + # residual: r_i = A*v_i - w_i*v_i + state.residual = lincomb([1, -folded_matrix.omega], + [Av, state.eigenvector], evaluate=True) + state.residual_norm = np.sqrt(state.residual @ state.residual) + # For DIIS, update the eigenvalue corresponding to the new eigenvector. + if diis_omegaUpdate: + state.eigenvalues[state.n_state] = Av @ state.eigenvector + folded_matrix.update_omega(state.eigenvalues[state.n_state]) + return state.residual_norm + + if conv_tol < matrix.shape[1] * np.finfo(float).eps: + warnings.warn(la.LinAlgWarning( + "Convergence tolerance (== {:5.2g}) lower than " + "estimated maximal numerical accuracy (== {:5.2g}). " + "Convergence might be hard to achieve." + "".format(conv_tol, matrix.shape[1] * np.finfo(float).eps) + )) + + state = DavidsonState(matrix, guesses) + state.eigenvalues = np.empty(n_ep) + state.eigenvectors = [] + state.residual_norms = np.empty(n_ep) + + folded_matrix = AdcMatrixFolded(matrix) + if preconditioner is not None and isinstance(preconditioner, type): + preconditioner = preconditioner(matrix) + preconditioner.diagonal = folded_matrix.diagonal() + + # Retain single part of guess vectors + guesses_i = [AmplitudeVector(ph=guess.__getitem__("ph")) for guess in guesses] + if omegas is None: + # Calculate the initial (guess) eigenvalue for state 0. + Avi = matrix.block_apply("ph_ph", guesses_i[0].ph) + state.eigenvalues[0] = Avi.dot(guesses_i[0].ph) + + state.timer.restart("folded iterations") + for n_state in range(n_ep): + # Initialize guess omega for excited states. + if omegas is None: + folded_matrix.update_omega(state.eigenvalues[n_state]) + else: + folded_matrix.update_omega(omegas[n_state]) + + state_i = FoldedDavidsonState(folded_matrix, guesses_i) + state_i.n_state = n_state + callback(state_i, "folded_start") + + # Macro iterations for state i + state_i.timer.restart("folded iterations") + while state_i.macro_iter < macro_max_iter: + state_i.macro_iter += 1 + # Micro davidson iteration for diagonalising A(w_i) + state_i = davidson_iterations( + folded_matrix, + state_i, + max_subspace, + max_iter, + n_ep=n_ep, + is_converged=convergence_micro, + callback=callback, + which=which, + preconditioner=preconditioner, + preconditioning_method=preconditioning_method, + debug_checks=debug_checks, + residual_min_norm=residual_min_norm, + explicit_symmetrisation=explicit_symmetrisation) + + state.n_iter += state_i.n_iter + # Update omega and calculate the residual_norm + # under the latest omega for state i. + folded_matrix.update_omega(state_i.eigenvalues[state_i.n_state]) + residualNorm_folded(state_i) + callback(state_i, "micro") + if state_i.residual_norm < macro_conv_tol: + state_i.converged_macro = True + break + if state_i.macro_iter == macro_max_iter: + warnings.warn(la.LinAlgWarning( + f"Maximum number of macro iterations ({macro_max_iter}) " + "reached in modified davidson procedure.")) + + callback(state_i, "macro_stop") + # DIIS to further converge + state_i.timer.restart("folded iterations") + diis = DIIS(num_diis_vecs=num_diis_vecs, start_iter=4) + if not state_i.converged_macro: + warnings.warn(la.LinAlgWarning( + "Macro iterations with Davidson diagonalization " + "is not converged yet.")) + + preconditioner.update_shifts(float(0)) + while diis.iter_idx < diis_max_iter: + b_i = state_i.eigenvector + preconditioner @ state_i.residual + # corrected vector: b_i = u_i + residual_i / D11 + state_i.eigenvector = diis.compute_new_vec(b_i, state_i.residual) + residualNorm_folded(state_i, diis_omegaUpdate=True) + callback(state_i, "DIIS_steps") + if state_i.residual_norm < conv_tol: + state_i.converged_diis = True + break + if diis.iter_idx == diis_max_iter: + warnings.warn(la.LinAlgWarning( + f"Maximum number of iterations (== {diis_max_iter}) " + "reached in DIIS procedure.")) + state_i.DIIS_iter = diis.iter_idx + state.DIIS_iter += state_i.DIIS_iter + callback(state_i, "DIIS_stop") + + guesses_i = state_i.subspace_vectors.copy() + # Orthonormalize and update guesses_i for the next state: + # for state 0, taking initial guesses as guess vectors; + # for the higher-excited states, taking all eigenvectors of + # current A(w_i) as guess vectors. + guess_vecs = guesses_i.copy() + del guess_vecs[n_state] + coefficient = np.hstack(([1], -(state_i.eigenvector @ guess_vecs))) + newVec = lincomb(coefficient, [state_i.eigenvector] + + guess_vecs, evaluate=True) + state_i.eigenvector = newVec / np.sqrt(newVec @ newVec) + guesses_i[n_state] = state_i.eigenvector + + # Collect results into the "DavidsonState" + state.n_applies += state_i.n_applies + state.eigenvalues[n_state:] = state_i.eigenvalues[n_state:] + state.residual_norms[n_state] = state_i.residual_norm + state_i.eigenvector = folded_matrix.unfold(state_i.eigenvector) + state_i.eigenvector /= np.sqrt(state_i.eigenvector @ state_i.eigenvector) + state.eigenvectors.append(state_i.eigenvector) + + if convergence_test(state): + callback(state, "sum_folded") + state.timer.stop("folded iterations") + return state + + def jacobi_davidson(*args, **kwargs): return eigsh(*args, preconditioner=JacobiPreconditioner, preconditioning_method="Davidson", **kwargs) @@ -409,3 +727,76 @@ def jacobi_davidson(*args, **kwargs): def davidson(*args, **kwargs): return eigsh(*args, preconditioner=None, **kwargs) + + +def davidson_folded_DIIS(*args, **kwargs): + return eigsh_folded(*args, preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson", **kwargs) + + +class DIIS: + """ + An implementation of DIIS acceleration, adapted from + https://github.com/edeprince3/pdaggerq/blob/master/examples/full_cc_codes/diis.py + """ + def __init__(self, num_diis_vecs: int, start_iter=4): + """ + Initialize DIIS updater + + :params num_diis_vecs: Integer number representing number of DIIS + vectors to keep + :param start_iter: optional (default=4) number to start DIIS iterations + """ + self.nvecs = num_diis_vecs + self.error_vecs = [] + self.prev_vecs = [] + self.start_iter = start_iter + self.iter_idx = 0 + + def compute_new_vec(self, iterate, error): + """ + Compute a DIIS update. Only perform diis update after start_iter + have been accumulated. + """ + # don't start DIIS until start_iter + if self.iter_idx < self.start_iter: + self.iter_idx += 1 + return iterate + + # add iterate and error to the list of error and iterates + self.prev_vecs.append(iterate) + self.error_vecs.append(error) + self.iter_idx += 1 + + # if prev_vecs is larger than the diis space size, then pop the oldest + if len(self.prev_vecs) > self.nvecs: + self.prev_vecs.pop(0) + self.error_vecs.pop(0) + + # construct bmat and solve ax=b diis problem + b_mat, rhs = self.get_bmat() + c = np.linalg.solve(b_mat, rhs) + c = c.flatten() + + # construct new iterate from solution to diis ax=b and previous vecs. + new_iterate = self.prev_vecs[0].zeros_like() + for ii in range(len(self.prev_vecs)): + new_iterate += c[ii] * self.prev_vecs[ii] + return new_iterate + + def get_bmat(self): + """ + Compute b-mat + """ + dim = len(self.prev_vecs) + b = np.zeros((dim, dim)) + for i, j in product(range(dim), repeat=2): + if i <= j: + b[i, j] = self.error_vecs[i].dot(self.error_vecs[j]) + b[j, i] = b[i, j] + b = np.hstack((b, -1 * np.ones((dim, 1)))) + b = np.vstack((b, -1 * np.ones((1, dim + 1)))) + b[-1, -1] = 0 + rhs = np.zeros((dim + 1, 1)) + rhs[-1, 0] = -1 + return b, rhs diff --git a/adcc/solver/test_davidson_folded.py b/adcc/solver/test_davidson_folded.py new file mode 100644 index 000000000..492e83db9 --- /dev/null +++ b/adcc/solver/test_davidson_folded.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +## --------------------------------------------------------------------- +## +## Copyright (C) 2018 by the adcc authors +## +## This file is part of adcc. +## +## adcc is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published +## by the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## adcc is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with adcc. If not, see . +## +## --------------------------------------------------------------------- +import adcc +import unittest + +from pytest import approx + +from adcc import LazyMp +from adcc.testdata.cache import cache +from adcc.solver.davidson import davidson_folded_DIIS + + +class TestSolverDavidson(unittest.TestCase): + def test_adc2_singlets_folded(self): + refdata = cache.reference_data["h2o_sto3g"] + matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) + + # Solve for singlets + guesses = adcc.guesses_singlet(matrix, n_guesses=8, block="ph") + res = davidson_folded_DIIS(matrix, guesses, n_ep=8) + + ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_singlets[:8]) + + def test_adc2_triplets_folded(self): + refdata = cache.reference_data["h2o_sto3g"] + matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) + + # Solve for triplets + guesses = adcc.guesses_triplet(matrix, n_guesses=8, block="ph") + res = davidson_folded_DIIS(matrix, guesses, n_ep=8) + + ref_triplets = refdata["adc2"]["triplet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_triplets[:8]) + + def test_adc2_singlets_folded_adc1Guesses(self): + from adcc.workflow import run_adc + + refdata = cache.reference_data["h2o_sto3g"] + matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) + + # run adc1 for initial guesses + matrix_adc1 = adcc.AdcMatrix("adc1", LazyMp(cache.refstate["h2o_sto3g"])) + adc1 = run_adc(matrix_adc1, method="adc1", n_singlets=8) + omegas = adc1.excitation_energy_uncorrected + guesses = adc1.excitation_vector + # Solve for singlets + res = davidson_folded_DIIS(matrix, guesses, omegas=omegas, n_ep=8) + + ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_singlets[:8]) + + def test_adc2_triplets_folded_adc1(self): + from adcc.workflow import run_adc + + refdata = cache.reference_data["h2o_sto3g"] + matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) + + # run adc1 for initial guesses + matrix_adc1 = adcc.AdcMatrix("adc1", LazyMp(cache.refstate["h2o_sto3g"])) + adc1 = run_adc(matrix_adc1, method="adc1", n_triplets=8) + omegas = adc1.excitation_energy_uncorrected + guesses = adc1.excitation_vector + # Solve for triplets + res = davidson_folded_DIIS(matrix, guesses=guesses, omegas=omegas, n_ep=8) + + ref_triplets = refdata["adc2"]["triplet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_triplets[:8]) diff --git a/adcc/test_workflow.py b/adcc/test_workflow.py index 8cb0f3069..1799b22e2 100644 --- a/adcc/test_workflow.py +++ b/adcc/test_workflow.py @@ -215,6 +215,26 @@ def test_diagonalise_adcmatrix(self): assert res.converged assert res.eigenvalues == approx(ref_singlets[:3]) + guesses = adcc.guesses_singlet(matrix, n_guesses=6, block="ph") + res = diagonalise_adcmatrix(matrix, n_states=3, kind="singlet", + guesses=guesses, fold=True) + ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_singlets[:3]) + + from adcc.workflow import run_adc + matrix_adc1 = adcc.AdcMatrix("adc1", adcc.LazyMp + (cache.refstate["h2o_sto3g"])) + adc1 = run_adc(matrix_adc1, method="adc1", kind="singlet", n_states=3) + omegas = adc1.excitation_energy_uncorrected + guesses = adc1.excitation_vector + res = diagonalise_adcmatrix(matrix, n_states=3, kind="singlet", + fold=True, guesses_fold="adc1", + guesses=guesses, omegas=omegas) + ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] + assert res.converged + assert res.eigenvalues == approx(ref_singlets[:3]) + with pytest.raises(InputError): # Too low tolerance res = diagonalise_adcmatrix(matrix, n_states=9, kind="singlet", eigensolver="davidson", diff --git a/adcc/workflow.py b/adcc/workflow.py index a290032e7..12f7f77e2 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -35,7 +35,7 @@ from .ExcitedStates import ExcitedStates from .ReferenceState import ReferenceState as adcc_ReferenceState from .solver.lanczos import lanczos -from .solver.davidson import jacobi_davidson +from .solver.davidson import jacobi_davidson, davidson_folded_DIIS from .solver.explicit_symmetrisation import (IndexSpinSymmetrisation, IndexSymmetrisation) @@ -47,7 +47,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, n_guesses_doubles=None, output=sys.stdout, core_orbitals=None, frozen_core=None, frozen_virtual=None, method=None, n_singlets=None, n_triplets=None, n_spin_flip=None, - environment=None, **solverargs): + environment=None, fold=None, guesses_fold=None, + omegas=None, **solverargs): """Run an ADC calculation. Main entry point to run an ADC calculation. The reference to build the ADC @@ -133,12 +134,29 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, The keywords to specify how coupling to an environment model, e.g. PE, is treated. For details see :ref:`environment`. + fold : bool, optional + Perfom modified Davidson algorithm when using doubles-folded ADC(2) matrix. + + guesses_fold: str, optional + There is one more option for the initial guess vectors of the modified + Davidson solver, i.e., "adc1" takes initial guess vectors and guess + eigenvalues as eigenvalues and eigenvectors of adc(1) matrix. + Other parameters ---------------- max_subspace : int, optional Maximal subspace size max_iter : int, optional Maximal number of iterations + macro_conv_tol : float, optional + Convergence tolerance on the l2 norm squared of residuals to consider + them converged during macro iterations. + macro_max_iter : int, optional + Maximal number of macro iterations + num_diis_vecs: int, optional (default=50) + Maximal number of DIIS vectors to keep + diis_max_iter : int, optional + Maximal number of DIIS iterations Returns ------- @@ -202,10 +220,17 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if env_matrix_term: matrix += env_matrix_term + if fold and guesses_fold == "adc1" and guesses is None: + fold_initial = run_adc(data_or_matrix, method="adc1", + kind=kind, n_states=n_states) + omegas = fold_initial.excitation_energy_uncorrected + guesses = fold_initial.excitation_vector + diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, n_guesses_doubles=n_guesses_doubles, conv_tol=conv_tol, output=output, - eigensolver=eigensolver, **solverargs) + eigensolver=eigensolver, fold=fold, guesses_fold=guesses_fold, + omegas=omegas, **solverargs) exstates = ExcitedStates(diagres) exstates.kind = kind exstates.spin_change = spin_change @@ -347,7 +372,8 @@ def validate_state_parameters(reference_state, n_states=None, n_singlets=None, def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", guesses=None, n_guesses=None, n_guesses_doubles=None, - conv_tol=None, output=sys.stdout, **solverargs): + conv_tol=None, output=sys.stdout, fold=None, + guesses_fold=None, omegas=None, **solverargs): """ This function seeks appropriate guesses and afterwards proceeds to diagonalise the ADC matrix using the specified eigensolver. @@ -375,10 +401,16 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", # Set some solver-specific parameters if eigensolver == "davidson": n_guesses_per_state = 2 - callback = setup_solver_printing( - "Jacobi-Davidson", matrix, kind, solver.davidson.default_print, - output=output) - run_eigensolver = jacobi_davidson + if fold: + callback = setup_solver_printing( + "Jacobi-Davidson using doubles-folding", matrix, kind, + solver.davidson.default_print, output=output) + run_eigensolver = davidson_folded_DIIS + else: + callback = setup_solver_printing( + "Jacobi-Davidson", matrix, kind, solver.davidson.default_print, + output=output) + run_eigensolver = jacobi_davidson elif eigensolver == "lanczos": n_guesses_per_state = 1 callback = setup_solver_printing( @@ -407,10 +439,16 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", "are explicitly provided.") solverargs.setdefault("which", "SA") - return run_eigensolver(matrix, guesses, n_ep=n_states, conv_tol=conv_tol, - callback=callback, - explicit_symmetrisation=explicit_symmetrisation, - **solverargs) + if fold and guesses_fold == "adc1": + return run_eigensolver(matrix, guesses=guesses, omegas=omegas, + n_ep=n_states, conv_tol=conv_tol, callback=callback, + explicit_symmetrisation=explicit_symmetrisation, + **solverargs) + else: + return run_eigensolver(matrix, guesses, n_ep=n_states, conv_tol=conv_tol, + callback=callback, + explicit_symmetrisation=explicit_symmetrisation, + **solverargs) def estimate_n_guesses(matrix, n_states, singles_only=True, From 4906cff365b6532323fb19b385b372906c910ca3 Mon Sep 17 00:00:00 2001 From: TianLin-228 Date: Tue, 8 Apr 2025 17:17:33 +0200 Subject: [PATCH 2/4] Merge branch 'master' into davidson_foldingMatrix --- adcc/solver/davidson.py | 21 ++++++- adcc/solver/test_davidson_folded.py | 92 ----------------------------- adcc/tests/solver/davidson_test.py | 78 +++++++++++++++++++++++- 3 files changed, 96 insertions(+), 95 deletions(-) delete mode 100644 adcc/solver/test_davidson_folded.py diff --git a/adcc/solver/davidson.py b/adcc/solver/davidson.py index 23abc9ecf..1010afbfd 100644 --- a/adcc/solver/davidson.py +++ b/adcc/solver/davidson.py @@ -555,7 +555,8 @@ def convergence_test(state): return state -def eigsh_folded(matrix, guesses, omegas=None, n_ep=None, max_subspace=None, +def eigsh_folded(matrix, guesses, omegas=None, n_ep=None, n_block=None, + max_subspace=None, conv_tol=1e-9, which="SA", max_iter=70, callback=None, preconditioner=None, preconditioning_method="Davidson", @@ -573,6 +574,9 @@ def eigsh_folded(matrix, guesses, omegas=None, n_ep=None, max_subspace=None, Guess vectors (fixes also the Davidson block size) n_ep : int or NoneType, optional Number of eigenpairs to be computed + n_block : int or NoneType, optional + The solver block size: the number of vectors that are added to the subspace + in each iteration max_subspace : int or NoneType, optional Maximal subspace size conv_tol : float, optional @@ -626,7 +630,18 @@ def callback(state, identifier): if n_ep is None: n_ep = len(guesses) elif n_ep > len(guesses): - raise ValueError("n_ep cannot exceed the number of guess vectors.") + raise ValueError(f"n_ep (= {n_ep}) cannot exceed the number of guess " + f"vectors (= {len(guesses)}).") + + if n_block is None: + n_block = n_ep + elif n_block < n_ep: + raise ValueError(f"n_block (= {n_block}) cannot be smaller than the number " + f"of states requested (= {n_ep}).") + elif n_block > len(guesses): + raise ValueError(f"n_block (= {n_block}) cannot exceed the number of guess " + f"vectors (= {len(guesses)}).") + if not max_subspace: # TODO Arnoldi uses this: # max_subspace = max(2 * n_ep + 1, 20) @@ -676,6 +691,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): # Retain single part of guess vectors guesses_i = [AmplitudeVector(ph=guess.__getitem__("ph")) for guess in guesses] + if omegas is None: # Calculate the initial (guess) eigenvalue for state 0. Avi = matrix.block_apply("ph_ph", guesses_i[0].ph) @@ -704,6 +720,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): max_subspace, max_iter, n_ep=n_ep, + n_block=n_block, is_converged=convergence_micro, callback=callback, which=which, diff --git a/adcc/solver/test_davidson_folded.py b/adcc/solver/test_davidson_folded.py deleted file mode 100644 index 492e83db9..000000000 --- a/adcc/solver/test_davidson_folded.py +++ /dev/null @@ -1,92 +0,0 @@ -#!/usr/bin/env python3 -## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab -## --------------------------------------------------------------------- -## -## Copyright (C) 2018 by the adcc authors -## -## This file is part of adcc. -## -## adcc is free software: you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published -## by the Free Software Foundation, either version 3 of the License, or -## (at your option) any later version. -## -## adcc is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with adcc. If not, see . -## -## --------------------------------------------------------------------- -import adcc -import unittest - -from pytest import approx - -from adcc import LazyMp -from adcc.testdata.cache import cache -from adcc.solver.davidson import davidson_folded_DIIS - - -class TestSolverDavidson(unittest.TestCase): - def test_adc2_singlets_folded(self): - refdata = cache.reference_data["h2o_sto3g"] - matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) - - # Solve for singlets - guesses = adcc.guesses_singlet(matrix, n_guesses=8, block="ph") - res = davidson_folded_DIIS(matrix, guesses, n_ep=8) - - ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] - assert res.converged - assert res.eigenvalues == approx(ref_singlets[:8]) - - def test_adc2_triplets_folded(self): - refdata = cache.reference_data["h2o_sto3g"] - matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) - - # Solve for triplets - guesses = adcc.guesses_triplet(matrix, n_guesses=8, block="ph") - res = davidson_folded_DIIS(matrix, guesses, n_ep=8) - - ref_triplets = refdata["adc2"]["triplet"]["eigenvalues"] - assert res.converged - assert res.eigenvalues == approx(ref_triplets[:8]) - - def test_adc2_singlets_folded_adc1Guesses(self): - from adcc.workflow import run_adc - - refdata = cache.reference_data["h2o_sto3g"] - matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) - - # run adc1 for initial guesses - matrix_adc1 = adcc.AdcMatrix("adc1", LazyMp(cache.refstate["h2o_sto3g"])) - adc1 = run_adc(matrix_adc1, method="adc1", n_singlets=8) - omegas = adc1.excitation_energy_uncorrected - guesses = adc1.excitation_vector - # Solve for singlets - res = davidson_folded_DIIS(matrix, guesses, omegas=omegas, n_ep=8) - - ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] - assert res.converged - assert res.eigenvalues == approx(ref_singlets[:8]) - - def test_adc2_triplets_folded_adc1(self): - from adcc.workflow import run_adc - - refdata = cache.reference_data["h2o_sto3g"] - matrix = adcc.AdcMatrix("adc2", LazyMp(cache.refstate["h2o_sto3g"])) - - # run adc1 for initial guesses - matrix_adc1 = adcc.AdcMatrix("adc1", LazyMp(cache.refstate["h2o_sto3g"])) - adc1 = run_adc(matrix_adc1, method="adc1", n_triplets=8) - omegas = adc1.excitation_energy_uncorrected - guesses = adc1.excitation_vector - # Solve for triplets - res = davidson_folded_DIIS(matrix, guesses=guesses, omegas=omegas, n_ep=8) - - ref_triplets = refdata["adc2"]["triplet"]["eigenvalues"] - assert res.converged - assert res.eigenvalues == approx(ref_triplets[:8]) diff --git a/adcc/tests/solver/davidson_test.py b/adcc/tests/solver/davidson_test.py index 2182d49d8..5c4e70866 100644 --- a/adcc/tests/solver/davidson_test.py +++ b/adcc/tests/solver/davidson_test.py @@ -25,7 +25,7 @@ import pytest from adcc import LazyMp -from adcc.solver.davidson import jacobi_davidson, eigsh +from adcc.solver.davidson import jacobi_davidson, eigsh, davidson_folded_DIIS from adcc.misc import cached_property from ..testdata_cache import testdata_cache @@ -127,3 +127,79 @@ def test_adc2_triplets(self): assert n_states > 1 assert res.converged assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states]) + + def test_adc2_singlets_folded(self): + refdata = testdata_cache.adcman_data( + system="h2o_sto3g", method="adc2", case="gen" + )["singlet"] + + # Solve for singlets + guesses = adcc.guesses_singlet(self.matrix, n_guesses=8, block="ph") + res = davidson_folded_DIIS(self.matrix, guesses, n_ep=8) + + ref_singlets = refdata["eigenvalues"] + n_states = min(len(ref_singlets), len(res.eigenvalues)) + assert n_states > 1 + assert res.converged + assert res.eigenvalues[:n_states] == pytest.approx(ref_singlets[:n_states]) + + def test_adc2_triplets_folded(self): + refdata = testdata_cache.adcman_data( + system="h2o_sto3g", method="adc2", case="gen" + )["triplet"] + # Solve for triplets + guesses = adcc.guesses_triplet(self.matrix, n_guesses=8, block="ph") + res = davidson_folded_DIIS(self.matrix, guesses, n_ep=8) + + ref_triplets = refdata["eigenvalues"] + n_states = min(len(ref_triplets), len(res.eigenvalues)) + assert n_states > 1 + assert res.converged + assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states]) + + def test_adc2_singlets_folded_adc1Guesses(self): + from adcc.workflow import run_adc + + refdata = testdata_cache.adcman_data( + system="h2o_sto3g", method="adc2", case="gen" + )["singlet"] + + # run adc1 as initial guesses + matrix_adc1 = adcc.AdcMatrix( + "adc1", LazyMp(testdata_cache.refstate("h2o_sto3g", case="gen")) + ) + adc1 = run_adc(matrix_adc1, method="adc1", n_singlets=8) + omegas = adc1.excitation_energy_uncorrected + guesses = adc1.excitation_vector + # Solve for singlets + res = davidson_folded_DIIS(self.matrix, guesses, omegas=omegas, n_ep=8) + + ref_singlets = refdata["eigenvalues"] + n_states = min(len(ref_singlets), len(res.eigenvalues)) + assert n_states > 1 + assert res.converged + assert res.eigenvalues[:n_states] == pytest.approx(ref_singlets[:n_states]) + + def test_adc2_triplets_folded_adc1Guesses(self): + from adcc.workflow import run_adc + + refdata = testdata_cache.adcman_data( + system="h2o_sto3g", method="adc2", case="gen" + )["triplet"] + + # run adc1 as initial guesses + matrix_adc1 = adcc.AdcMatrix( + "adc1", LazyMp(testdata_cache.refstate("h2o_sto3g", case="gen")) + ) + adc1 = run_adc(matrix_adc1, method="adc1", n_triplets=8) + omegas = adc1.excitation_energy_uncorrected + guesses = adc1.excitation_vector + # Solve for triplets + res = davidson_folded_DIIS(self.matrix, guesses=guesses, + omegas=omegas, n_ep=8) + + ref_triplets = refdata["eigenvalues"] + n_states = min(len(ref_triplets), len(res.eigenvalues)) + assert n_states > 1 + assert res.converged + assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states]) From bcc13f9a70c6b294fdc474ed922b7e8efa84c3af Mon Sep 17 00:00:00 2001 From: TianLin-228 Date: Tue, 15 Apr 2025 13:16:33 +0200 Subject: [PATCH 3/4] delete update_omega & test for AdcMatrixFolded --- adcc/AdcMatrix.py | 27 +++----- adcc/solver/davidson.py | 12 ++-- adcc/tests/AdcMatrix_test.py | 100 ++++++++++++++++++++++++++++- adcc/tests/solver/davidson_test.py | 1 + adcc/tests/workflow_test.py | 16 ++--- 5 files changed, 122 insertions(+), 34 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 14c76de83..2df69d4a3 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -666,9 +666,6 @@ def __init__(self, matrix): block_orders=matrix.block_orders, intermediates=matrix.intermediates) - def update_omega(self, omega): - self.omega = omega - def diagonal(self): """ Return the approximate diagonal of the ADC(2) matrix with doubles-folding. @@ -686,19 +683,11 @@ def matvec(self, other): return AmplitudeVector(ph=self.block_apply("ph_ph", other.ph) + self.block_apply("ph_pphh", u2)) - def __matmul__(self, other): - if isinstance(other, AmplitudeVector): - return self.matvec(other) - if isinstance(other, list): - if all(isinstance(elem, AmplitudeVector) for elem in other): - return [self.matvec(v) for v in other] - return NotImplemented - - def unfold(self, u1): - """ - recompute the doubles component and return the complete vector. - """ - diag = super().diagonal().pphh - e = diag.ones_like() - u2 = self.block_apply("pphh_ph", u1.ph) / (e * self.omega - diag) - return AmplitudeVector(ph=u1.ph, pphh=u2) + # def unfold(self, u1): + # """ + # recompute the doubles component and return the complete vector. + # """ + # diag = super().diagonal().pphh + # e = diag.ones_like() + # u2 = self.block_apply("pphh_ph", u1.ph) / (e * self.omega - diag) + # return AmplitudeVector(ph=u1.ph, pphh=u2) \ No newline at end of file diff --git a/adcc/solver/davidson.py b/adcc/solver/davidson.py index 1010afbfd..7dce0ad66 100644 --- a/adcc/solver/davidson.py +++ b/adcc/solver/davidson.py @@ -668,7 +668,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): # For DIIS, update the eigenvalue corresponding to the new eigenvector. if diis_omegaUpdate: state.eigenvalues[state.n_state] = Av @ state.eigenvector - folded_matrix.update_omega(state.eigenvalues[state.n_state]) + folded_matrix.omega = state.eigenvalues[state.n_state] return state.residual_norm if conv_tol < matrix.shape[1] * np.finfo(float).eps: @@ -691,7 +691,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): # Retain single part of guess vectors guesses_i = [AmplitudeVector(ph=guess.__getitem__("ph")) for guess in guesses] - + if omegas is None: # Calculate the initial (guess) eigenvalue for state 0. Avi = matrix.block_apply("ph_ph", guesses_i[0].ph) @@ -701,9 +701,9 @@ def residualNorm_folded(state, diis_omegaUpdate=False): for n_state in range(n_ep): # Initialize guess omega for excited states. if omegas is None: - folded_matrix.update_omega(state.eigenvalues[n_state]) + folded_matrix.omega = state.eigenvalues[n_state] else: - folded_matrix.update_omega(omegas[n_state]) + folded_matrix.omega = omegas[n_state] state_i = FoldedDavidsonState(folded_matrix, guesses_i) state_i.n_state = n_state @@ -733,7 +733,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): state.n_iter += state_i.n_iter # Update omega and calculate the residual_norm # under the latest omega for state i. - folded_matrix.update_omega(state_i.eigenvalues[state_i.n_state]) + folded_matrix.omega = state_i.eigenvalues[state_i.n_state] residualNorm_folded(state_i) callback(state_i, "micro") if state_i.residual_norm < macro_conv_tol: @@ -788,7 +788,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): state.n_applies += state_i.n_applies state.eigenvalues[n_state:] = state_i.eigenvalues[n_state:] state.residual_norms[n_state] = state_i.residual_norm - state_i.eigenvector = folded_matrix.unfold(state_i.eigenvector) + # state_i.eigenvector = folded_matrix.unfold(state_i.eigenvector) state_i.eigenvector /= np.sqrt(state_i.eigenvector @ state_i.eigenvector) state.eigenvectors.append(state_i.eigenvector) diff --git a/adcc/tests/AdcMatrix_test.py b/adcc/tests/AdcMatrix_test.py index d3ec88bdf..146b92293 100644 --- a/adcc/tests/AdcMatrix_test.py +++ b/adcc/tests/AdcMatrix_test.py @@ -25,7 +25,11 @@ import numpy as np from numpy.testing import assert_allclose -from adcc.AdcMatrix import AdcExtraTerm, AdcMatrixProjected, AdcMatrixShifted +from adcc.AdcMatrix import ( + AdcExtraTerm, + AdcMatrixProjected, + AdcMatrixShifted, + AdcMatrixFolded) from adcc.Intermediates import Intermediates from adcc.adc_pp.matrix import AdcBlock @@ -401,3 +405,97 @@ def test_matmul(self, system): assert_nonzero_blocks(ores.pphh, pres.pphh, nonzeros["pphh"], tol=1e-14) # TODO Test block_view, block_apply + + +@pytest.mark.parametrize("system", ["h2o_sto3g", "cn_sto3g"]) +class TestAdcMatrixFolded: + def construct_matrices(self, system: str): + reference_state = testdata_cache.refstate(system, case="gen") + ground_state = adcc.LazyMp(reference_state) + matrix = adcc.AdcMatrix("adc2", ground_state) + + folded = AdcMatrixFolded(matrix) + return matrix, folded + + def test_diagonal(self, system: str): + matrix, folded = self.construct_matrices(system) + + odiag = matrix.diagonal() + fdiag = folded.diagonal() + assert np.max(np.abs(fdiag["ph"].to_ndarray() + - odiag["ph"].to_ndarray())) < 1e-12 + assert not hasattr(fdiag, "pphh") + + def test_matmul(self, system: str): + omega = 0.5 + matrix, folded = self.construct_matrices(system) + folded.omega = omega + + vec = adcc.guess_zero(matrix) + vec.set_random() + + ores = matrix @ vec + fres = folded @ vec + + assert ores.ph.describe_symmetry() == fres.ph.describe_symmetry() + assert not hasattr(fres, "pphh") + # assert_equal_symmetry(ores.ph, fres.ph) + + # Consistency of the different part of ores.ph and freq.sh, i.e A_12 * V2 + diff_o = matrix.block_apply("ph_pphh", vec.pphh) + + diag = matrix.diagonal().pphh + e = diag.ones_like() + u2 = matrix.block_apply("pphh_ph", vec.ph) / (e * folded.omega - diag) + diff_f = matrix.block_apply("ph_pphh", u2) + + diff_matmul = ores.ph - fres.ph + diff_man = diff_o - diff_f + diff = diff_man - diff_matmul + assert np.max(np.abs(diff.to_ndarray())) < 1e-12 + + +""" + def test_unfold(self, system: str): + + matrix, folded = self.construct_matrices(system) + + # omega = 0.5 + + from adcc.workflow import run_adc + adc2 = run_adc(matrix, method="adc2", kind="singlet", n_states=1) + omega = adc2.excitation_energy_uncorrected + vec = adc2.excitation_vector[0] + + # vec = adcc.guess_zero(matrix) + # vec.set_random() + print(vec) + + # from adcc.solver.conjugate_gradient import conjugate_gradient + # matrix_shifted = AdcMatrixShifted(matrix, shift=-omega) + # rhs = vec.zeros_like() + # fres_unfold = conjugate_gradient(matrix_shifted, rhs=rhs, x0=vec) + # print(fres_unfold) + folded.omega = omega + fres_unfold = folded.unfold(vec) + + # # Consistency of manully recomputed vector and "unfold" + # diff = vec.pphh - fres_unfold.pphh + # assert np.max(np.abs(diff.to_ndarray())) < 1e-12 + + # Consistency of symmetry + assert vec.ph.describe_symmetry() == fres_unfold.ph.describe_symmetry() + print("vec.pphh", vec.pphh.describe_symmetry()) + print("fres_unfold.pphh", fres_unfold.pphh.describe_symmetry()) + assert vec.pphh.describe_symmetry() == fres_unfold.pphh.describe_symmetry() + # assert_equal_symmetry(vec.pphh, fres_unfold.pphh) + + + # diag = matrix.diagonal().pphh + # e = diag.ones_like() + # u2_man = matrix.block_apply("pphh_ph", vec.ph) / (e * folded.omega - diag) + # assert u2_man.describe_symmetry() == fres_unfold.pphh.describe_symmetry() + # # Consistency of manully recomputed vector and "unfold" + # diff = u2_man - fres_unfold.pphh + # assert np.max(np.abs(diff.to_ndarray())) < 1e-12 +""" diff --git a/adcc/tests/solver/davidson_test.py b/adcc/tests/solver/davidson_test.py index 5c4e70866..f6d353255 100644 --- a/adcc/tests/solver/davidson_test.py +++ b/adcc/tests/solver/davidson_test.py @@ -147,6 +147,7 @@ def test_adc2_triplets_folded(self): refdata = testdata_cache.adcman_data( system="h2o_sto3g", method="adc2", case="gen" )["triplet"] + # Solve for triplets guesses = adcc.guesses_triplet(self.matrix, n_guesses=8, block="ph") res = davidson_folded_DIIS(self.matrix, guesses, n_ep=8) diff --git a/adcc/tests/workflow_test.py b/adcc/tests/workflow_test.py index 45cf078e5..c78656457 100644 --- a/adcc/tests/workflow_test.py +++ b/adcc/tests/workflow_test.py @@ -220,25 +220,25 @@ def test_diagonalise_adcmatrix(self): assert res.converged assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) + # for ADC problems with doubles-folding guesses = adcc.guesses_singlet(matrix, n_guesses=6, block="ph") - res = diagonalise_adcmatrix(matrix, n_states=3, kind="singlet", + res = diagonalise_adcmatrix(matrix, n_states=3, kind=kind, guesses=guesses, fold=True) - ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] assert res.converged - assert res.eigenvalues == approx(ref_singlets[:3]) + assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) from adcc.workflow import run_adc - matrix_adc1 = adcc.AdcMatrix("adc1", adcc.LazyMp - (cache.refstate["h2o_sto3g"])) - adc1 = run_adc(matrix_adc1, method="adc1", kind="singlet", n_states=3) + matrix_adc1 = adcc.AdcMatrix( + "adc1", adcc.LazyMp(testdata_cache.refstate("h2o_sto3g", case="gen")) + ) + adc1 = run_adc(matrix_adc1, method="adc1", kind=kind, n_states=3) omegas = adc1.excitation_energy_uncorrected guesses = adc1.excitation_vector res = diagonalise_adcmatrix(matrix, n_states=3, kind="singlet", fold=True, guesses_fold="adc1", guesses=guesses, omegas=omegas) - ref_singlets = refdata["adc2"]["singlet"]["eigenvalues"] assert res.converged - assert res.eigenvalues == approx(ref_singlets[:3]) + assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) with pytest.raises(InputError): # Too low tolerance # SCF tolerance = 1e-14 currently From 1197c30a71f9ddaaa832af7841ab118d82904764 Mon Sep 17 00:00:00 2001 From: TianLin-228 Date: Mon, 27 Oct 2025 02:50:26 +0100 Subject: [PATCH 4/4] emplicit symmetry; workflow test --- adcc/AdcMatrix.py | 40 ++++--- adcc/solver/davidson.py | 168 ++++++++++------------------- adcc/tests/AdcMatrix_test.py | 123 +++++++++++---------- adcc/tests/solver/davidson_test.py | 73 ++++++++++++- adcc/tests/workflow_test.py | 15 ++- adcc/workflow.py | 54 +++++++--- 6 files changed, 266 insertions(+), 207 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 2df69d4a3..b04589a64 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -30,7 +30,7 @@ from .functions import ones_like from .Intermediates import Intermediates from .AmplitudeVector import AmplitudeVector - +# from .solver.explicit_symmetrisation import IndexSymmetrisation class AdcExtraTerm: def __init__(self, matrix, blocks): @@ -665,6 +665,11 @@ def __init__(self, matrix): super().__init__(matrix.method, matrix.ground_state, block_orders=matrix.block_orders, intermediates=matrix.intermediates) + if matrix.method.name != "adc2": + raise TypeError("The matrix should be ADC2 Matrix for doubles-folding.") + + from .solver.explicit_symmetrisation import IndexSymmetrisation + self.isymm = IndexSymmetrisation(matrix) def diagonal(self): """ @@ -672,22 +677,29 @@ def diagonal(self): """ return AmplitudeVector(ph=super().diagonal().ph) - def matvec(self, other): + def u2_fold(self, other): """ - Compute the doubles-folded matrix-vector product of the singles vector with - an effective ADC matrix which depends on the eigenvalue ω. + Return the doubles component of u in terms of its singles component. """ diag = super().diagonal().pphh e = diag.ones_like() u2 = self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag) + return self.isymm.symmetrise(AmplitudeVector(pphh=u2)) + + def matvec(self, other): + """ + Compute the doubles-folded matrix-vector product of the singles vector with + an effective ADC matrix which depends on the eigenvalue(omega). + """ + u2 = self.u2_fold(other) return AmplitudeVector(ph=self.block_apply("ph_ph", other.ph) - + self.block_apply("ph_pphh", u2)) - - # def unfold(self, u1): - # """ - # recompute the doubles component and return the complete vector. - # """ - # diag = super().diagonal().pphh - # e = diag.ones_like() - # u2 = self.block_apply("pphh_ph", u1.ph) / (e * self.omega - diag) - # return AmplitudeVector(ph=u1.ph, pphh=u2) \ No newline at end of file + + self.block_apply("ph_pphh", u2.pphh)) + + def unfold(self, other): + """ + Recompute the doubles component and return the complete vector. + """ + other_symm = self.isymm.symmetrise(other) + u2 = self.u2_fold(other_symm) + + return AmplitudeVector(ph=other_symm.ph, pphh=u2.pphh) \ No newline at end of file diff --git a/adcc/solver/davidson.py b/adcc/solver/davidson.py index 7dce0ad66..6f1fc4119 100644 --- a/adcc/solver/davidson.py +++ b/adcc/solver/davidson.py @@ -34,7 +34,7 @@ from .preconditioner import JacobiPreconditioner from .SolverStateBase import EigenSolverStateBase from .explicit_symmetrisation import IndexSymmetrisation -from itertools import product +from .DIIS import DIIS class DavidsonState(EigenSolverStateBase): @@ -289,6 +289,8 @@ def form_residual(rval, rvec): if i in epair_mask] callback(state, "is_converged") else: + print(state.n_iter, rvecs.shape, len(np.transpose(rvecs)[0]),len(SS)) + # update guesses vectors for next macro iteration state.subspace_vectors = [lincomb(v, SS, evaluate=True) for v in np.transpose(rvecs)] @@ -321,6 +323,7 @@ def form_residual(rval, rvec): if n_ss_vec + n_block > max_subspace: callback(state, "restart") with state.timer.record("projection"): + print("TRUE") # The addition of the preconditioned vectors goes beyond max. # subspace size => Collapse first, ie keep current Ritz vectors # as new subspace @@ -415,6 +418,7 @@ def form_residual(rval, rvec): else: # Compute all eigenvectors as guesses vectors # for the next macro iteration. + print(state.n_iter, rvecs.shape, len(np.transpose(rvecs)[0]),len(SS)) state.subspace_vectors = [lincomb(v, SS, evaluate=True) for v in np.transpose(rvecs)] assert len(state.subspace_vectors) == n_block @@ -572,6 +576,10 @@ def eigsh_folded(matrix, guesses, omegas=None, n_ep=None, n_block=None, ADC(2) matrix instance guesses : list Guess vectors (fixes also the Davidson block size) + omegas : numpy.ndarray or NoneType, optional + Initial omega values for the modified Davidson solver when using + doubles-folded ADC(2) matrix. If not provided, the initial omega values + will be derived from the provided guess vectors. n_ep : int or NoneType, optional Number of eigenpairs to be computed n_block : int or NoneType, optional @@ -619,6 +627,9 @@ def callback(state, identifier): pass if not isinstance(matrix, AdcMatrixlike): raise TypeError("matrix is not of type AdcMatrixlike") + elif matrix.method.name != "adc2": + raise TypeError("matrix is not ADC2 Matrix") + for guess in guesses: if not isinstance(guess, AmplitudeVector): raise TypeError("One of the guesses is not of type AmplitudeVector") @@ -646,18 +657,24 @@ def callback(state, identifier): # TODO Arnoldi uses this: # max_subspace = max(2 * n_ep + 1, 20) max_subspace = max(6 * n_ep, 20, 5 * len(guesses)) + elif max_subspace < 2 * n_block: + raise ValueError(f"max_subspace (= {max_subspace}) needs to be at least " + f"twice as large as n_block (n_block = {n_block}).") + elif max_subspace < len(guesses): + raise ValueError(f"max_subspace (= {max_subspace}) cannot be smaller than " + f"the number of guess vectors (= {len(guesses)}).") - def convergence_test(state): + def convergence_test_final(state): state.residuals_converged = state.residual_norms < conv_tol state.converged = np.all(state.residuals_converged) return state.converged - def convergence_micro(state): # really rough + def convergence_test_micro(state): # really rough state.converged = np.abs(state.history_rval[0] - state.history_rval[1]) > state.energy_diff return state.converged - def residualNorm_folded(state, diis_omegaUpdate=False): + def residual_norm_folded(state, diis_omegaUpdate=False): state.eigenvector /= np.sqrt(state.eigenvector @ state.eigenvector) Av = folded_matrix @ state.eigenvector state.n_applies += 1 @@ -665,10 +682,9 @@ def residualNorm_folded(state, diis_omegaUpdate=False): state.residual = lincomb([1, -folded_matrix.omega], [Av, state.eigenvector], evaluate=True) state.residual_norm = np.sqrt(state.residual @ state.residual) - # For DIIS, update the eigenvalue corresponding to the new eigenvector. + # For DIIS, update the omega corresponding to the new eigenvector. if diis_omegaUpdate: - state.eigenvalues[state.n_state] = Av @ state.eigenvector - folded_matrix.omega = state.eigenvalues[state.n_state] + folded_matrix.omega = Av @ state.eigenvector return state.residual_norm if conv_tol < matrix.shape[1] * np.finfo(float).eps: @@ -686,24 +702,31 @@ def residualNorm_folded(state, diis_omegaUpdate=False): folded_matrix = AdcMatrixFolded(matrix) if preconditioner is not None and isinstance(preconditioner, type): - preconditioner = preconditioner(matrix) - preconditioner.diagonal = folded_matrix.diagonal() + preconditioner = preconditioner(folded_matrix) # Retain single part of guess vectors guesses_i = [AmplitudeVector(ph=guess.__getitem__("ph")) for guess in guesses] - + # If not given, the omega value of first state is computed + # from the guess vectors; the other omegas are set to the eigenvalues + # of the folded matrix under the previous state. + # If given, the guesses and omegas should be the pairs of eigenvalues and + # eigenvectors of the adc1 matrix. if omegas is None: - # Calculate the initial (guess) eigenvalue for state 0. + omegas = np.full(n_ep, np.nan) + # Calculate the initial (guess) omega for state 0. Avi = matrix.block_apply("ph_ph", guesses_i[0].ph) - state.eigenvalues[0] = Avi.dot(guesses_i[0].ph) + omegas[0] = Avi.dot(guesses_i[0].ph) + elif len(omegas) < n_ep: + raise ValueError(f"The number of omegas (={len(omegas)}) should be equal " + f"to the number of requested state(={n_ep}).") + print("OMEGAS2:", type(omegas),type(omegas[0])) state.timer.restart("folded iterations") for n_state in range(n_ep): - # Initialize guess omega for excited states. - if omegas is None: - folded_matrix.omega = state.eigenvalues[n_state] - else: - folded_matrix.omega = omegas[n_state] + # Initialize guess omega for the target state. + if np.isnan(omegas[n_state]): + omegas[n_state] = state_i.eigenvalues[n_state] + folded_matrix.omega = omegas[n_state] state_i = FoldedDavidsonState(folded_matrix, guesses_i) state_i.n_state = n_state @@ -713,28 +736,22 @@ def residualNorm_folded(state, diis_omegaUpdate=False): state_i.timer.restart("folded iterations") while state_i.macro_iter < macro_max_iter: state_i.macro_iter += 1 - # Micro davidson iteration for diagonalising A(w_i) - state_i = davidson_iterations( - folded_matrix, - state_i, - max_subspace, - max_iter, - n_ep=n_ep, - n_block=n_block, - is_converged=convergence_micro, - callback=callback, - which=which, - preconditioner=preconditioner, - preconditioning_method=preconditioning_method, - debug_checks=debug_checks, - residual_min_norm=residual_min_norm, - explicit_symmetrisation=explicit_symmetrisation) + # Micro davidson iterations for diagonalising A(w_i) + davidson_iterations(folded_matrix, state_i, max_subspace, + max_iter, n_ep=n_ep, n_block=n_block, + is_converged=convergence_test_micro, + callback=callback, which=which, + preconditioner=preconditioner, + preconditioning_method=preconditioning_method, + debug_checks=debug_checks, + residual_min_norm=residual_min_norm, + explicit_symmetrisation=explicit_symmetrisation) state.n_iter += state_i.n_iter # Update omega and calculate the residual_norm # under the latest omega for state i. folded_matrix.omega = state_i.eigenvalues[state_i.n_state] - residualNorm_folded(state_i) + residual_norm_folded(state_i) callback(state_i, "micro") if state_i.residual_norm < macro_conv_tol: state_i.converged_macro = True @@ -747,7 +764,7 @@ def residualNorm_folded(state, diis_omegaUpdate=False): callback(state_i, "macro_stop") # DIIS to further converge state_i.timer.restart("folded iterations") - diis = DIIS(num_diis_vecs=num_diis_vecs, start_iter=4) + diis = DIIS(num_diis_vecs=num_diis_vecs, start_iter=0) if not state_i.converged_macro: warnings.warn(la.LinAlgWarning( "Macro iterations with Davidson diagonalization " @@ -755,10 +772,10 @@ def residualNorm_folded(state, diis_omegaUpdate=False): preconditioner.update_shifts(float(0)) while diis.iter_idx < diis_max_iter: - b_i = state_i.eigenvector + preconditioner @ state_i.residual # corrected vector: b_i = u_i + residual_i / D11 + b_i = state_i.eigenvector + preconditioner @ state_i.residual state_i.eigenvector = diis.compute_new_vec(b_i, state_i.residual) - residualNorm_folded(state_i, diis_omegaUpdate=True) + residual_norm_folded(state_i, diis_omegaUpdate=True) callback(state_i, "DIIS_steps") if state_i.residual_norm < conv_tol: state_i.converged_diis = True @@ -769,13 +786,14 @@ def residualNorm_folded(state, diis_omegaUpdate=False): "reached in DIIS procedure.")) state_i.DIIS_iter = diis.iter_idx state.DIIS_iter += state_i.DIIS_iter + state_i.eigenvalues[n_state] = folded_matrix.omega callback(state_i, "DIIS_stop") - guesses_i = state_i.subspace_vectors.copy() # Orthonormalize and update guesses_i for the next state: # for state 0, taking initial guesses as guess vectors; # for the higher-excited states, taking all eigenvectors of # current A(w_i) as guess vectors. + guesses_i = state_i.subspace_vectors.copy() guess_vecs = guesses_i.copy() del guess_vecs[n_state] coefficient = np.hstack(([1], -(state_i.eigenvector @ guess_vecs))) @@ -786,13 +804,13 @@ def residualNorm_folded(state, diis_omegaUpdate=False): # Collect results into the "DavidsonState" state.n_applies += state_i.n_applies - state.eigenvalues[n_state:] = state_i.eigenvalues[n_state:] state.residual_norms[n_state] = state_i.residual_norm - # state_i.eigenvector = folded_matrix.unfold(state_i.eigenvector) + state.eigenvalues[n_state] = state_i.eigenvalues[n_state] + state_i.eigenvector = folded_matrix.unfold(state_i.eigenvector) state_i.eigenvector /= np.sqrt(state_i.eigenvector @ state_i.eigenvector) state.eigenvectors.append(state_i.eigenvector) - if convergence_test(state): + if convergence_test_final(state): callback(state, "sum_folded") state.timer.stop("folded iterations") return state @@ -810,71 +828,3 @@ def davidson(*args, **kwargs): def davidson_folded_DIIS(*args, **kwargs): return eigsh_folded(*args, preconditioner=JacobiPreconditioner, preconditioning_method="Davidson", **kwargs) - - -class DIIS: - """ - An implementation of DIIS acceleration, adapted from - https://github.com/edeprince3/pdaggerq/blob/master/examples/full_cc_codes/diis.py - """ - def __init__(self, num_diis_vecs: int, start_iter=4): - """ - Initialize DIIS updater - - :params num_diis_vecs: Integer number representing number of DIIS - vectors to keep - :param start_iter: optional (default=4) number to start DIIS iterations - """ - self.nvecs = num_diis_vecs - self.error_vecs = [] - self.prev_vecs = [] - self.start_iter = start_iter - self.iter_idx = 0 - - def compute_new_vec(self, iterate, error): - """ - Compute a DIIS update. Only perform diis update after start_iter - have been accumulated. - """ - # don't start DIIS until start_iter - if self.iter_idx < self.start_iter: - self.iter_idx += 1 - return iterate - - # add iterate and error to the list of error and iterates - self.prev_vecs.append(iterate) - self.error_vecs.append(error) - self.iter_idx += 1 - - # if prev_vecs is larger than the diis space size, then pop the oldest - if len(self.prev_vecs) > self.nvecs: - self.prev_vecs.pop(0) - self.error_vecs.pop(0) - - # construct bmat and solve ax=b diis problem - b_mat, rhs = self.get_bmat() - c = np.linalg.solve(b_mat, rhs) - c = c.flatten() - - # construct new iterate from solution to diis ax=b and previous vecs. - new_iterate = self.prev_vecs[0].zeros_like() - for ii in range(len(self.prev_vecs)): - new_iterate += c[ii] * self.prev_vecs[ii] - return new_iterate - - def get_bmat(self): - """ - Compute b-mat - """ - dim = len(self.prev_vecs) - b = np.zeros((dim, dim)) - for i, j in product(range(dim), repeat=2): - if i <= j: - b[i, j] = self.error_vecs[i].dot(self.error_vecs[j]) - b[j, i] = b[i, j] - b = np.hstack((b, -1 * np.ones((dim, 1)))) - b = np.vstack((b, -1 * np.ones((1, dim + 1)))) - b[-1, -1] = 0 - rhs = np.zeros((dim + 1, 1)) - rhs[-1, 0] = -1 - return b, rhs diff --git a/adcc/tests/AdcMatrix_test.py b/adcc/tests/AdcMatrix_test.py index 146b92293..d1cfbaf73 100644 --- a/adcc/tests/AdcMatrix_test.py +++ b/adcc/tests/AdcMatrix_test.py @@ -30,6 +30,7 @@ AdcMatrixProjected, AdcMatrixShifted, AdcMatrixFolded) +# from adcc.solver.AdcMatrixFolded import AdcMatrixFolded from adcc.Intermediates import Intermediates from adcc.adc_pp.matrix import AdcBlock @@ -414,88 +415,86 @@ def construct_matrices(self, system: str): ground_state = adcc.LazyMp(reference_state) matrix = adcc.AdcMatrix("adc2", ground_state) - folded = AdcMatrixFolded(matrix) - return matrix, folded + folded_matrix = AdcMatrixFolded(matrix) + return matrix, folded_matrix + + def test_adc2_matrix(self, system: str): + reference_state = testdata_cache.refstate(system, case="gen") + ground_state = adcc.LazyMp(reference_state) + matrix = adcc.AdcMatrix("adc3", ground_state) + with pytest.raises(TypeError): + AdcMatrixFolded(matrix) def test_diagonal(self, system: str): - matrix, folded = self.construct_matrices(system) + matrix, folded_matrix = self.construct_matrices(system) odiag = matrix.diagonal() - fdiag = folded.diagonal() + fdiag = folded_matrix.diagonal() assert np.max(np.abs(fdiag["ph"].to_ndarray() - odiag["ph"].to_ndarray())) < 1e-12 assert not hasattr(fdiag, "pphh") - def test_matmul(self, system: str): - omega = 0.5 - matrix, folded = self.construct_matrices(system) - folded.omega = omega + def test_u2_fold(self, system: str): + matrix, folded_matrix = self.construct_matrices(system) + folded_matrix.omega = 0.5 vec = adcc.guess_zero(matrix) vec.set_random() - ores = matrix @ vec - fres = folded @ vec - - assert ores.ph.describe_symmetry() == fres.ph.describe_symmetry() - assert not hasattr(fres, "pphh") - # assert_equal_symmetry(ores.ph, fres.ph) - - # Consistency of the different part of ores.ph and freq.sh, i.e A_12 * V2 - diff_o = matrix.block_apply("ph_pphh", vec.pphh) - diag = matrix.diagonal().pphh e = diag.ones_like() - u2 = matrix.block_apply("pphh_ph", vec.ph) / (e * folded.omega - diag) - diff_f = matrix.block_apply("ph_pphh", u2) - - diff_matmul = ores.ph - fres.ph - diff_man = diff_o - diff_f - diff = diff_man - diff_matmul - assert np.max(np.abs(diff.to_ndarray())) < 1e-12 - - -""" - def test_unfold(self, system: str): + u2_manual = matrix.block_apply("pphh_ph", vec.ph) / (e * folded_matrix.omega - diag) - matrix, folded = self.construct_matrices(system) + u2_fold = folded_matrix.u2_fold(vec) - # omega = 0.5 + diff = u2_manual - u2_fold.pphh + assert np.max(np.abs(diff.to_ndarray())) < 1e-12 - from adcc.workflow import run_adc - adc2 = run_adc(matrix, method="adc2", kind="singlet", n_states=1) - omega = adc2.excitation_energy_uncorrected - vec = adc2.excitation_vector[0] + def test_matvec(self, system: str): + matrix, folded_matrix = self.construct_matrices(system) + folded_matrix.omega = 0.5 - # vec = adcc.guess_zero(matrix) - # vec.set_random() - print(vec) + vec = adcc.guess_zero(matrix) + vec.set_random() - # from adcc.solver.conjugate_gradient import conjugate_gradient - # matrix_shifted = AdcMatrixShifted(matrix, shift=-omega) - # rhs = vec.zeros_like() - # fres_unfold = conjugate_gradient(matrix_shifted, rhs=rhs, x0=vec) - # print(fres_unfold) - folded.omega = omega - fres_unfold = folded.unfold(vec) + ores = matrix @ vec + fres = folded_matrix @ vec - # # Consistency of manully recomputed vector and "unfold" - # diff = vec.pphh - fres_unfold.pphh - # assert np.max(np.abs(diff.to_ndarray())) < 1e-12 + assert ores.ph.describe_symmetry() == fres.ph.describe_symmetry() + assert not hasattr(fres, "pphh") - # Consistency of symmetry - assert vec.ph.describe_symmetry() == fres_unfold.ph.describe_symmetry() - print("vec.pphh", vec.pphh.describe_symmetry()) - print("fres_unfold.pphh", fres_unfold.pphh.describe_symmetry()) - assert vec.pphh.describe_symmetry() == fres_unfold.pphh.describe_symmetry() - # assert_equal_symmetry(vec.pphh, fres_unfold.pphh) + # Consistency of the different part + diag = matrix.diagonal().pphh + e = diag.ones_like() + u2 = matrix.block_apply("pphh_ph", vec.ph) / (e * folded_matrix.omega - diag) + diff_f = matrix.block_apply("ph_pphh", u2) #A_12 * u2_fold + + diff_o = matrix.block_apply("ph_pphh", vec.pphh) #A_12 * V2 + + diff_manual = diff_o - diff_f + diff_matvec = ores.ph - fres.ph + diff = diff_manual - diff_matvec + assert np.max(np.abs(diff.to_ndarray())) < 1e-12 + def test_unfold_random_vec(self,system:str): + matrix, folded_matrix = self.construct_matrices(system) + + vec = adcc.guess_zero(matrix) + vec.set_random() + + omega = 0.5 + folded_matrix.omega = omega + fres_unfold = folded_matrix.unfold(vec) - # diag = matrix.diagonal().pphh - # e = diag.ones_like() - # u2_man = matrix.block_apply("pphh_ph", vec.ph) / (e * folded.omega - diag) - # assert u2_man.describe_symmetry() == fres_unfold.pphh.describe_symmetry() - # # Consistency of manully recomputed vector and "unfold" - # diff = u2_man - fres_unfold.pphh - # assert np.max(np.abs(diff.to_ndarray())) < 1e-12 -""" + diag = matrix.diagonal().pphh + e = diag.ones_like() + u2_man = matrix.block_apply("pphh_ph", vec.ph) / (e * folded_matrix.omega - diag) + u2_man_symm = folded_matrix.isymm.symmetrise(adcc.AmplitudeVector(pphh=u2_man)) + mres_unfold = adcc.AmplitudeVector(ph=vec.ph, pphh=u2_man_symm.pphh) + + diff_s = mres_unfold.ph - fres_unfold.ph + diff_d = mres_unfold.pphh - fres_unfold.pphh + assert np.max(np.abs(diff_s.to_ndarray())) < 1e-12 + assert np.max(np.abs(diff_d.to_ndarray())) < 1e-12 + assert mres_unfold.ph.describe_symmetry() == fres_unfold.ph.describe_symmetry() + assert mres_unfold.pphh.describe_symmetry() == fres_unfold.pphh.describe_symmetry() \ No newline at end of file diff --git a/adcc/tests/solver/davidson_test.py b/adcc/tests/solver/davidson_test.py index f6d353255..5f9d0cf7a 100644 --- a/adcc/tests/solver/davidson_test.py +++ b/adcc/tests/solver/davidson_test.py @@ -25,7 +25,7 @@ import pytest from adcc import LazyMp -from adcc.solver.davidson import jacobi_davidson, eigsh, davidson_folded_DIIS +from adcc.solver.davidson import jacobi_davidson, eigsh, davidson_folded_DIIS, eigsh_folded from adcc.misc import cached_property from ..testdata_cache import testdata_cache @@ -127,6 +127,59 @@ def test_adc2_triplets(self): assert n_states > 1 assert res.converged assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states]) + +from adcc.solver.preconditioner import JacobiPreconditioner +class TestSolverDavidsonFolded(unittest.TestCase): + @cached_property + def matrix(self): + return adcc.AdcMatrix( + "adc2", LazyMp(testdata_cache.refstate("h2o_sto3g", case="gen")) + ) + + def test_n_guesses(self): + # we have to have a guess for each state + guesses = adcc.guesses_singlet(self.matrix, n_guesses=1, block="ph") + with pytest.raises(ValueError): + eigsh_folded(self.matrix, guesses, n_ep=2) + res = eigsh_folded(self.matrix, guesses, n_ep=1, max_iter=1, + preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson") + assert len(res.eigenvalues) == 1 + # by default: construct 1 state for each guess + res = eigsh_folded(self.matrix, guesses, max_iter=1, + preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson") + assert len(res.eigenvalues) == 1 + + def test_n_block(self): + # has to be: n_ep <= n_block <= n_guesses + guesses = adcc.guesses_singlet(self.matrix, n_guesses=3, block="ph") + with pytest.raises(ValueError): + eigsh_folded(self.matrix, guesses, n_ep=2, n_block=1) + with pytest.raises(ValueError): + eigsh_folded(self.matrix, guesses, n_ep=2, n_block=4) + # defaults to n_ep + res = eigsh_folded(self.matrix, guesses, n_ep=2, max_iter=2, + preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson") + assert len(res.eigenvalues) == 2 + + def test_max_subspace(self): + # max_subspace >= 2 * n_block + guesses = adcc.guesses_singlet(self.matrix, n_guesses=3, block="ph") + with pytest.raises(ValueError): + eigsh_folded(self.matrix, guesses, n_ep=1, n_block=2, max_subspace=3) + res = eigsh_folded(self.matrix, guesses, n_ep=1, n_block=2, max_subspace=4, + preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson") + assert len(res.eigenvalues) == 1 + # max_subspace >= n_guesses + with pytest.raises(ValueError): + eigsh_folded(self.matrix, guesses, n_ep=1, n_block=1, max_subspace=2) + res = eigsh_folded(self.matrix, guesses, n_ep=1, n_block=1, max_subspace=3, + preconditioner=JacobiPreconditioner, + preconditioning_method="Davidson") + assert len(res.eigenvalues) == 1 def test_adc2_singlets_folded(self): refdata = testdata_cache.adcman_data( @@ -143,6 +196,18 @@ def test_adc2_singlets_folded(self): assert res.converged assert res.eigenvalues[:n_states] == pytest.approx(ref_singlets[:n_states]) + ref_singlets_vec_s = refdata["eigenvectors_singles"] + ref_singlets_vec_d = refdata["eigenvectors_doubles"] + + from adcc.AmplitudeVector import AmplitudeVector + print(ref_singlets_vec_s,type(ref_singlets_vec_s),type(res.eigenvectors[0]),res.eigenvectors[0]) + ref_vec = AmplitudeVector(ph=ref_singlets_vec_s,pphh=ref_singlets_vec_d) + assert res.eigenvectors[0].ph.describe_symmetry() == ref_vec.ph.describe_symmetry() + assert res.eigenvectors[0].pphh.describe_symmetry() == ref_vec.pphh.describe_symmetry() + + # assert res.eigenvectors[0].ph.describe_symmetry() == guesses[0].ph.describe_symmetry() + # assert res.eigenvectors[0].pphh.describe_symmetry() == guesses[0].pphh.describe_symmetry() + def test_adc2_triplets_folded(self): refdata = testdata_cache.adcman_data( system="h2o_sto3g", method="adc2", case="gen" @@ -158,7 +223,7 @@ def test_adc2_triplets_folded(self): assert res.converged assert res.eigenvalues[:n_states] == pytest.approx(ref_triplets[:n_states]) - def test_adc2_singlets_folded_adc1Guesses(self): + def test_adc2_singlets_folded_guess_adc1(self): from adcc.workflow import run_adc refdata = testdata_cache.adcman_data( @@ -172,6 +237,8 @@ def test_adc2_singlets_folded_adc1Guesses(self): adc1 = run_adc(matrix_adc1, method="adc1", n_singlets=8) omegas = adc1.excitation_energy_uncorrected guesses = adc1.excitation_vector + print("OMEGAS:", type(omegas),type(omegas[0])) + # Solve for singlets res = davidson_folded_DIIS(self.matrix, guesses, omegas=omegas, n_ep=8) @@ -181,7 +248,7 @@ def test_adc2_singlets_folded_adc1Guesses(self): assert res.converged assert res.eigenvalues[:n_states] == pytest.approx(ref_singlets[:n_states]) - def test_adc2_triplets_folded_adc1Guesses(self): + def test_adc2_triplets_folded_guess_adc1(self): from adcc.workflow import run_adc refdata = testdata_cache.adcman_data( diff --git a/adcc/tests/workflow_test.py b/adcc/tests/workflow_test.py index c78656457..1e1a8ab29 100644 --- a/adcc/tests/workflow_test.py +++ b/adcc/tests/workflow_test.py @@ -223,10 +223,16 @@ def test_diagonalise_adcmatrix(self): # for ADC problems with doubles-folding guesses = adcc.guesses_singlet(matrix, n_guesses=6, block="ph") res = diagonalise_adcmatrix(matrix, n_states=3, kind=kind, - guesses=guesses, fold=True) + guesses=guesses, doubles_folding=True) assert res.converged assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) + + with pytest.warns(UserWarning) as record: + diagonalise_adcmatrix(matrix, n_states=3, kind=kind, + doubles_folding=True, guesses=guesses, omegas=None) + assert len(record) == 2 + from adcc.workflow import run_adc matrix_adc1 = adcc.AdcMatrix( "adc1", adcc.LazyMp(testdata_cache.refstate("h2o_sto3g", case="gen")) @@ -235,7 +241,7 @@ def test_diagonalise_adcmatrix(self): omegas = adc1.excitation_energy_uncorrected guesses = adc1.excitation_vector res = diagonalise_adcmatrix(matrix, n_states=3, kind="singlet", - fold=True, guesses_fold="adc1", + doubles_folding=True, guesses=guesses, omegas=omegas) assert res.converged assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) @@ -255,6 +261,11 @@ def test_diagonalise_adcmatrix(self): eigensolver="davidson", guesses=guesses) + with pytest.raises(InputError): # Wrong solver for doubles folding + res = diagonalise_adcmatrix(matrix, n_states=9, kind=kind, + eigensolver="blubber", + doubles_folding=True) + def test_estimate_n_guesses(self): from adcc.workflow import estimate_n_guesses diff --git a/adcc/workflow.py b/adcc/workflow.py index 4fa30d3fd..54a23a62c 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -47,7 +47,7 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, n_guesses_doubles=None, output=sys.stdout, core_orbitals=None, frozen_core=None, frozen_virtual=None, method=None, n_singlets=None, n_triplets=None, n_spin_flip=None, - environment=None, fold=None, guesses_fold=None, + environment=None, doubles_folding=None, folding_guesses_adc1=False, omegas=None, **solverargs): """Run an ADC calculation. @@ -106,8 +106,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, guesses : list, optional Provide the guess vectors to be employed for the ADC run. Takes - preference over `n_guesses` and `n_guesses_doubles`, such that these - parameters are ignored. + preference over `n_guesses` `n_guesses_doubles` and `folding_guesses_adc1`, + such that these parameters are ignored. output : stream, optional Python stream to which output will be written. If `None` all output @@ -134,13 +134,18 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, The keywords to specify how coupling to an environment model, e.g. PE, is treated. For details see :ref:`environment`. - fold : bool, optional + doubles_folding : bool, optional Perfom modified Davidson algorithm when using doubles-folded ADC(2) matrix. - guesses_fold: str, optional + folding_guesses_adc1 : bool, optional There is one more option for the initial guess vectors of the modified - Davidson solver, i.e., "adc1" takes initial guess vectors and guess - eigenvalues as eigenvalues and eigenvectors of adc(1) matrix. + Davidson solver, i.e., "adc1", which sets initial guess vectors and omega + values as eigenvalues and eigenvectors of ADC(1) matrix. + + omegas : numpy.ndarray or NoneType, optional + Initial omega values for the modified Davidson solver when using + doubles-folded ADC(2) matrix. If not provided, the initial omega values + will be derived from the provided guess vectors. Other parameters ---------------- @@ -220,16 +225,17 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if env_matrix_term: matrix += env_matrix_term - if fold and guesses_fold == "adc1" and guesses is None: - fold_initial = run_adc(data_or_matrix, method="adc1", - kind=kind, n_states=n_states) - omegas = fold_initial.excitation_energy_uncorrected - guesses = fold_initial.excitation_vector + if doubles_folding and folding_guesses_adc1 == True and guesses is None: + fold_initial = run_adc(data_or_matrix, method="adc1", + kind=kind, n_states=n_states) + omegas = fold_initial.excitation_energy_uncorrected + guesses = fold_initial.excitation_vector + # print("OMEGAS1:", type(omegas),type(omegas[0])) diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, n_guesses_doubles=n_guesses_doubles, conv_tol=conv_tol, output=output, - eigensolver=eigensolver, fold=fold, guesses_fold=guesses_fold, + eigensolver=eigensolver, doubles_folding=doubles_folding, omegas=omegas, **solverargs) exstates = ExcitedStates(diagres) exstates.kind = kind @@ -372,8 +378,8 @@ def validate_state_parameters(reference_state, n_states=None, n_singlets=None, def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", guesses=None, n_guesses=None, n_guesses_doubles=None, - conv_tol=None, output=sys.stdout, fold=None, - guesses_fold=None, omegas=None, **solverargs): + conv_tol=None, output=sys.stdout, doubles_folding=None, + omegas=None, **solverargs): """ This function seeks appropriate guesses and afterwards proceeds to diagonalise the ADC matrix using the specified eigensolver. @@ -401,7 +407,7 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", # Set some solver-specific parameters if eigensolver == "davidson": n_guesses_per_state = 2 - if fold: + if doubles_folding: callback = setup_solver_printing( "Jacobi-Davidson using doubles-folding", matrix, kind, solver.davidson.default_print, output=output) @@ -412,6 +418,9 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", output=output) run_eigensolver = jacobi_davidson elif eigensolver == "lanczos": + if doubles_folding: + raise InputError("Lanczos eigensolver not supported for " + "doubles-folded ADC(2) matrices.") n_guesses_per_state = 1 callback = setup_solver_printing( "Lanczos", matrix, kind, solver.lanczos.default_print, @@ -420,6 +429,17 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", else: raise InputError(f"Solver {eigensolver} unknown, try 'davidson'.") + # Doubles folding warnings + if doubles_folding and guesses is not None: + warnings.warn("The guess vectors are explicitly provided.\n" + "It is recommended to use ADC(1) transition densities and excitation " + "energies as the starting guesses for the eigenvectors and eigenvalues.") + if omegas is None: + warnings.warn("The guess vectors are explicitly provided, but no initial omegas.\n" + "In the modified Davidson algorithm for the doubles-folded ADC(2) matrix, " + "the initial omega value(s) are derived from the provided guess vector(s). " + "Ensure their values and ordering match the intended states.") + # Obtain or check guesses if guesses is None: if n_guesses is None: @@ -445,7 +465,7 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson", "are explicitly provided.") solverargs.setdefault("which", "SA") - if fold and guesses_fold == "adc1": + if doubles_folding: return run_eigensolver(matrix, guesses=guesses, omegas=omegas, n_ep=n_states, conv_tol=conv_tol, callback=callback, explicit_symmetrisation=explicit_symmetrisation,