#
# @BEGIN LICENSE
#
# Psi4: an open-source quantum chemistry software package
#
# Copyright (c) 2007-2018 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
#
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
import re
import sys
import math
import numpy as np
from psi4 import core
from psi4.driver import qcdb
from psi4.driver import p4util
from psi4.driver import driver_util
from psi4.driver import constants
from psi4.driver.p4util.exceptions import *
from psi4.driver.procrouting.interface_cfour import cfour_psivar_list
zeta_values = ['d', 't', 'q', '5', '6', '7', '8']
zeta_val2sym = {k + 2: v for k, v in zip(range(7), zeta_values)}
zeta_sym2val = {v: k for k, v in zeta_val2sym.items()}
def _expand_bracketed_basis(basisstring, molecule=None):
r"""Function to transform and validate basis series specification
*basisstring* for cbs(). A basis set with no paired square brackets is
passed through with zeta level 0 (e.g., '6-31+G(d,p)' is returned as
[6-31+G(d,p)] and [0]). A basis set with square brackets is checked
for sensible sequence and Dunning-ness and returned as separate basis
sets (e.g., 'cc-pV[Q5]Z' is returned as [cc-pVQZ, cc-pV5Z] and [4,
5]). This function checks that the basis is valid by trying to build
the qcdb.BasisSet object for *molecule* or for H2 if None. Allows
out-of-order zeta specification (e.g., [qtd]) and numeral for number
(e.g., [23]) but not skipped zetas (e.g., [dq]) or zetas outside [2,
8] or non-Dunning or non-Ahlrichs or non-Jensen sets or
non-findable .gbs sets.
"""
BSET = []
ZSET = []
legit_compound_basis = re.compile(r'^(?P<pre>.*cc-.*|def2-|.*pcs+eg-)\[(?P<zeta>[dtq2345678,s1]*)\](?P<post>.*z.*|)$', re.IGNORECASE)
pc_basis = re.compile(r'.*pcs+eg-$', re.IGNORECASE)
def2_basis = re.compile(r'def2-', re.IGNORECASE)
if legit_compound_basis.match(basisstring):
basisname = legit_compound_basis.match(basisstring)
# handle def2-svp* basis sets as double-zeta
if def2_basis.match(basisname.group('pre')):
bn_gz = basisname.group('zeta').replace("s","d")
# handle pc-n basis set polarisation -> zeta conversion
elif pc_basis.match(basisname.group('pre')):
bn_gz = basisname.group('zeta').replace("4","5").replace("3","4").replace("2","3").replace("1","2")
else:
bn_gz = basisname.group('zeta')
# filter out commas and be forgiving of e.g., t5q or 3q
zetas = [z for z in zeta_values if (z in bn_gz or str(zeta_values.index(z) + 2) in bn_gz)]
for b in zetas:
if ZSET and (int(ZSET[len(ZSET) - 1]) - zeta_values.index(b)) != 1:
raise ValidationError("""Basis set '%s' has skipped zeta level '%s'.""" % (basisstring, zeta_val2sym[zeta_sym2val[b] - 1]))
# reassemble def2-svp* properly instead of def2-dzvp*
if def2_basis.match(basisname.group('pre')) and b == "d":
BSET.append(basisname.group('pre') + "s" + basisname.group('post')[1:])
# reassemble pc-n basis sets properly
elif pc_basis.match(basisname.group('pre')):
BSET.append(basisname.group('pre') + "{0:d}".format(zeta_sym2val[b] - 1))
else:
BSET.append(basisname.group('pre') + b + basisname.group('post'))
ZSET.append(zeta_values.index(b) + 2)
elif re.match(r'.*\[.*\].*$', basisstring, flags=re.IGNORECASE):
raise ValidationError("""Basis series '%s' invalid. Specify a basis series matching"""
""" '*cc-*[dtq2345678,]*z*'. or 'def2-[sdtq]zvp*' or '*pcs[s]eg-[1234]'""" % (basisstring))
else:
BSET.append(basisstring)
ZSET.append(0)
if molecule is None:
molecule = """\nH\nH 1 1.00\n"""
for basis in BSET:
try:
qcdb.BasisSet.pyconstruct(molecule, "BASIS", basis)
except qcdb.BasisSetNotFound:
e=sys.exc_info()[1]
raise ValidationError("""Basis set '%s' not available for molecule.""" % (basis))
return (BSET, ZSET)
def _contract_bracketed_basis(basisarray):
r"""Function to reform a bracketed basis set string from a sequential series
of basis sets *basisarray* (e.g, form 'cc-pv[q5]z' from array [cc-pvqz, cc-pv5z]).
Used to print a nicely formatted basis set string in the results table.
"""
if len(basisarray) == 1:
return basisarray[0]
else:
zetaindx = [i for i in range(len(basisarray[0])) if basisarray[0][i] != basisarray[1][i]][0]
ZSET = [bas[zetaindx] for bas in basisarray]
pre = basisarray[1][:zetaindx]
post = basisarray[1][zetaindx + 1:]
basisstring = pre + '[' + ''.join(ZSET) + ']' + post
return basisstring
[docs]def xtpl_highest_1(functionname, zHI, valueHI, verbose=True, **kwargs):
r"""Scheme for total or correlation energies with a single basis or the highest
zeta-level among an array of bases. Used by :py:func:`~psi4.cbs`.
.. math:: E_{total}^X = E_{total}^X
"""
if isinstance(valueHI, float):
if verbose:
# Output string with extrapolation parameters
cbsscheme = ''
cbsscheme += """\n ==> %s <==\n\n""" % (functionname.upper())
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
core.print_out(cbsscheme)
return valueHI
elif isinstance(valueHI, (core.Matrix, core.Vector)):
if verbose > 2:
core.print_out(""" HI-zeta (%s) Total Energy:\n""" % (str(zHI)))
valueHI.print_out()
return valueHI
[docs]def scf_xtpl_helgaker_2(functionname, zLO, valueLO, zHI, valueHI, verbose=True, alpha=None):
r"""Extrapolation scheme using exponential form for reference energies with two adjacent zeta-level bases.
Used by :py:func:`~psi4.cbs`.
Halkier, Helgaker, Jorgensen, Klopper, & Olsen, Chem. Phys. Lett. 302 (1999) 437-446,
DOI: 10.1016/S0009-2614(99)00179-7
.. math:: E_{total}^X = E_{total}^{\infty} + \beta e^{-\alpha X}, \alpha = 1.63
"""
if type(valueLO) != type(valueHI):
raise ValidationError("scf_xtpl_helgaker_2: Inputs must be of the same datatype! (%s, %s)"
% (type(valueLO), type(valueHI)))
if alpha is None:
alpha = 1.63
beta_division = 1 / (math.exp(-1 * alpha * zLO) * (math.exp(-1 * alpha) - 1))
beta_mult = math.exp(-1 * alpha * zHI)
if isinstance(valueLO, float):
beta = (valueHI - valueLO) * beta_division
value = valueHI - beta * beta_mult
if verbose:
# Output string with extrapolation parameters
cbsscheme = ''
cbsscheme += """\n ==> Helgaker 2-point exponential SCF extrapolation for method: %s <==\n\n""" % (functionname.upper())
cbsscheme += """ LO-zeta (%s) Energy: % 16.12f\n""" % (str(zLO), valueLO)
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
cbsscheme += """ Alpha (exponent) Value: % 16.12f\n""" % (alpha)
cbsscheme += """ Beta (coefficient) Value: % 16.12f\n\n""" % (beta)
name_str = "%s/(%s,%s)" % (functionname.upper(), zeta_val2sym[zLO].upper(), zeta_val2sym[zHI].upper())
cbsscheme += """ @Extrapolated """
cbsscheme += name_str + ':'
cbsscheme += " " * (18 - len(name_str))
cbsscheme += """% 16.12f\n\n""" % value
core.print_out(cbsscheme)
return value
elif isinstance(valueLO, (core.Matrix, core.Vector)):
beta = valueHI.clone()
beta.name = 'Helgaker SCF (%s, %s) beta' % (zLO, zHI)
beta.subtract(valueLO)
beta.scale(beta_division)
beta.scale(beta_mult)
value = valueHI.clone()
value.subtract(beta)
value.name = 'Helgaker SCF (%s, %s) data' % (zLO, zHI)
if verbose > 2:
core.print_out("""\n ==> Helgaker 2-point exponential SCF extrapolation for method: %s <==\n\n""" % (functionname.upper()))
core.print_out(""" LO-zeta (%s)""" % str(zLO))
core.print_out(""" LO-zeta Data""")
valueLO.print_out()
core.print_out(""" HI-zeta (%s)""" % str(zHI))
core.print_out(""" HI-zeta Data""")
valueHI.print_out()
core.print_out(""" Extrapolated Data:\n""")
value.print_out()
core.print_out(""" Alpha (exponent) Value: %16.8f\n""" % (alpha))
core.print_out(""" Beta Data:\n""")
beta.print_out()
return value
else:
raise ValidationError("scf_xtpl_helgaker_2: datatype is not recognized '%s'." % type(valueLO))
def scf_xtpl_truhlar_2(functionname, zLO, valueLO, zHI, valueHI, verbose=True, alpha=None):
r"""Extrapolation scheme using power form for reference energies with two adjacent zeta-level bases.
Used by :py:func:`~psi4.cbs`.
Truhlar, Chem. Phys. Lett. 294 (1998) 45-48, DOI: 10.1016/S0009-2614(98)00866-5
.. math:: E_{total}^X = E_{total}^{\infty} + \beta X^{-\alpha}, \alpha = 3.4
"""
if type(valueLO) != type(valueHI):
raise ValidationError("scf_xtpl_truhlar_2: Inputs must be of the same datatype! (%s, %s)"
% (type(valueLO), type(valueHI)))
if alpha is None:
alpha = 3.40
beta_division = 1 / (zHI ** (-1 * alpha) - zLO ** (-1 * alpha))
beta_mult = zHI ** (-1 * alpha)
if isinstance(valueLO, float):
beta = (valueHI - valueLO) * beta_division
value = valueHI - beta * beta_mult
if verbose:
# Output string with extrapolation parameters
cbsscheme = ''
cbsscheme += """\n ==> Truhlar 2-point power form SCF extrapolation for method: %s <==\n\n""" % (functionname.upper())
cbsscheme += """ LO-zeta (%s) Energy: % 16.12f\n""" % (str(zLO), valueLO)
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
cbsscheme += """ Alpha (exponent) Value: % 16.12f\n""" % (alpha)
cbsscheme += """ Beta (coefficient) Value: % 16.12f\n\n""" % (beta)
name_str = "%s/(%s,%s)" % (functionname.upper(), zeta_val2sym[zLO].upper(), zeta_val2sym[zHI].upper())
cbsscheme += """ @Extrapolated """
cbsscheme += name_str + ':'
cbsscheme += " " * (18 - len(name_str))
cbsscheme += """% 16.12f\n\n""" % value
core.print_out(cbsscheme)
return value
elif isinstance(valueLO, (core.Matrix, core.Vector)):
beta = valueHI.clone()
beta.name = 'Truhlar SCF (%s, %s) beta' % (zLO, zHI)
beta.subtract(valueLO)
beta.scale(beta_division)
beta.scale(beta_mult)
value = valueHI.clone()
value.subtract(beta)
value.name = 'Truhlar SCF (%s, %s) data' % (zLO, zHI)
if verbose > 2:
core.print_out("""\n ==> Truhlar 2-point power from SCF extrapolation for method: %s <==\n\n""" % (functionname.upper()))
core.print_out(""" LO-zeta (%s)""" % str(zLO))
core.print_out(""" LO-zeta Data""")
valueLO.print_out()
core.print_out(""" HI-zeta (%s)""" % str(zHI))
core.print_out(""" HI-zeta Data""")
valueHI.print_out()
core.print_out(""" Extrapolated Data:\n""")
value.print_out()
core.print_out(""" Alpha (exponent) Value: %16.8f\n""" % (alpha))
core.print_out(""" Beta Data:\n""")
beta.print_out()
return value
else:
raise ValidationError("scf_xtpl_truhlar_2: datatype is not recognized '%s'." % type(valueLO))
def scf_xtpl_karton_2(functionname, zLO, valueLO, zHI, valueHI, verbose=True, alpha=None):
r"""Extrapolation scheme using root-power form for reference energies with two adjacent zeta-level bases.
Used by :py:func:`~psi4.cbs`.
Karton, Martin, Theor. Chem. Acc. 115 (2006) 330-333, DOI: 10.1007/s00214-005-0028-6
.. math:: E_{total}^X = E_{total}^{\infty} + \beta e^{-\alpha\sqrt{X}}, \alpha = 6.3
"""
if type(valueLO) != type(valueHI):
raise ValidationError("scf_xtpl_karton_2: Inputs must be of the same datatype! (%s, %s)"
% (type(valueLO), type(valueHI)))
if alpha is None:
alpha = 6.30
beta_division = 1 / (math.exp(-1 * alpha) * (math.exp(math.sqrt(zHI)) - math.exp(math.sqrt(zLO))))
beta_mult = math.exp(-1 * alpha * math.sqrt(zHI))
if isinstance(valueLO, float):
beta = (valueHI - valueLO) * beta_division
value = valueHI - beta * beta_mult
if verbose:
# Output string with extrapolation parameters
cbsscheme = ''
cbsscheme += """\n ==> Karton 2-point power form SCF extrapolation for method: %s <==\n\n""" % (functionname.upper())
cbsscheme += """ LO-zeta (%s) Energy: % 16.12f\n""" % (str(zLO), valueLO)
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
cbsscheme += """ Alpha (exponent) Value: % 16.12f\n""" % (alpha)
cbsscheme += """ Beta (coefficient) Value: % 16.12f\n\n""" % (beta)
name_str = "%s/(%s,%s)" % (functionname.upper(), zeta_val2sym[zLO].upper(), zeta_val2sym[zHI].upper())
cbsscheme += """ @Extrapolated """
cbsscheme += name_str + ':'
cbsscheme += " " * (18 - len(name_str))
cbsscheme += """% 16.12f\n\n""" % value
core.print_out(cbsscheme)
return value
elif isinstance(valueLO, (core.Matrix, core.Vector)):
beta = valueHI.clone()
beta.name = 'Karton SCF (%s, %s) beta' % (zLO, zHI)
beta.subtract(valueLO)
beta.scale(beta_division)
beta.scale(beta_mult)
value = valueHI.clone()
value.subtract(beta)
value.name = 'Karton SCF (%s, %s) data' % (zLO, zHI)
if verbose > 2:
core.print_out("""\n ==> Karton 2-point power from SCF extrapolation for method: %s <==\n\n""" % (functionname.upper()))
core.print_out(""" LO-zeta (%s)""" % str(zLO))
core.print_out(""" LO-zeta Data""")
valueLO.print_out()
core.print_out(""" HI-zeta (%s)""" % str(zHI))
core.print_out(""" HI-zeta Data""")
valueHI.print_out()
core.print_out(""" Extrapolated Data:\n""")
value.print_out()
core.print_out(""" Alpha (exponent) Value: %16.8f\n""" % (alpha))
core.print_out(""" Beta Data:\n""")
beta.print_out()
return value
else:
raise ValidationError("scf_xtpl_Karton_2: datatype is not recognized '%s'." % type(valueLO))
[docs]def scf_xtpl_helgaker_3(functionname, zLO, valueLO, zMD, valueMD, zHI, valueHI, verbose=True, alpha=None):
r"""Extrapolation scheme for reference energies with three adjacent zeta-level bases.
Used by :py:func:`~psi4.cbs`.
Halkier, Helgaker, Jorgensen, Klopper, & Olsen, Chem. Phys. Lett. 302 (1999) 437-446,
DOI: 10.1016/S0009-2614(99)00179-7
.. math:: E_{total}^X = E_{total}^{\infty} + \beta e^{-\alpha X}, \alpha = 3.0
"""
if (type(valueLO) != type(valueMD)) or (type(valueMD) != type(valueHI)):
raise ValidationError("scf_xtpl_helgaker_3: Inputs must be of the same datatype! (%s, %s, %s)"
% (type(valueLO), type(valueMD), type(valueHI)))
if isinstance(valueLO, float):
ratio = (valueHI - valueMD) / (valueMD - valueLO)
alpha = -1 * math.log(ratio)
beta = (valueHI - valueMD) / (math.exp(-1 * alpha * zMD) * (ratio - 1))
value = valueHI - beta * math.exp(-1 * alpha * zHI)
if verbose:
# Output string with extrapolation parameters
cbsscheme = ''
cbsscheme += """\n ==> Helgaker 3-point SCF extrapolation for method: %s <==\n\n""" % (functionname.upper())
cbsscheme += """ LO-zeta (%s) Energy: % 16.12f\n""" % (str(zLO), valueLO)
cbsscheme += """ MD-zeta (%s) Energy: % 16.12f\n""" % (str(zMD), valueMD)
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
cbsscheme += """ Alpha (exponent) Value: % 16.12f\n""" % (alpha)
cbsscheme += """ Beta (coefficient) Value: % 16.12f\n\n""" % (beta)
name_str = "%s/(%s,%s,%s)" % (functionname.upper(), zeta_val2sym[zLO].upper(), zeta_val2sym[zMD].upper(),
zeta_val2sym[zHI].upper())
cbsscheme += """ @Extrapolated """
cbsscheme += name_str + ':'
cbsscheme += " " * (18 - len(name_str))
cbsscheme += """% 16.12f\n\n""" % value
core.print_out(cbsscheme)
return value
elif isinstance(valueLO, (core.Matrix, core.Vector)):
valueLO = np.array(valueLO)
valueMD = np.array(valueMD)
valueHI = np.array(valueHI)
nonzero_mask = np.abs(valueHI) > 1.e-14
top = (valueHI - valueMD)[nonzero_mask]
bot = (valueMD - valueLO)[nonzero_mask]
ratio = top / bot
alpha = -1 * np.log(np.abs(ratio))
beta = top / (np.exp(-1 * alpha * zMD) * (ratio - 1))
np_value = valueHI.copy()
np_value[nonzero_mask] -= beta * np.exp(-1 * alpha * zHI)
np_value[~nonzero_mask] = 0.0
# Build and set from numpy routines
value = core.Matrix(*valueHI.shape)
value_view = np.asarray(value)
value_view[:] = np_value
return value
else:
raise ValidationError("scf_xtpl_helgaker_3: datatype is not recognized '%s'." % type(valueLO))
#def corl_xtpl_helgaker_2(functionname, valueSCF, zLO, valueLO, zHI, valueHI, verbose=True):
[docs]def corl_xtpl_helgaker_2(functionname, zLO, valueLO, zHI, valueHI, verbose=True, alpha=None):
r"""Extrapolation scheme for correlation energies with two adjacent zeta-level bases.
Used by :py:func:`~psi4.cbs`.
Halkier, Helgaker, Jorgensen, Klopper, Koch, Olsen, & Wilson, Chem. Phys. Lett. 286 (1998) 243-252,
DOI: 10.1016/S0009-2614(99)00179-7
.. math:: E_{corl}^X = E_{corl}^{\infty} + \beta X^{-alpha}
"""
if type(valueLO) != type(valueHI):
raise ValidationError("corl_xtpl_helgaker_2: Inputs must be of the same datatype! (%s, %s)"
% (type(valueLO), type(valueHI)))
if alpha is None:
alpha = 3.0
if isinstance(valueLO, float):
value = (valueHI * zHI ** alpha - valueLO * zLO ** alpha) / (zHI ** alpha - zLO ** alpha)
beta = (valueHI - valueLO) / (zHI ** (-alpha) - zLO ** (-alpha))
# final = valueSCF + value
final = value
if verbose:
# Output string with extrapolation parameters
cbsscheme = """\n\n ==> Helgaker 2-point correlated extrapolation for method: %s <==\n\n""" % (functionname.upper())
# cbsscheme += """ HI-zeta (%1s) SCF Energy: % 16.12f\n""" % (str(zHI), valueSCF)
cbsscheme += """ LO-zeta (%s) Energy: % 16.12f\n""" % (str(zLO), valueLO)
cbsscheme += """ HI-zeta (%s) Energy: % 16.12f\n""" % (str(zHI), valueHI)
cbsscheme += """ Alpha (exponent) Value: % 16.12f\n""" % alpha
cbsscheme += """ Extrapolated Energy: % 16.12f\n\n""" % value
#cbsscheme += """ LO-zeta (%s) Correlation Energy: % 16.12f\n""" % (str(zLO), valueLO)
#cbsscheme += """ HI-zeta (%s) Correlation Energy: % 16.12f\n""" % (str(zHI), valueHI)
#cbsscheme += """ Beta (coefficient) Value: % 16.12f\n""" % beta
#cbsscheme += """ Extrapolated Correlation Energy: % 16.12f\n\n""" % value
name_str = "%s/(%s,%s)" % (functionname.upper(), zeta_val2sym[zLO].upper(), zeta_val2sym[zHI].upper())
cbsscheme += """ @Extrapolated """
cbsscheme += name_str + ':'
cbsscheme += " " * (19 - len(name_str))
cbsscheme += """% 16.12f\n\n""" % final
core.print_out(cbsscheme)
return final
elif isinstance(valueLO, (core.Matrix, core.Vector)):
beta = valueHI.clone()
beta.subtract(valueLO)
beta.scale(1 / (zHI ** (-alpha) - zLO ** (-alpha)))
beta.name = 'Helgaker Corl (%s, %s) beta' % (zLO, zHI)
value = valueHI.clone()
value.scale(zHI ** alpha)
tmp = valueLO.clone()
tmp.scale(zLO ** alpha)
value.subtract(tmp)
value.scale(1 / (zHI ** alpha - zLO ** alpha))
value.name = 'Helgaker Corr (%s, %s) data' % (zLO, zHI)
if verbose > 2:
core.print_out("""\n ==> Helgaker 2-point correlated extrapolation for """
"""method: %s <==\n\n""" % (functionname.upper()))
core.print_out(""" LO-zeta (%s) Data\n""" % (str(zLO)))
valueLO.print_out()
core.print_out(""" HI-zeta (%s) Data\n""" % (str(zHI)))
valueHI.print_out()
core.print_out(""" Extrapolated Data:\n""")
value.print_out()
core.print_out(""" Alpha (exponent) Value: %16.8f\n""" % alpha)
core.print_out(""" Beta Data:\n""")
beta.print_out()
# value.add(valueSCF)
return value
else:
raise ValidationError("corl_xtpl_helgaker_2: datatype is not recognized '%s'." % type(valueLO))
def return_energy_components():
VARH = {}
VARH['scf'] = {
'scf': 'SCF TOTAL ENERGY'}
VARH['hf'] = {
'hf': 'HF TOTAL ENERGY'}
VARH['mp2'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY'}
VARH['mp2.5'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY'}
VARH['mp3'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY'}
VARH['mp4(sdq)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY',
'mp4(sdq)': 'MP4(SDQ) TOTAL ENERGY'}
VARH['mp4'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY',
'mp4(sdq)': 'MP4(SDQ) TOTAL ENERGY',
'mp4': 'MP4(SDTQ) TOTAL ENERGY'}
VARH['omp2'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'omp2': 'OMP2 TOTAL ENERGY'}
VARH['omp2.5'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'omp2.5': 'OMP2.5 TOTAL ENERGY'}
VARH['omp3'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY',
'omp3': 'OMP3 TOTAL ENERGY'}
VARH['olccd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'olccd': 'OLCCD TOTAL ENERGY'}
VARH['lccd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'lccd': 'LCCD TOTAL ENERGY'}
VARH['lccsd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'lccsd': 'LCCSD TOTAL ENERGY'}
VARH['cepa(0)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'cepa(0)': 'CEPA(0) TOTAL ENERGY'}
VARH['cepa(1)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'cepa(1)': 'CEPA(1) TOTAL ENERGY'}
VARH['cepa(3)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'cepa(3)': 'CEPA(3) TOTAL ENERGY'}
VARH['acpf'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'acpf': 'ACPF TOTAL ENERGY'}
VARH['aqcc'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'aqcc': 'AQCC TOTAL ENERGY'}
VARH['qcisd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY',
'mp4(sdq)': 'MP4(SDQ) TOTAL ENERGY',
'qcisd': 'QCISD TOTAL ENERGY'}
VARH['cc2'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'cc2': 'CC2 TOTAL ENERGY'}
VARH['ccsd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'ccsd': 'CCSD TOTAL ENERGY'}
VARH['bccd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'bccd': 'CCSD TOTAL ENERGY'}
VARH['cc3'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'cc3': 'CC3 TOTAL ENERGY'}
VARH['fno-ccsd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'fno-ccsd': 'CCSD TOTAL ENERGY'}
VARH['fno-ccsd(t)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'ccsd': 'CCSD TOTAL ENERGY',
'fno-ccsd(t)': 'CCSD(T) TOTAL ENERGY'}
VARH['qcisd(t)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mp2.5': 'MP2.5 TOTAL ENERGY',
'mp3': 'MP3 TOTAL ENERGY',
'mp4(sdq)': 'MP4(SDQ) TOTAL ENERGY',
'qcisd': 'QCISD TOTAL ENERGY',
'qcisd(t)': 'QCISD(T) TOTAL ENERGY'}
VARH['ccsd(t)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'ccsd': 'CCSD TOTAL ENERGY',
'ccsd(t)': 'CCSD(T) TOTAL ENERGY'}
VARH['bccd(t)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'ccsd': 'CCSD TOTAL ENERGY',
'bccd(t)': 'CCSD(T) TOTAL ENERGY'}
VARH['cisd'] = {
'hf': 'HF TOTAL ENERGY',
'cisd': 'CISD TOTAL ENERGY'}
VARH['cisdt'] = {
'hf': 'HF TOTAL ENERGY',
'cisdt': 'CISDT TOTAL ENERGY'}
VARH['cisdtq'] = {
'hf': 'HF TOTAL ENERGY',
'cisdtq': 'CISDTQ TOTAL ENERGY'}
VARH['fci'] = {
'hf': 'HF TOTAL ENERGY',
'fci': 'FCI TOTAL ENERGY'}
VARH['mrccsd'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mrccsd': 'CCSD TOTAL ENERGY'}
VARH['mrccsd(t)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mrccsd': 'CCSD TOTAL ENERGY',
'mrccsd(t)': 'CCSD(T) TOTAL ENERGY'}
VARH['mrccsdt'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mrccsdt': 'CCSDT TOTAL ENERGY'}
VARH['mrccsdt(q)'] = {
'hf': 'HF TOTAL ENERGY',
'mp2': 'MP2 TOTAL ENERGY',
'mrccsdt': 'CCSDT TOTAL ENERGY',
'mrccsdt(q)': 'CCSDT(Q) TOTAL ENERGY'}
for cilevel in range(2, 99):
VARH['ci%s' % (str(cilevel))] = {
'hf': 'HF TOTAL ENERGY',
'ci%s' % (str(cilevel)): 'CI TOTAL ENERGY'}
for mplevel in range(5, 99):
VARH['mp%s' % (str(mplevel))] = {
'hf': 'HF TOTAL ENERGY',
'mp%s' % (str(mplevel)): 'MP%s TOTAL ENERGY' % (str(mplevel))}
for mplevel2 in range(2, mplevel):
VARH['mp%s' % (str(mplevel))]['mp%s' % (str(mplevel2))] = \
'MP%s TOTAL ENERGY' % (str(mplevel2))
# Integrate CFOUR methods
VARH.update(cfour_psivar_list())
return VARH
VARH = return_energy_components()
###################################
## Start of Complete Basis Set ##
###################################
[docs]def cbs(func, label, **kwargs):
r"""Function to define a multistage energy method from combinations of
basis set extrapolations and delta corrections and condense the
components into a minimum number of calculations.
:aliases: complete_basis_set()
:returns: (*float*) -- Total electronic energy in Hartrees
:PSI variables:
.. hlist::
:columns: 1
* :psivar:`CBS TOTAL ENERGY <CBSTOTALENERGY>`
* :psivar:`CBS REFERENCE ENERGY <CBSREFERENCEENERGY>`
* :psivar:`CBS CORRELATION ENERGY <CBSCORRELATIONENERGY>`
* :psivar:`CURRENT ENERGY <CURRENTENERGY>`
* :psivar:`CURRENT REFERENCE ENERGY <CURRENTREFERENCEENERGY>`
* :psivar:`CURRENT CORRELATION ENERGY <CURRENTCORRELATIONENERGY>`
.. caution:: Some features are not yet implemented. Buy a developer a coffee.
- No way to tell function to boost fitting basis size for all calculations.
- No way to extrapolate def2 family basis sets
- Need to add more extrapolation schemes
As represented in the equation below, a CBS energy method is defined in several
sequential stages (scf, corl, delta, delta2, delta3, delta4, delta5) covering treatment
of the reference total energy, the correlation energy, a delta correction to the
correlation energy, and a second delta correction, etc.. Each is activated by its
stage_wfn keyword and is only allowed if all preceding stages are active.
.. include:: ../cbs_eqn.rst
* Energy Methods
The presence of a stage_wfn keyword is the indicator to incorporate
(and check for stage_basis and stage_scheme keywords) and compute
that stage in defining the CBS energy.
The cbs() function requires, at a minimum, ``name='scf'`` and ``scf_basis``
keywords to be specified for reference-step only jobs and ``name`` and
``corl_basis`` keywords for correlated jobs.
The following energy methods have been set up for cbs().
.. hlist::
:columns: 5
* scf
* hf
* mp2
* mp2.5
* mp3
* mp4(sdq)
* mp4
* mp\ *n*
* omp2
* omp2.5
* omp3
* olccd
* lccd
* lccsd
* cepa(0)
* cepa(1)
* cepa(3)
* acpf
* aqcc
* qcisd
* cc2
* ccsd
* fno-ccsd
* bccd
* cc3
* qcisd(t)
* ccsd(t)
* fno-ccsd(t)
* bccd(t)
* cisd
* cisdt
* cisdtq
* ci\ *n*
* fci
* mrccsd
* mrccsd(t)
* mrccsdt
* mrccsdt(q)
:type name: string
:param name: ``'scf'`` || ``'ccsd'`` || etc.
First argument, usually unlabeled. Indicates the computational method
for the correlation energy, unless only reference step to be performed,
in which case should be ``'scf'``. Overruled if stage_wfn keywords supplied.
:type scf_wfn: string
:param scf_wfn: |dl| ``'scf'`` |dr| || ``'c4-scf'`` || etc.
Indicates the energy method for which the reference energy is to be
obtained. Generally unnecessary, as 'scf' is *the* scf in |PSIfour| but
can be used to direct lone scf components to run in |PSIfour| or Cfour
in a mixed-program composite method.
:type corl_wfn: string
:param corl_wfn: ``'mp2'`` || ``'ccsd(t)'`` || etc.
Indicates the energy method for which the correlation energy is to be
obtained. Can also be specified with ``name`` or as the unlabeled
first argument to the function.
:type delta_wfn: string
:param delta_wfn: ``'ccsd'`` || ``'ccsd(t)'`` || etc.
Indicates the (superior) energy method for which a delta correction
to the correlation energy is to be obtained.
:type delta_wfn_lesser: string
:param delta_wfn_lesser: |dl| ``corl_wfn`` |dr| || ``'mp2'`` || etc.
Indicates the inferior energy method for which a delta correction
to the correlation energy is to be obtained.
:type delta2_wfn: string
:param delta2_wfn: ``'ccsd'`` || ``'ccsd(t)'`` || etc.
Indicates the (superior) energy method for which a second delta correction
to the correlation energy is to be obtained.
:type delta2_wfn_lesser: string
:param delta2_wfn_lesser: |dl| ``delta_wfn`` |dr| || ``'ccsd(t)'`` || etc.
Indicates the inferior energy method for which a second delta correction
to the correlation energy is to be obtained.
:type delta3_wfn: string
:param delta3_wfn: ``'ccsd'`` || ``'ccsd(t)'`` || etc.
Indicates the (superior) energy method for which a third delta correction
to the correlation energy is to be obtained.
:type delta3_wfn_lesser: string
:param delta3_wfn_lesser: |dl| ``delta2_wfn`` |dr| || ``'ccsd(t)'`` || etc.
Indicates the inferior energy method for which a third delta correction
to the correlation energy is to be obtained.
:type delta4_wfn: string
:param delta4_wfn: ``'ccsd'`` || ``'ccsd(t)'`` || etc.
Indicates the (superior) energy method for which a fourth delta correction
to the correlation energy is to be obtained.
:type delta4_wfn_lesser: string
:param delta4_wfn_lesser: |dl| ``delta3_wfn`` |dr| || ``'ccsd(t)'`` || etc.
Indicates the inferior energy method for which a fourth delta correction
to the correlation energy is to be obtained.
:type delta5_wfn: string
:param delta5_wfn: ``'ccsd'`` || ``'ccsd(t)'`` || etc.
Indicates the (superior) energy method for which a fifth delta correction
to the correlation energy is to be obtained.
:type delta5_wfn_lesser: string
:param delta5_wfn_lesser: |dl| ``delta4_wfn`` |dr| || ``'ccsd(t)'`` || etc.
Indicates the inferior energy method for which a fifth delta correction
to the correlation energy is to be obtained.
* Basis Sets
Currently, the basis set set through ``set`` commands have no influence
on a cbs calculation.
:type scf_basis: :ref:`basis string <apdx:basisElement>`
:param scf_basis: |dl| ``corl_basis`` |dr| || ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the reference energy.
If any correlation method is specified, ``scf_basis`` can default
to ``corl_basis``.
:type corl_basis: :ref:`basis string <apdx:basisElement>`
:param corl_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the correlation energy.
:type delta_basis: :ref:`basis string <apdx:basisElement>`
:param delta_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the delta correction
to the correlation energy.
:type delta2_basis: :ref:`basis string <apdx:basisElement>`
:param delta2_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the second delta correction
to the correlation energy.
:type delta3_basis: :ref:`basis string <apdx:basisElement>`
:param delta3_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the third delta correction
to the correlation energy.
:type delta4_basis: :ref:`basis string <apdx:basisElement>`
:param delta4_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the fourth delta correction
to the correlation energy.
:type delta5_basis: :ref:`basis string <apdx:basisElement>`
:param delta5_basis: ``'cc-pV[TQ]Z'`` || ``'jun-cc-pv[tq5]z'`` || ``'6-31G*'`` || etc.
Indicates the sequence of basis sets employed for the fifth delta correction
to the correlation energy.
* Schemes
Transformations of the energy through basis set extrapolation for each
stage of the CBS definition. A complaint is generated if number of basis
sets in stage_basis does not exactly satisfy requirements of stage_scheme.
An exception is the default, ``'xtpl_highest_1'``, which uses the best basis
set available. See :ref:`sec:cbs_xtpl` for all available schemes.
:type scf_scheme: function
:param scf_scheme: |dl| ``xtpl_highest_1`` |dr| || ``scf_xtpl_helgaker_3`` || etc.
Indicates the basis set extrapolation scheme to be applied to the reference energy.
Defaults to :py:func:`~scf_xtpl_helgaker_3` if three valid basis sets
present in ``scf_basis``, :py:func:`~scf_xtpl_helgaker_2` if two valid basis
sets present in ``scf_basis``, and :py:func:`~xtpl_highest_1` otherwise.
.. hlist::
:columns: 1
* xtpl_highest_1
* scf_xtpl_helgaker_3
* scf_xtpl_helgaker_2
* scf_xtpl_truhlar_2
* scf_xtpl_karton_2
:type corl_scheme: function
:param corl_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``corl_basis`` and :py:func:`~xtpl_highest_1` otherwise.
.. hlist::
:columns: 1
* xtpl_highest_1
* corl_xtpl_helgaker_2
:type delta_scheme: function
:param delta_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the delta correction
to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``delta_basis`` and :py:func:`~xtpl_highest_1` otherwise.
.. hlist::
:columns: 1
* xtpl_highest_1
* corl_xtpl_helgaker_2
:type delta2_scheme: function
:param delta2_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the second delta correction
to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``delta2_basis`` and :py:func:`~xtpl_highest_1` otherwise.
.. hlist::
:columns: 1
* xtpl_highest_1
* corl_xtpl_helgaker_2
:type delta3_scheme: function
:param delta3_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the third delta correction
to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``delta3_basis`` and :py:func:`~xtpl_highest_1` otherwise.
:type delta4_scheme: function
:param delta4_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the fourth delta correction
to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``delta4_basis`` and :py:func:`~xtpl_highest_1` otherwise.
:type delta5_scheme: function
:param delta5_scheme: |dl| ``xtpl_highest_1`` |dr| || ``corl_xtpl_helgaker_2`` || etc.
Indicates the basis set extrapolation scheme to be applied to the fifth delta correction
to the correlation energy.
Defaults to :py:func:`~corl_xtpl_helgaker_2` if two valid basis sets
present in ``delta5_basis`` and :py:func:`~xtpl_highest_1` otherwise.
:type scf_alpha: float
Overrides the default \alpha parameter used in the listed SCF extrapolation procedures.
Has no effect on others, including :py:func:`~xtpl_highest_1` and :py:func:`~scf_xtpl_helgaker_3`.
.. hlist::
:columns: 1
* :py:func:`scf_xtpl_helgaker_2`
* :py:func:`scf_xtpl_truhlar_2`
* :py:func:`scf_xtpl_karton_2`
:type corl_alpha: float
Overrides the default \alpha parameter used in the listed :py:func:`corl_xtpl_helgaker_2` correlation
extrapolation to the corl stage. The supplied \alpha does not impact delta or any further stages.
.. hlist::
:columns: 1
* :py:func:`corl_xtpl_helgaker_2`
:type delta_alpha: float
Overrides the default \alpha parameter used in the listed
:py:func:`corl_xtpl_helgaker_2` correlation extrapolation for the delta correction. Useful when
delta correction is performed using smaller basis sets for which a different \alpha might
be more appropriate.
.. hlist::
:columns: 1
* :py:func:`corl_xtpl_helgaker_2`
:type molecule: :ref:`molecule <op_py_molecule>`
:param molecule: ``h2o`` || etc.
The target molecule, if not the last molecule defined.
:examples:
>>> # [1] replicates with cbs() the simple model chemistry scf/cc-pVDZ: set basis cc-pVDZ energy('scf')
>>> cbs(name='scf', scf_basis='cc-pVDZ')
>>> # [2] replicates with cbs() the simple model chemistry mp2/jun-cc-pVDZ: set basis jun-cc-pVDZ energy('mp2')
>>> cbs(name='mp2', corl_basis='jun-cc-pVDZ')
>>> # [3] DTQ-zeta extrapolated scf reference energy
>>> cbs(name='scf', scf_basis='cc-pV[DTQ]Z', scf_scheme=scf_xtpl_helgaker_3)
>>> # [4] DT-zeta extrapolated mp2 correlation energy atop a T-zeta reference
>>> cbs(corl_wfn='mp2', corl_basis='cc-pv[dt]z', corl_scheme=corl_xtpl_helgaker_2)
>>> # [5] a DT-zeta extrapolated coupled-cluster correction atop a TQ-zeta extrapolated mp2 correlation energy atop a Q-zeta reference (both equivalent)
>>> cbs(corl_wfn='mp2', corl_basis='aug-cc-pv[tq]z', delta_wfn='ccsd(t)', delta_basis='aug-cc-pv[dt]z')
>>> cbs(energy, wfn='mp2', corl_basis='aug-cc-pv[tq]z', corl_scheme=corl_xtpl_helgaker_2, delta_wfn='ccsd(t)', delta_basis='aug-cc-pv[dt]z', delta_scheme=corl_xtpl_helgaker_2)
>>> # [6] a D-zeta ccsd(t) correction atop a DT-zeta extrapolated ccsd cluster correction atop a TQ-zeta extrapolated mp2 correlation energy atop a Q-zeta reference
>>> cbs(name='mp2', corl_basis='aug-cc-pv[tq]z', corl_scheme=corl_xtpl_helgaker_2, delta_wfn='ccsd', delta_basis='aug-cc-pv[dt]z', delta_scheme=corl_xtpl_helgaker_2, delta2_wfn='ccsd(t)', delta2_wfn_lesser='ccsd', delta2_basis='aug-cc-pvdz')
>>> # [7] cbs() coupled with database()
>>> TODO database('mp2', 'BASIC', subset=['h2o','nh3'], symm='on', func=cbs, corl_basis='cc-pV[tq]z', corl_scheme=corl_xtpl_helgaker_2, delta_wfn='ccsd(t)', delta_basis='sto-3g')
>>> # [8] cbs() coupled with optimize()
>>> TODO optimize('mp2', corl_basis='cc-pV[DT]Z', corl_scheme=corl_xtpl_helgaker_2, func=cbs)
"""
kwargs = p4util.kwargs_lower(kwargs)
return_wfn = kwargs.pop('return_wfn', False)
verbose = kwargs.pop('verbose', 0)
ptype = kwargs.pop('ptype')
cbs_alpha = {
'scf': kwargs.get('cbs_scf_alpha', kwargs.get('scf_alpha', None)),
'corl': kwargs.get('cbs_corl_alpha', kwargs.get('corl_alpha', None)),
'delta': kwargs.get('cbs_delta_alpha', kwargs.get('delta_alpha', None)),
'delta2': kwargs.get('cbs_delta2_alpha', kwargs.get('delta2_alpha', None)),
}
# Establish function to call (only energy makes sense for cbs)
if ptype not in ['energy', 'gradient', 'hessian']:
raise ValidationError("""Wrapper complete_basis_set is unhappy to be calling function '%s' instead of 'energy'.""" % ptype)
optstash = p4util.OptionsState(
['BASIS'],
['WFN'],
['WRITER_FILE_LABEL'])
# Define some quantum chemical knowledge, namely what methods are subsumed in others
do_scf = True
do_corl = False
do_delta = False
do_delta2 = False
do_delta3 = False
do_delta4 = False
do_delta5 = False
user_writer_file_label = core.get_global_option('WRITER_FILE_LABEL')
# Make sure the molecule the user provided is the active one
molecule = kwargs.pop('molecule', core.get_active_molecule())
molecule.update_geometry()
molstr = molecule.create_psi4_string_from_molecule()
natom = molecule.natom()
# Establish method for correlation energy
cbs_corl_wfn = kwargs.pop('corl_wfn', '').lower()
if cbs_corl_wfn:
do_corl = True
# Establish method for reference energy
if do_corl and cbs_corl_wfn.startswith('c4-'):
default_scf = 'c4-hf'
else:
default_scf = 'hf'
cbs_scf_wfn = kwargs.pop('scf_wfn', default_scf).lower()
if do_scf:
if cbs_scf_wfn not in VARH.keys():
raise ValidationError("""Requested SCF method '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_scf_wfn))
# ... resume correlation logic
if do_corl:
if cbs_corl_wfn not in VARH.keys():
raise ValidationError("""Requested CORL method '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_corl_wfn))
cbs_corl_wfn_lesser = kwargs.get('corl_wfn_lesser', cbs_scf_wfn).lower()
if cbs_corl_wfn_lesser not in VARH.keys():
raise ValidationError("""Requested CORL method lesser '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta_wfn_lesser))
# Establish method for delta correction energy
if 'delta_wfn' in kwargs:
do_delta = True
cbs_delta_wfn = kwargs['delta_wfn'].lower()
if cbs_delta_wfn not in VARH.keys():
raise ValidationError("""Requested DELTA method '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta_wfn))
cbs_delta_wfn_lesser = kwargs.get('delta_wfn_lesser', cbs_corl_wfn).lower()
if cbs_delta_wfn_lesser not in VARH.keys():
raise ValidationError("""Requested DELTA method lesser '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta_wfn_lesser))
# Establish method for second delta correction energy
if 'delta2_wfn' in kwargs:
do_delta2 = True
cbs_delta2_wfn = kwargs['delta2_wfn'].lower()
if cbs_delta2_wfn not in VARH.keys():
raise ValidationError("""Requested DELTA2 method '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta2_wfn))
cbs_delta2_wfn_lesser = kwargs.get('delta2_wfn_lesser', cbs_delta_wfn).lower()
if cbs_delta2_wfn_lesser not in VARH.keys():
raise ValidationError("""Requested DELTA2 method lesser '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta2_wfn_lesser))
# # Establish method for third delta correction energy
# if 'delta3_wfn' in kwargs:
# do_delta3 = True
# cbs_delta3_wfn = kwargs['delta3_wfn'].lower()
# if cbs_delta3_wfn not in VARH.keys():
# raise ValidationError("""Requested DELTA3 method '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta3_wfn))
#
# cbs_delta3_wfn_lesser = kwargs.get('delta3_wfn_lesser', cbs_delta2_wfn).lower()
# if not (cbs_delta3_wfn_lesser in VARH.keys()):
# raise ValidationError("""Requested DELTA3 method lesser '%s' is not recognized. Add it to VARH in wrapper.py to proceed.""" % (cbs_delta3_wfn_lesser))
#
# # Establish method for fourth delta correction energy
# if 'delta4_wfn' in kwargs:
# do_delta4 = True
# cbs_delta4_wfn = kwargs['delta4_wfn'].lower()
# if not (cbs_delta4_wfn in VARH.keys()):
# raise ValidationError('Requested DELTA4 method \'%s\' is not recognized. Add it to VARH in wrapper.py to proceed.' % (cbs_delta4_wfn))
#
# if 'delta4_wfn_lesser' in kwargs:
# cbs_delta4_wfn_lesser = kwargs['delta4_wfn_lesser'].lower()
# else:
# cbs_delta4_wfn_lesser = cbs_delta3_wfn
# if not (cbs_delta4_wfn_lesser in VARH.keys()):
# raise ValidationError('Requested DELTA4 method lesser \'%s\' is not recognized. Add it to VARH in wrapper.py to proceed.' % (cbs_delta4_wfn_lesser))
#
# # Establish method for fifth delta correction energy
# if 'delta5_wfn' in kwargs:
# do_delta5 = True
# cbs_delta5_wfn = kwargs['delta5_wfn'].lower()
# if not (cbs_delta5_wfn in VARH.keys()):
# raise ValidationError('Requested DELTA5 method \'%s\' is not recognized. Add it to VARH in wrapper.py to proceed.' % (cbs_delta5_wfn))
#
# if 'delta5_wfn_lesser' in kwargs:
# cbs_delta5_wfn_lesser = kwargs['delta5_wfn_lesser'].lower()
# else:
# cbs_delta5_wfn_lesser = cbs_delta4_wfn
# if not (cbs_delta5_wfn_lesser in VARH.keys()):
# raise ValidationError('Requested DELTA5 method lesser \'%s\' is not recognized. Add it to VARH in wrapper.py to proceed.' % (cbs_delta5_wfn_lesser))
# Check that user isn't skipping steps in scf + corl + delta + delta2 sequence
if do_scf and not do_corl and not do_delta and not do_delta2 and not do_delta3 and not do_delta4 and not do_delta5:
pass
elif do_scf and do_corl and not do_delta and not do_delta2 and not do_delta3 and not do_delta4 and not do_delta5:
pass
elif do_scf and do_corl and do_delta and not do_delta2 and not do_delta3 and not do_delta4 and not do_delta5:
pass
elif do_scf and do_corl and do_delta and do_delta2 and not do_delta3 and not do_delta4 and not do_delta5:
pass
#elif do_scf and do_corl and do_delta and do_delta2 and do_delta3 and not do_delta4 and not do_delta5:
# pass
#elif do_scf and do_corl and do_delta and do_delta2 and do_delta3 and do_delta4 and not do_delta5:
# pass
#elif do_scf and do_corl and do_delta and do_delta2 and do_delta3 and do_delta4 and do_delta5:
# pass
else:
raise ValidationError('Requested scf (%s) + corl (%s) + delta (%s) + delta2 (%s) + delta3 (%s) + delta4 (%s) + delta5 (%s) not valid. These steps are cummulative.' %
(do_scf, do_corl, do_delta, do_delta2, do_delta3, do_delta4, do_delta5))
# Establish list of valid basis sets for correlation energy
if do_corl:
if 'corl_basis' in kwargs:
BSTC, ZETC = _expand_bracketed_basis(kwargs['corl_basis'].lower(), molecule=molstr)
else:
raise ValidationError("""CORL basis sets through keyword '%s' are required.""" % ('corl_basis'))
# Establish list of valid basis sets for scf energy
if 'scf_basis' in kwargs:
BSTR, ZETR = _expand_bracketed_basis(kwargs['scf_basis'].lower(), molecule=molstr)
elif do_corl:
BSTR = BSTC[:]
ZETR = ZETC[:]
else:
raise ValidationError("""SCF basis sets through keyword '%s' are required. Or perhaps you forgot the '%s'.""" % ('scf_basis', 'corl_wfn'))
# Establish list of valid basis sets for delta correction energy
if do_delta:
if 'delta_basis' in kwargs:
BSTD, ZETD = _expand_bracketed_basis(kwargs['delta_basis'].lower(), molecule=molstr)
else:
raise ValidationError("""DELTA basis sets through keyword '%s' are required.""" % ('delta_basis'))
# Establish list of valid basis sets for second delta correction energy
if do_delta2:
if 'delta2_basis' in kwargs:
BSTD2, ZETD2 = _expand_bracketed_basis(kwargs['delta2_basis'].lower(), molecule=molstr)
else:
raise ValidationError("""DELTA2 basis sets through keyword '%s' are required.""" % ('delta2_basis'))
# # Establish list of valid basis sets for third delta correction energy
# if do_delta3:
# if 'delta3_basis' in kwargs:
# BSTD3, ZETD3 = validate_bracketed_basis(kwargs['delta3_basis'].lower())
# else:
# raise ValidationError('DELTA3 basis sets through keyword \'%s\' are required.' % ('delta3_basis'))
#
# # Establish list of valid basis sets for fourth delta correction energy
# if do_delta4:
# if 'delta4_basis' in kwargs:
# BSTD4, ZETD4 = validate_bracketed_basis(kwargs['delta4_basis'].lower())
# else:
# raise ValidationError('DELTA4 basis sets through keyword \'%s\' are required.' % ('delta4_basis'))
#
# # Establish list of valid basis sets for fifth delta correction energy
# if do_delta5:
# if 'delta5_basis' in kwargs:
# BSTD5, ZETD5 = validate_bracketed_basis(kwargs['delta5_basis'].lower())
# else:
# raise ValidationError('DELTA5 basis sets through keyword \'%s\' are required.' % ('delta5_basis'))
# Establish treatment for scf energy (validity check useless since python will catch it long before here)
if (len(BSTR) == 3) and ('scf_basis' in kwargs):
cbs_scf_scheme = scf_xtpl_helgaker_3
elif (len(BSTR) == 2) and ('scf_basis' in kwargs):
cbs_scf_scheme = scf_xtpl_helgaker_2
elif (len(BSTR) == 1) and ('scf_basis' in kwargs):
cbs_scf_scheme = xtpl_highest_1
elif 'scf_basis' in kwargs:
raise ValidationError("""SCF basis sets of number %d cannot be handled.""" % (len(BSTR)))
elif do_corl:
cbs_scf_scheme = xtpl_highest_1
BSTR = [BSTC[-1]]
ZETR = [ZETC[-1]]
if 'scf_scheme' in kwargs:
cbs_scf_scheme = kwargs['scf_scheme']
# Establish treatment for correlation energy
if do_corl:
if len(BSTC) == 2:
cbs_corl_scheme = corl_xtpl_helgaker_2
elif len(BSTC) > 2:
raise ValidationError("""Cannot extrapolate correlation with %d basis sets. Use highest 2.""" % (len(BSTC)))
else:
cbs_corl_scheme = xtpl_highest_1
if 'corl_scheme' in kwargs:
cbs_corl_scheme = kwargs['corl_scheme']
# Establish treatment for delta correction energy
if do_delta:
if len(BSTD) == 2:
cbs_delta_scheme = corl_xtpl_helgaker_2
else:
cbs_delta_scheme = xtpl_highest_1
if 'delta_scheme' in kwargs:
cbs_delta_scheme = kwargs['delta_scheme']
# Establish treatment for delta2 correction energy
if do_delta2:
if len(BSTD2) == 2:
cbs_delta2_scheme = corl_xtpl_helgaker_2
else:
cbs_delta2_scheme = xtpl_highest_1
if 'delta2_scheme' in kwargs:
cbs_delta2_scheme = kwargs['delta2_scheme']
# # Establish treatment for delta3 correction energy
# if do_delta3:
# if len(BSTD3) == 2:
# cbs_delta3_scheme = corl_xtpl_helgaker_2
# else:
# cbs_delta3_scheme = xtpl_highest_1
# if 'delta3_scheme' in kwargs:
# cbs_delta3_scheme = kwargs['delta3_scheme']
#
# # Establish treatment for delta4 correction energy
# if do_delta4:
# if len(BSTD4) == 2:
# cbs_delta4_scheme = corl_xtpl_helgaker_2
# else:
# cbs_delta4_scheme = xtpl_highest_1
# if 'delta4_scheme' in kwargs:
# cbs_delta4_scheme = kwargs['delta4_scheme']
#
# # Establish treatment for delta5 correction energy
# if do_delta5:
# if len(BSTD5) == 2:
# cbs_delta5_scheme = corl_xtpl_helgaker_2
# else:
# cbs_delta5_scheme = xtpl_highest_1
# if 'delta5_scheme' in kwargs:
# cbs_delta5_scheme = kwargs['delta5_scheme']
# Build string of title banner
cbsbanners = ''
cbsbanners += """core.print_out('\\n')\n"""
cbsbanners += """p4util.banner(' CBS Setup: %s ' % label)\n"""
cbsbanners += """core.print_out('\\n')\n\n"""
exec(cbsbanners)
# Call schemes for each portion of total energy to 'place orders' for calculations needed
d_fields = ['d_stage', 'd_scheme', 'd_basis', 'd_wfn', 'd_need', 'd_coef', 'd_energy', 'd_gradient', 'd_hessian']
f_fields = ['f_wfn', 'f_basis', 'f_zeta', 'f_energy', 'f_gradient', 'f_hessian']
GRAND_NEED = []
MODELCHEM = []
if do_scf:
NEED = _expand_scheme_orders(cbs_scf_scheme, BSTR, ZETR, cbs_scf_wfn, natom)
GRAND_NEED.append(dict(zip(d_fields, ['scf', cbs_scf_scheme,
_contract_bracketed_basis(BSTR), cbs_scf_wfn, NEED, +1, 0.0, None, None])))
if do_corl:
NEED = _expand_scheme_orders(cbs_corl_scheme, BSTC, ZETC, cbs_corl_wfn, natom)
GRAND_NEED.append(dict(zip(d_fields, ['corl', cbs_corl_scheme,
_contract_bracketed_basis(BSTC), cbs_corl_wfn, NEED, +1, 0.0, None, None])))
NEED = _expand_scheme_orders(cbs_corl_scheme, BSTC, ZETC, cbs_corl_wfn_lesser, natom)
GRAND_NEED.append(dict(zip(d_fields, ['corl', cbs_corl_scheme,
_contract_bracketed_basis(BSTC), cbs_corl_wfn_lesser, NEED, -1, 0.0, None, None])))
if do_delta:
NEED = _expand_scheme_orders(cbs_delta_scheme, BSTD, ZETD, cbs_delta_wfn, natom)
GRAND_NEED.append(dict(zip(d_fields, ['delta', cbs_delta_scheme,
_contract_bracketed_basis(BSTD), cbs_delta_wfn, NEED, +1, 0.0, None, None])))
NEED = _expand_scheme_orders(cbs_delta_scheme, BSTD, ZETD, cbs_delta_wfn_lesser, natom)
GRAND_NEED.append(dict(zip(d_fields, ['delta', cbs_delta_scheme,
_contract_bracketed_basis(BSTD), cbs_delta_wfn_lesser, NEED, -1, 0.0, None, None])))
if do_delta2:
NEED = _expand_scheme_orders(cbs_delta2_scheme, BSTD2, ZETD2, cbs_delta2_wfn, natom)
GRAND_NEED.append(dict(zip(d_fields, ['delta2', cbs_delta2_scheme,
_contract_bracketed_basis(BSTD2), cbs_delta2_wfn, NEED, +1, 0.0, None, None])))
NEED = _expand_scheme_orders(cbs_delta2_scheme, BSTD2, ZETD2, cbs_delta2_wfn_lesser, natom)
GRAND_NEED.append(dict(zip(d_fields, ['delta2', cbs_delta2_scheme,
_contract_bracketed_basis(BSTD2), cbs_delta2_wfn_lesser, NEED, -1, 0.0, None, None])))
# if do_delta3:
# NEED = call_function_in_1st_argument(cbs_delta3_scheme,
# mode='requisition', basisname=BSTD3, basiszeta=ZETD3, wfnname=cbs_delta3_wfn)
# GRAND_NEED.append(dict(zip(d_fields, ['delta3', cbs_delta3_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta3_wfn, NEED, +1, 0.0])))
#
# NEED = call_function_in_1st_argument(cbs_delta3_scheme,
# mode='requisition', basisname=BSTD3, basiszeta=ZETD3, wfnname=cbs_delta3_wfn_lesser)
# GRAND_NEED.append(dict(zip(d_fields, ['delta3', cbs_delta3_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta3_wfn_lesser, NEED, -1, 0.0])))
#
# if do_delta4:
# NEED = call_function_in_1st_argument(cbs_delta4_scheme,
# mode='requisition', basisname=BSTD4, basiszeta=ZETD4, wfnname=cbs_delta4_wfn)
# GRAND_NEED.append(dict(zip(d_fields, ['delta4', cbs_delta4_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta4_wfn, NEED, +1, 0.0])))
#
# NEED = call_function_in_1st_argument(cbs_delta4_scheme,
# mode='requisition', basisname=BSTD4, basiszeta=ZETD4, wfnname=cbs_delta4_wfn_lesser)
# GRAND_NEED.append(dict(zip(d_fields, ['delta4', cbs_delta4_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta4_wfn_lesser, NEED, -1, 0.0])))
#
# if do_delta5:
# NEED = call_function_in_1st_argument(cbs_delta5_scheme,
# mode='requisition', basisname=BSTD5, basiszeta=ZETD5, wfnname=cbs_delta5_wfn)
# GRAND_NEED.append(dict(zip(d_fields, ['delta5', cbs_delta5_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta5_wfn, NEED, +1, 0.0])))
#
# NEED = call_function_in_1st_argument(cbs_delta5_scheme,
# mode='requisition', basisname=BSTD5, basiszeta=ZETD5, wfnname=cbs_delta5_wfn_lesser)
# GRAND_NEED.append(dict(zip(d_fields, ['delta5', cbs_delta5_scheme,
# reconstitute_bracketed_basis(NEED), cbs_delta5_wfn_lesser, NEED, -1, 0.0])))
for stage in GRAND_NEED:
for lvl in stage['d_need'].items():
MODELCHEM.append(lvl[1])
# Apply chemical reasoning to choose the minimum computations to run
JOBS = MODELCHEM[:]
addlremark = {'energy': '', 'gradient': ', GRADIENT', 'hessian': ', HESSIAN'}
instructions = ''
instructions += """ Naive listing of computations required.\n"""
for mc in JOBS:
instructions += """ %12s / %-24s for %s%s\n""" % \
(mc['f_wfn'], mc['f_basis'], VARH[mc['f_wfn']][mc['f_wfn']], addlremark[ptype])
# Remove duplicate modelchem portion listings
for mc in MODELCHEM:
dups = -1
for indx_job, job in enumerate(JOBS):
if (job['f_wfn'] == mc['f_wfn']) and (job['f_basis'] == mc['f_basis']):
dups += 1
if dups >= 1:
del JOBS[indx_job]
# Remove chemically subsumed modelchem portion listings
if ptype == 'energy':
for mc in MODELCHEM:
for wfn in VARH[mc['f_wfn']]:
for indx_job, job in enumerate(JOBS):
if (VARH[mc['f_wfn']][wfn] == VARH[job['f_wfn']][job['f_wfn']]) and \
(mc['f_basis'] == job['f_basis']) and not \
(mc['f_wfn'] == job['f_wfn']):
del JOBS[indx_job]
instructions += """\n Enlightened listing of computations required.\n"""
for mc in JOBS:
instructions += """ %12s / %-24s for %s%s\n""" % \
(mc['f_wfn'], mc['f_basis'], VARH[mc['f_wfn']][mc['f_wfn']], addlremark[ptype])
# Expand listings to all that will be obtained
JOBS_EXT = []
for job in JOBS:
for wfn in VARH[job['f_wfn']]:
JOBS_EXT.append(dict(zip(f_fields, [wfn, job['f_basis'], job['f_zeta'],
0.0,
core.Matrix(natom, 3),
core.Matrix(3 * natom, 3 * natom)])))
instructions += """\n Full listing of computations to be obtained (required and bonus).\n"""
for mc in JOBS_EXT:
instructions += """ %12s / %-24s for %s%s\n""" % \
(mc['f_wfn'], mc['f_basis'], VARH[mc['f_wfn']][mc['f_wfn']], addlremark[ptype])
core.print_out(instructions)
psioh = core.IOManager.shared_object()
psioh.set_specific_retention(constants.PSIF_SCF_MOS, True)
# projection across point groups not allowed and cbs() usually a mix of symm-enabled and symm-tol calls
# needs to be communicated to optimize() so reset by that optstash
core.set_local_option('SCF', 'GUESS_PERSIST', True)
Njobs = 0
# Run necessary computations
for mc in JOBS:
kwargs['name'] = mc['f_wfn']
# Build string of title banner
cbsbanners = ''
cbsbanners += """core.print_out('\\n')\n"""
cbsbanners += """p4util.banner(' CBS Computation: %s / %s%s ')\n""" % \
(mc['f_wfn'].upper(), mc['f_basis'].upper(), addlremark[ptype])
cbsbanners += """core.print_out('\\n')\n\n"""
exec(cbsbanners)
# Build string of molecule and commands that are dependent on the database
commands = '\n'
commands += """\ncore.set_global_option('BASIS', '%s')\n""" % (mc['f_basis'])
commands += """core.set_global_option('WRITER_FILE_LABEL', '%s')\n""" % \
(user_writer_file_label + ('' if user_writer_file_label == '' else '-') + mc['f_wfn'].lower() + '-' + mc['f_basis'].lower())
exec(commands)
# Make energy(), etc. call
response = func(molecule=molecule, **kwargs)
if ptype == 'energy':
mc['f_energy'] = response
elif ptype == 'gradient':
mc['f_gradient'] = response
mc['f_energy'] = core.get_variable('CURRENT ENERGY')
if verbose > 1:
mc['f_gradient'].print_out()
elif ptype == 'hessian':
mc['f_hessian'] = response
mc['f_energy'] = core.get_variable('CURRENT ENERGY')
if verbose > 1:
mc['f_hessian'].print_out()
Njobs += 1
if verbose > 1:
core.print_out("\nCURRENT ENERGY: %14.16f\n" % mc['f_energy'])
# Fill in energies for subsumed methods
if ptype == 'energy':
for wfn in VARH[mc['f_wfn']]:
for job in JOBS_EXT:
if (wfn == job['f_wfn']) and (mc['f_basis'] == job['f_basis']):
job['f_energy'] = core.get_variable(VARH[wfn][wfn])
if verbose > 1:
core.print_variables()
core.clean_variables()
core.clean()
# Copy data from 'run' to 'obtained' table
for mce in JOBS_EXT:
if (mc['f_wfn'] == mce['f_wfn']) and (mc['f_basis'] == mce['f_basis']):
mce['f_energy'] = mc['f_energy']
mce['f_gradient'] = mc['f_gradient']
mce['f_hessian'] = mc['f_hessian']
psioh.set_specific_retention(constants.PSIF_SCF_MOS, False)
# Build string of title banner
cbsbanners = ''
cbsbanners += """core.print_out('\\n')\n"""
cbsbanners += """p4util.banner(' CBS Results: %s ' % label)\n"""
cbsbanners += """core.print_out('\\n')\n\n"""
exec(cbsbanners)
# Insert obtained energies into the array that stores the cbs stages
for stage in GRAND_NEED:
for lvl in stage['d_need'].items():
MODELCHEM.append(lvl[1])
for job in JOBS_EXT:
# Dont ask
if (((lvl[1]['f_wfn'] == job['f_wfn']) or
((lvl[1]['f_wfn'][3:] == job['f_wfn']) and lvl[1]['f_wfn'].startswith('c4-')) or
((lvl[1]['f_wfn'] == job['f_wfn'][3:]) and job['f_wfn'].startswith('c4-')) or
(('c4-' + lvl[1]['f_wfn']) == job['f_wfn']) or
(lvl[1]['f_wfn'] == ('c4-' + job['f_wfn']))) and
(lvl[1]['f_basis'] == job['f_basis'])):
lvl[1]['f_energy'] = job['f_energy']
lvl[1]['f_gradient'] = job['f_gradient']
lvl[1]['f_hessian'] = job['f_hessian']
# Make xtpl() call
finalenergy = 0.0
finalgradient = core.Matrix(natom, 3)
finalhessian = core.Matrix(3 * natom, 3 * natom)
for stage in GRAND_NEED:
hiloargs = {'alpha': cbs_alpha[stage['d_stage']]}
hiloargs.update(_contract_scheme_orders(stage['d_need'], 'f_energy'))
stage['d_energy'] = stage['d_scheme'](**hiloargs)
finalenergy += stage['d_energy'] * stage['d_coef']
if ptype == 'gradient':
hiloargs.update(_contract_scheme_orders(stage['d_need'], 'f_gradient'))
stage['d_gradient'] = stage['d_scheme'](**hiloargs)
work = stage['d_gradient'].clone()
work.scale(stage['d_coef'])
finalgradient.add(work)
elif ptype == 'hessian':
hiloargs.update(_contract_scheme_orders(stage['d_need'], 'f_hessian'))
stage['d_hessian'] = stage['d_scheme'](**hiloargs)
work = stage['d_hessian'].clone()
work.scale(stage['d_coef'])
finalhessian.add(work)
# Build string of results table
table_delimit = ' ' + '-' * 105 + '\n'
tables = ''
tables += """\n ==> %s <==\n\n""" % ('Components')
tables += table_delimit
tables += """ %6s %20s %1s %-26s %3s %16s %-s\n""" % ('', 'Method', '/', 'Basis', 'Rqd', 'Energy [Eh]', 'Variable')
tables += table_delimit
for job in JOBS_EXT:
star = ''
for mc in MODELCHEM:
if (job['f_wfn'] == mc['f_wfn']) and (job['f_basis'] == mc['f_basis']):
star = '*'
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % ('', job['f_wfn'],
'/', job['f_basis'], star, job['f_energy'], VARH[job['f_wfn']][job['f_wfn']])
tables += table_delimit
tables += """\n ==> %s <==\n\n""" % ('Stages')
tables += table_delimit
tables += """ %6s %20s %1s %-27s %2s %16s %-s\n""" % ('Stage', 'Method', '/', 'Basis', 'Wt', 'Energy [Eh]', 'Scheme')
tables += table_delimit
for stage in GRAND_NEED:
tables += """ %6s %20s %1s %-27s %2d %16.8f %-s\n""" % (stage['d_stage'], stage['d_wfn'],
'/', stage['d_basis'], stage['d_coef'], stage['d_energy'], stage['d_scheme'].__name__)
tables += table_delimit
tables += """\n ==> %s <==\n\n""" % ('CBS')
tables += table_delimit
tables += """ %6s %20s %1s %-27s %2s %16s %-s\n""" % ('Stage', 'Method', '/', 'Basis', '', 'Energy [Eh]', 'Scheme')
tables += table_delimit
if do_scf:
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[0]['d_stage'],
GRAND_NEED[0]['d_wfn'], '/', GRAND_NEED[0]['d_basis'], '',
GRAND_NEED[0]['d_energy'],
GRAND_NEED[0]['d_scheme'].__name__)
if do_corl:
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[1]['d_stage'],
GRAND_NEED[1]['d_wfn'], '/', GRAND_NEED[1]['d_basis'], '',
GRAND_NEED[1]['d_energy'] - GRAND_NEED[2]['d_energy'],
GRAND_NEED[1]['d_scheme'].__name__)
if do_delta:
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[3]['d_stage'],
GRAND_NEED[3]['d_wfn'] + ' - ' + GRAND_NEED[4]['d_wfn'], '/', GRAND_NEED[3]['d_basis'], '',
GRAND_NEED[3]['d_energy'] - GRAND_NEED[4]['d_energy'],
GRAND_NEED[3]['d_scheme'].__name__)
if do_delta2:
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[5]['d_stage'],
GRAND_NEED[5]['d_wfn'] + ' - ' + GRAND_NEED[6]['d_wfn'], '/', GRAND_NEED[5]['d_basis'], '',
GRAND_NEED[5]['d_energy'] - GRAND_NEED[6]['d_energy'],
GRAND_NEED[5]['d_scheme'].__name__)
# if do_delta3:
# tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[6]['d_stage'], GRAND_NEED[6]['d_wfn'] + ' - ' + GRAND_NEED[7]['d_wfn'],
# '/', GRAND_NEED[6]['d_basis'], '', GRAND_NEED[6]['d_energy'] - GRAND_NEED[7]['d_energy'], GRAND_NEED[6]['d_scheme'].__name__)
# if do_delta4:
# tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[8]['d_stage'], GRAND_NEED[8]['d_wfn'] + ' - ' + GRAND_NEED[9]['d_wfn'],
# '/', GRAND_NEED[8]['d_basis'], '', GRAND_NEED[8]['d_energy'] - GRAND_NEED[9]['d_energy'], GRAND_NEED[8]['d_scheme'].__name__)
# if do_delta5:
# tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % (GRAND_NEED[10]['d_stage'], GRAND_NEED[10]['d_wfn'] + ' - ' + GRAND_NEED[11]['d_wfn'],
# '/', GRAND_NEED[10]['d_basis'], '', GRAND_NEED[10]['d_energy'] - GRAND_NEED[11]['d_energy'], GRAND_NEED[10]['d_scheme'].__name__)
tables += """ %6s %20s %1s %-27s %2s %16.8f %-s\n""" % ('total', 'CBS', '', '', '', finalenergy, '')
tables += table_delimit
core.print_out(tables)
core.set_variable('CBS REFERENCE ENERGY', GRAND_NEED[0]['d_energy'])
core.set_variable('CBS CORRELATION ENERGY', finalenergy - GRAND_NEED[0]['d_energy'])
core.set_variable('CBS TOTAL ENERGY', finalenergy)
core.set_variable('CURRENT REFERENCE ENERGY', GRAND_NEED[0]['d_energy'])
core.set_variable('CURRENT CORRELATION ENERGY', finalenergy - GRAND_NEED[0]['d_energy'])
core.set_variable('CURRENT ENERGY', finalenergy)
core.set_variable('CBS NUMBER', Njobs)
# new skeleton wavefunction w/mol, highest-SCF basis (just to choose one), & not energy
basis = core.BasisSet.build(molecule, "ORBITAL", 'def2-svp')
wfn = core.Wavefunction(molecule, basis)
optstash.restore()
if ptype == 'energy':
finalquantity = finalenergy
elif ptype == 'gradient':
finalquantity = finalgradient
wfn.set_gradient(finalquantity)
if finalquantity.rows(0) < 20:
core.print_out('CURRENT GRADIENT')
finalquantity.print_out()
elif ptype == 'hessian':
finalquantity = finalhessian
wfn.set_hessian(finalquantity)
if finalquantity.rows(0) < 20:
core.print_out('CURRENT HESSIAN')
finalquantity.print_out()
if return_wfn:
return (finalquantity, wfn)
else:
return finalquantity
_lmh_labels = {1: ['HI'],
2: ['LO', 'HI'],
3: ['LO', 'MD', 'HI'],
4: ['LO', 'MD', 'M2', 'HI'],
5: ['LO', 'MD', 'M2', 'M3', 'HI']}
def _expand_scheme_orders(scheme, basisname, basiszeta, wfnname, natom):
"""Check that the length of *basiszeta* array matches the implied degree of
extrapolation in *scheme* name. Return a dictionary of same length as
basiszeta, with *basisname* and *basiszeta* distributed therein.
"""
Nxtpl = len(basiszeta)
if int(scheme.__name__.split('_')[-1]) != Nxtpl:
raise ValidationError("""Call to '%s' not valid with '%s' basis sets.""" % (scheme.__name__, len(basiszeta)))
f_fields = ['f_wfn', 'f_basis', 'f_zeta', 'f_energy', 'f_gradient', 'f_hessian']
NEED = {}
for idx in range(Nxtpl):
NEED[_lmh_labels[Nxtpl][idx]] = dict(zip(f_fields, [wfnname, basisname[idx], basiszeta[idx],
0.0,
core.Matrix(natom, 3),
core.Matrix(3 * natom, 3 * natom)]))
return NEED
def _contract_scheme_orders(needdict, datakey='f_energy'):
"""Prepared named arguments for extrapolation functions by
extracting zetas and values (which one determined by *datakey*) out
of *needdict* and returning a dictionary whose keys are contructed
from _lmh_labels.
"""
largs = {}
largs['functionname'] = needdict['HI']['f_wfn']
Nxtpl = len(needdict)
zlabels = _lmh_labels[Nxtpl] # e.g., ['LO', 'HI']
for zeta in range(Nxtpl):
zlab = zlabels[zeta] # e.g., LO
largs['z' + zlab] = needdict[zlab]['f_zeta']
largs['value' + zlab] = needdict[zlab][datakey]
return largs
## Aliases ##
complete_basis_set = cbs
def _cbs_wrapper_methods(**kwargs):
cbs_method_kwargs = ['scf_wfn', 'corl_wfn', 'delta_wfn']
cbs_method_kwargs += ['delta%d_wfn' % x for x in range(2, 6)]
cbs_methods = []
for method in cbs_method_kwargs:
if method in kwargs:
cbs_methods.append(kwargs[method])
return cbs_methods
def _parse_cbs_gufunc_string(method_name):
method_name_list = re.split( """\+(?=\s*[Dd]:)""", method_name)
if len(method_name_list) > 2:
raise ValidationError("CBS gufunc: Text parsing is only valid for a single delta, please use the CBS wrapper directly")
method_list = []
basis_list = []
for num, method_str in enumerate(method_name_list):
if (method_str.count("[") > 1) or (method_str.count("]") > 1):
raise ValidationError("""CBS gufunc: Too many brackets given! %s """ % method_str)
if method_str.count('/') != 1:
raise ValidationError("""CBS gufunc: All methods must specify a basis with '/'. %s""" % method_str)
if num > 0:
method_str = method_str.strip()
if method_str[:2].lower() != 'd:':
raise ValidationError("""CBS gufunc: Delta method must start with 'D:'.""")
else:
method_str = method_str[2:]
method, basis = method_str.split('/')
method_list.append(method)
basis_list.append(basis)
return method_list, basis_list
def _cbs_gufunc(func, total_method_name, **kwargs):
"""
Text based wrapper of the CBS function.
"""
# Catch kwarg issues for all methods
kwargs = p4util.kwargs_lower(kwargs)
return_wfn = kwargs.pop('return_wfn', False)
core.clean_variables()
ptype = kwargs.pop('ptype', None)
# Make sure the molecule the user provided is the active one
molecule = kwargs.pop('molecule', core.get_active_molecule())
molecule.update_geometry()
# Sanitize total_method_name
label = total_method_name
total_method_name = total_method_name.lower()
total_method_name = total_method_name.replace(' ', '')
# Split into components
method_list, basis_list = _parse_cbs_gufunc_string(total_method_name)
# Single energy call?
single_call = len(method_list) == 1
single_call &= '[' not in basis_list[0]
single_call &= ']' not in basis_list[0]
if single_call:
method_name = method_list[0]
basis = basis_list[0]
# Save some global variables so we can reset them later
optstash = p4util.OptionsState(['BASIS'])
core.set_global_option('BASIS', basis)
ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)
core.clean()
optstash.restore()
if return_wfn:
return (ptype_value, wfn)
else:
return ptype_value
# Drop out for props and freqs
if ptype in ["properties", "frequency"]:
raise ValidationError("%s: Cannot extrapolate or delta correct %s yet." % (ptype.title(), ptype))
# Catch kwarg issues for CBS methods only
user_dertype = kwargs.pop('dertype', None)
cbs_verbose = kwargs.pop('cbs_verbose', False)
# If we are not a single call, let CBS wrapper handle it!
cbs_kwargs = {}
cbs_kwargs['ptype'] = ptype
cbs_kwargs['return_wfn'] = True
cbs_kwargs['molecule'] = molecule
cbs_kwargs['verbose'] = cbs_verbose
# Find method and basis
if method_list[0] in ['scf', 'hf', 'c4-scf', 'c4-hf']:
cbs_kwargs['scf_wfn'] = method_list[0]
cbs_kwargs['scf_basis'] = basis_list[0]
if 'scf_scheme' in kwargs:
cbs_kwargs['scf_scheme'] = kwargs['scf_scheme']
else:
cbs_kwargs['corl_wfn'] = method_list[0]
cbs_kwargs['corl_basis'] = basis_list[0]
if 'corl_scheme' in kwargs:
cbs_kwargs['corl_scheme'] = kwargs['corl_scheme']
if len(method_list) > 1:
cbs_kwargs['delta_wfn'] = method_list[1]
cbs_kwargs['delta_basis'] = basis_list[1]
if 'delta_scheme' in kwargs:
cbs_kwargs['delta_scheme'] = kwargs['delta_scheme']
ptype_value, wfn = cbs(func, label, **cbs_kwargs)
if return_wfn:
return (ptype_value, wfn)
else:
return ptype_value
