- ONIOMTheory: linkatom changes, RCD option, ONIOM-2 gradient complet…

…e, ONIOM-3 gradient not ready

- QMMMTheory: linkatom changes, RCD shifting,
- calculate_RMSD function for 2 fragments, allows subset match and heavyatoms option
- get_linkatom_positions function extension: simple and ratio method
- grabatomcharges_ORCA: bugfix for BS job,
- write freqs and normal mode comps to file by default
- various cleanup

ash/__init__.py

@@ -31,7 +31,7 @@
     angle_between_atoms, dihedral_between_atoms, pdb_to_smiles, xyz_to_pdb_with_connectivity, writepdb_with_connectivity, mol_to_pdb, sdf_to_pdb
 from .modules.module_coords import remove_atoms_from_system_CHARMM, add_atoms_to_system_CHARMM, getwaterconstraintslist,\
     QMregionfragexpand, QMPC_fragexpand, read_xyzfiles, Reaction, define_XH_constraints, simple_get_water_constraints, print_internal_coordinate_table,\
-    flexible_align_pdb, flexible_align_xyz, flexible_align, insert_solute_into_solvent, nuc_nuc_repulsion
+    flexible_align_pdb, flexible_align_xyz, flexible_align, insert_solute_into_solvent, nuc_nuc_repulsion, calculate_RMSD
 
 # Singlepoint
 import ash.modules.module_singlepoint

ash/interfaces/interface_MLatom.py

@@ -9,7 +9,7 @@
 ##########################
 # MLatom Theory interface
 ##########################
-# NOTE: we only intend to support mostly ML functionality in MLatom (not all basic QM methods)
+# NOTE: we mostly intend to support - ML functionality in MLatom (not all basic QM methods)
 
 # MLatom has 3 types of models:
 # 1) methods: these are pre-trained or at least standalone electronic structure methods
@@ -56,8 +56,6 @@ def __init__(self, printlevel=2, numcores=1, label="mlatom",
             print("Neither a method or ml_model was selected for MLatomTheory interface. Exiting.")
             ashexit()
 
-
-
         # METHODS: pre-trained models
         # Note: useful method keywords in MLAatom below
         # AIQMx models require either MNDO or Sparrow. Also dftd4 (except AIQM1@DFT* but that one is bad anyway)
@@ -243,7 +241,8 @@ def train(self, molDB_xyzfile=None, molDB_scalarproperty_file=None,
                 if 'lambda' in hyperparameters:
                     print('Optimized lambda:', lmbd)
                 print('Optimized validation loss:', valloss)
-
+            else:
+                print("No hyperparameters provided (NOT recommended)")
             print("\nNow training...")
             self.model.train(molecular_database=molDB,
                             property_to_learn=property_to_learn,

ash/interfaces/interface_ORCA.py

@@ -1954,6 +1954,10 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
     charges=[]
     BS=False #if broken-symmetry job
     column=None
+
+    numatoms = int(pygrep('Number of atoms                             ...', outputfile)[-1])
+    print("numatoms:", numatoms)
+
     #if
     if len(pygrep2("WARNING: Broken symmetry calculations", outputfile)):
         BS=True
@@ -1970,7 +1974,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
                         charges.append(float(line.split()[2]))
                 if 'Atom No    Charge        Core      Valence    Rydberg      Total' in line:
                     grab=True
-    elif chargemodel=="CHELPG":
+    elif chargemodel.upper() =="CHELPG":
         with open(outputfile) as ofile:
             for line in ofile:
                 if grab==True:
@@ -1982,7 +1986,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
                     grab=True
                     #Setting charges list to zero in case of multiple charge-tables. Means we grab second table
                     charges=[]
-    elif chargemodel=="Hirshfeld":
+    elif chargemodel.upper() =="HIRSHFELD":
         with open(outputfile) as ofile:
             for line in ofile:
                 if grab==True:
@@ -1994,7 +1998,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
                     grab=True
                     #Setting charges list to zero in case of multiple charge-tables. Means we grab second table
                     charges=[]
-    elif chargemodel=="CM5":
+    elif chargemodel.upper()=="CM5":
         elems = []
         coords = []
         with open(outputfile) as ofile:
@@ -2023,7 +2027,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
         atomicnumbers=ash.modules.module_coords.elemstonuccharges(elems)
         charges = ash.functions.functions_elstructure.calc_cm5(atomicnumbers, coords, charges)
         print("CM5 charges :", list(charges))
-    elif chargemodel == "Mulliken":
+    elif chargemodel.upper() == "MULLIKEN":
         with open(outputfile) as ofile:
             for line in ofile:
                 if grab==True:
@@ -2038,7 +2042,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
                     else:
                         column=-1
 
-    elif chargemodel == "Loewdin":
+    elif chargemodel.upper() == "LOEWDIN":
         with open(outputfile) as ofile:
             for line in ofile:
                 if grab==True:
@@ -2054,7 +2058,7 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
                         column=-2
                     else:
                         column=-1
-    elif chargemodel == "IAO":
+    elif chargemodel.upper() == "IAO":
         with open(outputfile) as ofile:
             for line in ofile:
                 if grab==True:
@@ -2073,8 +2077,9 @@ def grabatomcharges_ORCA(chargemodel,outputfile):
     #If BS then we have grabbed charges for both high-spin and BS solution
     if BS is True:
         print("Broken-symmetry job detected. Only taking BS-state populations")
-        charges=charges[int(len(charges)/2):]
-
+        if len(charges) != numatoms:
+            charges=charges[numatoms:]
+        print("charges:", charges)
     return charges
 
 
@@ -2552,16 +2557,17 @@ def create_ASH_otool(basename=None, theoryfile=None, scriptlocation=None, charge
         #otool.write("theory=ZeroTheory()\n")
         #otool.write("theory=ZeroTheory()\n")
         otool.write("result=Singlepoint(theory=theory,fragment=frag,Grad=True, charge={}, mult={})\n".format(charge,mult))
-        otool.write("energy = result.energy")
-        otool.write("gradient = result.gradient")
+        otool.write("energy = result.energy\n")
+        otool.write("gradient = result.gradient\n")
         otool.write("print(gradient)\n")
         otool.write("ash.interfaces.interface_ORCA.print_gradient_in_ORCAformat(energy,gradient,\"{}\")\n".format(basename))
     st = os.stat(scriptlocation+"/otool_external")
     os.chmod(scriptlocation+"/otool_external", st.st_mode | stat.S_IEXEC)
 
 # Using ORCA as External Optimizer for ASH
 #Will only work for theories that can be pickled: not OpenMMTheory, probably not QMMMTheory
-def ORCA_External_Optimizer(fragment=None, theory=None, orcadir=None, charge=None, mult=None):
+def ORCA_External_Optimizer(fragment=None, theory=None, orcadir=None, charge=None, mult=None,
+                            ORCA_jobkeyword="Opt"):
     print_line_with_mainheader("ORCA_External_Optimizer")
     if fragment == None or theory == None:
         print("ORCA_External_Optimizer requires fragment and theory keywords")
@@ -2607,7 +2613,7 @@ def ORCA_External_Optimizer(fragment=None, theory=None, orcadir=None, charge=Non
 
     #ORCA input file
     with open(basename+".inp", 'w') as o:
-        o.write("! ExtOpt Opt\n")
+        o.write(f"! ExtOpt {ORCA_jobkeyword}\n")
         o.write("\n")
         o.write("*xyzfile {} {} {}\n".format(charge,mult,xyzfile))
 
@@ -2631,7 +2637,10 @@ def ORCA_External_Optimizer(fragment=None, theory=None, orcadir=None, charge=Non
     fragment.coords=coords
 
     #Grabbing final energy
-    energy = ORCAfinalenergygrab(basename+".out")
+    #energy = ORCAfinalenergygrab(basename+".out")
+    #energy = pygrep(,f"{basename}.out")
+    energylines = pygrep2("FINAL SINGLE POINT ENERGY (From external program)",f"{basename}.out", errors="ignore")
+    energy = float(energylines[-1].split()[-1])
     print("Final energy from external ORCA optimization:", energy)
 
     return energy

ash/interfaces/interface_OpenMM.py

@@ -150,7 +150,7 @@ def __init__(self, printlevel=2, platform='CPU', numcores=1, topoforce=False, fo
         # Platform (CPU, CUDA, OpenCL) and Parallelization
         self.platform_choice = platform
 
-        if properties == None:
+        if properties is None:
             self.properties = {}
         else:
             self.properties = properties
@@ -5331,16 +5331,16 @@ def trim_list_of_lists(k):
 def merge_pdb_files(pdbfile_1,pdbfile_2,outputname="merged.pdb"):
     import openmm.app
 
-    #Funciton to merge PDB-files (e.g. protein and ligand) while preserving and updating connectivity records
+    #Function to merge PDB-files (e.g. protein and ligand) while preserving and updating connectivity records
     #PDB inputfiles
     pdb1 = openmm.app.PDBFile(pdbfile_1)
     pdb2 = openmm.app.PDBFile(pdbfile_2)
 
     #Create modeller object
     modeller = openmm.app.Modeller(pdb1.topology, pdb1.positions) #Add pdbfile1
     modeller.add(pdb2.topology, pdb2.positions) #Add pdbfile2
-    mergedTopology = modeller.topology #Merging topology files
-    mergedPositions = modeller.positions #merging ositions
+    mergedTopology = modeller.topology #Merging topology
+    mergedPositions = modeller.positions #merging positions
 
     #Write merged topology and positions to new PDB file
     write_pdbfile_openMM(modeller.topology, mergedPositions, outputname)

ash/interfaces/interface_knarr.py

@@ -342,7 +342,7 @@ def NEB(reactant=None, product=None, theory=None, images=8, CI=True, free_end=Fa
     elif runmode == 'serial' and numcores == 1:
         print(BC.WARNING,"NEB runmode is serial, i.e. running one image after another.", BC.END)
         print("theory:", theory)
-        print("theory.numcors:", theory.numcores)
+        print("Theory.numcores:", theory.numcores)
         if theory.numcores > 1:
             print(BC.WARNING,f"Theory parallelization is active and will utilize: {theory.numcores} CPU cores per image.",BC.END)
         else: