Skip to content

Commit

Permalink
finished prepare_mols, small fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@svela-bs
svela-bs committed Feb 27, 2024
1 parent b46921d commit 127f67a
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 84 deletions.
3 changes: 2 additions & 1 deletion cell2mol/c2m_module.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
87 changes: 21 additions & 66 deletions cell2mol/charge_assignment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}")

Expand All @@ -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
Expand All @@ -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 = []
Expand All @@ -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

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

Expand Down
82 changes: 65 additions & 17 deletions cell2mol/classes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down Expand Up @@ -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
Expand All @@ -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

#######################################################
Expand Down Expand Up @@ -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 ######
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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")

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

0 comments on commit 127f67a

Please sign in to comment.