For Setting up Gromacs-Protein Ligand System

combineGro_prot_lig.py

## Usage: python combineGro_prot_lig.py protein_clean.pdb UNK.pdb 
import sys
import os 

def GetCharge(): 
    lines = open('LOG_CHARGE').readlines()
    charge = []
    for line in lines: 
        if 'System has non-zero total charge' in line: 
            charge.append(float(line.split()[-1]))
#    print(charge[0])
    return(charge[0])
     

def getNumProtResid(fl):
    nres = int(fl[0:5])
    return(nres+1)

def WriteComplexGro(pro_dat,lig_dat):
    ofile=open('complex.gro','w')
    ofile.write('%s\n'%pro_dat[0])
    ofile.write('%5d\n'%tot_num)
    for line in pro_dat[2:-1]: ofile.write('%s\n'%line)
    for line in lig_dat[2:-1]: ofile.write('%s\n'%('%5d'%res_num+line[5:]))
    ofile.write('%s\n'%pro_dat[-1])
    ofile.close()
    return None

def EditTopolFile(resid):
    topol = open('topol.top','r').readlines()
    olines = []
    for line in topol: 
        olines.append(line)
        if 'forcefield.itp' in line: olines.append('#include \"%s.itp\"'%resid) 
    olines.append('%-15s     1\n'%resid)
    ofile = open('topol_PLIG.top','w')
    for line in olines: ofile.write('%s'%line)
    return None 

def SetupPLIGSystemGMX(): 
    ofile = open('ions.mdp','w')
    ofile.write('''
; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator	= steep		; Algorithm (steep = steepest descent minimization)
emtol		= 1000.0  	; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep          = 0.01          ; Energy step size
nsteps		= 50000	  	; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist		= 1		; Frequency to update the neighbor list and long range forces
ns_type		= grid		; Method to determine neighbor list (simple, grid)
rlist		= 1.0		; Cut-off for making neighbor list (short range forces)
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 1.0		; Short-range electrostatic cut-off
rvdw		= 1.0		; Short-range Van der Waals cut-off
pbc		= xyz 		; Periodic Boundary Conditions (yes/no)
'''
)
    ofile.close()
    commands = [
                'gmx editconf -f complex.gro  -o complex_box.gro -c -d 1.2 -bt cubic', 
                'gmx solvate -cp complex_box.gro -cs tip4p.gro -o complex_solvated.pdb -p topol_PLIG.top', 
                'gmx grompp -f ions.mdp -c complex_solvated.pdb -p topol_PLIG.top -o ions.tpr &> LOG_CHARGE',
                'echo 15| gmx genion -s ions.tpr -o complex_solvated.pdb -p topol_PLIG.top -pname NA -nname CL -neutral'
                ]
    for com in commands: 
        os.system(com)
    return None

pro_gro = sys.argv[1]
lig_gro = sys.argv[2]
lig_rid = sys.argv[2]
pro = open(pro_gro,'r').readlines()
lig = open(lig_gro,'r').readlines()
pro_dat = [line.rstrip() for line in pro]
lig_dat = [line.rstrip() for line in lig]
res_num = getNumProtResid(pro_dat[-2])
tot_num = int(pro_dat[1])+int(lig_dat[1])
## Write Protein-Ligand Complex GRO file 
WriteComplexGro(pro_dat,lig_dat)
## Edit topol.top to add ligand part 
EditTopolFile(resid=lig_rid)
## Create Solvated Protein-Ligand System 
SetupPLIGSystemGMX()