From c11fe7f99510bb569f4ba654a25995767d2993ff Mon Sep 17 00:00:00 2001 From: allenjbTX Date: Fri, 8 May 2026 14:04:57 -0400 Subject: [PATCH] - Ensure 'gcp_gradient' file is made by removing existing 'gradient' file created by xTBTheory - Update WrapTheory to optionally return pc_gradient when Grad and PC flags are set --- ash/interfaces/interface_Grimme_corrections.py | 8 +++++++- ash/modules/module_hybridtheory.py | 5 +++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/ash/interfaces/interface_Grimme_corrections.py b/ash/interfaces/interface_Grimme_corrections.py index 9cf5cda4c..43dd36f18 100644 --- a/ash/interfaces/interface_Grimme_corrections.py +++ b/ash/interfaces/interface_Grimme_corrections.py @@ -129,6 +129,10 @@ def calc_gcp(fragment=None, xyzfile=None, current_coords=None, elems=None, funct 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'): + os.remove('gradient') print("command_list:", command_list) with open('gcp.out', 'w') as ofile: @@ -171,7 +175,9 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el elems=None, Grad=False, PC=False, numcores=None, restart=False, label=None, charge=None, mult=None): - frag = ash.Fragment(elems=elems, coords=current_coords,printlevel=0) + if qm_elems is None: + qm_elems = elems + frag = ash.Fragment(elems=qm_elems, coords=current_coords,printlevel=0) eg = calc_gcp(fragment=frag, functional=self.functional, Grad=Grad, printlevel=self.printlevel) if Grad: self.energy = eg[0] diff --git a/ash/modules/module_hybridtheory.py b/ash/modules/module_hybridtheory.py index a3ec0fae5..d41a461f8 100644 --- a/ash/modules/module_hybridtheory.py +++ b/ash/modules/module_hybridtheory.py @@ -433,6 +433,7 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el full_dimension=current_coords.shape[0] energies=[] gradients=[] + pc_gradient=None chosen_coords=current_coords chosen_elems=qm_elems for i,theory in enumerate(self.theories): @@ -468,6 +469,8 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el #print(f"Theory: {theory.theorynamelabel} gradient shape", eg_tuple[1].shape) energy = eg_tuple[0] tempgrad = eg_tuple[1] + if PC and len(eg_tuple) == 3: + pc_gradient = eg_tuple[2] # Assemble gradient of correct dimension if i+1 == 1 and self.theory1_atoms is not None: fullgrad=np.zeros((full_dimension,3)) @@ -539,6 +542,8 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el print("Gradient (Combined):", self.gradient) if Grad: + if PC and pc_gradient is not None: + return self.energy, self.gradient, pc_gradient return self.energy, self.gradient else: return self.energy