diff --git a/cell2mol/c2m_module.py b/cell2mol/c2m_module.py index 5ff3be27..9224a65a 100644 --- a/cell2mol/c2m_module.py +++ b/cell2mol/c2m_module.py @@ -31,7 +31,8 @@ def cell2mol(newcell: object, reconstruction: bool=True, charge_assignment: bool tini = time.time() if not newcell.is_fragmented: newcell.reset_charges() - newcell.determine_charge(debug=debug) + newcell.assign_charges(debug=debug) + newcell.create_bonds(debug=debug) tend = time.time() if debug >= 1: print(f"\nTotal execution time for Charge Assignment: {tend - tini:.2f} seconds") diff --git a/cell2mol/charge_assignment.py b/cell2mol/charge_assignment.py index 89067c3e..09d2366a 100644 --- a/cell2mol/charge_assignment.py +++ b/cell2mol/charge_assignment.py @@ -897,15 +897,17 @@ def prepare_mols(moleclist: list, unique_indices: list, unique_species: list, se else: if debug >= 2: print(f"PREPARE: Error doing molecule {idx}. Created Charge State is different than Target {new_cs.corr_total_charge} vs {cs.corr_total_charge}") return None + if not allocated: Warning = False ########################### ###### FOR LIGANDS ###### ########################### - elif mol.iscomplex: + elif mol.subtype == "molecule" and mol.iscomplex: if debug >= 2: print(f"PREPARE: Molecule {moleclist.index(mol)} has {len(mol.ligands)} ligands") for kdx, lig in enumerate(mol.ligands): spec = unique_species[lig.unique_index] # This is the reference specie + if debug >= 2: print("") if debug >= 2: print(f"PREPARE: Doing Ligand {kdx} with unique_index: {lig.unique_index}") if debug >= 2: print(f"PREPARE: Specie with poscharges: {spec.possible_cs}") @@ -922,19 +924,7 @@ def prepare_mols(moleclist: list, unique_indices: list, unique_species: list, se ref_data.append(a.label+str(a.mconnec)) dummy1, dummy2, map12 = reorder(ref_data, target_data, spec.coord, lig.coord) ############### - - ## Recomputes the protonation states, but I think it should be able to use the one in cs prot = cs.protonation.reorder(map12) - #list_of_protonations = get_protonation_states(lig, debug=debug) - #found_prot = False - #for p in list_of_protonations: - # if p.os == spec.protonation.os and p.added_atoms == spec.protonation.added_atoms and not found_prot: - # tmp_addedlist = list(np.array(p.addedlist)[map12]) - # if debug >= 2: print(f"PREPARE: tmp_addedlist={tmp_addedlist}") - # if all(tmp_addedlist[idx] == spec.protonation.addedlist[idx] for idx in range(len(p.addedlist))): - # if debug >= 2: print(f"PREPARE: found match in protonation with tmpsmiles:{p.tmpsmiles}") - # prot = p - # found_prot = True if lig.is_nitrosyl: if lig.NO_type == "Linear": NOcharge = 1 #NOcharge is the charge with which I need to run getcharge to make it work @@ -952,64 +942,27 @@ def prepare_mols(moleclist: list, unique_indices: list, unique_species: list, se else: lig.set_charges(new_cs.corr_total_charge, new_cs.corr_atom_charges, new_cs.smiles, new_cs.rdkit_mol) if debug >= 1: print(f"PREPARE: Success doing ligand {kdx}. Created Charge State with total_charge={new_cs.corr_total_charge}") + allocated = True - ### SERGI: DONE UNTIL HERE + if allocated: idxtoallocate += 1 + else: Warning = False - else: - if debug >= 1: print(f"PREPARE: WARNING, I Could not identify the protonation state. I'll try to obtain the desired result") - found_charge_state = False - for prot in list_of_protonations: - list_of_charge_states = [] - list_of_protonations_for_each_state = [] - - tmpobject = ["Ligand", lig, mol] - chargestried = get_list_of_charges_to_try(spec, prot) - for ich in chargestried: - ch_state = get_charge(prot.labels, prot.coords, prot.conmat, ich, prot.cov_factor) - list_of_charge_states.append(ch_state) - list_of_protonations_for_each_state.append(prot) - if debug >= 1: print(f" POSCHARGE: charge 0 with smiles {ch_state.smiles}") - - if len(list_of_charge_states) > 0: - best_charge_distr_idx = select_charge_distr(list_of_charge_states, debug=debug) - else: - if debug >= 1: print(f" POSCHARGE. found EMPTY best_charge_distr_idx for PROTONATION state") - best_charge_distr_idx = [] - - if debug >= 2: print(f" POSCHARGE. best_charge_distr_idx={best_charge_distr_idx}") - for idx in best_charge_distr_idx: - c = list_of_charge_states[idx] - p = list_of_protonations_for_each_state[idx] - if debug >= 2: print(f" POSCHARGE. {c.corr_total_charge}={ch}, {p.added_atoms}={tgt_protonation.added_atoms}") - if c.corr_total_charge == ch and p.added_atoms == tgt_protonation.added_atoms: - lig.charge(c.corr_atom_charges, c.corr_total_charge, c.rdkit_mol, c.smiles) - if debug >= 1: print(f"PREPARE: Success doing ligand {kdx}. Created Charge State with total_charge={c.corr_total_charge}") - found_charge_state = True - - if not found_charge_state: Warning = True - if allocated: - idxtoallocate += 1 - else: - idxtoallocate += 1 - if debug >= 1: print(f"PREPARE: Warning allocating molecule {idx} with {final_charge_distribution[idxtoallocate]} as target charge") - Warning = True - for kdx, met in enumerate(mol.metals): - specie = unique_indices[idxtoallocate] - spec_object = unique_species[specie][1] - allocated = False + spec = unique_species[met.unique_index] if debug >= 2: print("") if debug >= 2: print(f"PREPARE: Metal {kdx}, label {met.label} is specie {specie}") - if debug >= 2: print(f"PREPARE: Metal possible_css: {spec_object.possible_cs}") - for jdx, ch in enumerate(spec_object.possible_cs): - if final_charge_distribution[idxtoallocate] == ch and not allocated: + if debug >= 2: print(f"PREPARE: Metal possible_css: {spec.possible_cs}") + allocated = False + for jdx, cs in enumerate(spec.possible_cs): + if final_charge_distribution[idxtoallocate] == cs.corr_total_charge and not allocated: allocated = True - met.charge(ch) - if allocated: - idxtoallocate += 1 + met.set_charge(cs.corr_total_charge) + + if allocated: idxtoallocate += 1 + else: Warning = False if not Warning: - # Now builds the Charge Data for the final molecule. Smiles is a list with all ligand smiles separately. + # Now sets the charge for the final molecule, using metal and ligand data if debug >= 2: print(f"PREPARE: Building Molecule {idx} From Ligand&Metal Information") tmp_atcharge = np.zeros((mol.natoms)) tmp_smiles = [] @@ -1018,13 +971,13 @@ def prepare_mols(moleclist: list, unique_indices: list, unique_species: list, se for kdx, a in enumerate(lig.parent_indices): tmp_atcharge[a] = lig.atcharge[kdx] for met in mol.metals: - tmp_atcharge[met.atlist] = met.totcharge - - mol.charge(tmp_atcharge, int(sum(tmp_atcharge)), [], tmp_smiles) + tmp_atcharge[met.parent_index] = met.charge + mol.set_charges(int(sum(tmp_atcharge)), atomic_charges=tmp_atcharge, smiles=tmp_smiles) return moleclist, Warning ####################################################### +<<<<<<< HEAD def build_bonds(moleclist: list, debug: int=0) -> list: from cell2mol.classes import bond @@ -1172,6 +1125,8 @@ def build_bonds(moleclist: list, debug: int=0) -> list: return moleclist ####################################################### +======= +>>>>>>> 3c676b24 (finished prepare_mols, small fixes) def correct_smiles_ligand(ligand: object): ## Receives a ligand class object and constructs the smiles and the rdkit_mol object from scratch, using atoms and bond information diff --git a/cell2mol/classes.py b/cell2mol/classes.py index 0124b5e1..fdfd92b5 100644 --- a/cell2mol/classes.py +++ b/cell2mol/classes.py @@ -212,7 +212,7 @@ def create_bonds(self, debug: int=0): #if debug >= 2: print("BUILD BONDS: atom", idx, a.label) rdkitatom = self.rdkit_mol.GetAtomWithIdx(idx) tmp = rdkitatom.GetSymbol() - if atom.label != tmp: print("Error in Build_Bonds. Atom labels do not coincide. GMOL vs. MOL:", atom.label, tmp) + if atom.label != tmp: print("Error in Create Bonds. Atom labels do not coincide. GMOL vs. MOL:", atom.label, tmp) else: # First part. Creates bond information for b in rdkitatom.GetBonds(): @@ -369,6 +369,13 @@ def __init__(self, labels: list, coord: list, parent_indices: list=None, radii: # remove = group.check_coordination(debug=debug) # self.reset_adjacencies_lig_metalist_v2 (remove, debug=debug) + ####################################################### + def correction_smiles(self): + if not hasattr(self.parent,"smiles"): self.parent.smiles = [] + self.smiles, self.rdkit_mol = correct_smiles_ligand(self) + self.parent.smiles.append(self.smiles) + + ####################################################### def get_connected_metals(self, debug: int=0): # metal.groups will be used for the calculation of the relative metal radius # and define the coordination geometry of the metal /hapicitiy/ hapttype @@ -379,7 +386,7 @@ def get_connected_metals(self, debug: int=0): tmplabels.append(met.label) tmpcoord.append(met.coord) isgood, tmpadjmat, tmpadjnum = get_adjmatrix(tmplabels, tmpcoord, metal_only=True) - if isgood and any(tmpadjnum) > 0: self.groups.append(group) + if isgood and any(tmpadjnum) > 0: self.metals.append(met) return self.metals ####################################################### @@ -571,9 +578,23 @@ def get_denticity(self, debug: int=0): ############### class bond(object): def __init__(self, atom1: object, atom2: object, bond_order: int=1): - self.atom1 = atom1 - self.atom2 = atom2 - self.order = bond_order + self.type = "bond" + self.version = "0.1" + self.atom1 = atom1 + self.atom2 = atom2 + self.order = bond_order + + def __repr__(self): + to_print = f'---------------------------------------------------\n' + to_print += ' >>> Cell2mol BOND Object >>> \n' + to_print += f'---------------------------------------------------\n' + to_print += f' Version = {self.version}\n' + to_print += f' Type = {self.type}\n' + to_print += f' Atom 1 = {self.atom1.parent_index}\n' + to_print += f' Atom 2 = {self.atom2.parent_index}\n' + to_print += f' Bond Order = {self.order}\n' + to_print += '----------------------------------------------------\n' + return to_print ############### ### ATOM ###### @@ -605,9 +626,18 @@ def check_connectivity(self, other: object, debug: int=0): else: return False ####################################################### - def add_bond(self,newbond: object): + def add_bond(self, newbond: object, debug: int=0): if not hasattr(self,"bonds"): self.bonds = [] - self.bonds.append(newbond) + at1 = newbond.atom1 + at2 = newbond.atom2 + found = False + for b in self.bonds: + if (b.atom1 == at1 and b.atom2 == at2) or (b.atom1 == at2 and b.atom2 == at1): + if debug > 0: print(f"ATOM.ADD_BOND found the same bond with atoms:") + if debug > 0: print(f"atom1: {b.atom1}") + if debug > 0: print(f"atom2: {b.atom2}") + found = True  ### It means that the same bond has already been defined + if not found: self.bonds.append(newbond) ####################################################### def set_adjacency_parameters(self, cov_factor: float, metal_factor: float) -> None: @@ -634,6 +664,19 @@ def set_adjacencies(self, adjmat, madjmat, connectivity: int, metal_connectivity for idx, c in enumerate(madjmat): ## The atom only receives one row of madjmat, so this is not a matrix anymore if c > 0: self.metal_adjacency.append(idx) + + ####################################################### + def get_connected_metals(self, metalist: list, debug: int=0): + self.metals = [] + for met in metalist: + tmplabels = self.label.copy() + tmpcoord = self.coord.copy() + tmplabels.append(met.label) + tmpcoord.append(met.coord) + isgood, tmpadjmat, tmpadjnum = get_adjmatrix(tmplabels, tmpcoord, metal_only=True) + if isgood and any(tmpadjnum) > 0: self.metals.append(met) + return self.metals + ####################################################### def get_closest_metal(self, metalist: list, debug: int=0): ## Here, the list of metal atoms must be provided @@ -754,7 +797,7 @@ def get_spin(self): ############ def reset_charge(self): - atom.reset_charge(self) ## First uses the generic specie class function for itself and its atoms + atom.reset_charge(self) ## First uses the generic atom class function for itself if hasattr(self,"poscharges"): delattr(self,"poscharge") ############## @@ -1053,23 +1096,28 @@ def create_bonds(self, debug: int=0): # isconnected = check_connectivity(at, met) if isconnected: newbond = bond(at, met, 0.5) + at.add_bond(newbond) + met.add_bond(newbond) count += 1 if count != at.mconnec: print(f"CELL.CREATE_BONDS: error creating bonds for atom: \n{atom}\n of ligand: \n{lig}\n") print(f"CELL.CREATE_BONDS: count differs from atom.mconnec: {count}, {at.mconnec}") - # Adds Metal-Metal Bonds, with an arbitrary 0.5 order: + # Adds Metal-Metal Bonds, with an arbitrary 0.5 order: + if mol.iscomplex: + for idx, met1 in mol.metals: + for jdx, met2 in mol.metals: + if idx <= jdx: continue + isconnected = check_connectivity(met1, met2) + if isconnected: + newbond = bond(met1, met2, 0.5) + met1.add_bond(newbond) + met2.add_bond(newbond) + # Fourth part : correction smiles of ligands mol.smiles_with_H = [lig.smiles for lig in mol.ligands] - for lig in mol.ligands: - lig.correction_smiles() - - ####################################################### - def correction_smiles (self): - if not hasattr(self.parent,"smiles"): self.parent.smiles = [] - self.smiles, self.rdkit_mol = correct_smiles_ligand(self) - self.parent.smiles.append(self.smiles) + mol.smiles = [lig.correction_smiles() for lig in mol.ligands] ####################################################### def assign_spin(self, debug: int=0) -> object: