From 10636228af7e631c75b1f596acbd2777ca2e996d Mon Sep 17 00:00:00 2001 From: vgapsys Date: Wed, 14 Jul 2021 11:45:22 +0200 Subject: [PATCH] updated cyclic peptide builder script --- src/pmx/scripts/build_cyclic_top.py | 143 +++++++++++++++++++++++++++- 1 file changed, 140 insertions(+), 3 deletions(-) diff --git a/src/pmx/scripts/build_cyclic_top.py b/src/pmx/scripts/build_cyclic_top.py index 104f7a96..fffff797 100755 --- a/src/pmx/scripts/build_cyclic_top.py +++ b/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' )