diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 516a401d..b04589a6 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): @@ -650,3 +650,56 @@ 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) + 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): + """ + Return the approximate diagonal of the ADC(2) matrix with doubles-folding. + """ + return AmplitudeVector(ph=super().diagonal().ph) + + def u2_fold(self, other): + """ + 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.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 d0588c7c..6f1fc411 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 .DIIS import DIIS 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 @@ -130,6 +197,9 @@ def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep, n_block, 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] @@ -201,19 +271,33 @@ 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([np.sqrt(r @ r) - for r in state.residuals]) + 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: + 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)] + 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 @@ -221,9 +305,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 @@ -231,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 @@ -318,9 +411,19 @@ 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. + 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 + 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 " @@ -456,6 +559,263 @@ def convergence_test(state): return state +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", + 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) + 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 + 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 + 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") + 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") + + 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(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) + 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_final(state): + state.residuals_converged = state.residual_norms < conv_tol + state.converged = np.all(state.residuals_converged) + return state.converged + + 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 residual_norm_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 omega corresponding to the new eigenvector. + if diis_omegaUpdate: + folded_matrix.omega = Av @ state.eigenvector + 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(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: + 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) + 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 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 + 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 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] + residual_norm_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=0) + 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: + # 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) + 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 + 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 + state_i.eigenvalues[n_state] = folded_matrix.omega + callback(state_i, "DIIS_stop") + + # 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))) + 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.residual_norms[n_state] = state_i.residual_norm + 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_final(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) @@ -463,3 +823,8 @@ 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) diff --git a/adcc/tests/AdcMatrix_test.py b/adcc/tests/AdcMatrix_test.py index d3ec88bd..d1cfbaf7 100644 --- a/adcc/tests/AdcMatrix_test.py +++ b/adcc/tests/AdcMatrix_test.py @@ -25,7 +25,12 @@ 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.solver.AdcMatrixFolded import AdcMatrixFolded from adcc.Intermediates import Intermediates from adcc.adc_pp.matrix import AdcBlock @@ -401,3 +406,95 @@ 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_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_matrix = self.construct_matrices(system) + + odiag = matrix.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_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() + + diag = matrix.diagonal().pphh + e = diag.ones_like() + u2_manual = matrix.block_apply("pphh_ph", vec.ph) / (e * folded_matrix.omega - diag) + + u2_fold = folded_matrix.u2_fold(vec) + + diff = u2_manual - u2_fold.pphh + assert np.max(np.abs(diff.to_ndarray())) < 1e-12 + + 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() + + ores = matrix @ vec + fres = folded_matrix @ vec + + assert ores.ph.describe_symmetry() == fres.ph.describe_symmetry() + assert not hasattr(fres, "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_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 2182d49d..5f9d0cf7 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, eigsh_folded from adcc.misc import cached_property from ..testdata_cache import testdata_cache @@ -127,3 +127,147 @@ 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( + 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]) + + 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" + )["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_guess_adc1(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 + print("OMEGAS:", type(omegas),type(omegas[0])) + + # 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_guess_adc1(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]) diff --git a/adcc/tests/workflow_test.py b/adcc/tests/workflow_test.py index 6dbc8e6d..1e1a8ab2 100644 --- a/adcc/tests/workflow_test.py +++ b/adcc/tests/workflow_test.py @@ -220,6 +220,32 @@ 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=kind, + 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")) + ) + 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", + doubles_folding=True, + guesses=guesses, omegas=omegas) + assert res.converged + assert res.eigenvalues[:n_states] == approx(ref_singlets[:n_states]) + with pytest.raises(InputError): # Too low tolerance # SCF tolerance = 1e-14 currently res = diagonalise_adcmatrix(matrix, n_states=9, kind=kind, @@ -235,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 6c79fb4c..54a23a62 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, doubles_folding=None, folding_guesses_adc1=False, + omegas=None, **solverargs): """Run an ADC calculation. Main entry point to run an ADC calculation. The reference to build the ADC @@ -105,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 @@ -133,12 +134,34 @@ 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`. + doubles_folding : bool, optional + Perfom modified Davidson algorithm when using doubles-folded ADC(2) matrix. + + folding_guesses_adc1 : bool, optional + There is one more option for the initial guess vectors of the modified + 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 ---------------- 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 +225,18 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if env_matrix_term: matrix += env_matrix_term + 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, **solverargs) + eigensolver=eigensolver, doubles_folding=doubles_folding, + omegas=omegas, **solverargs) exstates = ExcitedStates(diagres) exstates.kind = kind exstates.spin_change = spin_change @@ -347,7 +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, **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. @@ -375,11 +407,20 @@ 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 doubles_folding: + 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": + 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, @@ -388,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: @@ -413,10 +465,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 doubles_folding: + 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,