Move MTMs, diffdms and transition dms to python

.travis.yml

@@ -18,8 +18,7 @@ addons:
       - libopenblas-base
 
 install:
-  - pip3 install pyscf
-  - pip3 install cppe
+  - pip3 install pyscf cppe
   - pip3 install --verbose -r requirements.txt
 script: python3 setup.py test -a '--cov=adcc'
 

adcc/ExcitedStates.py

@@ -23,13 +23,14 @@
 import warnings
 import numpy as np
 
+from . import adc_pp
 from .misc import cached_property
 from .timings import Timer, timed_member_call
 from .AdcMethod import AdcMethod
+from .Excitation import Excitation, mark_excitation_property
 from .FormatIndex import (FormatIndexAdcc, FormatIndexBase,
                           FormatIndexHfProvider, FormatIndexHomoLumo)
 from .visualisation import ExcitationSpectrum
-from .state_densities import compute_gs2state_optdm, compute_state_diffdm
 from .OneParticleOperator import product_trace
 from .FormatDominantElements import FormatDominantElements
 
@@ -38,7 +39,6 @@
 from scipy import constants
 from matplotlib import pyplot as plt
 from .solver.SolverStateBase import EigenSolverStateBase
-from .Excitation import Excitation, mark_excitation_property
 
 
 class FormatExcitationVector:
@@ -285,8 +285,8 @@ def excitation_vector(self):
     @timed_member_call(timer="_property_timer")
     def transition_dm(self):
         """List of transition density matrices of all computed states"""
-        return [compute_gs2state_optdm(self.property_method, self.ground_state,
-                                       evec, self.matrix.intermediates)
+        return [adc_pp.transition_dm(self.property_method, self.ground_state,
+                                     evec, self.matrix.intermediates)
                 for evec in self.excitation_vector]
 
     @cached_property
@@ -369,8 +369,8 @@ def rotatory_strength(self):
     @timed_member_call(timer="_property_timer")
     def state_diffdm(self):
         """List of difference density matrices of all computed states"""
-        return [compute_state_diffdm(self.property_method, self.ground_state,
-                                     evec, self.matrix.intermediates)
+        return [adc_pp.state_diffdm(self.property_method, self.ground_state,
+                                    evec, self.matrix.intermediates)
                 for evec in self.excitation_vector]
 
     @property

adcc/OneParticleOperator.py

@@ -26,15 +26,15 @@
 
 
 class OneParticleOperator(libadcc.OneParticleOperator):
-    def __init__(self, mospaces, is_symmetric=True, cartesian_transform="1"):
+    def __init__(self, spaces, is_symmetric=True, cartesian_transform="1"):
         """
         Construct an OneParticleOperator object. All blocks are initialised
         as zero blocks.
 
         Parameters
         ----------
 
-        mospaces : adcc.Mospaces
+        spaces : adcc.Mospaces or adcc.ReferenceState or adcc.LazyMp
             MoSpaces object
 
         is_symmetric : bool
@@ -43,7 +43,14 @@ def __init__(self, mospaces, is_symmetric=True, cartesian_transform="1"):
         cartesian_transform : str
             Symbol for cartesian transformation (see make_symmetry for details.)
         """
-        super().__init__(mospaces, is_symmetric, cartesian_transform)
+        if isinstance(spaces, libadcc.ReferenceState):
+            super().__init__(spaces.mospaces, is_symmetric, cartesian_transform)
+            self.reference_state = spaces
+        elif isinstance(spaces, libadcc.LazyMp):
+            super().__init__(spaces.mospaces, is_symmetric, cartesian_transform)
+            self.reference_state = spaces.reference_state
+        else:
+            super().__init__(spaces, is_symmetric, cartesian_transform)
 
     @classmethod
     def from_cpp(cls, cpp_operator):
@@ -170,6 +177,11 @@ def __mul__(self, other):
     def __rmul__(self, other):
         return self.copy().__imul__(other)
 
+    def evaluate(self):
+        for b in self.blocks_nonzero:
+            self.block(b).evaluate()
+        return self
+
 
 def product_trace(op1, op2):
     # TODO use blocks_nonzero and build the set intersection

adcc/ReferenceState.py

@@ -22,15 +22,15 @@
 ## ---------------------------------------------------------------------
 import numpy as np
 
-import libadcc
-
 from .misc import cached_property
 from .Tensor import Tensor
 from .MoSpaces import MoSpaces
 from .backends import import_scf_results
 from .OperatorIntegrals import OperatorIntegrals
 from .OneParticleOperator import OneParticleOperator, product_trace
 
+import libadcc
+
 __all__ = ["ReferenceState"]
 
 
@@ -159,6 +159,14 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
             if hasattr(hfdata, name):
                 setattr(self, name, getattr(hfdata, name))
 
+    def __getattr__(self, attr):
+        from . import block as b
+
+        if attr.startswith("f"):
+            return self.fock(b.__getattr__(attr[1:]))
+        else:
+            return self.eri(b.__getattr__(attr))
+
     @property
     def mospaces(self):
         return self._mospaces
@@ -187,9 +195,9 @@ def density(self):
         Return the Hartree-Fock density in the MO basis
         """
         density = OneParticleOperator(self.mospaces, is_symmetric=True)
-        for b in density.blocks:
-            sym = libadcc.make_symmetry_operator(self.mospaces, b, True, "1")
-            density.set_block(b, Tensor(sym))
+        for block in density.blocks:
+            sym = libadcc.make_symmetry_operator(self.mospaces, block, True, "1")
+            density.set_block(block, Tensor(sym))
         for ss in self.mospaces.subspaces_occupied:
             density[ss + ss].set_mask("ii", 1)
         density.reference_state = self

adcc/adc_pp/__init__.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from .state_diffdm import state_diffdm
+from .transition_dm import transition_dm
+from .state2state_transition_dm import state2state_transition_dm
+from .modified_transition_moments import modified_transition_moments
+
+"""
+Submodule, which contains rather lengthy low-level kernels
+(e.g. matrix-vector products or working equations), which are called
+from the high-level objects in the adcc main module.
+"""
+
+__all__ = ["state_diffdm", "state2state_transition_dm", "transition_dm",
+           "modified_transition_moments"]

adcc/adc_pp/modified_transition_moments.py

@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from math import sqrt
+
+from adcc import block as b
+from adcc.AdcMethod import AdcMethod
+from adcc.functions import einsum, evaluate
+from adcc.AmplitudeVector import AmplitudeVector
+
+import libadcc
+
+
+def mtm_adc0(mp, dipop, intermediates):
+    return AmplitudeVector(dipop[b.ov])
+
+
+def mtm_adc1(mp, dipop, intermediates):
+    f1 = dipop[b.ov] - einsum("ijab,jb->ia", mp.t2(b.oovv), dipop[b.ov])
+    return AmplitudeVector(f1)
+
+
+def mtm_adc2(mp, dipop, intermediates):
+    t2 = mp.t2(b.oovv)
+    p0 = mp.mp2_diffdm
+
+    f1 = (
+        + dipop[b.ov]
+        - einsum("ijab,jb->ia", t2,
+                 + dipop[b.ov] - 0.5 * einsum("jkbc,kc->jb", t2, dipop[b.ov]))
+        + 0.5 * einsum("ij,ja->ia", p0[b.oo], dipop[b.ov])
+        - 0.5 * einsum("ib,ab->ia", dipop[b.ov], p0[b.vv])
+        + einsum("ib,ab->ia", p0[b.ov], dipop[b.vv])
+        - einsum("ij,ja->ia", dipop[b.oo], p0[b.ov])
+        - einsum("ijab,jb->ia", mp.td2(b.oovv), dipop[b.ov])
+    )
+    f2 = (
+        + einsum("ijac,bc->ijab", t2, dipop[b.vv]).antisymmetrise(2, 3)
+        - einsum("ik,kjab->ijab", dipop[b.oo], t2).antisymmetrise(0, 1)
+    )
+    return AmplitudeVector(f1, f2)
+
+
+def mtm_cvs_adc0(mp, dipop, intermediates):
+    return AmplitudeVector(dipop[b.cv])
+
+
+def mtm_cvs_adc2(mp, dipop, intermediates):
+    f1 = (
+        + dipop[b.cv]
+        - einsum("Ib,ba->Ia", dipop[b.cv], intermediates.cv_p_vv)
+        - einsum("Ij,ja->Ia", dipop[b.co], intermediates.cv_p_ov)
+    )
+    f2 = (1 / sqrt(2)) * einsum("Ik,kjab->jIab", dipop[b.co], mp.t2(b.oovv))
+    return AmplitudeVector(f1, f2)
+
+
+DISPATCH = {
+    "adc0": mtm_adc0,
+    "adc1": mtm_adc1,
+    "adc2": mtm_adc2,
+    "cvs-adc0": mtm_cvs_adc0,
+    "cvs-adc1": mtm_cvs_adc0,  # Identical to CVS-ADC(0)
+    "cvs-adc2": mtm_cvs_adc2,
+}
+
+
+def modified_transition_moments(method, ground_state, dipole_operator=None,
+                                intermediates=None):
+    """Compute the modified transition moments (MTM) for the provided
+    ADC method with reference to the passed ground state.
+
+    Parameters
+    ----------
+    method: adc.Method
+        Provide a method at which to compute the MTMs
+    ground_state : adcc.LazyMp
+        The MP ground state
+    dipole_operator : adcc.OneParticleOperator or list, optional
+        Only required if different dipole operators than the standard
+        dipole operators in the MO basis should be used.
+    intermediates : AdcIntermediates
+        Intermediates from the ADC calculation to reuse
+
+    Returns
+    -------
+    adcc.AmplitudeVector or list of adcc.AmplitudeVector
+    """
+    if not isinstance(method, AdcMethod):
+        method = AdcMethod(method)
+    if not isinstance(ground_state, libadcc.LazyMp):
+        raise TypeError("ground_state should be a LazyMp object.")
+    if intermediates is None:
+        intermediates = libadcc.AdcIntermediates(ground_state)
+
+    unpack = False
+    if dipole_operator is None:
+        dipole_operator = ground_state.reference_state.operators.electric_dipole
+    elif not isinstance(dipole_operator, list):
+        unpack = True
+        dipole_operator = [dipole_operator]
+    if method.name not in DISPATCH:
+        raise NotImplementedError("modified_transition_moments is not "
+                                  f"implemented for {method.name}.")
+
+    ret = [DISPATCH[method.name](ground_state, dipop, intermediates)
+           for dipop in dipole_operator]
+    if unpack:
+        assert len(ret) == 1
+        ret = ret[0]
+    return evaluate(ret)

adcc/adc_pp/state2state_transition_dm.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2018 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+from adcc.AdcMethod import AdcMethod
+from adcc.AmplitudeVector import AmplitudeVector
+
+import libadcc
+
+DISPATCH = {}  # None implemented
+
+
+def state2state_transition_dm(method, ground_state, amplitude_from,
+                              amplitude_to, intermediates=None):
+    """
+    Compute the state to state transition density matrix
+    state in the MO basis using the intermediate-states representation.
+
+    Parameters
+    ----------
+    method : str, AdcMethod
+        The method to use for the computation (e.g. "adc2")
+    ground_state : LazyMp
+        The ground state upon which the excitation was based
+    amplitude_from : AmplitudeVector
+        The amplitude vector of the state to start from
+    amplitude_to : AmplitudeVector
+        The amplitude vector of the state to excite to
+    intermediates : AdcIntermediates
+        Intermediates from the ADC calculation to reuse
+    """
+    if not isinstance(method, AdcMethod):
+        method = AdcMethod(method)
+    if not isinstance(ground_state, libadcc.LazyMp):
+        raise TypeError("ground_state should be a LazyMp object.")
+    if not isinstance(amplitude_from, AmplitudeVector):
+        raise TypeError("amplitude_from should be an AmplitudeVector object.")
+    if not isinstance(amplitude_to, AmplitudeVector):
+        raise TypeError("amplitude_to should be an AmplitudeVector object.")
+    if intermediates is None:
+        intermediates = libadcc.AdcIntermediates(ground_state)
+
+    if method.name not in DISPATCH:
+        raise NotImplementedError("state2state_transition_dm is not implemented "
+                                  f"for {method.name}.")
+    else:
+        ret = DISPATCH[method.name](ground_state, amplitude_from, amplitude_to,
+                                    intermediates)
+        return ret.evaluate()