From 5da258363f39dd2ecf197eb5ff81bd856f5e1608 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Sat, 8 Aug 2026 16:12:32 +0200 Subject: [PATCH 01/20] xtb-gradient grab bugfix: failed for larger systems, now fixed --- ash/interfaces/interface_xtb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ash/interfaces/interface_xtb.py b/ash/interfaces/interface_xtb.py index 32773f9c3..09d254b4a 100644 --- a/ash/interfaces/interface_xtb.py +++ b/ash/interfaces/interface_xtb.py @@ -754,7 +754,7 @@ def xtbgradientgrab(numatoms): with open('gradient') as f: for line in reverse_lines(f): if ' cycle =' in line: - energy=float(line.split()[6]) + energy = float(line.split("SCF energy =")[1].split()[0]) return energy, gradient if count==numatoms: grab=False From df32a2a3d0beee81d77eec44ca244277a6061e84 Mon Sep 17 00:00:00 2001 From: wuls968 <201544050+wuls968@users.noreply.github.com> Date: Tue, 11 Aug 2026 13:43:45 +0800 Subject: [PATCH 02/20] fix: stop PySCF runs before using unconverged states --- ash/interfaces/interface_pyscf.py | 8 ++ ash/tests/test_pyscf_scf_convergence.py | 145 ++++++++++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 ash/tests/test_pyscf_scf_convergence.py diff --git a/ash/interfaces/interface_pyscf.py b/ash/interfaces/interface_pyscf.py index 10a4bb426..f065a4c6c 100644 --- a/ash/interfaces/interface_pyscf.py +++ b/ash/interfaces/interface_pyscf.py @@ -2569,6 +2569,14 @@ def run_SCF(self,mf=None, dm=None, max_cycle=None): #Grid printing scf_result = mf.run(dm) + # max_cycle<=0 is an intentional non-SCF evaluation used by scf_noiter. + if not scf_result.converged and mf.max_cycle > 0: + ashexit( + errormessage=( + f"PySCF SCF did not converge (max_cycle={mf.max_cycle}). " + "Refusing to continue with an unconverged electronic state." + ) + ) #Grid if 'KS' in self.scf_type: try: diff --git a/ash/tests/test_pyscf_scf_convergence.py b/ash/tests/test_pyscf_scf_convergence.py new file mode 100644 index 000000000..7987de124 --- /dev/null +++ b/ash/tests/test_pyscf_scf_convergence.py @@ -0,0 +1,145 @@ +import sys +from types import ModuleType, SimpleNamespace + +import numpy as np +import pytest + +from ash.interfaces.interface_pyscf import PySCFTheory + + +class FakeRHF: + def __init__(self, *, converged, max_cycle=50): + self.converged = converged + self.max_cycle = max_cycle + self.e_tot = -75.0 + self.mo_occ = np.array([2.0]) + self.make_rdm1_calls = 0 + self.nuc_grad_method_calls = 0 + self.grad_hcore_mm_calls = 0 + self.grad_nuc_mm_calls = 0 + + def run(self, dm): + return self + + def make_rdm1(self): + self.make_rdm1_calls += 1 + return np.array([[2.0]]) + + def nuc_grad_method(self): + self.nuc_grad_method_calls += 1 + return FakeGradient(self) + + +class FakeGradient: + def __init__(self, mf): + self.mf = mf + + def kernel(self): + return np.array([[0.0, 0.0, 0.0]]) + + def grad_hcore_mm(self, dm): + self.mf.grad_hcore_mm_calls += 1 + return np.array([[0.0, 0.0, 0.0]]) + + def grad_nuc_mm(self): + self.mf.grad_nuc_mm_calls += 1 + return np.array([[0.0, 0.0, 0.0]]) + + +def install_fake_pyscf(monkeypatch): + pyscf_module = ModuleType("pyscf") + pyscf_dft_module = ModuleType("pyscf.dft") + pyscf_dft_module.rks = SimpleNamespace(RKS=type("FakeRKS", (), {})) + pyscf_module.dft = pyscf_dft_module + pyscf_module.scf = SimpleNamespace(hf=SimpleNamespace(RHF=FakeRHF)) + monkeypatch.setitem(sys.modules, "pyscf", pyscf_module) + monkeypatch.setitem(sys.modules, "pyscf.dft", pyscf_dft_module) + + +def make_theory(mf): + theory = PySCFTheory.__new__(PySCFTheory) + theory.printlevel = 0 + theory.mf = mf + theory.scf_type = "RHF" + theory.functional = None + theory.platform = "CPU" + theory.periodic = False + return theory + + +def test_run_scf_aborts_before_density_for_unconverged_result(monkeypatch): + install_fake_pyscf(monkeypatch) + mf = FakeRHF(converged=False) + theory = make_theory(mf) + + with pytest.raises(SystemExit): + theory.run_SCF() + + assert mf.make_rdm1_calls == 0 + + +def test_run_scf_allows_explicit_zero_cycle_evaluation(monkeypatch): + install_fake_pyscf(monkeypatch) + mf = FakeRHF(converged=False) + theory = make_theory(mf) + + result = theory.run_SCF(max_cycle=0) + + assert result is mf + assert mf.max_cycle == 0 + assert mf.make_rdm1_calls == 1 + np.testing.assert_array_equal(theory.dm, np.array([[2.0]])) + + +def test_run_scf_allows_explicit_negative_cycle_evaluation(monkeypatch): + install_fake_pyscf(monkeypatch) + mf = FakeRHF(converged=False, max_cycle=-1) + theory = make_theory(mf) + + result = theory.run_SCF() + + assert result is mf + assert mf.make_rdm1_calls == 1 + + +def test_actualrun_aborts_before_gradient_and_pc_gradient_for_unconverged_scf(monkeypatch): + install_fake_pyscf(monkeypatch) + mf = FakeRHF(converged=False) + theory = make_theory(mf) + theory.SCF = True + theory.verbose_setting = 0 + theory.write_chkfile_name = None + theory.BS = False + theory.setup_guess = lambda: None + theory.run_stability_analysis = lambda: None + theory.do_pop_analysis = False + theory.get_dipole_moment = lambda: None + theory.dispersion = None + theory.NMF = False + theory.losc = False + theory.postSCF = False + theory.PC_gradient_code = "new" + + with pytest.raises(SystemExit): + theory.actualrun( + Grad=True, + PC=True, + current_MM_coords=np.array([[1.0, 0.0, 0.0]]), + MMcharges=np.array([1.0]), + ) + + assert mf.nuc_grad_method_calls == 0 + assert mf.grad_hcore_mm_calls == 0 + assert mf.grad_nuc_mm_calls == 0 + + +def test_run_scf_preserves_density_for_converged_result(monkeypatch): + install_fake_pyscf(monkeypatch) + mf = FakeRHF(converged=True) + theory = make_theory(mf) + + result = theory.run_SCF() + + assert result is mf + assert mf.make_rdm1_calls == 1 + np.testing.assert_array_equal(theory.dm, np.array([[2.0]])) From 1ead87acd94f3986d5a510bf9738a451d1ce7074 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 19 Aug 2026 17:49:23 +0200 Subject: [PATCH 03/20] cp2k interface update for cp2k 2026.2: fix for xtb-tblite calculations --- ash/interfaces/interface_CP2K.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index b2047de4c..9e7613583 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -743,6 +743,7 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, xtbcode = int(''.join(filter(str.isdigit, xtb_type))) inpfile.write(f' &XTB\n') if xtb_tblite is True: + inpfile.write(f' GFN_TYPE TBLITE\n') #NOTE inpfile.write(f' &TBLITE\n') inpfile.write(f' METHOD {xtb_type}\n') inpfile.write(f' &END\n') From 22e55b7db20b679f358c4998b63d3805f946adbd Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Thu, 20 Aug 2026 12:53:41 +0200 Subject: [PATCH 04/20] WrapTheory: basic support for PBC (untested) --- ash/modules/module_hybridtheory.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/ash/modules/module_hybridtheory.py b/ash/modules/module_hybridtheory.py index d41a461f8..80fd1d070 100644 --- a/ash/modules/module_hybridtheory.py +++ b/ash/modules/module_hybridtheory.py @@ -11,6 +11,7 @@ from ash.functions.functions_general import BC, ashexit, print_time_rel,print_line_with_mainheader,listdiff from ash.modules.module_theory import Theory +from ash.modules.module_coords_PBC import cell_params_to_vectors, cell_vectors_to_params from ash.interfaces.interface_ORCA import grabatomcharges_ORCA from ash.interfaces.interface_xtb import grabatomcharges_xTB,grabatomcharges_xTB_output import ash.constants @@ -373,7 +374,8 @@ class WrapTheory(Theory): def __init__(self, theory1=None, theory2=None, theories=None, printlevel=2, label=None, theory1_atoms=None, theory2_atoms=None, theory3_atoms=None, theory4_atoms=None, theory5_atoms=None, - theory_operators=None): + theory_operators=None, + periodic=False, periodic_cell_vectors=None, periodic_cell_dimensions=None): super().__init__() self.theorytype="QM" @@ -391,6 +393,11 @@ def __init__(self, theory1=None, theory2=None, theories=None, printlevel=2, labe self.theory4_atoms=theory4_atoms self.theory5_atoms=theory5_atoms + # PBC + self.periodic_cell_dimensions=periodic_cell_dimensions + self.periodic_cell_vectors=periodic_cell_vectors + self.periodic=periodic + print_line_with_mainheader(f"{self.theorynamelabel} initialization") print("Creating WrapTheory object") @@ -416,6 +423,27 @@ def __init__(self, theory1=None, theory2=None, theories=None, printlevel=2, labe print(f"Error: Number of theory-operators {len(self.theory_operators)} is not equal to number of theories {len(self.theories)}") ashexit() + # Return cell gradient if periodic is True + # Find first theory that has a cell gradient and return it + def get_cell_gradient(self): + self.cell_gradient = None + # Loop over theories and find the first theory that has a cell gradient + for theory in self.theories: + if hasattr(theory, 'cell_gradient') and theory.cell_gradient is not None: + self.cell_gradient = theory.cell_gradient + return self.cell_gradient + + # Update cell using either periodic_cell_vectors or periodic_cell_dimensions + def update_cell(self,periodic_cell_vectors=None, periodic_cell_dimensions=None): + print("Updating cell vectors") + if periodic_cell_vectors is not None: + self.periodic_cell_vectors = periodic_cell_vectors + + self.periodic_cell_dimensions = cell_vectors_to_params(periodic_cell_vectors) + elif periodic_cell_dimensions is not None: + self.periodic_cell_dimensions=periodic_cell_dimensions + + self.periodic_cell_vectors = cell_params_to_vectors(periodic_cell_dimensions) def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_elems=None, mm_elems=None, elems=None, Grad=False, PC=False, numcores=None, restart=False, label=None, From 6b8bfb5ab7c8cce912e3b7fe45d7d15bcc6fb505 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Thu, 20 Aug 2026 12:56:12 +0200 Subject: [PATCH 05/20] fix --- ash/modules/module_hybridtheory.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/ash/modules/module_hybridtheory.py b/ash/modules/module_hybridtheory.py index 80fd1d070..1b4299b8c 100644 --- a/ash/modules/module_hybridtheory.py +++ b/ash/modules/module_hybridtheory.py @@ -394,10 +394,16 @@ def __init__(self, theory1=None, theory2=None, theories=None, printlevel=2, labe self.theory5_atoms=theory5_atoms # PBC - self.periodic_cell_dimensions=periodic_cell_dimensions - self.periodic_cell_vectors=periodic_cell_vectors self.periodic=periodic - + if self.periodic: + if periodic_cell_dimensions is not None: + print("periodic_cell_dimensions:", periodic_cell_dimensions) + self.periodic_cell_dimensions = periodic_cell_dimensions + # Convert to cell vectors + self.periodic_cell_vectors = cell_params_to_vectors(periodic_cell_dimensions) + elif periodic_cell_vectors is not None: + self.periodic_cell_vectors = periodic_cell_vectors + self.periodic_cell_dimensions = cell_vectors_to_params(periodic_cell_vectors) print_line_with_mainheader(f"{self.theorynamelabel} initialization") print("Creating WrapTheory object") From 00255159a3abe387ae41b22b4ab2a7de35a1f3f0 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 15:53:52 +0200 Subject: [PATCH 06/20] CP2K: Skala Gauxc option --- ash/interfaces/interface_CP2K.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 9e7613583..ab860d2ea 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -53,7 +53,8 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel center_coords=False, scf_maxiter=50, outer_scf_maxiter=10, scf_convergence=1e-6, eps_default=1e-10, coupling='GAUSSIAN', GEEP_num_gauss=6, MM_radius_scaling=1, mm_radii=None, OT=True, OT_minimizer='DIIS', OT_preconditioner='FULL_ALL', - OT_linesearch='3PNT', outer_SCF=True, outer_SCF_optimizer='SD', OT_energy_gap=0.08): + OT_linesearch='3PNT', outer_SCF=True, outer_SCF_optimizer='SD', OT_energy_gap=0.08, + path_to_gauxc_model=None): self.theorytype="QM" self.theorynamelabel="CP2K" @@ -238,6 +239,10 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel self.outer_scf_maxiter=outer_scf_maxiter self.eps_default=eps_default + # Skala GauXC stuff + # NOte: functional needs to be Skala and then path_to_gauxc_model should be path to the .fun model file + self.path_to_gauxc_model=path_to_gauxc_model + # K-points self.kpoint_settings=kpoint_settings @@ -458,7 +463,8 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el scf_maxiter=self.scf_maxiter, outer_scf_maxiter=self.outer_scf_maxiter, ngrids=self.ngrids, xc_finer_grid=self.xc_finer_grid, cutoff=self.cutoff, rel_cutoff=self.rel_cutoff, printlevel=self.printlevel, OT=self.OT, OT_minimizer=self.OT_minimizer, OT_preconditioner=self.OT_preconditioner, OT_linesearch=self.OT_linesearch, - outer_SCF=self.outer_SCF, outer_SCF_optimizer=self.outer_SCF_optimizer, OT_energy_gap=self.OT_energy_gap) + outer_SCF=self.outer_SCF, outer_SCF_optimizer=self.outer_SCF_optimizer, OT_energy_gap=self.OT_energy_gap, + path_to_gauxc_model=self.path_to_gauxc_model) else: #No QM/MM #QM-CELL @@ -500,7 +506,8 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el basis_file=self.basis_file, potential_file=self.potential_file, psolver=self.psolver, OT=self.OT, OT_minimizer=self.OT_minimizer, OT_preconditioner=self.OT_preconditioner, OT_linesearch=self.OT_linesearch, - outer_SCF=self.outer_SCF, outer_SCF_optimizer=self.outer_SCF_optimizer, OT_energy_gap=self.OT_energy_gap) + outer_SCF=self.outer_SCF, outer_SCF_optimizer=self.outer_SCF_optimizer, OT_energy_gap=self.OT_energy_gap, + path_to_gauxc_model=self.path_to_gauxc_model) #Delete old forces file if present try: @@ -635,7 +642,8 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, qm_kind_dict=None, mm_kind_list=None, mm_ewald_type='NONE', mm_ewald_alpha=0.35, mm_ewald_gmax="21 21 21", printlevel=2, OT=False, OT_minimizer='DIIS', OT_preconditioner='FULL_ALL', OT_linesearch='3PNT', - outer_SCF=False, outer_SCF_optimizer='DIIS', OT_energy_gap=0.08): + outer_SCF=False, outer_SCF_optimizer='DIIS', OT_energy_gap=0.08, + path_to_gauxc_model=None): if method == 'QMMM': if mm_radii == None: print("No user MM radii provided. Will use default radii from internal dict (element_radii_for_cp2k).") @@ -797,6 +805,11 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, inpfile.write(f' &END PAIR_POTENTIAL\n') inpfile.write(f' &END VDW_POTENTIAL\n') inpfile.write(f' &XC_FUNCTIONAL {functional}\n') + # GauXC + if functional.upper() in ['SKALA']: + inpfile.write(f' &GauXC\n') + inpfile.write(f' MODEL {path_to_gauxc_model}\n') + inpfile.write(f' &END GauXC\n') inpfile.write(f' &END XC_FUNCTIONAL\n') inpfile.write(f' &END XC\n') From 31ca241b6b46ff0f1a43d4d9503872cf03f2be00 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 16:00:52 +0200 Subject: [PATCH 07/20] cp2k: fix for nonperiodic calcs --- ash/interfaces/interface_CP2K.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index ab860d2ea..0e4a3c173 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -116,6 +116,8 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel else: print("Periodic is False") self.periodic_type='NONE' + self.periodic_cell_dimensions=None + self.periodic_cell_vectors=None print("PERIODIC_TYPE:", self.periodic_type) # Parallelization From bef0a9833e75843693fc80d52369ea145cea64a9 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 16:04:30 +0200 Subject: [PATCH 08/20] cp2k: skala gauxc fix --- ash/interfaces/interface_CP2K.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 0e4a3c173..0394e609e 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -806,12 +806,16 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, inpfile.write(f' &END PRINT_DFTD\n') inpfile.write(f' &END PAIR_POTENTIAL\n') inpfile.write(f' &END VDW_POTENTIAL\n') - inpfile.write(f' &XC_FUNCTIONAL {functional}\n') + # GauXC if functional.upper() in ['SKALA']: + inpfile.write(f' &XC_FUNCTIONAL \n') inpfile.write(f' &GauXC\n') inpfile.write(f' MODEL {path_to_gauxc_model}\n') inpfile.write(f' &END GauXC\n') + else: + inpfile.write(f' &XC_FUNCTIONAL {functional}\n') + inpfile.write(f' &END XC_FUNCTIONAL\n') inpfile.write(f' &END XC\n') From e0c4df3f0f0a87ee45baa7704a0718ad9edfb965 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 16:06:39 +0200 Subject: [PATCH 09/20] fx --- ash/interfaces/interface_CP2K.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 0394e609e..a227018af 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -810,9 +810,9 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, # GauXC if functional.upper() in ['SKALA']: inpfile.write(f' &XC_FUNCTIONAL \n') - inpfile.write(f' &GauXC\n') + inpfile.write(f' &GAUXC\n') inpfile.write(f' MODEL {path_to_gauxc_model}\n') - inpfile.write(f' &END GauXC\n') + inpfile.write(f' &END GAUXC\n') else: inpfile.write(f' &XC_FUNCTIONAL {functional}\n') From 49872a7826eb64688e08bb5841548628855dbd30 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 16:27:26 +0200 Subject: [PATCH 10/20] cp2k: fix for nonperiodic calculations (was broken) --- ash/interfaces/interface_CP2K.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index a227018af..68092e32c 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -485,6 +485,8 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el qm_box_dims=np.around(qm_box_dims + self.qm_cell_shift_par,1) print(f"Setting Cell box dimensions to {qm_box_dims} Angstrom") self.cell_dimensions=list(qm_box_dims)+[90.0,90.0,90.0] + self.cell_vectors=cell_params_to_vectors(self.cell_dimensions) + self.periodic_cell_vectors=self.cell_vectors #Write xyz-file with coordinates system_xyzfile="system_cp2k" From 9da43235fda27b2b9a87e09fdff1be203d2fc41e Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Mon, 24 Aug 2026 16:30:15 +0200 Subject: [PATCH 11/20] fix --- ash/interfaces/interface_CP2K.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 68092e32c..de929d772 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -116,8 +116,9 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel else: print("Periodic is False") self.periodic_type='NONE' - self.periodic_cell_dimensions=None - self.periodic_cell_vectors=None + # Note: still using these keywords for simpler logic + self.periodic_cell_dimensions=cell_dimensions + self.periodic_cell_vectors=cell_vectors print("PERIODIC_TYPE:", self.periodic_type) # Parallelization From fee0452666c6f73da1c763e2b4d6bd25f22684a9 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 11:18:11 +0200 Subject: [PATCH 12/20] localized_orbital_analysis function --- ash/__init__.py | 2 +- ash/functions/functions_elstructure.py | 416 ++++++++++++++++++++++++- ash/interfaces/interface_ORCA.py | 6 +- 3 files changed, 421 insertions(+), 3 deletions(-) diff --git a/ash/__init__.py b/ash/__init__.py index a3ec98dd1..ac5abb93a 100644 --- a/ash/__init__.py +++ b/ash/__init__.py @@ -72,7 +72,7 @@ from .functions.functions_elstructure import read_cube, write_cube, write_cube_diff, diffdens_tool, create_cubefile_from_orbfile, diffdens_of_cubefiles, \ NOCV_density_ORCA, difference_density_ORCA, NOCV_Multiwfn,write_cube_sum,write_cube_product,create_density_from_orb, make_molden_file, \ diagonalize_DM_AO, diagonalize_DM, DM_AO_to_MO, DM_AO_to_MO, DM_MO_to_AO, select_space_from_occupations,select_indices_from_occupations, ASH_write_integralfile, \ - density_sensitivity_metric + density_sensitivity_metric, localized_orbital_analysis #multiwfn interface import ash.interfaces.interface_multiwfn diff --git a/ash/functions/functions_elstructure.py b/ash/functions/functions_elstructure.py index 8391faf10..b87ea1b63 100644 --- a/ash/functions/functions_elstructure.py +++ b/ash/functions/functions_elstructure.py @@ -1,3 +1,5 @@ +import sys + import numpy as np import math import shutil @@ -16,6 +18,7 @@ from ash.dictionaries_lists import eldict from ash.constants import hartokcal from ash.interfaces.interface_multiwfn import multiwfn_run +from ash.modules.module_singlepoint import Singlepoint #CM5. from https://github.com/patrickmelix/CM5-calculator/blob/master/cm5calculator.py @@ -2610,4 +2613,415 @@ def boltzmann_populations(energies, temperature=298.15): boltzmann_factors = np.exp(-1*rel_energies * beta) populations = boltzmann_factors / np.sum(boltzmann_factors) print("Boltzmann populations at", temperature, "K:", populations) - return populations \ No newline at end of file + return populations + +################################## +# localized_orbital_analysis() +################################### +# Processing of Loewdin reduced MO analysis from ORCA output +# Functions get_loewdin_table, convert_loewdin_text_to_dataframe, select_contr and plot_dataframe_table +# were written by Nico Spiller (taken from https://gitlab.gwdg.de/orca-helpers/orca-helpers/-/tree/master/orca_helpers/output_files/loewdin_plotter) +def get_loewdin_table(file): + '''get LOEWDIN REDUCED ORBITAL POPULATIONS PER MO data from output file + + returns list of blocks + if unrestricted, list[0] contains spin up, list[1] contains spin down + if restricted, list[0] contains restricted orbitals + + + example unrestricted: + + ------------------------------------------ + LOEWDIN REDUCED ORBITAL POPULATIONS PER MO + ------------------------------------------- + THRESHOLD FOR PRINTING IS 0.1% + SPIN UP + 0 1 2 3 4 5 + -144.18793 -18.63385 -14.84937 -12.05575 -12.05213 -12.05213 + 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 + -------- -------- -------- -------- -------- -------- + 0 Ca s 100.0 0.0 99.7 0.0 0.0 0.0 + 0 Ca pz 0.0 0.0 0.0 100.0 0.0 0.0 + + + example restricted: + + ------------------------------------------ + LOEWDIN REDUCED ORBITAL POPULATIONS PER MO + ------------------------------------------- + THRESHOLD FOR PRINTING IS 0.1% + 0 1 2 3 4 5 + -144.18789 -18.63406 -14.84933 -12.05572 -12.05209 -12.05209 + 2.00000 2.00000 2.00000 2.00000 2.00000 2.00000 + -------- -------- -------- -------- -------- -------- + 0 Ca s 100.0 0.0 99.7 0.0 0.0 0.0 + 0 Ca pz 0.0 0.0 0.0 100.0 0.0 0.0 + 0 Ca px 0.0 0.0 0.0 0.0 0.0 100.0 + ''' + + empty_lines = 0 # counter for empty lines + match = 'LOEWDIN REDUCED ORBITAL POPULATIONS PER MO' # start parsing here + match_casscf = 'LOEWDIN ORBITAL-COMPOSITIONS' # CASSCF prints different string + break_casscf = '----------------------------' # break early for CASSCF + parse = False # switch parsing + spin = 2 # 0: spin up, 1: spin down, 2: restricted + data = [ [], [], [] ] + + with open(file) as f: + for line in f: + + # turn on parsing + if not parse and ( match in line or match_casscf in line ): + parse = True + next(f) # skip next two lines + next(f) + continue + + # stop reading + elif parse and ( not len(line.split()) or break_casscf in line ): + empty_lines += 1 + # stop if two consecutive empty lines found in spin = 1 or 2 blocks + if empty_lines == 2: + if spin: + break + else: + continue + + else: + empty_lines = 0 + + # read data + if parse: + if 'SPIN UP' in line: + spin = 0 + continue + elif 'SPIN DOWN' in line: + spin = 1 + continue + + # writes to: data[0] for alpha, data[1] for beta, data[2] for restricted + data[spin].append(line) + + # remove empty blocks + data = [ i for i in data if i ] + + return data + +def convert_loewdin_text_to_dataframe(block, nmax=None, only_occ=False): + '''convert plain text data into pandas dataframes + returns two dataframes, the first containing the contributions, the second orbital information + + takes a list containing the plain text data from ORCA output file + + optional arguments + nmax: integer or None, index of highest MO to consider + only_occ: True or False, stop at first occupation 0.0 + ''' + import pandas as pd + # get reduced ao basis for df: + # dict of dict of set: {a_no1: {l1: {ml11, ml12 ,...}, l2: {ml21, ml22}, }, ... } + basis = dict() + # dict with atoms to index list: {'Fe': [0, 1], 'S': [2, 3]} + atom2index = dict() + + for line in block: + + l = line[:12].split() # look at fixed width column only + + if len(l): # only lines with data in first couple of characters + a_no = int(l[0]) + a = l[1] + ao_l = l[-1][0] + ao_ml = l[-1][1:] + + # fill nested dict, or create subdict + if a_no in basis: + + if ao_l in basis[a_no]: + basis[a_no][ao_l].add(ao_ml) + else: + basis[a_no][ao_l] = { ao_ml } + + else: + basis[a_no] = {ao_l: { ao_ml } } + + # fill atom2index + if a in atom2index: + atom2index[a].add(a_no) + else: + atom2index[a] = { a_no } + + # ... dict with basis is complete now + + # convert to list of tuples for df index + columns1 = [] + for k1, v1 in basis.items(): + for k2, v2 in v1.items(): + for v in v2: + columns1.append((k1, k2, v)) + + # determine # of mos + empty_lines = [ i for i, line in enumerate(block) if not len(line.split())] + last_empty = empty_lines[-2] + last_mo = int( block[last_empty + 1].split()[-1] ) + + # basis for pandas dataframe + columns1 = pd.MultiIndex.from_tuples(columns1) + index = range( last_mo + 1) + + # empty numpy array for contributions + arr1 = np.zeros((len(index), len(columns1)) ) + # flatten multiindex: simple dict with mapping of each index to integer + columns2col_idx = { '{}{}{}'.format(*c): i for i, c in enumerate(columns1) } + + # empty numpy array for orbital information + arr2 = np.empty((len(index), 3)) # no, erg, occ + arr2[:] = np.nan + + # start parsing + parse_rows = False + header = True + parse_finished = False + + block_iterator = iter(block) + for line in block_iterator: + l = line.split() + + # check for new header (necessary to capture other than first) + if len(l) == 0: + header = True + continue + + # parse header lines: every 4 lines after empy line + if header: + # stop, if premature end + if parse_finished: + last_mo = row_idx[-1] + break + + row_idx = [int(i) for i in l] + + l = next(block_iterator).split() + + es = [ i for i in l ] + l = next(block_iterator).split() + + occs = [ float(i) for i in l ] + next(block_iterator) + + # check if reached nmax + if nmax and nmax in row_idx: + final_idx = row_idx.index(nmax) + parse_finished = True + + # check if still occ orbitals + elif only_occ and 0.0 in occs: + final_idx = occs.index(0.0) + parse_finished = True + + if parse_finished: + row_idx = row_idx[:final_idx+1] + es = es[:final_idx+1] + occs = occs[:final_idx+1] + + arr2[row_idx, 0] = row_idx + arr2[row_idx, 1] = es + arr2[row_idx, 2] = occs + + header = False + continue + + # parse data lines + contr = l[3:len(row_idx)+3] + a_no = l[0] + ele = l[1] + ao = l[2] + ao_l = ao[0] + ao_ml = ao[1:] + + # fill numpy array + col_idx = columns2col_idx['{}{}{}'.format(a_no, ao_l, ao_ml)] + arr1[row_idx, col_idx] = contr + + + # df1 is a multiindex dataframe + df1 = pd.DataFrame(data=arr1, index=index, columns=columns1) + df1 = df1.sort_index(axis=1, level=0, sort_remaining=False) # sort columns by atom numbers + df1 = df1.loc[:last_mo,] + + # df2 is a regular 2D dataframe + columns2 = [ 'mo', 'erg', 'occ' ] + df2 = pd.DataFrame(data=arr2, index=index, columns=columns2) + df2.loc[:, 'mo'] = pd.to_numeric(df2.loc[:, 'mo'], downcast='integer', errors='coerce') + df2 = df2.loc[:last_mo,] + + return df1, df2, atom2index + +def select_contr(df, atom2index, mos=[], atoms=[], subshells=[], collapse=0): + '''Select only part of the dataframe and add atom names + mos list of orbital indice (default: all) + atoms list of atom indices to + return smaller df''' + import pandas as pd + if not atoms: + atoms = slice(None) + print('INFO selecting atoms: all') + else: + try: + atoms = [ atom2index[i] if i in atom2index else int(i) for i in atoms ] # substitute atom names by indices + except ValueError: + print('ERROR: One of {} a valid atom'.format(' '.join(atoms))) + sys.exit(1) + + atoms = list(set(list(pd.core.common.flatten(atoms)))) # flatten possible sets and make sure only unique + print('INFO selecting atoms: {}'.format(' '.join([ str(i) for i in atoms]))) + + if not subshells: + subshells = slice(None) + print('INFO selecting l: all') + else: + print('INFO selecting l: {}'.format(' '.join(subshells))) + + if mos: + mos = range(mos[0], mos[1] + 1) + else: + mos = slice(None) + + + df = df.loc[mos, ( atoms, subshells, ) ] + + if collapse == 1: + #df = df.sum(level=(0, 1,), axis=1) + df = df.T.groupby(level=(0, 1)).sum().T + print('INFO collapsing data: ml') + elif collapse == 2: + #df = df.sum(level=(0, ), axis=1) + df = df.T.groupby(level=(0, )).sum().T + print('INFO collapsing data: l and ml') + + # rename df columns to include atom names + renamedict = dict() + for k, v in atom2index.items(): + for i in v: + renamedict[i] = '{}{}'.format(k, i) + df = df.rename(columns=renamedict) + + return df + +def plot_dataframe_table(df, path=None, annot=True, s=1.0): + '''Plot a single DataFrame''' + import matplotlib.pylab as plt + import seaborn as sns + y, x = [i/2/s for i in df.shape] + + fig, ax = plt.subplots(figsize=(x, y)) + + sns.heatmap( + data=df, + ax=ax, + annot=annot, + square=True, + linewidth=0.01, + vmin=0, + vmax=100, + cbar=not annot, + fmt='g', + cmap='magma' + ) + + ax.set_ylabel('MO #') + ax.set_xlabel('AO contr (Atom #-l-ml)') + + fig.tight_layout() + + if path: + fig.savefig(path, transparent=False) + + return fig, ax + +# Function to automate localized orbital analysis from ORCA output +def localized_orbital_analysis(theory=None, fragment=None, + angmom='d', threshold=70, atomlist=None, + plot_table=True, create_cubefiles=True, gridval=60, + debug=False): + + if theory is None or fragment is None: + print("Error: please provide a valid ASH Theory and Fragment") + ashexit() + if atomlist is None: + print("Error: please provide a list of atoms to analyze via atomlist keyword") + ashexit() + + # Run theory to get orbitals + print("Running theory to get orbitals") + + Singlepoint(theory=theory, fragment=fragment) + + + # ORCA noiter job to get Löwdin decomposition of localized orbitals + dummytheory=ORCATheory(orcasimpleinput=theory.orcasimpleinput+' noiter normalprint', + moreadfile=theory.filename+'.loc', filename="orca_noiter") + Singlepoint(theory=dummytheory, fragment=fragment) + + #Loewdin analysis + blocks = get_loewdin_table(dummytheory.filename+'.out') + nmax=None #highest mo index to consider + only_occ=True + dfarr = [] + for b in blocks: + df, _, atom2index = convert_loewdin_text_to_dataframe(b, nmax=nmax, only_occ=only_occ) + dfarr.append(df) + + # reduce dataframes + print("Selecting localized orbitals...") + for i, df in enumerate(dfarr): + print('INFO processing block: {}'.format(i)) + dfarr[i] = select_contr(df, atom2index, mos=[], atoms=atomlist, subshells=[], collapse=1) + dfarr_a= dfarr[0] + dfarr_b= dfarr[1] + if debug: + print("dfarr_a:", dfarr_a) + print("dfarr_b:", dfarr_b) + + # Getting column + col_a = next(c for c in dfarr_a.columns if c[1] == angmom) + col_b = next(c for c in dfarr_b.columns if c[1] == angmom) + + # Select d-orbitals for alpha and beta + indices_a = dfarr_a.index[dfarr_a[col_a] > threshold].tolist() + indices_b = dfarr_b.index[dfarr_b[col_b] > threshold].tolist() + + print() + print() + # Print relevant MOs + print("Found these localized MOs:") + print("alpha:") + for i in indices_a: + perc = dfarr_a.loc[i, col_a] + print(f"MO{i} {perc:7.1f}% {angmom}-character") + print("beta:") + for i in indices_b: + perc = dfarr_b.loc[i, col_b] + print(f"MO{i} {perc:7.1f}% {angmom}-character") + + if plot_table: + print("\n\nNow plotting tables for alpha and beta orbitals") + output="table_a.png" + print('INFO plotting: {}'.format(output)) + plot_dataframe_table(dfarr_a, path=output, annot=True, s=1.0) + output="table_b.png" + print('INFO plotting: {}'.format(output)) + plot_dataframe_table(dfarr_b, path=output, annot=True, s=1.0) + + #Create cubefiles for selected MO-numbers + if create_cubefiles: + print("\n\n Now creating cubefiles for selected MO-numbers") + for MONUMBER in indices_a: + print(f"Now plotting alpha MO {MONUMBER}") + run_orca_plot("orca.loc", "mo", gridvalue=gridval, mo_operator=0, mo_number=MONUMBER) + os.rename(f"orca.mo{MONUMBER}a.cube", f"orcaloccalc_alpha_mo{MONUMBER}.cube") #Renaming cube file + print() + for MONUMBER in indices_b: + print(f"Now plotting beta MO {MONUMBER}") + run_orca_plot("orca.loc", "mo", gridvalue=gridval, mo_operator=1, mo_number=MONUMBER) + os.rename(f"orca.mo{MONUMBER}b.cube", f"orcaloccalc_beta_mo{MONUMBER}.cube") #Renaming cube file \ No newline at end of file diff --git a/ash/interfaces/interface_ORCA.py b/ash/interfaces/interface_ORCA.py index 9f9d9d387..312b2b58d 100644 --- a/ash/interfaces/interface_ORCA.py +++ b/ash/interfaces/interface_ORCA.py @@ -2941,11 +2941,15 @@ def orblocfind(outputfile, atomindex_strings=None, popthreshold=0.1): stronggrab=True if stronggrab == True: for atindex in atomindex_strings: + #print("atindex:", atindex) if str(atindex) in line: - #print(line) + #print("if. line:", line) if operator == 'alpha': + #print("alpha") atom=line.split()[2] + #print("atom is", atom) monumber=line.split()[1][:-1] + #print("monumber is", monumber) dict_alpha.setdefault(atom, []).append(int(monumber)) elif operator == 'beta': atom=line.split()[2] From 33dfdba05dc2c2bb83aecbc5987bcd59503545f4 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 11:53:23 +0200 Subject: [PATCH 13/20] added tm ox assignment to localized_orbital analysis --- ash/functions/functions_elstructure.py | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/ash/functions/functions_elstructure.py b/ash/functions/functions_elstructure.py index b87ea1b63..95289fec0 100644 --- a/ash/functions/functions_elstructure.py +++ b/ash/functions/functions_elstructure.py @@ -2816,7 +2816,7 @@ def convert_loewdin_text_to_dataframe(block, nmax=None, only_occ=False): # check if still occ orbitals elif only_occ and 0.0 in occs: - final_idx = occs.index(0.0) + final_idx = occs.index(0.0)-1 parse_finished = True if parse_finished: @@ -2966,10 +2966,9 @@ def localized_orbital_analysis(theory=None, fragment=None, #Loewdin analysis blocks = get_loewdin_table(dummytheory.filename+'.out') nmax=None #highest mo index to consider - only_occ=True dfarr = [] for b in blocks: - df, _, atom2index = convert_loewdin_text_to_dataframe(b, nmax=nmax, only_occ=only_occ) + df, _, atom2index = convert_loewdin_text_to_dataframe(b, nmax=nmax, only_occ=True) dfarr.append(df) # reduce dataframes @@ -3004,6 +3003,26 @@ def localized_orbital_analysis(theory=None, fragment=None, perc = dfarr_b.loc[i, col_b] print(f"MO{i} {perc:7.1f}% {angmom}-character") + tmlist=['Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg'] + elements=['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr'] + if len(atomlist) == 1 and atomlist[0] in tmlist: + print("\n\nFound a single transition metal atom in atomlist, trying to determine oxidation state") + metal = atomlist[0] + + oxnumbers={-1:'-I',0:'0',-2:'-II',-3:'-III',-4:'-IV',-5:'-V',-6:'-VI',-7:'-VII',-8:'-VIII',-9:'-IX',-10:'-X',1:'I',2:'II',3:'III',4:'IV',5:'V',6:'VI',7:'VII',8:'VIII',9:'IX',10:'X'} + #D-valence electrons for some elements (counting atomic s-electrons as d) + delectrondict={26:['Fe',8],42:['Mo',6],23:['V',5]} + + alphanum=len(indices_a) + betanum=len(indices_b) + delectrons=alphanum+betanum + print("Whole d-electrons:",delectrons) + print(metal, str(alphanum)+'α', str(betanum)+'β') + oxstate=delectrondict[elements.index(metal)+1][1] - delectrons + print("Oxidation state: ", metal, oxnumbers[oxstate]) + + + if plot_table: print("\n\nNow plotting tables for alpha and beta orbitals") output="table_a.png" From 716f1c84833412c5ffe76252a5d52eca3ab793e3 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 12:00:08 +0200 Subject: [PATCH 14/20] ox state changes --- ash/functions/functions_elstructure.py | 51 +++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/ash/functions/functions_elstructure.py b/ash/functions/functions_elstructure.py index 95289fec0..dac686c4b 100644 --- a/ash/functions/functions_elstructure.py +++ b/ash/functions/functions_elstructure.py @@ -3010,8 +3010,55 @@ def localized_orbital_analysis(theory=None, fragment=None, metal = atomlist[0] oxnumbers={-1:'-I',0:'0',-2:'-II',-3:'-III',-4:'-IV',-5:'-V',-6:'-VI',-7:'-VII',-8:'-VIII',-9:'-IX',-10:'-X',1:'I',2:'II',3:'III',4:'IV',5:'V',6:'VI',7:'VII',8:'VIII',9:'IX',10:'X'} - #D-valence electrons for some elements (counting atomic s-electrons as d) - delectrondict={26:['Fe',8],42:['Mo',6],23:['V',5]} + # D-valence electrons for transition metal ions (NOTE: counting atomic s-electrons as d) + delectrondict = { + # 3d series + 21: ['Sc', 3], + 22: ['Ti', 4], + 23: ['V', 5], + 24: ['Cr', 6], + 25: ['Mn', 7], + 26: ['Fe', 8], + 27: ['Co', 9], + 28: ['Ni', 10], + 29: ['Cu', 11], + 30: ['Zn', 12], + + # 4d series + 39: ['Y', 3], + 40: ['Zr', 4], + 41: ['Nb', 5], + 42: ['Mo', 6], + 43: ['Tc', 7], + 44: ['Ru', 8], + 45: ['Rh', 9], + 46: ['Pd', 10], + 47: ['Ag', 11], + 48: ['Cd', 12], + + # 5d series + 57: ['La', 3], + 72: ['Hf', 4], + 73: ['Ta', 5], + 74: ['W', 6], + 75: ['Re', 7], + 76: ['Os', 8], + 77: ['Ir', 9], + 78: ['Pt', 10], + 79: ['Au', 11], + 80: ['Hg', 12], + + # 6d series + 89: ['Ac', 3], + 104: ['Rf', 4], + 105: ['Db', 5], + 106: ['Sg', 6], + 107: ['Bh', 7], + 108: ['Hs', 8], + 109: ['Mt', 9], + 110: ['Ds', 10], + 111: ['Rg', 11], + 112: ['Cn', 12]} alphanum=len(indices_a) betanum=len(indices_b) From 869c91f2146023d8a3b5c8c98cbe95a4708e80bb Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 12:02:39 +0200 Subject: [PATCH 15/20] linter fix --- ash/functions/functions_elstructure.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ash/functions/functions_elstructure.py b/ash/functions/functions_elstructure.py index dac686c4b..889debdd6 100644 --- a/ash/functions/functions_elstructure.py +++ b/ash/functions/functions_elstructure.py @@ -2784,6 +2784,7 @@ def convert_loewdin_text_to_dataframe(block, nmax=None, only_occ=False): parse_finished = False block_iterator = iter(block) + row_idx=None for line in block_iterator: l = line.split() From c04d349c4c9c64f83d24b7ef071386415be53f59 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 12:12:00 +0200 Subject: [PATCH 16/20] cp2k: fix for skala gauxc addition --- ash/interfaces/interface_CP2K.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index de929d772..83d99969c 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -130,7 +130,6 @@ def __init__(self, cp2kdir=None, cp2k_bin_name=None, filename='cp2k', printlevel #Finding CP2K dir and binaries self.cp2kdir, self.cp2k_bin_name = find_cp2k(cp2kdir,cp2k_bin_name) - #Checking OpenMPI if self.numcores != 1: print(f"Parallel job requested with numcores: {self.numcores}") @@ -810,17 +809,18 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, inpfile.write(f' &END PAIR_POTENTIAL\n') inpfile.write(f' &END VDW_POTENTIAL\n') - # GauXC - if functional.upper() in ['SKALA']: - inpfile.write(f' &XC_FUNCTIONAL \n') - inpfile.write(f' &GAUXC\n') - inpfile.write(f' MODEL {path_to_gauxc_model}\n') - inpfile.write(f' &END GAUXC\n') - else: - inpfile.write(f' &XC_FUNCTIONAL {functional}\n') + # XC and GauXC + if basis_method.upper() != "XTB": + if functional.upper() in ['SKALA']: + inpfile.write(f' &XC_FUNCTIONAL \n') + inpfile.write(f' &GAUXC\n') + inpfile.write(f' MODEL {path_to_gauxc_model}\n') + inpfile.write(f' &END GAUXC\n') + else: + inpfile.write(f' &XC_FUNCTIONAL {functional}\n') - inpfile.write(f' &END XC_FUNCTIONAL\n') - inpfile.write(f' &END XC\n') + inpfile.write(f' &END XC_FUNCTIONAL\n') + inpfile.write(f' &END XC\n') inpfile.write(f' &END DFT\n\n') From b0c09ed168e01e812cbef8a1bf4bcc020bbdff3f Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 12:13:33 +0200 Subject: [PATCH 17/20] fix --- ash/interfaces/interface_CP2K.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/ash/interfaces/interface_CP2K.py b/ash/interfaces/interface_CP2K.py index 83d99969c..682ccc5f1 100644 --- a/ash/interfaces/interface_CP2K.py +++ b/ash/interfaces/interface_CP2K.py @@ -820,8 +820,7 @@ def write_CP2K_input(method='QUICKSTEP', jobname='ash-CP2K', center_coords=True, inpfile.write(f' &XC_FUNCTIONAL {functional}\n') inpfile.write(f' &END XC_FUNCTIONAL\n') - inpfile.write(f' &END XC\n') - + inpfile.write(f' &END XC\n') inpfile.write(f' &END DFT\n\n') #QM/MM From f09a0cdd57c6f96f91b3aec039a310405c4f5ca8 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 14:35:27 +0200 Subject: [PATCH 18/20] gcp: debug printing --- ash/interfaces/interface_Grimme_corrections.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ash/interfaces/interface_Grimme_corrections.py b/ash/interfaces/interface_Grimme_corrections.py index 43dd36f18..c4f5a1d0d 100644 --- a/ash/interfaces/interface_Grimme_corrections.py +++ b/ash/interfaces/interface_Grimme_corrections.py @@ -126,12 +126,15 @@ def calc_gcp(fragment=None, xyzfile=None, current_coords=None, elems=None, funct xyzfile="gcpgeo.xyz" numatoms=len(current_coords) + + command_list=['mctc-gcp', xyzfile, '-l', functional] if Grad: command_list.append('--grad') # mctc-gcp appends into an existing 'gradient' file (Turbomole format) instead of # writing 'gcp_gradient'; remove it so the standalone file is always created. if os.path.exists('gradient'): + print("Removing old gradient file") os.remove('gradient') print("command_list:", command_list) From 3787806441ee41f93e58f18c254d54287d121452 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 16:19:53 +0200 Subject: [PATCH 19/20] gcp-interface: accounting for different behaviour of gcp-correction (2.3.2 vs 2.4.0), filename of gradient behaviour --- ash/interfaces/interface_Grimme_corrections.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/ash/interfaces/interface_Grimme_corrections.py b/ash/interfaces/interface_Grimme_corrections.py index c4f5a1d0d..d11112f2f 100644 --- a/ash/interfaces/interface_Grimme_corrections.py +++ b/ash/interfaces/interface_Grimme_corrections.py @@ -145,8 +145,21 @@ def calc_gcp(fragment=None, xyzfile=None, current_coords=None, elems=None, funct energy = float(result[1]) print("gcp energy:", energy) + # get gcp version number via mctc-gcp --version + version_result = pygrep("mctc-gcp", "gcp.out") + print("version_result:", version_result) + gcp_version = version_result[-1] if version_result else "unknown" + print("gcp_version:", gcp_version) + if Grad: - gradient = gcpgradientgrab(numatoms,"gcp_gradient") + # if gcp version >= 2.4.0 the file is called gradient + if gcp_version >= "2.4.0": + # gcp gradient is written to a file called 'gradient' + from ash.interfaces.interface_Turbomole import grab_gradient + gradient = grab_gradient(numatoms,file="gradient") + elif gcp_version < "2.4.0": + # if gcp version < 2.4.0 the gcp gradient is written to a file called 'gcp_gradient' + gradient = gcpgradientgrab(numatoms,"gcp_gradient") if printlevel > 2: print("gcp gradient:", gradient) return energy, gradient From 072801b91f180e69179a668ae08bae51271c7b99 Mon Sep 17 00:00:00 2001 From: RagnarB83 Date: Wed, 26 Aug 2026 16:24:36 +0200 Subject: [PATCH 20/20] grab_gradient turbomole routine now less sensitive to encoding issues --- ash/interfaces/interface_Turbomole.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ash/interfaces/interface_Turbomole.py b/ash/interfaces/interface_Turbomole.py index ec5995665..f2b0973ec 100644 --- a/ash/interfaces/interface_Turbomole.py +++ b/ash/interfaces/interface_Turbomole.py @@ -634,7 +634,7 @@ def grab_gradient_old(numatoms,file="gradient"): def grab_gradient(numatoms,file="gradient"): gradient = np.zeros((numatoms,3)) - with open(file, 'r') as gradfile: + with open(file, 'r', encoding="utf-8", errors="replace") as gradfile: gradlines = gradfile.readlines() #Reverse lines gradlines.reverse()