insert_dbonds.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

def WriteLinePDBform(atom_d, atom_id, fout, new_lipid_name, tar_lipid_name):
    if atom_d[1] == tar_lipid_name:
        if atom_id < 100000:
            fout.write("%-6s%5d%5s %-5s%-8d%8.3f%8.3f%8.3f%6.2f%6.2f      %s\n"%("ATOM", atom_id, atom_d[0], new_lipid_name, atom_d[2], atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]))
        else:
            fout.write("%-6s%5s%5s %-5s%-8d%8.3f%8.3f%8.3f%6.2f%6.2f      %s\n"%("ATOM", "*****", atom_d[0], new_lipid_name, atom_d[2], atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]))
    else:
        if atom_id < 100000:
            fout.write("%-6s%5d%5s %-5s%-8d%8.3f%8.3f%8.3f%6.2f%6.2f      %s\n"%("ATOM", atom_id, atom_d[0], atom_d[1], atom_d[2], atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]))
        else:
            fout.write("%-6s%5s%5s %-5s%-8d%8.3f%8.3f%8.3f%6.2f%6.2f      %s\n"%("ATOM", "*****", atom_d[0], atom_d[1], atom_d[2], atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]))

def WriteLineGROform(atom_d, atom_id, fout, new_lipid_name, tar_lipid_name):
    atom_id = atom_id % 100000
    assert atom_d[2] < 100000, "Max # of molecules is 99999."
    if atom_d[1] == tar_lipid_name:
        fout.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n" % (atom_d[2], new_lipid_name, atom_d[0], atom_id, atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]) )
    else:
        fout.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n" % (atom_d[2], atom_d[1], atom_d[0], atom_id, atom_d[3], atom_d[4], atom_d[5], atom_d[6], atom_d[7], atom_d[8]) )

#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]) + "Y", "H" + str(dbond_id[i] + 1) + "X", "H" + str(dbond_id[i] + 1) + "Y", "H" + str(dbond_id[i]) + "R", "H" + str(dbond_id[i]) + "S", "H" + str(dbond_id[i] + 1) + "R", "H" + str(dbond_id[i] + 1) + "S")
        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(8):
            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(lines):
    atoms_d = []
    for line in lines[2:-1]:
        atoms_d.append(ReadLineGROform(line) )
    return atoms_d

def ReadBoxGRO(lines):
    return lines[-1]

def ReviseHydroName(name):
    if name[-1] == "Y":
        new_name = name[0:-1] + "X"
        return new_name
    elif name[-1] == "S":
        new_name = name[0:-1] + "R"
        return new_name
    else:
        print "Unknwon hydro name " + name + "."
        sys.exit()
    
def WriteAll(fout, atoms_d, dbond_id, new_lipid_name, tar_lipid_name, writefunc):
    atom_id = 1
    for i in range(len(atoms_d)):
        flag = IsTargetCarbon(atoms_d[i], dbond_id, tar_lipid_name)
        
        if flag == RetFlag.Target:
            #phi1034, phi1035, phi2034, phi2035
            phi = CalcDihedral(atoms_d[i+1], atoms_d[i], atoms_d[i+3], atoms_d[i+4]), CalcDihedral(atoms_d[i+1], atoms_d[i], atoms_d[i+3], atoms_d[i+5]), CalcDihedral(atoms_d[i+2], atoms_d[i], atoms_d[i+3], atoms_d[i+4]), CalcDihedral(atoms_d[i+2], atoms_d[i], atoms_d[i+3], atoms_d[i+5])

            min_phi = min(phi)
            print phi
            assert min_phi < (0.5 * math.pi)
            
            min_idx = phi.index(min_phi)
            writefunc(atoms_d[i], atom_id, fout, new_lipid_name, tar_lipid_name)
            atom_id += 1
            
            if (min_idx / 2) == 0:
                writefunc(atoms_d[i+1], atom_id, fout, new_lipid_name, tar_lipid_name)
                atom_id += 1
            else:
                atoms_d[i+2][0] = ReviseHydroName(atoms_d[i+2][0])
                writefunc(atoms_d[i+2], atom_id, fout, new_lipid_name, tar_lipid_name)
                atom_id += 1
            
            writefunc(atoms_d[i+3], atom_id, fout, new_lipid_name, tar_lipid_name)
            atom_id += 1
            
            if (min_idx % 2) == 0:
                writefunc(atoms_d[i+4], atom_id, fout, new_lipid_name, tar_lipid_name)
                atom_id += 1
            else:
                atoms_d[i+5][0] = ReviseHydroName(atoms_d[i+5][0])
                writefunc(atoms_d[i+5], atom_id, fout, new_lipid_name, tar_lipid_name)
                atom_id += 1
        elif flag == RetFlag.Write:
            writefunc(atoms_d[i], atom_id, fout, new_lipid_name, tar_lipid_name)
            atom_id += 1
        elif flag == RetFlag.SkipObject:
            print "skip atom"
            
    return (atom_id, atoms_d[-1][1], atoms_d[-1][2])
    
def ConverStrs2Int(array):
    for i in range(len(array)):
        array[i] = int(array[i])
    return array

def WriteRemarkPDB(fout):
    fout.write("REMARK  CREATED BY Koh NAKAGAWA(http://noguchi.issp.u-tokyo.ac.jp/nakagawa/index.html)\n")
    fout.write("REMAKR  DATE: %s\n" % datetime.datetime.today())

def WriteEndPDB(fout, last_atom_id, last_mol_id, last_mol_name):
    fout.write("%-6s%5d      %-5s%4d\n" % ("TER", last_atom_id, last_mol_name, last_mol_id))
    fout.write("END")

def WriteTitleGRO(fout, new_lipid_name):
    fout.write("%s in Water\n" % (new_lipid_name) )

def WriteEndGRO(fout, box):
    fout.write(box)

def WriteAtomNum(fout, atom_num):
    fout.write("%d\n", atom_num)

if __name__ == "__main__":
    argvs = sys.argv
    if(len(argvs) < 6):
        print "Usage $ python %s input_filename output_filename converted_lipid_name new_lipid_name positions_of_double_bond" % argvs[0]
        quit()
    print "Remove some H atoms and introduce some double bonds."
    print "Convert %s to %s and remove some bonds." % (argvs[3], argvs[4])

    dbond_id       = ConverStrs2Int(argvs[5:len(argvs)])
    fin            = open(argvs[1], "r")
    fout           = open(argvs[2], "w")
    tar_lipid_name = argvs[3]
    new_lipid_name = argvs[4]

    extent = argvs[1][len(argvs[1]) - 3:]
    if extent == "pdb":
        #PDB format
        atoms_d = ReadAllPDB(fin)
        WriteRemarkPDB(fout)
        last_atom_id, last_mol_name, last_mol_id = WriteAll(fout, atoms_d, dbond_id, new_lipid_name, tar_lipid_name, WriteLinePDBform)
        WriteEndPDB(fout, last_atom_id, last_mol_id, last_mol_name)
    elif extent == "gro":
        #GRO format
        lines   = fin.readlines()
        atoms_d = ReadAllGRO(lines)
        box     = ReadBoxGRO(lines)
        WriteTitleGRO(fout, new_lipid_name)
        last_atom_id, last_mol_name, last_mol_id = WriteAll(fout, atoms_d, dbond_id, new_lipid_name, tar_lipid_name, WriteLineGROform)
        WriteEndGRO(fout, box)
        fout.close() #flush
        
        lines = open(argvs[2], "r").readlines()
        lines.insert(1, "%d\n" % (last_atom_id - 1))
        lines = "".join(lines)
        open(argvs[2], "w").write(lines)
    else:
        print "Unknown file extension %s" % extent