Add files via upload

update the supporting for Pb

autosolvate/autoMCPB.py

@@ -6,6 +6,7 @@
 from molSimplify.Classes.ligand import ligand_breakdown
 import io
 from contextlib import redirect_stdout, redirect_stderr
+import networkx
 
 import os
 import re
@@ -17,15 +18,18 @@
 atom_zoo = ['H','C','N','O','CL','BR','I','P','S','F']
 
 
-metals  = [
-    "Li", "Na", "K", "Rb", "Cs", "Fr", "Pb",
-    "Be", "Mg", "Ca", "Sr", "Ba", "Ra",  "Sn","La"
-    "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",  
-    "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd",  
-    "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",  
-    "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn" 
+metals = [
+    "Li", "Be", "Na", "Mg", "Al", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", 
+    "Co", "Ni", "Cu", "Zn", "Ga", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", 
+    "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", 
+    "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", 
+    "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "Fr", "Ra", "Ac", 
+    "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", 
+    "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc", 
+    "Lv", "Ts"
 ]
 
+
 cheat_metal = {1:'Cu',2:'Fe',3:'Fe',4:'Mn'}
 
 valence_electrons = {
@@ -102,12 +106,33 @@ def check_ligand(sdf):
     return check, ligand_charge
 
 
-def ligandbreakdown_graph(xyzfile):
+def ligandbreakdown_BFGS(xyzfile,metalID):
     prefix = os.path.splitext(xyzfile)[0]
     sdf = prefix + '.sdf'
-    cmd = 'obabel ' + xyzfile + ' -O ' + prefix
-    pass
-
+    print(sdf)
+    cmd = 'obabel ' + xyzfile + ' -O ' + sdf
+    subprocess.call(cmd,shell=True)
+    with open(sdf,'r') as f:
+        data = f.readlines()
+        linker_dic = {}
+        type_dic = {}
+        connected = []
+        atomnumbers = int(data[3].split()[0])
+        bondsnumber = int(data[3].split()[1])
+        connectioninfo = data[4+atomnumbers:bondsnumber+4+atomnumbers]
+        atomtypes = data[4:4+atomnumbers]
+        G = networkx.Graph()
+        for line in connectioninfo:
+            atomx = int(line.split()[0]) - 1
+            atomy = int(line.split()[1]) - 1
+            if metalID not in [atomx,atomy]:
+                G.add_edge(atomx, atomy)
+        molecules = list(networkx.connected_components(G))
+        molecules_all = []
+        for i in molecules:
+            molecules_all.append(list(i))
+        return molecules_all
+
 
 class AtomInfo():
     def __init__(self,atomnum,atomtype,coordinate):
@@ -286,13 +311,16 @@ def ligand_identify(self):
         liglist = ligands_breakdown_infos[0]
        # denticity = ligands_breakdown_infos[1]
         ligcons =  ligands_breakdown_infos[2]
-        self.liglist = liglist
-        self.ligcons = ligcons
+        if len(liglist) > 0:
+            self.liglist = liglist
+            self.ligcons = ligcons
         #self.denticity = denticity
-        if len(liglist) == 0:
-            print('Warning: No ligand is identified by molsimplyfy, try to use graph search to break ligand')
+        elif len(liglist) == 0:
+            print('Warning: No ligand is identified by molsimplify, try to use graph search to break ligand')
+            liglist = ligandbreakdown_BFGS(self.xyzfile,int(self.metal_ID))
+            self.ligcons = 'NA'
+            self.liglist = liglist
 
-          #  sys.exit()
 
     def modifiy_pdb(self,inputpdb,mol2,outputpdb):
         r"""
@@ -671,9 +699,10 @@ def combine_pdbs(self):
             ligandname = 'LG' + str(ID)
             cmd = cmd + ' ' + ligandname + '.pdb '
         cmd = cmd + '> ' + temppdb
+        print('combine ' + cmd)
         subprocess.call(cmd,shell=True)
         cmd =self.amberhome + 'pdb4amber -i ' + temppdb + ' -o ' +  self.filename + '_final.pdb'
-        with open(ligandname + '_pdb4amber.log', 'w') as f:
+        with open(temppdb + '_pdb4amber.log', 'w') as f:
             subprocess.call(cmd, shell=True, stdout=f, stderr=subprocess.STDOUT)
 
     def get_bonded_pairs(self):
@@ -686,25 +715,36 @@ def get_bonded_pairs(self):
         -------
         None
         """
-        pdbin = self.filename + '_final.pdb'
-        xyzout = self.filename + '_final.xyz'
-        write(xyzout, read(pdbin))        
-        new_complex_mol = mol3D()
-        new_complex_mol.readfromxyz(self.filename + '_final.xyz',)
-        metalID = new_complex_mol.findMetal(transition_metals_only=False)
-        self.new_complex_mol = new_complex_mol
-        pairs = ligand_breakdown(new_complex_mol,transition_metals_only=False)[2]
+        if self.ligcons != 'NA':
+            pdbin = self.filename + '_final.pdb'
+            xyzout = self.filename + '_final.xyz'
+            write(xyzout, read(pdbin))        
+            new_complex_mol = mol3D()
+            new_complex_mol.readfromxyz(self.filename + '_final.xyz',)
+            metalID = new_complex_mol.findMetal(transition_metals_only=False)
+            self.new_complex_mol = new_complex_mol
+            pairs = ligand_breakdown(new_complex_mol,transition_metals_only=False)[2]
 
-        add_bonded_pairs = 'add_bonded_pairs '
-        for denticitys in pairs:
-            for denticity in denticitys:
-                metal_denticity = str(metalID[0]+1) + '-' + str(denticity+1) + ' '
-                add_bonded_pairs = add_bonded_pairs + metal_denticity
-
-        self.add_bonded_pairs = add_bonded_pairs
+            add_bonded_pairs = 'add_bonded_pairs '
+            for denticitys in pairs:
+                for denticity in denticitys:
+                    metal_denticity = str(metalID[0]+1) + '-' + str(denticity+1) + ' '
+                    add_bonded_pairs = add_bonded_pairs + metal_denticity
+
+            self.add_bonded_pairs = add_bonded_pairs
+        else:
+            self.add_bonded_pairs = ''
 
     def generate_MCPB_input(self):
         metalname = self.atomsinfo[self.metal_ID].atomtype
+        if self.ligcons == 'NA':
+            pdbin = self.filename + '_final.pdb'
+            xyzout = self.filename + '_final.xyz'
+            write(xyzout, read(pdbin))        
+            new_complex_mol = mol3D()
+            new_complex_mol.readfromxyz(self.filename + '_final.xyz',)
+            metalID = new_complex_mol.findMetal(transition_metals_only=False)
+            self.new_complex_mol = new_complex_mol
         metalID = self.new_complex_mol.findMetal(transition_metals_only=False)
         ion_ids = 'ion_ids ' + str(metalID[0]+1)
         self.ion_ids = ion_ids
@@ -727,7 +767,8 @@ def generate_MCPB_input(self):
             f.write('group_name ' + self.filename + '\n')
             f.write('cutoff 2.8\n')
             f.write(ion_ids + '\n')
-            f.write(self.add_bonded_pairs + '\n')
+            if self.ligcons != 'NA':
+                f.write(self.add_bonded_pairs + '\n')
             f.write(ion_mol2files+ '\n')
             f.write(naa_mol2files+ '\n')
             f.write(frcmod_files + '\n')
@@ -1273,7 +1314,12 @@ def build(self):
         elif self.round in ['4','4b','4n1','4n2']:
             self.run_MCPB_4()
         elif self.round in ['5','check','CHECK']:
-            self.checkingFF()
+            self.coordinates_reader_xyz()
+            self.metal_list()
+            self.check_atoms()
+            self.ligand_identify()
+            if self.ligcons != 'NA':
+                self.checkingFF()
 
 def startautoMCPB(argumentList):
     options = "hn:c:u:m:f:s:x:A:"