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)