Fixes to gcp_gradient and WrapTheory pc_gradient

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]

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