+#
+# @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
+#
+
+"""Module with functions that encode the sequence of PSI module
+calls for each of the *name* values of the energy(), optimize(),
+response(), and frequency() function. *name* can be assumed lowercase by here.
+
+"""
+import os
+import re
+import shutil
+import subprocess
+import sys
+import warnings
+from typing import Dict, List, Union
+
+import numpy as np
+from qcelemental.util import parse_version, which
+
+from psi4 import core, extras
+
+from .. import p4util
+from .. import psifiles as psif
+from .. import qcdb
+from ..constants import constants
+from ..p4util.exceptions import (
+ ManagedMethodError,
+ PastureRequiredError,
+ UpgradeHelper,
+ ValidationError,
+ docs_table_link,
+)
+
+#from psi4.driver.molutil import *
+from ..qcdb.basislist import corresponding_basis
+from . import dft, empirical_dispersion, mcscf, proc_util, response, solvent
+from .proc_data import method_algorithm_type
+from .roa import run_roa
+
+# never import driver, wrappers, or aliases into this file
+
+
+
+# ADVICE on new additions:
+# * two choices: basic `def run` or managed `def select`
+# * consult http://psicode.org/psi4manual/master/proc_py.html --or-- <psi4-repo>/doc/sphinxman/source/proc_py.rst
+
+
+def select_scf_gradient(name, **kwargs):
+ """Function selecting the algorithm for an SCF gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type("scf") # `"scf"` instead of `name` avoids adding every functional to governing dict in proc_data.py
+ module = core.get_global_option('QC_MODULE')
+
+ if mtd_type == 'CD':
+ # manifestation of `"""No analytic derivatives for SCF_TYPE CD."""`.
+ # here, only hits upon `gradient("scf")` so above message also present in driver.py to catch e.g., mp2 gradient atop a cd reference.
+ func = None
+ else:
+ func = run_scf_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp2(name, **kwargs):
+ """Function selecting the algorithm for a MP2 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ func = run_detci
+ elif module == 'FNOCC':
+ func = run_fnocc
+ elif module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'DFMP2']:
+ func = run_dfmp2
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'DFMP2']:
+ func = run_dfmp2
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ raise UpgradeHelper("energy('mp2')", "energy('zapt2')", 1.7,
+ " Replace method MP with method ZAPT for ROHF reference. DETCI is orders-of-magnitude inefficient for perturbation theory.")
+ elif module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'DFMP2']:
+ func = run_dfmp2
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference in ['RKS', 'UKS']:
+ if mtd_type == 'DF':
+ if module in ['', 'DFMP2']:
+ func = run_dfmp2
+
+ if module == 'DETCI':
+ core.print_out("""\nDETCI is ill-advised for method MP2 as it is available inefficiently as a """
+ """byproduct of a CISD computation.\n""")
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp2_gradient(name, **kwargs):
+ """Function selecting the algorithm for a MP2 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+ all_electron = (core.get_global_option('FREEZE_CORE') == "FALSE")
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc_gradient
+ elif module in ['', 'DFMP2']:
+ func = run_dfmp2_gradient
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module, all_electron])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp2_property(name, **kwargs):
+ """Function selecting the algorithm for a MP2 property call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'DF':
+ #if module == 'OCC':
+ # func = run_dfocc_property
+ if module in ['', 'DFMP2']:
+ func = run_dfmp2_property
+ #elif reference == 'UHF':
+ # if mtd_type == 'DF':
+ # if module in ['', 'OCC']:
+ # func = run_dfocc_property
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2(name, **kwargs):
+ """Function selecting the algorithm for an OMP2 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2_gradient(name, **kwargs):
+ """Function selecting the algorithm for an OMP2 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2_property(name, **kwargs):
+ """Function selecting the algorithm for an OMP2 property call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_property
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2p5_property(name, **kwargs):
+ """Function selecting the algorithm for an OMP2.5 property call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_property
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp3_property(name, **kwargs):
+ """Function selecting the algorithm for an OMP3 property call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_property
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_olccd_property(name, **kwargs):
+ """Function selecting the algorithm for an OLCCD property call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_property
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp3(name, **kwargs):
+ """Function selecting the algorithm for a MP3 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ func = run_detci
+ elif module == 'FNOCC':
+ func = run_fnocc
+ elif module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ raise UpgradeHelper("energy('mp3')", "energy('zapt3')", 1.7,
+ " Replace method MP with method ZAPT for ROHF reference. DETCI is orders-of-magnitude inefficient for perturbation theory.")
+
+ if module == 'DETCI':
+ core.print_out("""\nDETCI is ill-advised for method MP3 as it is available inefficiently as a """
+ """byproduct of a CISD computation.\n""")
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp3_gradient(name, **kwargs):
+ """Function selecting the algorithm for a MP3 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+ all_electron = (core.get_global_option('FREEZE_CORE') == "FALSE")
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module, all_electron])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp3(name, **kwargs):
+ """Function selecting the algorithm for an OMP3 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp3_gradient(name, **kwargs):
+ """Function selecting the algorithm for an OMP3 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp2p5(name, **kwargs):
+ """Function selecting the algorithm for a MP2.5 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp2p5_gradient(name, **kwargs):
+ """Function selecting the algorithm for a MP2.5 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+ all_electron = (core.get_global_option('FREEZE_CORE') == "FALSE")
+
+ func = None
+ if reference in ['RHF', 'UHF']:
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module, all_electron])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2p5(name, **kwargs):
+ """Function selecting the algorithm for an OMP2.5 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_omp2p5_gradient(name, **kwargs):
+ """Function selecting the algorithm for an OMP2.5 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_lccd(name, **kwargs):
+ """Function selecting the algorithm for a LCCD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'OCC':
+ func = run_occ
+ elif module in ['', 'FNOCC']:
+ func = run_cepa
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_lccd_gradient(name, **kwargs):
+ """Function selecting the algorithm for a LCCD gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+ all_electron = (core.get_global_option('FREEZE_CORE') == "FALSE")
+
+ func = None
+ if reference in ['RHF', 'UHF']:
+ if mtd_type == 'CONV':
+ if all_electron:
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module, all_electron])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_olccd(name, **kwargs):
+ """Function selecting the algorithm for an OLCCD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_olccd_gradient(name, **kwargs):
+ """Function selecting the algorithm for an OLCCD gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF', 'UHF', 'ROHF', 'RKS', 'UKS']:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_fnoccsd(name, **kwargs):
+ """Function selecting the algorithm for a FNO-CCSD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'FNOCC']:
+ func = run_fnocc
+ elif mtd_type == 'DF':
+ if module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif mtd_type == 'CD':
+ if module in ['', 'FNOCC']:
+ func = run_fnodfcc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccsd(name, **kwargs):
+ """Function selecting the algorithm for a CCSD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ # [Aug 2022] DF CCSD through CCENERGY for (RHF|ROHF) not enabled here since not advertised. It does run, though, see #2710
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'FNOCC':
+ func = run_fnocc
+ elif module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module == 'CCT3' and extras.addons("cct3"):
+ import cct3
+ func = cct3.run_cct3
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif mtd_type == 'CD':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type in ["DF", "CD"]:
+ if module in ["", "OCC"]:
+ func = run_dfocc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'CCT3' and extras.addons("cct3"):
+ import cct3
+ func = cct3.run_cct3
+ elif module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccsd_gradient(name, **kwargs):
+ """Function selecting the algorithm for a CCSD gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'CCENERGY']:
+ func = run_ccenergy_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'CCENERGY']:
+ func = run_ccenergy_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'CCENERGY']:
+ func = run_ccenergy_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_fnoccsd_t_(name, **kwargs):
+ """Function selecting the algorithm for a FNO-CCSD(T) energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'FNOCC']:
+ func = run_fnocc
+ elif mtd_type == 'DF':
+ if module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif mtd_type == 'CD':
+ if module in ['', 'FNOCC']:
+ func = run_fnodfcc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccsd_t_(name, **kwargs):
+ """Function selecting the algorithm for a CCSD(T) energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'FNOCC':
+ func = run_fnocc
+ elif module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type == 'DF':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif mtd_type == 'CD':
+ if module == 'OCC':
+ func = run_dfocc
+ elif module in ['', 'FNOCC']:
+ func = run_fnodfcc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type in ["DF", "CD"]:
+ if module in ["OCC"]: # SOON "",
+ func = run_dfocc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccsd_t__gradient(name, **kwargs):
+ """Function selecting the algorithm for a CCSD(T) gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ['RHF']:
+ if mtd_type == 'CONV':
+ if module in ['CCENERGY']: # FORMERLY ""
+ func = run_ccenergy_gradient
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc_gradient
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['CCENERGY']: # FORMERLY ""
+ func = run_ccenergy_gradient
+ elif mtd_type == 'DF':
+ if module in ['OCC']: # SOON "",
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccsd_at_(name, **kwargs):
+ """Function selecting the algorithm for a a-CCSD(T) energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ if name.lower() == "a-ccsd(t)":
+ pass
+ elif name.lower() in ["ccsd(at)", "lambda-ccsd(t)", "ccsd(t)_l"]:
+ core.print_out(f"""\nMethod "{name.lower()}" has been regularized to "a-ccsd(t)" for QCVariables.""")
+ name = "a-ccsd(t)"
+
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif reference == "UHF":
+ if mtd_type == 'CONV':
+ if module in ['', 'MRCC'] and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif mtd_type == "DF":
+ if module in ["OCC"]: # SOON "",
+ func = run_dfocc
+ elif mtd_type == "CD":
+ if module in ["OCC"]: # SOON "",
+ func = run_dfocc
+ elif reference == "ROHF":
+ if mtd_type == 'CONV':
+ if module in ['', 'MRCC'] and which("dmrcc", return_bool=True):
+ func = run_mrcc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_cisd(name, **kwargs):
+ """Function selecting the algorithm for a CISD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ func = run_detci
+ elif module in ['', 'FNOCC']:
+ func = run_cepa
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'DETCI']:
+ func = run_detci
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mp4(name, **kwargs):
+ """Function selecting the algorithm for a MP4 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ func = run_detci
+ elif module in ['', 'FNOCC']:
+ func = run_fnocc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'DETCI':
+ raise UpgradeHelper("energy('mp4')", "energy('zapt4')", 1.7,
+ " Replace method MP with method ZAPT for ROHF reference. DETCI is orders-of-magnitude inefficient for perturbation theory.")
+
+ if module == 'DETCI':
+ core.print_out("""\nDETCI is ill-advised for method MP4 as it is available inefficiently as a """
+ """byproduct of a CISDT computation.\n""")
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_remp2(name, **kwargs):
+ """Function selecting the algorithm for a REMP2 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference in ["RHF", "UHF"]:
+ if mtd_type == 'CONV':
+ if module in ['', 'OCC']:
+ func = run_occ
+ elif mtd_type == 'DF':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+ elif mtd_type == 'CD':
+ if module in ['', 'OCC']:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccd(name, **kwargs):
+ """Function selecting the algorithm for a CCD energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option("SCF", "REFERENCE")
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option("QC_MODULE")
+
+ func = None
+ if reference in ["RHF", "UHF"]:
+ if mtd_type == "CONV":
+ if module in [""]:
+ core.print_out("""\nThis method is not available with conventional integrals. Add "set """
+ """cc_type df" or "set cc_type cd" to input to access this method.\n""")
+ elif mtd_type == "DF":
+ if module in ["", "OCC"]:
+ func = run_dfocc
+ elif mtd_type == "CD":
+ if module in ["", "OCC"]:
+ func = run_dfocc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_ccd_gradient(name, **kwargs):
+ """Function selecting the algorithm for a CCD gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option("SCF", "REFERENCE")
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option("QC_MODULE")
+
+ func = None
+ if reference in ["RHF", "UHF"]:
+ if mtd_type == "CONV":
+ if module in [""]:
+ core.print_out("""\nThis method is not available with conventional integrals. Add "set """
+ """cc_type df" or "set cc_type cd" to input to access this method.\n""")
+ elif mtd_type == "DF":
+ if module in ["", "OCC"]:
+ func = run_dfocc_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop("probe", False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_cc2(name, **kwargs):
+ """Function selecting the algorithm for a CC2 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ # [LAB Aug 2022] I'm leaving MRCC CC2 in as a route, but my c.2014 MRCC consistently yields:
+ # "Approximate CC methods are not implemented for excitation level 2!"
+ # [LAB Aug 2022] DF CC2 enabled for test_gradient but only by deliberate `set qc_module ccenergy`
+ # since not advertised. See #2710.
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif mtd_type == 'DF':
+ if module in ['CCENERGY']:
+ func = run_ccenergy
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_cc2_gradient(name, **kwargs):
+ """Function selecting the algorithm for a CC2 gradient call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ # [LAB Aug 2022] Both UHF and ROHF gradients run in ccenergy but ROHF is slightly off (1.e-5)
+ # and UHF is more off (1.e-4). Moreover, manual only claims RHF are working, so restricting here.
+ # [LAB Aug 2022] DF CC2 enabled for test_gradient but only by deliberate `set qc_module ccenergy`
+ # since not advertised. See #2710.
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'CCENERGY']:
+ func = run_ccenergy_gradient
+ elif mtd_type == 'DF':
+ if module in ['CCENERGY']:
+ func = run_ccenergy_gradient
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_cc3(name, **kwargs):
+ """Function selecting the algorithm for a CC3 energy call
+ and directing to specified or best-performance default modules.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module == 'MRCC' and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif module in ['', 'CCENERGY']:
+ func = run_ccenergy
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ # ROHF MRCC CC3 CCn methods are not implemented for ROHF reference!
+ if module in ['', 'CCENERGY']:
+ func = run_ccenergy
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def select_mrcc(name, **kwargs):
+ """Function selecting the algorithm for a CC* energy call
+ and directing to specified MRCC module.
+
+ This function is unusual among "select" functions in that it services multiple methods and a
+ single module. This function could have been skipped and the methods associated directly with
+ run_rmcc; however, routing through this function screens for conv only while
+ providing uniform error messages with other select functions.
+
+ """
+ reference = core.get_option('SCF', 'REFERENCE')
+ type_var, _, mtd_type = method_algorithm_type(name)
+ module = core.get_global_option('QC_MODULE')
+
+ # todo fix table link anchor
+
+ func = None
+ if reference == 'RHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'MRCC'] and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif reference == 'UHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'MRCC'] and which("dmrcc", return_bool=True):
+ func = run_mrcc
+ elif reference == 'ROHF':
+ if mtd_type == 'CONV':
+ if module in ['', 'MRCC'] and which("dmrcc", return_bool=True):
+ func = run_mrcc
+
+ if func is None:
+ raise ManagedMethodError([__name__, name, type_var, mtd_type, reference, module])
+
+ if kwargs.pop('probe', False):
+ return
+ else:
+ return func(name, **kwargs)
+
+
+def build_functional_and_disp(name, restricted, save_pairwise_disp=False, **kwargs):
+
+ if core.has_option_changed("SCF", "DFT_DISPERSION_PARAMETERS"):
+ modified_disp_params = core.get_option("SCF", "DFT_DISPERSION_PARAMETERS")
+ else:
+ modified_disp_params = None
+
+ # Figure out functional
+ superfunc, disp_type = dft.build_superfunctional(name, restricted)
+
+ if disp_type:
+ if isinstance(name, dict):
+ # user dft_functional={} spec - type for lookup, dict val for param defs,
+ # name & citation discarded so only param matches to existing defs will print labels
+ _disp_functor = empirical_dispersion.EmpiricalDispersion(name_hint='',
+ level_hint=disp_type["type"],
+ param_tweaks=disp_type["params"],
+ save_pairwise_disp=save_pairwise_disp,
+ engine=kwargs.get('engine', None))
+ else:
+ # dft/*functionals.py spec - name & type for lookup, option val for param tweaks
+ _disp_functor = empirical_dispersion.EmpiricalDispersion(name_hint=superfunc.name(),
+ level_hint=disp_type["type"],
+ param_tweaks=modified_disp_params,
+ save_pairwise_disp=save_pairwise_disp,
+ engine=kwargs.get('engine', None))
+
+ # [Aug 2018] there once was a breed of `disp_type` that quacked
+ # like a list rather than the more common dict handled above. if
+ # ever again sighted, make an issue so this code can accommodate.
+
+ _disp_functor.print_out()
+ return superfunc, _disp_functor
+
+ else:
+ return superfunc, None
+
+
+
[docs]
+
def scf_wavefunction_factory(name, ref_wfn, reference, **kwargs):
+
"""Builds the correct (R/U/RO/CU HF/KS) wavefunction from the
+
provided information, sets relevant auxiliary basis sets on it,
+
and prepares any empirical dispersion.
+
+
"""
+
# Figure out functional and dispersion
+
superfunc, _disp_functor = build_functional_and_disp(name, restricted=(reference in ["RKS", "RHF"]), **kwargs)
+
+
# Build the wavefunction
+
core.prepare_options_for_module("SCF")
+
if reference in ["RHF", "RKS"]:
+
wfn = core.RHF(ref_wfn, superfunc)
+
elif reference == "ROHF":
+
wfn = core.ROHF(ref_wfn, superfunc)
+
elif reference in ["UHF", "UKS"]:
+
wfn = core.UHF(ref_wfn, superfunc)
+
elif reference == "CUHF":
+
wfn = core.CUHF(ref_wfn, superfunc)
+
else:
+
raise ValidationError("SCF: Unknown reference (%s) when building the Wavefunction." % reference)
+
+
if _disp_functor and _disp_functor.engine != 'nl':
+
wfn._disp_functor = _disp_functor
+
+
# Set the DF basis sets
+
df_needed = core.get_global_option("SCF_TYPE") in ["DF", "MEM_DF", "DISK_DF" ]
+
df_needed |= "DFDIRJ" in core.get_global_option("SCF_TYPE")
+
df_needed |= (core.get_global_option("SCF_TYPE") == "DIRECT" and core.get_option("SCF", "DF_SCF_GUESS"))
+
if df_needed:
+
aux_basis = core.BasisSet.build(wfn.molecule(), "DF_BASIS_SCF",
+
core.get_option("SCF", "DF_BASIS_SCF"),
+
"JKFIT", core.get_global_option('BASIS'),
+
puream=wfn.basisset().has_puream())
+
wfn.set_basisset("DF_BASIS_SCF", aux_basis)
+
else:
+
wfn.set_basisset("DF_BASIS_SCF", core.BasisSet.zero_ao_basis_set())
+
+
# Set the relativistic basis sets
+
if core.get_global_option("RELATIVISTIC") in ["X2C", "DKH"]:
+
decon_basis = core.BasisSet.build(wfn.molecule(), "BASIS_RELATIVISTIC",
+
core.get_option("SCF", "BASIS_RELATIVISTIC"),
+
"DECON", core.get_global_option('BASIS'),
+
puream=wfn.basisset().has_puream())
+
wfn.set_basisset("BASIS_RELATIVISTIC", decon_basis)
+
+
# Set the multitude of SAD basis sets
+
if (core.get_option("SCF", "GUESS") in ["SAD", "SADNO", "HUCKEL", "MODHUCKEL"]):
+
sad_basis_list = core.BasisSet.build(wfn.molecule(), "ORBITAL",
+
core.get_global_option("BASIS"),
+
puream=wfn.basisset().has_puream(),
+
return_atomlist=True)
+
wfn.set_sad_basissets(sad_basis_list)
+
+
if ("DF" in core.get_option("SCF", "SAD_SCF_TYPE")):
+
# We need to force this to spherical regardless of any user or other demands.
+
optstash = p4util.OptionsState(['PUREAM'])
+
core.set_global_option('PUREAM', True)
+
sad_fitting_list = core.BasisSet.build(wfn.molecule(), "DF_BASIS_SAD",
+
core.get_option("SCF", "DF_BASIS_SAD"),
+
puream=True,
+
return_atomlist=True)
+
wfn.set_sad_fitting_basissets(sad_fitting_list)
+
optstash.restore()
+
+
if core.get_option("SCF", "GUESS") == "SAPGAU":
+
# Populate sapgau basis
+
sapgau = core.BasisSet.build(wfn.molecule(), "SAPGAU_BASIS", core.get_global_option("SAPGAU_BASIS"))
+
wfn.set_basisset("SAPGAU", sapgau)
+
+
if hasattr(core, "EXTERN") and 'external_potentials' in kwargs:
+
core.print_out("\n Warning! Both an external potential EXTERN object and the external_potential" +
+
" keyword argument are specified. The external_potentials keyword argument will be ignored.\n")
+
raise ValidationError("double extern")
+
+
ep = kwargs.get("external_potentials", None)
+
if ep is not None:
+
_set_external_potentials_to_wavefunction(ep, wfn)
+
+
return wfn
+
+
+
+def _set_external_potentials_to_wavefunction(external_potential: Union[List, Dict[str, List]], wfn: "core.Wavefunction"):
+ """Initialize :py:class:`psi4.core.ExternalPotential` object(s) from charges and locations and set on **wfn**.
+
+ Parameters
+ ----------
+ external_potential
+ List-like structure where each row corresponds to a charge. Lines can be composed of ``q, [x, y, z]`` or
+ ``q, x, y, z``. Locations are in [a0].
+ Or, dictionary where keys are FI-SAPT fragments A, B, or C and values are as above.
+
+ """
+ from psi4.driver.qmmm import QMMMbohr
+
+ def validate_qxyz(qxyz):
+ if len(qxyz) == 2:
+ return qxyz[0], qxyz[1][0], qxyz[1][1], qxyz[1][2]
+ elif len(qxyz) == 4:
+ return qxyz[0], qxyz[1], qxyz[2], qxyz[3]
+ else:
+ raise ValidationError(f"Point charge '{qxyz}' not mapping into 'chg, [x, y, z]' or 'chg, x, y, z'")
+
+ if isinstance(external_potential, dict):
+ # For FSAPT, we can take a dictionary of external potentials, e.g.,
+ # external_potentials={'A': potA, 'B': potB, 'C': potC} (any optional)
+ # For the dimer SAPT calculation, we need to account for the external potential
+ # in all of the subsystems A, B, C. So we add them all in total_external_potential
+ # and set the external potential to the dimer wave function
+
+ total_external_potential = core.ExternalPotential()
+
+ for frag, frag_qxyz in external_potential.items():
+ if frag.upper() in "ABC":
+ chrgfield = QMMMbohr()
+ for qxyz in frag_qxyz:
+ chrgfield.extern.addCharge(*validate_qxyz(qxyz))
+
+ wfn.set_potential_variable(frag.upper(), chrgfield.extern)
+ total_external_potential.appendCharges(chrgfield.extern.getCharges())
+
+ else:
+ core.print_out("\n Warning! Unknown key for the external_potentials argument: %s" % frag)
+
+ wfn.set_external_potential(total_external_potential)
+
+ else:
+ chrgfield = QMMMbohr()
+ for qxyz in external_potential:
+ chrgfield.extern.addCharge(*validate_qxyz(qxyz))
+ wfn.set_potential_variable("C", chrgfield.extern) # This is for the FSAPT procedure
+ wfn.set_external_potential(chrgfield.extern)
+
+
+
+
[docs]
+
def scf_helper(name, post_scf=True, **kwargs):
+
"""Function serving as helper to SCF, choosing whether to cast
+
up or just run SCF with a standard guess. This preserves
+
previous SCF options set by other procedures (e.g., SAPT
+
output file types for SCF). Most run_* functions should call
+
this function, common exceptions being when multireference
+
SCF is needed or when restarting from converged SCF.
+
+
"""
+
+
if post_scf:
+
name = "scf"
+
+
optstash = p4util.OptionsState(
+
['PUREAM'],
+
['BASIS'],
+
['QMEFP'],
+
['INTS_TOLERANCE'],
+
['DF_BASIS_SCF'],
+
['SCF', 'GUESS'],
+
['SCF', 'DF_INTS_IO'],
+
['SCF', 'ORBITALS_WRITE'],
+
['SCF_TYPE'], # Hack: scope gets changed internally with the Andy trick
+
)
+
+
optstash2 = p4util.OptionsState(
+
['BASIS'],
+
['DF_BASIS_SCF'],
+
['SCF_TYPE'],
+
['SCF', 'DF_INTS_IO'],
+
)
+
+
# Make sure we grab the correctly scoped integral threshold for SCF
+
core.set_global_option('INTS_TOLERANCE', core.get_option('SCF', 'INTS_TOLERANCE'))
+
+
# Grab a few kwargs
+
use_c1 = kwargs.get('use_c1', False)
+
scf_molecule = kwargs.get('molecule', core.get_active_molecule())
+
read_orbitals = core.get_option('SCF', 'GUESS') == "READ"
+
do_timer = kwargs.pop("do_timer", True)
+
ref_wfn = kwargs.pop('ref_wfn', None)
+
if ref_wfn is not None:
+
raise ValidationError("Cannot seed an SCF calculation with a reference wavefunction ('ref_wfn' kwarg).")
+
+
# decide if we keep the checkpoint file
+
_chkfile = kwargs.get('write_orbitals', True)
+
write_checkpoint_file = False
+
if isinstance(_chkfile, str):
+
write_checkpoint_file = True
+
filename = kwargs.get('write_orbitals')
+
core.set_local_option("SCF", "ORBITALS_WRITE", filename)
+
elif _chkfile is True:
+
write_checkpoint_file = True
+
+
# Continuum solvation needs to be run w/o symmetry
+
if core.get_option("SCF", "PCM") or core.get_option("SCF", "DDX"):
+
c1_molecule = scf_molecule.clone()
+
c1_molecule.reset_point_group('c1')
+
c1_molecule.update_geometry()
+
+
scf_molecule = c1_molecule
+
core.print_out(""" PCM or DDX continuum solvation does not make use of molecular symmetry: """
+
"""further calculations in C1 point group.\n""")
+
+
# PE needs to use exactly input orientation to correspond to potfile
+
if core.get_option("SCF", "PE") or "external_potentials" in kwargs:
+
c1_molecule = scf_molecule.clone()
+
if getattr(scf_molecule, "_initial_cartesian", None) is not None:
+
c1_molecule._initial_cartesian = scf_molecule._initial_cartesian.clone()
+
c1_molecule.set_geometry(c1_molecule._initial_cartesian)
+
c1_molecule.reset_point_group("c1")
+
c1_molecule.fix_orientation(True)
+
c1_molecule.fix_com(True)
+
c1_molecule.update_geometry()
+
else:
+
raise ValidationError("Set no_com/no_reorient/symmetry c1 by hand for PE or external potentials on non-Cartesian molecules.")
+
+
scf_molecule = c1_molecule
+
core.print_out(""" PE does not make use of molecular symmetry: """
+
"""further calculations in C1 point group.\n""")
+
core.print_out(""" PE geometry must align with POTFILE keyword: """
+
"""resetting coordinates with fixed origin and orientation.\n""")
+
+
# SCF Banner data
+
banner = kwargs.pop('banner', None)
+
bannername = name
+
+
# Did we pass in a DFT functional?
+
dft_func = kwargs.pop('dft_functional', None)
+
if dft_func is not None:
+
if name.lower() != "scf":
+
raise ValidationError("dft_functional was supplied to SCF, but method name was not SCF ('%s')" % name)
+
name = dft_func
+
bannername = name
+
if isinstance(name, dict):
+
bannername = name.get("name", "custom functional")
+
+
p4util.libint2_print_out()
+
+
# Setup the timer
+
if do_timer:
+
core.tstart()
+
+
# Second-order SCF requires non-symmetric density matrix support
+
if core.get_option('SCF', 'SOSCF'):
+
proc_util.check_non_symmetric_jk_density("Second-order SCF")
+
+
# sort out cast_up settings. no need to stash these since only read, never reset
+
cast = False
+
if core.has_option_changed('SCF', 'BASIS_GUESS'):
+
cast = core.get_option('SCF', 'BASIS_GUESS')
+
if p4util.yes.match(str(cast)):
+
cast = True
+
elif p4util.no.match(str(cast)):
+
cast = False
+
+
if cast:
+
+
# A user can set "BASIS_GUESS" to True and we default to 3-21G
+
if cast is True:
+
guessbasis = corresponding_basis(core.get_global_option('BASIS'), 'GUESS')[0]
+
if guessbasis is None:
+
guessbasis = '3-21G' # guess of last resort
+
else:
+
guessbasis = cast
+
core.set_global_option('BASIS', guessbasis)
+
+
castdf = 'DF' in core.get_global_option('SCF_TYPE')
+
+
if core.has_option_changed('SCF', 'DF_BASIS_GUESS'):
+
castdf = core.get_option('SCF', 'DF_BASIS_GUESS')
+
if p4util.yes.match(str(castdf)):
+
castdf = True
+
elif p4util.no.match(str(castdf)):
+
castdf = False
+
+
if castdf:
+
core.set_global_option('SCF_TYPE', 'DF')
+
core.set_local_option('SCF', 'DF_INTS_IO', 'none')
+
+
# Figure out the fitting basis set
+
if castdf is True:
+
core.set_global_option('DF_BASIS_SCF', '')
+
elif isinstance(castdf, str):
+
core.set_global_option('DF_BASIS_SCF', castdf)
+
else:
+
raise ValidationError("Unexpected castdf option (%s)." % castdf)
+
+
+
# Switch to the guess namespace
+
namespace = core.IO.get_default_namespace()
+
guesspace = namespace + '.guess'
+
if namespace == '':
+
guesspace = 'guess'
+
core.IO.set_default_namespace(guesspace)
+
+
# Print some info about the guess
+
core.print_out('\n')
+
p4util.banner('Guess SCF, %s Basis' % (guessbasis))
+
core.print_out('\n')
+
+
if cast and read_orbitals:
+
raise ValidationError("""Detected options to both cast and read orbitals""")
+
+
if (core.get_option('SCF', 'STABILITY_ANALYSIS') == 'FOLLOW') and (core.get_option('SCF', 'REFERENCE') != 'UHF'):
+
raise ValidationError(f"""Stability analysis root following is only available for unrestricted Hartree--Fock, not present {core.get_option('SCF', 'REFERENCE')}""")
+
+
# If GUESS is auto guess what it should be
+
if core.get_option('SCF', 'GUESS') == "AUTO":
+
if (scf_molecule.natom() > 1):
+
core.set_local_option('SCF', 'GUESS', 'SAD')
+
else:
+
core.set_local_option('SCF', 'GUESS', 'CORE')
+
+
if core.get_global_option('BASIS') in ['', '(AUTO)']:
+
if name in ['hf3c', 'hf-3c']:
+
core.set_global_option('BASIS', 'minix')
+
elif name in ['pbeh3c', 'pbeh-3c']:
+
core.set_global_option('BASIS', 'def2-msvp')
+
elif name in ['r2scan3c', 'r2scan-3c']:
+
core.set_global_option('BASIS', 'def2-mTZVPP')
+
elif name in ['b973c', 'b97-3c']:
+
core.set_global_option('BASIS', 'def2-mTZVP')
+
elif name in ['wb97x3c', 'wb97x-3c']:
+
core.set_global_option('BASIS', 'vdzp')
+
+
# the FIRST scf call
+
if cast:
+
# Cast is a special case
+
base_wfn = core.Wavefunction.build(scf_molecule, core.get_global_option('BASIS'))
+
core.print_out("\n ---------------------------------------------------------\n")
+
if banner:
+
core.print_out(" " + banner.center(58))
+
if cast:
+
core.print_out(" " + "SCF Castup computation".center(58))
+
ref_wfn = scf_wavefunction_factory(name, base_wfn, core.get_option('SCF', 'REFERENCE'), **kwargs)
+
+
# Compute additive correction: dftd3, mp2d, dftd4, etc.
+
if hasattr(ref_wfn, "_disp_functor"):
+
disp_energy = ref_wfn._disp_functor.compute_energy(ref_wfn.molecule())
+
ref_wfn.set_variable("-D Energy", disp_energy)
+
ref_wfn.compute_energy()
+
+
# cast clean-up
+
if cast:
+
+
# Move files to proper namespace
+
core.IO.change_file_namespace(180, guesspace, namespace)
+
core.IO.set_default_namespace(namespace)
+
+
optstash2.restore()
+
+
# Print the banner for the standard operation
+
core.print_out('\n')
+
p4util.banner(bannername.upper())
+
core.print_out('\n')
+
+
# the SECOND scf call
+
base_wfn = core.Wavefunction.build(scf_molecule, core.get_global_option('BASIS'))
+
if banner:
+
core.print_out("\n ---------------------------------------------------------\n")
+
core.print_out(" " + banner.center(58))
+
+
scf_wfn = scf_wavefunction_factory(name, base_wfn, core.get_option('SCF', 'REFERENCE'), **kwargs)
+
+
# The wfn from_file routine adds the npy suffix if needed, but we add it here so that
+
# we can use os.path.isfile to query whether the file exists before attempting to read
+
read_filename = scf_wfn.get_scratch_filename(180) + '.npy'
+
if ((core.get_option('SCF', 'GUESS') == 'READ') and os.path.isfile(read_filename)):
+
old_wfn = core.Wavefunction.from_file(read_filename)
+
+
Ca_occ = old_wfn.Ca_subset("SO", "OCC")
+
Cb_occ = old_wfn.Cb_subset("SO", "OCC")
+
+
if old_wfn.molecule().schoenflies_symbol() != scf_molecule.schoenflies_symbol():
+
raise ValidationError("Cannot compute projection of different symmetries.")
+
+
if old_wfn.basisset().name() == scf_wfn.basisset().name():
+
core.print_out(f" Reading orbitals from file {read_filename}, no projection.\n\n")
+
scf_wfn.guess_Ca(Ca_occ)
+
scf_wfn.guess_Cb(Cb_occ)
+
else:
+
core.print_out(f" Reading orbitals from file {read_filename}, projecting to new basis.\n\n")
+
core.print_out(" Computing basis projection from %s to %s\n\n" % (old_wfn.basisset().name(), scf_wfn.basisset().name()))
+
+
pCa = scf_wfn.basis_projection(Ca_occ, old_wfn.nalphapi(), old_wfn.basisset(), scf_wfn.basisset())
+
pCb = scf_wfn.basis_projection(Cb_occ, old_wfn.nbetapi(), old_wfn.basisset(), scf_wfn.basisset())
+
scf_wfn.guess_Ca(pCa)
+
scf_wfn.guess_Cb(pCb)
+
+
# Strip off headers to only get R, RO, U, CU
+
old_ref = old_wfn.name().replace("KS", "").replace("HF", "")
+
new_ref = scf_wfn.name().replace("KS", "").replace("HF", "")
+
if old_ref != new_ref:
+
scf_wfn.reset_occ_ = True
+
+
elif (core.get_option('SCF', 'GUESS') == 'READ') and not os.path.isfile(read_filename):
+
core.print_out(f"\n !!! Unable to find file {read_filename}, defaulting to SAD guess. !!!\n\n")
+
core.set_local_option('SCF', 'GUESS', 'SAD')
+
sad_basis_list = core.BasisSet.build(scf_wfn.molecule(), "ORBITAL",
+
core.get_global_option("BASIS"),
+
puream=scf_wfn.basisset().has_puream(),
+
return_atomlist=True)
+
scf_wfn.set_sad_basissets(sad_basis_list)
+
+
if ("DF" in core.get_option("SCF", "SAD_SCF_TYPE")):
+
sad_fitting_list = core.BasisSet.build(scf_wfn.molecule(), "DF_BASIS_SAD",
+
core.get_option("SCF", "DF_BASIS_SAD"),
+
puream=scf_wfn.basisset().has_puream(),
+
return_atomlist=True)
+
scf_wfn.set_sad_fitting_basissets(sad_fitting_list)
+
+
+
if cast:
+
core.print_out("\n Computing basis projection from %s to %s\n\n" % (ref_wfn.basisset().name(), base_wfn.basisset().name()))
+
if ref_wfn.basisset().n_ecp_core() != base_wfn.basisset().n_ecp_core():
+
raise ValidationError("Projecting from basis ({}) with ({}) ECP electrons to basis ({}) with ({}) ECP electrons will be a disaster. Select a compatible cast-up basis with `set guess_basis YOUR_BASIS_HERE`.".format(
+
ref_wfn.basisset().name(), ref_wfn.basisset().n_ecp_core(), base_wfn.basisset().name(), base_wfn.basisset().n_ecp_core()))
+
+
pCa = ref_wfn.basis_projection(ref_wfn.Ca(), ref_wfn.nalphapi(), ref_wfn.basisset(), scf_wfn.basisset())
+
pCb = ref_wfn.basis_projection(ref_wfn.Cb(), ref_wfn.nbetapi(), ref_wfn.basisset(), scf_wfn.basisset())
+
scf_wfn.guess_Ca(pCa)
+
scf_wfn.guess_Cb(pCb)
+
+
+
# Print basis set info
+
if core.get_option("SCF", "PRINT_BASIS"):
+
scf_wfn.basisset().print_detail_out()
+
+
# Compute additive correction: dftd3, mp2d, dftd4, etc.
+
if hasattr(scf_wfn, "_disp_functor"):
+
disp_energy = scf_wfn._disp_functor.compute_energy(scf_wfn.molecule(), scf_wfn)
+
scf_wfn.set_variable("-D Energy", disp_energy)
+
+
# PCM preparation
+
if core.get_option('SCF', 'PCM'):
+
if core.get_option('SCF', 'PE'):
+
raise ValidationError("""Error: 3-layer QM/MM/PCM not implemented.\n""")
+
pcmsolver_parsed_fname = core.get_local_option('PCM', 'PCMSOLVER_PARSED_FNAME')
+
pcm_print_level = core.get_option('SCF', "PRINT")
+
scf_wfn.set_PCM(core.PCM(pcmsolver_parsed_fname, pcm_print_level, scf_wfn.basisset()))
+
+
# DDPCM preparation
+
if core.get_option('SCF', 'DDX'):
+
if not solvent._have_ddx:
+
raise ModuleNotFoundError('Python module ddx not found. Solve by installing it: `conda install -c conda-forge pyddx` or `pip install pyddx`')
+
ddx_options = solvent.ddx.get_ddx_options(scf_molecule)
+
scf_wfn.ddx = solvent.ddx.DdxInterface(
+
molecule=scf_molecule, options=ddx_options,
+
basisset=scf_wfn.basisset()
+
)
+
scf_wfn.ddx_state = None
+
+
# PE preparation
+
if core.get_option('SCF', 'PE'):
+
if not solvent._have_pe:
+
raise ModuleNotFoundError('Python module cppe not found. Solve by installing it: `conda install -c conda-forge cppe`')
+
# PE needs information about molecule and basis set
+
pol_embed_options = solvent.pol_embed.get_pe_options()
+
core.print_out(f""" Using potential file
+
{pol_embed_options["potfile"]}
+
for Polarizable Embedding calculation.\n""")
+
scf_wfn.pe_state = solvent.pol_embed.CppeInterface(
+
molecule=scf_molecule, options=pol_embed_options,
+
basisset=scf_wfn.basisset()
+
)
+
+
e_scf = scf_wfn.compute_energy()
+
for obj in [core, scf_wfn]:
+
# set_variable("SCF TOTAL ENERGY") # P::e SCF
+
for pv in ["SCF TOTAL ENERGY", "CURRENT ENERGY", "CURRENT REFERENCE ENERGY"]:
+
obj.set_variable(pv, e_scf)
+
+
# We always would like to print a little property information
+
if kwargs.get('scf_do_properties', True):
+
oeprop = core.OEProp(scf_wfn)
+
oeprop.set_title("SCF")
+
+
# Figure our properties, if empty do dipole
+
props = [x.upper() for x in core.get_option("SCF", "SCF_PROPERTIES")]
+
if "DIPOLE" not in props:
+
props.append("DIPOLE")
+
+
proc_util.oeprop_validator(props)
+
for x in props:
+
oeprop.add(x)
+
+
# Populate free-atom volumes
+
# if we're doing MBIS
+
if 'MBIS_VOLUME_RATIOS' in props:
+
p4util.free_atom_volumes(scf_wfn)
+
+
# Compute properties
+
oeprop.compute()
+
for obj in [core, scf_wfn]:
+
with warnings.catch_warnings():
+
warnings.simplefilter("ignore")
+
obj.set_variable("CURRENT DIPOLE", obj.variable("SCF DIPOLE")) # P::e SCF
+
+
# Write out MO's
+
if core.get_option("SCF", "PRINT_MOS"):
+
mowriter = core.MOWriter(scf_wfn)
+
mowriter.write()
+
+
# Write out a molden file
+
if core.get_option("SCF", "MOLDEN_WRITE"):
+
filename = core.get_writer_file_prefix(scf_molecule.name()) + ".molden"
+
dovirt = bool(core.get_option("SCF", "MOLDEN_WITH_VIRTUAL"))
+
+
occa = scf_wfn.occupation_a()
+
occb = scf_wfn.occupation_a()
+
+
mw = core.MoldenWriter(scf_wfn)
+
mw.write(filename, scf_wfn.Ca(), scf_wfn.Cb(), scf_wfn.epsilon_a(),
+
scf_wfn.epsilon_b(), scf_wfn.occupation_a(),
+
scf_wfn.occupation_b(), dovirt)
+
+
# Write checkpoint file (orbitals and basis); Can be disabled, e.g., for findif displacements
+
if write_checkpoint_file and isinstance(_chkfile, str):
+
filename = kwargs['write_orbitals']
+
scf_wfn.to_file(filename)
+
# core.set_local_option("SCF", "ORBITALS_WRITE", filename)
+
elif write_checkpoint_file:
+
filename = scf_wfn.get_scratch_filename(180)
+
scf_wfn.to_file(filename)
+
extras.register_numpy_file(filename) # retain with -m (messy) option
+
+
if do_timer:
+
core.tstop()
+
+
optstash.restore()
+
+
if (not use_c1) or (scf_molecule.schoenflies_symbol() == 'c1'):
+
return scf_wfn
+
else:
+
# C1 copy quietly
+
c1_optstash = p4util.OptionsState(['PRINT'])
+
core.set_global_option("PRINT", 0)
+
+
# If we force c1 copy the active molecule
+
scf_molecule.update_geometry()
+
core.print_out("""\n A requested method does not make use of molecular symmetry: """
+
"""further calculations in C1 point group.\n\n""")
+
c1_molecule = scf_molecule.clone()
+
c1_molecule.reset_point_group('c1')
+
c1_molecule.fix_orientation(True)
+
c1_molecule.fix_com(True)
+
c1_molecule.update_geometry()
+
c1_basis = core.BasisSet.build(c1_molecule, "ORBITAL", core.get_global_option('BASIS'), quiet=True)
+
tmp = scf_wfn.c1_deep_copy(c1_basis)
+
if not scf_wfn.has_variable("-D ENERGY"):
+
tmp.del_variable("-D ENERGY")
+
c1_jkbasis = core.BasisSet.build(c1_molecule, "DF_BASIS_SCF",
+
core.get_global_option("DF_BASIS_SCF"),
+
"JKFIT", core.get_global_option('BASIS'), quiet=True)
+
tmp.set_basisset("DF_BASIS_SCF", c1_jkbasis)
+
c1_optstash.restore()
+
return tmp
+
+
+
+def run_dct(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density cumulant theory calculation.
+
+ """
+
+ if (core.get_global_option('FREEZE_CORE') == 'TRUE'):
+ raise ValidationError('Frozen core is not available for DCT.')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs)
+
+ if (core.get_global_option("DCT_TYPE") == "DF"):
+ core.print_out(" Constructing Basis Sets for DCT...\n\n")
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_DCT",
+ core.get_global_option("DF_BASIS_DCT"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_DCT", aux_basis)
+
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+ dct_wfn = core.dct(ref_wfn)
+
+ else:
+ # Ensure IWL files have been written for non DF-DCT
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+ dct_wfn = core.dct(ref_wfn)
+
+ for k, v in dct_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return dct_wfn
+
+
+def run_dct_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ DCT gradient calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['GLOBALS', 'DERTYPE'])
+
+
+ core.set_global_option('DERTYPE', 'FIRST')
+ dct_wfn = run_dct_property(name, **kwargs)
+
+ derivobj = core.Deriv(dct_wfn)
+ derivobj.set_tpdm_presorted(True)
+ if core.get_option('DCT', 'DCT_TYPE') == 'CONV':
+ grad = derivobj.compute()
+ else:
+ grad = derivobj.compute_df('DF_BASIS_SCF', 'DF_BASIS_DCT')
+
+ dct_wfn.set_gradient(grad)
+
+ optstash.restore()
+ return dct_wfn
+
+
+def run_dct_property(name, **kwargs):
+ """ Function encoding sequence of PSI module calls for
+ DCT property calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['DCT', 'OPDM'])
+
+ core.set_local_option('DCT', 'OPDM', 'true')
+ dct_wfn = run_dct(name, **kwargs)
+
+ # Run OEProp
+ oe = core.OEProp(dct_wfn)
+ oe.set_title("DCT")
+ for prop in kwargs.get("properties", []):
+ prop = prop.upper()
+ if prop in core.OEProp.valid_methods or "MULTIPOLE(" in prop:
+ oe.add(prop)
+ oe.compute()
+ dct_wfn.oeprop = oe
+
+ for k, v in dct_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return dct_wfn
+
+
+def run_dfocc(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density-fitted or Cholesky-decomposed
+ (non-)orbital-optimized MPN or CC computation.
+
+ """
+ dtl = docs_table_link("dummy", "occ_nonoo")
+
+ optstash = p4util.OptionsState(
+ ['SCF', 'DF_INTS_IO'],
+ ['DFOCC', 'WFN_TYPE'],
+ ['DFOCC', 'ORB_OPT'],
+ ['DFOCC', 'DO_SCS'],
+ ['DFOCC', 'DO_SOS'],
+ ['DFOCC', 'READ_SCF_3INDEX'],
+ ['DFOCC', 'CHOLESKY'],
+ ['DFOCC', 'CC_LAMBDA'])
+
+ def set_cholesky_from(corl_type):
+ if corl_type == 'DF':
+ core.set_local_option('DFOCC', 'CHOLESKY', 'FALSE')
+ proc_util.check_disk_df(name.upper(), optstash)
+
+ elif corl_type == 'CD':
+ core.set_local_option('DFOCC', 'CHOLESKY', 'TRUE')
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ optstash.add_option(['SCF_TYPE'])
+ core.set_global_option('SCF_TYPE', 'CD')
+ core.print_out(""" SCF Algorithm Type (re)set to CD.\n""")
+ if core.get_global_option('SCF_TYPE') != 'CD':
+ core.set_local_option('DFOCC', 'READ_SCF_3INDEX', 'FALSE')
+ else:
+ raise ValidationError(f"""Invalid type '{corl_type}' for DFOCC. See Capabilities Table at {dtl}""")
+
+ if name in ["mp2.5", "mp3"] and not core.has_global_option_changed("MP_TYPE"):
+ core.print_out(f" Information: {name.upper()} default algorithm changed to DF in August 2020. Use `set mp_type conv` for previous behavior.\n")
+
+ director = {
+ "mp2": {"wfn_type": "DF-OMP2", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp2": {"wfn_type": "DF-OMP2", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "mp2.5": {"wfn_type": "DF-OMP2.5", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp2.5": {"wfn_type": "DF-OMP2.5", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "mp3": {"wfn_type": "DF-OMP3", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp3": {"wfn_type": "DF-OMP3", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "remp2": {"wfn_type": "DF-OREMP", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "oremp2": {"wfn_type": "DF-OREMP", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "lccd": {"wfn_type": "DF-OLCCD", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "olccd": {"wfn_type": "DF-OLCCD", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "ccd": {"wfn_type": "DF-CCD", "orb_opt": "FALSE", "nat_orbs": "FALSE",}, # changes to DF-OCCD
+
+ "ccsd": {"wfn_type": "DF-CCSD", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+
+ "ccsd(t)": {"wfn_type": "DF-CCSD(T)", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+
+ "a-ccsd(t)": {"wfn_type": "DF-CCSD(AT)", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for DFOCC energy")
+
+ # throw exception for CONV (approximately). run reference defaulting logic
+ set_cholesky_from(method_algorithm_type(name).now)
+
+ for k, v in director[name].items():
+ core.set_local_option("DFOCC", k.upper(), v)
+
+ core.set_local_option('DFOCC', 'DO_SCS', 'FALSE')
+ core.set_local_option('DFOCC', 'DO_SOS', 'FALSE')
+ core.set_local_option('SCF', 'DF_INTS_IO', 'SAVE')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ else:
+ if ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" DFOCC does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ if not core.get_local_option("DFOCC", "CHOLESKY"):
+ core.print_out(" Constructing Basis Sets for DFOCC...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_CC",
+ core.get_global_option("DF_BASIS_CC"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_CC", aux_basis)
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+
+ dfocc_wfn = core.dfocc(ref_wfn)
+
+ # Shove variables into global space
+ for k, v in dfocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return dfocc_wfn
+
+
+def run_dfocc_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density-fitted (non-)orbital-optimized MPN or CC computation.
+
+ """
+ dtl = docs_table_link("dummy", "occ_nonoo")
+
+ optstash = p4util.OptionsState(
+ ['SCF', 'DF_INTS_IO'],
+ ['REFERENCE'],
+ ['DFOCC', 'WFN_TYPE'],
+ ['DFOCC', 'ORB_OPT'],
+ ['DFOCC', 'CC_LAMBDA'],
+ ['GLOBALS', 'DERTYPE'])
+
+ if name in ["mp2.5", "mp3"] and not core.has_global_option_changed("MP_TYPE"):
+ core.print_out(f" Information: {name.upper()} default algorithm changed to DF in August 2020. Use `set mp_type conv` for previous behavior.\n")
+
+ # CC_LAMBDA keyword was being set TRUE sporadically, but that's covered c-side
+
+ director = {
+ "mp2": {"wfn_type": "DF-OMP2", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp2": {"wfn_type": "DF-OMP2", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "mp2.5": {"wfn_type": "DF-OMP2.5", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp2.5": {"wfn_type": "DF-OMP2.5", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "mp3": {"wfn_type": "DF-OMP3", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "omp3": {"wfn_type": "DF-OMP3", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "remp2": {"wfn_type": "DF-OREMP", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "oremp2": {"wfn_type": "DF-OREMP", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "lccd": {"wfn_type": "DF-OLCCD", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ "olccd": {"wfn_type": "DF-OLCCD", "orb_opt": "TRUE", "nat_orbs": "FALSE",},
+
+ "ccd": {"wfn_type": "DF-CCD", "orb_opt": "FALSE", "nat_orbs": "FALSE",}, # changes to DF-OCCD
+
+ "ccsd": {"wfn_type": "DF-CCSD", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+
+ "ccsd(t)": {"wfn_type": "DF-CCSD(T)", "orb_opt": "FALSE", "nat_orbs": "FALSE",},
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for DFOCC gradient")
+
+ # throw exception for CONV (approximately)
+ if (corl_type := method_algorithm_type(name).now) not in ["DF", "CD"]:
+ raise ValidationError(f"Invalid type {corl_type} for DFOCC gradient. See Capabilities Table at {dtl}")
+
+ proc_util.check_disk_df(name.upper(), optstash)
+
+ # throw exception for SCF_TYPE
+ if core.get_global_option('SCF_TYPE') != 'DISK_DF':
+ raise ValidationError('DFOCC gradients need DF-SCF reference.')
+
+ for k, v in director[name].items():
+ core.set_local_option("DFOCC", k.upper(), v)
+
+ core.set_global_option('DERTYPE', 'FIRST')
+ core.set_local_option('DFOCC', 'DO_SCS', 'FALSE')
+ core.set_local_option('DFOCC', 'DO_SOS', 'FALSE')
+ core.set_local_option('SCF', 'DF_INTS_IO', 'SAVE')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ else:
+ if ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" DFOCC does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ core.print_out(" Constructing Basis Sets for DFOCC...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_CC",
+ core.get_global_option("DF_BASIS_CC"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_CC", aux_basis)
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+ dfocc_wfn = core.dfocc(ref_wfn)
+
+ derivobj = core.Deriv(dfocc_wfn)
+ derivobj.compute_df("DF_BASIS_SCF", "DF_BASIS_CC")
+
+ dfocc_wfn.set_variable(f"{name.upper()} TOTAL GRADIENT", dfocc_wfn.gradient())
+
+ # Shove variables into global space
+ for k, v in dfocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return dfocc_wfn
+
+
+def run_dfocc_property(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density-fitted (non-)orbital-optimized MPN or CC computation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['SCF', 'DF_INTS_IO'],
+ ['DFOCC', 'WFN_TYPE'],
+ ['DFOCC', 'ORB_OPT'],
+ ['DFOCC', 'OEPROP'])
+
+ if name in ['mp2', 'omp2']:
+ core.set_local_option('DFOCC', 'WFN_TYPE', 'DF-OMP2')
+ elif name in ['omp3']:
+ core.set_local_option('DFOCC', 'WFN_TYPE', 'DF-OMP3')
+ elif name in ['omp2.5']:
+ core.set_local_option('DFOCC', 'WFN_TYPE', 'DF-OMP2.5')
+ elif name in ['olccd']:
+ core.set_local_option('DFOCC', 'WFN_TYPE', 'DF-OLCCD')
+ else:
+ raise ValidationError('Unidentified method ' % (name))
+
+ proc_util.check_disk_df(name.upper(), optstash)
+
+ if name in ['mp2']:
+ core.set_local_option('DFOCC', 'ORB_OPT', 'FALSE')
+ elif name in ['omp2', 'omp3', 'omp2.5', 'olccd']:
+ core.set_local_option('DFOCC', 'ORB_OPT', 'TRUE')
+
+ core.set_local_option('DFOCC', 'OEPROP', 'TRUE')
+ core.set_local_option('DFOCC', 'DO_SCS', 'FALSE')
+ core.set_local_option('DFOCC', 'DO_SOS', 'FALSE')
+ core.set_local_option('SCF', 'DF_INTS_IO', 'SAVE')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ else:
+ if ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" DFOCC does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ core.print_out(" Constructing Basis Sets for DFOCC...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_CC",
+ core.get_global_option("DF_BASIS_CC"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_CC", aux_basis)
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+ dfocc_wfn = core.dfocc(ref_wfn)
+
+ # Shove variables into global space
+ # TODO: Make other methods in DFOCC update all variables, then add them to the list. Adding now, risks setting outdated information.
+ if name in ['mp2', 'omp2']:
+ for k, v in dfocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return dfocc_wfn
+
+
+def run_qchf(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an quadratically-convergent SCF computation.
+
+ """
+ dtl = docs_table_link("dummy", "occ_nonoo")
+
+ optstash = p4util.OptionsState(
+ ['SCF', 'DF_INTS_IO'],
+ ['DF_BASIS_SCF'],
+ ['SCF', 'FAIL_ON_MAXITER'],
+ ['MAXITER'],
+ ['DFOCC', 'ORB_OPT'],
+ ['DFOCC', 'WFN_TYPE'],
+ ['DFOCC', 'QCHF'],
+ ['DFOCC', 'E_CONVERGENCE'])
+
+ # throw exception for CONV
+ if (corl_type := method_algorithm_type(name).now) not in ["DISK_DF", "DF", "CD"]:
+ raise ValidationError(f"Invalid type {corl_type} for QCHF energy through `run_qchf`. See Capabilities Table at {dtl}")
+
+ core.set_local_option('DFOCC', 'ORB_OPT', 'TRUE')
+ core.set_local_option('DFOCC', 'WFN_TYPE', 'QCHF')
+ core.set_local_option('DFOCC', 'QCHF', 'TRUE')
+ core.set_local_option('DFOCC', 'E_CONVERGENCE', 8)
+
+ core.set_local_option('SCF', 'DF_INTS_IO', 'SAVE')
+ core.set_local_option('SCF', 'FAIL_ON_MAXITER', False)
+ core.set_local_option('SCF', 'MAXITER', 1)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ else:
+ if ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" QCHF does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+ dfocc_wfn = core.dfocc(ref_wfn)
+
+ # Shove variables into global space
+ for k, v in dfocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return dfocc_wfn
+
+
+def run_occ(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a conventional integral (O)MPN computation
+
+ """
+ # Stash these options so we can reload them at computation end.
+ optstash = p4util.OptionsState(
+ ['OCC', 'SPIN_SCALE_TYPE'],
+ ['OCC', 'ORB_OPT'],
+ ['OCC', 'WFN_TYPE'])
+
+ director = {
+ "mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "NONE", },
+ "scs-mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "SCS", },
+ "scs(n)-mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "SCSN", },
+ "scs-mp2-vdw": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "SCSVDW",},
+ "sos-mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "SOS", },
+ "sos-pi-mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "SOSPI", },
+ "custom-scs-mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE", "spin_scale_type": "CUSTOM",},
+
+ "omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", "spin_scale_type": "NONE", },
+ "scs-omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", "spin_scale_type": "SCS", },
+ "sos-omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", "spin_scale_type": "SOS", },
+ "custom-scs-omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", "spin_scale_type": "CUSTOM",},
+
+ "mp2.5": {"wfn_type": "OMP2.5", "orb_opt": "FALSE", "spin_scale_type": "NONE", },
+ "custom-scs-mp2.5": {"wfn_type": "OMP2.5", "orb_opt": "FALSE", "spin_scale_type": "CUSTOM",},
+
+ "omp2.5": {"wfn_type": "OMP2.5", "orb_opt": "TRUE", "spin_scale_type": "NONE", },
+ "custom-scs-omp2.5": {"wfn_type": "OMP2.5", "orb_opt": "TRUE", "spin_scale_type": "CUSTOM",},
+
+ "mp3": {"wfn_type": "OMP3", "orb_opt": "FALSE", "spin_scale_type": "NONE", },
+ "scs-mp3": {"wfn_type": "OMP3", "orb_opt": "FALSE", "spin_scale_type": "SCS", },
+ "custom-scs-mp3": {"wfn_type": "OMP3", "orb_opt": "FALSE", "spin_scale_type": "CUSTOM",},
+
+ "omp3": {"wfn_type": "OMP3", "orb_opt": "TRUE", "spin_scale_type": "NONE", },
+ "scs-omp3": {"wfn_type": "OMP3", "orb_opt": "TRUE", "spin_scale_type": "SCS", },
+ "sos-omp3": {"wfn_type": "OMP3", "orb_opt": "TRUE", "spin_scale_type": "SOS", },
+ "custom-scs-omp3": {"wfn_type": "OMP3", "orb_opt": "TRUE", "spin_scale_type": "CUSTOM",},
+
+ "remp2": {"wfn_type": "REMP", "orb_opt": "FALSE", "spin_scale_type": "NONE", },
+
+ "oremp2": {"wfn_type": "OREMP", "orb_opt": "TRUE", "spin_scale_type": "NONE", },
+
+ "lccd": {"wfn_type": "OCEPA", "orb_opt": "FALSE", "spin_scale_type": "NONE", },
+ "custom-scs-lccd": {"wfn_type": "OCEPA", "orb_opt": "FALSE", "spin_scale_type": "CUSTOM",},
+
+ "olccd": {"wfn_type": "OCEPA", "orb_opt": "TRUE", "spin_scale_type": "NONE", },
+ "custom-scs-olccd": {"wfn_type": "OCEPA", "orb_opt": "TRUE", "spin_scale_type": "CUSTOM",},
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for OCC energy")
+
+ for k, v in director[name].items():
+ core.set_local_option("OCC", k.upper(), v)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+
+ occ_wfn = core.occ(ref_wfn)
+
+ # Shove variables into global space
+ keep_custom_spin_scaling = core.has_option_changed("OCC", "SS_SCALE") or core.has_option_changed("OCC", "OS_SCALE")
+ for k, v in occ_wfn.variables().items():
+ # Custom spin component scaling variables are meaningless if custom scalings hasn't been set. Delete them.
+ if k.startswith("CUSTOM SCS") and not keep_custom_spin_scaling:
+ occ_wfn.del_variable(k)
+ else:
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return occ_wfn
+
+
+def run_occ_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a conventional integral (O)MPN computation
+ """
+ optstash = p4util.OptionsState(
+ ['OCC', 'ORB_OPT'],
+ ['OCC', 'WFN_TYPE'],
+ ['OCC', 'DO_SCS'],
+ ['OCC', 'DO_SOS'],
+ ['GLOBALS', 'DERTYPE'])
+
+ if core.get_global_option('SCF_TYPE') in ['CD', 'DF', 'MEM_DF', 'DISK_DF']:
+ raise ValidationError('OCC gradients need conventional SCF reference.')
+
+ director = {
+ "mp2": {"wfn_type": "OMP2", "orb_opt": "FALSE",},
+ "omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", },
+ "conv-omp2": {"wfn_type": "OMP2", "orb_opt": "TRUE", },
+
+ "mp2.5": {"wfn_type": "OMP2.5", "orb_opt": "FALSE",},
+ "omp2.5": {"wfn_type": "OMP2.5", "orb_opt": "TRUE", },
+
+ "mp3": {"wfn_type": "OMP3", "orb_opt": "FALSE",},
+ "omp3": {"wfn_type": "OMP3", "orb_opt": "TRUE", },
+
+ "oremp2": {"wfn_type": "OREMP", "orb_opt": "TRUE", },
+
+ "lccd": {"wfn_type": "OCEPA", "orb_opt": "FALSE",},
+ "olccd": {"wfn_type": "OCEPA", "orb_opt": "TRUE", },
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for OCC gradient")
+
+ for k, v in director[name].items():
+ core.set_local_option("OCC", k.upper(), v)
+
+ core.set_global_option('DERTYPE', 'FIRST')
+
+ # locking out SCS through explicit keyword setting
+ # * so that current energy must match call
+ # * since grads not avail for scs
+ core.set_local_option('OCC', 'SPIN_SCALE_TYPE', 'NONE')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ ref_wfn.semicanonicalize()
+
+ occ_wfn = core.occ(ref_wfn)
+
+ derivobj = core.Deriv(occ_wfn)
+ grad = derivobj.compute()
+
+ occ_wfn.set_gradient(grad)
+ occ_wfn.set_variable(f"{name.upper()} TOTAL GRADIENT", grad)
+
+ # Shove variables into global space
+ keep_custom_spin_scaling = core.has_option_changed("OCC", "SS_SCALE") or core.has_option_changed("OCC", "OS_SCALE")
+ for k, v in occ_wfn.variables().items():
+ # Custom spin component scaling variables are meaningless if custom scalings hasn't been set. Delete them.
+ if k.startswith("CUSTOM SCS") and not keep_custom_spin_scaling:
+ occ_wfn.del_variable(k)
+ else:
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return occ_wfn
+
+
+def run_scf(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a self-consistent-field theory (HF & DFT) calculation.
+ """
+ optstash_mp2 = p4util.OptionsState(
+ ['DF_BASIS_MP2'],
+ ['DFMP2', 'MP2_OS_SCALE'],
+ ['DFMP2', 'MP2_SS_SCALE'])
+
+ dft_func = False
+ if "dft_functional" in kwargs:
+ dft_func = True
+
+ optstash_scf = proc_util.scf_set_reference_local(name, is_dft=dft_func)
+
+ # See if we're doing TDSCF after, keep JK if so
+ if sum(core.get_option("SCF", "TDSCF_STATES")) > 0:
+ core.set_local_option("SCF", "SAVE_JK", True)
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ scf_wfn = scf_helper(name, post_scf=False, **kwargs)
+ returnvalue = scf_wfn.energy()
+
+ ssuper = scf_wfn.functional()
+
+ if ssuper.is_c_hybrid():
+
+ # throw exception for CONV/CD MP2
+ if (mp2_type := core.get_global_option("MP2_TYPE")) != "DF":
+ raise ValidationError(f"Invalid MP2 type {mp2_type} for DF-DFT energy. See capabilities Table.")
+
+ core.tstart()
+ aux_basis = core.BasisSet.build(scf_wfn.molecule(), "DF_BASIS_MP2",
+ core.get_option("DFMP2", "DF_BASIS_MP2"),
+ "RIFIT", core.get_global_option('BASIS'),
+ puream=-1)
+ scf_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
+ if ssuper.is_c_scs_hybrid():
+ core.set_local_option('DFMP2', 'MP2_OS_SCALE', ssuper.c_os_alpha())
+ core.set_local_option('DFMP2', 'MP2_SS_SCALE', ssuper.c_ss_alpha())
+ dfmp2_wfn = core.dfmp2(scf_wfn)
+ dfmp2_wfn.compute_energy()
+
+ vdh = dfmp2_wfn.variable('CUSTOM SCS-MP2 CORRELATION ENERGY')
+
+ else:
+ dfmp2_wfn = core.dfmp2(scf_wfn)
+ dfmp2_wfn.compute_energy()
+ vdh = ssuper.c_alpha() * dfmp2_wfn.variable('MP2 CORRELATION ENERGY')
+
+ # remove misleading MP2 psivars computed with DFT, not HF, reference
+ for var in dfmp2_wfn.variables():
+ if var.startswith('MP2 ') and ssuper.name() not in ['MP2D']:
+ scf_wfn.del_variable(var)
+
+ scf_wfn.set_variable("DOUBLE-HYBRID CORRECTION ENERGY", vdh) # P::e SCF
+ scf_wfn.set_variable("{} DOUBLE-HYBRID CORRECTION ENERGY".format(ssuper.name()), vdh)
+ returnvalue += vdh
+ scf_wfn.set_variable("DFT TOTAL ENERGY", returnvalue) # P::e SCF
+ for pv, pvv in scf_wfn.variables().items():
+ if pv.endswith('DISPERSION CORRECTION ENERGY') and pv.startswith(ssuper.name()):
+ fctl_plus_disp_name = pv.split()[0]
+ scf_wfn.set_variable(fctl_plus_disp_name + ' TOTAL ENERGY', returnvalue)
+ break
+ else:
+ scf_wfn.set_variable('{} TOTAL ENERGY'.format(ssuper.name()), returnvalue)
+
+ scf_wfn.set_variable('CURRENT ENERGY', returnvalue)
+ scf_wfn.set_energy(returnvalue)
+ core.print_out('\n\n')
+ core.print_out(' %s Energy Summary\n' % (name.upper()))
+ core.print_out(' ' + '-' * (15 + len(name)) + '\n')
+ core.print_out(' DFT Reference Energy = %22.16lf\n' % (returnvalue - vdh))
+ core.print_out(' Scaled MP2 Correlation = %22.16lf\n' % (vdh))
+ core.print_out(' @Final double-hybrid DFT total energy = %22.16lf\n\n' % (returnvalue))
+ core.tstop()
+
+ if ssuper.name() == 'MP2D':
+ for pv, pvv in dfmp2_wfn.variables().items():
+ scf_wfn.set_variable(pv, pvv)
+
+ # Conversely, remove DFT qcvars from MP2D
+ for var in scf_wfn.variables():
+ if 'DFT ' in var or 'DOUBLE-HYBRID ' in var:
+ scf_wfn.del_variable(var)
+
+ # DFT groups dispersion with SCF. Reshuffle so dispersion with MP2 for MP2D.
+ for pv in ['SCF TOTAL ENERGY', 'SCF ITERATION ENERGY', 'MP2 TOTAL ENERGY']:
+ scf_wfn.set_variable(pv, scf_wfn.variable(pv) - scf_wfn.variable('DISPERSION CORRECTION ENERGY'))
+
+ scf_wfn.set_variable('MP2D CORRELATION ENERGY', scf_wfn.variable('MP2 CORRELATION ENERGY') + scf_wfn.variable('DISPERSION CORRECTION ENERGY'))
+ scf_wfn.set_variable('MP2D TOTAL ENERGY', scf_wfn.variable('MP2D CORRELATION ENERGY') + scf_wfn.variable('HF TOTAL ENERGY'))
+ scf_wfn.set_variable('CURRENT ENERGY', scf_wfn.variable('MP2D TOTAL ENERGY'))
+ scf_wfn.set_variable('CURRENT CORRELATION ENERGY', scf_wfn.variable('MP2D CORRELATION ENERGY'))
+ scf_wfn.set_variable('CURRENT REFERENCE ENERGY', scf_wfn.variable('SCF TOTAL ENERGY'))
+
+ # Shove variables into global space
+ for k, v in scf_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash_scf.restore()
+ optstash_mp2.restore()
+ return scf_wfn
+
+
+def run_scf_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a SCF gradient calculation.
+
+ """
+
+ dft_func = False
+ if "dft_functional" in kwargs:
+ dft_func = True
+
+ optstash = proc_util.scf_set_reference_local(name, is_dft=dft_func)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = run_scf(name, **kwargs)
+
+ if core.get_option('SCF', 'REFERENCE') in ['ROHF', 'CUHF']:
+ ref_wfn.semicanonicalize()
+
+ if hasattr(ref_wfn, "_disp_functor"):
+ disp_grad = ref_wfn._disp_functor.compute_gradient(ref_wfn.molecule(), ref_wfn)
+ ref_wfn.set_variable("-D Gradient", disp_grad)
+
+ grad = core.scfgrad(ref_wfn)
+
+ ref_wfn.set_gradient(grad)
+
+ ref_wfn.set_variable("SCF TOTAL GRADIENT", grad) # P::e SCF
+ if ref_wfn.functional().needs_xc():
+ ref_wfn.set_variable("DFT TOTAL GRADIENT", grad) # overwritten later for DH -- TODO when DH gradients # P::e SCF
+ else:
+ ref_wfn.set_variable("HF TOTAL GRADIENT", grad) # P::e SCF
+
+ # Shove variables into global space
+ for k, v in ref_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return ref_wfn
+
+
+def run_scf_hessian(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an SCF hessian calculation.
+
+ """
+ optstash = proc_util.scf_set_reference_local(name)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = run_scf(name, **kwargs)
+
+ badref = core.get_option('SCF', 'REFERENCE') in ['ROHF', 'CUHF']
+ badint = core.get_global_option('SCF_TYPE') in [ 'CD', 'OUT_OF_CORE']
+ if badref or badint:
+ raise ValidationError("Only RHF/UHF/RKS/UKS Hessians are currently implemented. SCF_TYPE either CD or OUT_OF_CORE not supported")
+
+ if hasattr(ref_wfn, "_disp_functor"):
+ disp_hess = ref_wfn._disp_functor.compute_hessian(ref_wfn.molecule(), ref_wfn)
+ ref_wfn.set_variable("-D Hessian", disp_hess)
+
+ H = core.scfhess(ref_wfn)
+ ref_wfn.set_hessian(H)
+
+ ref_wfn.set_variable("SCF TOTAL HESSIAN", H) # P::e SCF
+ if ref_wfn.functional().needs_xc():
+ ref_wfn.set_variable("DFT TOTAL HESSIAN", H) # overwritten later for DH -- TODO when DH Hessians # P::e SCF
+ else:
+ ref_wfn.set_variable("HF TOTAL HESSIAN", H) # P::e SCF
+
+ # Shove variables into global space
+ for k, v in ref_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return ref_wfn
+
+
+def run_mcscf(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a multiconfigurational self-consistent-field calculation.
+
+ """
+ # Make sure the molecule the user provided is the active one
+ mcscf_molecule = kwargs.get('molecule', core.get_active_molecule())
+ mcscf_molecule.update_geometry()
+ if 'ref_wfn' in kwargs:
+ raise ValidationError("It is not possible to pass run_mcscf a reference wavefunction")
+ new_wfn = core.Wavefunction.build(mcscf_molecule, core.get_global_option('BASIS'))
+
+ return core.mcscf(new_wfn)
+
+
+def run_dfmp2_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a DFMP2 gradient calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['DF_BASIS_SCF'],
+ ['DF_BASIS_MP2'],
+ ['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+ core.print_out(""" SCF Algorithm Type (re)set to DF.\n""")
+
+ if "DF" not in core.get_global_option('SCF_TYPE'):
+ raise ValidationError('DF-MP2 gradients need DF-SCF reference.')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ core.tstart()
+ core.print_out('\n')
+ p4util.banner('DFMP2')
+ core.print_out('\n')
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_MP2",
+ core.get_option("DFMP2", "DF_BASIS_MP2"),
+ "RIFIT", core.get_global_option('BASIS'))
+ ref_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
+
+ dfmp2_wfn = core.dfmp2(ref_wfn)
+ grad = dfmp2_wfn.compute_gradient()
+ dfmp2_wfn.set_gradient(grad)
+
+ # Shove variables into global space
+ dfmp2_wfn.set_variable("MP2 TOTAL GRADIENT", grad) # P::e DFMP2
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('MP2 TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('MP2 CORRELATION ENERGY'))
+ for k, v in dfmp2_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ core.tstop()
+ return dfmp2_wfn
+
+
+def run_dfmp2d_gradient(name, **kwargs):
+ """Encode MP2-D method."""
+
+ dfmp2_wfn = run_dfmp2_gradient('mp2', **kwargs)
+ wfn_grad = dfmp2_wfn.gradient().clone()
+
+ _, _disp_functor = build_functional_and_disp('MP2D', restricted=True)
+ disp_grad = _disp_functor.compute_gradient(dfmp2_wfn.molecule(), dfmp2_wfn)
+ wfn_grad.add(disp_grad)
+ dfmp2_wfn.set_gradient(wfn_grad)
+
+ dfmp2_wfn.set_variable('MP2D CORRELATION ENERGY', dfmp2_wfn.variable('MP2 CORRELATION ENERGY') + dfmp2_wfn.variable('DISPERSION CORRECTION ENERGY'))
+ dfmp2_wfn.set_variable('MP2D TOTAL ENERGY', dfmp2_wfn.variable('MP2D CORRELATION ENERGY') + dfmp2_wfn.variable('HF TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('MP2D TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('MP2D CORRELATION ENERGY'))
+
+ # Shove variables into global space
+ for k, v in dfmp2_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return dfmp2_wfn
+
+
+def run_ccenergy(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a CCSD, CC2, and CC3 calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['TRANSQT2', 'WFN'],
+ ['CCSORT', 'WFN'],
+ ['CCENERGY', 'WFN'])
+
+ if name == 'ccsd':
+ core.set_local_option('TRANSQT2', 'WFN', 'CCSD')
+ core.set_local_option('CCSORT', 'WFN', 'CCSD')
+ core.set_local_option('CCTRANSORT', 'WFN', 'CCSD')
+ core.set_local_option('CCENERGY', 'WFN', 'CCSD')
+ elif name == 'ccsd(t)':
+ core.set_local_option('TRANSQT2', 'WFN', 'CCSD_T')
+ core.set_local_option('CCSORT', 'WFN', 'CCSD_T')
+ core.set_local_option('CCTRANSORT', 'WFN', 'CCSD_T')
+ core.set_local_option('CCENERGY', 'WFN', 'CCSD_T')
+ elif name == 'a-ccsd(t)':
+ core.set_local_option('TRANSQT2', 'WFN', 'CCSD_AT')
+ core.set_local_option('CCSORT', 'WFN', 'CCSD_AT')
+ core.set_local_option('CCTRANSORT', 'WFN', 'CCSD_AT')
+ core.set_local_option('CCENERGY', 'WFN', 'CCSD_AT')
+ core.set_local_option('CCHBAR', 'WFN', 'CCSD_AT')
+ core.set_local_option('CCLAMBDA', 'WFN', 'CCSD_AT')
+ elif name == 'cc2':
+ core.set_local_option('TRANSQT2', 'WFN', 'CC2')
+ core.set_local_option('CCSORT', 'WFN', 'CC2')
+ core.set_local_option('CCTRANSORT', 'WFN', 'CC2')
+ core.set_local_option('CCENERGY', 'WFN', 'CC2')
+ elif name == 'cc3':
+ core.set_local_option('TRANSQT2', 'WFN', 'CC3')
+ core.set_local_option('CCSORT', 'WFN', 'CC3')
+ core.set_local_option('CCTRANSORT', 'WFN', 'CC3')
+ core.set_local_option('CCENERGY', 'WFN', 'CC3')
+ elif name == 'eom-cc2':
+ core.set_local_option('TRANSQT2', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCSORT', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCTRANSORT', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCENERGY', 'WFN', 'EOM_CC2')
+ elif name == 'eom-ccsd':
+ core.set_local_option('TRANSQT2', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCSORT', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCTRANSORT', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCENERGY', 'WFN', 'EOM_CCSD')
+ # Call a plain energy('ccenergy') and have full control over options, incl. wfn
+ elif name == 'ccenergy':
+ pass
+
+ # Bypass routine scf if user did something special to get it to converge
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ if core.get_global_option("CC_TYPE") == "DF":
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_CC",
+ core.get_global_option("DF_BASIS_CC"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_CC", aux_basis)
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ # Obtain semicanonical orbitals
+ if ((core.get_option('SCF', 'REFERENCE') == 'ROHF')
+ and ((name in ['ccsd(t)', 'a-ccsd(t)', 'cc2', 'cc3', 'eom-cc2', 'eom-cc3'])
+ or core.get_option('CCTRANSORT', 'SEMICANONICAL'))):
+ ref_wfn.semicanonicalize()
+
+ if core.get_global_option('RUN_CCTRANSORT'):
+ core.cctransort(ref_wfn)
+ else:
+ try:
+ from psi4.driver.pasture import addins
+ addins.ccsort_transqt2(ref_wfn)
+ except Exception:
+ raise PastureRequiredError("RUN_CCTRANSORT")
+
+
+ ccwfn = core.ccenergy(ref_wfn)
+ if core.get_global_option('PE'):
+ ccwfn.pe_state = ref_wfn.pe_state
+
+ if name == 'a-ccsd(t)':
+ core.cchbar(ref_wfn)
+ lambdawfn = core.cclambda(ref_wfn)
+ for k, v in lambdawfn.variables().items():
+ ccwfn.set_variable(k, v)
+
+ optstash.restore()
+ return ccwfn
+
+
+def run_ccenergy_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a CCSD and CCSD(T) gradient calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['GLOBALS', 'DERTYPE'],
+ ['CCLAMBDA', 'WFN'],
+ ['CCDENSITY', 'WFN'])
+
+ core.set_global_option('DERTYPE', 'FIRST')
+
+ if core.get_global_option('FREEZE_CORE') not in ["FALSE", "0"]:
+ raise ValidationError('Frozen core is not available for the CC gradients.')
+
+ ccwfn = run_ccenergy(name, **kwargs)
+
+ if name == 'cc2':
+ core.set_local_option('CCHBAR', 'WFN', 'CC2')
+ core.set_local_option('CCLAMBDA', 'WFN', 'CC2')
+ core.set_local_option('CCDENSITY', 'WFN', 'CC2')
+ if name == 'ccsd':
+ core.set_local_option('CCLAMBDA', 'WFN', 'CCSD')
+ core.set_local_option('CCDENSITY', 'WFN', 'CCSD')
+ elif name == 'ccsd(t)':
+ core.set_local_option('CCLAMBDA', 'WFN', 'CCSD_T')
+ core.set_local_option('CCDENSITY', 'WFN', 'CCSD_T')
+
+ core.cchbar(ccwfn)
+ core.cclambda(ccwfn)
+ core.ccdensity(ccwfn)
+
+ derivobj = core.Deriv(ccwfn)
+ grad = derivobj.compute()
+ del derivobj
+
+ ccwfn.set_gradient(grad)
+ ccwfn.set_variable(f"{name.upper()} TOTAL GRADIENT", grad)
+ core.set_variable(f"{name.upper()} TOTAL GRADIENT", grad)
+ core.set_variable("CURRENT GRADIENT", grad)
+
+ optstash.restore()
+ return ccwfn
+
+
+def run_bccd(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a Brueckner CCD calculation.
+
+ """
+ dtl = docs_table_link("dummy", "ccenergy")
+
+ optstash = p4util.OptionsState(
+ ['TRANSQT2', 'WFN'],
+ ['CCSORT', 'WFN'],
+ ['CCENERGY', 'WFN'])
+
+ if name == 'bccd':
+ core.set_local_option('TRANSQT2', 'WFN', 'BCCD')
+ core.set_local_option('CCSORT', 'WFN', 'BCCD')
+ core.set_local_option('CCTRANSORT', 'WFN', 'BCCD')
+ core.set_local_option('CCENERGY', 'WFN', 'BCCD')
+
+ elif name == 'bccd(t)':
+ core.set_local_option('TRANSQT2', 'WFN', 'BCCD_T')
+ core.set_local_option('CCSORT', 'WFN', 'BCCD_T')
+ core.set_local_option('CCENERGY', 'WFN', 'BCCD_T')
+ core.set_local_option('CCTRANSORT', 'WFN', 'BCCD_T')
+ core.set_local_option('CCTRIPLES', 'WFN', 'BCCD_T')
+ else:
+ raise ValidationError("proc.py:run_bccd name %s not recognized" % name)
+
+ if (corl_type := method_algorithm_type(name).now) != "CONV":
+ raise ValidationError(f"Invalid type {corl_type} for CCENERGY energy through `run_bccd`. See Capabilities Table at {dtl}")
+
+ # Bypass routine scf if user did something special to get it to converge
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ # Needed for (T).
+ if (core.get_option('SCF', 'REFERENCE') == 'ROHF'):
+ ref_wfn.semicanonicalize()
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ core.set_local_option('CCTRANSORT', 'DELETE_TEI', 'false')
+
+ bcc_iter_cnt = 0
+ if (core.get_global_option("RUN_CCTRANSORT")):
+ sort_func = core.cctransort
+ else:
+ try:
+ from psi4.driver.pasture import addins
+ core.set_local_option('TRANSQT2', 'DELETE_TEI', 'false')
+ sort_func = addins.ccsort_transqt2
+ except Exception:
+ raise PastureRequiredError("RUN_CCTRANSORT")
+
+ hold_qcvars = {
+ "MP2 TOTAL ENERGY": None,
+ "MP2 CORRELATION ENERGY": None,
+ "MP2 SAME-SPIN CORRELATION ENERGY": None,
+ "MP2 OPPOSITE-SPIN CORRELATION ENERGY": None,
+ "MP2 SINGLES ENERGY": None,
+ "MP2 DOUBLES ENERGY": None,
+ "CCSD TOTAL ENERGY": None,
+ "CCSD CORRELATION ENERGY": None,
+ "CCSD SAME-SPIN CORRELATION ENERGY": None,
+ "CCSD OPPOSITE-SPIN CORRELATION ENERGY": None,
+ "CCSD SINGLES ENERGY": None,
+ "CCSD DOUBLES ENERGY": None,
+ }
+
+ while True:
+ sort_func(ref_wfn)
+
+ ref_wfn = core.ccenergy(ref_wfn)
+ core.print_out('Brueckner convergence check: %s\n' % bool(core.variable('BRUECKNER CONVERGED')))
+ if core.variable('BRUECKNER CONVERGED'):
+ break
+
+ if bcc_iter_cnt >= core.get_option('CCENERGY', 'BCCD_MAXITER'):
+ core.print_out("\n\nWarning! BCCD did not converge within the maximum number of iterations.")
+ core.print_out("You can increase the number of BCCD iterations by changing BCCD_MAXITER.\n\n")
+ break
+ bcc_iter_cnt += 1
+
+ if bcc_iter_cnt == 1:
+ for pv in hold_qcvars:
+ hold_qcvars[pv] = ref_wfn.variable(pv)
+
+ ref_wfn.set_variable("BCCD TOTAL ENERGY", ref_wfn.variable("CCSD TOTAL ENERGY"))
+ ref_wfn.set_variable("BCCD CORRELATION ENERGY", ref_wfn.variable("BCCD TOTAL ENERGY") - ref_wfn.variable("SCF TOTAL ENERGY"))
+ ref_wfn.set_variable("CURRENT CORRELATION ENERGY", ref_wfn.variable("BCCD CORRELATION ENERGY"))
+
+ # copy back canonical MP2 and CCSD from initial iteration
+ for pv, v in hold_qcvars.items():
+ if v is not None:
+ ref_wfn.set_variable(pv, v)
+ core.set_variable(pv, v)
+
+ if name == 'bccd(t)':
+ core.cctriples(ref_wfn)
+ ref_wfn.set_variable("B(T) CORRECTION ENERGY", ref_wfn.variable("(T) CORRECTION ENERGY"))
+ ref_wfn.set_variable("BCCD(T) TOTAL ENERGY", ref_wfn.variable("CCSD(T) TOTAL ENERGY"))
+ ref_wfn.set_variable("BCCD(T) CORRELATION ENERGY", ref_wfn.variable("BCCD(T) TOTAL ENERGY") - ref_wfn.variable("SCF TOTAL ENERGY")) # note != CCSD(T) CORRELATION ENERGY
+ ref_wfn.set_variable("CURRENT CORRELATION ENERGY", ref_wfn.variable("BCCD(T) CORRELATION ENERGY"))
+
+ for pv in ["(T) CORRECTION ENERGY", "CCSD(T) TOTAL ENERGY", "CCSD(T) CORRELATION ENERGY"]:
+ ref_wfn.del_variable(pv)
+ core.del_variable(pv)
+
+ for pv in [
+ "BCCD TOTAL ENERGY",
+ "BCCD CORRELATION ENERGY",
+ "B(T) CORRECTION ENERGY",
+ "BCCD(T) TOTAL ENERGY",
+ "BCCD(T) CORRELATION ENERGY",
+ "CURRENT CORRELATION ENERGY",
+ ]:
+ if ref_wfn.has_variable(pv):
+ core.set_variable(pv, ref_wfn.variable(pv))
+
+ # Notes
+ # * BCCD or BCCD(T) correlation energy is total energy of last Brueckner iteration minus HF energy of first Brueckner iteration
+
+ optstash.restore()
+ return ref_wfn
+
+
+def run_tdscf_excitations(wfn,**kwargs):
+
+ states = core.get_option("SCF","TDSCF_STATES")
+
+ # some sanity checks
+ if sum(states) == 0:
+ raise ValidationError("TDSCF: No states requested in TDSCF_STATES")
+
+ # unwrap 1-membered list of states, regardless of symmetry
+ # we will apportion states per irrep later on
+ if len(states) == 1:
+ states = states[0]
+
+ # Tie TDSCF_R_CONVERGENCE to D_CONVERGENCE in SCF reference
+ if core.has_option_changed('SCF', 'TDSCF_R_CONVERGENCE'):
+ r_convergence = core.get_option('SCF', 'TDSCF_R_CONVERGENCE')
+ else:
+ r_convergence = min(1.e-4, core.get_option('SCF', 'D_CONVERGENCE') * 1.e2)
+
+ # "anonymous" return value, as we stash observables in the passed Wavefunction object internally
+ _ = response.scf_response.tdscf_excitations(wfn,
+ states=states,
+ triplets=core.get_option("SCF", "TDSCF_TRIPLETS"),
+ tda=core.get_option("SCF", "TDSCF_TDA"),
+ r_convergence=r_convergence,
+ maxiter=core.get_option("SCF", "TDSCF_MAXITER"),
+ guess=core.get_option("SCF", "TDSCF_GUESS"),
+ verbose=core.get_option("SCF", "TDSCF_PRINT"),
+ coeff_cutoff=core.get_option("SCF", "TDSCF_COEFF_CUTOFF"),
+ tdm_print=core.get_option("SCF", "TDSCF_TDM_PRINT"))
+
+ # Shove variables into global space
+ for k, v in wfn.variables().items():
+ core.set_variable(k, v)
+
+ return wfn
+
+
+def run_tdscf_energy(name, **kwargs):
+
+ # Get a wfn in case we aren't given one
+ ref_wfn = kwargs.get('ref_wfn', None)
+
+ if ref_wfn is None:
+ if name is None:
+ raise ValidationError("TDSCF: No reference wave function!")
+ else:
+ ref_wfn = run_scf(name.strip('td-'), **kwargs)
+
+ return run_tdscf_excitations(ref_wfn, **kwargs)
+
+
+def run_scf_property(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ SCF calculations. This is a simple alias to :py:func:`~proc.run_scf`
+ since SCF properties all handled through oeprop.
+
+ """
+
+ core.tstart()
+ optstash = proc_util.scf_set_reference_local(name)
+
+ properties = kwargs.pop('properties')
+
+ # What response do we need?
+ response_list_vals = list(response.scf_response.property_dicts)
+ oeprop_list_vals = core.OEProp.valid_methods
+
+ oe_properties = []
+ linear_response = []
+ unknown_property = []
+ for prop in properties:
+
+ prop = prop.upper()
+ if prop in response_list_vals:
+ linear_response.append(prop)
+ elif (prop in oeprop_list_vals) or ("MULTIPOLE(" in prop):
+ oe_properties.append(prop)
+ else:
+ unknown_property.append(prop)
+
+ if "DIPOLE" not in oe_properties:
+ oe_properties.append("DIPOLE")
+
+ # Throw if we dont know what something is
+ if len(unknown_property):
+ complete_options = oeprop_list_vals + response_list_vals
+ alt_method_name = p4util.text.find_approximate_string_matches(unknown_property[0],
+ complete_options, 2)
+ alternatives = ""
+ if len(alt_method_name) > 0:
+ alternatives = " Did you mean? %s" % (" ".join(alt_method_name))
+
+ raise ValidationError("SCF Property: Feature '%s' is not recognized. %s" % (unknown_property[0], alternatives))
+
+ # Validate OEProp
+ if len(oe_properties):
+ proc_util.oeprop_validator(oe_properties)
+
+ if len(linear_response):
+ optstash_jk = p4util.OptionsState(["SAVE_JK"])
+ core.set_global_option("SAVE_JK", True)
+
+ # Compute the Wavefunction
+ scf_wfn = run_scf(name, scf_do_properties=False, do_timer=False, **kwargs)
+
+ # Run OEProp
+ oe = core.OEProp(scf_wfn)
+ oe.set_title(name.upper())
+ for prop in oe_properties:
+ oe.add(prop.upper())
+ oe.compute()
+ scf_wfn.oeprop = oe
+
+ core.set_variable("SCF DIPOLE", core.variable(name + " DIPOLE")) # P::e SCF
+
+ # Run Linear Respsonse
+ if len(linear_response):
+ core.prepare_options_for_module("SCF")
+ ret = response.scf_response.cpscf_linear_response(scf_wfn, *linear_response,
+ conv_tol = core.get_global_option("SOLVER_CONVERGENCE"),
+ max_iter = core.get_global_option("SOLVER_MAXITER"),
+ print_lvl = (core.get_global_option("PRINT") + 1))
+ optstash_jk.restore()
+
+ core.tstop()
+ optstash.restore()
+ return scf_wfn
+
+
+def run_cc_property(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ all CC property calculations.
+
+ """
+ optstash = p4util.OptionsState(
+ ['WFN'],
+ ['DERTYPE'],
+ ['ONEPDM'],
+ ['PROPERTY'],
+ ['CCLAMBDA', 'R_CONVERGENCE'],
+ ['CCEOM', 'R_CONVERGENCE'],
+ ['CCEOM', 'E_CONVERGENCE']) # yapf:disable
+
+ oneel_properties = core.OEProp.valid_methods
+ twoel_properties = []
+ response_properties = ['POLARIZABILITY', 'ROTATION', 'ROA', 'ROA_TENSOR']
+ excited_properties = ['OSCILLATOR_STRENGTH', 'ROTATIONAL_STRENGTH']
+
+ one = []
+ two = []
+ response = []
+ excited = []
+ invalid = []
+
+ if 'properties' in kwargs:
+ properties = kwargs['properties']
+
+ for prop in properties:
+
+ prop = prop.upper()
+ if prop in oneel_properties:
+ one.append(prop)
+ elif prop in twoel_properties:
+ two.append(prop)
+ elif prop in response_properties:
+ response.append(prop)
+ elif prop in excited_properties:
+ excited.append(prop)
+ else:
+ invalid.append(prop)
+ else:
+ raise ValidationError("""The "properties" keyword is required with the property() function.""")
+
+ # People are used to requesting dipole/quadrupole and getting dipole,quadrupole,mulliken_charges and NO_occupations
+ if ('DIPOLE' in one) or ('QUADRUPOLE' in one):
+ one = list(set(one + ['DIPOLE', 'QUADRUPOLE', 'MULLIKEN_CHARGES', 'NO_OCCUPATIONS']))
+
+ n_one = len(one)
+ n_two = len(two)
+ n_response = len(response)
+ n_excited = len(excited)
+ n_invalid = len(invalid)
+
+ if n_invalid > 0:
+ print("""The following properties are not currently supported: %s""" % invalid)
+
+ if n_excited > 0 and (name not in ['eom-ccsd', 'eom-cc2']):
+ raise ValidationError("""Excited state CC properties require EOM-CC2 or EOM-CCSD.""")
+
+ if (name in ['eom-ccsd', 'eom-cc2']) and n_response > 0:
+ raise ValidationError("""Cannot (yet) compute response properties for excited states.""")
+
+ if 'roa' in response:
+ # Perform distributed roa job
+ run_roa(name, **kwargs)
+ return # Don't do anything further
+
+ if (n_one > 0 or n_two > 0) and (n_response > 0):
+ print("""Computing both density- and response-based properties.""")
+
+ if n_response > 0:
+ if ("ref_wfn" in kwargs and not kwargs["ref_wfn"].same_a_b_orbs()) or core.get_option('SCF', 'REFERENCE') != 'RHF':
+ raise ValidationError(f"Non-RHF CC response properties are not implemented.")
+
+ if name in ['ccsd', 'cc2', 'eom-ccsd', 'eom-cc2']:
+ this_name = name.upper().replace('-', '_')
+ core.set_global_option('WFN', this_name)
+ ccwfn = run_ccenergy(name, **kwargs)
+ core.set_global_option('WFN', this_name)
+ else:
+ raise ValidationError(f"CC property name {name.upper()} not recognized")
+
+ # Need cchbar for everything
+ core.cchbar(ccwfn)
+
+ # Need ccdensity at this point only for density-based props
+ if n_one > 0 or n_two > 0:
+ if name == 'eom-ccsd':
+ core.set_global_option('WFN', 'EOM_CCSD')
+ core.set_global_option('DERTYPE', 'NONE')
+ core.cceom(ccwfn)
+ elif name == 'eom-cc2':
+ core.set_global_option('WFN', 'EOM_CC2')
+ core.set_global_option('DERTYPE', 'NONE')
+ core.cceom(ccwfn)
+ core.set_global_option('DERTYPE', 'NONE')
+ if core.get_option('CCDENSITY', 'OPDM_RELAX') or n_two > 0:
+ # WARNING!!! A one-particle property computed _with_ a two-particle property will differ
+ # from a one-particle property computed by itself. There are no two-particle properties at
+ # present, so we can kick the issue further down the road.
+ core.set_global_option('OPDM_ONLY', 'FALSE')
+ else:
+ core.set_global_option('OPDM_ONLY', 'TRUE')
+ core.cclambda(ccwfn)
+ core.ccdensity(ccwfn)
+
+ # Need ccresponse only for response-type props
+ if n_response > 0:
+ core.set_global_option('DERTYPE', 'RESPONSE')
+ core.cclambda(ccwfn)
+ for prop in response:
+ core.set_global_option('PROPERTY', prop)
+ core.ccresponse(ccwfn)
+
+ # Excited-state transition properties
+ if n_excited > 0:
+ if name == 'eom-ccsd':
+ core.set_global_option('WFN', 'EOM_CCSD')
+ elif name == 'eom-cc2':
+ core.set_global_option('WFN', 'EOM_CC2')
+ else:
+ raise ValidationError("""Unknown excited-state CC wave function.""")
+ core.set_global_option('DERTYPE', 'NONE')
+ if core.get_option('CCDENSITY', 'OPDM_RELAX'):
+ core.set_global_option('OPDM_ONLY', 'FALSE')
+ else:
+ core.set_global_option('OPDM_ONLY', 'TRUE')
+ # Tight convergence unnecessary for transition properties
+ core.set_local_option('CCLAMBDA', 'R_CONVERGENCE', 1e-4)
+ core.set_local_option('CCEOM', 'R_CONVERGENCE', 1e-4)
+ core.set_local_option('CCEOM', 'E_CONVERGENCE', 1e-5)
+ core.cceom(ccwfn)
+ core.cclambda(ccwfn)
+ core.ccdensity(ccwfn)
+
+ # => Make OEProp calls <=
+ if n_one > 0:
+ # ==> Initialize OEProp <==
+ oe = core.OEProp(ccwfn)
+ for oe_prop_name in one:
+ oe.add(oe_prop_name.upper())
+ # ==> OEProp for the ground state <==
+ # TODO: When Psi is Py 3.9+, transition to the removeprefix version.
+ title = name.upper().replace("EOM-", "")
+ #title = name.upper().removeprefix("EOM-")
+ oe.set_title(title)
+ set_of_names = {title + " {}", "CC {}"}
+ if name.startswith("eom"):
+ gs_h = 0
+ for h, i in enumerate(ccwfn.soccpi()):
+ if i % 2:
+ gs_h = gs_h ^ h
+ ct = ccwfn.molecule().point_group().char_table()
+ total_h_lbl = ct.gamma(0).symbol()
+ gs_h_lbl = ct.gamma(gs_h).symbol()
+ set_of_names.update({title + " ROOT 0 {}", "CC ROOT 0 {}",
+ f"{title} ROOT 0 {{}} - {total_h_lbl} TRANSITION",
+ f"CC ROOT 0 {{}} - {total_h_lbl} TRANSITION",
+ f"{title} ROOT 0 ({gs_h_lbl}) {{}}", f"CC ROOT 0 ({gs_h_lbl}) {{}}",
+ f"{title} ROOT 0 (IN {gs_h_lbl}) {{}}", f"CC ROOT 0 (IN {gs_h_lbl}) {{}}"})
+ oe.set_names(set_of_names)
+ oe.compute()
+
+ # ==> OEProp for Excited States <==
+ if name.startswith('eom'):
+ n_root_pi = core.get_global_option("ROOTS_PER_IRREP")
+ for h in range(ccwfn.nirrep()):
+ root_h_lbl = ct.gamma(h).symbol()
+ trans_h_lbl = ct.gamma(h ^ gs_h).symbol()
+ # Don't forget to count the ground state!
+ for i in range(n_root_pi[h]):
+ if h == gs_h: i += 1
+ root_title = title + f" ROOT {i} (IN {root_h_lbl})"
+ oe.set_title(root_title)
+ total_idx = ccwfn.total_index(i, h)
+ set_of_names = {f"{title} ROOT {total_idx} {{}}", f"CC ROOT {total_idx} {{}}",
+ f"{title} ROOT {total_idx} {{}} - {trans_h_lbl} TRANSITION",
+ f"CC ROOT {total_idx} {{}} - {trans_h_lbl} TRANSITION",
+ f"{title} ROOT {total_idx} ({root_h_lbl}) {{}}", f"CC ROOT {total_idx} ({root_h_lbl}) {{}}",
+ f"{title} ROOT {i} (IN {root_h_lbl}) {{}}", f"CC ROOT {i} (IN {root_h_lbl}) {{}}"}
+ oe.set_names(set_of_names)
+ Da = ccwfn.get_density(root_title + " ALPHA")
+ oe.set_Da_so(Da)
+ if not ccwfn.same_a_b_dens():
+ Db = ccwfn.get_density(root_title + " BETA")
+ oe.set_Db_so(Db)
+ oe.compute()
+
+ core.set_global_option('WFN', 'SCF')
+ core.revoke_global_option_changed('WFN')
+ core.set_global_option('DERTYPE', 'NONE')
+ core.revoke_global_option_changed('DERTYPE')
+
+ optstash.restore()
+ return ccwfn
+
+
+def run_dfmp2_property(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a DFMP2 property calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['DF_BASIS_SCF'],
+ ['DF_BASIS_MP2'],
+ ['ONEPDM'],
+ ['OPDM_RELAX'],
+ ['SCF_TYPE'])
+
+ core.set_global_option('ONEPDM', 'TRUE')
+ core.set_global_option('OPDM_RELAX', 'TRUE')
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF') # local set insufficient b/c SCF option read in DFMP2
+ core.print_out(""" SCF Algorithm Type (re)set to DF.\n""")
+
+ if 'DF' not in core.get_global_option('SCF_TYPE'):
+ raise ValidationError('DF-MP2 properties need DF-SCF reference.')
+
+ properties = kwargs.pop('properties')
+ proc_util.oeprop_validator(properties)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, scf_do_properties=False, use_c1=True, **kwargs) # C1 certified
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_MP2",
+ core.get_option("DFMP2", "DF_BASIS_MP2"),
+ "RIFIT", core.get_global_option('BASIS'))
+ ref_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
+
+ core.tstart()
+ core.print_out('\n')
+ p4util.banner('DFMP2')
+ core.print_out('\n')
+
+ dfmp2_wfn = core.dfmp2(ref_wfn)
+ grad = dfmp2_wfn.compute_gradient()
+
+ if name == 'scs-mp2':
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('SCS-MP2 TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('SCS-MP2 CORRELATION ENERGY'))
+ elif name == 'mp2':
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('MP2 TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('MP2 CORRELATION ENERGY'))
+
+ # Run OEProp
+ oe = core.OEProp(dfmp2_wfn)
+ oe.set_title(name.upper())
+ for prop in properties:
+ oe.add(prop.upper())
+ oe.compute()
+ dfmp2_wfn.oeprop = oe
+
+ # Shove variables into global space
+ for k, v in dfmp2_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ core.tstop()
+ return dfmp2_wfn
+
+
+def _clean_detci(keep: bool=True):
+ psioh = core.IOManager.shared_object()
+ psio = core.IO.shared_object()
+ cifl = core.get_option("DETCI", "CI_FILE_START")
+ for fl in range(cifl, cifl + 4):
+ if psio.open_check(fl):
+ psio.close(fl, keep)
+
+
+def run_detci_property(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a configuration interaction calculation, namely FCI,
+ CIn, MPn, and ZAPTn, computing properties.
+
+ """
+ optstash = p4util.OptionsState(
+ ['OPDM'],
+ ['TDM'])
+
+ # Find valid properties
+ valid_transition = ['TRANSITION_DIPOLE', 'TRANSITION_QUADRUPOLE']
+
+ ci_prop = []
+ ci_trans = []
+ properties = kwargs.pop('properties')
+ for prop in properties:
+ if prop.upper() in valid_transition:
+ ci_trans.append(prop)
+ else:
+ ci_prop.append(prop)
+
+ proc_util.oeprop_validator(ci_prop)
+
+ core.set_global_option('OPDM', 'TRUE')
+ if len(ci_trans):
+ core.set_global_option('TDM', 'TRUE')
+
+ # Compute
+ if name in ['mcscf', 'rasscf', 'casscf']:
+ ciwfn = run_detcas(name, **kwargs)
+ else:
+ ciwfn = run_detci(name, **kwargs)
+
+ # All property names are just CI
+ if 'CI' in name.upper():
+ name = 'CI'
+
+ states = core.get_global_option('avg_states')
+ nroots = core.get_global_option('num_roots')
+ if len(states) != nroots:
+ states = range(nroots)
+
+ # Run OEProp
+ oe = core.OEProp(ciwfn)
+ oe.set_title(name.upper())
+ for prop in ci_prop:
+ oe.add(prop.upper())
+
+ # Compute "the" CI density
+ oe.compute()
+ ciwfn.oeprop = oe
+
+ # If we have more than one root, compute all data
+ if nroots > 1:
+ core.print_out("\n ===> %s properties for all CI roots <=== \n\n" % name.upper())
+ for root in states:
+ oe.set_title("%s ROOT %d" % (name.upper(), root))
+ if ciwfn.same_a_b_dens():
+ oe.set_Da_mo(ciwfn.get_opdm(root, root, "A", True))
+ else:
+ oe.set_Da_mo(ciwfn.get_opdm(root, root, "A", True))
+ oe.set_Db_mo(ciwfn.get_opdm(root, root, "B", True))
+ oe.compute()
+
+ # Transition density matrices
+ if (nroots > 1) and len(ci_trans):
+ oe.clear()
+ for tprop in ci_trans:
+ oe.add(tprop.upper())
+
+ core.print_out("\n ===> %s properties for all CI transition density matrices <=== \n\n" % name.upper())
+ for root in states[1:]:
+ oe.set_title("%s ROOT %d -> ROOT %d" % (name.upper(), 0, root))
+ if ciwfn.same_a_b_dens():
+ oe.set_Da_mo(ciwfn.get_opdm(0, root, "A", True))
+ else:
+ oe.set_Da_mo(ciwfn.get_opdm(0, root, "A", True))
+ oe.set_Db_mo(ciwfn.get_opdm(0, root, "B", True))
+ oe.compute()
+
+ _clean_detci()
+ optstash.restore()
+ return ciwfn
+
+
+def run_eom_cc(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an EOM-CC calculation, namely EOM-CC2, EOM-CCSD, and EOM-CC3.
+
+ """
+ optstash = p4util.OptionsState(
+ ['TRANSQT2', 'WFN'],
+ ['CCSORT', 'WFN'],
+ ['CCENERGY', 'WFN'],
+ ['CCHBAR', 'WFN'],
+ ['CCEOM', 'WFN'])
+
+ if name == 'eom-ccsd':
+ core.set_local_option('TRANSQT2', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCSORT', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCENERGY', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCHBAR', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCEOM', 'WFN', 'EOM_CCSD')
+ ref_wfn = run_ccenergy('ccsd', **kwargs)
+
+ elif name == 'eom-cc2':
+ user_ref = core.get_option('CCENERGY', 'REFERENCE')
+ if (user_ref != 'RHF') and (user_ref != 'UHF'):
+ raise ValidationError('Reference %s for EOM-CC2 is not available.' % user_ref)
+ core.set_local_option('TRANSQT2', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCSORT', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCENERGY', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCHBAR', 'WFN', 'EOM_CC2')
+ core.set_local_option('CCEOM', 'WFN', 'EOM_CC2')
+ ref_wfn = run_ccenergy('cc2', **kwargs)
+
+ elif name == 'eom-cc3':
+ core.set_local_option('TRANSQT2', 'WFN', 'EOM_CC3')
+ core.set_local_option('CCSORT', 'WFN', 'EOM_CC3')
+ core.set_local_option('CCENERGY', 'WFN', 'EOM_CC3')
+ core.set_local_option('CCHBAR', 'WFN', 'EOM_CC3')
+ core.set_local_option('CCEOM', 'WFN', 'EOM_CC3')
+ ref_wfn = run_ccenergy('cc3', **kwargs)
+
+ core.cchbar(ref_wfn)
+ core.cceom(ref_wfn)
+
+ optstash.restore()
+ return ref_wfn
+ # TODO ask if all these cc modules not actually changing wfn
+
+
+def run_eom_cc_gradient(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an EOM-CCSD gradient calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['CCDENSITY', 'XI'],
+ ['CCDENSITY', 'ZETA'],
+ ['CCLAMBDA', 'ZETA'],
+ ['DERTYPE'],
+ ['CCDENSITY', 'WFN'],
+ ['CCLAMBDA', 'WFN'])
+
+ core.set_global_option('DERTYPE', 'FIRST')
+
+ if name == 'eom-ccsd':
+ core.set_local_option('CCLAMBDA', 'WFN', 'EOM_CCSD')
+ core.set_local_option('CCDENSITY', 'WFN', 'EOM_CCSD')
+ ref_wfn = run_eom_cc(name, **kwargs)
+ else:
+ core.print_out('DGAS: proc.py:1599 hitting an undefined sequence')
+ core.clean()
+ raise ValueError('Hit a wall in proc.py:1599')
+
+ core.set_local_option('CCLAMBDA', 'ZETA', 'FALSE')
+ core.set_local_option('CCDENSITY', 'ZETA', 'FALSE')
+ core.set_local_option('CCDENSITY', 'XI', 'TRUE')
+ core.cclambda(ref_wfn)
+ core.ccdensity(ref_wfn)
+ core.set_local_option('CCLAMBDA', 'ZETA', 'TRUE')
+ core.set_local_option('CCDENSITY', 'ZETA', 'TRUE')
+ core.set_local_option('CCDENSITY', 'XI', 'FALSE')
+ core.cclambda(ref_wfn)
+ core.ccdensity(ref_wfn)
+
+ derivobj = core.Deriv(ref_wfn)
+ grad = derivobj.compute()
+
+ ref_wfn.set_gradient(grad)
+
+ optstash.restore()
+ return ref_wfn
+
+
+def run_adcc(name, **kwargs):
+ """Prepare and run an ADC calculation in adcc, interpret the result and return
+ as a wavefunction.
+
+ """
+ # TODO Maybe it would improve readability if this function was spilt
+ # up and the whole thing went to a separate file (like for sapt,
+ # interface_cfour.py, ...
+
+ try:
+ import adcc
+ from adcc.exceptions import InvalidReference
+ except ModuleNotFoundError:
+ raise ValidationError("adcc extras qc_module not available. Try installing "
+ "via 'pip install adcc' or 'conda install -c conda-forge adcc'.")
+
+ min_version = "0.15.16"
+ if parse_version(adcc.__version__) < parse_version(min_version):
+ raise ModuleNotFoundError("adcc version {} is required at least. "
+ "Version {}"
+ " was found.".format(min_version,
+ adcc.__version__))
+
+ if core.get_option('ADC', 'REFERENCE') not in ["RHF", "UHF"]:
+ raise ValidationError('adcc requires reference RHF or UHF')
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.pop('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs)
+
+ # Start timer
+ do_timer = kwargs.pop("do_timer", True)
+ if do_timer:
+ core.tstart()
+
+ #
+ # Build kwargs for adcc
+ #
+ kwargs.pop("molecule", None)
+
+ if ref_wfn.frzcpi()[0] > 0:
+ kwargs["frozen_core"] = ref_wfn.frzcpi()[0]
+ if ref_wfn.frzvpi()[0] > 0:
+ kwargs["frozen_virtual"] = ref_wfn.frzvpi()[0]
+ if core.get_option("ADC", "NUM_CORE_ORBITALS"):
+ kwargs["core_orbitals"] = core.get_option("ADC", "NUM_CORE_ORBITALS")
+
+ scf_accuracy = max(core.get_option("SCF", "E_CONVERGENCE"),
+ core.get_option("SCF", "D_CONVERGENCE"))
+ if core.get_option("ADC", "R_CONVERGENCE") < 0:
+ kwargs["conv_tol"] = max(100 * scf_accuracy, 1e-6)
+ else:
+ kwargs["conv_tol"] = core.get_option("ADC", "R_CONVERGENCE")
+
+ n_roots = core.get_option('ADC', 'ROOTS_PER_IRREP')
+ if len(n_roots) > 1:
+ raise ValidationError("adcc can only deal with a single irrep.")
+ kwargs["n_states"] = n_roots[0]
+
+ if core.get_option("ADC", "NUM_GUESSES") > 0:
+ kwargs["n_guesses"] = core.get_option("ADC", "NUM_GUESSES")
+ if core.get_option("ADC", "MAX_NUM_VECS") > 0:
+ kwargs["max_subspace"] = core.get_option("ADC", "MAX_NUM_VECS")
+
+ kind = core.get_option("ADC", "KIND").lower()
+ if isinstance(ref_wfn, core.UHF):
+ if not core.has_option_changed("ADC", "KIND"):
+ kind = "any"
+ elif kind not in ["any", "spin_flip"]:
+ raise ValidationError("For UHF references the only valid values for 'KIND' are "
+ "'SPIN_FLIP' or 'ANY' and not '{}.".format(kind.upper()))
+ elif kind not in ["singlet", "triplet", "any"]:
+ raise ValidationError("For RHF references the value '{}' for 'KIND' is "
+ "not supported.".format(kind.upper()))
+ kwargs["kind"] = kind
+ kwargs["max_iter"] = core.get_option("ADC", "MAXITER")
+
+ #
+ # Determine ADC function method from adcc to run ADC
+ #
+ adcrunner = {
+ "cvs-adc(1)": adcc.cvs_adc1, "cvs-adc(2)": adcc.cvs_adc2,
+ "cvs-adc(2)-x": adcc.cvs_adc2x, "cvs-adc(3)": adcc.cvs_adc3,
+ "adc(1)": adcc.adc1, "adc(2)": adcc.adc2,
+ "adc(2)-x": adcc.adc2x, "adc(3)": adcc.adc3,
+ }
+ if name not in adcrunner:
+ raise ValidationError(f"Unsupported ADC method: {name}")
+ if "cvs" in name and "core_orbitals" not in kwargs:
+ raise ValidationError("If a CVS-ADC method is requested, the NUM_CORE_ORBITALS option "
+ "needs to be set.")
+ if "core_orbitals" in kwargs and "cvs" not in name:
+ raise ValidationError("The NUM_CORE_ORBITALS option needs to be set to '0' or absent "
+ "unless a CVS ADC method is requested.")
+ if "cvs" in name and kwargs["kind"] in ["spin_flip"]:
+ raise ValidationError("Spin-flip for CVS-ADC variants is not available.")
+
+ #
+ # Launch the rocket
+ #
+ # Copy thread setup from psi4
+ adcc.set_n_threads(core.get_num_threads())
+
+ # Hack to direct the stream-like interface adcc expects to the string interface of Psi4 core
+ class CoreStream:
+ def write(self, text):
+ core.print_out(text)
+
+ core.print_out("\n" + adcc.banner(colour=False) + "\n")
+ try:
+ state = adcrunner[name](ref_wfn, **kwargs, output=CoreStream())
+ except InvalidReference as ex:
+ raise ValidationError("Cannot run adcc because the passed reference wavefunction is "
+ "not supported in adcc. Check Psi4 SCF parameters. adcc reports: "
+ f"{ex}")
+ except Exception as ex:
+ raise ValidationError("Unknown exception occured while "
+ f"running adcc: '{ex}' ({type(ex).__name__})")
+ core.print_out("\n")
+
+ # TODO Should a non-converged calculation throw?
+
+ #
+ # Interpret results
+ #
+ # Note: This wavefunction is not consistent ... the density
+ # is e.g. not the proper one (i.e. not the MP(n) one)
+ adc_wfn = core.Wavefunction(ref_wfn.molecule(), ref_wfn.basisset())
+ adc_wfn.shallow_copy(ref_wfn)
+ adc_wfn.set_reference_wavefunction(ref_wfn)
+ adc_wfn.set_name(name)
+ adc_wfn.set_module("adcc")
+
+ # MP(3) energy for CVS-ADC(3) calculations is still a missing feature in adcc
+ # ... we store this variant here to be able to fall back to MP(2) energies.
+ is_cvs_adc3 = state.method.level >= 3 and state.ground_state.has_core_occupied_space
+
+ # Ground-state energies
+ mp = state.ground_state
+ mp_energy = mp.energy(state.method.level if not is_cvs_adc3 else 2)
+ mp_corr = 0.0
+ if state.method.level > 1:
+ core.print_out("Ground state energy breakdown:\n")
+ core.print_out(" Energy SCF {0:15.8g} [Eh]\n".format(ref_wfn.energy()))
+ for level in range(2, state.method.level + 1):
+ if level >= 3 and is_cvs_adc3:
+ continue
+ energy = mp.energy_correction(level)
+ mp_corr += energy
+ adc_wfn.set_variable(f"MP{level} CORRELATION ENERGY", energy)
+ adc_wfn.set_variable(f"MP{level} TOTAL ENERGY", mp.energy(level))
+ core.print_out(f" Energy correlation MP{level} {energy:15.8g} [Eh]\n")
+ core.print_out(" Energy total {0:15.8g} [Eh]\n".format(mp_energy))
+ adc_wfn.set_variable("CURRENT CORRELATION ENERGY", mp_corr) # P::e ADC
+ adc_wfn.set_variable("CURRENT ENERGY", mp_energy) # P::e ADC
+
+ # Set results of excited-states computation
+ # TODO Does not work: Can't use strings
+ # adc_wfn.set_variable("excitation kind", state.kind)
+ adc_wfn.set_variable("number of excited states", len(state.excitation_energy))
+ adc_wfn.set_variable("ADC ITERATIONS", state.n_iter) # P::e ADC
+ methods = [name.upper(), 'ADC']
+ for excitation in state.excitations:
+ root_index = excitation.index + 1
+ for method in methods:
+ adc_wfn.set_variable(f"{method} ROOT 0 (A) -> ROOT {root_index} (A) EXCITATION ENERGY",
+ excitation.excitation_energy)
+ adc_wfn.set_variable(f"{method} ROOT 0 (IN A) -> ROOT {root_index} (IN A) EXCITATION ENERGY",
+ excitation.excitation_energy)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} EXCITATION ENERGY",
+ excitation.excitation_energy)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} EXCITATION ENERGY - A TRANSITION",
+ excitation.excitation_energy)
+
+ core.print_out("\n\n ==> Excited states summary <== \n")
+ core.print_out("\n" + state.describe(oscillator_strengths=False) + "\n")
+
+ # TODO Setting the excitation amplitude elements inside the wavefunction is a little
+ # challenging, since for each excitation vector one needs to extract the elements
+ # and map the indices from the adcc to the Psi4 convention. For this reason it
+ # is not yet done.
+
+ core.print_out("\n ==> Dominant amplitudes per state <== \n\n")
+ tol_ampl = core.get_option("ADC", "CUTOFF_AMPS_PRINT")
+ core.print_out(state.describe_amplitudes(tolerance=tol_ampl) + "\n\n")
+
+ # Shove variables into global space
+ for k, v in adc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ if do_timer:
+ core.tstop()
+ adc_wfn.adcc_state = state
+ return adc_wfn
+
+
+def run_adcc_property(name, **kwargs):
+ """Run a ADC excited-states property calculation in adcc
+ and return the resulting properties.
+
+ """
+ # TODO Things available in ADCC, but not yet implemented here:
+ # Export of difference and transition density matrices for all states
+
+ properties = [prop.upper() for prop in kwargs.pop('properties')]
+ valid_properties = ['DIPOLE', 'OSCILLATOR_STRENGTH', 'TRANSITION_DIPOLE',
+ 'ROTATIONAL_STRENGTH']
+ unknown_properties = [prop for prop in properties if prop not in valid_properties]
+
+ if unknown_properties:
+ alternatives = ""
+ alt_method_name = p4util.text.find_approximate_string_matches(unknown_properties[0],
+ valid_properties, 2)
+ if alt_method_name:
+ alternatives = " Did you mean? " + " ".join(alt_method_name)
+
+ raise ValidationError("ADC property: Feature '{}' is not recognized. {}"
+ "".format(unknown_properties[0], alternatives))
+
+ # Start timer
+ do_timer = kwargs.pop("do_timer", True)
+ if do_timer:
+ core.tstart()
+ adc_wfn = run_adcc(name, do_timer=False, **kwargs)
+ state = adc_wfn.adcc_state
+ hf = state.reference_state
+ mp = state.ground_state
+
+ # Formats and indention
+ ind = " "
+ def format_vector(label, data):
+ assert data.ndim == 1
+ return f"{label:<40s} " + " ".join(f"{d:12.6g}" for d in data)
+
+ if "DIPOLE" in properties:
+ lines = ["\nGround state properties"]
+ lines += [ind + "Hartree-Fock (HF)"]
+ lines += [ind + ind + format_vector("Dipole moment (in a.u.)", hf.dipole_moment)]
+
+ if state.method.level > 1:
+ lines += [ind + "Møller Plesset 2nd order (MP2)"]
+ lines += [ind + ind + format_vector("Dipole moment (in a.u.)", mp.dipole_moment(2))]
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ for i, cart in enumerate(["X", "Y", "Z"]):
+ # retire components at v1.5
+ adc_wfn.set_variable("MP2 dipole " + cart, mp.dipole_moment(2)[i])
+ adc_wfn.set_variable("current dipole " + cart, mp.dipole_moment(2)[i])
+ adc_wfn.set_variable("MP2 dipole", mp.dipole_moment(2))
+ adc_wfn.set_variable("current dipole", mp.dipole_moment(2))
+ lines += [""]
+ core.print_out("\n".join(lines) + "\n")
+
+ gauge = core.get_option("ADC", "GAUGE").lower()
+ if gauge == "velocity":
+ gauge_short = "VEL"
+ elif gauge == "length":
+ gauge_short = "LEN"
+ else:
+ raise ValidationError(f"Gauge {gauge} not recognised for ADC calculations.")
+
+ computed = []
+ methods = [name.upper(), 'ADC']
+ for excitation in state.excitations:
+ root_index = excitation.index + 1
+ props = {}
+ if any(prop in properties for prop in ("TRANSITION_DIPOLE", "OSCILLATOR_STRENGTH")):
+ data = excitation.transition_dipole_moment
+ props["Transition dipole moment (in a.u.)"] = data
+ data_mat = data.reshape(1, 3)
+ for method in methods:
+ adc_wfn.set_variable(f"{method} ROOT 0 (A) -> ROOT {root_index} (A) "
+ "ELECTRIC TRANSITION DIPOLE MOMENT (LEN)", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT 0 (IN A) -> ROOT {root_index} (IN A) "
+ "ELECTRIC TRANSITION DIPOLE MOMENT (LEN)", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ "ELECTRIC TRANSITION DIPOLE MOMENT (LEN)", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ "ELECTRIC TRANSITION DIPOLE MOMENT (LEN) - A TRANSITION", data_mat)
+
+ if "OSCILLATOR_STRENGTH" in properties:
+ if gauge == "velocity":
+ data = excitation.oscillator_strength_velocity
+ else:
+ data = excitation.oscillator_strength
+ props[f"Oscillator strength ({gauge} gauge)"] = data.reshape(-1)
+ for method in methods:
+ adc_wfn.set_variable(f"{method} ROOT 0 (A) -> ROOT {root_index} (A) "
+ f"OSCILLATOR STRENGTH ({gauge_short})", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 (IN A) -> ROOT {root_index} (IN A) "
+ f"OSCILLATOR STRENGTH ({gauge_short})", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ f"OSCILLATOR STRENGTH ({gauge_short})", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ f"OSCILLATOR STRENGTH ({gauge_short}) - A TRANSITION", data)
+
+ if "ROTATIONAL_STRENGTH" in properties:
+ data = excitation.rotatory_strength
+ props["Rotational strength (velocity gauge)"] = data.reshape(-1)
+ for method in methods:
+ adc_wfn.set_variable(f"{method} ROOT 0 (A) -> ROOT {root_index} (A) "
+ "ROTATORY STRENGTH (VEL)", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 (IN A) -> ROOT {root_index} (IN A) "
+ "ROTATORY STRENGTH (VEL)", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ "ROTATORY STRENGTH (VEL)", data)
+ adc_wfn.set_variable(f"{method} ROOT 0 -> ROOT {root_index} "
+ "ROTATORY STRENGTH (VEL) - A TRANSITION", data)
+
+ if "DIPOLE" in properties:
+ data = excitation.state_dipole_moment
+ props["State dipole moment (in a.u.)"] = data
+ data_mat = data.reshape(1, 3)
+ for method in methods:
+ adc_wfn.set_variable(f"{method} ROOT {root_index} DIPOLE MOMENT", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT {root_index} DIPOLE MOMENT - A TRANSITION", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT {root_index} (A) DIPOLE MOMENT", data_mat)
+ adc_wfn.set_variable(f"{method} ROOT {root_index} (IN A) DIPOLE MOMENT", data_mat)
+
+ computed.append(props)
+
+ # for Psivar scraper
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n EXCITATION ENERGY") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (IN h) -> ROOT n (IN i) EXCITATION ENERGY") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (h) -> ROOT n (i) EXCITATION ENERGY") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n EXCITATION ENERGY - h TRANSITION") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT 0 -> ROOT n ELECTRIC TRANSITION DIPOLE MOMENT (LEN)") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT 0 (IN h) -> ROOT n (IN i) ELECTRIC TRANSITION DIPOLE MOMENT (LEN)") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT 0 (h) -> ROOT n (i) ELECTRIC TRANSITION DIPOLE MOMENT (LEN)") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT 0 -> ROOT n ELECTRIC TRANSITION DIPOLE MOMENT (LEN) - h TRANSITION") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT n DIPOLE MOMENT") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT n (IN i) DIPOLE MOMENT") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT n (i) DIPOLE MOMENT") # P::e ADC
+ # wfn.set_array_variable("ADC ROOT n DIPOLE MOMENT - h TRANSITION") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n OSCILLATOR STRENGTH (LEN)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (IN h) -> ROOT n (IN i) OSCILLATOR STRENGTH (LEN)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (h) -> ROOT n (i) OSCILLATOR STRENGTH (LEN)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n OSCILLATOR STRENGTH (LEN) - h TRANSITION") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n OSCILLATOR STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (IN h) -> ROOT n (IN i) OSCILLATOR STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (h) -> ROOT n (i) OSCILLATOR STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n OSCILLATOR STRENGTH (VEL) - h TRANSITION") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n ROTATORY STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (IN h) -> ROOT n (IN i) ROTATORY STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 (h) -> ROOT n (i) ROTATORY STRENGTH (VEL)") # P::e ADC
+ # wfn.set_variable("ADC ROOT 0 -> ROOT n ROTATORY STRENGTH (VEL) - h TRANSITION") # P::e ADC
+
+ core.print_out("\nExcited state properties:\n")
+ for i, props in enumerate(computed):
+ lines = [ind + f"Excited state {i}"]
+ for prop, data in sorted(props.items()):
+ lines += [ind + ind + format_vector(prop, data)]
+ core.print_out("\n".join(lines) + "\n")
+
+ # Shove variables into global space
+ for k, v in adc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ if do_timer:
+ core.tstop()
+ return adc_wfn
+
+
+def run_detci(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a configuration interaction calculation, namely FCI,
+ CIn, MPn, and ZAPTn.
+
+ """
+ # dtl = docs_table_link("dummy", "detci")
+
+ optstash = p4util.OptionsState(
+ ['DETCI', 'WFN'],
+ ['DETCI', 'MAX_NUM_VECS'],
+ ['DETCI', 'MPN_ORDER_SAVE'],
+ ['DETCI', 'MPN'],
+ ['DETCI', 'FCI'],
+ ['DETCI', 'EX_LEVEL'])
+
+ # throw exception for UHF
+ if core.get_option('DETCI', 'REFERENCE') not in ['RHF', 'ROHF']:
+ raise ValidationError('Reference %s for DETCI is not available.' %
+ core.get_option('DETCI', 'REFERENCE'))
+
+ # throw exception for DF/CD. many of these pre-trapped by select_* functions but some escape, incl. zapt
+ if (corl_type := method_algorithm_type(name).now) != "CONV":
+ raise ValidationError(f"Invalid type {corl_type} for DETCI energy through `run_detci`.") # See Capabilities Table")
+
+ mtdlvl_mobj = re.match(r"""\A(?P<method>[a-z]+)(?P<level>\d+)\Z""", name.lower())
+
+ if mtdlvl_mobj and mtdlvl_mobj.group("method") == "zapt":
+ level = int(mtdlvl_mobj.group("level"))
+
+ # throw exception for non-ROHF
+ if (ref := core.get_option("SCF", "REFERENCE")) != "ROHF":
+ raise UpgradeHelper(f"energy('zapt{level}')", f"energy('mp{level}')", 1.7,
+ " Replace method ZAPT with method MP for RHF reference. DETCI is orders-of-magnitude inefficient for perturbation theory.")
+
+ core.set_local_option('DETCI', 'WFN', 'ZAPTN')
+ maxnvect = int((level + 1) / 2) + (level + 1) % 2
+ core.set_local_option('DETCI', 'MAX_NUM_VECS', maxnvect)
+ if (level + 1) % 2:
+ core.set_local_option('DETCI', 'MPN_ORDER_SAVE', 2)
+ else:
+ core.set_local_option('DETCI', 'MPN_ORDER_SAVE', 1)
+ elif mtdlvl_mobj and mtdlvl_mobj.group("method") == "mp":
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'MPN', 'TRUE')
+ level = int(mtdlvl_mobj.group("level"))
+ maxnvect = int((level + 1) / 2) + (level + 1) % 2
+ core.set_local_option('DETCI', 'MAX_NUM_VECS', maxnvect)
+ if (level + 1) % 2:
+ core.set_local_option('DETCI', 'MPN_ORDER_SAVE', 2)
+ else:
+ core.set_local_option('DETCI', 'MPN_ORDER_SAVE', 1)
+ elif name == 'ccsd':
+ # untested
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'CC', 'TRUE')
+ core.set_local_option('DETCI', 'CC_EX_LEVEL', 2)
+ elif name == 'fci':
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'FCI', 'TRUE')
+ elif name == 'cisd':
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'EX_LEVEL', 2)
+ elif name == 'cisdt':
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'EX_LEVEL', 3)
+ elif name == 'cisdtq':
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ core.set_local_option('DETCI', 'EX_LEVEL', 4)
+ elif mtdlvl_mobj and mtdlvl_mobj.group("method") == "ci":
+ core.set_local_option('DETCI', 'WFN', 'DETCI')
+ level = int(mtdlvl_mobj.group("level"))
+ core.set_local_option('DETCI', 'EX_LEVEL', level)
+ elif name == 'detci':
+ pass
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ ciwfn = core.detci(ref_wfn)
+
+ # Shove variables into global space
+ for k, v in ciwfn.variables().items():
+ core.set_variable(k, v)
+
+ print_nos = False
+ if core.get_option("DETCI", "NAT_ORBS"):
+ ciwfn.ci_nat_orbs()
+ print_nos = True
+
+ proc_util.print_ci_results(ciwfn, name.upper(), ciwfn.variable("HF TOTAL ENERGY"), ciwfn.variable("CURRENT ENERGY"), print_nos)
+
+ core.print_out("\t\t \"A good bug is a dead bug\" \n\n")
+ core.print_out("\t\t\t - Starship Troopers\n\n")
+ core.print_out("\t\t \"I didn't write FORTRAN. That's the problem.\"\n\n")
+ core.print_out("\t\t\t - Edward Valeev\n")
+
+ if core.get_global_option("DIPMOM") and ("mp" not in name.lower()):
+ # We always would like to print a little dipole information
+ oeprop = core.OEProp(ciwfn)
+ oeprop.set_title(name.upper())
+ oeprop.add("DIPOLE")
+ oeprop.compute()
+ ciwfn.oeprop = oeprop
+ core.set_variable("CURRENT DIPOLE", core.variable(name.upper() + " DIPOLE"))
+
+ ciwfn.cleanup_ci()
+ ciwfn.cleanup_dpd()
+ _clean_detci()
+
+ for lvl in range(4, 11):
+ if ciwfn.has_variable(f"MP{lvl} CORRELATION ENERGY") and ciwfn.has_variable(f"MP{lvl-1} CORRELATION ENERGY"):
+ ciwfn.set_variable(f"MP{lvl} CORRECTION ENERGY", ciwfn.variable(f"MP{lvl} CORRELATION ENERGY") - ciwfn.variable(f"MP{lvl-1} CORRELATION ENERGY"))
+ core.set_variable(f"MP{lvl} CORRECTION ENERGY", ciwfn.variable(f"MP{lvl} CORRELATION ENERGY") - ciwfn.variable(f"MP{lvl-1} CORRELATION ENERGY"))
+
+ optstash.restore()
+ return ciwfn
+
+
+def run_dfmp2(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density-fitted MP2 calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['DF_BASIS_MP2'],
+ ['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+ core.print_out(""" SCF Algorithm Type (re)set to DF.\n""")
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ core.tstart()
+ core.print_out('\n')
+ p4util.banner('DFMP2')
+ core.print_out('\n')
+
+ if core.get_global_option('REFERENCE') == "ROHF":
+ ref_wfn.semicanonicalize()
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_MP2",
+ core.get_option("DFMP2", "DF_BASIS_MP2"),
+ "RIFIT", core.get_global_option('BASIS'))
+ ref_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
+
+ dfmp2_wfn = core.dfmp2(ref_wfn)
+ dfmp2_wfn.compute_energy()
+
+ if name == 'scs-mp2':
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('SCS-MP2 TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('SCS-MP2 CORRELATION ENERGY'))
+
+ elif name == 'mp2':
+ dfmp2_wfn.set_variable('CURRENT ENERGY', dfmp2_wfn.variable('MP2 TOTAL ENERGY'))
+ dfmp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dfmp2_wfn.variable('MP2 CORRELATION ENERGY'))
+
+ # Shove variables into global space
+ for k, v in dfmp2_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ core.tstop()
+ return dfmp2_wfn
+
+
+def run_dfep2(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a density-fitted MP2 calculation.
+
+ """
+ core.tstart()
+ optstash = p4util.OptionsState(
+ ['DF_BASIS_MP2'],
+ ['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+ core.print_out(""" SCF Algorithm Type (re)set to DF.\n""")
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ if core.get_global_option('REFERENCE') != "RHF":
+ raise ValidationError("DF-EP2 is not available for %s references.",
+ core.get_global_option('REFERENCE'))
+
+
+ # Build the wavefunction
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_EP2",
+ core.get_option("DFEP2", "DF_BASIS_EP2"),
+ "RIFIT", core.get_global_option('BASIS'))
+ ref_wfn.set_basisset("DF_BASIS_EP2", aux_basis)
+
+ dfep2_wfn = core.DFEP2Wavefunction(ref_wfn)
+
+ # Figure out what were doing
+ if core.has_option_changed('DFEP2', 'EP2_ORBITALS'):
+ ep2_input = core.get_global_option("EP2_ORBITALS")
+
+ else:
+ n_ip = core.get_global_option("EP2_NUM_IP")
+ n_ea = core.get_global_option("EP2_NUM_EA")
+
+ eps = np.hstack(dfep2_wfn.epsilon_a().nph)
+ irrep_map = np.hstack([np.ones_like(dfep2_wfn.epsilon_a().nph[x]) * x for x in range(dfep2_wfn.nirrep())])
+ sort = np.argsort(eps)
+
+ ip_map = sort[dfep2_wfn.nalpha() - n_ip:dfep2_wfn.nalpha()]
+ ea_map = sort[dfep2_wfn.nalpha():dfep2_wfn.nalpha() + n_ea]
+
+ ep2_input = [[] for x in range(dfep2_wfn.nirrep())]
+ nalphapi = tuple(dfep2_wfn.nalphapi())
+
+ # Add IP info
+ ip_info = np.unique(irrep_map[ip_map], return_counts=True)
+ for irrep, cnt in zip(*ip_info):
+ irrep = int(irrep)
+ ep2_input[irrep].extend(range(nalphapi[irrep] - cnt, nalphapi[irrep]))
+
+ # Add EA info
+ ea_info = np.unique(irrep_map[ea_map], return_counts=True)
+ for irrep, cnt in zip(*ea_info):
+ irrep = int(irrep)
+ ep2_input[irrep].extend(range(nalphapi[irrep], nalphapi[irrep] + cnt))
+
+ # Compute
+ ret = dfep2_wfn.compute(ep2_input)
+
+ # Resort it...
+ ret_eps = []
+ for h in range(dfep2_wfn.nirrep()):
+ ep2_data = ret[h]
+ inp_data = ep2_input[h]
+
+ for i in range(len(ep2_data)):
+ tmp = [h, ep2_data[i][0], ep2_data[i][1], dfep2_wfn.epsilon_a().get(h, inp_data[i]), inp_data[i]]
+ ret_eps.append(tmp)
+
+ ret_eps.sort(key=lambda x: x[3])
+
+ h2ev = constants.hartree2ev
+ irrep_labels = dfep2_wfn.molecule().irrep_labels()
+
+ core.print_out(" ==> Results <==\n\n")
+ core.print_out(" %8s %12s %12s %8s\n" % ("Orbital", "Koopmans (eV)", "EP2 (eV)", "EP2 PS"))
+ core.print_out(" ----------------------------------------------\n")
+ for irrep, ep2, ep2_ps, kt, pos in ret_eps:
+ label = str(pos + 1) + irrep_labels[irrep]
+ core.print_out(" %8s % 12.3f % 12.3f % 6.3f\n" % (label, (kt * h2ev), (ep2 * h2ev), ep2_ps))
+ core.set_variable("EP2 " + label.upper() + " ENERGY", ep2)
+ core.print_out(" ----------------------------------------------\n\n")
+
+ # Figure out the IP and EA
+ sorted_vals = np.array([x[1] for x in ret_eps])
+ ip_vals = sorted_vals[sorted_vals < 0]
+ ea_vals = sorted_vals[sorted_vals > 0]
+
+ ip_value = None
+ ea_value = None
+ if len(ip_vals):
+ core.set_variable("EP2 IONIZATION POTENTIAL", ip_vals[-1])
+ core.set_variable("CURRENT ENERGY", ip_vals[-1])
+ if len(ea_vals):
+ core.set_variable("EP2 ELECTRON AFFINITY", ea_vals[0])
+ if core.variable("EP2 IONIZATION POTENTIAL") == 0.0:
+ core.set_variable("CURRENT ENERGY", ea_vals[0])
+
+ core.print_out(" EP2 has completed successfully!\n\n")
+
+ core.tstop()
+ return dfep2_wfn
+
+
+def run_dlpnomp2(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a DLPNO-MP2 calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['DF_BASIS_MP2'],
+ ['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+ core.print_out(""" SCF Algorithm Type (re)set to DF.\n""")
+
+ # DLPNO-MP2 is only DF
+ if core.get_global_option('MP2_TYPE') != "DF":
+ raise ValidationError(""" DLPNO-MP2 is only implemented with density fitting.\n"""
+ """ 'mp2_type' must be set to 'DF'.\n""")
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ elif ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" DLPNO-MP2 does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ if core.get_global_option('REFERENCE') != "RHF":
+ raise ValidationError("DLPNO-MP2 is not available for %s references.",
+ core.get_global_option('REFERENCE'))
+
+ core.tstart()
+ core.print_out('\n')
+ p4util.banner('DLPNO-MP2')
+ core.print_out('\n')
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_MP2",
+ core.get_option("DLPNO", "DF_BASIS_MP2"),
+ "RIFIT", core.get_global_option('BASIS'))
+ ref_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
+
+ dlpnomp2_wfn = core.dlpno(ref_wfn)
+ dlpnomp2_wfn.compute_energy()
+
+ if name == 'scs-dlpno-mp2':
+ dlpnomp2_wfn.set_variable('CURRENT ENERGY', dlpnomp2_wfn.variable('SCS-MP2 TOTAL ENERGY'))
+ dlpnomp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dlpnomp2_wfn.variable('SCS-MP2 CORRELATION ENERGY'))
+
+ elif name == 'dlpno-mp2':
+ dlpnomp2_wfn.set_variable('CURRENT ENERGY', dlpnomp2_wfn.variable('MP2 TOTAL ENERGY'))
+ dlpnomp2_wfn.set_variable('CURRENT CORRELATION ENERGY', dlpnomp2_wfn.variable('MP2 CORRELATION ENERGY'))
+
+ # Shove variables into global space
+ for k, v in dlpnomp2_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ core.tstop()
+ return dlpnomp2_wfn
+
+
+def run_dmrgscf(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an DMRG calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['SCF_TYPE'],
+ ['DMRG', 'DMRG_CASPT2_CALC'])
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs)
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ if 'CASPT2' in name.upper():
+ core.set_local_option("DMRG", "DMRG_CASPT2_CALC", True)
+
+ dmrg_wfn = core.dmrg(ref_wfn)
+ optstash.restore()
+
+ # Shove variables into global space
+ for k, v in dmrg_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return dmrg_wfn
+
+
+def run_dmrgci(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an DMRG calculation.
+
+ """
+ optstash = p4util.OptionsState(
+ ['SCF_TYPE'],
+ ['DMRG', 'DMRG_SCF_MAX_ITER'])
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs)
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ core.set_local_option('DMRG', 'DMRG_SCF_MAX_ITER', 1)
+
+ dmrg_wfn = core.dmrg(ref_wfn)
+ optstash.restore()
+
+ # Shove variables into global space
+ for k, v in dmrg_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return dmrg_wfn
+
+
+def run_psimrcc(name, **kwargs):
+ """Function encoding sequence of PSI module calls for a PSIMRCC computation
+ using a reference from the MCSCF module
+
+ """
+ mcscf_wfn = run_mcscf(name, **kwargs)
+ psimrcc_wfn = core.psimrcc(mcscf_wfn)
+
+ # Shove variables into global space
+ for k, v in psimrcc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return psimrcc_wfn
+
+
+def run_psimrcc_scf(name, **kwargs):
+ """Function encoding sequence of PSI module calls for a PSIMRCC computation
+ using a reference from the SCF module
+
+ """
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs)
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+
+ psimrcc_wfn = core.psimrcc(ref_wfn)
+
+ # Shove variables into global space
+ for k, v in psimrcc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ return psimrcc_wfn
+
+
+def run_sapt(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a SAPT calculation of any level.
+
+ """
+ optstash = p4util.OptionsState(['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ # Get the molecule of interest
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ sapt_dimer = kwargs.pop('molecule', core.get_active_molecule())
+ else:
+ core.print_out('Warning! SAPT argument "ref_wfn" is only able to use molecule information.')
+ sapt_dimer = ref_wfn.molecule()
+
+ sapt_basis = kwargs.pop('sapt_basis', 'dimer')
+
+ sapt_dimer, monomerA, monomerB = proc_util.prepare_sapt_molecule(sapt_dimer, sapt_basis)
+
+ # Need to ensure consistent orbital freezing
+ # between monomer and dimer computations
+ monomerA_basis = core.BasisSet.build(monomerA, "BASIS", core.get_global_option("BASIS"))
+ monomerB_basis = core.BasisSet.build(monomerB, "BASIS", core.get_global_option("BASIS"))
+ nfc_ab = monomerA_basis.n_frozen_core() + monomerB_basis.n_frozen_core()
+
+ if (core.get_option('SCF', 'REFERENCE') != 'RHF') and (name.upper() != "SAPT0"):
+ raise ValidationError('Only SAPT0 supports a reference different from \"reference rhf\".')
+
+ do_delta_mp2 = True if name.endswith('dmp2') else False
+ do_empirical_disp = True if '-d' in name.lower() else False
+
+ if do_empirical_disp:
+ ## Make sure we are turning SAPT0 dispersion off
+ core.set_local_option('SAPT', 'SAPT0_E10', True)
+ core.set_local_option('SAPT', 'SAPT0_E20IND', True)
+ core.set_local_option('SAPT', 'SAPT0_E20Disp', False)
+
+ ri = core.get_global_option('SCF_TYPE')
+ df_ints_io = core.get_option('SCF', 'DF_INTS_IO')
+ # inquire if above at all applies to dfmp2
+
+ core.IO.set_default_namespace('dimer')
+ core.print_out('\n')
+ p4util.banner('Dimer HF')
+ core.print_out('\n')
+
+ # Compute dimer wavefunction
+
+ if (sapt_basis == 'dimer') and (ri == 'DF'):
+ core.set_global_option('DF_INTS_IO', 'SAVE')
+
+ optstash2 = p4util.OptionsState(['NUM_FROZEN_DOCC'])
+ core.set_global_option("NUM_FROZEN_DOCC", nfc_ab)
+ core.timer_on("SAPT: Dimer SCF")
+ dimer_wfn = scf_helper('RHF', molecule=sapt_dimer, **kwargs)
+ core.timer_off("SAPT: Dimer SCF")
+
+ necp_ab = dimer_wfn.basisset().n_ecp_core()
+ share_df_ints = ((sapt_basis == 'dimer') and (ri == 'DF') and not (necp_ab and (os.name == 'nt')))
+ if (sapt_basis == 'dimer') and (ri == 'DF') and not share_df_ints:
+ core.print_out(f"\n Turning off SAPT DF integrals sharing because of ECP: {necp_ab}\n\n")
+
+ if do_delta_mp2:
+ select_mp2("mp2", ref_wfn=dimer_wfn, **kwargs)
+ mp2_corl_interaction_e = core.variable('MP2 CORRELATION ENERGY')
+
+ optstash2.restore()
+ if share_df_ints:
+ core.set_global_option('DF_INTS_IO', 'LOAD')
+
+ # Compute Monomer A wavefunction
+ if share_df_ints:
+ core.IO.change_file_namespace(97, 'dimer', 'monomerA')
+
+ core.IO.set_default_namespace('monomerA')
+ core.print_out('\n')
+ p4util.banner('Monomer A HF')
+ core.print_out('\n')
+
+ core.timer_on("SAPT: Monomer A SCF")
+ monomerA_wfn = scf_helper('RHF', molecule=monomerA, **kwargs)
+ core.timer_off("SAPT: Monomer A SCF")
+
+ if do_delta_mp2:
+ select_mp2("mp2", ref_wfn=monomerA_wfn, **kwargs)
+ mp2_corl_interaction_e -= core.variable('MP2 CORRELATION ENERGY')
+
+ # Compute Monomer B wavefunction
+ if share_df_ints:
+ core.IO.change_file_namespace(97, 'monomerA', 'monomerB')
+ core.IO.set_default_namespace('monomerB')
+ core.print_out('\n')
+ p4util.banner('Monomer B HF')
+ core.print_out('\n')
+
+ core.timer_on("SAPT: Monomer B SCF")
+ monomerB_wfn = scf_helper('RHF', molecule=monomerB, **kwargs)
+ core.timer_off("SAPT: Monomer B SCF")
+
+ # Delta MP2
+ if do_delta_mp2:
+ select_mp2("mp2", ref_wfn=monomerB_wfn, **kwargs)
+ mp2_corl_interaction_e -= core.variable('MP2 CORRELATION ENERGY')
+ core.set_variable("SAPT MP2 CORRELATION ENERGY", mp2_corl_interaction_e) # P::e SAPT
+ core.set_global_option('DF_INTS_IO', df_ints_io)
+
+ if core.get_option('SCF', 'REFERENCE') == 'RHF':
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERA, 'monomerA', 'dimer')
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERB, 'monomerB', 'dimer')
+
+ core.IO.set_default_namespace('dimer')
+ core.set_local_option('SAPT', 'E_CONVERGENCE', 10e-10)
+ core.set_local_option('SAPT', 'D_CONVERGENCE', 10e-10)
+ if name in ['sapt0', 'ssapt0']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT0')
+ elif name == 'sapt2':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2')
+ elif name in ['sapt2+', 'sapt2+dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+')
+ core.set_local_option('SAPT', 'DO_CCD_DISP', False)
+ elif name in ['sapt2+(3)', 'sapt2+(3)dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', False)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', False)
+ elif name in ['sapt2+3', 'sapt2+3dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', True)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', False)
+ elif name in ['sapt2+(ccd)', 'sapt2+(ccd)dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+')
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+ elif name in ['sapt2+(3)(ccd)', 'sapt2+(3)(ccd)dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', False)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+ elif name in ['sapt2+3(ccd)', 'sapt2+3(ccd)dmp2']:
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', True)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+
+ # Make sure we are not going to run CPHF on ROHF, since its MO Hessian
+ # is not SPD
+ if core.get_option('SCF', 'REFERENCE') == 'ROHF':
+ core.set_local_option('SAPT', 'COUPLED_INDUCTION', False)
+ core.print_out(' Coupled induction not available for ROHF.\n')
+ core.print_out(' Proceeding with uncoupled induction only.\n')
+
+ core.print_out(" Constructing Basis Sets for SAPT...\n\n")
+ aux_basis = core.BasisSet.build(dimer_wfn.molecule(), "DF_BASIS_SAPT", core.get_global_option("DF_BASIS_SAPT"),
+ "RIFIT", core.get_global_option("BASIS"))
+ dimer_wfn.set_basisset("DF_BASIS_SAPT", aux_basis)
+
+ aux_basis = core.BasisSet.build(dimer_wfn.molecule(), "DF_BASIS_ELST", core.get_global_option("DF_BASIS_ELST"),
+ "JKFIT", core.get_global_option("BASIS"))
+ dimer_wfn.set_basisset("DF_BASIS_ELST", aux_basis)
+
+ core.print_out('\n')
+ p4util.banner(name.upper())
+ core.print_out('\n')
+ e_sapt = core.sapt(dimer_wfn, monomerA_wfn, monomerB_wfn)
+ dimer_wfn.set_module("sapt")
+
+ from psi4.driver.qcdb.psivardefs import sapt_psivars
+ p4util.expand_psivars(sapt_psivars())
+ optstash.restore()
+
+ # Get the SAPT name right if doing empirical dispersion
+ if do_empirical_disp:
+ sapt_name = "sapt0"
+ else:
+ sapt_name = name
+
+ # Make sure we got induction, otherwise replace it with uncoupled induction
+ which_ind = 'IND'
+ target_ind = 'IND'
+ if not core.has_variable(' '.join((sapt_name.upper(), which_ind, 'ENERGY'))):
+ which_ind = 'IND,U'
+
+ for term in ['ELST', 'EXCH', 'DISP', 'TOTAL']:
+ core.set_variable(' '.join(['SAPT', term, 'ENERGY']),
+ core.variable(' '.join([sapt_name.upper(), term, 'ENERGY'])))
+ # Special induction case
+ core.set_variable(' '.join(['SAPT', target_ind, 'ENERGY']),
+ core.variable(' '.join([sapt_name.upper(), which_ind, 'ENERGY'])))
+ core.set_variable('CURRENT ENERGY', core.variable('SAPT TOTAL ENERGY'))
+
+ # Empirical dispersion
+ if do_empirical_disp:
+ proc_util.sapt_empirical_dispersion(name, dimer_wfn)
+
+ return dimer_wfn
+
+
+def run_sapt_ct(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a charge-transfer SAPT calcuation of any level.
+
+ """
+ optstash = p4util.OptionsState(
+ ['SCF_TYPE'])
+
+ if 'ref_wfn' in kwargs:
+ core.print_out('\nWarning! Argument ref_wfn is not valid for sapt computations\n')
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ # Get the molecule of interest
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ sapt_dimer = kwargs.pop('molecule', core.get_active_molecule())
+ else:
+ core.print_out('Warning! SAPT argument "ref_wfn" is only able to use molecule information.')
+ sapt_dimer = ref_wfn.molecule()
+
+ sapt_dimer, monomerA, monomerB = proc_util.prepare_sapt_molecule(sapt_dimer, "dimer")
+ monomerAm = sapt_dimer.extract_subsets(1)
+ monomerAm.set_name('monomerAm')
+ monomerBm = sapt_dimer.extract_subsets(2)
+ monomerBm.set_name('monomerBm')
+
+ if core.get_option('SCF', 'REFERENCE') != 'RHF':
+ raise ValidationError('SAPT requires requires \"reference rhf\".')
+
+ ri = core.get_global_option('SCF_TYPE')
+ df_ints_io = core.get_option('SCF', 'DF_INTS_IO')
+ # inquire if above at all applies to dfmp2
+
+ core.IO.set_default_namespace('dimer')
+ core.print_out('\n')
+ p4util.banner('Dimer HF')
+ core.print_out('\n')
+ core.set_global_option('DF_INTS_IO', 'SAVE')
+ dimer_wfn = scf_helper('RHF', molecule=sapt_dimer, **kwargs)
+ core.set_global_option('DF_INTS_IO', 'LOAD')
+
+ if (ri == 'DF'):
+ core.IO.change_file_namespace(97, 'dimer', 'monomerA')
+ core.IO.set_default_namespace('monomerA')
+ core.print_out('\n')
+ p4util.banner('Monomer A HF (Dimer Basis)')
+ core.print_out('\n')
+ monomerA_wfn = scf_helper('RHF', molecule=monomerA, **kwargs)
+
+ if (ri == 'DF'):
+ core.IO.change_file_namespace(97, 'monomerA', 'monomerB')
+ core.IO.set_default_namespace('monomerB')
+ core.print_out('\n')
+ p4util.banner('Monomer B HF (Dimer Basis)')
+ core.print_out('\n')
+ monomerB_wfn = scf_helper('RHF', molecule=monomerB, **kwargs)
+ core.set_global_option('DF_INTS_IO', df_ints_io)
+
+ core.IO.set_default_namespace('monomerAm')
+ core.print_out('\n')
+ p4util.banner('Monomer A HF (Monomer Basis)')
+ core.print_out('\n')
+ monomerAm_wfn = scf_helper('RHF', molecule=monomerAm, **kwargs)
+
+ core.IO.set_default_namespace('monomerBm')
+ core.print_out('\n')
+ p4util.banner('Monomer B HF (Monomer Basis)')
+ core.print_out('\n')
+ monomerBm_wfn = scf_helper('RHF', molecule=monomerBm, **kwargs)
+
+ core.IO.set_default_namespace('dimer')
+ core.set_local_option('SAPT', 'E_CONVERGENCE', 10e-10)
+ core.set_local_option('SAPT', 'D_CONVERGENCE', 10e-10)
+ if name == 'sapt0-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT0')
+ elif name == 'sapt2-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2')
+ elif name == 'sapt2+-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+')
+ elif name == 'sapt2+(3)-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', False)
+ elif name == 'sapt2+3-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', True)
+ elif name == 'sapt2+(ccd)-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+')
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+ elif name == 'sapt2+(3)(ccd)-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', False)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+ elif name == 'sapt2+3(ccd)-ct':
+ core.set_local_option('SAPT', 'SAPT_LEVEL', 'SAPT2+3')
+ core.set_local_option('SAPT', 'DO_THIRD_ORDER', True)
+ core.set_local_option('SAPT', 'DO_CCD_DISP', True)
+
+ core.print_out('\n')
+
+ aux_basis = core.BasisSet.build(dimer_wfn.molecule(), "DF_BASIS_SAPT",
+ core.get_global_option("DF_BASIS_SAPT"),
+ "RIFIT", core.get_global_option("BASIS"))
+ dimer_wfn.set_basisset("DF_BASIS_SAPT", aux_basis)
+ if core.get_global_option("DF_BASIS_ELST") == "":
+ dimer_wfn.set_basisset("DF_BASIS_ELST", aux_basis)
+ else:
+ aux_basis = core.BasisSet.build(dimer_wfn.molecule(), "DF_BASIS_ELST",
+ core.get_global_option("DF_BASIS_ELST"),
+ "RIFIT", core.get_global_option("BASIS"))
+ dimer_wfn.set_basisset("DF_BASIS_ELST", aux_basis)
+
+
+ core.print_out('\n')
+ p4util.banner('SAPT Charge Transfer')
+ core.print_out('\n')
+
+ core.print_out('\n')
+ p4util.banner('Dimer Basis SAPT')
+ core.print_out('\n')
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERA, 'monomerA', 'dimer')
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERB, 'monomerB', 'dimer')
+ e_sapt = core.sapt(dimer_wfn, monomerA_wfn, monomerB_wfn)
+ CTd = core.variable('SAPT CT ENERGY')
+ dimer_wfn.set_module("sapt")
+
+ core.print_out('\n')
+ p4util.banner('Monomer Basis SAPT')
+ core.print_out('\n')
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERA, 'monomerAm', 'dimer')
+ core.IO.change_file_namespace(psif.PSIF_SAPT_MONOMERB, 'monomerBm', 'dimer')
+ e_sapt = core.sapt(dimer_wfn, monomerAm_wfn, monomerBm_wfn)
+ CTm = core.variable('SAPT CT ENERGY')
+ CT = CTd - CTm
+
+ units = (1000.0, constants.hartree2kcalmol, constants.hartree2kJmol)
+ core.print_out('\n\n')
+ core.print_out(' SAPT Charge Transfer Analysis\n')
+ core.print_out(' ------------------------------------------------------------------------------------------------\n')
+ core.print_out(' SAPT Induction (Dimer Basis) %12.4lf [mEh] %12.4lf [kcal/mol] %12.4lf [kJ/mol]\n' %
+ tuple(CTd * u for u in units))
+ core.print_out(' SAPT Induction (Monomer Basis)%12.4lf [mEh] %12.4lf [kcal/mol] %12.4lf [kJ/mol]\n' %
+ tuple(CTm * u for u in units))
+ core.print_out(' SAPT Charge Transfer %12.4lf [mEh] %12.4lf [kcal/mol] %12.4lf [kJ/mol]\n\n' %
+ tuple(CT * u for u in units))
+ core.set_variable("SAPT CT ENERGY", CT) # P::e SAPT
+
+ optstash.restore()
+ return dimer_wfn
+
+def run_fisapt(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ an F/ISAPT0 computation
+
+ """
+ optstash = p4util.OptionsState(['SCF_TYPE'])
+
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ # Get the molecule of interest
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ sapt_dimer = kwargs.pop('molecule', core.get_active_molecule())
+ else:
+ core.print_out('Warning! FISAPT argument "ref_wfn" is only able to use molecule information.')
+ sapt_dimer = ref_wfn.molecule()
+ sapt_dimer.update_geometry() # make sure since mol from wfn, kwarg, or P::e
+
+ # Shifting to C1 so we need to copy the active molecule
+ if sapt_dimer.schoenflies_symbol() != 'c1':
+ core.print_out(' FISAPT does not make use of molecular symmetry, further calculations in C1 point group.\n')
+ _ini_cart = getattr(sapt_dimer, "_initial_cartesian", None)
+ sapt_dimer = sapt_dimer.clone()
+ sapt_dimer.reset_point_group('c1')
+ sapt_dimer.fix_orientation(True)
+ sapt_dimer.fix_com(True)
+ sapt_dimer.update_geometry()
+ if _ini_cart:
+ sapt_dimer._initial_cartesian = _ini_cart
+
+ if core.get_option('SCF', 'REFERENCE') != 'RHF':
+ raise ValidationError('FISAPT requires requires \"reference rhf\".')
+
+ if ref_wfn is None:
+ core.timer_on("FISAPT: Dimer SCF")
+ ref_wfn = scf_helper('RHF', molecule=sapt_dimer, **kwargs)
+ core.timer_off("FISAPT: Dimer SCF")
+
+ core.print_out(" Constructing Basis Sets for FISAPT...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(),
+ "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT",
+ core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ sapt_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SAPT", core.get_global_option("DF_BASIS_SAPT"),
+ "RIFIT", core.get_global_option("BASIS"),
+ ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SAPT", sapt_basis)
+
+ minao = core.BasisSet.build(ref_wfn.molecule(), "BASIS", core.get_global_option("MINAO_BASIS"))
+ ref_wfn.set_basisset("MINAO", minao)
+
+ # Turn of dispersion for -d
+ if "-d" in name.lower():
+ core.set_local_option("FISAPT", "FISAPT_DO_FSAPT_DISP", False)
+
+ fisapt_wfn = core.FISAPT(ref_wfn)
+ from .sapt import fisapt_proc
+ fisapt_wfn.compute_energy(external_potentials=kwargs.get("external_potentials", None))
+
+ # Compute -D dispersion
+ if "-d" in name.lower():
+ proc_util.sapt_empirical_dispersion(name, ref_wfn)
+
+ optstash.restore()
+ return ref_wfn
+
+def run_mrcc(name, **kwargs):
+ """Function that prepares environment and input files
+ for a calculation calling Kallay's MRCC code.
+
+ """
+ from .proc_data import mrcc_methods
+
+ # level is a dictionary of settings to be passed to core.mrcc
+ try:
+ level = mrcc_methods[name.lower()]
+ except KeyError:
+ if name.lower() == "a-ccsd(t)":
+ level = mrcc_methods["ccsd(t)_l"]
+ else:
+ raise ValidationError(f"""MRCC method '{name}' invalid.""")
+
+ # Check to see if we really need to run the SCF code.
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs)
+ vscf = ref_wfn.variable('SCF TOTAL ENERGY')
+
+ # Fullname is the string we need to search for in iface
+ fullname = level['fullname']
+
+ # User can provide 'keep' to the method.
+ # When provided, do not delete the MRCC scratch directory.
+ keep = False
+ if 'keep' in kwargs:
+ keep = kwargs['keep']
+
+ # Save current directory location
+ current_directory = os.getcwd()
+
+ # Find environment by merging PSIPATH and PATH environment variables
+ lenv = {
+ 'PATH': ':'.join([os.path.abspath(x) for x in os.environ.get('PSIPATH', '').split(':') if x != '']) +
+ ':' + os.environ.get('PATH'),
+ 'LD_LIBRARY_PATH': os.environ.get('LD_LIBRARY_PATH')
+ }
+ # Filter out None values as subprocess will fault on them
+ lenv = {k: v for k, v in lenv.items() if v is not None}
+
+ # Need to move to the scratch directory, perferrably into a separate directory in that location
+ psi_io = core.IOManager.shared_object()
+ os.chdir(psi_io.get_default_path())
+
+ # Make new directory specifically for mrcc
+ mrcc_tmpdir = 'mrcc_' + str(os.getpid())
+ if 'path' in kwargs:
+ mrcc_tmpdir = kwargs['path']
+
+ # Check to see if directory already exists, if not, create.
+ if os.path.exists(mrcc_tmpdir) is False:
+ os.mkdir(mrcc_tmpdir)
+
+ # Move into the new directory
+ os.chdir(mrcc_tmpdir)
+
+ # Generate integrals and input file (dumps files to the current directory)
+ if core.get_option('FNOCC', 'NAT_ORBS'):
+ mints = core.MintsHelper(ref_wfn.basisset())
+ mints.set_print(1)
+ mints.integrals()
+
+ core.mrcc_generate_input(ref_wfn, level)
+ ref_wfn.set_module("mrcc")
+
+ # Load the fort.56 file
+ # and dump a copy into the outfile
+ core.print_out('\n===== Begin fort.56 input for MRCC ======\n')
+ core.print_out(open('fort.56', 'r').read())
+ core.print_out('===== End fort.56 input for MRCC ======\n')
+
+ # Modify the environment:
+ # PGI Fortan prints warning to screen if STOP is used
+ lenv['NO_STOP_MESSAGE'] = '1'
+
+ # Obtain the number of threads MRCC should use
+ lenv['OMP_NUM_THREADS'] = str(core.get_num_threads())
+
+ # If the user provided MRCC_OMP_NUM_THREADS set the environ to it
+ if core.has_option_changed('MRCC', 'MRCC_OMP_NUM_THREADS'):
+ lenv['OMP_NUM_THREADS'] = str(core.get_option('MRCC', 'MRCC_OMP_NUM_THREADS'))
+
+ # Call dmrcc, directing all screen output to the output file
+ external_exe = 'dmrcc'
+ try:
+ retcode = subprocess.Popen([external_exe], bufsize=0, stdout=subprocess.PIPE, env=lenv)
+ except OSError as e:
+ sys.stderr.write('Program %s not found in path or execution failed: %s\n' % (external_exe, e.strerror))
+ core.print_out('Program %s not found in path or execution failed: %s\n' % (external_exe, e.strerror))
+ message = ("Program %s not found in path or execution failed: %s\n" % (external_exe, e.strerror))
+ raise ValidationError(message)
+
+ c4out = ''
+ while True:
+ data = retcode.stdout.readline()
+ if not data:
+ break
+ core.print_out(data.decode('utf-8'))
+ c4out += data.decode('utf-8')
+
+ # Scan iface file and grab the file energy.
+ ene = 0.0
+ for line in open('iface'):
+ fields = line.split()
+ m = fields[1]
+ try:
+ ene = float(fields[5])
+ if m == "MP(2)":
+ m = "MP2"
+ elif m == "CCSD[T]":
+ m = "CCSD+T(CCSD)"
+ elif m == "CCSD(T)_L":
+ m = "A-CCSD(T)"
+ ref_wfn.set_variable(m + ' TOTAL ENERGY', ene)
+ ref_wfn.set_variable(m + ' CORRELATION ENERGY', ene - vscf)
+ except ValueError:
+ continue
+
+ # The last 'ene' in iface is the one the user requested.
+ ref_wfn.set_variable('CURRENT ENERGY', ene)
+ ref_wfn.set_variable('CURRENT CORRELATION ENERGY', ene - vscf)
+
+ for subm in ["MP2", "CCSD"]:
+ if ref_wfn.has_variable(f"{subm} TOTAL ENERGY") and core.get_option("SCF", "REFERENCE") in ["RHF", "UHF"]:
+ ref_wfn.set_variable(f"{subm} SINGLES ENERGY", 0.0)
+ ref_wfn.set_variable(f"{subm} DOUBLES ENERGY", ref_wfn.variable(f"{subm} CORRELATION ENERGY"))
+
+ if ref_wfn.has_variable("CCSD(T) CORRELATION ENERGY") and ref_wfn.has_variable("CCSD CORRELATION ENERGY"):
+ ref_wfn.set_variable("(T) CORRECTION ENERGY", ref_wfn.variable("CCSD(T) CORRELATION ENERGY") - ref_wfn.variable("CCSD CORRELATION ENERGY"))
+
+ if ref_wfn.has_variable("A-CCSD(T) CORRELATION ENERGY") and ref_wfn.has_variable("CCSD CORRELATION ENERGY"):
+ ref_wfn.set_variable("A-(T) CORRECTION ENERGY", ref_wfn.variable("A-CCSD(T) CORRELATION ENERGY") - ref_wfn.variable("CCSD CORRELATION ENERGY"))
+
+ # Load the iface file
+ iface = open('iface', 'r')
+ iface_contents = iface.read()
+
+ # Delete mrcc tempdir
+ os.chdir('..')
+ try:
+ # Delete unless we're told not to
+ if (keep is False and not('path' in kwargs)):
+ shutil.rmtree(mrcc_tmpdir)
+ except OSError as e:
+ print('Unable to remove MRCC temporary directory %s' % e, file=sys.stderr)
+ exit(1)
+
+ # Return to submission directory
+ os.chdir(current_directory)
+
+ # If we're told to keep the files or the user provided a path, do nothing.
+ if keep or ('path' in kwargs):
+ core.print_out('\nMRCC scratch files have been kept.\n')
+ core.print_out('They can be found in ' + mrcc_tmpdir)
+
+ # Dump iface contents to output
+ core.print_out('\n')
+ p4util.banner('Full results from MRCC')
+ core.print_out('\n')
+ core.print_out(iface_contents)
+
+ # Shove variables into global space
+ for k, v in ref_wfn.variables().items():
+ core.set_variable(k, v)
+
+ core.print_variables()
+ return ref_wfn
+
+
+def run_fnodfcc(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a [FNO-](DF|CD)-CCSD[(T)] computation.
+
+ >>> set cc_type df
+ >>> energy('fno-ccsd(t)')
+
+ """
+ kwargs = p4util.kwargs_lower(kwargs)
+ dtl = docs_table_link("dummy", "fnocc")
+
+ # stash user options
+ optstash = p4util.OptionsState(
+ ['FNOCC', 'COMPUTE_TRIPLES'],
+ ['FNOCC', 'DFCC'],
+ ['FNOCC', 'NAT_ORBS'],
+ ['FNOCC', 'RUN_CEPA'],
+ ['FNOCC', 'DF_BASIS_CC'],
+ ['SCF', 'DF_BASIS_SCF'],
+ ['SCF', 'DF_INTS_IO'])
+
+ def set_cholesky_from(type_val):
+ if type_val == 'CD':
+ core.set_local_option('FNOCC', 'DF_BASIS_CC', 'CHOLESKY')
+ # Alter default algorithm
+ if not core.has_global_option_changed('SCF_TYPE'):
+ optstash.add_option(['SCF_TYPE'])
+ core.set_global_option('SCF_TYPE', 'CD')
+ core.print_out(""" SCF Algorithm Type (re)set to CD.\n""")
+
+ elif type_val in ['DISK_DF', 'DF']:
+ if core.get_option('FNOCC', 'DF_BASIS_CC') == 'CHOLESKY':
+ core.set_local_option('FNOCC', 'DF_BASIS_CC', '')
+
+ proc_util.check_disk_df(name.upper(), optstash)
+ else:
+ raise ValidationError(f"Invalid type {type_val} for FNOCC energy through `run_fnodfcc`. See Capabilities Table at {dtl}")
+
+ director = {
+ # Note "nat_orbs" not set defensively False for non-"fno" calls
+ "ccsd": { "dfcc": True, "run_cepa": False, "compute_triples": False,},
+ "fno-ccsd": {"nat_orbs": True, "dfcc": True, "run_cepa": False, "compute_triples": False,},
+
+ "ccsd(t)": { "dfcc": True, "run_cepa": False, "compute_triples": True, },
+ "fno-ccsd(t)": {"nat_orbs": True, "dfcc": True, "run_cepa": False, "compute_triples": True, },
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for FNOCC energy")
+
+ # throw exception for open-shells
+ if (ref := core.get_option("SCF", "REFERENCE")) != "RHF":
+ raise ValidationError(f"Invalid reference type {ref} != RHF for FNOCC energy. See Capabilities Table at {dtl}.")
+
+ # throw exception for CONV (approximately). after defaulting logic, throw exception for SCF_TYPE CONV (approximately)
+ set_cholesky_from(method_algorithm_type(name).now)
+ if (scf_type := core.get_global_option("SCF_TYPE")) not in ["CD", "DISK_DF"]:
+ raise ValidationError(f"Invalid {scf_type=} for FNOCC energy through `run_fnodfcc`. See Capabilities Table at {dtl}")
+
+ for k, v in director[name].items():
+ core.set_local_option("FNOCC", k.upper(), v)
+
+ # save DF or CD ints generated by SCF for use in CC
+ core.set_local_option('SCF', 'DF_INTS_IO', 'SAVE')
+
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, use_c1=True, **kwargs) # C1 certified
+ else:
+ if ref_wfn.molecule().schoenflies_symbol() != 'c1':
+ raise ValidationError(""" FNOCC does not make use of molecular symmetry: """
+ """reference wavefunction must be C1.\n""")
+
+ core.print_out(" Constructing Basis Sets for FNOCC...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_CC",
+ core.get_global_option("DF_BASIS_CC"),
+ "RIFIT", core.get_global_option("BASIS"))
+ ref_wfn.set_basisset("DF_BASIS_CC", aux_basis)
+
+ if core.get_global_option("RELATIVISTIC") in ["X2C", "DKH"]:
+ rel_bas = core.BasisSet.build(ref_wfn.molecule(), "BASIS_RELATIVISTIC",
+ core.get_option("SCF", "BASIS_RELATIVISTIC"),
+ "DECON", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset('BASIS_RELATIVISTIC',rel_bas)
+
+ fnocc_wfn = core.fnocc(ref_wfn)
+
+ # Shove variables into global space
+ for k, v in fnocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return fnocc_wfn
+
+
+def run_fnocc(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a QCISD(T), CCSD(T), MP2.5, MP3, and MP4 computation.
+
+ >>> energy('fno-ccsd(t)')
+
+ """
+ kwargs = p4util.kwargs_lower(kwargs)
+ dtl = docs_table_link("dummy", "fnocc")
+
+ # stash user options:
+ optstash = p4util.OptionsState(
+ ['TRANSQT2', 'WFN'],
+ ['FNOCC', 'RUN_MP2'],
+ ['FNOCC', 'RUN_MP3'],
+ ['FNOCC', 'RUN_MP4'],
+ ['FNOCC', 'RUN_CCSD'],
+ ['FNOCC', 'COMPUTE_TRIPLES'],
+ ['FNOCC', 'COMPUTE_MP4_TRIPLES'],
+ ['FNOCC', 'DFCC'],
+ ['FNOCC', 'RUN_CEPA'],
+ ['FNOCC', 'USE_DF_INTS'],
+ ['FNOCC', 'NAT_ORBS'])
+
+ # AED reinforces default
+ core.set_local_option('FNOCC', 'USE_DF_INTS', False)
+
+ if name in ["mp3", "fno-mp3"] and not core.has_global_option_changed("MP_TYPE"):
+ core.print_out(f" Information: {name.upper()} default algorithm changed to DF in August 2020. Use `set mp_type conv` for previous behavior.\n")
+
+ director = {
+ # Note "nat_orbs" not set defensively False for non-"fno" calls
+ "mp2": { "dfcc": False, "run_cepa": False, "run_mp2": True, },
+
+ "mp3": { "dfcc": False, "run_cepa": False, "run_mp3": True, },
+ "fno-mp3": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "run_mp3": True, },
+
+ "mp4(sdq)": { "dfcc": False, "run_cepa": False, "compute_triples": False, "run_mp4": True, "compute_mp4_triples": False,},
+ "fno-mp4(sdq)": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "compute_triples": False, "run_mp4": True, "compute_mp4_triples": False,},
+
+ "mp4": { "dfcc": False, "run_cepa": False, "compute_triples": True, "run_mp4": True, "compute_mp4_triples": True, },
+ "fno-mp4": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "compute_triples": True, "run_mp4": True, "compute_mp4_triples": True, },
+
+ "qcisd": { "dfcc": False, "run_cepa": False, "run_ccsd": False, "compute_triples": False, },
+ "fno-qcisd": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "run_ccsd": False, "compute_triples": False, },
+
+ "qcisd(t)": { "dfcc": False, "run_cepa": False, "run_ccsd": False, "compute_triples": True, },
+ "fno-qcisd(t)": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "run_ccsd": False, "compute_triples": True, },
+
+ "ccsd": { "dfcc": False, "run_cepa": False, "run_ccsd": True, "compute_triples": False, },
+ "fno-ccsd": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "run_ccsd": True, "compute_triples": False, },
+
+ "ccsd(t)": { "dfcc": False, "run_cepa": False, "run_ccsd": True, "compute_triples": True, },
+ "fno-ccsd(t)": {"nat_orbs": True, "dfcc": False, "run_cepa": False, "run_ccsd": True, "compute_triples": True, },
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for FNOCC energy")
+
+ # throw exception for open-shells
+ if (ref := core.get_option("SCF", "REFERENCE")) != "RHF":
+ raise ValidationError(f"Invalid reference type {ref} != RHF for FNOCC energy. See Capabilities Table at {dtl}")
+
+ # throw exception for DF/CD. most of these pre-trapped by select_* functions but some escape, incl. mp4(sdq) and qcisd variants
+ if (corl_type := method_algorithm_type(name).now) != "CONV":
+ raise ValidationError(f"Invalid type {corl_type} for FNOCC energy through `run_fnocc`. See Capabilities Table at {dtl}")
+
+ for k, v in director[name].items():
+ core.set_local_option("FNOCC", k.upper(), v)
+
+ # Bypass the scf call if a reference wavefunction is given
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ if not core.get_option('FNOCC', 'USE_DF_INTS'):
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+ else:
+ core.print_out(" Constructing Basis Sets for FNOCC...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ if core.get_global_option("RELATIVISTIC") in ["X2C", "DKH"]:
+ rel_bas = core.BasisSet.build(ref_wfn.molecule(), "BASIS_RELATIVISTIC",
+ core.get_option("SCF", "BASIS_RELATIVISTIC"),
+ "DECON", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset('BASIS_RELATIVISTIC',rel_bas)
+
+ fnocc_wfn = core.fnocc(ref_wfn)
+
+ # set current correlation energy and total energy. only need to treat mpn here.
+ if (lbl := name.replace("fno-", "")) in ["mp3", "mp4(sdq)", "mp4"]:
+ fnocc_wfn.set_variable("CURRENT ENERGY", fnocc_wfn.variable(f"{lbl.upper()} TOTAL ENERGY"))
+ fnocc_wfn.set_variable("CURRENT CORRELATION ENERGY", fnocc_wfn.variable(f"{lbl.upper()} CORRELATION ENERGY"))
+ if lbl == "mp4":
+ fnocc_wfn.set_variable("MP4 CORRECTION ENERGY", fnocc_wfn.variable("MP4 CORRELATION ENERGY") - fnocc_wfn.variable("MP3 CORRELATION ENERGY"))
+
+ # Shove variables into global space
+ for k, v in fnocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return fnocc_wfn
+
+
+def run_cepa(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ a cepa-like calculation.
+
+ >>> energy('cepa(1)')
+
+ """
+ kwargs = p4util.kwargs_lower(kwargs)
+ dtl = docs_table_link("dummy", "fnocc")
+
+ # save user options
+ optstash = p4util.OptionsState(
+ ['TRANSQT2', 'WFN'],
+ ['FNOCC', 'NAT_ORBS'],
+ ['FNOCC', 'DFCC'],
+ ['FNOCC', 'RUN_CEPA'],
+ ['FNOCC', 'USE_DF_INTS'],
+ ['FNOCC', 'CEPA_LEVEL'],
+ ['FNOCC', 'CEPA_NO_SINGLES'])
+
+ # AED reinforces default
+ core.set_local_option('FNOCC', 'USE_DF_INTS', False)
+
+ director = {
+ # Note "nat_orbs" not set defensively False for non-"fno" calls
+ "lccd": { "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": True, },
+ "fno-lccd": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": True, },
+
+ "lccsd": { "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": False,},
+ "fno-lccsd": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": False,},
+ "cepa(0)": { "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": False,},
+ "fno-cepa(0)": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cepa(0)", "cepa_no_singles": False,},
+
+ "cepa(1)": { "dfcc": False, "run_cepa": True, "cepa_level": "cepa(1)", "cepa_no_singles": False,},
+ "fno-cepa(1)": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cepa(1)", "cepa_no_singles": False,},
+
+ "cepa(3)": { "dfcc": False, "run_cepa": True, "cepa_level": "cepa(3)", "cepa_no_singles": False,},
+ "fno-cepa(3)": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cepa(3)", "cepa_no_singles": False,},
+
+ "acpf": { "dfcc": False, "run_cepa": True, "cepa_level": "acpf", "cepa_no_singles": False,},
+ "fno-acpf": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "acpf", "cepa_no_singles": False,},
+
+ "aqcc": { "dfcc": False, "run_cepa": True, "cepa_level": "aqcc", "cepa_no_singles": False,},
+ "fno-aqcc": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "aqcc", "cepa_no_singles": False,},
+
+ "cisd": { "dfcc": False, "run_cepa": True, "cepa_level": "cisd", "cepa_no_singles": False,},
+ "fno-cisd": {"nat_orbs": True, "dfcc": False, "run_cepa": True, "cepa_level": "cisd", "cepa_no_singles": False,},
+ }
+
+ if name not in director:
+ raise ValidationError(f"Invalid method {name} for FNOCC energy")
+
+ # throw exception for open-shells
+ if (ref := core.get_option("SCF", "REFERENCE")) != "RHF":
+ raise ValidationError(f"Invalid reference type {ref} != RHF for FNOCC energy. See Capabilities Table at {dtl}")
+
+ # throw exception for DF/CD. some of these pre-trapped by select_* functions but others escape, incl. cepa variants
+ if (corl_type := method_algorithm_type(name).now) != "CONV":
+ raise ValidationError(f"Invalid type {corl_type} for FNOCC energy through `run_cepa`. See Capabilities Table at {dtl}")
+
+ for k, v in director[name].items():
+ core.set_local_option("FNOCC", k.upper(), v)
+
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+ ref_wfn = scf_helper(name, **kwargs) # C1 certified
+
+ if not core.get_option('FNOCC', 'USE_DF_INTS'):
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+ else:
+ core.print_out(" Constructing Basis Sets for FISAPT...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ fnocc_wfn = core.fnocc(ref_wfn)
+
+ # one-electron properties
+ cepa_level = director[name]["cepa_level"]
+ if core.get_option('FNOCC', 'DIPMOM'):
+ if cepa_level in ['cepa(1)', 'cepa(3)']:
+ core.print_out("""\n Error: one-electron properties not implemented for %s\n\n""" % name)
+ elif core.get_option('FNOCC', 'NAT_ORBS'):
+ core.print_out("""\n Error: one-electron properties not implemented for %s\n\n""" % name)
+ else:
+ p4util.oeprop(fnocc_wfn, 'DIPOLE', 'QUADRUPOLE', 'MULLIKEN_CHARGES', 'NO_OCCUPATIONS', title=cepa_level.upper())
+
+ # Shove variables into global space
+ for k, v in fnocc_wfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return fnocc_wfn
+
+
+def run_detcas(name, **kwargs):
+ """Function encoding sequence of PSI module calls for
+ determinant-based multireference wavefuncations,
+ namely CASSCF and RASSCF.
+ """
+ optstash = p4util.OptionsState(
+ ['DETCI', 'WFN'],
+ ['SCF_TYPE'],
+ ['ONEPDM'],
+ ['OPDM_RELAX']
+ )
+
+ user_ref = core.get_option('DETCI', 'REFERENCE')
+ if user_ref not in ['RHF', 'ROHF']:
+ raise ValidationError('Reference %s for DETCI is not available.' % user_ref)
+
+ if name == 'rasscf':
+ core.set_local_option('DETCI', 'WFN', 'RASSCF')
+ elif name == 'casscf':
+ core.set_local_option('DETCI', 'WFN', 'CASSCF')
+ else:
+ raise ValidationError("Run DETCAS: Name %s not understood" % name)
+
+ ref_wfn = kwargs.get('ref_wfn', None)
+ if ref_wfn is None:
+
+ ref_optstash = p4util.OptionsState(
+ ['SCF_TYPE'],
+ ['DF_BASIS_SCF'],
+ ['DF_BASIS_MP2'],
+ ['ONEPDM'],
+ ['OPDM_RELAX']
+ )
+
+ # No real reason to do a conventional guess
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ # If RHF get MP2 NO's
+ # Why doesnt this work for conv?
+ if (('DF' in core.get_global_option('SCF_TYPE')) and (user_ref == 'RHF') and
+ (core.get_option('DETCI', 'MCSCF_TYPE') in ['DF', 'AO']) and
+ (core.get_option("DETCI", "MCSCF_GUESS") == "MP2")):
+ core.set_global_option('ONEPDM', True)
+ core.set_global_option('OPDM_RELAX', False)
+ ref_wfn = run_dfmp2_gradient(name, **kwargs)
+ else:
+ ref_wfn = scf_helper(name, **kwargs)
+
+ # Ensure IWL files have been written
+ if (core.get_option('DETCI', 'MCSCF_TYPE') == 'CONV'):
+ mints = core.MintsHelper(ref_wfn.basisset())
+ mints.set_print(1)
+ mints.integrals()
+
+ ref_optstash.restore()
+
+ # The DF case
+ if core.get_option('DETCI', 'MCSCF_TYPE') == 'DF':
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DF')
+
+ core.print_out(" Constructing Basis Sets for MCSCF...\n\n")
+ scf_aux_basis = core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_SCF",
+ core.get_option("SCF", "DF_BASIS_SCF"),
+ "JKFIT", core.get_global_option('BASIS'),
+ puream=ref_wfn.basisset().has_puream())
+ ref_wfn.set_basisset("DF_BASIS_SCF", scf_aux_basis)
+
+ # The AO case
+ elif core.get_option('DETCI', 'MCSCF_TYPE') == 'AO':
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'DIRECT')
+
+ # The conventional case
+ elif core.get_option('DETCI', 'MCSCF_TYPE') == 'CONV':
+ if not core.has_global_option_changed('SCF_TYPE'):
+ core.set_global_option('SCF_TYPE', 'PK')
+
+ # Ensure IWL files have been written
+ proc_util.check_iwl_file_from_scf_type(core.get_global_option('SCF_TYPE'), ref_wfn)
+ else:
+ raise ValidationError("Run DETCAS: MCSCF_TYPE %s not understood." % str(core.get_option('DETCI', 'MCSCF_TYPE')))
+
+
+ # Second-order SCF requires non-symmetric density matrix support
+ if core.get_option('DETCI', 'MCSCF_ALGORITHM') in ['AH', 'OS']:
+ proc_util.check_non_symmetric_jk_density("Second-order MCSCF")
+
+ ciwfn = mcscf.mcscf_solver(ref_wfn)
+
+ # We always would like to print a little dipole information
+ oeprop = core.OEProp(ciwfn)
+ oeprop.set_title(name.upper())
+ oeprop.add("DIPOLE")
+ oeprop.compute()
+ ciwfn.oeprop = oeprop
+ core.set_variable("CURRENT DIPOLE", core.variable(name.upper() + " DIPOLE"))
+
+ # Shove variables into global space
+ for k, v in ciwfn.variables().items():
+ core.set_variable(k, v)
+
+ optstash.restore()
+ return ciwfn
+
+def run_efp(name, **kwargs):
+ """Function encoding sequence of module calls for a pure EFP
+ computation (ignore any QM atoms).
+
+ """
+
+ efp_molecule = kwargs.get('molecule', core.get_active_molecule())
+ try:
+ efpobj = efp_molecule.EFP
+ except AttributeError:
+ raise ValidationError("""Method 'efp' not available without EFP fragments in molecule""")
+
+ # print efp geom in [A]
+ core.print_out(efpobj.banner())
+ core.print_out(efpobj.geometry_summary(units_to_bohr=constants.bohr2angstroms))
+
+ # set options
+ # * 'chtr', 'qm_exch', 'qm_disp', 'qm_chtr' may be enabled in a future libefp release
+ efpopts = {}
+ for opt in ['elst', 'exch', 'ind', 'disp',
+ 'elst_damping', 'ind_damping', 'disp_damping']:
+ psiopt = 'EFP_' + opt.upper()
+ if core.has_option_changed('EFP', psiopt):
+ efpopts[opt] = core.get_option('EFP', psiopt)
+ efpopts['qm_elst'] = False
+ efpopts['qm_ind'] = False
+ efpobj.set_opts(efpopts, label='psi', append='psi')
+ do_gradient = core.get_option('EFP', 'DERTYPE') == 'FIRST'
+
+ # compute and report
+ efpobj.compute(do_gradient=do_gradient)
+ core.print_out(efpobj.energy_summary(label='psi'))
+
+ ene = efpobj.get_energy(label='psi')
+ core.set_variable('EFP ELST ENERGY', ene['electrostatic'] + ene['charge_penetration'] + ene['electrostatic_point_charges'])
+ core.set_variable('EFP IND ENERGY', ene['polarization'])
+ core.set_variable('EFP DISP ENERGY', ene['dispersion'])
+ core.set_variable('EFP EXCH ENERGY', ene['exchange_repulsion'])
+ core.set_variable('EFP TOTAL ENERGY', ene['total'])
+ core.set_variable('CURRENT ENERGY', ene['total'])
+
+ if do_gradient:
+ core.print_out(efpobj.gradient_summary())
+
+ torq = efpobj.get_gradient()
+ torq = core.Matrix.from_array(np.asarray(torq).reshape(-1, 6))
+ core.set_variable("EFP TORQUE", torq) # P::e EFP
+
+ return ene['total']
+