Source code for psi4.driver.procrouting.empirical_dispersion

#
# @BEGIN LICENSE
#
# Psi4: an open-source quantum chemistry software package
#
# Copyright (c) 2007-2023 The Psi4 Developers.
#
# The copyrights for code used from other parties are included in
# the corresponding files.
#
# This file is part of Psi4.
#
# Psi4 is free software; you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, version 3.
#
# Psi4 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 Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Psi4; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# @END LICENSE
#

import collections
from typing import Dict, List, Tuple, Union

import numpy as np
from qcelemental.models import AtomicInput
import qcengine as qcng

from psi4 import core
from psi4.driver import p4util
from psi4.driver import driver_findif
from psi4.driver.p4util.exceptions import ValidationError

_engine_can_do = collections.OrderedDict([
    # engine order establishes default for each disp
    ("libdisp",  ["d1", "d2",                                                                                                               "chg", "das2009", "das2010",]),
    ("s-dftd3",  [            "d3zero2b", "d3bj2b", "d3mzero2b", "d3mbj2b", "d3zeroatm", "d3bjatm", "d3mzeroatm", "d3mbjatm",                                           ]),
    ("dftd3",    [      "d2", "d3zero2b", "d3bj2b", "d3mzero2b", "d3mbj2b",                                                                                             ]),
    ("nl",       [                                                                                                                          "nl",                       ]),
    ("mp2d",     [                                                                                                                          "dmp2",                     ]),
    ("dftd4",    [                                                                                                            "d4bjeeqatm",                             ]),
    ("mctc-gcp", [                                                                                                                          "3c",                       ]),
    ("gcp",      [                                                                                                                          "3c",                       ]),
]) # yapf: disable


def _capable_engines_for_disp()-> Dict[str, List[str]]:
    """Invert _engine_can_do dictionary and check program detection.

    Returns a dictionary with keys all dispersion levels and values a list of all
    capable engines, where the engine in the first element is available, if any are.

    """
    try:
        from qcengine.testing import _programs as _programs_qcng
    except ModuleNotFoundError:
        # _programs_qcng is up-to-date with current harnesses but it requires pytest present, so let's provide a workaround
        from qcelemental.util import which, which_import
        _programs_qcng = {
            "dftd3": which("dftd3", return_bool=True),
            "dftd4": which_import("dftd4", return_bool=True),
            "s-dftd3": which_import("dftd3", return_bool=True),
            "mctc-gcp": which("mctc-gcp", return_bool=True),
            "gcp": which("gcp", return_bool=True),
            "mp2d": which("mp2d", return_bool=True),
        }

    programs_disp = {k: v for k, v in _programs_qcng.items() if k in _engine_can_do}
    programs_disp["libdisp"] = True
    programs_disp["nl"] = True

    capable = collections.defaultdict(list)
    capable_sorted_by_available = collections.defaultdict(list)
    for eng, disps in _engine_can_do.items():
        for disp in disps:
            capable[disp].append(eng)
    for disp, engines in capable.items():
        capable_sorted_by_available[disp] = sorted(engines, key=lambda x: (not programs_disp[x], x))

    return capable_sorted_by_available


[docs] class EmpiricalDispersion(): """Lightweight unification of empirical dispersion calculation modes. Attributes ---------- dashlevel : str {"d1", "d2", "chg", "das2009", "das2010", "nl", "dmp2", "d3zero2b", "d3bj2b", "d3mzero2b", "d3mbj2b", "d3zeroatm", "d3bjatm", "d3mzeroatm", "d3mbjatm", "d4bjeeqatm"} Name of dispersion correction to be applied. Resolved from `name_hint` and/or `level_hint` into a key of `empirical_dispersion_resources.dashcoeff`. dashparams : dict Complete set of parameter values defining the flexible parts of :py:attr:`dashlevel`. Number and parameter names vary by :py:attr:`dashlevel`. Resolved into a complete set (keys of dashcoeff[dashlevel]['default']) from `name_hint` and/or `dashcoeff_supplement` and/or user `param_tweaks`. fctldash : str If :py:attr:`dashparams` for :py:attr:`dashlevel` corresponds to a defined, named, untweaked "functional-dashlevel" set, then that functional. Otherwise, empty string. description : str Tagline for dispersion :py:attr:`dashlevel`. dashlevel_citation : str Literature reference for dispersion :py:attr:`dashlevel` in general, *not necessarily* for :py:attr:`dashparams`. dashparams_citation : str Literature reference for dispersion parameters, if :py:attr:`dashparams` corresponds to a defined, named, untweaked "functional-dashlevel" set with a citation. Otherwise, empty string. dashcoeff_supplement : dict See description in `qcengine.programs.empirical_dispersion_resources.from_arrays`. Used here to "bless" the dispersion definitions attached to the procedures/dft/<rung>_functionals-defined dictionaries as legit, non-custom, and of equal validity to `qcengine.programs.empirical_dispersion_resources.dashcoeff` itself for purposes of validating :py:attr:`fctldash`. engine : str {'libdisp', "s-dftd3", 'dftd3', 'nl', 'mp2d', "dftd4"} Compute engine for dispersion. One of Psi4's internal libdisp library, external Grimme or Beran projects, or nl. disp : Dispersion Only present for :py:attr:`engine` `=libdisp`. Psi4 class instance prepared to compute dispersion. ordered_params : list Fixed-order list of relevant parameters for :py:attr:`dashlevel`. Matches :rst:psivar:`DISPERSION CORRECTION ENERGY` ordering. Used for printing. Parameters ---------- name_hint Name of functional (func only, func & disp, or disp only) for which to compute dispersion (e.g., blyp, BLYP-D2, blyp-d3bj, blyp-d3(bj), hf+d). Any or all parameters initialized from ``dashcoeff[dashlevel][functional-without-dashlevel]`` or ``dashcoeff_supplement[dashlevel][functional-with-dashlevel]`` can be overwritten via `param_tweaks`. level_hint Name of dispersion correction to be applied (e.g., d, D2, d3(bj), das2010). Must be key in `dashcoeff` or "alias" or "formal" to one. param_tweaks Values for the same keys as `dashcoeff[dashlevel]['default']` (and same order if list) used to override any or all values initialized by `name_hint`. Extra parameters will error. engine Override which code computes dispersion. See above for allowed values. Formerly (pre Nov 2022) only relevant for -D2, which can be computed by libdisp or dftd3. Now (post Nov 2022) also relevant for -D3 variants, which can be computed by dftd3 executable or simple-dftd3 Python module. gcp_engine Override which code computes the gcp correction. Now can use classic gcp or mctc-gcp executables. save_pairwise_disp Whether to request atomic pairwise analysis. """ def __init__(self, *, name_hint: str = None, level_hint: str = None, param_tweaks: Union[Dict, List] = None, engine: str = None, gcp_engine: str = None, save_pairwise_disp: bool = False): from .dft import dashcoeff_supplement self.dashcoeff_supplement = dashcoeff_supplement self.save_pairwise_disp = save_pairwise_disp resolved = qcng.programs.empirical_dispersion_resources.from_arrays( name_hint=name_hint, level_hint=level_hint, param_tweaks=param_tweaks, dashcoeff_supplement=self.dashcoeff_supplement) self.fctldash = resolved['fctldash'] self.dashlevel = resolved['dashlevel'] self.dashparams = resolved['dashparams'] self.description = qcng.programs.empirical_dispersion_resources.dashcoeff[self.dashlevel]['description'] self.ordered_params = qcng.programs.empirical_dispersion_resources.dashcoeff[self.dashlevel]['default'].keys() self.dashlevel_citation = qcng.programs.empirical_dispersion_resources.dashcoeff[self.dashlevel]['citation'] self.dashparams_citation = resolved['dashparams_citation'] capable_engines_for_disp = _capable_engines_for_disp() if engine is None: self.engine = capable_engines_for_disp[self.dashlevel][0] else: if self.dashlevel in _engine_can_do[engine]: self.engine = engine else: raise ValidationError(f"This little engine ({engine}) can't ({self.dashlevel})") if self.engine == 'libdisp': self.disp = core.Dispersion.build(self.dashlevel, **resolved['dashparams']) if gcp_engine is None: self.gcp_engine = capable_engines_for_disp["3c"][0] else: if "3c" in _engine_can_do[gcp_engine]: self.gcp_engine = gcp_engine else: raise ValidationError(f"This little engine ({engine}) can't (3c)")
[docs] def print_out(self): """Format dispersion parameters of `self` for output file.""" text = [] text.append(" => {}: Empirical Dispersion <=".format( (self.fctldash.upper() if self.fctldash.upper() else 'Custom'))) text.append('') text.append(self.description) text.append(self.dashlevel_citation.rstrip()) if self.dashparams_citation: text.append(" Parametrisation from:{}".format(self.dashparams_citation.rstrip())) text.append('') for op in self.ordered_params: text.append(" %6s = %14.6f" % (op, self.dashparams[op])) text.append('\n') core.print_out('\n'.join(text))
[docs] def compute_energy(self, molecule: core.Molecule, wfn: core.Wavefunction = None) -> float: """Compute dispersion energy based on engine, dispersion level, and parameters in `self`. Parameters ---------- molecule System for which to compute empirical dispersion correction. wfn Location to set QCVariables Returns ------- float Dispersion energy [Eh]. Notes ----- :psivar:`DISPERSION CORRECTION ENERGY` Disp always set. Overridden in SCF finalization, but that only changes for "-3C" methods. :psivar:`fctl DISPERSION CORRECTION ENERGY` Set if :py:attr:`fctldash` nonempty. """ if self.engine in ["s-dftd3", 'dftd3', 'mp2d', "dftd4"]: resi = AtomicInput( **{ 'driver': 'energy', 'model': { 'method': self.fctldash, 'basis': '(auto)', }, 'keywords': { 'level_hint': self.dashlevel, 'params_tweaks': self.dashparams, 'dashcoeff_supplement': self.dashcoeff_supplement, 'pair_resolved': self.save_pairwise_disp, 'apply_qcengine_aliases': True, # for s-dftd3 'verbose': 1, }, 'molecule': molecule.to_schema(dtype=2), 'provenance': p4util.provenance_stamp(__name__), }) jobrec = qcng.compute( resi, self.engine, raise_error=True, task_config={"scratch_directory": core.IOManager.shared_object().get_default_path(), "ncores": core.get_num_threads()}) dashd_part = float(jobrec.extras['qcvars']['DISPERSION CORRECTION ENERGY']) if wfn is not None: for k, qca in jobrec.extras['qcvars'].items(): if ("CURRENT" not in k) and ("PAIRWISE" not in k): wfn.set_variable(k, float(qca) if isinstance(qca, str) else qca) # Pass along the pairwise dispersion decomposition if we need it if self.save_pairwise_disp is True: wfn.set_variable("PAIRWISE DISPERSION CORRECTION ANALYSIS", jobrec.extras['qcvars']["2-BODY PAIRWISE DISPERSION CORRECTION ANALYSIS"]) if self.fctldash in ['hf3c', 'pbeh3c']: jobrec = qcng.compute( resi, self.gcp_engine, raise_error=True, task_config={"scratch_directory": core.IOManager.shared_object().get_default_path(), "ncores": core.get_num_threads()}) gcp_part = jobrec.return_result dashd_part += gcp_part return dashd_part else: ene = self.disp.compute_energy(molecule) core.set_variable('DISPERSION CORRECTION ENERGY', ene) if self.fctldash: core.set_variable(f"{self.fctldash} DISPERSION CORRECTION ENERGY", ene) return ene
[docs] def compute_gradient(self, molecule: core.Molecule, wfn: core.Wavefunction = None) -> core.Matrix: """Compute dispersion gradient based on engine, dispersion level, and parameters in `self`. Parameters ---------- molecule System for which to compute empirical dispersion correction. wfn Location to set QCVariables Returns ------- Matrix (nat, 3) dispersion gradient [Eh/a0]. """ if self.engine in ["s-dftd3", 'dftd3', 'mp2d', "dftd4"]: resi = AtomicInput( **{ 'driver': 'gradient', 'model': { 'method': self.fctldash, 'basis': '(auto)', }, 'keywords': { 'level_hint': self.dashlevel, 'params_tweaks': self.dashparams, 'dashcoeff_supplement': self.dashcoeff_supplement, 'apply_qcengine_aliases': True, # for s-dftd3 'verbose': 1, }, 'molecule': molecule.to_schema(dtype=2), 'provenance': p4util.provenance_stamp(__name__), }) jobrec = qcng.compute( resi, self.engine, raise_error=True, task_config={"scratch_directory": core.IOManager.shared_object().get_default_path(), "ncores": core.get_num_threads()}) dashd_part = core.Matrix.from_array(jobrec.extras['qcvars']['DISPERSION CORRECTION GRADIENT']) if wfn is not None: for k, qca in jobrec.extras['qcvars'].items(): if "CURRENT" not in k: wfn.set_variable(k, float(qca) if isinstance(qca, str) else qca) if self.fctldash in ['hf3c', 'pbeh3c']: jobrec = qcng.compute( resi, self.gcp_engine, raise_error=True, task_config={"scratch_directory": core.IOManager.shared_object().get_default_path(), "ncores": core.get_num_threads()}) gcp_part = core.Matrix.from_array(jobrec.return_result) dashd_part.add(gcp_part) return dashd_part else: return self.disp.compute_gradient(molecule)
[docs] def compute_hessian(self, molecule: core.Molecule, wfn: core.Wavefunction = None) -> core.Matrix: """Compute dispersion Hessian based on engine, dispersion level, and parameters in `self`. Uses finite difference, as no dispersion engine has analytic second derivatives. Parameters ---------- molecule System for which to compute empirical dispersion correction. wfn Location to set QCVariables Returns ------- Matrix (3*nat, 3*nat) dispersion Hessian [Eh/a0/a0]. """ optstash = p4util.OptionsState(['PRINT'], ['PARENT_SYMMETRY']) core.set_global_option('PRINT', 0) core.print_out("\n\n Analytical Dispersion Hessians are not supported by any engine.\n") core.print_out(" Computing the Hessian through finite difference of gradients.\n\n") # Setup the molecule molclone = molecule.clone() molclone.reinterpret_coordentry(False) molclone.fix_orientation(True) molclone.fix_com(True) # Record undisplaced symmetry for projection of diplaced point groups core.set_global_option("PARENT_SYMMETRY", molecule.schoenflies_symbol()) findif_meta_dict = driver_findif.hessian_from_gradients_geometries(molclone, -1) for displacement in findif_meta_dict["displacements"].values(): geom_array = np.reshape(displacement["geometry"], (-1, 3)) molclone.set_geometry(core.Matrix.from_array(geom_array)) molclone.update_geometry() displacement["gradient"] = self.compute_gradient(molclone).np.ravel().tolist() H = driver_findif.assemble_hessian_from_gradients(findif_meta_dict, -1) if wfn is not None: wfn.set_variable('DISPERSION CORRECTION HESSIAN', H) optstash.restore() return core.Matrix.from_array(H)