diff --git a/pysisyphus/calculators/Calculator.py b/pysisyphus/calculators/Calculator.py index 212a56d4cc..e8157ca7c6 100644 --- a/pysisyphus/calculators/Calculator.py +++ b/pysisyphus/calculators/Calculator.py @@ -17,6 +17,7 @@ from pysisyphus.constants import BOHR2ANG from pysisyphus.helpers import geom_loader from pysisyphus.linalg import finite_difference_hessian +from pysisyphus.wavefunction import Wavefunction KeepKind = Enum("KeepKind", ["ALL", "LATEST", "NONE"]) @@ -37,7 +38,6 @@ def __post_init__(self): class Calculator: - conf_key = None _set_plans = [] @@ -238,6 +238,12 @@ def restore_org_hessian(self): self.get_hessian = self._org_get_hessian self.hessian_kind = HessKind["ORG"] + def load_wavefunction_from_file(self, fn, **kwargs): + return Wavefunction.from_file(fn, **kwargs) + + def get_stored_wavefunction(self, **kwargs): + raise Exception("Not implemented!") + def make_fn(self, name, counter=None, return_str=False): """Make a full filename. diff --git a/pysisyphus/calculators/ORCA.py b/pysisyphus/calculators/ORCA.py index 9477432c34..cf920701b6 100644 --- a/pysisyphus/calculators/ORCA.py +++ b/pysisyphus/calculators/ORCA.py @@ -367,7 +367,6 @@ def get_exc_ens_fosc(wf_fn, cis_fn, log_fn): class ORCA(OverlapCalculator): - conf_key = "orca" _set_plans = ( "gbw", @@ -375,6 +374,7 @@ class ORCA(OverlapCalculator): "cis", "densities", ("molden", "mwfn_wf"), + "json", ) def __init__( @@ -384,7 +384,8 @@ def __init__( gbw=None, do_stable=False, numfreq=False, - json_dump=True, + json_dump=None, + wavefunction_dump=True, **kwargs, ): """ORCA calculator. @@ -412,7 +413,9 @@ def __init__( before every calculation. numfreq : bool, optional Use numerical frequencies instead of analytical ones. - json_dump : bool, optional + json_dump : bool, optional, deprecated + Use 'wavefunction_dump' instead. + wavefunction_dump : bool, optional Whether to dump the wavefunction to JSON via orca_2json. The JSON can become very large in calculations comprising many basis functions. """ @@ -423,7 +426,13 @@ def __init__( self.gbw = gbw self.do_stable = bool(do_stable) self.freq_keyword = "numfreq" if numfreq else "freq" - self.json_dump = bool(json_dump) + if json_dump is not None: + warnings.warn( + "Use of 'json_dump' is deprecated! Use 'wavefunction_dump' instead!", + DeprecationWarning, + ) + wavefunction_dump = json_dump + self.wavefunction_dump = bool(wavefunction_dump) assert ("pal" not in keywords) and ("nprocs" not in blocks), ( "PALn/nprocs not " "allowed! Use 'pal: n' in the 'calc' section instead." @@ -610,11 +619,14 @@ def get_hessian(self, atoms, coords, **prepare_kwargs): inp = self.prepare_input(atoms, coords, calc_type, **prepare_kwargs) results = self.run(inp, calc="hessian") - # results = self.store_and_track( - # results, self.get_hessian, atoms, coords, **prepare_kwargs - # ) + results = self.store_and_track( + results, self.get_hessian, atoms, coords, **prepare_kwargs + ) return results + def get_stored_wavefunction(self, **kwargs): + return self.load_wavefunction_from_file(self.json, **kwargs) + def run_calculation(self, atoms, coords, **prepare_kwargs): """Basically some kind of dummy method that can be called to execute ORCA with the stored cmd of this calculator.""" @@ -634,7 +646,7 @@ def run_after(self, path): self.popen(cmd, cwd=path) shutil.copy(path / "orca.molden.input", path / "orca.molden") - if self.json_dump: + if self.wavefunction_dump: # Will silently fail with ECPs cmd = "orca_2json orca" proc = self.popen(cmd, cwd=path) diff --git a/pysisyphus/calculators/XTB.py b/pysisyphus/calculators/XTB.py index 7e59501827..4af1fd165e 100644 --- a/pysisyphus/calculators/XTB.py +++ b/pysisyphus/calculators/XTB.py @@ -21,13 +21,8 @@ class XTB(Calculator): - conf_key = "xtb" - _set_plans = ( - "charges", - "json", - "xtbrestart", - ) + _set_plans = ("charges", "json", "xtbrestart", "molden") def __init__( self, @@ -42,6 +37,7 @@ def __init__( topo=None, topo_update=None, quiet=False, + wavefunction_dump=False, **kwargs, ): """XTB calculator. @@ -74,6 +70,8 @@ def __init__( Mememory per core in MB. quiet : bool, optional Suppress creation of log files. + wavefunction_dump : bool + Whether to dump a molden file. """ super().__init__(**kwargs) @@ -92,6 +90,7 @@ def __init__( self.topo = topo self.topo_update = topo_update self.quiet = quiet + self.wavefunction_dump = wavefunction_dump self.topo_used = 0 self.xtbrestart = None @@ -109,12 +108,12 @@ def __init__( "xtbopt.xyz", "g98.out", "xtb.trj", - # "json:xtbout.json", "charges:charges", "xcontrol", + "molden:molden.input", ) if self.restart: - self.to_keep += ("xtbrestart", ) + self.to_keep += ("xtbrestart",) if self.quiet: self.to_keep = () @@ -193,6 +192,8 @@ def prepare_add_args(self, xcontrol=None): elif self.alpb: alpb = f"--alpb {self.alpb}".split() add_args = add_args + alpb + if self.wavefunction_dump: + add_args = add_args + ["--molden"] # Select parametrization gfn = ["--gfnff"] if self.gfn == "ff" else f"--gfn {self.gfn}".split() add_args = add_args + gfn @@ -238,6 +239,9 @@ def get_hessian(self, atoms, coords, **prepare_kwargs): results = self.run(inp, **kwargs) return results + def get_stored_wavefunction(self, **kwargs): + return self.load_wavefunction_from_file(self.molden, xtb_nuc_charges=True) + def run_calculation(self, atoms, coords, **prepare_kwargs): self.prepare_input(atoms, coords, "calculation", **prepare_kwargs) inp = self.prepare_coords(atoms, coords) @@ -377,7 +381,7 @@ def parse_hessian(self, path): with open(path / "hessian") as handle: text = handle.read() hessian = np.array(text.split()[1:], dtype=float) - coord_num = int(hessian.size ** 0.5) + coord_num = int(hessian.size**0.5) hessian = hessian.reshape(coord_num, coord_num) energy = self.parse_energy(path) results = { diff --git a/pysisyphus/io/molden.py b/pysisyphus/io/molden.py index ad8550f3e6..ace668de5e 100644 --- a/pysisyphus/io/molden.py +++ b/pysisyphus/io/molden.py @@ -2,12 +2,19 @@ import numpy as np import pyparsing as pp +from scipy.special import gamma -from pysisyphus.elem_data import ATOMIC_NUMBERS +from pysisyphus.elem_data import ATOMIC_NUMBERS, nuc_charges_for_atoms from pysisyphus.Geometry import Geometry from pysisyphus.io.xyz import parse_xyz from pysisyphus.helpers_pure import file_or_str -from pysisyphus.wavefunction import MoldenShells, Shell, Wavefunction +from pysisyphus.wavefunction import ( + get_l, + MoldenShells, + ORCAMoldenShells, + Shell, + Wavefunction, +) from pysisyphus.wavefunction.helpers import BFType @@ -182,6 +189,56 @@ def pgto_exps_coeffs(s, loc, toks): return as_dict +def get_xtb_nuc_charges(atoms, as_ecp_electrons=False): + """Modified nuclear charges w/o core electrons. + + Adapated from Multiwfn 3.8 + """ + + org_nuc_charges = nuc_charges_for_atoms(atoms) + ecp_electrons = np.zeros_like(org_nuc_charges, dtype=int) + for i, Z in enumerate(org_nuc_charges): + # H, He + if 1 <= Z <= 2: + ecpel = 0 + # Li - Ne + elif 3 <= Z <= 10: + ecpel = 2 + # Na - Ar + elif 11 <= Z <= 18: + ecpel = 10 + # K - Cu + elif 19 <= Z <= 29: + ecpel = 18 + # Zn - Kr + elif 30 <= Z <= 36: + ecpel = 28 + # Rb - Ag + elif 37 <= Z <= 47: + ecpel = 36 + # Cd - Xe + elif 48 <= Z <= 54: + ecpel = 46 + # Cs - Ba + elif 55 <= Z <= 56: + ecpel = 54 + # La - Lu, crazy + elif 57 <= Z <= 71: + ecpel = Z - 3 # Resulting charge is always 3 + # Hf - Au + elif 72 <= Z <= 79: + ecpel = 68 + # Hg - Rn + elif 80 <= Z <= 86: + ecpel = 78 + ecp_electrons[i] = ecpel + if as_ecp_electrons: + result = ecp_electrons + else: + result = org_nuc_charges - ecp_electrons + return result + + def parse_molden_atoms(data): atoms = list() coords = list() @@ -226,46 +283,147 @@ def shells_from_molden(text): return shells +def radial_integral(l, exponent): + """ + Integrates + (r r**l * exp(-exponent * r**2))**2 dr from r=0 to r=oo + as described in the SI of the JANPA paper [1] (see top of page 8, + second integral in the square root. + + In my opinion, the integrals lacks a factor 'r'. Below, some sympy code + can be found to solve this integral (including 1*r). + + import sympy as sym + r, z = sym.symbols("r z", positive=True) + l = sym.symbols("l", integer=True, positive=True) + sym.integrate((r * r**l * sym.exp(-z*r**2))**2, (r, 0, sym.oo)) + + The 'solved' integral on page 8 is correct again. + + ∞ + ⌠ + ⎮ 2 + ⎮ 2 2⋅l -2⋅r ⋅z + ⎮ r ⋅r ⋅ℯ dr = (2*z)**(-l - 1/2)*gamma(l + 3/2)/(4*z) + ⌡ + 0 + """ + return (2 * exponent) ** (-l - 1 / 2) * gamma(l + 3 / 2) / (4 * exponent) + + +@file_or_str(".molden", ".input") +def shells_from_orca_molden(text): + molden_shells = shells_from_molden(text) + + dividers = { + 0: 1, + 1: 1, + 2: 3, + 3: 15, + 4: 35, + } + + def fix_contr_coeffs(l, coeffs, exponents): + """Fix contraction coefficients. Based on equations found in the SI + of the JANPA paper [1].""" + l = get_l(l) + divider = dividers[l] + rad_ints = radial_integral(l, exponents) + norms2 = rad_ints * 4 * np.pi / (2 * l + 1) / divider + norms = np.sqrt(norms2) + normed_coeffs = coeffs * norms + return normed_coeffs + + _shells = list() + for shell in molden_shells.shells: + L, center, _, exps, *_ = shell.as_tuple() + fixed_coeffs = fix_contr_coeffs(L, shell.coeffs_org, exps) + fixed_shell = Shell( + L=L, + center=center, + coeffs=fixed_coeffs, + exps=exps, + center_ind=shell.center_ind, + atomic_num=shell.atomic_num, + ) + _shells.append(fixed_shell) + shells = ORCAMoldenShells(_shells) + return shells + + @file_or_str(".molden", ".input") -def wavefunction_from_molden(text, charge=None, shells_func=None, **wf_kwargs): +def wavefunction_from_molden( + text, + charge=None, + orca_contr_renorm=False, + xtb_nuc_charges: bool = False, + **wf_kwargs, +): """Construct Wavefunction object from .molden file. - shells_func is used for ORCA, to modify the contraction coefficients. + orca_contr_renorm fixes contraction coefficients, as required for + ORCA and XTB. """ data = parse_molden(text) atoms, coords, nuc_charges = parse_molden_atoms(data) nuc_charge = sum(nuc_charges) + # XTB only describes valence electrons. The missing electrons will be + # treated as ECP electrons. + if xtb_nuc_charges: + ecp_electrons = get_xtb_nuc_charges(atoms, as_ecp_electrons=True) + wf_kwargs["ecp_electrons"] = ecp_electrons + nuc_charge = nuc_charge - sum(ecp_electrons) spins = list() occ_a = 0.0 occ_b = 0.0 + occ_a_unpaired = 0.0 + occ_b_unpaired = 0.0 Ca = list() Cb = list() for mo in data["mos"]: spin = mo["spin"].lower() occ = mo["occup"] + assert occ >= 0.0 + is_unpaired = 0.0 < occ < 2.0 spins.append(spin) coeffs = mo["coeffs"] if spin == "alpha": occ_a += occ Ca.append(coeffs) + if is_unpaired: + occ_a_unpaired += occ elif spin == "beta": occ_b += occ Cb.append(coeffs) + if is_unpaired: + occ_b_unpaired += occ else: raise Exception(f"Spin can only be 'alpha' or 'beta', but got '{spin}'!") assert occ_a.is_integer assert occ_b.is_integer occ_a = int(occ_a) occ_b = int(occ_b) + occ_a_unpaired = int(occ_a_unpaired) + occ_b_unpaired = int(occ_b_unpaired) + unpaired_electrons = occ_a_unpaired or occ_b_unpaired + unrestricted = set(spins) == {"Alpha", "Beta"} or unpaired_electrons + restricted = not unrestricted # MOs must be in columns Ca = np.array(Ca, dtype=float).T Cb = np.array(Cb, dtype=float).T # Restricted calculation - if Cb.size == 0: + if restricted: Cb = Ca.copy() occ_a = occ_b = occ_a // 2 + # The second case is hit w/ XTB molden files and restriced-open-shell calculatons + elif unrestricted and Cb.size == 0: + assert occ_a >= occ_b + common_electrons = (occ_a - occ_a_unpaired) // 2 + Cb = Ca.copy() + occ_a -= common_electrons + occ_b = common_electrons C = np.stack((Ca, Cb)) molden_charge = nuc_charge - (occ_a + occ_b) @@ -273,10 +431,12 @@ def wavefunction_from_molden(text, charge=None, shells_func=None, **wf_kwargs): charge = molden_charge # Multiplicity = (2S + 1), S = number of unpaired elecs. * 0.5 mult = (occ_a - occ_b) + 1 - unrestricted = set(spins) == {"Alpha", "Beta"} - if shells_func is None: + if orca_contr_renorm: + shells_func = shells_from_orca_molden + else: shells_func = shells_from_molden + shells = shells_func(text) return Wavefunction( diff --git a/pysisyphus/io/orca.py b/pysisyphus/io/orca.py index dcf7878fbd..c3cec94604 100644 --- a/pysisyphus/io/orca.py +++ b/pysisyphus/io/orca.py @@ -6,7 +6,6 @@ import json import numpy as np -from scipy.special import gamma from pysisyphus.constants import BOHR2ANG from pysisyphus.helpers_pure import file_or_str @@ -15,7 +14,6 @@ get_l, Shell, ORCAShells, - ORCAMoldenShells, Wavefunction, ) from pysisyphus.wavefunction.helpers import BFType @@ -132,77 +130,8 @@ def get_occ_and_mo_coeffs(mos): ) -def radial_integral(l, exponent): - """ - Integrates - (r r**l * exp(-exponent * r**2))**2 dr from r=0 to r=oo - as described in the SI of the JANPA paper [1] (see top of page 8, - second integral in the square root. - - In my opinion, the integrals lacks a factor 'r'. Below, some sympy code - can be found to solve this integral (including 1*r). - - import sympy as sym - r, z = sym.symbols("r z", positive=True) - l = sym.symbols("l", integer=True, positive=True) - sym.integrate((r * r**l * sym.exp(-z*r**2))**2, (r, 0, sym.oo)) - - The 'solved' integral on page 8 is correct again. - - ∞ - ⌠ - ⎮ 2 - ⎮ 2 2⋅l -2⋅r ⋅z - ⎮ r ⋅r ⋅ℯ dr = (2*z)**(-l - 1/2)*gamma(l + 3/2)/(4*z) - ⌡ - 0 - """ - return (2 * exponent) ** (-l - 1 / 2) * gamma(l + 3 / 2) / (4 * exponent) - - -@file_or_str(".molden", ".input") -def shells_from_molden(text): - molden_shells = molden.shells_from_molden(text) - - dividers = { - 0: 1, - 1: 1, - 2: 3, - 3: 15, - 4: 35, - } - - def fix_contr_coeffs(l, coeffs, exponents): - """Fix contraction coefficients. Based on equations found in the SI - of the JANPA paper [1].""" - l = get_l(l) - divider = dividers[l] - rad_ints = radial_integral(l, exponents) - norms2 = rad_ints * 4 * np.pi / (2 * l + 1) / divider - norms = np.sqrt(norms2) - normed_coeffs = coeffs * norms - return normed_coeffs - - _shells = list() - for shell in molden_shells.shells: - L, center, _, exps, *_ = shell.as_tuple() - fixed_coeffs = fix_contr_coeffs(L, shell.coeffs_org, exps) - fixed_shell = Shell( - L=L, - center=center, - coeffs=fixed_coeffs, - exps=exps, - center_ind=shell.center_ind, - atomic_num=shell.atomic_num, - ) - _shells.append(fixed_shell) - shells = ORCAMoldenShells(_shells) - return shells - - @file_or_str(".molden", ".input") -def wavefunction_from_molden(text, charge=None, **wf_kwargs): - shells_func = shells_from_molden +def wavefunction_from_orca_molden(text, charge=None, **wf_kwargs): return molden.wavefunction_from_molden( - text, charge=charge, shells_func=shells_func, **wf_kwargs + text, charge=charge, orca_contr_renorm=True, **wf_kwargs ) diff --git a/pysisyphus/wavefunction/pop_analysis.py b/pysisyphus/wavefunction/pop_analysis.py index 70326c8fdb..3d25cfee35 100644 --- a/pysisyphus/wavefunction/pop_analysis.py +++ b/pysisyphus/wavefunction/pop_analysis.py @@ -11,8 +11,11 @@ @dataclass class PopAnalysis: + # Also store atoms and coordinates?! # atoms: List[str] # coords3d: NDArray[float] + + # Population of alpha and beta electrons pop_a: NDArray[float] pop_b: NDArray[float] nuc_charges: NDArray[float] diff --git a/pysisyphus/wavefunction/wavefunction.py b/pysisyphus/wavefunction/wavefunction.py index ffa60b0929..5e14ab78c4 100644 --- a/pysisyphus/wavefunction/wavefunction.py +++ b/pysisyphus/wavefunction/wavefunction.py @@ -135,6 +135,9 @@ def mo_num(self): def from_file(fn, **kwargs): path = Path(fn) + if not path.exists(): + raise FileNotFoundError(path) + from_funcs = { ".json": Wavefunction.from_orca_json, ".fchk": Wavefunction.from_fchk, @@ -142,10 +145,16 @@ def from_file(fn, **kwargs): from_funcs_for_str = ( # ORCA ("Molden file created by orca_2mkl", Wavefunction.from_orca_molden), + # XTB detection, as done in Multiwfn. Doesn't appear too stable/specific ... + ("[Atoms] AU", Wavefunction.from_orca_molden), + # AOMix, e.g. from Turbomole ("[AOMix Format", Wavefunction.from_aomix), - # OpenMolcas - # ("[N_Atoms]", Wavefunction.from_molden), # seems buggy right now + # OpenMolcas; seems buggy. Maybe also related to messed up contr. coeffs? + # ("[N_Atoms]", Wavefunction.from_molden), + # # General Molden fallback + # ("[Molden Format]", Wavefunction.from_molden), ) + # If possible I would advise to stay away from .molden files :) try: from_func = from_funcs[path.suffix.lower().strip()] except KeyError: @@ -194,9 +203,9 @@ def from_orca_molden(text, **kwargs): wavefunction_from_molden will come up with an absurdly high charge. """ - from pysisyphus.io.orca import wavefunction_from_molden + from pysisyphus.io.orca import wavefunction_from_orca_molden - wf = wavefunction_from_molden(text, **kwargs) + wf = wavefunction_from_orca_molden(text, **kwargs) return wf @staticmethod diff --git a/tests/test_orca/test_orca.py b/tests/test_orca/test_orca.py index a0711abe64..20b5c732f6 100644 --- a/tests/test_orca/test_orca.py +++ b/tests/test_orca/test_orca.py @@ -136,9 +136,9 @@ def test_parse_orca_cis(method, tda, triplets, this_dir): cis_fn = this_dir / "ref_cis" / f"{base_name}_000.000.orca.cis" Xa, Ya, Xb, Yb = parse_orca_cis(cis_fn) - ac = Xa ** 2 - Ya ** 2 + ac = Xa**2 - Ya**2 a_sum = np.sum(ac, axis=(1, 2)) - bc = Xb ** 2 - Yb ** 2 + bc = Xb**2 - Yb**2 b_sum = np.sum(bc, axis=(1, 2)) Y0 = np.zeros_like(Ya) @@ -170,7 +170,7 @@ def print_summary(X, Y, aorb): print(X) print(f"Y{aorb}") print(Y) - norms = X ** 2 - Y ** 2 + norms = X**2 - Y**2 norms_sum = np.sum(norms, axis=(1, 2)) print(f"@{method}, {tda=}, {aorb} coeffs, {norms_sum}") print(norms) @@ -199,7 +199,6 @@ def assert_dens_mats(dens_dict, json_fn): ), ) def test_orca_gs_densities(dens_fn, json_fn): - dens_dict = parse_orca_densities(WF_LIB_DIR / dens_fn) _ = assert_dens_mats(dens_dict, WF_LIB_DIR / json_fn) @@ -230,7 +229,17 @@ def test_orca_es_densities(dens_fn, json_fn, ref_dpm): @using("orca") def test_orca_hf(): geom = geom_loader("lib:h2o.xyz") - calc = ORCA("hf sto-3g", json_dump=True) + calc = ORCA("hf sto-3g") geom.set_calculator(calc) energy = geom.energy assert energy == pytest.approx(-74.960702484) + + +@using("orca") +def test_orca_stored_wavefunction(): + geom = geom_loader("lib:h2o.xyz") + calc = ORCA("hf sto-3g", wavefunction_dump=True) + geom.set_calculator(calc) + geom.energy + # Wavefunction already does some internal sanity checking + calc.get_stored_wavefunction() diff --git a/tests/test_xtb/test_xtb.py b/tests/test_xtb/test_xtb.py index 73f9bc4e24..486755b532 100644 --- a/tests/test_xtb/test_xtb.py +++ b/tests/test_xtb/test_xtb.py @@ -2,10 +2,13 @@ import pytest from pysisyphus.calculators import XTB +from pysisyphus.io.molden import get_xtb_nuc_charges from pysisyphus.dynamics.helpers import get_mb_velocities_for_geom +from pysisyphus.elem_data import INV_ATOMIC_NUMBERS from pysisyphus.helpers import geom_loader from pysisyphus.helpers_pure import eigval_to_wavenumber from pysisyphus.testing import using +from pysisyphus.wavefunction.pop_analysis import mulliken_charges_from_wf @pytest.fixture @@ -81,3 +84,26 @@ def test_xtb_retry_calc(this_dir): geom.set_calculator(XTB(pal=6, retry_etemp=1000.0, retry_calc=1)) en = geom.energy assert en == pytest.approx(-108.70462911348) + + +@using("xtb") +def test_xtb_stored_wavefunction(): + geom = geom_loader("lib:h2o.xyz") + calc = XTB(pal=1, wavefunction_dump=True) + geom.set_calculator(calc) + geom.energy + # Wavefunction already does some internal sanity checking + wf = calc.get_stored_wavefunction() + assert wf.charge == 0 + pa = mulliken_charges_from_wf(wf) + ref_charges = (-0.56619393, 0.28309696, 0.28309696) + np.testing.assert_allclose(pa.charges, ref_charges, atol=1e-6) + + +def test_xtb_nuc_charges(): + # Parametrized up to Z = 86 + atomic_nums = np.arange(1, 87, dtype=int) + all_atoms = [INV_ATOMIC_NUMBERS[Z] for Z in atomic_nums] + mod_charges = get_xtb_nuc_charges(all_atoms) + ecp_electrons = get_xtb_nuc_charges(all_atoms, as_ecp_electrons=True) + np.testing.assert_allclose(mod_charges + ecp_electrons, atomic_nums)