Skip to content

Commit

Permalink
Add distance plotting tutorial.
Browse files Browse the repository at this point in the history
  • Loading branch information
mgt16-LANL committed Dec 19, 2023
1 parent 24c87a1 commit d9b62ba
Show file tree
Hide file tree
Showing 5 changed files with 1,464 additions and 26 deletions.
7 changes: 4 additions & 3 deletions architector/complex_construction.py
Original file line number Diff line number Diff line change
Expand Up @@ -280,9 +280,10 @@ def final_eval(self,single_point=False):
fix_m_neighbors=True,
ff_preopt_run=True
)
self.calculator = CalcExecutor(tmp_relax.mol, parameters=self.parameters,
final_sanity_check=self.parameters['full_sanity_checks'],
relax=(not single_point)
self.calculator = CalcExecutor(tmp_relax.mol,
parameters=self.parameters,
final_sanity_check=self.parameters['full_sanity_checks'],
relax=(not single_point)
)
else: # Ensure calculation object at least exists
self.calculator = CalcExecutor(self.complexMol,method='UFF',fix_m_neighbors=False,relax=False)
Expand Down
49 changes: 42 additions & 7 deletions architector/io_molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -815,12 +815,7 @@ def append_ligand(self, ligand=None, non_coordinating=False):
newligbodict.update({newkey:val})
if lig_constraints is not None:
for ind,dist in lig_constraints.items():
if non_coordinating:
newind = natoms + ind
elif ind > 1:
newind = natoms + ind - 1
else:
newind = 1
newind = natoms + ind
self.ase_constraints.update({(0,newind):dist})
self.BO_dict.update(newligbodict)
self.ase_atoms += lig_ase_atoms
Expand Down Expand Up @@ -1295,6 +1290,7 @@ def classify_metal_geo_type(self,return_result=False):
def get_lig_dists(self,
calc_nonbonded_dists=True,
skin=0.3,
radius=None,
ref_ind='metals'):
"""Calculate metal-ligand distances and tabulate for a given structure.
Expand All @@ -1305,6 +1301,8 @@ def get_lig_dists(self,
skin : int, optional
Cutoff for nonbonded distances in Angstroms (distance considered "close" by within this many angstroms
of another coordinating atom to the metal), by default 0.3
radius : float, optional
Cutoff for radius of interaction to visualize around the molecule, by default None.
ref_ind : int, optional
Index of the atom to reference to, None will reference to all metals.
If integer passed, the distances will be calculated from this index
Expand All @@ -1321,18 +1319,24 @@ def get_lig_dists(self,
if isinstance(ref_ind,str):
if ref_ind == 'metals':
atoms = np.array(self.find_metals())
metals = atoms
else:
raise NotImplementedError('Not yet implemented for other keywords - only "metals".')
elif isinstance(ref_ind,bool):
atoms = np.array(self.find_metals())
metals = atoms
elif isinstance(ref_ind,(int,float)):
atoms = np.array([int(ref_ind)])
metals = np.array(self.find_metals())
elif isinstance(ref_ind,list):
atoms = np.array(ref_ind)
metals = np.array(self.find_metals())
elif isinstance(ref_ind,np.ndarray):
atoms = ref_ind
metals = np.array(self.find_metals())
elif ref_ind is None:
atoms = []
metals = np.array(self.find_metals())
symbols = self.ase_atoms.get_chemical_symbols()
ml_dist_dicts = []
index = 0
Expand All @@ -1347,6 +1351,7 @@ def get_lig_dists(self,
for met in atoms:
con_atoms = np.nonzero(self.graph[met])[0]
con_atom_dists = distmat[met][con_atoms]
m_visited = False
for j,c in enumerate(con_atoms):
for i,ind_set in enumerate(info_dict['original_lig_inds']):
if c in ind_set: # Find ligand this atom belongs to.
Expand All @@ -1362,13 +1367,30 @@ def get_lig_dists(self,
'atom_symbols':'{}-{}'.format(symbols[met],symbols[c])
})
index += 1
elif (c in metals) and (c != met) and (not m_visited):
m_visited=True
ml_dist_dicts.append({
'atom_pair':(met,c),
'bond_type':'explicit_bond',
'smiles':None,
'smiles_index':None,
'distance':con_atom_dists[j],
'sum_cov_radii':io_ptable.rcov1[io_ptable.elements.index(symbols[met])] + \
io_ptable.rcov1[io_ptable.elements.index(symbols[c])],
'atom_symbols':'{}-{}'.format(symbols[met],symbols[c])
})
index += 1
if calc_nonbonded_dists:
other_close_atoms = np.where(distmat[met] < (np.max(con_atom_dists)+skin))[0]
if radius is None:
other_close_atoms = np.where(distmat[met] < (np.max(con_atom_dists)+skin))[0]
else:
other_close_atoms = np.where(distmat[met] < (radius))[0]
other_close_atoms = np.array([x for x in other_close_atoms if x not in (con_atoms.tolist() + \
[int(met)])])
if len(other_close_atoms) > 0:
other_close_atom_dists = distmat[met][other_close_atoms]
for j,c in enumerate(other_close_atoms):
m_visited = False
for i,ind_set in enumerate(info_dict['original_lig_inds']):
if c in ind_set: # Find ligand this atom belongs to.
ind_in_ligand = np.where(ind_set == c)[0][0]
Expand All @@ -1383,4 +1405,17 @@ def get_lig_dists(self,
'atom_symbols':'{}-{}'.format(symbols[met],symbols[c])
})
index += 1
elif (c in metals) and (c != met) and (not m_visited):
m_visited=True
ml_dist_dicts.append({
'atom_pair':(met,c),
'bond_type':'implicit_bond',
'smiles':None,
'smiles_index':None,
'distance':other_close_atom_dists[j],
'sum_cov_radii':io_ptable.rcov1[io_ptable.elements.index(symbols[met])] + \
io_ptable.rcov1[io_ptable.elements.index(symbols[c])],
'atom_symbols':'{}-{}'.format(symbols[met],symbols[c])
})
index += 1
return pd.DataFrame(ml_dist_dicts)
3 changes: 2 additions & 1 deletion architector/io_obabel.py
Original file line number Diff line number Diff line change
Expand Up @@ -986,7 +986,8 @@ def obmol_lig_split(mol2string,
if calc_all:
for i,lig_obmol in enumerate(lig_obmols):
_,anums,_ = get_OBMol_coords_anums_graph(lig_obmol,get_types=False)
lig_coord_ats.append(','.join([io_ptable.elements[x] for x in np.array(anums)[np.array(coord_atom_lists[i])]]))
if len(coord_atom_lists[i]) > 0:
lig_coord_ats.append(','.join([io_ptable.elements[x] for x in np.array(anums)[np.array(coord_atom_lists[i])]]))
info_dict['lig_coord_ats'] = lig_coord_ats
info_dict['original_lig_inds'] = ligs_inds
info_dict['mapped_smiles_inds'] = mapped_smiles_inds
Expand Down
61 changes: 46 additions & 15 deletions architector/visualization.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,11 +58,13 @@ def type_convert(structures):
def add_bonds(view_ats,
mol,
labelsize=12,
distradius=0.3,
distvisradius=0.3,
distcolor='black',
distskin=0.3,
distopacity=0.85,
vis_distances=None,
vis_distances=None,
distradius=None,
distlabelposit=1.0,
viewer=None):
"""Add bonds to visualization displayer?
Expand All @@ -80,22 +82,38 @@ def add_bonds(view_ats,
vis_distances='metals' will do the same
vis_distances=0 will add arrows and distance labels from the atom 0 to nearby atoms.
vis_distances=[0,1] will add arrows and distances labesl from both atoms 0 and 1 to nearby atoms.
distradius : float,
distvisradius : float,
radius of drawn distance vectors, by default 0.3
distcolor : str,
color of drawn distance vectors, by default 'black'
distopacity: float,
opacity (from 0(transparent) to 1(no transparent) for drawn distance vectors, by default 0.85
distskin : float,
skin around given atom to flag "nearby" neighbors, by default 0.3
distradius : float,
Radius around a given atom to flag "nearby" neighbors, by default None.
distlabelposit : float,
Fraction of the distance (towards the ending atom) that the distance label should be placed, by default 1.0
viewer : None,
which viewer to add the arrows to, by default None
"""
if vis_distances is not None:
bondsdf = mol.get_lig_dists(calc_nonbonded_dists=True,
skin=distskin,
ref_ind=vis_distances)
ref_ind=vis_distances,
radius=distradius)
visited = list()
count = 0
for i,row in bondsdf.iterrows():
# Allow for multiple different colors of interatomic distances.
if (row['atom_pair'][0] in visited) and (isinstance(distcolor,(list,np.ndarray))):
tcolor = distcolor[visited.index(row['atom_pair'][0])]
elif (isinstance(distcolor,(list,np.ndarray))):
tcolor = distcolor[count]
visited.append(row['atom_pair'][0])
count += 1
else:
tcolor = distcolor
starting = mol.ase_atoms.get_positions()[row['atom_pair'][0]] # Should be metal.
ending = mol.ase_atoms.get_positions()[row['atom_pair'][1]]
sx= starting[0]
Expand All @@ -104,19 +122,21 @@ def add_bonds(view_ats,
ex= ending[0]
ey= ending[1]
ez= ending[2]
dxyz = np.array([ex-sx,ey-sy,ez-sz])
vector = {'start': {'x':sx, 'y':sy, 'z':sz}, 'end': {'x':ex, 'y':ey, 'z':ez},
'radius':distradius,'color':distcolor,'opacity':distopacity}
'radius':distvisradius,'color':tcolor,'opacity':distopacity}
lposit = starting + distlabelposit * dxyz
if viewer is None:
view_ats.addArrow(vector)
view_ats.addLabel("{0:.2f}".format(row['distance']), {'position':{'x':'{}'.format(ex),
'y':'{}'.format(ey),'z':'{}'.format(ez)},
view_ats.addLabel("{0:.2f}".format(row['distance']), {'position':{'x':'{}'.format(lposit[0]),
'y':'{}'.format(lposit[1]),'z':'{}'.format(lposit[2])},
'backgroundColor':"'black'",'backgroundOpacity':'0.3',
'fontOpacity':'1', 'fontSize':'{}'.format(labelsize),
'fontColor':"white",'inFront':'true'})
else:
view_ats.addArrow(vector,viewer=viewer)
view_ats.addLabel("{0:.2f}".format(row['distance']), {'position':{'x':'{}'.format(ex),
'y':'{}'.format(ey),'z':'{}'.format(ez)},
view_ats.addLabel("{0:.2f}".format(row['distance']), {'position':{'x':'{}'.format(lposit[0]),
'y':'{}'.format(lposit[1]),'z':'{}'.format(lposit[2])},
'backgroundColor':"'black'",'backgroundOpacity':'0.3',
'fontOpacity':'1', 'fontSize':'{}'.format(labelsize),
'fontColor':"white",'inFront':'true'},viewer=viewer)
Expand All @@ -126,7 +146,8 @@ def add_bonds(view_ats,
def view_structures(structures, w=200, h=200, columns=4, representation='ball_stick', labelsize=12,
labels=False, labelinds=None, vector=None, sphere_scale=0.3, stick_scale=0.25,
metal_scale=0.75, modes=None, trajectory=False, interval=200, vis_distances=None,
distradius=0.3, distcolor='black', distopacity=0.85, distskin=0.3):
distvisradius=0.3, distcolor='black', distopacity=0.85, distskin=0.3, distradius=None,
distlabelposit=1.0):
"""view_structures
Jupyter-notebook-based visualization of molecular structures.
Expand Down Expand Up @@ -190,14 +211,18 @@ def view_structures(structures, w=200, h=200, columns=4, representation='ball_st
vis_distances='metals' will do the same
vis_distances=0 will add arrows and distance labels from the atom 0 to nearby atoms.
vis_distances=[0,1] will add arrows and distances labesl from both atoms 0 and 1 to nearby atoms.
distradius : float,
distvisradius : float,
radius of drawn distance vectors, by default 0.3
distcolor : str,
color of drawn distance vectors, by default 'black'
distopacity: float,
opacity (from 0(transparent) to 1(no transparent) for drawn distance vectors, by default 0.85
distskin : float,
skin around given atom to flag "nearby" neighbors, by default 0.3
"Skin" on top of sum of cov radii around given atom to flag "nearby" neighbors, by default 0.3
distradius : float,
Radius around a given atom to flag "nearby" neighbors, by default None.
distlabelposit : float,
Fraction of the distance (towards the ending atom) that the distance label should be placed, by default 1.0
"""
mols = type_convert(structures)
if len(mols) == 1:
Expand Down Expand Up @@ -278,10 +303,12 @@ def view_structures(structures, w=200, h=200, columns=4, representation='ball_st
view_ats.addArrow(vector)
add_bonds(view_ats, mol,
labelsize=labelsize,
distradius=distradius,
distvisradius=distvisradius,
distcolor=distcolor,
distskin=distskin,
distopacity=distopacity,
distradius=distradius,
distlabelposit=distlabelposit,
vis_distances=vis_distances)
view_ats.zoomTo()
view_ats.show()
Expand Down Expand Up @@ -385,10 +412,12 @@ def view_structures(structures, w=200, h=200, columns=4, representation='ball_st
if vector:
view_ats.addArrow(vector, viewer=(x,y))
add_bonds(view_ats, mol,
distradius=distradius,
distvisradius=distvisradius,
distcolor=distcolor,
distskin=distskin,
distopacity=distopacity,
distradius=distradius,
distlabelposit=distlabelposit,
labelsize=labelsize, vis_distances=vis_distances, viewer=(x,y))
view_ats.zoomTo(viewer=(x,y))
if y+1 < columns: # Fill in columns
Expand Down Expand Up @@ -424,10 +453,12 @@ def view_structures(structures, w=200, h=200, columns=4, representation='ball_st
if vector:
view_ats.addArrow(vector)
add_bonds(view_ats, mol,
distradius=distradius,
distvisradius=distvisradius,
distcolor=distcolor,
distskin=distskin,
distopacity=distopacity,
distradius=distradius,
distlabelposit=distlabelposit,
labelsize=labelsize, vis_distances=vis_distances)
view_ats.zoomTo()
view_ats.animate({'interval':interval,'loop':'forward'}) # Infinite repetition
Expand Down
Loading

0 comments on commit d9b62ba

Please sign in to comment.