Merge pull request #21 from adc-connect/integration-psi4

Changes for integration with Psi4

LICENSE_adccore

@@ -2,7 +2,7 @@
 #--  Introduction  --#
 ######################
 adcc and the adcc source code is released under the terms of the GNU General Public
- License Version 3 (LGPLv3). This terms of release, however, do not apply to the
+ License Version 3 (GPLv3). This terms of release, however, do not apply to the
 adccore binary file `libadccore.so` (or a simlilarly named dylib file for MacOS),
 contained in the `adcc/lib` directory of the python code tarball or the python
 package directory.
@@ -119,19 +119,6 @@ them, and agree to be bound by their terms and conditions.
 This section lists the third-party components integrated into adccore in binary
 form, provides the copyright holders and their terms of use.
 
-#
-# ctx
-#
-- **Description:** Library providing key-value C++ datastructures
-                   for organised hierarchical storage.
-- **Copyright:**   Copyright (c) 2019 Michael F. Herbst
-- **License:**     Apache License 2.0
-- **Website:**     https://github.com/mfherbst/ctx
-- **Publication:** http://doi.org/10.5281/zenodo.2590706
-
-To view the license statement and copyright notices, see the files
-in the folder adcc/lib/libadccore_thirdparty/ctx
-
 #
 # libadc
 #

MANIFEST.in

@@ -14,11 +14,9 @@ include extension/AdcCore.py
 # (i.e. include files and compiled binary)
 include extension/adccore/adccore_config.json
 recursive-include extension/adccore/include/adcc/ *.hh
-recursive-include extension/adccore/include/ctx/ *.hh
 include adcc/lib/libadccore.so
 include adcc/lib/libadccore.dylib
 include adcc/lib/libadccore.*.dylib
 include adcc/lib/libstdc++.so.*
 include adcc/lib/libc++.so.*
 include adcc/lib/libadccore_LICENSE
-recursive-include adcc/lib/libadccore_thirdparty/ctx/ *

adcc/ExcitedStates.py

@@ -41,7 +41,7 @@
 
 
 class FormatExcitationVector:
-    def __init__(self, matrix, tolerance=1e-3, index_format=None):
+    def __init__(self, matrix, tolerance=0.01, index_format=None):
         """
         Set up a formatter class for formatting excitation vectors.
 
@@ -471,7 +471,7 @@ def _repr_pretty_(self, pp, cycle):
         else:
             pp.text(self.describe())
 
-    def describe_amplitudes(self, tolerance=1e-3, index_format=None):
+    def describe_amplitudes(self, tolerance=0.01, index_format=None):
         """
         Return a string describing the dominant amplitudes of each
         excitation vector in human-readable form. The ``kwargs``

adcc/FormatDominantElements.py

@@ -29,9 +29,9 @@
 
 
 class FormatDominantElements:
-    def __init__(self, mospaces, tolerance=1e-3, index_format=FormatIndexAdcc):
+    def __init__(self, mospaces, tolerance=0.01, index_format=FormatIndexAdcc):
         self.mospaces = mospaces
-        self.tolerance = 1e-3
+        self.tolerance = tolerance
         self.value_format = "{:+8.3g}"  # Formatting used for the values
 
         if isinstance(index_format, type):

adcc/FormatIndex.py

@@ -245,8 +245,7 @@ def optimise_formatting(self, space_index_pairs):
         """
         maxlen = max(len(self._translate_index(space, idx)[1])
                      for space, idx in space_index_pairs)
-        log_max_idx = int(np.log(max(1, maxlen)))
-        self.maxlen_offset = max(log_max_idx, self.maxlen_offset, 2)
+        self.maxlen_offset = max(maxlen, self.maxlen_offset, 2)
 
     def format(self, space, idx, concat_spin=True):
         word, offset, spin = self._translate_index(space, idx)

adcc/LazyMp.py

@@ -81,6 +81,21 @@ def dipole_moment(self, level=2):
             raise NotImplementedError("Only dipole moments for level 1 and 2"
                                       " are implemented.")
 
+    def energy(self, level=2):
+        """
+        Obtain the total energy (SCF energy plus all corrections)
+        at a particular level of perturbation theory.
+        """
+        if level == 0:
+            # Sum of orbital energies ...
+            raise NotImplementedError("Total MP(0) energy not implemented.")
+
+        # Accumulator for all energy terms
+        energies = [self.reference_state.energy_scf]
+        for il in range(2, level + 1):
+            energies.append(self.energy_correction(il))
+        return sum(energies)
+
     @property
     def mp2_density(self):
         return self.density(2)

adcc/backends/available_backends.py → adcc/backends/InvalidReference.py

@@ -22,33 +22,10 @@
 ## ---------------------------------------------------------------------
 
 
-def is_module_available(module):
-    import importlib
-
-    try:
-        importlib.import_module(module)
-        return True
-    except ImportError:
-        return False
-
-
-status = {
-    "pyscf": is_module_available("pyscf"),
-    "psi4": is_module_available("psi4") and is_module_available("psi4.core"),
-    "veloxchem": is_module_available("veloxchem"),
-    "molsturm": is_module_available("molsturm"),
-}
-
-
-available = sorted([b for b in status if status[b]])
-
-
-def first_available():
-    if len(available) == 0:
-        raise RuntimeError("No backend available.")
-    else:
-        return available[0]
-
-
-def have_backend(backend):
-    return status.get(backend, False)
+class InvalidReference(ValueError):
+    """
+    Exception thrown if a passed SCF reference is invalid, e.g. because
+    a feature like density-fitting has been applied, which is inconsistent
+    with the current capabilities of adcc.
+    """
+    pass

adcc/backends/__init__.py

@@ -23,11 +23,66 @@
 import os
 import warnings
 
-from .available_backends import available, first_available, have_backend
+from pkg_resources import parse_version
 
 import h5py
 
-__all__ = ["import_scf_results", "run_hf", "have_backend", "available"]
+from .InvalidReference import InvalidReference
+
+__all__ = ["import_scf_results", "run_hf", "have_backend", "available",
+           "InvalidReference"]
+
+
+def is_module_available(module, min_version=None):
+    """Check using importlib if a module is available."""
+    import importlib
+
+    try:
+        mod = importlib.import_module(module)
+    except ImportError:
+        return False
+
+    if not min_version:  # No version check
+        return True
+
+    if not hasattr(mod, "__version__"):
+        warnings.warn(
+            "Could not check host program {} minimal version, "
+            "since __version__ tag not found. Proceeding anyway."
+            "".format(module)
+        )
+        return True
+
+    if parse_version(mod.__version__) < parse_version(min_version):
+        warnings.warn(
+            "Found host program module {}, but its version {} is below "
+            "the least required (== {}). This host program will be ignored."
+            "".format(module, mod.__version__, min_version)
+        )
+        return False
+    return True
+
+
+# Cache for the list of available backends ... cannot be filled right now,
+# since this can lead to import loops when adcc is e.g. used from Psi4
+__status = dict()
+
+
+def available():
+    global __status
+    if not __status:
+        status = {
+            "pyscf": is_module_available("pyscf", "1.5.0"),
+            "psi4": is_module_available("psi4", "1.2.1") and is_module_available("psi4.core"),
+            "veloxchem": is_module_available("veloxchem"),  # Exports no version info
+            "molsturm": is_module_available("molsturm"),    # Exports no version info
+        }
+    return sorted([b for b in status if status[b]])
+
+
+def have_backend(backend):
+    """Is a particular backend available?"""
+    return backend in available()
 
 
 def import_scf_results(res):
@@ -100,10 +155,20 @@ def run_hf(backend=None, xyz=None, basis="sto-3g", charge=0, multiplicity=1,
         conv_tol:       energy convergence tolerance
         conv_tol_grad:  convergence tolerance of the electronic gradient
         max_iter:       maximum number of SCF iterations
+
+        Note: This function only exists for testing purposes and should
+        not be used in production calculations.
     """
 
     if not backend:
-        backend = first_available()
+        if len(available()) == 0:
+            raise RuntimeError(
+                "No supported host-program available as SCF backend. "
+                "See https://adc-connect.org/installation.html#install-hostprogram "
+                "for installation instructions."
+            )
+        else:
+            backend = available()[0]
         warnings.warn("No backend specified. Using {}.".format(backend))
 
     if not have_backend(backend):

adcc/backends/psi4.py

@@ -23,7 +23,9 @@
 import psi4
 import numpy as np
 
+from .InvalidReference import InvalidReference
 from .eri_build_helper import EriBuilder
+
 from adcc.misc import cached_property
 
 from libadcc import HartreeFockProvider
@@ -33,7 +35,7 @@ class Psi4OperatorIntegralProvider:
     def __init__(self, wfn):
         self.wfn = wfn
         self.backend = "psi4"
-        self.mints = psi4.core.MintsHelper(self.wfn.basisset())
+        self.mints = psi4.core.MintsHelper(self.wfn)
 
     @cached_property
     def electric_dipole(self):
@@ -43,7 +45,7 @@ def electric_dipole(self):
 class Psi4EriBuilder(EriBuilder):
     def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta):
         self.wfn = wfn
-        self.mints = psi4.core.MintsHelper(self.wfn.basisset())
+        self.mints = psi4.core.MintsHelper(self.wfn)
         super().__init__(n_orbs, n_orbs_alpha, n_alpha, n_beta)
 
     @property
@@ -72,7 +74,7 @@ def __init__(self, wfn):
         super().__init__()
 
         if not isinstance(wfn, psi4.core.RHF):
-            raise TypeError("Only restricted references (RHF) are supported.")
+            raise InvalidReference("Only restricted references (RHF) are supported.")
 
         self.wfn = wfn
         self.eri_ffff = None
@@ -157,33 +159,46 @@ def flush_cache(self):
         self.eri_cache = None
 
 
-def import_scf(scfdrv):
-    if not isinstance(scfdrv, psi4.core.HF):
-        raise TypeError("Unsupported type for backends.psi4.import_scf.")
-
-    if not isinstance(scfdrv, psi4.core.RHF):
-        raise TypeError("Only restricted references (RHF) are supported.")
-
-    # TODO: Psi4 throws an exception if SCF is not converged
-    # and there is, to the best of my knowledge, no `is_converged` property
-    # of the psi4 wavefunction
-    # if not scfdrv.is_converged:
-    #     raise ValueError("Cannot start an adc calculation on top of an SCF, "
-    #                      "which is not converged.")
+def import_scf(wfn):
+    if not isinstance(wfn, psi4.core.HF):
+        raise InvalidReference(
+            "Only psi4.core.HF and its subtypes are supported references in "
+            "backends.psi4.import_scf. This indicates that you passed an unsupported "
+            "SCF reference. Make sure you did a restricted or unrestricted HF "
+            "calculation."
+        )
 
-    provider = Psi4HFProvider(scfdrv)
+    if not isinstance(wfn, psi4.core.RHF):
+        raise InvalidReference("Right now only restricted references (RHF) are "
+                               "supported for Psi4.")
+
+    # TODO This is not fully correct, because the core.Wavefunction object
+    #      has an internal, but py-invisible Options structure, which contains
+    #      the actual set of options ... theoretically they could differ
+    scf_type = psi4.core.get_global_option('SCF_TYPE')
+    unsupported_scf_types = ["CD", "DISK_DF", "MEM_DF"]  # Choleski or density-fitting
+    if scf_type in unsupported_scf_types:
+        raise InvalidReference(f"Unsupported Psi4 SCF_TYPE, should not be one "
+                               "of {unsupported_scf_types}")
+
+    if wfn.nirrep() > 1:
+        raise InvalidReference("The passed Psi4 wave function object needs to have "
+                               "exactly one irrep, i.e. be of C1 symmetry.")
+
+    # Psi4 throws an exception if SCF is not converged, so there is no need
+    # to assert that here.
+    provider = Psi4HFProvider(wfn)
     return provider
 
 
-basissets = {
-    "sto3g": "sto-3g",
-    "def2tzvp": "def2-tzvp",
-    "ccpvdz": "cc-pvdz",
-}
-
-
 def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
            conv_tol_grad=1e-8, max_iter=150):
+    basissets = {
+        "sto3g": "sto-3g",
+        "def2tzvp": "def2-tzvp",
+        "ccpvdz": "cc-pvdz",
+    }
+
     mol = psi4.geometry("""
         {charge} {multiplicity}
         {xyz}