diff --git a/molSimplify/Classes/mol3D.py b/molSimplify/Classes/mol3D.py index 5e1e1da2..ea105983 100644 --- a/molSimplify/Classes/mol3D.py +++ b/molSimplify/Classes/mol3D.py @@ -10,6 +10,7 @@ import sys import tempfile import time +import copy import xml.etree.ElementTree as ET import numpy as np import networkx as nx @@ -6471,3 +6472,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 diff --git a/molSimplify/Scripts/generator.py b/molSimplify/Scripts/generator.py index 6429892f..d9290f01 100644 --- a/molSimplify/Scripts/generator.py +++ b/molSimplify/Scripts/generator.py @@ -236,7 +236,7 @@ def startgen(argv, flag, gui, inputfile_str=None, write_files=True): print(('adding ' + str(args.ligadd) + ' to ligand database with name ' + args.ligname + ' and connection atom(s) ' + str(args.ligcon))) addtoldb(smimol=args.ligadd, sminame=args.ligname, smident=len(args.ligcon), - smicat=str(args.ligcon).strip('[]'), smigrps="custom", smictg="custom", ffopt=args.ligffopt) + smicat=str(args.ligcon).strip('[]'), smigrps="custom", smictg="custom", ffopt=args.ligffopt, overwrite=args.overwrite) # normal structure generation or transition state building else: diff --git a/molSimplify/Scripts/inparse.py b/molSimplify/Scripts/inparse.py index 6a58fc14..5895fbea 100644 --- a/molSimplify/Scripts/inparse.py +++ b/molSimplify/Scripts/inparse.py @@ -529,6 +529,7 @@ def parseinputfile(args, inputfile_str=None): args.dbvlinks = False args.rprompt = False set_rundir = False + args.overwrite = False # (we should remove these where posible) # THIS NEEDS CLEANING UP TO MINIMIZE DUPLICATION WITH parsecommandline @@ -930,6 +931,8 @@ def parseinputfile(args, inputfile_str=None): args.ligcon = list_to_add[0] if (l[0] == '-ligffopt'): args.ffopt = l[1] + if (l[0] == '-overwrite'): + args.overwrite = l[1] # parse postprocessing arguments if (l[0] == '-postp'): args.postp = True diff --git a/molSimplify/Scripts/structgen.py b/molSimplify/Scripts/structgen.py index 1d516e84..e5201fa2 100644 --- a/molSimplify/Scripts/structgen.py +++ b/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, @@ -2259,7 +2261,7 @@ def align_dent3_lig(args, cpoint, batoms, m3D, core3D, coreref, ligand, lig3D, return lig3D_aligned, frozenats, MLoptbds -def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] +def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int], smart_generation: bool ) -> Tuple[mol3D, List[mol3D], str, run_diag, List[int], List[int]]: """Main ligand placement routine @@ -2306,6 +2308,7 @@ def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] 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 @@ -2319,6 +2322,7 @@ def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] batslist: List[List[int]] = [] bats: List[int] = [] ffoption_list = [] # for each ligand, keeps track of what the forcefield option is. + copied = False # for determinining if core3D needs to be copied or not # load bond data MLbonds = loaddata('/Data/ML.dat') # calculate occurrences, denticities etc for all ligands @@ -2507,6 +2511,8 @@ def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] 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 +2716,39 @@ def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] 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 and copied == False: + core3D_copy = mol3D() + core3D_copy.copymol3D(core3D) + copied = True + elif copied == False: + core3D_copy = mol3D() + core3D_copy.copymol3D(core3D) + copied = True + 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 +2821,9 @@ def mcomplex(args: Namespace, ligs: List[str], ligoc: List[int] core3D.writexyz('complex_after_final_ff.xyz') # core3D,enc = ffopt(args.ff,core3D,connected,1,frozenats,freezeangles,MLoptbds,'Adaptive',args.debug) + if smart_generation == True: + core3D.bo_dict = core3D_copy.bo_dict + return core3D, complex3D, emsg, this_diag, subcatoms_ext, mligcatoms_ext @@ -2924,7 +2956,7 @@ def generate_report(args: Namespace, ligands: List[str], ligoc: List[int] 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 = False) -> Tuple[List[str], str, run_diag]: """Main structure generation routine - multiple structures Parameters @@ -2967,7 +2999,7 @@ def structgen(args: Namespace, rootdir: str, ligands: List[str], ligoc: List[int return strfiles, emsg, this_diag else: core3D, complex3D, emsg, this_diag, subcatoms_ext, mligcatoms_ext = mcomplex( - args, ligands, ligoc) + args, ligands, ligoc, smart_generation = smart_generation) if args.debug: print(('subcatoms_ext are ' + str(subcatoms_ext))) if emsg: @@ -3016,6 +3048,128 @@ def structgen(args: Namespace, rootdir: str, ligands: List[str], ligoc: List[int # 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.0 + 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