Testing issue with many atom for which big frag multipole is necessary

PoltypeModules/lFragmenterForDMA.py

@@ -55,7 +55,7 @@ def saveFragment(mol, atomsToKeep, fragname, bondsToDelete=[]):
 
 def growFragment(atomidx, sdffile):
   mol = Chem.MolFromMolFile(sdffile,removeHs=False)
-
+  print(mol) 
   # 1. Atom in three rings
   # 2. Atom in two rings
   # 3. Atom in one ring
@@ -103,6 +103,16 @@ def growFragment(atomidx, sdffile):
     if len(match) == 1:
       atomsInSixMemRing.append(match[0])
 
+  atomsInTripleBond = []
+  pattern = Chem.MolFromSmarts('[$(*#*)]')
+  matches = mol.GetSubstructMatches(pattern)
+  for match in matches:
+    if len(match) == 1:
+      atomsInTripleBond.append(match[0])
+
+  print('Atoms in triple bond',atomsInTripleBond)
+
+
   # Case 1: the atom is in three rings
   if (atomidx - 1) in atomsInThreeRing:
     atomsToKeep = []
@@ -207,9 +217,165 @@ def growFragment(atomidx, sdffile):
         shutil.move(f"Frag_Atom{atomidx:03d}_0.mol", f"Frag_Atom{atomidx:03d}.mol")
       else:
         os.system(f"rm -f Frag_Atom{atomidx:03d}_0.mol")
-
+
+  # Case when atom is part of triple bond
+  if (atomidx -1) in atomsInTripleBond:
+    atomsToAdd = []
+    atomsToKeep = [atomidx-1]
+    t1 = mol.GetAtomWithIdx(atomidx-1)
+    for neig in t1.GetNeighbors():
+      neig_idx = neig.GetIdx()
+      atomsToAdd.append(neig_idx)
+      atomsToKeep.append(neig_idx)
+    print('Atoms to add ',atomsToAdd)
+    print('Atoms in ring',atomsInRing) 
+    # Loop through the atoms to add and check:
+    # if the atom is in an aromatic ring
+    # if not, then try to find coonected atoms with specific rules
+    for neig in atomsToAdd:
+
+      if neig in atomsInRing:
+        print('Found atom in ring as first neighboor')
+        connectedAtoms = findConnectedAtoms(mol, atomsToAdd, atomsInRing)
+        atomsToKeep += connectedAtoms
+
+      elif neig in atomsInTripleBond:
+        print(f'Neig {neig} is in atomsinTripleBond')
+        atomsToKeep += [neig]
+      else:
+        print(f'Neig {neig} is neither in atomsinTripleBOnd nor in atomsinRing')
+        connectedAtoms = findConnectedAtomsGeneral(mol, atomsToAdd, atomidx, neig)
+        atomsToKeep += connectedAtoms
+
+    print('Final atoms to keep',atomsToKeep)
+    saveFragment(mol, atomsToKeep, f"Frag_Atom{atomidx:03d}.mol")
+    if os.path.isfile(f"Frag_Atom{atomidx:03d}_0.mol"):
+      testmol = Chem.MolFromMolFile(f"Frag_Atom{atomidx:03d}.mol",removeHs=False)
+      if len(testmol.GetAtoms()) == len(mol.GetAtoms()):
+        shutil.move(f"Frag_Atom{atomidx:03d}_0.mol", f"Frag_Atom{atomidx:03d}.mol")
+      else:
+        os.system(f"rm -f Frag_Atom{atomidx:03d}_0.mol")
+
   return 
 
+
+def find_patern(mol,atm_to_check,pat_to_check):
+
+  """
+  This function find the atoms present in a specific pattern
+  and check if a given atom is part of the atom present in the 
+  pattern.
+
+  Inputs:
+    -   mol: rdkit mol object
+    -   atm_to_check: rdkit atom id to check
+    -   pat_to_check: SMART string to check as pattern
+
+  Outputs:
+    -   is_present: if atom_to_check exist or not in the SMART pattern (bool)
+    -   Search_pattern: list of rdkit atom ID present in the pattern
+  """
+
+
+  Search_pattern = []
+  pattern = Chem.MolFromSmarts(pat_to_check)
+  matches = mol.GetSubstructMatches(pattern)
+  for match in matches:
+    if len(match) == 1:
+      Search_pattern.append(match[0])
+
+  if atm_to_check in Search_pattern:
+    is_present = True
+  else:
+    is_present = False
+
+  return is_present, Search_pattern
+
+
+def findConnectedAtomsGeneral(mol, atomsToadd, atomidx, neig): 
+
+  """
+  This function attempt to find atoms connected to triple bond.
+  By default, it will cut just after an atom that does not have triple bond:
+  if the atom is part of triple bond, then stop after max two iterations.
+
+  TODO:
+    - If atoms neig as more than 1 neighboor, then need to add more options
+
+  Inputs:
+    -   mol: RDKIT mol object
+    -   atomsToadd: list of atom neighboor to main atom to check
+    -   atomidx: main atom to check ID
+    -   neig: neighboor atom of the main atom to check
+
+  Outputs:
+    -   Co_atoms: list of connected atoms
+  """
+
+  Break = False
+  prev_atm_idx = atomidx - 1
+  current_atm_idx = neig
+  Co_atoms = []
+  cc = 0
+
+  # While loop until the check for connected atoms stops
+  while not Break:
+
+    # Grab the neighboor atoms of neig
+    Neig = get_neigh_atoms(mol, current_atm_idx, prev_atm_idx)
+    print(f'Neig of {current_atm_idx}: {Neig}')
+    if len(Neig) == 1 and list(Neig.keys())[0] not in Co_atoms:
+        prev_atm_idx = current_atm_idx
+        current_atm_idx = list(Neig.keys())[0]
+
+        #if find_patern(mol,current_atm_idx,'[*;R]')[0]:
+        #  atomsToAdd = list(Neig.keys())
+        #  connectedAtoms = findConnectedAtoms(mol, atomsToAdd, find_patern(mol,current_atm_idx,'[*;R]')[1])
+        #  atomsToAdd += connectedAtoms
+        #  Co_atoms += atomsToAdd
+        #  Break = True
+        #else:
+        Co_atoms += list(Neig.keys())
+        Break = False
+
+    else:
+        Co_atoms += list(Neig.keys())
+        Break = True
+
+    cc += 1
+    if cc == 1:
+        Break = True
+
+  return Co_atoms
+
+
+def get_neigh_atoms(mol, Current_Idx, Prev_Idx):
+
+  """
+  This function create a dictionnary such as:
+  {Atom ID: Atomic Num}
+  Make sure that none of previous atom are selected
+
+  Inputs:
+    -   mol: RDKIT mol object
+    -   Current_Idx: Atom id to check
+    -   Prev_Idx: Atom id of the previous atom (to skip)
+
+  Outputs:
+    -   Neighboor: dictionnary with the neighboored atoms 
+
+  """
+  C_atm = mol.GetAtomWithIdx(Current_Idx)
+
+  Neighboor = {}
+  for neib in C_atm.GetNeighbors():
+    if neib.GetIdx() != Prev_Idx:
+      IDx = mol.GetAtomWithIdx(neib.GetIdx())
+      Neighboor[neib.GetIdx()] = IDx.GetAtomicNum()
+
+  return Neighboor
+
+
 # helper function to determine where to cut
 def findConnectedAtoms(rdkitmol, atomList, atomsInRing):
   # C-C single bond 
@@ -267,5 +433,7 @@ def specialCaseSixFusedFiveMemRing(fragmol):
 
   global sdffile
   sdffile = sys.argv[1]
-  atom_id = int(sys.argv[2]) 
+  atom_id = int(sys.argv[2])
+  print(sdffile)
+  print(atom_id) 
   growFragment(atom_id, sdffile)

PoltypeModules/lmodifytinkerkey.py

@@ -241,6 +241,8 @@ def modkey2_fragmpole(poltype):
           pt.write(f"numproc={poltype.numproc}\n")
           pt.write(f"maxmem={poltype.maxmem}\n")
           pt.write(f"maxdisk={poltype.maxdisk}\n")
+          if check_if_pattern_exist(f'./{fname}/{f}','[$(*#*)]'):
+            pt.write(f"new_gdma=True\n")
 
         # Run Poltype Job
         os.chdir(os.path.join(homedir, 'Fragments_DMA', fname))
@@ -269,6 +271,23 @@ def modkey2_fragmpole(poltype):
     shutil.move('tmpkeyforfrag.key', key2)
   return
 
+
+def check_if_pattern_exist(sdffile,pattern_to_check):
+
+    mol = Chem.MolFromMolFile(sdffile,removeHs=False)
+
+    atomsInPattern = []
+    pattern = Chem.MolFromSmarts(pattern_to_check)
+    matches = mol.GetSubstructMatches(pattern)
+    for match in matches:
+        if len(match) == 1:
+            atomsInPattern.append(match[0])
+
+    if len(atomsInPattern) > 0:
+        return True
+    else:
+        return False 
+
 # helper function to transfer multipole back to parent
 def transferMultipoleToParent(poltype, frag_sdf, frag_xyz, frag_key, parent_sdf, parent_xyz, parent_key):