-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPrepare.py
executable file
·76 lines (72 loc) · 2.51 KB
/
Prepare.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
import subprocess
import os
from simtk.openmm import app
## makeLib
# parametrizes a pdb or mol2 file to generate a .lib library file for tleap / nab
def makeLib(file_path, residue_name, connect0=None, connect1=None, charges='bcc', atom_type='gaff', force_field='leaprc.ff12SB', parameterized=False):
name, extension = file_path.split('/')[-1].split(".")
lib_path = '/'.join(file_path.split('/')[:-1]+[residue_name])
tleap_input = """
source %s
source leaprc.gaff
%s = loadmol2 %s.mol2
loadamberparams %s.frcmod"""%(force_field,
residue_name, name,
lib_path)
if connect0 and connect1:
tleap_input += """
set %s head %s.1.%s
set %s tail %s.1.%s
set %s.1 connect0 %s.head
set %s.1 connect1 %s.tail"""%(residue_name, residue_name, connect0,
residue_name, residue_name, connect1,
residue_name, residue_name,
residue_name, residue_name)
tleap_input +="""
check %s
saveoff %s %s.lib
savepdb %s %s_tmp.pdb
quit
"""%(residue_name, residue_name, lib_path, residue_name, lib_path)
if parameterized:
tleap_input = """
source leaprc.gaff
source %s
%s = loadpdb %s.pdb
saveoff %s %s.lib
savepdb %s %s_tmp.pdb
quit
"""%(force_field, residue_name, name, residue_name, lib_path,
residue_name, lib_path)
#Write LEaP infile
with open("%s.in"%name,"w") as fil:
fil.write(tleap_input)
if not parameterized:
#Execute
subprocess.call("conda run -n maws_p2 && antechamber -i %s -fi %s -o %s.mol2 -fo mol2 -c %s -rn %s -at %s"%(file_path,
extension,
name,
charges,
residue_name,
atom_type),
shell=True)
subprocess.call("conda run -n maws_p2 && parmchk2 -i %s.mol2 -f mol2 -o %s.frcmod"%(name, lib_path), shell=True)
subprocess.call("conda run -n maws_p2 && tleap -f %s.in"%name, shell=True)
PDB = app.PDBFile(lib_path + "_tmp.pdb")
length = sum([1 for atom in PDB.topology.atoms()])
#Cleanup
if not parameterized:
os.remove(name+".mol2")
#os.remove(name+".frcmod")
os.remove("%s.in"%name)
os.remove("%s_tmp.pdb"%(lib_path)) # lib_path = PDB
os.remove("leap.log")
return length
## toggleHydrogens
# Toggles presence/abscence of Hydrogens as indicated by
# @param boolean
# true: hydrogens or false: hydrogens
def toggleHydrogens(path, boolean=True):
subprocess.call("conda run -n maws_p2 && reduce -Trim %s > %s"%(path,path), shell=True)
if boolean:
subprocess.call("conda run -n maws_p2 && reduce -Build %s > %s"%(path,path), shell=True)