Add files via upload

1_post_process_traj_v1.py

@@ -0,0 +1,38 @@
+#!/home/magarkar/anaconda3/envs/simba/bin/python
+
+import MDAnalysis as mda
+import subprocess,os
+
+gro='prod.gro'
+xtc='prod.xtc'
+
+u=mda.Universe(gro,xtc)
+
+def postprocess_gromacs(xtc):
+    if os.path.isfile(xtc) and os.path.isfile('prod.tpr'):
+        subprocess.getoutput("echo 0 | gmx trjconv -s prod.tpr -f %s -o whole.xtc -pbc whole"%xtc)
+        subprocess.getoutput("echo 0 | gmx trjconv -s prod.tpr -f whole.xtc -o nojump.xtc -pbc nojump")
+        subprocess.getoutput("echo Protein System | gmx trjconv -s prod.tpr -f nojump.xtc -o center.xtc -pbc mol -center")
+        subprocess.getoutput("echo C-alpha System | gmx trjconv -s prod.tpr -f center.xtc -o traj.xtc -fit rot+trans")
+        subprocess.getoutput("echo C-alpha System | gmx trjconv -s prod.tpr -f center.xtc -o traj.gro -e 0 -fit rot+trans")
+    else:
+        print("file not found")
+
+###############################################################################################################################
+
+
+total_frames=u.trajectory.n_frames
+number_of_frames=50
+
+if total_frames>number_of_frames:
+    interval=int(total_frames/number_of_frames)
+    fiter = u.trajectory[::interval]
+    frames = [ts.frame for ts in fiter]
+    u.atoms.write('short.xtc', frames=frames)
+    postprocess_gromacs("short.xtc")
+else:
+    postprocess_gromacs("prod.xtc")
+
+if os.path.isfile("traj.xtc") and os.path.isfile("traj.gro"):
+    print("Cleaning UP")
+    subprocess.getoutput("rm center.xtc whole.xtc short.xtc nojump.xtc")

2_dg_calc_v2.py

@@ -0,0 +1,51 @@
+#!/home/magarkar/anaconda3/bin/python
+import parmed as pmd
+import subprocess, os
+import MDAnalysis as mda
+
+#1 Process gromacs trajectory
+#print("Step 1/6 Postprocessing trajectory")
+
+#subprocess.getoutput("echo 0 | gmx trjconv -s prod.tpr -f prod.xtc -o whole.xtc -pbc whole -e 1000")
+#subprocess.getoutput("echo 0 | gmx trjconv -s prod.tpr -f whole.xtc -o nojump.xtc -pbc nojump")
+#subprocess.getoutput("echo Protein System | gmx trjconv -s prod.tpr -f nojump.xtc -o energy.xtc -pbc mol -center")
+#subprocess.getoutput("echo Protein System | gmx trjconv -s prod.tpr -f nojump.xtc -o energy.gro -e 0 -pbc mol -center")
+
+#2 parmed to convert gromacs to AMBER
+print("Step 2/6 Convering to AMBER topology")
+
+top = pmd.load_file("../topol2.top",defines=dict(DEFINE='DUMMY'))
+struct = pmd.load_file("traj.gro")
+
+top.save('complex.parm7',overwrite=True)
+struct.save('complex.rst7',overwrite=True)
+
+#3 MDAnalysis xtc to netcdf
+print("Step 3/6 Convering to AMBER trajectory")
+
+u = mda.Universe('traj.gro', 'traj.xtc')
+u.atoms.write('energy.ncdf', frames='all')
+
+#4 Prepare trajectory with waters
+print("Step 4/6 Preparing trajectory with waters")
+
+subprocess.getoutput("cpptraj -p complex.parm7 -i /data/projects/SIMBA_PROJECT/scripts/dg_calc/closest20.in")
+
+#5 Separate components
+print("Step 5/6 Separating topology components")
+
+subprocess.getoutput('ante-MMPBSA.py -p N20.complex.parm7 -c complex_nowat.prmtop -r protein.prmtop -l ligand.prmtop -s ":NA,CL,Na+,Cl-" -n ":MOL" --radii=mbondi2') 
+
+#6 Calculate DG
+print("Step 6/6 Running MMGBSA")
+
+subprocess.getoutput("MMPBSA.py -O -i /data/projects/SIMBA_PROJECT/scripts/dg_calc/mmbpsa.in -o result.dat -sp N20.complex.parm7 -cp complex_nowat.prmtop -rp protein.prmtop -lp ligand.prmtop -y N20.nc")
+
+subprocess.getoutput("rm _MM* *#*")
+#subprocess.getoutput("rm whole.xtc nojump.xtc energy.ncdf N20.nc")
+
+u2 = mda.Universe("N20.complex.parm7","N20.nc")
+u2.atoms.write('for_simba.xtc', frames='all')
+u2.atoms.write('for_simba.gro', frames=[0])
+
+subprocess.getoutput("/home/magarkar/anaconda3/envs/gromacs/bin/gmx editconf -f for_simba.gro -o for_simba.pdb")