Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

smart generation added to structgen.py #290

Closed
wants to merge 5 commits into from
Closed

smart generation added to structgen.py #290

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e53af75
structgen upload... scary
rolan701 Nov 15, 2024
dc29fa7
Update structgen.py
rolan701 Nov 15, 2024
249141f
Update structgen.py
rolan701 Nov 15, 2024
efaa28c
[pre-commit.ci] auto fixes from pre-commit hooks
pre-commit-ci[bot] Nov 15, 2024
d172d30
Update structgen.py
rolan701 Nov 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 148 additions & 1 deletion molSimplify/Scripts/structgen.py
View file
Open in desktop
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 @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Check warning

Code scanning / CodeQL

Variable defined multiple times Warning

This assignment to 'dist' is unnecessary as it is
redefined
Loading
before this value is used.
This assignment to 'dist' is unnecessary as it is
redefined
Loading
before this value is used.
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
Loading