diff --git a/PoltypeModules/lFragmenterForDMA.py b/PoltypeModules/lFragmenterForDMA.py index d69e1e10..6e4cbc0a 100644 --- a/PoltypeModules/lFragmenterForDMA.py +++ b/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) diff --git a/PoltypeModules/lmodifytinkerkey.py b/PoltypeModules/lmodifytinkerkey.py index f4b074a6..7180cffb 100644 --- a/PoltypeModules/lmodifytinkerkey.py +++ b/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):