Test run 4

cell2mol/charge_assignment.py

@@ -33,7 +33,7 @@ def get_possible_charge_state(spec: object, debug: int=0):
     if not hasattr(spec,"protonation_states"): spec.get_protonation_states(debug=debug)
     if spec.protonation_states is None:                                             return None
     if spec.subtype == "group" or (spec.subtype == 'molecule' and spec.iscomplex):  return None
-
+    print(f"GET_POSSIBLE_CHARGE_STATE: {spec.formula} ({spec.subtype}) {spec.cov_factor=}")
     charge_states = []
     ### Evaluates possible charges for each protonation state ###
     for prot in spec.protonation_states:
@@ -392,7 +392,9 @@ def get_protonation_states_specie(specie: object, debug: int=0) -> list:
                         elemlist[idx] = "Cl"
                         addedlist[idx] = 1
                     else:
-                        block[idx] = 1
+                        # block[idx] = 1
+                        elemlist[idx] = "Cl"
+                        addedlist[idx] = 1
                 # Nitrogen
                 elif a.label == "N":
                     # Nitrosyl
@@ -601,15 +603,15 @@ def get_charge(charge: int, prot: object, allow: bool=True, embed_chiral: bool=T
 
     natoms = prot.natoms
     atnums = prot.atnums
-    print(f"\nGET_CHARGE. Starting get_charge with charge {charge} and {prot.formula} ({prot.added_atoms=}")
+    print(f"\nGET_CHARGE. Starting get_charge with charge {charge} and {prot.formula} {prot.added_atoms=}")
     # prot.coords and prot.cov_factor will not be used
     mols = xyz2mol(atnums, prot.coords, prot.adjmat, prot.cov_factor, charge=charge, allow_charged_fragments=allow)
-    print(f"\tGET_CHARGE.{len(mols)=} received from xyz2mol with charge {charge}")
+    print(f"GET_CHARGE.{len(mols)=} received from xyz2mol with charge {charge}")
 
     if len(mols) > 1: 
-        if debug >=1 : print(f"\tGET_CHARGE. WARNING: More than 1 mol received from xyz2mol for initcharge: {charge}")
+        if debug >=1 : print(f"GET_CHARGE. WARNING: More than 1 mol received from xyz2mol for initcharge: {charge}")
     elif len(mols) == 0:
-        if debug >=1 : print(f"\tGET_CHARGE. WARNING: No mol received from xyz2mol for initcharge: {charge}")
+        if debug >=1 : print(f"GET_CHARGE. WARNING: No mol received from xyz2mol for initcharge: {charge}")
         return None
     else :
         pass
@@ -626,7 +628,9 @@ def get_charge(charge: int, prot: object, allow: bool=True, embed_chiral: bool=T
             for mol in mols:
                 for i in range(natoms):
                     a = mol.GetAtomWithIdx(i)
-                    if debug >=2 : print(f"GET_CHARGE. {i} {a.GetSymbol()=}, {a.GetFormalCharge()=}, {a.GetImplicitValence()=}, {a.GetExplicitValence()=} {a.GetTotalValence()=}")
+                    if a.GetExplicitValence() != a.GetTotalValence():
+                        if debug >=2 : print(f"GET_CHARGE. {i} {a.GetSymbol()=}, {a.GetFormalCharge()=}, \
+                                             {a.GetImplicitValence()=}, {a.GetExplicitValence()=} {a.GetTotalValence()=}")
             return None
 
     # Smiles are generated with rdkit

cell2mol/classes.py

@@ -1524,7 +1524,7 @@ def reset_charge_assignment(self, debug: int=0):
     #######################################################
     def get_selected_cs(self, debug: int=0):
         if not hasattr(self, "unique_species"): self.get_unique_species(debug=debug)  
-        
+
         self.selected_cs = []
         for specie in self.unique_species:
             tmp = specie.get_possible_cs(debug=debug)

cell2mol/new_c2m_driver.py

@@ -9,7 +9,7 @@
 from cell2mol.other import handle_error
 from cell2mol.cell_operations import frac2cart_fromparam
 from cell2mol.new_charge_assignment import assign_charge_state_for_unique_species, balance_charge
-from cell2mol.new_cell_reconstruction import determine_wrap_keywords_pbc
+from cell2mol.new_cell_reconstruction import determine_wrap_keywords_pbc, modify_cov_factor_due_to_H, modify_cov_factor_due_to_possible_charges
 
 if __name__ == "__main__" or __name__ == "cell2mol.new_c2m_driver":
 
@@ -106,40 +106,17 @@
     refcell.get_subtype("reference")
     refcell.get_reference_molecules(ref_labels, ref_fracs, cov_factor=cov_factor, debug=debug)
 
-    if not refcell.has_isolated_H:  
-        refcell.check_missing_H(debug=debug)                                     
-    else:
-        while refcell.has_isolated_H :
-            # Increase covalent factor for H atoms
-            cov_factor += 0.05
-            refcell.get_reference_molecules(ref_labels, ref_fracs, cov_factor=cov_factor, debug=debug)
-        if debug >= 1: print(f"Covalent factor increases: {cov_factor=}")
-        refcell.check_missing_H(debug=debug)
-    refcell.assess_errors(mode="hydrogens")
+    refcell = modify_cov_factor_due_to_H(refcell, debug=debug)
     refcell.save(ref_cell_fname)
 
     if refcell.error_case == 0:
-        # Get unique species for the reference cell
-        refcell.get_unique_species(debug=debug) # Get unique_species, unique_indices, and species_list
+        refcell.get_unique_species(debug=debug) # Get unique_species, unique_indices, and species_list of the reference cell
         if debug >= 1:
             print(f"refcell.unique_species {[specie.formula for specie in refcell.unique_species]} {refcell.unique_indices=}")
             print(f"refcell.species_list {[specie.formula for specie in refcell.species_list]}\n")
-        # Get possible charge states for the unique species in the reference cell
-        refcell.get_selected_cs(debug=debug)
-        refcell.assess_errors(mode="possible_charges")
-
-        if refcell.error_get_poscharges:
-            if debug >= 1: print(f"{refcell.selected_cs=}")
-            while refcell.error_get_poscharges and cov_factor > 1.15:
-                # Decrease covalent factor for H atoms
-                cov_factor -= 0.05
-                refcell.get_reference_molecules(ref_labels, ref_fracs, cov_factor=cov_factor, debug=debug)
-                refcell.check_missing_H(debug=debug)
-                if not refcell.has_isolated_H:
-                    refcell.get_unique_species(debug=debug)
-                    refcell.get_selected_cs(debug=debug)
-
-            if debug >= 1: print(f"Covalent factor decreases: {cov_factor=}")
+
+        refcell = modify_cov_factor_due_to_possible_charges(refcell, debug=debug)
+        refcell.get_selected_cs(debug=debug) # for the last time for unique species
         refcell.assess_errors(mode="possible_charges")
     # Save reference cell object
     refcell.save(ref_cell_fname)
@@ -155,7 +132,7 @@
         reconstruction = True
         charge_assignment = False
         spin_assignment = False
-
+        cov_factor = refcell.refmoleclist[0].cov_factor
         # Get reference molecules
         newcell.get_reference_molecules(refcell.labels, refcell.frac_coord, cov_factor=cov_factor, debug=-1)
         if not newcell.has_isolated_H:  

cell2mol/new_cell_reconstruction.py

@@ -9,6 +9,54 @@
 # from cell2mol.read_write import writexyz
 # import pickle
 
+######################################################
+def modify_cov_factor_due_to_H (refcell, debug: int=0):
+    cov_factor = refcell.refmoleclist[0].cov_factor
+    if not refcell.has_isolated_H:  
+        refcell.check_missing_H(debug=debug)                                     
+    else:
+        if debug >= 1: print(f"Initial covalent factor: {cov_factor=} before increasing")
+        while refcell.has_isolated_H :
+            # Increase covalent factor for H atoms
+            cov_factor += 0.05
+            refcell.get_reference_molecules(refcell.labels, refcell.frac_coord, cov_factor=cov_factor, debug=0)
+        if debug >= 1: print(f"Covalent factor increases: {cov_factor=}")
+        refcell.check_missing_H(debug=debug)
+    refcell.assess_errors(mode="hydrogens")
+    return refcell
+
+######################################################
+def modify_cov_factor_due_to_possible_charges (refcell, debug: int=0):
+
+    cov_factor = refcell.refmoleclist[0].cov_factor
+    print(f"Initial ovalent factor: {cov_factor=}")
+
+    temp_selection = []
+    while (len(temp_selection) != len(refcell.species_list) or (refcell.has_isolated_H or refcell.has_missing_H) ) and (cov_factor > 1.15):
+        for specie in refcell.species_list:
+            tmp = specie.get_possible_cs(debug=debug)
+            if tmp is None:
+                cov_factor -= 0.05
+                refcell.get_reference_molecules(refcell.labels, refcell.frac_coord, cov_factor=cov_factor, debug=0)
+                if not refcell.has_isolated_H : refcell.check_missing_H(debug=debug)
+                if not refcell.has_missing_H :  refcell.get_unique_species(debug=debug)
+                temp_selection = []
+                break
+            elif specie.subtype != "metal":
+                temp_selection.append(list([cs.corr_total_charge for cs in specie.possible_cs]))
+            else :
+                temp_selection.append(specie.possible_cs)    
+
+    if debug >= 1: print(f"Covalent factor decreases: {cov_factor=}")
+    refcell.assess_errors(mode="hydrogens")
+    if refcell.error_case == 0:
+        if debug >= 1: print(f"OK with decreasing cov_factor {cov_factor=}")
+        return refcell
+    else:
+        if debug >= 1: print(f"Error with decreasing cov_factor {cov_factor=}")
+        refcell = modify_cov_factor_due_to_H(refcell, debug=debug)
+        return refcell
+
 ######################################################
 def apply_symmetry_operations_reference (refcell, cell_vector, sym_ops, normalize: bool=True, pbc: bool=True):
 

cell2mol/test/check_Cell_object.ipynb

@@ -372,19 +372,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [],
    "source": [
     "# name = \"ACEYOW\"\n",
-    "name = \"BILQUI\"\n",
+    "name = \"FIYFEB\"\n",
     "refcell = np.load(f\"{name}/Ref_Cell_{name}.cell\", allow_pickle=True)\n",
-    "# cell = np.load(f\"{name}/Cell_{name}.cell\", allow_pickle=True)"
+    "cell = np.load(f\"{name}/Cell_{name}.cell\", allow_pickle=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 40,
    "metadata": {
     "scrolled": true
    },
@@ -396,18 +396,18 @@
        " Version               = 2.0\n",
        " Type                  = cell\n",
        " Sub-Type              = reference\n",
-       " Name (Refcode)        = BILQUI\n",
-       " Num Atoms             = 119\n",
-       " Cell Parameters a:c   = [13.459 13.517 13.843]\n",
-       " Cell Parameters al:ga = [76.75 74.45 78.76]\n",
+       " Name (Refcode)        = FIYFEB\n",
+       " Num Atoms             = 70\n",
+       " Cell Parameters a:c   = [ 9.1546 20.9176 15.0072]\n",
+       " Cell Parameters al:ga = [90.    91.405 90.   ]\n",
        "---------------------------------------------------\n",
        " # of Ref Molecules:   = 2\n",
        " With Formulae:                                  \n",
-       "    0: H20-B-C24 \n",
-       "    1: H41-C28-P4-Fe "
+       "    0: H43-C19-N2-Si4-Fe \n",
+       "    1: K "
       ]
      },
-     "execution_count": 4,
+     "execution_count": 40,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -416,6 +416,73 @@
     "refcell"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.25\n",
+      "1.25\n",
+      "1.25\n",
+      "1.25\n"
+     ]
+    }
+   ],
+   "source": [
+    "for ref in refcell.refmoleclist:\n",
+    "    if ref.iscomplex:\n",
+    "        for lig in ref.ligands:\n",
+    "            print(lig.cov_factor)\n",
+    "    else:\n",
+    "        print(ref.cov_factor)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 50,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "H7-C7 1.25\n",
+      "H18-C6-N-Si2 1.25\n",
+      "H18-C6-N-Si2 1.25\n",
+      "K 1.25\n",
+      "H7-C7 1.25\n",
+      "H18-C6-N-Si2 1.25\n",
+      "K 1.25\n"
+     ]
+    }
+   ],
+   "source": [
+    "for specie in refcell.species_list:\n",
+    "    if specie.subtype != \"metal\":\n",
+    "        print(specie.formula, specie.cov_factor)\n",
+    "for specie in refcell.unique_species:\n",
+    "    if specie.subtype != \"metal\":\n",
+    "        print(specie.formula, specie.cov_factor)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for ref in refcell.refmoleclist:\n",
+    "    if ref.iscomplex:\n",
+    "        for lig in ref.ligands:\n",
+    "            print(lig.cov_factor)\n",
+    "    else:\n",
+    "        print(ref.cov_factor)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 5,

cell2mol/xyz2mol.py

@@ -485,7 +485,7 @@ def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
     # make a list of valences, e.g. for CO: [[4],[2,1]]
     valences_list_of_lists = []
     AC_valence = list(AC.sum(axis=1))
-    print(f"{AC_valence=}")
+    #print(f"{AC_valence=}")
     wrong = 0
 
     for i, (atomicNum, valence) in enumerate(zip(atoms, AC_valence)):
@@ -511,12 +511,12 @@ def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
             wrong += 1
             # sys.exit()
         valences_list_of_lists.append(possible_valence)
-    print(f"{wrong=}")
+    #print(f"{wrong=}")
     if wrong > 0:
         # print(f"AC2BO: {wrong=}")
         return None, atomic_valence_electrons
 
-    print(f"\tAC2BO: {valences_list_of_lists=}")
+    #print(f"\tAC2BO: {valences_list_of_lists=}")
 
     # convert [[4],[2,1]] to [[4,2],[4,1]]
     valences_list = []
@@ -553,7 +553,7 @@ def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
             check_bo = None
 
         if check_len and check_bo:
-            print(f"\tAC2BO: return AC", check_len, check_bo)
+            #print(f"\tAC2BO: return AC", check_len, check_bo)
             return AC, atomic_valence_electrons
 
         UA_pairs_list = get_UA_pairs(UA, AC, use_graph=use_graph)
@@ -581,7 +581,7 @@ def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
             )
 
             if status:
-                print(f"\tAC2BO: status", status)
+                #print(f"\tAC2BO: status", status)
                 return BO, atomic_valence_electrons
             elif (
                 BO.sum() >= best_BO.sum()
@@ -601,7 +601,7 @@ def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
             #         best_BO = BO.copy()
             #         print("AC2BO: best bo", best_BO)
         # print("best bo", best_BO)
-    print(f"\tAC2BO: return best bo")
+    #print(f"\tAC2BO: return best bo")
     #print("AC2BO: return best bo", best_BO)
     return best_BO, atomic_valence_electrons