Skip to content

Commit

Permalink
Debugging get_coordination_geometry
Browse files Browse the repository at this point in the history
  • Loading branch information
@choglass
choglass committed Jan 2, 2025
1 parent b124483 commit 253c2ff
Show file tree
Hide file tree
Showing 7 changed files with 178 additions and 32 deletions.
6 changes: 6 additions & 0 deletions cell2mol/charge_assignment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -546,6 +546,7 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:
else:
combinations = [0,1]
count = 0
print(f"{combinations=}")
for com in combinations:
newlab = local_labels.copy()
newcoord = local_coords.copy()
Expand All @@ -561,6 +562,11 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:

o_s = np.sum(com)
toallocate = int(0)
print(f"{non_local_groups=}")
for jdx, a in enumerate(ligand.atoms):
if a.mconnec >= 1 and a.label not in avoid and block[jdx] == 0:
print(jdx, a.label, a.mconnec)
print("====")
for jdx, a in enumerate(ligand.atoms):
if a.mconnec >= 1 and a.label not in avoid and block[jdx] == 0:
print(a.label)
Expand Down
37 changes: 33 additions & 4 deletions cell2mol/classes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -654,6 +654,9 @@ def split_ligand(self, debug: int=0):
newgroup.get_denticity(debug=debug)
# Top-down hierarchy
self.groups.append(newgroup)
elif len(conn_idx) == 0:
if debug > 1 : print(f"\tLIGAND.SPLIT_LIGAND: no group is found")
continue
else:
if debug > 1 : print(f"\tenterting SPLIT_GROUP for the GROUP {newgroup.formula} with {conn_idx=}")
splitted_groups = split_group(newgroup, conn_idx, final_ligand_indices, debug=debug)
Expand Down Expand Up @@ -1117,6 +1120,7 @@ def get_coord_sphere(self):
def get_coord_sphere_formula(self):
if not hasattr(self,"coord_sphere"): self.get_coord_sphere()
self.coord_sphere_formula = labels2formula(list([at.label for at in self.coord_sphere]))
print(f"METAL.Get_coord_sphere_formula: {self.get_parent_index('molecule')} {self.label} {self.coord_sphere_formula}")
return self.coord_sphere_formula

#######################################################
Expand Down Expand Up @@ -1161,7 +1165,6 @@ def get_connected_groups(self, debug: int=2):
#######################################################
def get_relative_metal_radius(self, debug: int=0):
if not hasattr(self,"groups"): self.get_connected_groups(debug=debug)

diff_list = []
for group in self.groups:
if group.is_haptic == False :
Expand All @@ -1182,14 +1185,39 @@ def get_relative_metal_radius(self, debug: int=0):
self.rel_metal_radius = round(average/elemdatabase.CovalentRadius3[self.label], 3)

return self.rel_metal_radius
#######################################################
def get_connected_metals(self, debug: int=2):
self.metals = []
mol = self.get_parent("molecule")
for met in mol.metals:
if met == self : continue
tmplabels = []
tmpcoord = []
tmplabels.append(self.label)
tmpcoord.append(self.coord)
tmplabels.append(met.label)
tmpcoord.append(met.coord)

if debug > 1: print(tmplabels, tmpcoord)

isgood, tmpadjmat, tmpadjnum = get_adjmatrix(tmplabels, tmpcoord, metal_only=True)
if isgood:
if debug > 1: print(met.label, tmpadjmat, tmpadjnum)
if all(tmpadjnum[1:]):
self.metals.append(met)
else:
if debug > 1: print(f"Metal {self.label} is not connected to {met.label}")
if debug >= 2 : print(f"METAL.Get_connected_metals: {self.label} connected to {len(self.metals)} metals {[m.label for m in self.metals]}")
return self.metals

#######################################################
def get_coordination_geometry(self: object, debug: int = 0):
coord_group = self.get_connected_groups()
self.coord_nr = len(coord_group)

if debug >= 1: print(f"\nMETAL.Get_coord_geometry: {self.label}")
if debug > 2 : print(f"METAL.Get_coord_geometry:\n{coord_group=}")
if debug >= 1: print(f"METAL.Get_coord_geometry: coord_nr={self.coord_nr}")
if debug >= 2 : print(f"METAL.Get_coord_geometry:\n{coord_group=}")
if debug >= 1: print(f"METAL.Get_ccoord_geometry: coord_nr={self.coord_nr}")

self.coord_geometry, self.geom_deviation = define_coordination_geometry(self, coord_group, debug = debug)
if debug >= 2: print(f"METAL.Get_coord_geometry: {self.coord_geometry=} {self.geom_deviation=}")
Expand Down Expand Up @@ -1390,7 +1418,8 @@ def get_reference_molecules(self, ref_labels: list, ref_fracs: list, cov_factor:
ref.get_hapticity(debug=debug)
for lig in ref.ligands:
lig.get_denticity(debug=debug)
for met in ref.metals:
for met in ref.metals:
met.get_connected_metals(debug=debug)
met.get_coordination_geometry(debug=debug)
met.get_coord_sphere_formula()
else:
Expand Down
2 changes: 1 addition & 1 deletion cell2mol/connectivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ def add_atom(labels: list, coords: list, site: int, ligand: object, metalist: li
# if debug >= 2: print(f"ADD_ATOM: received {newlab=}")
# if debug >= 2: print(f"ADD_ATOM: received {tmpconmat=}")
# if debug >= 2: print(f"ADD_ATOM: received {tmpconnec=}")
# if debug >= 2: print(f"ADD_ATOM: received {tmpconnec[posadded]=}")
if debug >= 2: print(f"ADD_ATOM: received {tmpconnec[posadded]=}")
newlab_with_metal = newlab.copy()
newcoord_with_metal = newcoord.copy()
newlab_with_metal.append(tgt.label)
Expand Down
21 changes: 12 additions & 9 deletions cell2mol/coordination_sphere.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,22 +61,25 @@ def shape_measure (symbols: list, positions: list, debug: int=0) -> dict:
if debug >= 2:print(f"SHAPE_MEASURE: {positions=}")

cn = len(symbols)-1 # coordination number of metal center

connectivity= [[1, i] for i in range(2, cn+2)]
if debug >= 2: print(f"SHAPE_MEASURE: coordination number of metal center {cn}")
if debug >= 2: print(f"SHAPE_MEASURE: connectivity of metal center(1) {connectivity}")
geometry = Geometry(positions=positions,
symbols=symbols,
connectivity=connectivity)

# ref_geom = np.array(shape_structure_references['{} Vertices'.format(cn)])
ref_geom = np.array(shape_structure_references_simplified['{} Vertices'.format(cn)])

posgeom_dev={}

for idx, rg in enumerate(ref_geom[:,0]):
shp_measure = geometry.get_shape_measure(rg, central_atom=1)
geom = ref_geom[:,3][idx]
posgeom_dev[geom]=round(shp_measure, 3)
if cn == 0 :
posgeom_dev = {}
elif cn == 1 :
posgeom_dev = {'Linear' : 0.0}
else :
posgeom_dev={}
ref_geom = np.array(shape_structure_references_simplified['{} Vertices'.format(cn)])
for idx, rg in enumerate(ref_geom[:,0]):
shp_measure = geometry.get_shape_measure(rg, central_atom=1)
geom = ref_geom[:,3][idx]
posgeom_dev[geom]=round(shp_measure, 3)

return posgeom_dev

Expand Down
99 changes: 83 additions & 16 deletions cell2mol/test/check_Cell_object.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 6,
"metadata": {},
"outputs": [
{
Expand All @@ -11,7 +11,7 @@
"'2023.09.6'"
]
},
"execution_count": 2,
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -24,14 +24,16 @@
"from rdkit.Chem.Draw import IPythonConsole\n",
"from rdkit.Chem import rdDetermineBonds\n",
"import matplotlib.pyplot as plt \n",
"from cell2mol.read_write import writexyz\n",
"\n",
"IPythonConsole.ipython_3d = False\n",
"import rdkit\n",
"rdkit.__version__"
]
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -44,18 +46,18 @@
},
{
"cell_type": "code",
"execution_count": 36,
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"refcode = \"AWEVON\"\n",
"refcode = \"YAGYIP01\"\n",
"refcell = np.load(f\"{refcode}/Ref_Cell_{refcode}.cell\", allow_pickle=True)\n",
"unitcell = np.load(f\"{refcode}/Cell_{refcode}.cell\", allow_pickle=True)\n"
"# unitcell = np.load(f\"{refcode}/Cell_{refcode}.cell\", allow_pickle=True)\n"
]
},
{
"cell_type": "code",
"execution_count": 33,
"execution_count": 9,
"metadata": {},
"outputs": [
{
Expand All @@ -65,19 +67,18 @@
" Version = 2.0\n",
" Type = cell\n",
" Sub-Type = reference\n",
" Name (Refcode) = AWEVON\n",
" Num Atoms = 67\n",
" Cell Parameters a:c = [ 9.3985 9.5346 11.7627]\n",
" Cell Parameters al:ga = [78.137 86.56 70.568]\n",
" Name (Refcode) = YAGYIP01\n",
" Num Atoms = 84\n",
" Cell Parameters a:c = [12.2832 18.0705 10.4767]\n",
" Cell Parameters al:ga = [90. 90.106 90. ]\n",
"---------------------------------------------------\n",
" # of Ref Molecules: = 3\n",
" # of Ref Molecules: = 2\n",
" With Formulae: \n",
" 0: I3 \n",
" 1: H2-C-Cl2 \n",
" 2: H32-C20-Fe-As4-I2 "
" 0: H30-C18-N24-Fe-Br6 \n",
" 1: B-F4 "
]
},
"execution_count": 33,
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -86,6 +87,72 @@
"refcell"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Reference Molecule : H30-C18-N24-Fe-Br6\n",
"\t Fe N6\n",
"\t H5-C3-N4-Br\n",
"\t H5-C3-N4-Br\n",
"\t H5-C3-N4-Br\n",
"\t H5-C3-N4-Br\n",
"\t H5-C3-N4-Br\n",
"\t H5-C3-N4-Br\n",
"Reference Molecule : B-F4\n"
]
}
],
"source": [
"for i, ref in enumerate(refcell.refmoleclist):\n",
" if hasattr(ref, \"totcharge\"):\n",
" print(f\"Reference Molecule :\", ref.formula, ref.totcharge)\n",
" else:\n",
" print(f\"Reference Molecule :\", ref.formula)\n",
" if ref.iscomplex:\n",
" for met in ref.metals:\n",
" if hasattr(met, \"charge\"):\n",
" print(\"\\t\", met.label, met.charge, met.coord_sphere_formula)\n",
" else:\n",
" print(\"\\t\", met.label, met.coord_sphere_formula)\n",
" for lig in ref.ligands:\n",
" if hasattr(lig, \"totcharge\"):\n",
" print(\"\\t\", lig.formula, lig.totcharge, lig.smiles)\n",
" else:\n",
" if hasattr(lig, \"smiles\"):\n",
" print(\"\\t\", lig.formula, lig.smiles)\n",
" else:\n",
" print(\"\\t\", lig.formula)\n",
"\n",
"# print(\"Unique Species in Reference Cell:\")\n",
"# for specie in refcell.unique_species:\n",
"# if specie.subtype == \"metal\":\n",
"# if hasattr(specie, \"charge\"):\n",
"# print(\"\\t\", specie.formula, specie.charge, specie.coord_sphere_formula)\n",
"# else:\n",
"# print(\"\\t\", specie.formula, specie.coord_sphere_formula)\n",
"# else:\n",
"# if hasattr(specie, \"totcharge\"):\n",
"# print(\"\\t\", specie.formula, specie.totcharge, specie.smiles)\n",
"# else:\n",
"# print(\"\\t\", specie.formula, specie.smiles) "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"for idx, ref in enumerate(refcell.refmoleclist):\n",
" writexyz(f\"{refcode}/\", f\"{refcode}_ref_{idx}.xyz\", ref.labels, ref.coord)"
]
},
{
"cell_type": "code",
"execution_count": 37,
Expand Down
43 changes: 42 additions & 1 deletion cell2mol/test/check_radius.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,55 @@
"cells": [
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from cell2mol.elementdata import ElementData\n",
"elemdatabase = ElementData()\n",
"from cell2mol.coordination_sphere import covalent_factor_for_metal_v2"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.575"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1.3*(elemdatabase.CovalentRadius3[\"Au\"] + elemdatabase.CovalentRadius3[\"Mn\"])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6390000000000002"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1.3*(elemdatabase.CovalentRadius3[\"Fe\"] + elemdatabase.CovalentRadius3[\"N\"])"
]
},
{
"cell_type": "code",
"execution_count": 12,
Expand Down
2 changes: 1 addition & 1 deletion cell2mol/test/cif_file_BOFFOS.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2239,7 +2239,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "cell2mol",
"language": "python",
"name": "python3"
},
Expand Down

0 comments on commit 253c2ff

Please sign in to comment.