Skip to content

Commit

Permalink
Refine Debugging mode 1
Browse files Browse the repository at this point in the history
  • Loading branch information
@choglass
choglass committed Apr 19, 2024
1 parent 55396f0 commit 70a99e9
Show file tree
Hide file tree
Showing 28 changed files with 1,393 additions and 24,545 deletions.
8 changes: 5 additions & 3 deletions cell2mol/c2m_driver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,16 +60,18 @@
################################
### PREPARES THE CELL OBJECT ###
################################
print(f"INITIATING cell object from input")

version = "2.0"
print(f"cell2mol version {version}")
print(f"INITIATING cell object from info path: {infopath}")
print(f"Debug level: {debug}")
# Reads reference molecules from info file, as well as labels and coordinates
labels, pos, ref_labels, ref_fracs, cellvec, cellparam = readinfo(infopath)
# Initiates cell
newcell = cell(name, labels, pos, cellvec, cellparam)
# Loads the reference molecules and checks_missing_H
# TODO : reconstruct the unit cell without using reference molecules
# TODO : reconstruct the unit cell using (only reconstruction of) reference molecules and Space group
newcell.get_reference_molecules(ref_labels, ref_fracs, debug=2)
newcell.get_reference_molecules(ref_labels, ref_fracs, debug=debug)
newcell.assess_errors()
newcell.save(ref_cell_fname)
######################
Expand Down
32 changes: 15 additions & 17 deletions cell2mol/cell_reconstruction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,15 +204,15 @@ def fragments_reconstruct(moleclist: list, fraglist: list, Hlist: list, refmolec

# Reconstruct Hydrogens with remaining Fragments
if len(remfrag) > 0 and len(Hlist) > 0:
print("FRAG_RECONSTRUCT.", len(fraglist), "molecules submitted to sequential with All")
for frag in fraglist:
print(frag.formula, frag.subtype, frag.labels)
print("FRAG_RECONSTRUCT.", len(fraglist), "fragments submitted to sequential with All")
if debug >= 2:
for frag in fraglist:
print("FRAG_RECONSTRUCT.", frag.formula, frag.subtype, frag.labels)

finalmols, remfrag = sequential(fraglist, refmoleclist, cellvec, factor, metal_factor, "All", debug)
moleclist.extend(finalmols)
print(f"{finalmols=}")
print(f"{remfrag=}")
for i, g in enumerate(finalmols):
if debug == 1: writexyz("/Users/ycho/cell2mol/cell2mol/test/YOBCUO/", f"reorder_molec_{i}.xyz", g.labels, g.coord)

if len(remfrag) > 0: Warning = True; print("FRAG_RECONSTRUCT. Remaining after Hydrogen reconstruction",remfrag)
else: Warning = False; print("FRAG_RECONSTRUCT. No remaining Molecules after Hydrogen reconstruction")
Expand Down Expand Up @@ -353,13 +353,13 @@ def sequential(fragmentlist: list, refmoleclist: list, cellvec: list, factor: fl
for i in range(0, len(list2)):
if i == Frag2_toallocate: sublist.append(list2[i])
elif i != Frag2_toallocate: keeplist2.append(list2[i])

print(f"sublist", len(sublist), [s.formula for s in sublist] )
print("list1", len(list1), [s.formula for s in list1])
print("list2", len(list2),[s.formula for s in list2])
print(f"keeplist1", len(keeplist1), [s.formula for s in keeplist1])
print(f"keeplist2", len(keeplist2), [s.formula for s in keeplist2])
print("")
if debug >= 2:
print(f"sublist", len(sublist), [s.formula for s in sublist] )
print("list1", len(list1), [s.formula for s in list1])
print("list2", len(list2),[s.formula for s in list2])
print(f"keeplist1", len(keeplist1), [s.formula for s in keeplist1])
print(f"keeplist2", len(keeplist2), [s.formula for s in keeplist2])
print("")
#################
# This part evaluates that the fragments that are going to be combined, can form one of the reference molecules. The resulting number of atoms is used.
#################
Expand All @@ -373,9 +373,6 @@ def sequential(fragmentlist: list, refmoleclist: list, cellvec: list, factor: fl
# Here, the function "combine" is called. It will try cell translations of one fragment, and check whether it eventually combines with the second fragment into either a bigger fragment or a molecule
#################
goodlist, avglist, badlist = combine(sublist, refmoleclist, cellvec, threshold_tmat, factor, metal_factor, debug=debug)
print(f"{goodlist=}")
print(f"{avglist=}")
print(f"{badlist=}")

#################
# This part handles the results of combine
Expand Down Expand Up @@ -516,8 +513,9 @@ def combine(tobemerged: list, references: list, cellvec: list, threshold_tmat: f
lig.get_denticity(debug=debug)
for met in reordered_newmolec.metals:
met.get_coordination_geometry(debug=debug)
print(f"{reordered_newmolec=}")
issame = compare_species(reordered_newmolec, ref, debug=2)
if debug >= 1: print(f"COMBINE: {reordered_newmolec.fomula=}")
if debug >= 1: print(f"COMBINE: {reordered_newmolec=}")
issame = compare_species(reordered_newmolec, ref, debug=debug)
if issame: ## Then is a molecule that appears in the reference list
found = True
reordered_newmolec.subtype = ref.subtype
Expand Down
51 changes: 25 additions & 26 deletions cell2mol/charge_assignment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:
for i in range(len(specie.labels)):
empty_list.append(int(0))
empty_protonation = protonation(specie.labels, specie.coord, specie.cov_factor, int(0), empty_list, empty_list, empty_list, empty_list, typ="Empty", parent=specie)
print("CREATED EMPTY PROTONATION", empty_protonation)
if debug >= 2: print(" CREATED EMPTY PROTONATION", empty_protonation)
return list([empty_protonation])

## If specie.subtype == "ligand":
Expand Down Expand Up @@ -211,8 +211,6 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:
if debug >= 2: print(" GET_PROTONATION_STATES: addressing group with hapticity:", g.haptic_type)
if debug >= 2: print(" GET_PROTONATION_STATES: and parent indices:", parent_indices)



if "h5-Cp" in g.haptic_type and not Selected_Hapticity:
Selected_Hapticity = True
tobeadded = 1
Expand Down Expand Up @@ -555,7 +553,7 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:
if debug >= 2: print(f" GET_PROTONATION_STATES: Protonation SAVED with {added_atoms} atoms added to ligand. status={new_prot.status}")
else:
if debug >= 2: print(f" GET_PROTONATION_STATES: Protonation DISCARDED. Steric Clashes found when adding atoms. status={new_prot.status}")
print(f"{protonation_states=}")
if debug >= 2: print(f" GET_PROTONATION_STATES:{protonation_states=}")
return protonation_states

#######################################################
Expand All @@ -567,7 +565,6 @@ def get_charge(ich: int, prot: object, allow: bool=True, debug: int=0):

natoms = prot.natoms
atnums = prot.atnums
if debug >= 2: print(f"*****get_charge******")

##########################
# xyz2mol is called here #
Expand All @@ -583,7 +580,7 @@ def get_charge(ich: int, prot: object, allow: bool=True, debug: int=0):

# Smiles are generated with rdkit
smiles = Chem.MolToSmiles(mols[0])
print(f"{smiles=}")
if debug >= 2: print(f"GET_CHARGE. {smiles=}")
# Gets the resulting charges
atom_charge = []
total_charge = 0
Expand Down Expand Up @@ -838,7 +835,7 @@ def balance_charge(unique_indices: list, unique_species: list, debug: int=0) ->

final_charge_distribution = []
for idx, d in enumerate(alldistr):
print(f"{d=}")
if debug >= 2: print(f"BALANCE: distribution={d}")
final_charge = np.sum(d)
if final_charge == 0:
final_charge_distribution.append(d)
Expand Down Expand Up @@ -875,12 +872,12 @@ def prepare_unresolved(unique_indices: list, unique_species: list, distributions
return list_molecules, list_indices, list_options

#######################################################
def set_charges_create_bonds (specie, unique_indices, unique_species, final_charge_distribution, debug):
def set_target_charge (specie, unique_indices, unique_species, final_charge_distribution, debug):

spec = unique_species[specie.unique_index]
indices = [index for index, value in enumerate(unique_indices) if value == specie.unique_index]
target_charge = [final_charge_distribution[i] for i in indices][0]
if debug > 1: print(spec, indices, target_charge)
if debug > 1: print("SET_TARGET_CHARGE:", spec, indices, target_charge)

if (specie.subtype == "molecule" and specie.iscomplex == False) or (specie.subtype == "ligand"):
formula = specie.formula
Expand All @@ -891,57 +888,55 @@ def set_charges_create_bonds (specie, unique_indices, unique_species, final_char
charge_list = spec.possible_cs

if target_charge in charge_list:
if debug > 1: print(f"Target charge {target_charge} of {formula} exists in {charge_list}." )
if debug > 1: print(f"SET_TARGET_CHARGE: Target charge {target_charge} of {formula} exists in {charge_list}." )
else:
if debug > 1: print(f"ERROR: Target charge {target_charge} of {formula} does not exist in {charge_list}." )
if debug >= 1: print(f"SET_TARGET_CHARGE: ERROR!! Target charge {target_charge} of {formula} does not exist in {charge_list}." )
return None

if (specie.subtype == "molecule" and specie.iscomplex == False) or (specie.subtype == "ligand"):
print(specie.formula)
if debug > 1: print("SET_TARGET_CHARGE:", specie.formula)
specie.get_protonation_states(debug=debug)
specie.get_possible_cs(debug=debug)
formula = specie.formula
charge_list = [cs.corr_total_charge for cs in specie.possible_cs]

if target_charge in charge_list:
if debug > 1: print(f"Target charge {target_charge} of {formula} exists in {charge_list}.")
if debug > 1: print(f"SET_TARGET_CHARGE: Target charge {target_charge} of {formula} exists in {charge_list}.")
idx = charge_list.index(target_charge)
cs = specie.possible_cs[idx]
#prot = cs.protonation
specie.set_charges(cs.corr_total_charge, cs.corr_atom_charges, cs.smiles, cs.rdkit_obj)

else:
if debug > 1: print(f"ERROR: Target charge {target_charge} of {formula} does not exist in {charge_list}." )
if debug >= 1: print(f"SET_TARGET_CHARGE: ERROR!! Target charge {target_charge} of {formula} does not exist in {charge_list}." )
return None

elif specie.subtype == "metal":
print(specie.label)
if debug > 1: print("SET_TARGET_CHARGE:", specie.label)
specie.get_possible_cs(debug=debug)
formula = specie.label
charge_list = spec.possible_cs

if target_charge in charge_list:
if debug > 1: print(f"Target charge {target_charge} of {formula} exists in {charge_list}." )
if debug > 1: print(f"SET_TARGET_CHARGE: Target charge {target_charge} of {formula} exists in {charge_list}." )
idx = charge_list.index(target_charge)
cs = specie.possible_cs[idx]
specie.set_charge(cs)
else:
if debug > 1: print(f"ERROR: Target charge {target_charge} of {formula} does not exist in {charge_list}." )
if debug >= 1: print(f"SET_TARGET_CHARGE: ERROR!! Target charge {target_charge} of {formula} does not exist in {charge_list}." )
return None
#######################################################
def prepare_mols_v4 (moleclist: list, unique_indices: list, unique_species: list, final_charge_distribution: list, debug: int=0):
count = 0
for mol in moleclist:
if mol.iscomplex == False:
set_charges_create_bonds(mol, unique_indices, unique_species, final_charge_distribution, debug)
set_target_charge(mol, unique_indices, unique_species, final_charge_distribution, debug)
count += 1

elif mol.iscomplex:
tmp_atcharge = np.zeros((mol.natoms))
tmp_smiles = []

for lig in mol.ligands:
set_charges_create_bonds(lig, unique_indices, unique_species, final_charge_distribution, debug)
set_target_charge(lig, unique_indices, unique_species, final_charge_distribution, debug)
count += 1

tmp_smiles.append(lig.smiles)
Expand All @@ -950,7 +945,7 @@ def prepare_mols_v4 (moleclist: list, unique_indices: list, unique_species: list
tmp_atcharge[a] = lig.atomic_charges[kdx]

for met in mol.metals:
set_charges_create_bonds(met, unique_indices, unique_species, final_charge_distribution, debug)
set_target_charge(met, unique_indices, unique_species, final_charge_distribution, debug)
count += 1
parent_index = met.get_parent_index("molecule")
tmp_atcharge[parent_index] = met.charge
Expand Down Expand Up @@ -1119,18 +1114,17 @@ def correct_smiles_ligand(ligand: object, debug: int=0) -> Tuple[str, object]:

# Sets bond information and hybridization
for jdx, atom in enumerate(ligand.atoms):
if debug >=2: print(jdx, atom.label)
nbonds = 0
for b in atom.bonds:
if debug >=1: print(b.atom1.label, b.atom2.label, b.order)
# if debug >=2 : print(b.atom1.label, b.atom2.label, b.order)
ismetal_1 = elemdatabase.elementblock[b.atom1.label] == "d" or elemdatabase.elementblock[b.atom1.label] == "f"
ismetal_2 = elemdatabase.elementblock[b.atom2.label] == "d" or elemdatabase.elementblock[b.atom2.label] == "f"
if ismetal_1 or ismetal_2:
pass
else:
begin_idx = b.atom1.get_parent_index("ligand")
end_idx = b.atom2.get_parent_index("ligand")
if debug >=2: print(begin_idx, end_idx)
# if debug >=2: print(begin_idx, end_idx)
nbonds += 1
if b.order== 1.0:
btype = Chem.BondType.SINGLE
Expand Down Expand Up @@ -1163,6 +1157,11 @@ def correct_smiles_ligand(ligand: object, debug: int=0) -> Tuple[str, object]:
## visulize a corrected rdkit object
if debug >=2:
from IPython.display import display
from rdkit.Chem.Draw import IPythonConsole
IPythonConsole.drawOptions.addAtomIndices = True
IPythonConsole.molSize = 300,300

print(f"{ligand.formula=} {smiles=}")
display(mol_with_atom_index(obj))

return smiles, obj
Expand Down Expand Up @@ -1229,7 +1228,7 @@ def get_smiles_complex (mol: object, debug: int=0) -> Tuple[str, object]:
IPythonConsole.drawOptions.addAtomIndices = True
IPythonConsole.molSize = 300,300

print(f"{smiles=}")
print(f"{mol.formula=} {smiles=}")
display(obj)

return smiles, obj
Expand Down
Loading

0 comments on commit 70a99e9

Please sign in to comment.