Skip to content

Commit

Permalink
Merge pull request #296 from rolan701/main
Browse files Browse the repository at this point in the history 
structure generation overhaul
  • Loading branch information
@rolan701
rolan701 authored Nov 21, 2024
2 parents 7765d71 + 7a86b90 commit 1c8fb7e
Show file tree
Hide file tree
Showing 5 changed files with 257 additions and 5 deletions.
94 changes: 94 additions & 0 deletions molSimplify/Classes/mol3D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
2 changes: 1 addition & 1 deletion molSimplify/Scripts/generator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
3 changes: 3 additions & 0 deletions molSimplify/Scripts/inparse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
160 changes: 157 additions & 3 deletions molSimplify/Scripts/structgen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -25,6 +26,7 @@
getPointu,
kabsch,
midpt,
move_point,
norm,
PointTranslateSph,
reflect_through_plane,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion tests/test_inparse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@ def test_parseinputfile_empty():
defaults = {'skipANN': False, 'oldANN': False,
'dbvdent': False, 'dbvconns': False,
'dbvhyb': False, 'dbvlinks': False,
'rprompt': False, 'rundir': f'{os.getcwd()}/Runs'}
'rprompt': False, 'rundir': f'{os.getcwd()}/Runs',
'overwrite': False,}

args = Namespace()
parseinputfile(args, inputfile_str=' ')
Expand Down

0 comments on commit 1c8fb7e

Please sign in to comment.