Source code for psi4.driver.driver_nbody

#
# @BEGIN LICENSE
#
# Psi4: an open-source quantum chemistry software package
#
# Copyright (c) 2007-2021 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 itertools
import math
from typing import Callable, Union

import numpy as np

from psi4 import core
from psi4.driver import p4util
from psi4.driver import constants
from psi4.driver.p4util.exceptions import *
from psi4.driver import driver_nbody_helper

### Math helper functions


def nCr(n, r):
    f = math.factorial
    return f(n) // f(r) // f(n - r)


### Begin CBS gufunc data


def _sum_cluster_ptype_data(ptype,
                            ptype_dict,
                            compute_list,
                            fragment_slice_dict,
                            fragment_size_dict,
                            ret,
                            vmfc=False,
                            n=0):
    """
    Sums gradient and hessian data from compute_list.

    compute_list comes in as a tuple(frag, basis)
    """

    if len(compute_list) == 0:
        return

    sign = 1

    # Do ptype
    if ptype == 'gradient':
        for fragn, basisn in compute_list:
            start = 0
            grad = np.asarray(ptype_dict[(fragn, basisn)])

            if vmfc:
                sign = ((-1)**(n - len(fragn)))

            for bas in basisn:
                end = start + fragment_size_dict[bas]
                ret[fragment_slice_dict[bas]] += sign * grad[start:end]
                start += fragment_size_dict[bas]

    elif ptype == 'hessian':
        for fragn, basisn in compute_list:
            hess = np.asarray(ptype_dict[(fragn, basisn)])

            if vmfc:
                sign = ((-1)**(n - len(fragn)))

            # Build up start and end slices
            abs_start, rel_start = 0, 0
            abs_slices, rel_slices = [], []
            for bas in basisn:
                rel_end = rel_start + 3 * fragment_size_dict[bas]
                rel_slices.append(slice(rel_start, rel_end))
                rel_start += 3 * fragment_size_dict[bas]

                tmp_slice = fragment_slice_dict[bas]
                abs_slices.append(slice(tmp_slice.start * 3, tmp_slice.stop * 3))

            for abs_sl1, rel_sl1 in zip(abs_slices, rel_slices):
                for abs_sl2, rel_sl2 in zip(abs_slices, rel_slices):
                    ret[abs_sl1, abs_sl2] += hess[rel_sl1, rel_sl2]

    else:
        raise KeyError("ptype can only be gradient or hessian How did you end up here?")


def _print_nbody_energy(energy_body_dict, header, embedding=False):
    core.print_out("""\n   ==> N-Body: %s energies <==\n\n""" % header)
    core.print_out("""   n-Body     Total Energy [Eh]       I.E. [kcal/mol]      Delta [kcal/mol]\n""")
    previous_e = energy_body_dict[1]
    if previous_e == 0.0:
        tot_e = False
    else:
        tot_e = True
    nbody_range = list(energy_body_dict)
    nbody_range.sort()
    for n in nbody_range:
        delta_e = (energy_body_dict[n] - previous_e)
        delta_e_kcal = delta_e * constants.hartree2kcalmol
        int_e_kcal = (
            energy_body_dict[n] - energy_body_dict[1]) * constants.hartree2kcalmol if not embedding else np.nan
        if tot_e:
            core.print_out("""     %4s  %20.12f  %20.12f  %20.12f\n""" % (n, energy_body_dict[n], int_e_kcal,
                                                                       delta_e_kcal))
        else:
            core.print_out("""     %4s  %20s  %20.12f  %20.12f\n""" % (n, "N/A", int_e_kcal,
                                                                       delta_e_kcal))
        previous_e = energy_body_dict[n]
    core.print_out("\n")


[docs]def nbody_gufunc(func: Union[str, Callable], method_string: str, **kwargs): """ Computes the nbody interaction energy, gradient, or Hessian depending on input. This is a generalized universal function for computing interaction and total quantities. :returns: *return type of func* |w--w| The data. :returns: (*float*, :py:class:`~psi4.core.Wavefunction`) |w--w| data and wavefunction with energy/gradient/hessian set appropriately when **return_wfn** specified. :type func: Callable :param func: ``energy`` || etc. Python function that accepts method_string and a molecule. Returns a energy, gradient, or Hessian as requested. :type method_string: str :param method_string: ``'scf'`` || ``'mp2'`` || ``'ci5'`` || etc. First argument, lowercase and usually unlabeled. Indicates the computational method to be passed to func. :type molecule: :ref:`molecule <op_py_molecule>` :param molecule: ``h2o`` || etc. The target molecule, if not the last molecule defined. :type return_wfn: :ref:`boolean <op_py_boolean>` :param return_wfn: ``'on'`` || |dl| ``'off'`` |dr| Indicate to additionally return the :py:class:`~psi4.core.Wavefunction` calculation result as the second element of a tuple. :type bsse_type: str or list :param bsse_type: ``'cp'`` || ``['nocp', 'vmfc']`` || |dl| ``None`` |dr| || etc. Type of BSSE correction to compute: CP for counterpoise correction, NoCP for plain supramolecular interaction energy, or VMFC for Valiron-Mayer Function Counterpoise correction. If a list is provided, the first string in the list determines which interaction or total energies/gradients/Hessians are returned by this function. By default, this function is not called. :type max_nbody: int :param max_nbody: ``3`` || etc. Maximum n-body to compute, cannot exceed the number of fragments in the molecule. :type ptype: str :param ptype: ``'energy'`` || ``'gradient'`` || ``'hessian'`` Type of the procedure passed in. :type return_total_data: :ref:`boolean <op_py_boolean>` :param return_total_data: ``'on'`` || |dl| ``'off'`` |dr| If True returns the total data (energy/gradient/Hessian) of the system, otherwise returns interaction data. Default is ``'off'`` for energies, ``'on'`` for gradients and Hessians. Note that the calculation of total counterpoise corrected energies implies the calculation of the energies of monomers in the monomer basis, hence specifying ``return_total_data = True`` may carry out more computations than ``return_total_data = False``. :type levels: dict :param levels: ``{1: 'ccsd(t)', 2: 'mp2', 'supersystem': 'scf'}`` || ``{1: 2, 2: 'ccsd(t)', 3: 'mp2'}`` || etc Dictionary of different levels of theory for different levels of expansion Note that method_string is not used in this case. ``supersystem`` computes all higher order n-body effects up to the number of fragments. :type embedding_charges: dict :param embedding_charges: ``{1: [-0.834, 0.417, 0.417], ..}`` Dictionary of atom-centered point charges. keys: 1-based index of fragment, values: list of charges for each fragment. :type charge_method: str :param charge_method: ``scf/6-31g`` || ``b3lyp/6-31g*`` || etc Method to compute point charges for monomers. Overridden by ``embedding_charges`` if both are provided. :type charge_type: str :param charge_type: ``MULLIKEN_CHARGES`` || ``LOWDIN_CHARGES`` Default is ``MULLIKEN_CHARGES`` """ # Initialize dictionaries for easy data passing metadata, component_results, nbody_results = {}, {}, {} # Parse some kwargs kwargs = p4util.kwargs_lower(kwargs) if kwargs.get('levels', False): return driver_nbody_helper.multi_level(func, **kwargs) metadata['ptype'] = kwargs.pop('ptype', None) metadata['return_wfn'] = kwargs.pop('return_wfn', False) metadata['return_total_data'] = kwargs.pop('return_total_data', None) metadata['molecule'] = kwargs.pop('molecule', core.get_active_molecule()) metadata['molecule'].update_geometry() metadata['molecule'].fix_com(True) metadata['molecule'].fix_orientation(True) metadata['embedding_charges'] = kwargs.get('embedding_charges', False) metadata['kwargs'] = kwargs core.clean_variables() if metadata['ptype'] not in ['energy', 'gradient', 'hessian']: raise ValidationError("""N-Body driver: The ptype '%s' is not regonized.""" % metadata['ptype']) if metadata['return_total_data'] is None: if metadata['ptype'] in ['gradient', 'hessian']: metadata['return_total_data'] = True else: metadata['return_total_data'] = False # Parse bsse_type, raise exception if not provided or unrecognized metadata['bsse_type_list'] = kwargs.pop('bsse_type') if metadata['bsse_type_list'] is None: raise ValidationError("N-Body GUFunc: Must pass a bsse_type") if not isinstance(metadata['bsse_type_list'], list): metadata['bsse_type_list'] = [metadata['bsse_type_list']] for num, btype in enumerate(metadata['bsse_type_list']): metadata['bsse_type_list'][num] = btype.lower() if btype.lower() not in ['cp', 'nocp', 'vmfc']: raise ValidationError("N-Body GUFunc: bsse_type '%s' is not recognized" % btype.lower()) metadata['max_nbody'] = kwargs.get('max_nbody', -1) if metadata['molecule'].nfragments() == 1: raise ValidationError("N-Body requires active molecule to have more than 1 fragment.") metadata['max_frag'] = metadata['molecule'].nfragments() if metadata['max_nbody'] == -1: metadata['max_nbody'] = metadata['molecule'].nfragments() else: metadata['max_nbody'] = min(metadata['max_nbody'], metadata['max_frag']) # Flip this off for now, needs more testing # If we are doing CP lets save them integrals #if 'cp' in bsse_type_list and (len(bsse_type_list) == 1): # # Set to save RI integrals for repeated full-basis computations # ri_ints_io = core.get_global_option('DF_INTS_IO') # # inquire if above at all applies to dfmp2 or just scf # core.set_global_option('DF_INTS_IO', 'SAVE') # psioh = core.IOManager.shared_object() # psioh.set_specific_retention(97, True) bsse_str = metadata['bsse_type_list'][0] if len(metadata['bsse_type_list']) > 1: bsse_str = str(metadata['bsse_type_list']) core.print_out("\n\n") core.print_out(" ===> N-Body Interaction Abacus <===\n") core.print_out(" BSSE Treatment: %s\n" % bsse_str) # Get compute list metadata = build_nbody_compute_list(metadata) # Compute N-Body components component_results = compute_nbody_components(func, method_string, metadata) # Assemble N-Body quantities nbody_results = assemble_nbody_components(metadata, component_results) # Build wfn and bind variables wfn = core.Wavefunction.build(metadata['molecule'], 'def2-svp') dicts = [ 'energies', 'ptype', 'intermediates', 'energy_body_dict', 'gradient_body_dict', 'hessian_body_dict', 'nbody', 'cp_energy_body_dict', 'nocp_energy_body_dict', 'vmfc_energy_body_dict' ] if metadata['ptype'] == 'gradient': wfn.set_gradient(nbody_results['ret_ptype']) nbody_results['gradient_body_dict'] = nbody_results['ptype_body_dict'] elif metadata['ptype'] == 'hessian': nbody_results['hessian_body_dict'] = nbody_results['ptype_body_dict'] wfn.set_hessian(nbody_results['ret_ptype']) component_results_gradient = component_results.copy() component_results_gradient['ptype'] = component_results_gradient['gradients'] metadata['ptype'] = 'gradient' nbody_results_gradient = assemble_nbody_components(metadata, component_results_gradient) wfn.set_gradient(nbody_results_gradient['ret_ptype']) nbody_results['gradient_body_dict'] = nbody_results_gradient['ptype_body_dict'] for r in [component_results, nbody_results]: for d in r: if d in dicts: for var, value in r[d].items(): try: wfn.set_scalar_variable(str(var), value) core.set_scalar_variable(str(var), value) except: wfn.set_array_variable(d.split('_')[0].upper() + ' ' + str(var), core.Matrix.from_array(value)) core.set_variable("CURRENT ENERGY", nbody_results['ret_energy']) wfn.set_variable("CURRENT ENERGY", nbody_results['ret_energy']) if metadata['ptype'] == 'gradient': core.set_variable("CURRENT GRADIENT", nbody_results['ret_ptype']) elif metadata['ptype'] == 'hessian': core.set_variable("CURRENT HESSIAN", nbody_results['ret_ptype']) if metadata['return_wfn']: return (nbody_results['ret_ptype'], wfn) else: return nbody_results['ret_ptype']
def build_nbody_compute_list(metadata): """Generates the list of N-Body computations to be performed for a given BSSE type. Parameters ---------- metadata : dict of str Dictionary containing N-body metadata. Required ``'key': value`` pairs: ``'bsse_type_list'``: list of str List of requested BSSE treatments. Possible values include lowercase ``'cp'``, ``'nocp'``, and ``'vmfc'``. ``'max_nbody'``: int Maximum number of bodies to include in the N-Body treatment. Possible: `max_nbody` <= `max_frag` Default: `max_nbody` = `max_frag` ``'max_frag'``: int Number of distinct fragments comprising full molecular supersystem. Returns ------- metadata : dict of str Dictionary containing N-body metadata. New ``'key': value`` pair: ``'compute_dict'`` : dict of str: dict Dictionary containing subdicts enumerating compute lists for each possible BSSE treatment. Contents: ``'all'``: dict of int: set Set containing full list of computations required ``'cp'``: dict of int: set Set containing list of computations required for CP procedure ``'nocp'``: dict of int: set Set containing list of computations required for non-CP procedure ``'vmfc_compute'``: dict of int: set Set containing list of computations required for VMFC procedure ``'vmfc_levels'``: dict of int: set Set containing list of levels required for VMFC procedure """ # What levels do we need? nbody_range = range(1, metadata['max_nbody'] + 1) fragment_range = range(1, metadata['max_frag'] + 1) cp_compute_list = {x: set() for x in nbody_range} nocp_compute_list = {x: set() for x in nbody_range} vmfc_compute_list = {x: set() for x in nbody_range} vmfc_level_list = {x: set() for x in nbody_range} # Need to sum something slightly different # Verify proper passing of bsse_type_list bsse_type_remainder = set(metadata['bsse_type_list']) - {'cp', 'nocp', 'vmfc'} if bsse_type_remainder: raise ValidationError("""Unrecognized BSSE type(s): %s Possible values are 'cp', 'nocp', and 'vmfc'.""" % ', '.join(str(i) for i in bsse_type_remainder)) # Build up compute sets if 'cp' in metadata['bsse_type_list']: # Everything is in dimer basis basis_tuple = tuple(fragment_range) for nbody in nbody_range: for x in itertools.combinations(fragment_range, nbody): if metadata['max_nbody'] == 1: break cp_compute_list[nbody].add((x, basis_tuple)) if 'nocp' in metadata['bsse_type_list'] or metadata['return_total_data']: # Everything in monomer basis for nbody in nbody_range: for x in itertools.combinations(fragment_range, nbody): nocp_compute_list[nbody].add((x, x)) if 'vmfc' in metadata['bsse_type_list']: # Like a CP for all combinations of pairs or greater for nbody in nbody_range: for cp_combos in itertools.combinations(fragment_range, nbody): basis_tuple = tuple(cp_combos) for interior_nbody in nbody_range: for x in itertools.combinations(cp_combos, interior_nbody): combo_tuple = (x, basis_tuple) vmfc_compute_list[interior_nbody].add(combo_tuple) vmfc_level_list[len(basis_tuple)].add(combo_tuple) # Build a comprehensive compute_range compute_list = {x: set() for x in nbody_range} for n in nbody_range: compute_list[n] |= cp_compute_list[n] compute_list[n] |= nocp_compute_list[n] compute_list[n] |= vmfc_compute_list[n] core.print_out(" Number of %d-body computations: %d\n" % (n, len(compute_list[n]))) metadata['compute_dict'] = { 'all': compute_list, 'cp': cp_compute_list, 'nocp': nocp_compute_list, 'vmfc_compute': vmfc_compute_list, 'vmfc_levels': vmfc_level_list } return metadata def compute_nbody_components(func, method_string, metadata): """Computes requested N-body components. Performs requested computations for psi4::Molecule object `molecule` according to `compute_list` with function `func` at `method_string` level of theory. Parameters ---------- func : str {'energy', 'gradient', 'hessian'} Function object to be called within N-Body procedure. method_string : str Indicates level of theory to be passed to function `func`. metadata : dict of str Dictionary of N-body metadata. Required ``'key': value`` pairs: ``'compute_list'``: dict of int: set List of computations to perform. Keys indicate body-levels, e.g,. `compute_list[2]` is the list of all 2-body computations required. ``'kwargs'``: dict Arbitrary keyword arguments to be passed to function `func`. Returns ------- dict of str: dict Dictionary containing computed N-body components. Contents: ``'energies'``: dict of set: float64 Dictionary containing all energy components required for given N-body procedure. ``'ptype'``: dict of set: float64 or dict of set: psi4.Matrix Dictionary of returned quantities from calls of function `func` during N-body computations ``'intermediates'``: dict of str: float64 Dictionary of psivars for intermediate N-body computations to be set at the end of the N-body procedure. """ # Get required metadata kwargs = metadata['kwargs'] molecule = metadata['molecule'] #molecule = core.get_active_molecule() compute_list = metadata['compute_dict']['all'] # Now compute the energies energies_dict = {} gradients_dict = {} ptype_dict = {} intermediates_dict = {} if kwargs.get('charge_method', False) and not metadata['embedding_charges']: metadata['embedding_charges'] = driver_nbody_helper.compute_charges(kwargs['charge_method'], kwargs.get('charge_type', 'MULLIKEN_CHARGES').upper(), molecule) for count, n in enumerate(compute_list.keys()): core.print_out("\n ==> N-Body: Now computing %d-body complexes <==\n\n" % n) total = len(compute_list[n]) for num, pair in enumerate(compute_list[n]): core.print_out( "\n N-Body: Computing complex (%d/%d) with fragments %s in the basis of fragments %s.\n\n" % (num + 1, total, str(pair[0]), str(pair[1]))) ghost = list(set(pair[1]) - set(pair[0])) current_mol = molecule.extract_subsets(list(pair[0]), ghost) current_mol.set_name("%s_%i_%i" % (current_mol.name(), count, num)) if metadata['embedding_charges']: driver_nbody_helper.electrostatic_embedding(metadata, pair=pair) # Save energies info ptype_dict[pair], wfn = func(method_string, molecule=current_mol, return_wfn=True, **kwargs) core.set_global_option_python('EXTERN', None) energies_dict[pair] = core.variable("CURRENT ENERGY") gradients_dict[pair] = wfn.gradient() var_key = "N-BODY (%s)@(%s) TOTAL ENERGY" % (', '.join([str(i) for i in pair[0]]), ', '.join( [str(i) for i in pair[1]])) intermediates_dict[var_key] = core.variable("CURRENT ENERGY") core.print_out("\n N-Body: Complex Energy (fragments = %s, basis = %s: %20.14f)\n" % (str( pair[0]), str(pair[1]), energies_dict[pair])) # Flip this off for now, needs more testing #if 'cp' in bsse_type_list and (len(bsse_type_list) == 1): # core.set_global_option('DF_INTS_IO', 'LOAD') core.clean() return { 'energies': energies_dict, 'gradients': gradients_dict, 'ptype': ptype_dict, 'intermediates': intermediates_dict } def assemble_nbody_components(metadata, component_results): """Assembles N-body components into interaction quantities according to requested BSSE procedure(s). Parameters ----------- metadata : dict of str Dictionary of N-body metadata. Required ``'key': value`` pairs: ``'ptype'``: {'energy', 'gradient', 'hessian'} Procedure which has generated the N-body components to be combined. ``'bsse_type_list'``: list of str List of requested BSSE treatments. Possible values include lowercase ``'cp'``, ``'nocp'``, and ``'vmfc'``. ``'max_nbody'``: int Maximum number of bodies to include in the N-Body treatment. Possible: `max_nbody` <= `max_frag` Default: `max_nbody` = `max_frag` ``'max_frag'``: int Number of distinct fragments comprising full molecular supersystem. ``'energies_dict'``: dict of set: float64 Dictionary containing all energy components required for given N-body procedure. ``'ptype_dict'``: dict of set: float64 or dict of set: psi4.Matrix Dictionary of returned quantities from calls of function `func` during N-body computations ``'compute_dict'``: dict of str: dict Dictionary containing {int: set} subdicts enumerating compute lists for each possible BSSE treatment. ``'kwargs'``: dict Arbitrary keyword arguments. component_results : dict of str: dict Dictionary containing computed N-body components. Required ``'key': value`` pairs: ``'energies'``: dict of set: float64 Dictionary containing all energy components required for given N-body procedure. ``'ptype'``: dict of set: float64 or dict of set: psi4.Matrix Dictionary of returned quantities from calls of function `func` during N-body computations ``'intermediates'``: dict of str: float64 Dictionary of psivars for intermediate N-body computations to be set at the end of the N-body procedure. Returns ------- results : dict of str Dictionary of all N-body results. Contents: ``'ret_energy'``: float64 Interaction data requested. If multiple BSSE types requested in `bsse_type_list`, the interaction data associated with the *first* BSSE type in the list is returned. ``'nbody_dict'``: dict of str: float64 Dictionary of relevant N-body psivars to be set ``'energy_body_dict'``: dict of int: float64 Dictionary of total energies at each N-body level, i.e., ``results['energy_body_dict'][2]`` is the sum of all 2-body total energies for the supersystem. May be empty if ``return_total_data`` is ``False``. ``'ptype_body_dict'``: dict or dict of int: array_like Empty dictionary if `ptype is ``'energy'``, or dictionary of total ptype arrays at each N-body level; i.e., ``results['ptype_body_dict'][2]`` for `ptype` ``'gradient'``is the total 2-body gradient. """ # Unpack metadata kwargs = metadata['kwargs'] nbody_range = range(1, metadata['max_nbody'] + 1) # Unpack compute list metadata compute_list = metadata['compute_dict']['all'] cp_compute_list = metadata['compute_dict']['cp'] nocp_compute_list = metadata['compute_dict']['nocp'] vmfc_compute_list = metadata['compute_dict']['vmfc_compute'] vmfc_level_list = metadata['compute_dict']['vmfc_levels'] # Build size and slices dictionaries fragment_size_dict = { frag: metadata['molecule'].extract_subsets(frag).natom() for frag in range(1, metadata['max_frag'] + 1) } start = 0 fragment_slice_dict = {} for k, v in fragment_size_dict.items(): fragment_slice_dict[k] = slice(start, start + v) start += v molecule_total_atoms = sum(fragment_size_dict.values()) # Final dictionaries cp_energy_by_level = {n: 0.0 for n in nbody_range} nocp_energy_by_level = {n: 0.0 for n in nbody_range} cp_energy_body_dict = {n: 0.0 for n in nbody_range} nocp_energy_body_dict = {n: 0.0 for n in nbody_range} vmfc_energy_body_dict = {n: 0.0 for n in nbody_range} # Build out ptype dictionaries if needed if metadata['ptype'] != 'energy': if metadata['ptype'] == 'gradient': arr_shape = (molecule_total_atoms, 3) elif metadata['ptype'] == 'hessian': arr_shape = (molecule_total_atoms * 3, molecule_total_atoms * 3) else: raise KeyError("N-Body: ptype '%s' not recognized" % ptype) cp_ptype_by_level = {n: np.zeros(arr_shape) for n in nbody_range} nocp_ptype_by_level = {n: np.zeros(arr_shape) for n in nbody_range} vmfc_ptype_by_level = {n: np.zeros(arr_shape) for n in nbody_range} cp_ptype_body_dict = {n: np.zeros(arr_shape) for n in nbody_range} nocp_ptype_body_dict = {n: np.zeros(arr_shape) for n in nbody_range} vmfc_ptype_body_dict = {n: np.zeros(arr_shape) for n in nbody_range} else: cp_ptype_by_level, cp_ptype_body_dict = {}, {} nocp_ptype_by_level, nocp_ptype_body_dict = {}, {} vmfc_ptype_body_dict = {} # Sum up all of the levels nbody_dict = {} for n in nbody_range: # Energy # Extract energies for monomers in monomer basis for CP total data if n == 1: monomers_in_monomer_basis = [v for v in nocp_compute_list[1] if len(v[1]) == 1] monomer_energies = 0.0 monomer_energy_list = [] for i in monomers_in_monomer_basis: monomer_energy_list.append(component_results['energies'][i]) monomer_energies += component_results['energies'][i] cp_energy_by_level[n] = sum(component_results['energies'][v] for v in cp_compute_list[n]) nocp_energy_by_level[n] = sum(component_results['energies'][v] for v in nocp_compute_list[n]) # Special vmfc case if n > 1: vmfc_energy_body_dict[n] = vmfc_energy_body_dict[n - 1] for tup in vmfc_level_list[n]: vmfc_energy_body_dict[n] += ((-1)**(n - len(tup[0]))) * component_results['energies'][tup] # Do ptype if metadata['ptype'] != 'energy': _sum_cluster_ptype_data(metadata['ptype'], component_results['ptype'], cp_compute_list[n], fragment_slice_dict, fragment_size_dict, cp_ptype_by_level[n]) _sum_cluster_ptype_data(metadata['ptype'], component_results['ptype'], nocp_compute_list[n], fragment_slice_dict, fragment_size_dict, nocp_ptype_by_level[n]) _sum_cluster_ptype_data( metadata['ptype'], component_results['ptype'], vmfc_level_list[n], fragment_slice_dict, fragment_size_dict, vmfc_ptype_by_level[n], vmfc=True, n=n) if metadata['ptype'] != 'energy': # Extract ptype data for monomers in monomer basis for CP total data monomer_ptype = np.zeros(arr_shape) _sum_cluster_ptype_data(metadata['ptype'], component_results['ptype'], monomers_in_monomer_basis, fragment_slice_dict, fragment_size_dict, monomer_ptype) # Compute cp energy and ptype if 'cp' in metadata['bsse_type_list']: for n in nbody_range: if n == metadata['max_frag']: cp_energy_body_dict[n] = cp_energy_by_level[n] - bsse if metadata['ptype'] != 'energy': cp_ptype_body_dict[n][:] = cp_ptype_by_level[n] - bsse_ptype continue for k in range(1, n + 1): take_nk = nCr(metadata['max_frag'] - k - 1, n - k) sign = ((-1)**(n - k)) value = cp_energy_by_level[k] cp_energy_body_dict[n] += take_nk * sign * value if metadata['ptype'] != 'energy': value = cp_ptype_by_level[k] cp_ptype_body_dict[n] += take_nk * sign * value if n == 1: bsse = cp_energy_body_dict[n] - monomer_energies cp_energy_body_dict[n] = monomer_energies if metadata['ptype'] != 'energy': bsse_ptype = cp_ptype_body_dict[n] - monomer_ptype cp_ptype_body_dict[n] = monomer_ptype.copy() else: cp_energy_body_dict[n] -= bsse if metadata['ptype'] != 'energy': cp_ptype_body_dict[n] -= bsse_ptype cp_interaction_energy = cp_energy_body_dict[metadata['max_nbody']] - cp_energy_body_dict[1] nbody_dict['Counterpoise Corrected Interaction Energy'] = cp_interaction_energy for n in nbody_range[1:]: var_key = 'CP-CORRECTED %d-BODY INTERACTION ENERGY' % n nbody_dict[var_key] = cp_energy_body_dict[n] - cp_energy_body_dict[1] _print_nbody_energy(cp_energy_body_dict, "Counterpoise Corrected (CP)", metadata['embedding_charges']) cp_interaction_energy = cp_energy_body_dict[metadata['max_nbody']] - cp_energy_body_dict[1] if monomer_energies != 0.0: nbody_dict['Counterpoise Corrected Total Energy'] = cp_energy_body_dict[metadata['max_nbody']] nbody_dict['Counterpoise Corrected Interaction Energy'] = cp_interaction_energy # Compute nocp energy and ptype if 'nocp' in metadata['bsse_type_list']: for n in nbody_range: if n == metadata['max_frag']: nocp_energy_body_dict[n] = nocp_energy_by_level[n] if metadata['ptype'] != 'energy': nocp_ptype_body_dict[n][:] = nocp_ptype_by_level[n] continue for k in range(1, n + 1): take_nk = nCr(metadata['max_frag'] - k - 1, n - k) sign = ((-1)**(n - k)) value = nocp_energy_by_level[k] nocp_energy_body_dict[n] += take_nk * sign * value if metadata['ptype'] != 'energy': value = nocp_ptype_by_level[k] nocp_ptype_body_dict[n] += take_nk * sign * value _print_nbody_energy(nocp_energy_body_dict, "Non-Counterpoise Corrected (NoCP)", metadata['embedding_charges']) nocp_interaction_energy = nocp_energy_body_dict[metadata['max_nbody']] - nocp_energy_body_dict[1] nbody_dict['Non-Counterpoise Corrected Total Energy'] = nocp_energy_body_dict[metadata['max_nbody']] nbody_dict['Non-Counterpoise Corrected Interaction Energy'] = nocp_interaction_energy for n in nbody_range[1:]: var_key = 'NOCP-CORRECTED %d-BODY INTERACTION ENERGY' % n nbody_dict[var_key] = nocp_energy_body_dict[n] - nocp_energy_body_dict[1] # Compute vmfc ptype if 'vmfc' in metadata['bsse_type_list']: if metadata['ptype'] != 'energy': for n in nbody_range: if n > 1: vmfc_ptype_body_dict[n] = vmfc_ptype_by_level[n - 1] vmfc_ptype_body_dict[n] += vmfc_ptype_by_level[n] _print_nbody_energy(vmfc_energy_body_dict, "Valiron-Mayer Function Couterpoise (VMFC)", metadata['embedding_charges']) vmfc_interaction_energy = vmfc_energy_body_dict[metadata['max_nbody']] - vmfc_energy_body_dict[1] nbody_dict['Valiron-Mayer Function Couterpoise Total Energy'] = vmfc_energy_body_dict[metadata['max_nbody']] nbody_dict['Valiron-Mayer Function Couterpoise Interaction Energy'] = vmfc_interaction_energy for n in nbody_range[1:]: var_key = 'VMFC-CORRECTED %d-BODY INTERACTION ENERGY' % n nbody_dict[var_key] = vmfc_energy_body_dict[n] - vmfc_energy_body_dict[1] # Returns results = {} results['nbody'] = nbody_dict for b in ['cp', 'nocp', 'vmfc']: if monomer_energies != 0.0: results['%s_energy_body_dict' % b] = eval('%s_energy_body_dict' % b) results['%s_energy_body_dict' % b] = {str(i) + b: j for i, j in results['%s_energy_body_dict' % b].items()} else: results['%s_energy_body_dict' % b] = {} # Figure out and build return types return_method = metadata['bsse_type_list'][0] if return_method == 'cp': results['ptype_body_dict'] = cp_ptype_body_dict results['energy_body_dict'] = cp_energy_body_dict elif return_method == 'nocp': results['ptype_body_dict'] = nocp_ptype_body_dict results['energy_body_dict'] = nocp_energy_body_dict elif return_method == 'vmfc': results['ptype_body_dict'] = vmfc_ptype_body_dict results['energy_body_dict'] = vmfc_energy_body_dict else: raise ValidationError( "N-Body Wrapper: Invalid return type. Should never be here, please post this error on github.") if metadata['return_total_data']: results['ret_energy'] = results['energy_body_dict'][metadata['max_nbody']] else: results['ret_energy'] = results['energy_body_dict'][metadata['max_nbody']] results['ret_energy'] -= results['energy_body_dict'][1] if metadata['ptype'] != 'energy': if metadata['return_total_data']: np_final_ptype = results['ptype_body_dict'][metadata['max_nbody']].copy() else: np_final_ptype = results['ptype_body_dict'][metadata['max_nbody']].copy() np_final_ptype -= results['ptype_body_dict'][1] results['ret_ptype'] = core.Matrix.from_array(np_final_ptype) else: results['ret_ptype'] = results['ret_energy'] if monomer_energies == 0.0: del results['energy_body_dict'] return results