check_dihedral.py

#!/usr/bin/python
import sys, datetime, math
import numpy as np

def ReadLinePDBform(line):
    ##01234567890123456789012345678901234567890123456789012345678901234567890123456789
    ##ATOM      1  N   DOPC    1      -5.364  -3.265  20.254  1.00  0.00      MEMB
    line_type = str(line[0:6].strip())
    ret = []    
    if line_type == "ATOM" or line_type == "HETATM":
        ret.append(str(line[12:16].strip())) #atom name   
        ret.append(str(line[16:21].strip())) #molecular name
        ret.append(int(line[22:30]))         #molecular number
        ret.append(float(line[30:38]))       #x coordinate
        ret.append(float(line[38:46]))       #y coordinate
        ret.append(float(line[46:54]))       #z coordinate
        ret.append(float(line[54:60]))
        ret.append(float(line[60:66]))
        ret.append(str(line[66:76].strip()))
        return ret
    return ret

def ReadLineGROform(line):
    #012345678901234567890123456789012345678901234567890
    #    1DOPC     N    1   6.661   5.287   5.129  0.1628  0.8558 -0.1599
    ret = [str(line[10:15]).strip(), str(line[5:10]).strip(), int(line[:5]), float(line[20:28]), float(line[28:36]), float(line[36:44]), float(line[44:52]), float(line[52:60]), float(line[60:68])]
    return ret

#calculate dihedral (radian) 
def CalcDihedral(atom_i, atom_j, atom_k, atom_l):
    posi = np.array( [atom_i[3], atom_i[4], atom_i[5]] )
    posj = np.array( [atom_j[3], atom_j[4], atom_j[5]] )
    posk = np.array( [atom_k[3], atom_k[4], atom_k[5]] )
    posl = np.array( [atom_l[3], atom_l[4], atom_l[5]] )
    
    i2j = posj - posi
    j2k = posk - posj
    k2l = posl - posk
    
    ijk_v = np.cross(i2j, j2k)
    jkl_v = np.cross(j2k, k2l)

    ijk_v = ijk_v / np.linalg.norm(ijk_v)
    jkl_v = jkl_v / np.linalg.norm(jkl_v)
    
    return math.acos(np.dot(ijk_v, jkl_v))
    
class RetFlag:
    Target = 0
    SkipObject = 1
    Write = 2

def IsTargetCarbon(atom_d, dbond_id, tar_lipname):
    for i in range(len(dbond_id)):
        car_atoms = ("C2" + str(dbond_id[i]), "C3" + str(dbond_id[i]), "C2" + str(dbond_id[i] + 1), "C3" + str(dbond_id[i] + 1))
        hyd_atoms = ("H" + str(dbond_id[i]) + "X", "H" + str(dbond_id[i] + 1) + "X", "H" + str(dbond_id[i]) + "R", "H" + str(dbond_id[i] + 1) + "R")
        flag_tar  = ( (atom_d[0] == car_atoms[0]) or (atom_d[0] == car_atoms[1]) ) and (atom_d[1] == tar_lipname)
        flag_skip = ( (atom_d[0] == car_atoms[2]) or (atom_d[0] == car_atoms[3]) )
        for j in range(4):
            flag_skip = flag_skip or (atom_d[0] == hyd_atoms[j])
        
        flag_skip = flag_skip and (atom_d[1] == tar_lipname)
        
        if flag_tar:
            return RetFlag.Target
        elif flag_skip:
            return RetFlag.SkipObject
        
    return RetFlag.Write

def ReadAllPDB(fin):
    atoms_d = []
    for line in fin:
        atom_d = ReadLinePDBform(line)
        if(len(atom_d) != 0):
            atoms_d.append(atom_d)
    return atoms_d

def ReadAllGRO(fin):
    lines = fin.readlines()
    atoms_d = []
    for line in lines[2:-1]:
        atoms_d.append(ReadLineGROform(line) )
    return atoms_d

def CheckAllAtoms(atoms_d, dbond_id, tar_lipid_name):
    atom_id = 1
    for i in range(len(atoms_d)):
        flag = IsTargetCarbon(atoms_d[i], dbond_id, tar_lipid_name)
        if flag == RetFlag.Target:
            phi = CalcDihedral(atoms_d[i+1], atoms_d[i], atoms_d[i+2], atoms_d[i+3])
            assert phi < (0.5 * math.pi)
            atom_id += 4
        elif flag == RetFlag.Write:
            atom_id += 1
    
def ConverStrs2Int(array):
    for i in range(len(array)):
        array[i] = int(array[i])
    return array

if __name__ == "__main__":
    argvs = sys.argv
    if(len(argvs) < 4):
        print "Usage $ python %s input_filename target_lipid_name positions_of_double_bond" % argvs[0]
        quit()
    print "Check dihedral angle"

    dbond_id       = ConverStrs2Int(argvs[3:len(argvs)])
    fin            = open(argvs[1], "r")
    tar_lipid_name = argvs[2]

    extent = argvs[1][len(argvs[1]) - 3:]

    if extent == "pdb":
        atoms_d = ReadAllPDB(fin)    
    elif extent == "gro":
        atoms_d = ReadAllGRO(fin)
    else:
        print "Unknown file extension %s" % extent 

    CheckAllAtoms(atoms_d, dbond_id, tar_lipid_name)
    print "Success"