Reattempt at committing

molSimplify/Classes/mol3D.py

@@ -6471,3 +6471,96 @@ def get_molecular_mass(self):
             self.mass = mol_mass
 
         return mol_mass
+
+    def mol3D_to_networkx(self,get_symbols:bool=True,get_bond_order:bool=True,get_bond_distance:bool=False):
+        g = nx.Graph()
+        # get every index of atoms in mol3D object
+        for atom_ind in range(0,self.natoms):
+            # set each atom as a node in the graph, and add symbol information if wanted
+            data={}
+            if get_symbols:
+                data['Symbol']=self.getAtom(atom_ind).symbol()
+                data['atom3D']=self.getAtom(atom_ind)
+            g.add_node(atom_ind,**data)
+        # get every bond in mol3D object
+        bond_info=self.bo_dict
+        for bond in bond_info:
+            # set each bond as an edge in the graph, and add symbol information if wanted
+            data={}
+            if get_bond_order:
+                data['bond_order']=bond_info[bond]
+            if get_bond_distance:
+                distance=self.getAtom(bond[0]).distance(self.getAtom(bond[1]))
+                data['bond_distance']=distance
+            g.add_edge(bond[0],bond[1],**data)
+        return g
+
+    def roland_combine(self, mol, catoms, bond_to_add=[], dirty=False):
+        """
+        Combines two molecules. Each atom in the second molecule
+        is appended to the first while preserving orders. Assumes
+        operation with a given mol3D instance, when handed a second mol3D instance.
+
+        Parameters
+        ----------
+            mol : mol3D
+                mol3D class instance containing molecule to be added.
+            bond_to_add : list, optional
+                List of tuples (ind1,ind2,order) bonds to add. Default is empty.
+            dirty : bool, optional
+                Add atoms without worrying about bond orders. Default is False.
+
+        Returns
+        -------
+            cmol : mol3D
+                New mol3D class containing the two molecules combined.
+        """
+
+        cmol = self
+        bo_dict = cmol.bo_dict
+
+        print('lig_dict')
+        print(mol.bo_dict)
+
+        if cmol.bo_dict == False:
+            # only central metal
+            bo_dict = {}
+        new_bo_dict = copy.deepcopy(bo_dict)
+
+        # add ligand connections
+        for bo in mol.bo_dict:
+            ind1 = bo[0] + len(cmol.atoms)
+            ind2 = bo[1] + len(cmol.atoms)
+            new_bo_dict[(ind1,ind2)]=mol.bo_dict[(bo[0],bo[1])]
+
+        # connect metal to ligand
+        metal_ind = cmol.findMetal()[0]
+        for atom in catoms:
+            ind1 = metal_ind
+            ind2 = atom + len(cmol.atoms)
+            new_bo_dict[(ind1,ind2)]=1
+
+        for atom in mol.atoms:
+            cmol.addAtom(atom, auto_populate_BO_dict = False)
+
+        cmol.bo_dict = new_bo_dict
+        return cmol
+
+    def graph_from_bodict(self, bo_dict):
+        g = []
+        for i, atom in enumerate(self.atoms):
+            sub_g = []
+            connected = []
+            for tup in bo_dict:
+                if i in tup:
+                    if i != tup[0]:
+                        connected.append(tup[0])
+                    else:
+                        connected.append(tup[1])
+            for j in range(0,len(self.atoms)):
+                if j in connected:
+                    sub_g.append(1.)
+                else:
+                    sub_g.append(0.)
+            g.append(sub_g)
+        return g

molSimplify/Scripts/structgen.py

@@ -12,6 +12,7 @@
     from openbabel import openbabel  # version 3 style import
 except ImportError:
     import openbabel  # fallback to version 2
+from openbabel import pybel
 import random
 import itertools
 import numpy as np
@@ -25,6 +26,7 @@
                                           getPointu,
                                           kabsch,
                                           midpt,
+                                          move_point,
                                           norm,
                                           PointTranslateSph,
                                           reflect_through_plane,
@@ -2306,6 +2308,7 @@
     emsg = ''
     complex3D: List[mol3D] = []
     occs0 = []      # occurrences of each ligand
+    complex2D = []
     toccs = 0       # total occurrence count (number of ligands)
     smilesligs = 0  # count how many smiles strings
     cats0: List[List[Union[int, str]]] = []  # connection atoms for ligands
@@ -2507,6 +2510,8 @@
                         lig = decorate_ligand(
                             lig, args.decoration[i], args.decoration_index[i], args.debug)
             lig.convert2mol3D()
+
+
             # initialize ligand
             lig3D, rempi, ligpiatoms = init_ligand(args, lig, tcats, keepHs, i)
             if emsg:
@@ -2710,16 +2715,34 @@
                 auxm = mol3D()
                 auxm.copymol3D(lig3D)
                 complex3D.append(auxm)
+
+                lig3D_copy = mol3D()
+                lig3D_copy.copymol3D(lig3D)
+                lig3D_copy.populateBOMatrix(bonddict=True)
+                lig2D = lig3D_copy.mol3D_to_networkx()
+                complex2D.append(lig2D)
+
                 if 'a' not in lig.ffopt.lower():
                     for latdix in range(0, lig3D.natoms):
                         if args.debug:
                             print(('a is not ff.lower, so adding atom:  ' +
                                    str(latdix+core3D.natoms) + ' to freeze'))
                         frozenats.append(latdix+core3D.natoms)
+
+
+
+                # combine molecules
+                if len(core3D.atoms) == 1:
+                    core3D_copy = mol3D()
+                    core3D_copy.copymol3D(core3D)
+                core3D_copy = core3D_copy.roland_combine(lig3D_copy, catoms)
+
+
                 # combine molecules
                 core3D = core3D.combine(lig3D)
                 core3D.convert2OBMol()
                 core3D.convert2mol3D()
+
                 # remove dummy cm atom if requested
                 if rempi:
                     # remove the fictitious center atom, for aromatic-bonding ligands like benzene
@@ -2792,6 +2815,8 @@
             core3D.writexyz('complex_after_final_ff.xyz')
         # core3D,enc = ffopt(args.ff,core3D,connected,1,frozenats,freezeangles,MLoptbds,'Adaptive',args.debug)
 
+    core3D.bo_dict = core3D_copy.bo_dict
+
     return core3D, complex3D, emsg, this_diag, subcatoms_ext, mligcatoms_ext
 
 
@@ -2924,7 +2949,7 @@
 
 
 def structgen(args: Namespace, rootdir: str, ligands: List[str], ligoc: List[int],
-              sernum: int, write_files: bool = True) -> Tuple[List[str], str, run_diag]:
+              sernum: int, write_files: bool = True, smart_generation: bool = True) -> Tuple[List[str], str, run_diag]:
     """Main structure generation routine - multiple structures
 
     Parameters
@@ -3016,6 +3041,128 @@
     # write xyz file
     if (not args.reportonly) and (write_files):
         core3D.writexyz(fname, no_tabs=args.no_tabs)
+        core3D.writemol2(fname)
+
+        if smart_generation == True:
+            # optimize
+            metal_ind = core3D.findMetal()[0]
+            freeze_inds = []
+            for bond in core3D.bo_dict:
+                if metal_ind in bond:
+                    freeze_inds.append(bond[0]+1)
+                    freeze_inds.append(bond[1]+1)
+            freeze_inds = list(set(list(freeze_inds)))
+
+            obConversion = openbabel.OBConversion()
+            OBMol = openbabel.OBMol()
+            constraints = openbabel.OBFFConstraints()
+            obConversion.SetInAndOutFormats("mol2", "mol2")
+            obConversion.ReadString(OBMol, core3D.writemol2('',writestring = True))
+            for atom in freeze_inds:
+                constraints.AddAtomConstraint(atom)
+            ff = pybel._forcefields["uff"]
+            success = ff.Setup(OBMol, constraints)
+            ff.SetConstraints(constraints)
+            if success:
+                ff.ConjugateGradients(10000)
+                ff.GetCoordinates(OBMol)
+            obConversion.WriteFile(OBMol,fname+'BABEL.mol2')
+            obConversion.SetOutFormat("xyz")
+            obConversion.WriteFile(OBMol,fname+'BABEL.xyz')
+
+            # check if bad
+            mol = mol3D()
+            mol.readfrommol2(fname+'BABEL.mol2')
+            overlap, mind = mol.sanitycheck(silence = True)
+            if overlap:
+                mind = 1000
+                errors_dict = {}
+                for ii, atom1 in enumerate(mol.atoms):
+                    for jj, atom0 in enumerate(mol.atoms):
+                        if jj > ii:
+                            if atom1.ismetal() or atom0.ismetal():
+                                cutoff = 0.6
+                            elif (atom0.sym in ['N', 'O'] and atom1.sym == 'H') or (atom1.sym in ['N', 'O'] and atom0.sym == 'H'):
+                                cutoff = 0.6
+                            else:
+                                cutoff = 0.65
+                            if distance(atom1.coords(), atom0.coords()) < cutoff * (atom1.rad + atom0.rad):
+                                norm = distance(
+                                    atom1.coords(), atom0.coords())/(atom1.rad+atom0.rad)
+                                errors_dict.update(
+                                    {atom1.sym + str(ii)+'-'+atom0.sym+str(jj)+'_normdist': norm})
+                                if distance(atom1.coords(), atom0.coords()) < mind:
+                                    mind = distance(atom1.coords(), atom0.coords())
+                                    if mind == 0.0:
+                                        # move atom0 over a little bit
+                                        atom0.setcoords(np.array(atom1.coords())+0.02)
+                                        obConversion.SetInAndOutFormats("mol2", "mol2")
+                                        OBMol = openbabel.OBMol()
+                                        obConversion.ReadString(OBMol, mol.writemol2('',writestring = True))
+
+                ff = pybel._forcefields["gaff"]
+                success = ff.Setup(OBMol, constraints)
+                ff.SetConstraints(constraints)
+                if success:
+                    ff.ConjugateGradients(10000)
+                    ff.GetCoordinates(OBMol)
+                    ff = pybel._forcefields["uff"]
+                    success = ff.Setup(OBMol, constraints)
+                    ff.SetConstraints(constraints)
+                    if success:
+                        ff.ConjugateGradients(10000)
+                        ff.GetCoordinates(OBMol)
+
+
+            obConversion.SetOutFormat("mol2")
+            obConversion.WriteFile(OBMol,fname+'BABEL.mol2')
+            obConversion.SetOutFormat("xyz")
+            obConversion.WriteFile(OBMol,fname+'BABEL.xyz')
+
+            # check if overextended H:
+            mol = mol3D()
+            mol.readfrommol2(fname+'BABEL.mol2')
+            changed = False
+            for bond in mol.bo_dict:
+                atom0 = mol.atoms[bond[0]]
+                atom1 = mol.atoms[bond[1]]
+                dist = -100000000000
+                if atom0.sym == 'C' and atom1.sym == 'H':
+                    dist = atom0.distance(atom1)
+                    L1 = np.array(tuple(atom0.coords()))
+                    L2 = np.array(tuple(atom1.coords()))
+                    if dist > 1.15:
+                        vector = L1 - L2
+                        required_dist = dist - 1.15
+                        new_point = move_point(atom1.coords(), vector, required_dist)
+                        atom1.setcoords(new_point)
+                        changed = True
+                elif atom0.sym == 'H' and atom1.sym == 'C':
+                    dist = atom0.distance(atom1)
+                    if dist > 1.15:
+                        L1 = np.array(tuple(atom0.coords()))
+                        L2 = np.array(tuple(atom1.coords()))
+                        vector = L2 - L1
+                        required_dist = dist - 1.15
+                        new_point = move_point(atom0.coords(), vector, required_dist)
+                        atom0.setcoords(new_point)
+                        changed = True
+            if changed:
+                mol.writemol2(fname+'BABEL.mol2')
+                obConversion = openbabel.OBConversion()
+                OBMol = openbabel.OBMol()
+                obConversion.SetInAndOutFormats("mol2", "mol2")
+                obConversion.ReadFile(OBMol, fname+'BABEL.mol2')
+                ff = pybel._forcefields["uff"]
+                success = ff.Setup(OBMol, constraints)
+                ff.SetConstraints(constraints)
+                if success:
+                    ff.ConjugateGradients(10000)
+                    ff.GetCoordinates(OBMol)
+                obConversion.WriteFile(OBMol,fname+'BABEL.mol2')
+                obConversion.SetOutFormat("xyz")
+                obConversion.WriteFile(OBMol,fname+'BABEL.xyz')
+
     strfiles.append(fname)
     if write_files:
         # write report file