From 4d12f12db97e7a9a894307dfb5b9f8337f8819a3 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Tue, 17 May 2022 15:28:10 +0200 Subject: [PATCH] Projectors and projected ADC matrix (#149) --- adcc/AdcMatrix.py | 90 ++++- adcc/MoSpaces.py | 12 +- adcc/guess/__init__.py | 25 +- adcc/guess/guesses_from_diagonal.py | 30 +- adcc/projection.py | 355 ++++++++++++++++++ adcc/test_AdcMatrix.py | 64 +++- adcc/test_projection.py | 243 ++++++++++++ adcc/testdata/cache.py | 2 +- .../water_core_by_projection.py | 46 +++ 9 files changed, 847 insertions(+), 20 deletions(-) create mode 100644 adcc/projection.py create mode 100644 adcc/test_projection.py create mode 100755 examples/2_core_excitations/water_core_by_projection.py diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 2028252e..016475b4 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -28,6 +28,7 @@ from .adc_pp import matrix as ppmatrix from .timings import Timer, timed_member_call from .AdcMethod import AdcMethod +from .functions import ones_like from .Intermediates import Intermediates from .AmplitudeVector import AmplitudeVector @@ -355,7 +356,7 @@ def block_view(self, block): # data structure as the AdcMatrix class as these block # are inherently different than an AdcMatrix (no Hermiticity # for example) and basically they only need to support some - # form of matrix-vector product and some stastics like + # form of matrix-vector product and some statistics like # spaces and sizes etc. block_orders = {bl: None for bl in self.block_orders.keys()} block_orders[block] = self.block_orders[block] @@ -628,3 +629,90 @@ def block_view(self, block): "shifted ADC matrices.") # TODO The way to implement this is to ask the inner matrix to # a block_view and then wrap that in an AdcMatrixShifted. + + +class AdcMatrixProjected(AdcMatrix): + def __init__(self, matrix, excitation_blocks, core_orbitals=None, + outer_virtuals=None): + """ + Initialise a projected ADC matrix, i.e. represents the expression + ``P @ M @ P`` where ``P`` is a projector onto a subset of + ``excitation_blocks``. + + The ``excitation_blocks`` are defined by partitioning the ``o1`` occupied + and ``v1`` virtual space of the ``matrix.mospaces`` into a core-occupied + ``c``, valence-occupied ``o``, inner-virtual ``v`` and outer-virtual ``w``. + This matrix will only keep selected blocks in the amplitudes non-zero, which + are selected in the ``excitation_blocks`` list + (e.g. ``["cv", "ccvv", "ocvv"]``). + + For details on the option how to select the spaces, see the documentation + in :py:`adcc.ReferenceState.__init__` (``outer_virtuals`` follows the same + rules as ``frozen_virtuals``). + + Parameters + ---------- + matrix : AdcMatrix + Matrix which is projected + excitation_blocks : list + Excitation blocks to keep in the Amplitudes. + core_orbitals : int or list or tuple, optional + The orbitals to be put into the ``c`` core space. + outer_virtuals : int or list or tuple, optional + The virtuals to be put into the ``w`` outer-virtual space. + """ + from .projection import Projector, SubspacePartitioning + + for sp in excitation_blocks: + if not any(len(sp) == len(ax_sp) + for ax_sp in matrix.axis_spaces.values()): + raise ValueError(f"Invalid partition block {sp}.") + + super().__init__(matrix.method, matrix.ground_state, + block_orders=matrix.block_orders, + intermediates=matrix.intermediates) + partitioning = SubspacePartitioning(matrix.mospaces, core_orbitals, + outer_virtuals) + + projectors = {} + for block in matrix.axis_spaces.keys(): + block_partitions = [sp for sp in excitation_blocks + if len(sp) == len(matrix.axis_spaces[block])] + projectors[block] = Projector(matrix.axis_spaces[block], + partitioning, block_partitions) + self.projectors = projectors + + def apply_projection(self, in_ampl): + return AmplitudeVector(**{ + block: self.projectors[block] @ in_ampl[block] + for block in in_ampl.keys() + }) + + def matvec(self, in_ampl): + in_proj = self.apply_projection(in_ampl) + out = super().matvec(in_proj) + return self.apply_projection(out) + + def block_apply(self, block, in_vec): + inblock, outblock = block.split("_") + in_proj = self.projectors[inblock].apply(in_vec) + ret = super().block_apply(block, in_proj) + return self.projectors[outblock].apply(ret) + + def diagonal(self): + blocks = {} + for (block, diagblock) in super().diagonal().items(): + # On the diagonal don't set the ignored amplitudes to zero, + # but instead set them to a very large value to (a) avoid these being + # selected for the initial guess and (b) ensure the preconditioning + # naturally reduces components along this direction. + P = self.projectors[block] + one = ones_like(diagblock) + blocks[block] = P @ diagblock + 100000 * (one - P @ one) + return AmplitudeVector(**blocks) + + def block_view(self, block): + raise NotImplementedError("Block-view not yet implemented for " + "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. diff --git a/adcc/MoSpaces.py b/adcc/MoSpaces.py index 8c8ce668..6782e050 100644 --- a/adcc/MoSpaces.py +++ b/adcc/MoSpaces.py @@ -48,7 +48,14 @@ def split_spaces(spacestr): raise ValueError(f"Invalid space string '{spacestr}'.") -def expand_spaceargs(hfdata, **spaceargs): +def expand_spaceargs(hfdata_or_n_orbs, **spaceargs): + if isinstance(hfdata_or_n_orbs, tuple): + (noa, nob) = hfdata_or_n_orbs + else: + hfdata = hfdata_or_n_orbs + noa = hfdata.n_orbs_alpha + nob = hfdata.n_orbs_beta + if isinstance(spaceargs.get("frozen_core", None), bool) \ and spaceargs.get("frozen_core", None): # Determine number of frozen core electrons automatically @@ -93,7 +100,6 @@ def expand_to_list(space, entry, from_min=None, from_max=None): "core_orbitals, frozen_virtual is an " "iterable, all must be.") - noa = hfdata.n_orbs_alpha n_orbs = [0, 0] for key in ["frozen_core", "core_orbitals"]: if key not in spaceargs or not spaceargs[key]: @@ -110,7 +116,7 @@ def expand_to_list(space, entry, from_min=None, from_max=None): if key in spaceargs and spaceargs[key]: spaceargs[key] = np.concatenate(( expand_to_list(key, spaceargs[key][0], from_max=noa), - expand_to_list(key, spaceargs[key][1], from_max=noa) + noa + expand_to_list(key, spaceargs[key][1], from_max=nob) + noa )).tolist() return spaceargs diff --git a/adcc/guess/__init__.py b/adcc/guess/__init__.py index 0cea7ebb..b2690876 100644 --- a/adcc/guess/__init__.py +++ b/adcc/guess/__init__.py @@ -28,6 +28,23 @@ "guesses_spin_flip", "guess_symmetries"] +def guess_kwargs_kind(kind): + """ + Return the kwargs required to be passed to `guesses_from_diagonal` to + computed states of the passed excitation `kind`. + """ + kwargsmap = dict( + singlet=dict(spin_block_symmetrisation="symmetric", spin_change=0), + triplet=dict(spin_block_symmetrisation="antisymmetric", spin_change=0), + spin_flip=dict(spin_block_symmetrisation="none", spin_change=-1), + any=dict(spin_block_symmetrisation="none", spin_change=0), + ) + try: + return kwargsmap[kind] + except KeyError: + raise ValueError(f"Kind not known: {kind}") + + def guesses_singlet(matrix, n_guesses, block="ph", **kwargs): """ Obtain guesses for computing singlet states by inspecting the passed @@ -41,8 +58,7 @@ def guesses_singlet(matrix, n_guesses, block="ph", **kwargs): kwargs Any other argument understood by guesses_from_diagonal. """ return guesses_from_diagonal(matrix, n_guesses, block=block, - spin_block_symmetrisation="symmetric", - spin_change=0, **kwargs) + **guess_kwargs_kind("singlet"), **kwargs) def guesses_triplet(matrix, n_guesses, block="ph", **kwargs): @@ -58,8 +74,7 @@ def guesses_triplet(matrix, n_guesses, block="ph", **kwargs): kwargs Any other argument understood by guesses_from_diagonal. """ return guesses_from_diagonal(matrix, n_guesses, block=block, - spin_block_symmetrisation="antisymmetric", - spin_change=0, **kwargs) + **guess_kwargs_kind("triplet"), **kwargs) # guesses for computing any state (singlet or triplet) @@ -79,4 +94,4 @@ def guesses_spin_flip(matrix, n_guesses, block="ph", **kwargs): kwargs Any other argument understood by guesses_from_diagonal. """ return guesses_from_diagonal(matrix, n_guesses, block=block, - spin_change=-1, **kwargs) + **guess_kwargs_kind("spin_flip"), **kwargs) diff --git a/adcc/guess/guesses_from_diagonal.py b/adcc/guess/guesses_from_diagonal.py index 8a53a137..2d94f972 100644 --- a/adcc/guess/guesses_from_diagonal.py +++ b/adcc/guess/guesses_from_diagonal.py @@ -33,7 +33,7 @@ def guesses_from_diagonal(matrix, n_guesses, block="ph", spin_change=0, spin_block_symmetrisation="none", - degeneracy_tolerance=1e-14): + degeneracy_tolerance=1e-14, max_diagonal_value=1000): """ Obtain guesses by inspecting a block of the diagonal of the passed ADC matrix. The symmetry of the returned vectors is already set-up properly. @@ -59,6 +59,9 @@ def guesses_from_diagonal(matrix, n_guesses, block="ph", spin_change=0, degeneracy_tolerance Tolerance for two entries of the diagonal to be considered degenerate, i.e. identical. + max_diagonal_value + Maximal diagonal value, which is considered as a valid candidate + to form a guess. """ if not isinstance(matrix, AdcMatrixlike): raise TypeError("matrix needs to be of type AdcMatrixlike") @@ -88,7 +91,8 @@ def guesses_from_diagonal(matrix, n_guesses, block="ph", spin_change=0, raise ValueError(f"Don't know how to generate guesses for block {block}") return guessfunction(matrix, n_guesses, spin_change, - spin_block_symmetrisation, degeneracy_tolerance) + spin_block_symmetrisation, degeneracy_tolerance, + max_diagonal_value) class TensorElement: @@ -195,7 +199,8 @@ def telem_nospin(telem): def guesses_from_diagonal_singles(matrix, n_guesses, spin_change=0, spin_block_symmetrisation="none", - degeneracy_tolerance=1e-14): + degeneracy_tolerance=1e-14, + max_diagonal_value=1000): motrans = MoIndexTranslation(matrix.mospaces, matrix.axis_spaces["ph"]) if n_guesses == 0: return [] @@ -207,10 +212,13 @@ def guesses_from_diagonal_singles(matrix, n_guesses, spin_change=0, # Search of the smallest elements # This predicate checks an index is an allowed element for the singles - # part of the guess vectors and has the requested spin-change + # part of the guess vectors and has the requested spin-change. + # Also it filters out too large diagonal entries (which are essentially + # hopeless to give useful excitations) def pred_singles(telem): return (ret[0].ph.is_allowed(telem.index) - and telem.spin_change == spin_change) + and telem.spin_change == spin_change + and abs(telem.value) <= max_diagonal_value) elements = find_smallest_matching_elements( pred_singles, matrix.diagonal().ph, motrans, n_guesses, @@ -265,7 +273,8 @@ def telem_nospin(telem): def guesses_from_diagonal_doubles(matrix, n_guesses, spin_change=0, spin_block_symmetrisation="none", - degeneracy_tolerance=1e-14): + degeneracy_tolerance=1e-14, + max_diagonal_value=1000): if n_guesses == 0: return [] @@ -286,9 +295,14 @@ def guesses_from_diagonal_doubles(matrix, n_guesses, spin_change=0, guesses_d, matrix.mospaces, df02, df13, spin_change_twice, degeneracy_tolerance ) - # Resize in case less guesses found than requested - return ret[:n_found] + ret = ret[:n_found] + + # Filter out elements above the noted diagonal value + diagonal_elements = [ret_d.pphh.dot(matrix.diagonal().pphh * ret_d.pphh) + for ret_d in ret] + return [ret[i] for (i, elem) in enumerate(diagonal_elements) + if elem <= max_diagonal_value] # TODO Generic algorithm for forming orthogonal spin components of arbitrary # size. Could be useful for the doubles guesses later on. diff --git a/adcc/projection.py b/adcc/projection.py new file mode 100644 index 00000000..d6eee258 --- /dev/null +++ b/adcc/projection.py @@ -0,0 +1,355 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +## --------------------------------------------------------------------- +## +## Copyright (C) 2022 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 libadcc +import warnings +import itertools + +import numpy as np + +from .guess import guess_kwargs_kind, guess_symmetries +from .Tensor import Tensor +from .MoSpaces import expand_spaceargs +from .Symmetry import Symmetry +from .AdcMatrix import AdcMatrixlike +from .AmplitudeVector import AmplitudeVector + +__all__ = ["SubspacePartitioning", "Projector"] + + +class SubspacePartitioning: + def __init__(self, mospaces, core_orbitals=None, outer_virtuals=None, + partitions=None, aliases=None): + """ + Partition the "o1" and "v1" occupied and virtual orbital subspace into two + (or more) parts. In combination with the :py:`adcc.Projector` + class this allows to disable some excitations in an amplitude vector. + + Parameters + ---------- + mospaces : MoSpaces + MoSpaces to base partitioning on. + core_orbitals : int or list or tuple, optional + The orbitals to be put into the ``c`` core space. + For details on the option how to select these spaces, see the + documentation in :py:`adcc.ReferenceState.__init__` (``outer_virtuals`` + follows the same rules as ``frozen_virtuals``). + outer_virtuals : int or list or tuple, optional + The virtuals to be put into the ``w`` outer-virtual space. + partitions : dict, optional + Partitioning of the spaces: Mapping from the partitioned spaces to + the indices of the orbitals which belong to this partition. The + indices are given relative to the corresponding full orbital subspace + (advanced argument). + aliases : dict, optional + Aliase mapping to give partitions a more user-friendly identifier + (advanced argument). + """ + if mospaces.has_core_occupied_space or "v2" in mospaces.subspaces_virtual: + raise ValueError("Cannot use MoSpaces with core-valence separation " + "or frozen virtuals.") + assert mospaces.subspaces_occupied == ["o1"] + + if not partitions or not aliases: + if core_orbitals is None and outer_virtuals is None: + raise ValueError("One of core_orbitals or outer_virtuals should " + "be passed.") + + # core / outer virtual part of the "normal" occupied / virtual subspace + noa = mospaces.n_orbs_alpha("o1") + nob = mospaces.n_orbs_beta("o1") + nva = mospaces.n_orbs_alpha("v1") + nvb = mospaces.n_orbs_beta("v1") + core = expand_spaceargs((noa, nob), core_orbitals=core_orbitals) + core = core["core_orbitals"] + fv = expand_spaceargs((nva, nvb), frozen_virtual=outer_virtuals) + fv = fv["frozen_virtual"] + + # The active core and virtual parts: + occ = [i for i in range(mospaces.n_orbs("o1")) if i not in core] + av = [i for i in range(mospaces.n_orbs("v1")) if i not in fv] + assert occ and av + + # Setup the space partitionings and aliases + partitions = {"o1.1": occ, "v1.1": av, } + aliases = {"o": "o1.1", "v": "v1.1", } + if core: + partitions["o1.2"] = core + aliases["c"] = "o1.2" + if fv: + partitions["v1.2"] = fv + aliases["w"] = "v1.2" + self.__init__(mospaces, partitions=partitions, aliases=aliases) + else: + assert core_orbitals is None + assert outer_virtuals is None + + self.partitions = partitions + self.aliases = aliases + self.mospaces = mospaces + + # Determine whether the partitioning keeps spin symmetry. + self.keeps_spin_symmetry = mospaces.restricted + for (label, indices) in partitions.items(): + if not self.keeps_spin_symmetry: + break + nalpha = mospaces.n_orbs_alpha(label[:2]) + alpha_part = [i for i in indices if i < nalpha] + beta_part = [i - nalpha for i in indices if i >= nalpha] + self.keeps_spin_symmetry = alpha_part == beta_part + + if mospaces.restricted and not self.keeps_spin_symmetry: + warnings.warn( + "For restricted references the case of a space " + "partitioning, which breaks spin (i.e. where excitations " + "differing only in spin are placed in different " + "partitions) has not been tested. You are on your own." + ) + + def list_space_partitions(self, space): + """ + Return the partitions into which the passed space is split. + Resolves aliases, i.e. does not return not the list of strings + such as ['o1.1', 'o1.2], but the aliased equivalent ["o", "c"]. + """ + partitions = [s for s in self.partitions if s.startswith(space + ".")] + partitions = [key for (key, value) in self.aliases.items() + if value in partitions] + return sorted(partitions) + + def get_partition(self, label): + """Return the index list corresponding to a partition label""" + return self.partitions[self.aliases[label]] + + +class Projector: + def __init__(self, subspaces, partitioning, blocks_to_keep): + """ + Initialise a projector, which upon multiplication with a tensor sets + a number of partitioning blocks of this tensor to zero. These are defined + by the passed `partitioning` and the list of `blocks_to_keep`. + + Parameters + ---------- + subspaces : list + List of subspaces to act upon (e.g. ``["o1", "v1", "v1"]``) + partitioning : SubspacePartitioning + Partitioning of the ``o1`` occupied and ``v1`` virtual space. + blocks_to_keep : list + Blocks which are kept as non-zero (e.g. ``["cv", "ccvv", "ocvv"]``). + """ + if len(subspaces) == 3 or len(subspaces) > 4 \ + or (len(subspaces) > 2 and subspaces[1] == subspaces[2]): + raise NotImplementedError( + f"Projector not implemented for subspaces {subspaces}" + ) + + def normalise_block(labels): + labels = list(labels) + if len(subspaces) <= 1: + return labels + if subspaces[0] == subspaces[1] and labels[0] > labels[1]: + labels[0], labels[1] = labels[1], labels[0] + if len(subspaces) <= 2: + return labels + if len(subspaces) == 3 or subspaces[1] == subspaces[2] \ + or len(subspaces) > 4: + raise NotImplementedError("not implemented") + if subspaces[2] == subspaces[3] and labels[2] > labels[3]: + labels[2], labels[3] = labels[3], labels[2] + return labels + + blocks_to_keep = ["".join(normalise_block(block)) + for block in blocks_to_keep] + partitions = [partitioning.list_space_partitions(space) + for space in subspaces] + all_partition_blocks = ["".join(normalise_block(block)) + for block in itertools.product(*partitions)] + for block in blocks_to_keep: + if block not in all_partition_blocks: + raise ValueError(f"Partition block {block} not known.") + + sym = Symmetry(partitioning.mospaces, "".join(subspaces)) + sym.irreps_allowed = ["A"] + if partitioning.keeps_spin_symmetry: + if len(subspaces) == 1: + sym.spin_block_maps = [("a", "b", 1)] + elif len(subspaces) == 2: + sym.spin_block_maps = [("aa", "bb", 1), ("ab", "ba", 1)] + elif len(subspaces) == 4: + sym.spin_block_maps = [ + ("aaaa", "bbbb", 1), ("aaab", "bbba", 1), + ("aaba", "bbab", 1), ("abaa", "babb", 1), + ("baaa", "abbb", 1), ("aabb", "bbaa", 1), + ("abba", "baab", 1), ("abab", "baba", 1) + ] + else: + raise NotImplementedError("not implemented") + if len(subspaces) == 4: + assert subspaces[0] == subspaces[1] + assert subspaces[2] == subspaces[3] + sym.permutations = ["ijkl", "jikl", "ijlk"] + + # Allocate zero tensor to hold projection kernel: + kernel = Tensor(sym) + + # Generic implementation not using symmetry + for block in blocks_to_keep: + ranges = [partitioning.get_partition(b) for b in block] + for index in itertools.product(*ranges): + kernel[index] = 1.0 + + self.subspaces = subspaces + self.partitioning = partitioning + self.blocks_to_keep = blocks_to_keep + self.kernel = kernel + + def matvec(self, v): + return self.kernel * v + + def __matmul__(self, other): + if isinstance(other, libadcc.Tensor): + return self.matvec(other) + if isinstance(other, list): + if all(isinstance(elem, libadcc.Tensor) for elem in other): + return [self.matvec(ov) for ov in other] + return NotImplemented + + +def transfer_amplitude_cvs_to_full(mospaces_cvs, cvs, symmetry_full, tol=1e-12): + """ + Take a CVS vector `cvs` taken from the MoSpace `mospaces_cvs` and copy + its data into a "full" ADC amplitude constructed according to the symmetry + `symmetry_full`. Notice that this will fail if the symmetry of `cvs` is + not compatible (e.g. differs in targeted spin manifold, restricted versus + unrestricted etc.). + """ + if cvs.keys() != symmetry_full.keys(): + raise ValueError("Blocks present in CVS and full vector should agree.") + if not mospaces_cvs.has_core_occupied_space: + raise ValueError("Should have CVS enabled in mospaces_cvs") + if not isinstance(cvs, AmplitudeVector): + raise TypeError("Expected cvs to be an AmplitudeVector.") + mospaces_full = symmetry_full["ph"].mospaces + n_c_a = mospaces_cvs.n_orbs_alpha("o2") + n_o_a = mospaces_cvs.n_orbs_alpha("o1") + n_of_a = mospaces_full.n_orbs_alpha("o1") + + cvs_contiguous_in_full = ( + mospaces_cvs.core_orbitals[:n_c_a] + == mospaces_full.occupied_orbitals[:n_c_a] + and mospaces_cvs.core_orbitals[n_c_a:] + == mospaces_full.occupied_orbitals[n_of_a:n_of_a + n_c_a] + ) + if not cvs_contiguous_in_full: + raise NotImplementedError( + "For transfer CVS to full only the case where the CVS space is " + "contiguous in the full space has been implemented, i.e. when the " + "`core_orbitals` are chosen as an integer in the run_adc routines." + ) + full = AmplitudeVector(**{block: Tensor(sym) + for block, sym in symmetry_full.items()}) + + # Singles block + assert cvs.ph.subspaces == ["o2", "v1"] + assert full.ph.subspaces == ["o1", "v1"] + assert cvs.ph.shape == (mospaces_cvs.n_orbs("o2"), mospaces_cvs.n_orbs("v1")) + + cvs_ph = cvs.ph.to_ndarray() + full_ph = np.zeros(full.ph.shape) + full_ph[0:n_c_a, :] = cvs_ph[:n_c_a, :] + full_ph[n_of_a:n_of_a + n_c_a, :] = cvs_ph[n_c_a:, :] + full.ph.set_from_ndarray(full_ph, tol) + + # Doubles block + if "pphh" in cvs: + assert cvs.pphh.subspaces == ["o1", "o2", "v1", "v1"] + assert full.pphh.subspaces == ["o1", "o1", "v1", "v1"] + + # Note 1/sqrt(2) factor is there because the ocvv ocurrs twice + # in the full amplitude, but only once in the CVS amplitude. + cvs_pphh = cvs.pphh.to_ndarray() / np.sqrt(2) + cvs_pphh_aa = cvs_pphh[:n_o_a, :n_c_a, :, :] + cvs_pphh_ab = cvs_pphh[:n_o_a, n_c_a:, :, :] + cvs_pphh_ba = cvs_pphh[n_o_a:, :n_c_a, :, :] + cvs_pphh_bb = cvs_pphh[n_o_a:, n_c_a:, :, :] + + full_pphh = np.zeros(full.pphh.shape) + full_pphh[n_c_a:n_of_a, 0:n_c_a, :, :] = cvs_pphh_aa + full_pphh[0:n_c_a, n_c_a:n_of_a, :, :] = ( + -cvs_pphh_aa.transpose((1, 0, 2, 3)) + ) + full_pphh[n_c_a:n_of_a, n_of_a:n_of_a + n_c_a, :, :] = cvs_pphh_ab + full_pphh[n_of_a:n_of_a + n_c_a, n_c_a:n_of_a, :, :] = ( + -cvs_pphh_ab.transpose((1, 0, 2, 3)) + ) + full_pphh[n_of_a + n_c_a:, 0:n_c_a, :, :] = cvs_pphh_ba + full_pphh[0:n_c_a, n_of_a + n_c_a:, :, :] = ( + -cvs_pphh_ba.transpose((1, 0, 2, 3)) + ) + full_pphh[n_of_a + n_c_a:, n_of_a:n_of_a + n_c_a, :, :] = cvs_pphh_bb + full_pphh[n_of_a:n_of_a + n_c_a, n_of_a + n_c_a:, :, :] = ( + -cvs_pphh_bb.transpose((1, 0, 2, 3)) + ) + full.pphh.set_from_ndarray(full_pphh, tol) + return full + + +def transfer_cvs_to_full(state_matrix_cvs, matrix_full, vector=None, kind=None, + spin_change=0, spin_block_symmetrisation="none"): + """ + Transfer `vector` (or a list of vectors) from the CVS space to the full space + defined by the ADC matrix `matrix_full`. `state_matrix_cvs` is either the + CVS matrix or the state containing the solved CVS excitations. To properly + set up the symmetry of the returned vectors, set the `kind`, `spin_change` + and `spin_block_symmetrisation` parameters. If the first argument is an + `ExcitedStates` object, the symmetry is automatically detected. + """ + if kind is None and hasattr(state_matrix_cvs, "kind"): + kind = state_matrix_cvs.kind + + if spin_change is None and spin_block_symmetrisation is None: + if kind is None: + raise ValueError("kind needs to be given if first argument is not an " + "ExcitedStates object and spin symmetry setup is not " + "explicitly given.") + return transfer_cvs_to_full(state_matrix_cvs, matrix_full, vector, kind, + **guess_kwargs_kind(kind)) + + if vector is None: + if hasattr(state_matrix_cvs, "excitation_vector"): + vector = state_matrix_cvs.excitation_vector + else: + raise ValueError("vector needs to be given if first argument is not an " + "ExcitedStates object.") + if isinstance(vector, list): + return [transfer_cvs_to_full(state_matrix_cvs, matrix_full, v, kind, + **guess_kwargs_kind(kind)) for v in vector] + + if isinstance(state_matrix_cvs, AdcMatrixlike): + mospaces_cvs = state_matrix_cvs.mospaces + else: + mospaces_cvs = state_matrix_cvs.matrix.mospaces + + sym = guess_symmetries(matrix_full, spin_change=spin_change, + spin_block_symmetrisation=spin_block_symmetrisation) + return transfer_amplitude_cvs_to_full(mospaces_cvs, vector, sym) diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index ea5d0b92..5d0a1e69 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -28,7 +28,7 @@ from numpy.testing import assert_allclose -from adcc.AdcMatrix import AdcMatrixShifted, AdcExtraTerm +from adcc.AdcMatrix import AdcExtraTerm, AdcMatrixProjected, AdcMatrixShifted from adcc.adc_pp.matrix import AdcBlock from adcc.testdata.cache import cache @@ -319,4 +319,64 @@ def template_matmul(self, case): assert np.max(np.abs(diff_s.to_ndarray())) < 1e-12 assert np.max(np.abs(diff_d.to_ndarray())) < 1e-12 - # TODO Test to_dense_matrix, compute_apply + # TODO Test block_view, block_apply + + +@expand_test_templates(testcases) +class TestAdcMatrixProjected(unittest.TestCase): + def construct_matrices(self, case, n_core, n_virt): + from .test_projection import construct_nonzero_blocks + + reference_state = cache.refstate[case] + ground_state = adcc.LazyMp(reference_state) + matrix = adcc.AdcMatrix("adc3", ground_state) + + out = construct_nonzero_blocks(reference_state.mospaces, n_core, n_virt) + spaces, nonzero_blocks = out + + excitation_blocks = spaces["ph"] + spaces["pphh"] + projected = AdcMatrixProjected(matrix, excitation_blocks, + core_orbitals=n_core, + outer_virtuals=n_virt) + return matrix, projected, nonzero_blocks + + def template_diagonal(self, case): + from .test_projection import assert_nonzero_blocks + + out = self.construct_matrices(case, n_core=2, n_virt=1) + matrix, projected, nonzeros = out + + for block in ("ph", "pphh"): + odiag = matrix.diagonal()[block] + pdiag = projected.diagonal()[block] + assert_nonzero_blocks(odiag, pdiag, nonzeros[block], zero_value=100000) + # TODO Manually verified to be identical, however, string parsing + # of the describe_symmetry output is not super reliable and so this + # test does not pass in CI. + # assert_equal_symmetry(odiag, pdiag) + + def template_matmul(self, case): + from .test_projection import (assert_equal_symmetry, + assert_nonzero_blocks) + + out = self.construct_matrices(case, n_core=1, n_virt=1) + matrix, projected, nonzeros = out + + spin_block_symmetrisation = "none" + if "h2o" in case: + spin_block_symmetrisation = "symmetric" + vec = adcc.guess_zero(matrix, + spin_block_symmetrisation=spin_block_symmetrisation) + vec.set_random() + pvec = projected.apply_projection(vec.copy()) # only apply projection + + pres = projected @ vec + ores = matrix @ pvec + res_for_sym = matrix @ vec + + assert_equal_symmetry(res_for_sym.ph, pres.ph) + assert_equal_symmetry(res_for_sym.pphh, pres.pphh) + assert_nonzero_blocks(ores.ph, pres.ph, nonzeros["ph"], tol=1e-14) + assert_nonzero_blocks(ores.pphh, pres.pphh, nonzeros["pphh"], tol=1e-14) + + # TODO Test block_view, block_apply diff --git a/adcc/test_projection.py b/adcc/test_projection.py new file mode 100644 index 00000000..38aa5e42 --- /dev/null +++ b/adcc/test_projection.py @@ -0,0 +1,243 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +## --------------------------------------------------------------------- +## +## Copyright (C) 2022 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 re +import adcc +import unittest +import itertools + +from adcc.projection import (Projector, SubspacePartitioning, + transfer_cvs_to_full) +from adcc.testdata.cache import cache + +import numpy as np + +from numpy.testing import assert_allclose + +from .misc import expand_test_templates +from .HfCounterData import HfCounterData + + +class TestSubspacePartitioning(unittest.TestCase): + def test_restricted(self): + n_alpha = 5 + n_beta = 5 + n_bas = 20 + n_orbs_alpha = 10 + restricted = True + data = HfCounterData(n_alpha, n_beta, n_bas, n_orbs_alpha, restricted) + refstate = adcc.ReferenceState(data, import_all_below_n_orbs=None) + + partitioning = SubspacePartitioning(refstate.mospaces, core_orbitals=2, + outer_virtuals=3) + assert "v" in partitioning.aliases + assert "w" in partitioning.aliases + assert "o" in partitioning.aliases + assert "c" in partitioning.aliases + assert partitioning.keeps_spin_symmetry + + assert partitioning.list_space_partitions("o1") == ["c", "o"] + assert partitioning.list_space_partitions("v1") == ["v", "w"] + + assert partitioning.get_partition("c") == [0, 1, 5, 6] + assert partitioning.get_partition("o") == [2, 3, 4, 7, 8, 9] + assert partitioning.get_partition("v") == [0, 1, 5, 6] + assert partitioning.get_partition("w") == [2, 3, 4, 7, 8, 9] + + def test_unrestricted(self): + n_alpha = 6 + n_beta = 5 + n_bas = 20 + n_orbs_alpha = 10 + restricted = False + data = HfCounterData(n_alpha, n_beta, n_bas, n_orbs_alpha, restricted) + refstate = adcc.ReferenceState(data, import_all_below_n_orbs=None) + + partitioning = SubspacePartitioning(refstate.mospaces, core_orbitals=2, + outer_virtuals=2) + assert "v" in partitioning.aliases + assert "w" in partitioning.aliases + assert "o" in partitioning.aliases + assert "c" in partitioning.aliases + assert not partitioning.keeps_spin_symmetry + + assert partitioning.list_space_partitions("o1") == ["c", "o"] + assert partitioning.list_space_partitions("v1") == ["v", "w"] + + assert partitioning.get_partition("c") == [0, 1, 6, 7] + assert partitioning.get_partition("o") == [2, 3, 4, 5, 8, 9, 10] + assert partitioning.get_partition("v") == [0, 1, 4, 5, 6] + assert partitioning.get_partition("w") == [2, 3, 7, 8] + + +def construct_nonzero_blocks(mospaces, n_core, n_virt): + nva = mospaces.n_orbs_alpha("v1") + nvb = mospaces.n_orbs_beta("v1") + noa = mospaces.n_orbs_alpha("o1") + nob = mospaces.n_orbs_beta("o1") + c = dict(a=slice(0, n_core), + b=slice(noa, noa + n_core)) + o = dict(a=slice(n_core, noa), + b=slice(n_core + noa, noa + nob)) + v = dict(a=slice(0, nva - n_virt), + b=slice(nva, nvb + nva - n_virt)) + w = dict(a=slice(nva - n_virt, nva), + b=slice(nva - n_virt + nvb, nvb + nva)) + + spaces_ph = ["cv", "ow"] + nonzero_blocks_ph = [] + for (s1, s2) in itertools.product(["a", "b"], ["a", "b"]): + nonzero_blocks_ph.append((c[s1], v[s2])) + nonzero_blocks_ph.append((o[s1], w[s2])) + + spaces_pphh = ["covv", "ccvw", "covw"] + nonzero_blocks_pphh = [] + for (s1, s2, s3, s4) in itertools.product(["a", "b"], ["a", "b"], + ["a", "b"], ["a", "b"]): + nonzero_blocks_pphh.append((c[s1], o[s2], v[s3], v[s4])) + nonzero_blocks_pphh.append((o[s1], c[s2], v[s3], v[s4])) + # + nonzero_blocks_pphh.append((c[s1], c[s2], v[s3], w[s4])) + nonzero_blocks_pphh.append((c[s1], c[s2], w[s3], v[s4])) + # + nonzero_blocks_pphh.append((c[s1], o[s2], v[s3], w[s4])) + nonzero_blocks_pphh.append((o[s1], c[s2], v[s3], w[s4])) + nonzero_blocks_pphh.append((c[s1], o[s2], w[s3], v[s4])) + nonzero_blocks_pphh.append((o[s1], c[s2], w[s3], v[s4])) + + spaces = dict(ph=spaces_ph, pphh=spaces_pphh) + nonzeros = dict(ph=nonzero_blocks_ph, pphh=nonzero_blocks_pphh) + return spaces, nonzeros + + +def assert_nonzero_blocks(orig, projected, nonzero_blocks, zero_value=0, tol=0): + diff = (orig - projected).to_ndarray() + reset = projected.to_ndarray() + for block in nonzero_blocks: + assert np.max(np.abs(diff[block])) <= tol # Values unchanged + reset[block] = zero_value + cond = np.max(np.abs(np.abs(reset) - zero_value)) <= tol + if not cond: + print(reset) + # All zero blocks are equal to zero_value or -zero_value + assert np.max(np.abs(np.abs(reset) - zero_value)) <= tol + + +def assert_equal_symmetry(sym1, sym2): + def split_blocks(sym): + str_sym = sym.describe_symmetry() + + # Perform some cleanup (effects only point group symmetry) + str_sym = str_sym.replace("0(1)", "0(0)").replace("1(1)", "1(0)") + + blocks = [] + for line in str_sym.split("\n"): + if re.match(r"^ [0-9]\.", line): + blocks.append([]) + line = line[4:] + blocks[-1].append(line) + str_blocks = ["\n".join(b).strip() for b in blocks] + return [b for b in str_blocks if not b.endswith("Mappings:")] + + blocks_sym1 = split_blocks(sym1) + blocks_sym2 = split_blocks(sym2) + assert len(blocks_sym1) == len(blocks_sym2) + for block in blocks_sym1: + assert block in blocks_sym2 + + +class TestProjector(unittest.TestCase): + def base_test(self, case, kind, n_core, n_virt): + state = cache.adc_states[case]["adc3"][kind] + mospaces = state.reference_state.mospaces + + partitioning = SubspacePartitioning(mospaces, + core_orbitals=n_core, + outer_virtuals=n_virt) + spaces, nonzero_blocks = construct_nonzero_blocks(mospaces, n_core, n_virt) + + # Singles + proj = Projector(["o1", "v1"], partitioning, spaces["ph"]) + rhs = state.excitation_vector[0].ph.copy().set_random() + out = proj @ rhs + assert_equal_symmetry(rhs, out) + assert_nonzero_blocks(rhs, out, nonzero_blocks["ph"]) + + # + # Doubles + partitioning = SubspacePartitioning(mospaces, + core_orbitals=n_core, + outer_virtuals=n_virt) + proj = Projector(["o1", "o1", "v1", "v1"], partitioning, spaces["pphh"]) + rhs = state.excitation_vector[0].pphh.copy().set_random() + out = proj @ rhs + assert_equal_symmetry(rhs, out) + assert_nonzero_blocks(rhs, out, nonzero_blocks["pphh"]) + + def test_h2o_sto3g_singlet(self): + self.base_test("h2o_sto3g", "singlet", n_core=2, n_virt=1) + + def test_h2o_sto3g_triplet(self): + self.base_test("h2o_sto3g", "triplet", n_core=1, n_virt=1) + + def test_cn_sto3g(self): + self.base_test("cn_sto3g", "state", n_core=2, n_virt=1) + + def test_h2o_def2tzvp_singlet(self): + self.base_test("h2o_def2tzvp", "singlet", n_core=2, n_virt=5) + + def test_h2o_def2tzvp_triplet(self): + self.base_test("h2o_def2tzvp", "triplet", n_core=1, n_virt=3) + + def test_cn_ccpvdz(self): + self.base_test("cn_ccpvdz", "state", n_core=2, n_virt=4) + + +testcases = [("h2o_sto3g", "singlet"), ("h2o_sto3g", "triplet"), + ("cn_sto3g", "any"), ("h2o_def2tzvp", "singlet"), + ("h2o_def2tzvp", "triplet"), ("cn_ccpvdz", "any")] + + +@expand_test_templates(testcases) +class TestCvsTransfer(unittest.TestCase): + def template_high_level(self, case, kind): + kindkey = kind if kind != "any" else "state" + print(cache.adc_states[case]["cvs-adc2x"].keys()) + state_cvs = cache.adc_states[case]["cvs-adc2x"][kindkey] + matrix = adcc.AdcMatrix("adc2x", cache.refstate[case]) + + orth = np.array([[v @ w for v in state_cvs.excitation_vector] + for w in state_cvs.excitation_vector]) + fullvecs = transfer_cvs_to_full(state_cvs, matrix) + orthfull = np.array([[v @ w for v in fullvecs] for w in fullvecs]) + assert_allclose(orth, orthfull, atol=1e-16) + + def template_random(self, case, kind): + kindkey = kind if kind != "any" else "state" + state_cvs = cache.adc_states[case]["cvs-adc2x"][kindkey] + vectors = [v.copy().set_random() for v in state_cvs.excitation_vector] + matrix = adcc.AdcMatrix("adc2x", cache.refstate[case]) + + orth = np.array([[v @ w for v in vectors] for w in vectors]) + fullvecs = transfer_cvs_to_full(state_cvs.matrix, matrix, vectors, kind) + orthfull = np.array([[v @ w for v in fullvecs] for w in fullvecs]) + assert_allclose(orth, orthfull, atol=1e-16) diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py index 95887f32..1f0e37f1 100644 --- a/adcc/testdata/cache.py +++ b/adcc/testdata/cache.py @@ -47,7 +47,7 @@ def make_mock_adc_state(refstate, matmethod, kind, reference): state.method = matrix.method state.ground_state = ground_state state.reference_state = refstate - state.kind = kind + state.kind = kind if kind != "state" else "any" state.eigenvalues = reference[kind]["eigenvalues"][:n_full] spin_change = 0 diff --git a/examples/2_core_excitations/water_core_by_projection.py b/examples/2_core_excitations/water_core_by_projection.py new file mode 100755 index 00000000..38bf86a7 --- /dev/null +++ b/examples/2_core_excitations/water_core_by_projection.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +import adcc +import adcc.projection + +from pyscf import gto, scf +from adcc.AdcMatrix import AdcMatrixProjected + +# Aim of this script is to compute core excitations of water +# using a projection approach, i.e. where the non-CVS ADC matrix +# is projected into parts of the orbital subspaces to target core +# excitations specifically. + +# Run SCF in pyscf +mol = gto.M( + atom='O 0 0 0;' + 'H 0 0 1.795239827225189;' + 'H 1.693194615993441 0 -0.599043184453037', + basis='cc-pvtz', + unit="Bohr" +) +scfres = scf.RHF(mol) +scfres.conv_tol = 1e-12 +scfres.conv_tol_grad = 1e-9 +scfres.kernel() +print(adcc.banner()) + +# Construct the standard ADC(2)-x matrix and project it, such that +# only cv, ccvv and ocvv excitations are allowed. Choose one orbital (oxygen 1s) +# to make up the core space. +matrix = adcc.AdcMatrix("adc2x", adcc.ReferenceState(scfres)) +pmatrix = AdcMatrixProjected(matrix, ["cv", "ccvv", "ocvv"], core_orbitals=1) + +# Run the ADC calculation. Note that doubles guesses are not yet supported in this +# setup and like lead to numerical issues (zero excitations) if selected. So we +# explicitly disable them here. +state_proj = adcc.run_adc(pmatrix, n_singlets=7, conv_tol=1e-8, n_guesses_doubles=0) +print(state_proj.describe()) + +# For comparison we also run a standard CVS-ADC(2)-x calculation: +state_cvs = adcc.cvs_adc2x(scfres, n_singlets=7, conv_tol=1e-8, core_orbitals=1) +print(state_cvs.describe()) + +# Note that this calculation could also be used as an initial guess: +guesses = adcc.projection.transfer_cvs_to_full(state_cvs, pmatrix) +state_proj2 = adcc.run_adc(pmatrix, n_singlets=7, conv_tol=1e-8, guesses=guesses)