Merge pull request #25 from dptech-corp/topology

unigbsa-traj result bug

unigbsa/CLI.py

@@ -266,10 +266,8 @@ def traj_pipeline(args=None):
     pbsa = GBSA()
     pbsa.set_paras(complexfile=complexFile, trajectoryfile=trajFile, topolfile=topolFile, indexfile=indexFile, mmpbsafile=mmpbsafile, pbsaParas=pbsaParas, nt=nt)
     pbsa.run(verbose=debug)
-    detal_G = pbsa.extract_result()
-    print("mode    detal_G(kcal/mole)    Std. Dev.")
-    for k, v in detal_G.items():
-        print('%4s    %18.4f    %9.4f'%(k, v[0], v[1]))
+    delta_G = pbsa.extract_result()
+    print(delta_G)
 
 def mmpbsa_plot():
     from unigbsa.gbsa import plots

unigbsa/simulation/topology.py

@@ -55,7 +55,10 @@ def build_lignad(ligandfile, forcefield="gaff2", charge_method="bcc", engine="ac
         'net_charge': charge,
         'sqm_key': "grms_tol=0.005,qm_theory='AM1',scfconv=1.d-8,ndiis_attempts=700,maxcyc=0",
     }
-    cmd = '''export OMP_NUM_THREADS={thread};acpype -i {ligandfile} -b {molname} -a {forcefield} -c {method} -n {net_charge} -k "{sqm_key}" -f >acpype.log 2>&1 '''.format(**paras)
+    if charge is None:
+        cmd = '''export OMP_NUM_THREADS={thread};acpype -i {ligandfile} -b {molname} -a {forcefield} -c {method} -k "{sqm_key}" -f >acpype.log 2>&1 '''.format(**paras)
+    else:
+        cmd = '''export OMP_NUM_THREADS={thread};acpype -i {ligandfile} -b {molname} -a {forcefield} -c {method} -n {net_charge} -k "{sqm_key}" -f >acpype.log 2>&1 '''.format(**paras)
     RC = os.system(cmd)
     if RC != 0:
         print(cmd)
@@ -76,6 +79,67 @@ def build_lignad(ligandfile, forcefield="gaff2", charge_method="bcc", engine="ac
         shutil.rmtree(ligandName)
     return moltop, molgro
 
+
+def pdb2amber(pdbfile, outfile=None):
+    structure = pmd.load_structure(pdbfile)
+    atomlist = []
+    for c, residue in enumerate(structure.residues, start=1):
+        # change atoms name from GROMACS to AMBER
+        for atom in residue.atoms:
+            if atom.name == 'OC1':
+                atom.name = 'O'
+            elif atom.name == 'OC2':
+                atom.name = 'OXT'
+                residue.ter = True  # parmed terminal
+        # change residues name according to AMBER
+        if residue.name == 'ILE':
+            for atom in residue.atoms:
+                if atom.name == 'CD':
+                    atom.name = 'CD1'
+                    break
+        elif residue.name == 'LYS':
+            atoms = [atom.name for atom in residue.atoms]
+            if 'HZ3' not in atoms:
+                residue.name = 'LYN'
+        elif residue.name == 'ASP':
+            atoms = [atom.name for atom in residue.atoms]
+            if 'HD2' in atoms:
+                residue.name = 'ASH'
+        elif residue.name == 'GLU':
+            atoms = [atom.name for atom in residue.atoms]
+            if 'HE2' in atoms:
+                residue.name = 'GLH'
+        elif residue.name in ['HIS', 'HIE', 'HID', 'HIP']:
+            atoms = [atom.name for atom in residue.atoms if atom.atomic_number == 1]
+            if 'HD1' in atoms and 'HE2' in atoms:
+                residue.name = 'HIP'
+            elif 'HD1' in atoms:
+                residue.name = 'HID'
+            elif 'HE2' in atoms:
+                residue.name = 'HIE'
+        for atom in residue.atoms:
+            if atom.element > 2:
+                atomlist.append(atom)
+    structure.atoms = pmd.structure.AtomList(atomlist)
+    if outfile is None:
+        outfile = pdbfile
+    structure.write_pdb(outfile)
+
+
+def build_protein_tleap(pdbfile, forcefield='', outtop=None, outcoor=None):
+    struc = pmd.load_structure(pdbfile)
+    leaprc = '''
+    source leaprc.protein.ff19SB
+    source leaprc.gaff
+    loadOff atomic_ions.lib
+    loadamberparams frcmod.ions234lm_126_tip3p
+    set default PBRadii mbondi2
+    REC = loadpdb %s
+    REC_OUT = combine { REC }
+    saveamberparm REC_OUT REC.prmtop REC.inpcrd
+    ''' % pdbfile
+
+
 def build_protein(pdbfile, forcefield='amber99sb-ildn', outtop=None, outcoord=None):
     """
     Build a protein topology and coordinate file from a PDB file

unigbsa/simulation/utils.py

@@ -298,14 +298,14 @@ def obtain_net_charge(molfile):
     mol2file = str(uuid.uuid4()) + '.mol2'
     cmd = f'obabel -i {ftype} {molfile} -omol2 -O {mol2file} --partialcharge gasteiger >/dev/null 2>&1 '
     RC = os.system(cmd)
-    if RC!=0:
-        logging.warning('Failed to obtain mol charge, use guess')
+    if RC != 0:
+        print('Failed to obtain mol charge, use guess')
         return None
     charge = 0
     with open(mol2file) as fr:
         for line in fr:
-            if line.startswith('charge'):
-                charge += int(line.split()[1].strip())
-    return charge
-# def obtain_net_charge(molfile):
+            llist = line.split()
+            if len(llist) >= 8:
+                charge += float(llist[-1])
+    return int(round(charge))
 

unigbsa/version.py

@@ -1,2 +1,2 @@
-__version__="0.1.1"
+__version__="0.1.2"