updated cyclic peptide builder script

src/pmx/scripts/build_cyclic_top.py

@@ -90,6 +90,113 @@ def renumber_atoms( topol, offset ):
 def remove_residues( topol ):
     del topol.residues[0]
     del topol.residues[-1]
+
+# add angles between the first and last residues
+def add_angles( topol, a1, a2, partners ):
+    # get 1-2-3 connections
+    list13_a1 = []
+    list13_a2 = []
+    get13( list13_a1, a1, partners, topol )
+    get13( list13_a2, a2, partners, topol )
+
+    # add angles
+    for a in list13_a1:
+        ang = [a[0],a[1],a[2],1]
+        ang_inv = [a[2],a[1],a[0],1]
+        if (ang in topol.angles) or (ang_inv in topol.angles):
+            continue
+        else:
+            topol.angles.append(ang)
+    for a in list13_a2:
+        ang = [a[0],a[1],a[2],1]
+        ang_inv = [a[2],a[1],a[0],1]
+        if (ang in topol.angles) or (ang_inv in topol.angles):
+            continue
+        else:
+            topol.angles.append(ang)
+
+# add dihedrals between the first and last residues
+def add_dihedrals( topol, a1, a2, partners ):
+    # get 1-2-3-4 connections
+    list14_a1 = []
+    list14_a2 = []
+    get14( list14_a1, a1, partners, topol )
+    get14( list14_a2, a2, partners, topol )
+    # also collect for partners of a1 and a2
+    for a in partners[a1]:
+        if a!=a2:
+            get14( list14_a1, a, partners, topol )
+    for a in partners[a2]:
+        if a!=a1:
+            get14( list14_a2, a, partners, topol )
+
+    # add dihedrals
+    for a in list14_a1:
+        dih = [a[0],a[1],a[2],a[3],9,'']
+        dih_inv = [a[3],a[2],a[1],a[0],9,'']
+        if (dih in topol.dihedrals) or (dih_inv in topol.dihedrals):
+            continue
+        else:
+            topol.dihedrals.append(dih)
+    for a in list14_a2:
+        dih = [a[0],a[1],a[2],a[3],9,'']
+        dih_inv = [a[3],a[2],a[1],a[0],9,'']
+        if (dih in topol.dihedrals) or (dih_inv in topol.dihedrals):
+            continue
+        else:
+            topol.dihedrals.append(dih)
+
+# add impropers (type 4 for amber)
+def add_improper_dihedrals( topol, a1, a2, partners, ff='amber99sb-star-ildn'):
+    resnameLast = a2.resname
+    resnameFirst = a1.resname
+    firstResNum = 1
+    lastResNum = len(topol.residues)
+
+    a_N_first = a1
+    a_C_last = a2
+
+    # first dihedral: Ca_last - N_first - C_last - O_last
+    a_Ca_last = ''
+    a_O_last = ''
+    for a in topol.atoms:
+        if a.resnr==lastResNum and a.name=='CA':
+            a_Ca_last = a
+        if a.resnr==lastResNum and a.name=='O':
+            a_O_last = a
+    dih = [a_Ca_last,a_N_first,a_C_last,a_O_last,4,'']
+    topol.dihedrals.append(dih)
+
+    # second dihedral: C_last - Ca_first - N_first - H_first
+    # (if Pro: C_last - CD_first - N_first - Ca_first  )
+    dih = []
+    a_Ca_first = ''
+    a_H_first = ''
+    a_CD_first = ''
+    for a in topol.atoms:
+        if a.resnr==firstResNum and a.name=='CA':
+            a_Ca_first = a
+        if a.resnr==firstResNum and a.name=='H':
+            a_H_first = a
+        if a.resnr==firstResNum and a.name=='CD':
+            a_CD_first = a
+    if a_N_first.resname!='PRO':
+        dih = [a_C_last,a_Ca_first,a_N_first,a_H_first,4,'']
+    else:
+        dih = [a_C_last,a_CD_first,a_N_first,a_Ca_first,4,'']
+
+    topol.dihedrals.append(dih)
+
+    # for amber-star: N_last - Ca_last - C_last - N_first
+    if ('star' in ff) or ('*' in ff):
+        if a_C_last.resname!='GLY':
+            a_N_last = ''
+            for a in topol.atoms:
+                if a.resnr==lastResNum and a.name=='N':
+                    a_N_last = a
+            dih = [a_N_last,a_Ca_last,a_C_last,a_N_first,4,'105.4 0.75 1']
+            topol.dihedrals.append(dih)
+
 
 # add a bond to close the cycle
 def add_bond( topol ):
@@ -139,7 +246,10 @@ def add_bond( topol ):
     for p in pairs:
         if [p[0],p[1]] not in list13:
             topol.pairs.append(p)
+
+    return(a1,a2,partners)
 
+# parse 1-4 interactions to construct pairs
 def parse14( pairs, list13, a, partners, topol ):
     for a12 in partners[a]:
         for a13 in partners[a12]:
@@ -158,11 +268,32 @@ def parse14( pairs, list13, a, partners, topol ):
                     continue
                 pairs.append(p)
 
+# gets the list of 1-2-3 connections (for angles)
+def get13( list13, a, partners, topol ):
+    for a12 in partners[a]:
+        for a13 in partners[a12]:
+            if a13==a:
+                continue
+            list13.append( [a,a12,a13] )
+#            list13.append( [a13,a12,a] )
+
+# gets the list of 1-4 connections (for dihedrals)
+def get14( list14, a, partners, topol ):
+    for a12 in partners[a]:
+        for a13 in partners[a12]:
+            if a13==a:
+                continue
+            for a14 in partners[a13]:
+                if a12==a14 or a==a14:
+                    continue
+                list14.append( [a,a12,a13,a14] )
+
 def main(argv):
 
    options=[
    Option( "-seq", "string", "PVWLVVV" , "peptide sequence"),
    Option( "-chir", "string", "LDLDLDL" , "peptide sequence"),
+   Option( "-ff", "string", "amber99sb-star-ildn" , "for now only amber is supported"),
         ]
 
    files = [
@@ -184,7 +315,7 @@ def main(argv):
 #   iname = cmdl['-f']
    oname = cmdl['-o']
    otopname = cmdl['-p']
-#   ff = cmdl['-ft']
+   ff = cmdl['-ff']
 
    # 1. construct the peptide, e.g XYZ
    m = build_chain(cmdl['-seq'],bCyclic=True,chirality=cmdl['-chir'])
@@ -215,8 +346,14 @@ def main(argv):
    remove_atoms( topol, removeID )   
    renumber_atoms( topol, offset )
    remove_residues( topol )
-   add_bond( topol )
-
+   a1,a2,partners = add_bond( topol ) # connecting bond and pairs
+   # a1 is N in residue 1
+   # a2 is C in the last residue
+   # partners is a dictionary of bonded partners for each atom: partners[atom] = [a1,a2,a3,...]
+   add_angles( topol, a1, a2, partners ) # connecting angles
+   add_dihedrals( topol, a1, a2, partners ) # connecting dihedrals
+   add_improper_dihedrals( topol, a1, a2, partners, ff=ff) # impropers
+
    # 5. modify the structure: remove first and last residues
    # read the structure from pdb2gmx
    m = Model( 'out.pdb' )