Merge branch 'master' into NEW

Author: RagnarB83

ash/interfaces/interface_DFTB.py

@@ -266,6 +266,7 @@ def write_DFTB_input(hamiltonian,xtbmethod,xyzfilename, elems,coords,charge,mult
     #############
     method1=f"Hamiltonian = {hamiltonian}" +"{"+"\n"
     inputlines.append(method1)
+    inputlines.append(f"Charge = {charge}\n")
     if 'XTB' in hamiltonian.upper():
         method2=f"Method = '{xtbmethod}'"+'\n\n'
         inputlines.append(method2)
@@ -417,9 +418,12 @@ def grab_energy_gradient_DFTB(file, numatoms, Grad=False, PC=False, numpc=None):
                 if 'Dipole moment:' in line:
                     pcgrab=False
 
-    if Grad is True and PC is True:
-        if np.any(pcgradient) is False:
-            print("Error: PCgradient array from DFTB output is zero. Something went wrong.")
+    if Grad is True:
+        if not np.any(gradient):
+            print("Error: gradient array from DFTB output is zero. Something went wrong (e.g. SCC not converged or DFTB+ crashed before writing forces). Check DFTB outputfile.")
+            ashexit()
+        if PC is True and not np.any(pcgradient):
+            print("Error: PCgradient array from DFTB output is zero. Something went wrong (e.g. SCC not converged or DFTB+ crashed before writing forces). Check DFTB outputfile.")
             ashexit()
 
     return energy,gradient,pcgradient

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/interfaces/interface_OpenMM.py

@@ -2695,8 +2695,6 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi
             xmlfile = "amber10.xml"
         elif forcefield == 'Amber14':
             xmlfile = "amber14-all.xml"
-        elif forcefield == 'Amber96':
-            xmlfile = "amber96.xml"
         elif forcefield == 'CHARMM36':
             xmlfile = "charmm36.xml"
         elif forcefield == 'CHARMM2013':
@@ -2708,53 +2706,45 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi
         else:
             print("Unknown forcefield")
             exit()
-    elif xmlfile is not None:
-        print("Using xmlfile:", xmlfile)
-    elif forcefield_object is not None:
-        pass
-    else:
-        print("You must provide a forcefield name, forcefieldobject or xmlfile keywords!")
-        ashexit()
-
-    #Water model selection for CHARMM forcefields
-    if 'CHARMM' in forcefield:
-        # Using specific CHARMM36 version of TIP3P
-        if watermodel is None:
-            print("No watermodel selected.")
-            if waterxmlfile is None:
-                print("No waterxmlfile selected either")
-                print("Selecting automatically recommended CHARMM-style TIP3P")
-                watermodel='tip3p'
-
-        print("watermodel:", watermodel)
-        if watermodel.lower() == "tip3p":
-            modeller_solvent_name="tip3p" #Used when adding solvent
-            waterxmlfile = "charmm36/water.xml"
-        print("Waterxmlfile selected:", waterxmlfile)
-
-
-    #Water model selection for AMber forcefields
-    if 'Amber' in forcefield:
-        if watermodel is None:
-            print("No watermodel selected.")
-            if waterxmlfile is None:
-                print("No waterxmlfile selected either")
-                print("Selecting automatically recommended TIP3P-4B (watermodel='tip3pfb')")
-                print("This is a reparameterized version of TIP3P")
-                watermodel='tip3pfb'
-        print("watermodel:", watermodel)
-        # Using specific Amber FB version of TIP3P
-        if watermodel.lower() == "tip3pfb" or watermodel.lower() == "tip3p-fb":
-            modeller_solvent_name="tip3p" #Used when adding solvent
-            waterxmlfile = "amber14/tip3pfb.xml" #NOTE: this is not actually TIP3P but a reparaterized version
-        elif watermodel.lower () == 'tip3p':
-                modeller_solvent_name="tip3p"
-                if forcefield == 'Amber14':
-                    waterxmlfile = "amber14/tip3p.xml"
-                else:
-                    waterxmlfile = "tip3p.xml"
-        print("Waterxmlfile selected:", waterxmlfile)
 
+        #Water model selection for CHARMM forcefields
+        if 'CHARMM' in forcefield:
+            # Using specific CHARMM36 version of TIP3P
+            if watermodel is None:
+                print("No watermodel selected.")
+                if waterxmlfile is None:
+                    print("No waterxmlfile selected either")
+                    print("Selecting automatically recommended CHARMM-style TIP3P")
+                    watermodel='tip3p'
+
+            print("watermodel:", watermodel)
+            if watermodel.lower() == "tip3p":
+                modeller_solvent_name="tip3p" #Used when adding solvent
+                waterxmlfile = "charmm36/water.xml"
+            print("Waterxmlfile selected:", waterxmlfile)
+
+
+        #Water model selection for AMber forcefields
+        if 'Amber' in forcefield:
+            if watermodel is None:
+                print("No watermodel selected.")
+                if waterxmlfile is None:
+                    print("No waterxmlfile selected either")
+                    print("Selecting automatically recommended TIP3P-4B (watermodel='tip3pfb')")
+                    print("This is a reparameterized version of TIP3P")
+                    watermodel='tip3pfb'
+            print("watermodel:", watermodel)
+            # Using specific Amber FB version of TIP3P
+            if watermodel.lower() == "tip3pfb" or watermodel.lower() == "tip3p-fb":
+                modeller_solvent_name="tip3p" #Used when adding solvent
+                waterxmlfile = "amber14/tip3pfb.xml" #NOTE: this is not actually TIP3P but a reparaterized version
+            elif watermodel.lower () == 'tip3p':
+                    modeller_solvent_name="tip3p"
+                    if forcefield == 'Amber14':
+                        waterxmlfile = "amber14/tip3p.xml"
+                    else:
+                        waterxmlfile = "tip3p.xml"
+            print("Waterxmlfile selected:", waterxmlfile)
 
     ############
     # Define a forcefield if using XML-files
@@ -2778,9 +2768,18 @@ def OpenMM_Modeller(pdbfile=None, forcefield_object=None, forcefield=None, xmlfi
             forcefield_obj = openmm_app.forcefield.ForceField(xmlfile,waterxmlfile)
         elif extraxmlfile is not None and waterxmlfile is not None:
             forcefield_obj = openmm_app.forcefield.ForceField(xmlfile,extraxmlfile,waterxmlfile)
-    else:
+
+    elif forcefield_object is not None:
         print("Using forcefield object provided")
         forcefield_obj= forcefield_object
+        
+        if watermodel is not None or waterxmlfile is not None:
+            print("Warning: watermodel/waterxmlfile ignored when forcefield_object is supplied")
+
+
+    else:
+        print("You must provide a forcefield name, forcefieldobject or xmlfile keywords!")
+        ashexit()
 
     print("PDBfile:", pdbfile)
     print("pH:", pH)

ash/interfaces/interface_dlfind.py

@@ -526,7 +526,7 @@ def setup_constraints_act_frozen(self):
                 self.frozenatoms = self.frozenatoms+frozenatoms_xyz
                 print("frozenatoms_z:", frozenatoms_z)
         if self.actatoms is not None:
-            print_if_level("Actatoms provided: {self.actatoms}", self.printlevel,2)
+            print_if_level(f"Actatoms provided: {self.actatoms}", self.printlevel, 2)
 
             if self.PBC:
                 print("PBC detected. Adding 4 dummy atoms to actatoms if not already present")
@@ -937,4 +937,4 @@ def define_residues(fragment=None, min_size=5, max_size=15):
         else:
             residues.append(current_res)
 
-    return residues
+    return residues

ash/interfaces/interface_xtb.py

@@ -121,7 +121,8 @@ def __init__(self, xtbdir=None, xtbmethod='GFN1', runmode='inputfile', numcores=
         self.analytic_hessian=False
         print_line_with_mainheader(f"{self.theorynamelabel}Theory initialization")
 
-        # Hardness of pointcharge. GAM factor. Big number means PC behaviour
+        # Hardness of pointcharge. GAM factor. Big number means PC behaviour. 
+        # If hardness is set to 'elements', xtb's element-specific hardness is used (see https://xtb-docs.readthedocs.io/en/latest/pcem.html).
         self.hardness=hardness_PC
 
         # Accuracy (0.1 it quite tight)
@@ -521,7 +522,11 @@ def run(self, current_coords=None, current_MM_coords=None, MMcharges=None, qm_el
                 print("...")
             # Create pcharge file if PC
             if PC:
-                create_xtb_pcfile_general(current_MM_coords, MMcharges, hardness=self.hardness)
+                if self.hardness == 'elements':
+                    hardness = mm_elems
+                else:
+                    hardness = [self.hardness] * len(MMcharges)
+                create_xtb_pcfile_general(current_MM_coords, MMcharges, hardness=hardness)
 
             # Run xTB (note: passing PC and Grad Booleans)
             run_xtb_SP(self.xtbdir, self.xtbmethod, coordfile, charge, mult, printlevel=self.printlevel, PC=PC, solvent=self.solvent,
@@ -1059,8 +1064,8 @@ def create_xtb_pcfile_general(coords,pchargelist,hardness=1000):
     #https://xtb-docs.readthedocs.io/en/latest/pcem.html
     with open('pcharge', 'w') as pcfile:
         pcfile.write(str(len(pchargelist))+'\n')
-        for p,c in zip(pchargelist,coords):
-            line = "{} {} {} {} {}".format(p, c[0], c[1], c[2], hardness)
+        for p,c,h in zip(pchargelist,coords,hardness):
+            line = "{} {} {} {} {}".format(p, c[0], c[1], c[2], h)
             pcfile.write(line+'\n')
 
 

ash/modules/module_QMMM.py

@@ -985,7 +985,8 @@ def mech_run(self, current_coords=None, elems=None, Grad=False, numcores=1, exit
 
                 self.MMenergy, self.MMgradient = self.mm_theory.run(current_coords=current_coords, qmatoms=self.qmatoms, Grad=True)
             else:
-                print("QM/MM Grad is false")
+                if self.printlevel >= 2:
+                    print("QM/MM Grad is false")
                 self.MMenergy = self.mm_theory.run(current_coords=current_coords, qmatoms=self.qmatoms)
         else:
             self.MMenergy=0
@@ -1404,7 +1405,8 @@ def elstat_run(self, current_coords=None, elems=None, Grad=False, numcores=1, ex
 
                 self.MMenergy, self.MMgradient= self.mm_theory.run(current_coords=current_coords, qmatoms=self.qmatoms, Grad=True)
             else:
-                print("QM/MM Grad is false")
+                if self.printlevel >= 2:
+                    print("QM/MM Grad is false")
                 self.MMenergy= self.mm_theory.run(current_coords=current_coords, qmatoms=self.qmatoms)
         else:
             self.MMenergy=0

ash/modules/module_coords.py

@@ -3940,13 +3940,15 @@ def check_charge_mult(charge, mult, theorytype, fragment, jobtype, theory=None,
 
         #Note: theory needs to be set
         if charge is None or mult is None:
-            print(BC.WARNING,f"Warning: Charge/mult was not provided to {jobtype}",BC.END)
-            print("Checking if present in QM/MM object")
+            if printlevel >= 1:
+                print(BC.WARNING,f"Warning: Charge/mult was not provided to {jobtype}",BC.END)
+                print("Checking if present in QM/MM object")
             if theory.qm_charge != None and theory.qm_mult != None:
-                print("Found qm_charge and qm_mult attributes.")
                 charge=theory.qm_charge
                 mult=theory.qm_mult
-                print(f"Using charge={charge} and mult={mult}")
+                if printlevel >= 1:
+                    print("Found qm_charge and qm_mult attributes.")
+                    print(f"Using charge={charge} and mult={mult}")
             elif fragment.charge != None and fragment.mult != None:
                 print(BC.WARNING,"Fragment contains charge/mult information: Charge: {} Mult: {} Using this instead".format(fragment.charge,fragment.mult), BC.END)
                 print(BC.WARNING,"Make sure this is what you want!", BC.END)

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

ash/modules/module_surface_new.py

@@ -91,6 +91,7 @@ def calc_surface(
                     'subfrctor': subfrctor,
                     'force_noPBC': force_noPBC, 'PBC_format_option': PBC_format_option,
                     'ActiveRegion': ActiveRegion, 
+                    'actatoms': actatoms, 
                     'result_write_to_disk':False,
                     'printlevel':printlevel,
                     }
@@ -154,7 +155,14 @@ def calc_surface(
         ashexit()
 
     # Build connectivity once
-    conn = _build_connectivity(fragment.coords, fragment.elems)
+    # Build connectivity once
+    if ActiveRegion is True and actatoms is not None and len(actatoms) > 0:
+        print(f"Building connectivity only inside active region: {len(actatoms)} atoms")
+        conn = _build_connectivity(fragment.coords, fragment.elems,
+                               atom_indices=actatoms)
+    else:
+        print(f"Building connectivity for full fragment: {len(fragment.coords)} atoms")
+        conn = _build_connectivity(fragment.coords, fragment.elems)
 
     # Changing printlevel of fragment
     fragment.printlevel=printlevel
@@ -406,7 +414,7 @@ def calc_surface(
                         optimizerobj.Hess_inv = None
 
                 #Running optimizer object, passing theory, fragment, constraints and possible extra kws
-                result = optimizerobj.run(theory=theory,fragment=fragment, constraints=allconstraints, **extraoopt_run_kws)
+                result = optimizerobj.run(theory=theory,fragment=fragment, constraints=allconstraints, charge=charge, mult=mult, **extraoopt_run_kws)
 
                 #if pointcount == 2:
                 #    print("2nd point optimization result:", result)
@@ -1150,21 +1158,41 @@ def _measure_dihedral(coords, i, j, k, l):
 # Connectivity
 # ---------------------------------------------------------------------------
 
-def _build_connectivity(coords, elems):
+def _build_connectivity(coords, elems, atom_indices=None):
     coords = np.asarray(coords)
     n = len(elems)
+
     radii = np.array([
-        _COVALENT_RADII.get(e.capitalize(), _DEFAULT_RADIUS) for e in elems
+        _COVALENT_RADII.get(e.capitalize(), _DEFAULT_RADIUS)
+        for e in elems
     ])
+
+    # Keep full-length connectivity list so downstream code
+    # can continue using global atom indices
     conn = [set() for _ in range(n)]
-    for i in range(n):
-        for j in range(i + 1, n):
+
+    # Default behaviour: full-system connectivity
+    if atom_indices is None:
+        atom_indices = range(n)
+    else:
+        atom_indices = list(atom_indices)
+
+    nsel = len(atom_indices)
+
+    for a in range(nsel):
+        i = atom_indices[a]
+
+        for b in range(a + 1, nsel):
+            j = atom_indices[b]
+
             dist = np.linalg.norm(coords[i] - coords[j])
             threshold = radii[i] + radii[j] + _CONNECTIVITY_TOLERANCE
-            # Ignore very short distances (e.g. same atom or ghost atoms)
+
+            # Ignore very short distances
             if 0.4 < dist < threshold:
                 conn[i].add(j)
                 conn[j].add(i)
+
     return conn
 
 def _atoms_on_side(start, fixed, conn):
@@ -1488,7 +1516,7 @@ def _preset_geometry_restraint(fragment, RC_list, rc_values, optimizerobj,
     preset_args = {k: v for k, v in opt_arguments.items()
                    if k not in ('constraints', 'constrainvalue')}
     # Optimizing with restraint theory, passing extraconstraints as contraints if present
-    optimizerobj.run(theory=restraint_theory,fragment=fragment, constraints=extraconstraints, **extraoopt_run_kws)
+    optimizerobj.run(theory=restraint_theory,fragment=fragment, constraints=extraconstraints, charge=charge, mult=mult, **extraoopt_run_kws)
 
     #optimizer(
     #    fragment=fragment, theory=restraint_theory, constraints=extraconstraints,