In process

uibcdf · Apr 30, 2024 · 3a7d7dd · 3a7d7dd
1 parent bccc30a
commit 3a7d7dd
Show file tree

Hide file tree

Showing 3 changed files with 425 additions and 32 deletions.
diff --git a/molsysmt/build/get_missing_bonds.py b/molsysmt/build/get_missing_bonds.py
@@ -40,6 +40,8 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
             h_atoms_not_water = [ii for ii in atoms_not_water if ii not in heavy_atoms_not_water] 
             h_atoms_water = [ii for ii in atoms_water if ii not in o_atoms_water] 
 
+
+            # bonds in water molecules are defined by template and independently
             if len(atoms_water):
                 if len(h_atoms_water):
                     print('System with hydrogens and waters... to be implemented')
@@ -72,8 +74,20 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
                 bonds_distances += contacts_heavy_atoms[0]
 
                 if len(h_atoms_not_water):
-                    print('System with hydrogens... to be implemented')
-                    raise NotImplementedError()
+
+                    aux_contacts = get_contacts(molecular_system, selection=heavy_atoms_not_water,
+                                                selection_2=h_atoms_not_water, threshold=threshold,
+                                                output_type='pairs', output_indices='atom',
+                                                pbc=True, skip_digestion=True)
+
+                    # Hydrogen bonds need to be removed: let's work with those bonds intra group only.
+                    contacts_heavy_atoms_with_h_atoms = []
+                    group_indices = get(molecular_system, element='atom', group_index=True)
+                    for ii,jj in aux_contacts[0]:
+                        if group_indices[ii]==group_indices[jj]:
+                            contacts_heavy_atoms_with_h_atoms.append([ii,jj])
+
+                    bonds_distances += contacts_heavy_atoms_with_h_atoms
 
             if with_templates:
 
@@ -202,11 +216,6 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
 
                 raise NotImplementedError
 
-            #if sorted:
-            #    output = _sorted(bonds)
-            #else:
-            #    output = bonds
-
             if with_distances and with_templates:
 
                 distances_in_templates = True
@@ -223,8 +232,11 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
                         break
 
                 if distances_in_templates and templates_in_distances:
-                   output = bonds_distances 
-
+                   output = bonds_distances
+                else:
+                    #raise ValueError('templates and distances do not match')
+                    print('templates and distances do not match')
+                    return bonds_distances, bonds_templates
             else:
 
                 raise NotImplementedError

diff --git a/molsysmt/structure/get_contacts.py b/molsysmt/structure/get_contacts.py
@@ -57,10 +57,17 @@ def get_contacts(molecular_system, selection=None, center_of_atoms=False, weight
         pairs = []
         n_contact_maps = contact_map.shape[0]
 
-        for ii in range(n_contact_maps):
-            aux_pairs = np.nonzero(np.triu(contact_map[ii],k=1)==True)
-            aux_pairs = np.column_stack(aux_pairs).tolist()
-            pairs.append(aux_pairs)
+
+        if selection_2 is None:
+            for ii in range(n_contact_maps):
+                aux_pairs = np.nonzero(np.triu(contact_map[ii],k=1)==True)
+                aux_pairs = np.column_stack(aux_pairs).tolist()
+                pairs.append(aux_pairs)
+        else:
+            for ii in range(n_contact_maps):
+                aux_pairs = np.nonzero(contact_map[ii])
+                aux_pairs = np.column_stack(aux_pairs).tolist()
+                pairs.append(aux_pairs)
 
         if output_indices=='atom':
 
@@ -73,8 +80,8 @@ def get_contacts(molecular_system, selection=None, center_of_atoms=False, weight
                 atom_indices_2 = np.array(atom_indices_2)
                 for ii in range(n_contact_maps):
                     aux_pairs = np.array(pairs[ii])
-                    pairs[ii]=np.column_stack(atom_indices[aux_pairs[:,0]],
-                            atom_indices_2[aux_pairs[:,1]]).tolist()
+                    pairs[ii]=np.column_stack([atom_indices[aux_pairs[:,0]],
+                            atom_indices_2[aux_pairs[:,1]]]).tolist()
 
         if output_type=='sorted pairs':
             for ii in range(n_contact_maps):