Source code for psi4.driver.frac

#
# @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
#

__all__ = [
    "frac_nuke",
    "frac_traverse",
    "ip_fitting",
]

from typing import Callable, Dict, Union

from psi4 import core

from . import driver, p4util
from .p4util.exceptions import *


[docs] def frac_traverse(name: Union[str, Callable], **kwargs) -> Dict[float, float]: """Scan electron occupancy from +1 electron to -1. Parameters ---------- name DFT functional string name or function defining functional whose omega is to be optimized. molecule : :ref:`molecule <op_py_molecule>`, optional Target molecule (neutral) for which omega is to be tuned, if not last defined. cation_mult : Optional[int] Multiplicity of cation, if not neutral multiplicity + 1. anion_mult : Optional[int] Multiplicity of anion, if not neutral multiplicity + 1. frac_start : Optional[int] Iteration at which to start frac procedure when not reading previous guess. Defaults to 25. HOMO_occs : Optional[List] Occupations to step through for cation, by default `[1 - 0.1 * x for x in range(11)]`. LUMO_occs : Optional[List] Occupations to step through for anion, by default `[1 - 0.1 * x for x in range(11)]`. HOMO : Optional[int] Index of HOMO. LUMO : Optional[int] Index of LUMO. frac_diis : Optional[bool] Do use DIIS for non-1.0-occupied points? neutral_guess : Optional[bool] Do use neutral orbitals as guess for the anion? hf_guess: Optional[bool] Do use UHF guess before UKS? continuous_guess : Optional[bool] Do carry along guess rather than reguessing at each occupation? filename : Optional[str] Result filename, if not name of molecule. Returns ------- ~typing.Dict[float, float] Dictionary associating SCF energies with occupations. """ optstash = p4util.OptionsState( ['SCF', 'GUESS'], ['SCF', 'DF_INTS_IO'], ['SCF', 'REFERENCE'], ["SCF", "FRAC_START"], ["SCF", "FRAC_RENORMALIZE"], #["SCF", "FRAC_LOAD"], ["SCF", "FRAC_OCC"], ["SCF", "FRAC_VAL"], ["SCF", "FRAC_DIIS"]) kwargs = p4util.kwargs_lower(kwargs) # Make sure the molecule the user provided is the active one, and neutral molecule = kwargs.pop('molecule', core.get_active_molecule()) molecule.update_geometry() if molecule.molecular_charge() != 0: raise ValidationError("""frac_traverse requires neutral molecule to start.""") if molecule.schoenflies_symbol() != 'c1': core.print_out(""" Requested procedure `frac_traverse` does not make use of molecular symmetry: """ """further calculations in C1 point group.\n""") molecule = molecule.clone() molecule.reset_point_group('c1') molecule.update_geometry() charge0 = molecule.molecular_charge() mult0 = molecule.multiplicity() chargep = charge0 + 1 chargem = charge0 - 1 multp = kwargs.get('cation_mult', mult0 + 1) multm = kwargs.get('anion_mult', mult0 + 1) # By default, we start the frac procedure on the 25th iteration # when not reading a previous guess frac_start = kwargs.get('frac_start', 25) # By default, we occupy by tenths of electrons HOMO_occs = kwargs.get('HOMO_occs', [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]) LUMO_occs = kwargs.get('LUMO_occs', [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]) # By default, HOMO and LUMO are both in alpha Z = 0 for A in range(molecule.natom()): Z += molecule.Z(A) Z -= charge0 HOMO = kwargs.get('HOMO', (Z / 2 + 1 if (Z % 2) else Z / 2)) LUMO = kwargs.get('LUMO', HOMO + 1) # By default, DIIS in FRAC (1.0 occupation is always DIIS'd) frac_diis = kwargs.get('frac_diis', True) # By default, use the neutral orbitals as a guess for the anion neutral_guess = kwargs.get('neutral_guess', True) # By default, burn-in with UHF first, if UKS hf_guess = False if core.get_local_option('SCF', 'REFERENCE') == 'UKS': hf_guess = kwargs.get('hf_guess', True) # By default, re-guess at each N continuous_guess = kwargs.get('continuous_guess', False) # By default, drop the files to the molecule's name root = kwargs.get('filename', molecule.name()) traverse_filename = root + '.traverse.dat' # => Traverse <= # occs = [] energies = [] potentials = [] convs = [] # => Run the neutral for its orbitals, if requested <= # core.set_local_option("SCF", "DF_INTS_IO", "SAVE") old_guess = core.get_local_option("SCF", "GUESS") if (neutral_guess): if (hf_guess): core.set_local_option("SCF", "REFERENCE", "UHF") driver.energy('scf', dft_functional=name, molecule=molecule, **kwargs) core.set_local_option("SCF", "GUESS", "READ") core.set_local_option("SCF", "DF_INTS_IO", "LOAD") # => Run the anion first <= # molecule.set_molecular_charge(chargem) molecule.set_multiplicity(multm) # => Burn the anion in with hf, if requested <= # if hf_guess: core.set_local_option("SCF", "REFERENCE","UHF") driver.energy('scf', dft_functional=name, molecule=molecule, **kwargs) core.set_local_option("SCF", "REFERENCE", "UKS") core.set_local_option("SCF", "GUESS", "READ") core.set_local_option("SCF", "DF_INTS_IO", "SAVE") core.set_local_option("SCF", "FRAC_START", frac_start) core.set_local_option("SCF", "FRAC_RENORMALIZE", True) # NYI core.set_local_option("SCF", "FRAC_LOAD", False) for occ in LUMO_occs: core.set_local_option("SCF", "FRAC_OCC", [LUMO]) core.set_local_option("SCF", "FRAC_VAL", [occ]) E, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, **kwargs) C = 1 if E == 0.0: E = core.variable('SCF ITERATION ENERGY') C = 0 if LUMO > 0: eps = wfn.epsilon_a() potentials.append(eps.get(int(LUMO) - 1)) else: eps = wfn.epsilon_b() potentials.append(eps.get(-int(LUMO) - 1)) occs.append(occ) energies.append(E) convs.append(C) core.set_local_option("SCF", "FRAC_START", 2) #core.set_local_option("SCF", "FRAC_LOAD", True) core.set_local_option("SCF", "GUESS", "READ") core.set_local_option("SCF", "FRAC_DIIS", frac_diis) core.set_local_option("SCF", "DF_INTS_IO", "LOAD") # => Run the neutral next <= # molecule.set_molecular_charge(charge0) molecule.set_multiplicity(mult0) # Burn the neutral in with hf, if requested <= # if not continuous_guess: core.set_local_option("SCF", "GUESS", old_guess) if hf_guess: core.set_local_option("SCF", "FRAC_START", 0) core.set_local_option("SCF", "REFERENCE", "UHF") driver.energy('scf', dft_functional=name, molecule=molecule, **kwargs) core.set_local_option("SCF", "REFERENCE", "UKS") core.set_local_option("SCF", "GUESS", "READ") # NYI core.set_local_option("SCF", "FRAC_LOAD", False) core.set_local_option("SCF", "FRAC_START", frac_start) core.set_local_option("SCF", "FRAC_RENORMALIZE", True) for occ in HOMO_occs: core.set_local_option("SCF", "FRAC_OCC", [HOMO]) core.set_local_option("SCF", "FRAC_VAL", [occ]) E, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, **kwargs) C = 1 if E == 0.0: E = core.variable('SCF ITERATION ENERGY') C = 0 if LUMO > 0: eps = wfn.epsilon_a() potentials.append(eps.get(int(HOMO) - 1)) else: eps = wfn.epsilon_b() potentials.append(eps.get(-int(HOMO) - 1)) occs.append(occ - 1.0) energies.append(E) convs.append(C) core.set_local_option("SCF", "FRAC_START", 2) # NYI core.set_local_option("SCF", "FRAC_LOAD", True) core.set_local_option("SCF", "GUESS", "READ") core.set_local_option("SCF", "FRAC_DIIS", frac_diis) core.set_local_option("SCF", "DF_INTS_IO", "LOAD") # => Print the results out <= # E = {} core.print_out("""\n ==> Fractional Occupation Traverse Results <==\n\n""") core.print_out(""" %-11s %-24s %-24s %11s\n""" % ('N', 'Energy', 'HOMO Energy', 'Converged')) for k in range(len(occs)): core.print_out(""" %11.3E %24.16E %24.16E %11d\n""" % (occs[k], energies[k], potentials[k], convs[k])) E[occs[k]] = energies[k] core.print_out(""" You trying to be a hero Watkins? Just trying to kill some bugs sir! -Starship Troopers""") # Drop the files out with open(traverse_filename, 'w') as fh: fh.write(""" %-11s %-24s %-24s %11s\n""" % ('N', 'Energy', 'HOMO Energy', 'Converged')) for k in range(len(occs)): fh.write(""" %11.3E %24.16E %24.16E %11d\n""" % (occs[k], energies[k], potentials[k], convs[k])) optstash.restore() return E
[docs] def frac_nuke(name: Union[str, Callable], **kwargs) -> Dict[float, float]: """Pull all the electrons out, one at a time""" optstash = p4util.OptionsState( ['SCF', 'GUESS'], ['SCF', 'DF_INTS_IO'], ["SCF", "FRAC_START"], ["SCF", "FRAC_RENORMALIZE"], # NYI ["SCF", "FRAC_LOAD"], ["SCF", "FRAC_OCC"], ["SCF", "FRAC_VAL"], ["SCF", "FRAC_DIIS"]) kwargs = p4util.kwargs_lower(kwargs) # Make sure the molecule the user provided is the active one, and neutral molecule = kwargs.pop('molecule', core.get_active_molecule()) molecule.update_geometry() if molecule.molecular_charge() != 0: raise ValidationError("""frac_nuke requires neutral molecule to start.""") if molecule.schoenflies_symbol() != 'c1': core.print_out(""" Requested procedure `frac_nuke` does not make use of molecular symmetry: """ """further calculations in C1 point group.\n""") molecule = molecule.clone() molecule.reset_point_group('c1') molecule.update_geometry() charge0 = molecule.molecular_charge() mult0 = molecule.multiplicity() # By default, we start the frac procedure on the 25th iteration # when not reading a previous guess frac_start = kwargs.get('frac_start', 25) # By default, we occupy by tenths of electrons foccs = kwargs.get('foccs', [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]) # By default, HOMO and LUMO are both in alpha N = 0 for A in range(molecule.natom()): N += molecule.Z(A) N -= charge0 N = int(N) Nb = int((N - mult0 + 1) / 2) Na = int(N - Nb) charge = charge0 mult = mult0 # By default, nuke all the electrons Nmin = 0 if 'nmax' in kwargs: Nmin = N - int(kwargs['nmax']) # By default, DIIS in FRAC (1.0 occupation is always DIIS'd) frac_diis = kwargs.get('frac_diis', True) # By default, drop the files to the molecule's name root = kwargs.get('filename', molecule.name()) traverse_filename = root + '.traverse.dat' stats_filename = root + '.stats.dat' # => Traverse <= # core.set_local_option("SCF", "DF_INTS_IO", "SAVE") Ns = [] energies = [] potentials = [] convs = [] stats = [] # Run one SCF to burn things in E, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, **kwargs) # Determine HOMO eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() eps_a.print_out() if Na == Nb: HOMO = -Nb elif Nb == 0: HOMO = Na else: E_a = eps_a.get(int(Na - 1)) E_b = eps_b.get(int(Nb - 1)) if E_a >= E_b: HOMO = Na else: HOMO = -Nb stats.append(""" %6d %6d %6d %6d %6d %6d\n""" % (N, Na, Nb, charge, mult, HOMO)) if HOMO > 0: Na -= 1 else: Nb -= 1 charge += 1 mult = Na - Nb + 1 core.set_local_option("SCF", "DF_INTS_IO", "LOAD") core.set_local_option("SCF", "FRAC_START", frac_start) core.set_local_option("SCF", "FRAC_RENORMALIZE", True) # Nuke 'em Rico! for Nintegral in range(N, Nmin, -1): # Nuke the current HOMO for occ in foccs: core.set_local_option("SCF", "FRAC_OCC", [HOMO]) core.set_local_option("SCF", "FRAC_VAL", [occ]) E, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, **kwargs) C = 1 if E == 0.0: E = core.variable('SCF ITERATION ENERGY') C = 0 if HOMO > 0: eps = wfn.epsilon_a() potentials.append(eps.np[HOMO - 1]) else: eps = wfn.epsilon_b() potentials.append(eps.np[-HOMO - 1]) Ns.append(Nintegral + occ - 1.0) energies.append(E) convs.append(C) core.set_local_option("SCF", "FRAC_START", 2) # NYI core.set_local_option("SCF", "FRAC_LOAD", True) core.set_local_option("SCF", "FRAC_DIIS", frac_diis) core.set_local_option("SCF", "GUESS", "READ") # Set the next charge/mult molecule.set_molecular_charge(charge) molecule.set_multiplicity(mult) # Determine HOMO print('DGAS: What ref should this point to?') eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() if Na == Nb: HOMO = -Nb elif Nb == 0: HOMO = Na else: E_a = eps_a.np[int(Na - 1)] E_b = eps_b.np[int(Nb - 1)] if E_a >= E_b: HOMO = Na else: HOMO = -Nb stats.append(""" %6d %6d %6d %6d %6d %6d\n""" % (Nintegral-1, Na, Nb, charge, mult, HOMO)) if HOMO > 0: Na -= 1 else: Nb -= 1 charge += 1 mult = Na - Nb + 1 core.set_local_option("SCF", "DF_INTS_IO", "NONE") # => Print the results out <= # E = {} core.print_out("""\n ==> Fractional Occupation Nuke Results <==\n\n""") core.print_out(""" %-11s %-24s %-24s %11s\n""" % ('N', 'Energy', 'HOMO Energy', 'Converged')) for k in range(len(Ns)): core.print_out(""" %11.3E %24.16E %24.16E %11d\n""" % (Ns[k], energies[k], potentials[k], convs[k])) E[Ns[k]] = energies[k] core.print_out('\n') core.print_out(""" %6s %6s %6s %6s %6s %6s\n""" % ('N', 'Na', 'Nb', 'Charge', 'Mult', 'HOMO')) for line in stats: core.print_out(line) core.print_out('\n "You shoot a nuke down a bug hole, you got a lot of dead bugs"\n') core.print_out(' -Starship Troopers\n') # Drop the files out with open(traverse_filename, 'w') as fh: fh.write(""" %-11s %-24s %-24s %11s\n""" % ('N', 'Energy', 'HOMO Energy', 'Converged')) for k in range(len(Ns)): fh.write(""" %11.3E %24.16E %24.16E %11d\n""" % (Ns[k], energies[k], potentials[k], convs[k])) with open(stats_filename, 'w') as fh: fh.write(""" %6s %6s %6s %6s %6s %6s\n""" % ('N', 'Na', 'Nb', 'Charge', 'Mult', 'HOMO')) for line in stats: fh.write(line) optstash.restore() return E
[docs] def ip_fitting(name: Union[str, Callable], omega_l: float = 0.05, omega_r: float = 2.5, omega_convergence: float = 1.0e-3, maxiter: int = 20, **kwargs) -> float: """Optimize DFT omega parameter for molecular system. Parameters ---------- name DFT functional string name or function defining functional whose omega is to be optimized. omega_l Minimum omega to be considered during fitting. omega_r Maximum omega to be considered during fitting. molecule : :ref:`molecule <op_py_molecule>`, optional Target molecule (neutral) for which omega is to be tuned, if not last defined. omega_convergence Threshold below which to consider omega converged. (formerly omega_tolerance) maxiter Maximum number of iterations towards omega convergence. Returns ------- float Optimal omega parameter. """ optstash = p4util.OptionsState( ['SCF', 'REFERENCE'], ['SCF', 'GUESS'], ['SCF', 'DF_INTS_IO'], ['SCF', 'DFT_OMEGA'], ['DOCC'], ['SOCC']) kwargs = p4util.kwargs_lower(kwargs) # By default, do not read previous 180 orbitals file read = False read180 = '' if 'read' in kwargs: read = True read180 = kwargs['read'] if core.get_option('SCF', 'REFERENCE') != 'UKS': core.print_out(""" Requested procedure `ip_fitting` runs further calculations with UKS reference.\n""") core.set_local_option('SCF', 'REFERENCE', 'UKS') # Make sure the molecule the user provided is the active one, and neutral molecule = kwargs.pop('molecule', core.get_active_molecule()) molecule.update_geometry() if molecule.molecular_charge() != 0: raise ValidationError("""IP Fitting requires neutral molecule to start.""") if molecule.schoenflies_symbol() != 'c1': core.print_out(""" Requested procedure `ip_fitting` does not make use of molecular symmetry: """ """further calculations in C1 point group.\n""") molecule = molecule.clone() molecule.reset_point_group('c1') molecule.update_geometry() charge0 = molecule.molecular_charge() mult0 = molecule.multiplicity() # How many electrons are there? N = 0 for A in range(molecule.natom()): N += molecule.Z(A) N -= charge0 N = int(N) Nb = int((N - mult0 + 1) / 2) Na = int(N - Nb) # Work in the ot namespace for this procedure core.IO.set_default_namespace("ot") # Burn in to determine orbital eigenvalues if read: core.set_local_option("SCF", "GUESS", "READ") copy_file_to_scratch(read180, 'psi', 'ot', 180) core.set_local_option("SCF", "DF_INTS_IO", "SAVE") E, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, banner='IP Fitting SCF: Burn-in', **kwargs) core.set_local_option("SCF", "DF_INTS_IO", "LOAD") if not wfn.functional().is_x_lrc(): raise ValidationError("""Not sensible to optimize omega for non-long-range-correction functional.""") # Determine HOMO, to determine mult1 eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() if Na == Nb: HOMO = -Nb elif Nb == 0: HOMO = Na else: E_a = eps_a.np[int(Na - 1)] E_b = eps_b.np[int(Nb - 1)] if E_a >= E_b: HOMO = Na else: HOMO = -Nb Na1 = Na Nb1 = Nb if HOMO > 0: Na1 -= 1 else: Nb1 -= 1 charge1 = charge0 + 1 mult1 = Na1 - Nb1 + 1 omegas = [] E0s = [] E1s = [] kIPs = [] IPs = [] types = [] # Right endpoint core.set_local_option('SCF', 'DFT_OMEGA', omega_r) # Neutral if read: core.set_local_option("SCF", "GUESS", "READ") p4util.copy_file_to_scratch(read180, 'psi', 'ot', 180) molecule.set_molecular_charge(charge0) molecule.set_multiplicity(mult0) E0r, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, banner='IP Fitting SCF: Neutral, Right Endpoint', **kwargs) eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() if Nb == 0: E_HOMO = eps_a.np[int(Na - 1)] else: E_a = eps_a.np[int(Na - 1)] E_b = eps_b.np[int(Nb - 1)] E_HOMO = max(E_a, E_b) E_HOMOr = E_HOMO core.IO.change_file_namespace(180, "ot", "neutral") # Cation if read: core.set_local_option("SCF", "GUESS", "READ") p4util.copy_file_to_scratch(read180, 'psi', 'ot', 180) molecule.set_molecular_charge(charge1) molecule.set_multiplicity(mult1) E1r = driver.energy('scf', dft_functional=name, molecule=molecule, banner='IP Fitting SCF: Cation, Right Endpoint', **kwargs) core.IO.change_file_namespace(180, "ot", "cation") IPr = E1r - E0r kIPr = -E_HOMOr delta_r = IPr - kIPr if IPr > kIPr: raise ValidationError("""\n***IP Fitting Error: Right Omega limit should have kIP > IP: {} !> {}""".format(kIPr, IPr)) omegas.append(omega_r) types.append('Right Limit') E0s.append(E0r) E1s.append(E1r) IPs.append(IPr) kIPs.append(kIPr) # Use previous orbitals from here out core.set_local_option("SCF", "GUESS", "READ") # Left endpoint core.set_local_option('SCF', 'DFT_OMEGA', omega_l) # Neutral core.IO.change_file_namespace(180, "neutral", "ot") molecule.set_molecular_charge(charge0) molecule.set_multiplicity(mult0) core.set_global_option("DOCC", [Nb]) core.set_global_option("SOCC", [Na - Nb]) E0l, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, banner='IP Fitting SCF: Neutral, Left Endpoint', **kwargs) eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() if Nb == 0: E_HOMO = eps_a.np[int(Na - 1)] else: E_a = eps_a.np[int(Na - 1)] E_b = eps_b.np[int(Nb - 1)] E_HOMO = max(E_a, E_b) E_HOMOl = E_HOMO core.IO.change_file_namespace(180, "ot", "neutral") # Cation core.IO.change_file_namespace(180, "cation", "ot") molecule.set_molecular_charge(charge1) molecule.set_multiplicity(mult1) core.set_global_option("DOCC", [Nb1]) core.set_global_option("SOCC", [Na1 - Nb1]) E1l = driver.energy('scf', dft_functional=name, molecule=molecule, banner='IP Fitting SCF: Cation, Left Endpoint', **kwargs) core.IO.change_file_namespace(180, "ot", "cation") IPl = E1l - E0l kIPl = -E_HOMOl delta_l = IPl - kIPl if IPl < kIPl: raise ValidationError("""\n***IP Fitting Error: Left Omega limit should have kIP < IP: {} !< {}""".format(kIPl, IPl)) omegas.append(omega_l) types.append('Left Limit') E0s.append(E0l) E1s.append(E1l) IPs.append(IPl) kIPs.append(kIPl) converged = False repeat_l = 0 repeat_r = 0 for step in range(maxiter): # Regula Falsi (modified) if repeat_l > 1: delta_l /= 2.0 if repeat_r > 1: delta_r /= 2.0 omega = - (omega_r - omega_l) / (delta_r - delta_l) * delta_l + omega_l core.set_local_option('SCF', 'DFT_OMEGA', omega) # Neutral core.IO.change_file_namespace(180, "neutral", "ot") molecule.set_molecular_charge(charge0) molecule.set_multiplicity(mult0) core.set_global_option("DOCC", [Nb]) core.set_global_option("SOCC", [Na - Nb]) E0, wfn = driver.energy('scf', dft_functional=name, return_wfn=True, molecule=molecule, banner='IP Fitting SCF: Neutral, Omega = {:11.3E}'.format(omega), **kwargs) eps_a = wfn.epsilon_a() eps_b = wfn.epsilon_b() if Nb == 0: E_HOMO = eps_a.np[int(Na - 1)] else: E_a = eps_a.np[int(Na - 1)] E_b = eps_b.np[int(Nb - 1)] E_HOMO = max(E_a, E_b) core.IO.change_file_namespace(180, "ot", "neutral") # Cation core.IO.change_file_namespace(180, "cation", "ot") molecule.set_molecular_charge(charge1) molecule.set_multiplicity(mult1) core.set_global_option("DOCC", [Nb1]) core.set_global_option("SOCC", [Na1 - Nb1]) E1 = driver.energy('scf', dft_functional=name, molecule=molecule, banner='IP Fitting SCF: Cation, Omega = {:11.3E}'.format(omega), **kwargs) core.IO.change_file_namespace(180, "ot", "cation") IP = E1 - E0 kIP = -E_HOMO delta = IP - kIP if kIP > IP: omega_r = omega E0r = E0 E1r = E1 IPr = IP kIPr = kIP delta_r = delta repeat_r = 0 repeat_l += 1 else: omega_l = omega E0l = E0 E1l = E1 IPl = IP kIPl = kIP delta_l = delta repeat_l = 0 repeat_r += 1 omegas.append(omega) types.append('Regula-Falsi') E0s.append(E0) E1s.append(E1) IPs.append(IP) kIPs.append(kIP) # Termination if abs(omega_l - omega_r) < omega_convergence: converged = True break core.IO.set_default_namespace("") core.print_out("""\n ==> IP Fitting Results <==\n\n""") core.print_out(""" => Occupation Determination <= \n\n""") core.print_out(""" %6s %6s %6s %6s %6s %6s\n""" % ('N', 'Na', 'Nb', 'Charge', 'Mult', 'HOMO')) core.print_out(""" Neutral: %6d %6d %6d %6d %6d %6d\n""" % (N, Na, Nb, charge0, mult0, HOMO)) core.print_out(""" Cation: %6d %6d %6d %6d %6d\n\n""" % (N - 1, Na1, Nb1, charge1, mult1)) core.print_out(""" => Regula Falsi Iterations <=\n\n""") core.print_out(""" %3s %11s %14s %14s %14s %s\n""" % ('N','Omega','IP','kIP','Delta','Type')) for k in range(len(omegas)): core.print_out(""" %3d %11.3E %14.6E %14.6E %14.6E %s\n""" % (k + 1, omegas[k], IPs[k], kIPs[k], IPs[k] - kIPs[k], types[k])) optstash.restore() if converged: core.print_out("""\n IP Fitting Converged\n""") core.print_out(""" Final omega = %14.6E\n""" % ((omega_l + omega_r) / 2)) core.print_out("""\n "M,I. does the dying. Fleet just does the flying."\n""") core.print_out(""" -Starship Troopers\n""") else: raise ConvergenceError("""IP Fitting """, step + 1) return ((omega_l + omega_r) / 2)