From 193df73a39f2c9569f18973c25e8445b9bd56a04 Mon Sep 17 00:00:00 2001 From: recombinatrix Date: Mon, 17 Feb 2025 22:49:46 +1100 Subject: [PATCH 1/2] extended polyconf examples --- .DS_Store | Bin 0 -> 6148 bytes .gitignore | 6 + polyconf/polyconf/Polymer.py | 19 +- polyconf_examples/01_build_PEI_linear.py | 298 ------------------ .../01a_build_PEI_with_extend.py | 80 +++++ .../01b_build_PEI_with_linear_extend.py | 141 +++++++++ polyconf_examples/01c_build_PEI_then_solve.py | 76 +++++ .../01d_build_PEI_conformation.py | 88 ++++++ .../01e_build_PEI_conformation_ensemble.py | 70 ++++ polyconf_examples/02a_build_PMMA_isotactic.py | 28 +- polyconf_examples/03_build_branched_PEI.py | 138 ++++++++ polyconf_examples/README.md | 14 +- 12 files changed, 631 insertions(+), 327 deletions(-) create mode 100644 .DS_Store delete mode 100644 polyconf_examples/01_build_PEI_linear.py create mode 100644 polyconf_examples/01a_build_PEI_with_extend.py create mode 100644 polyconf_examples/01b_build_PEI_with_linear_extend.py create mode 100644 polyconf_examples/01c_build_PEI_then_solve.py create mode 100644 polyconf_examples/01d_build_PEI_conformation.py create mode 100644 polyconf_examples/01e_build_PEI_conformation_ensemble.py create mode 100644 polyconf_examples/03_build_branched_PEI.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..9379249c38c97b981e7daf4702875d0f9ff82b62 GIT binary patch literal 6148 zcmeHKyGjH>5Ukb<4ouF>aDKr*SdQ}x?gLbOCF`DeD4K2Lcll{nKMr{`Jw*X{Gr?}mrG>F`mCN&zV#1*Cu!kOIF@z-upUd6KA63P=Gd@U4J<9~#}U zD;yK!)4?HH0OEq-FwUcwAT|#WyTUP%5t=2Hm{hA4!;;Q;tGupoOiVhgnh&d+tvVEo z+j)MAbXb?DQ3^NtxGt&i$@%Obj~XK_}{Gz;%&HfxlMZ3m&`}ivR!s literal 0 HcmV?d00001 diff --git a/.gitignore b/.gitignore index d556509..71a04bf 100644 --- a/.gitignore +++ b/.gitignore @@ -134,3 +134,9 @@ tests/output/ # virtual environments .venv + +# Output of polyconf tutorial + +polyconf_examples/PEI_linear_* +polyconf_examples/PEI_branched_* +polyconf_examples/PMMA_* diff --git a/polyconf/polyconf/Polymer.py b/polyconf/polyconf/Polymer.py index 1eb0454..5c6d925 100644 --- a/polyconf/polyconf/Polymer.py +++ b/polyconf/polyconf/Polymer.py @@ -75,9 +75,10 @@ def copy(self) -> Polymer: def renamer(self, resid, namein, nameout='X'): """ - Change selected atom names to a new name like X1, X2, X3. Intended to - flag dummy atoms for removal. Selected atoms are given a basename, - e.g. 'X' defined by the nameout argument, as well as a number. + Change selected atom names to a new name. + Intended to flag dummy atoms for removal. + Selected atoms are given a basename defined by the nameout argument, + e.g. 'X' , and a number. :param resid: which residue number to select and rename atoms from :type resid: int @@ -101,8 +102,8 @@ def newresid(self): :return: the polymer's current highest resid plus one. :rtype: int """ - n = max(self.polymer.residues.resids) + 1 - return n + n = int(max(self.polymer.residues.resids) + 1) + return (n) def maxresid(self): @@ -112,8 +113,8 @@ def maxresid(self): :return: the polymer's current highest resid :rtype: int """ - n = max(self.polymer.residues.resids) - return n + n = int(max(self.polymer.residues.resids)) + return (n) def extend(self, monomer, n, nn, names, joins, ortho=[1,1,1], @@ -195,8 +196,8 @@ def extend(self, monomer, n, nn, names, joins, ortho=[1,1,1], # cool. # love that for me. - Q = self.polymer.select_atoms(f'resid {n} and name {names["Q"]}').positions[-1] - S = self.polymer.select_atoms(f'resid {n} and name {names["S"]}').positions[-1] + Q = self.polymer.select_atoms(f"resid {n} and name {names['Q']}").positions[-1] + S = self.polymer.select_atoms(f"resid {n} and name {names['S']}").positions[-1] u_ = monomer # monomer = mdict['path'][monomer] u_.atoms.tempfactors=float(beta) diff --git a/polyconf_examples/01_build_PEI_linear.py b/polyconf_examples/01_build_PEI_linear.py deleted file mode 100644 index 3a18aa5..0000000 --- a/polyconf_examples/01_build_PEI_linear.py +++ /dev/null @@ -1,298 +0,0 @@ -#!/usr/bin/env python - -from polyconf.Monomer import Monomer -from polyconf.Polymer import Polymer -from polyconf.PDB import PDB - - -############################################################################### -# # -# Polymer one: Building a polymer with extend() # -# # -############################################################################### - -# this polymer will be 128 monomers long. -# there will be an initial monomer ( PEI_start.pdb ), 126 middle monomers (PEI_monomer.pdb) , and a final monomer ('PEI_end.pdb') - -# first, we initialise a new polymer with the starting monomer - -polymer1=Polymer(Monomer('PEI_start.pdb')) - -# then, we will add the middle monomers -# we will do this by adding one middle monomer, and repeating 126 times - -adds=126 - -for i in range (0,adds): - polymer1.extend( # extend with one monomer, aligned along this step's linearization vector - Monomer('PEI_monomer.pdb'), # extend with this monomer - n=polymer1.maxresid(), # we will allways add onto the existing monomer with the highest resid - nn=polymer1.newresid(), # the incoming monomer needs a new resid - names=dict(Q='CX',P='C1',S='N1',R='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i - joins=[('N1','C1')],# new connection between N1_i and C1_i+1 - ) - - -# now we add the final monomer -# the process is the same, but we are using a different pdb file - -polymer1.extend( - Monomer('PEI_end.pdb'), - n=polymer1.maxresid(), - nn=polymer1.newresid(), - names=dict(Q='CX',P='C1',S='N1',R='NX'), - joins=[('N1','C1')], - ) - -# Now we can save our polymer as a pdb - -# we start by making a PDB object, which will save our file - -Saver = PDB(polymer1) - -Saver.cleanup() # this puts our polymer in the center of the box - -# saver can automatically remove dummy atoms -# here we are saying that every atom with the name CX or NX is a dummy atom -# this works because we set up our files so those dummy atoms had consistent names in every single monomer - -Saver.save(dummyAtoms='CX NX',fname='polymer_01_vanilla-extend') - -# When you examine polymer one , you can see that the resulting strucure is a tightly coiled helix. -# This structure is highly ordered, and the turns of the helix are very close. -# It's not a very likely conformation - -# Lets use some other methods to make different conformations - -############################################################################### -# # -# Polymer two: Building a linear polymer with linear extend() # -# # -############################################################################### - -# this polymer will also be 128 monomers long, and we will use the same input files as before -# however, this time we will extend the polymer by along a vector [1,0,0]. - -polymer2=Polymer(Monomer('PEI_start.pdb')) # initialise the polymer the same way - -adds=126 -aligner=[1,0,0] # linearization vector - -for i in range (0,adds): - - polymer2.extend( - Monomer('PEI_monomer.pdb'), - n=polymer2.maxresid(), - nn=polymer2.newresid(), - names=dict(Q='CX',P='C1',S='N1',R='NX', V1='C1',V2='N1'), # we are aligning each monomer so the vector A1 A2 lies along aligner - linearise=True, # linearize incoming monomer so monomer vector lies along along ortho - ortho=aligner, - joins=[('N1','C1')], - ) - -# now we add the final monomer -# the process is the same, but we are using a different pdb file -polymer2.extend( - Monomer('PEI_end.pdb'), - n=polymer2.maxresid(), - nn=polymer2.newresid(), - names=dict(Q='CX',P='C1',S='N1',R='NX', V1='C1',V2='N1'), # we are aligning each monomer so the vector A1 A2 lies along aligner - linearise=True, # linearize incoming monomer so monomer vector lies along along ortho - ortho=aligner, - joins=[('N1','C1')], - ) - - -# save the polymer as before - -Saver = PDB(polymer2) -Saver.cleanup() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_02_simple-linear-extend') # save, excluding dummy atoms - -# When you examine polymer two , you can see that the resulting structure is extended in a long linear chain. -# The linearization process has introduced some unrealistic bond angles, the angle [ N1_n , C1_n+1 , C2_n+1] is completely linear. -# Additionally, the structure is still highly ordered. - -# These pseudolinear structures can be useful for inspecting your polymer to ensure you have built it correctly. -# They may also be useful starting points for generating valid conformations with shuffle() and dihedral_solver() -# If your simualtion pramaters are correct, the unrealistic angles will be corrected during energy minimzation and equilibration. - - -############################################################################### -# # -# Polymer three: Building a linear polymer with alternating extend() steps # -# # -############################################################################### - -# this polymer will also be 128 monomers long, and we will use the same input files as before -# however, this time we will align to alternating vectors -# even monomers will be aligned to the vector [1,0,0] -# odd monomers will be aligned to the vector [0,1,0] - -# they will continue to switch between the two until they reach the end - -polymer3=Polymer(Monomer('PEI_start.pdb')) # initialise -adds=126 # lay 24 additional monomers -alternator=True # this is a boolean we will use to alternate between our two alignment vectors - -vector1=[1,0,0] # linearization vector -vector2=[0,1,0] # linearization vector - -for i in range (0,adds): - # first, choose which alignment vector to use - - if alternator: - aligner=vector1 - else: - aligner=vector2 - - alternator=not alternator #flip the alternator - - polymer3.extend( - Monomer('PEI_monomer.pdb'), - n=polymer3.maxresid(), - nn=polymer3.newresid(), - names=dict(Q='CX',P='C1',S='N1',R='NX', V1='C1',V2='N1'), - linearise=True, - ortho=aligner, - joins=[('N1','C1')], - ) - -# then add the final monomer - - if alternator: - aligner=vector1 - else: - aligner=vector2 - -polymer3.extend( - Monomer('PEI_end.pdb'), - n=polymer3.maxresid(), - nn=polymer3.newresid(), - names=dict(Q='CX',P='C1',S='N1',R='NX', V1='C1',V2='N1'), - linearise=True, - ortho=aligner, - joins=[('N1','C1')], - ) - -# save this polymer - -Saver = PDB(polymer3) -Saver.cleanup() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_03_alternating-linear-extend') # save, excluding dummy atoms - -# When you examine polymer three , you can see that the resulting structure is again extended in a long linear chain. -# The conformation of this chain is different, and the bond angles are less unrealistic. Hopefully, energy minimzation would correct this -# however there are still some problems. -# This conformation has ammonium groupos positioned close to each other, and may be high energy. -# Additionally, the structure is still highly ordered. - -# for our next polymer, we will address the problem of highly ordered starting structures. - -############################################################################### -# # -# Polymer four: Building a reasonable conformation with dihedral_solver() # -# # -############################################################################### - -# in this example, we will make a copy of polymer 1, then use dihedral_solver() to generate a more reasonable starting conformation - -polymer4=polymer1.copy() # make a copy of our first polymer, which was a tightly packed helix - -# we need to choose a set of deihdrals to shuffle -# I am going to choose every dihedral centered on the bond between atoms C2 and N1 - -# we use gen pairlist to generate a list of all of these dihedrals -# note that we only specify the two atoms in the center of the dihedral - -CN_dihedrals=polymer4.gen_pairlist(J='C2',K='N1',first_resid=1,last_resid=128,mult=3) - -# next, we use dihedral solver to try to solve this conformation -# solver will rotate each dihedral to reach a conformation where no atoms on one side of the dihedral are within a specified distance of any atom on the other side of the dihedral - -# by default, this only looks at atoms from residues with resid up to the current monomer -# we are using a cutoff of 0.7 angrstroms -# this means that when you attempt to solve dihedral in monomer 5, it splits monomer 5 into two halves, separated by the dihedral -# and it tries to find a conformation where two things are true: -# no atom on one side of the dihedral is within 0.7 A of any atom on the other side of the dihedral -# no atom in monomer 5 is within 0.7 A of any atom on monomers with resid between zero and 5 - -# it does not look beyond residue 5. -# the reason for this is that we want to solve the polymer starting from one end, and working along the chain - -cutoff = 0.7 # we the generated conformation to have no two atoms within 0.7 A of each other -# this sort of distance is typically far enough to remove very high energy from overlapping atoms so that the structure can be further refined with energy minimzation - -polymer4.dihedral_solver(CN_dihedrals,dummies='CX,NX',cutoff=cutoff) #TODO why have I got a comma in there -# this will iterate through every bond in C2-N1 bond, rotating them so that atoms one either side of the bond are always more than 0.7 A apart. -# when it checks interatom distances, it ignores dummy atoms -# we have to specify the names of dummy atoms - -Saver = PDB(polymer4) - -Saver.cleanup() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_04_solved') # save, excluding dummy atoms - -# When you examine polymer four, you can see that the resulting structure is much more reasonable than polymer 1. -# The polymer is arranged in a twisting helix. The loops are much less tightly packed, and the distances between adjacent loops are much more reasonable. -# Furthermore, because we did not use a linear extend, the bonds and angles in this structure are consistent with those in the original monomers -# However, dihedal_solver() is algorithmic, and so it has still resulted in a highly ordered structure. -# This polymer may be reasonable for illustrative purposes, but we may wish to generate an ensemble of less ordered structure for production simulations. - -############################################################################### -# # -# Polymer five: Generating a randomised starting conformation shuffler() # -# # -############################################################################### - -# in this example, we will make a copy of polymer 1, then use shuffler() to generate a random conformation -# we will then use dihedral_solver() to convert the random conformation into something reasonable - -#first, make a copy of polymer one - -polymer5=polymer1.copy() - - -import random - -random.seed(8) # by setting the random seed we can ensure the random conformations are replicable. - -# we are going to generate a conformation in a two step process -# first, we will are going to randomly rotate the dihedrals centered on the bond between atoms N1_n and C1_n+1 -# then, we are going to solve the polymer by rotating each dihedral centered on a C1-C2 bond - -# in principle, you can shuffle as many dihedrals as you like -# however, the process of randomly rotating dihedrals does not consider interatomic distances -# in my experience, randomising more dihedrals results in very tight compacted conformations with may overlapping atoms -# it is best to choose only one ot two bonds to shuffle -# if you are not having success shuffling one bond, it may be useful to try shuffling a different bond instead -# or shuffling a second time to find a more amenable starting conformation - -NC_dihedrals=polymer5.gen_pairlist(J='N1',K='C1',first_resid=1,last_resid=127,mult=6,same_res=False) -# these dihedrals are centered on the bond that links two monomers, so we need to set same_res to False -# we have also set the multiplicity of these dihedrals to 6, to increase the search space and improve the chance of finding a valid conformation - -alkane_dhedrals=polymer5.gen_pairlist(J='C1',K='C2',first_resid=1,last_resid=128,mult=3) - -polymer5.shuffler(NC_dihedrals,) # this generates a random conformation - -# lets look at the intermediate conformation we got by random shuffling - -Saver = PDB(polymer5) -Saver.cleanup() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_05_shuffled') # save, excluding dummy atoms - -# now solvce the conformation -polymer5.dihedral_solver(alkane_dhedrals,dummies='CX NX',cutoff=1) # this converts the shuffled conformation into one without overlapping atoms - -Saver = PDB(polymer5) -Saver.cleanup() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_05_shuffled_then_solved') # save, excluding dummy atoms - -# Polymer five has two output files, shuffled, and shuffled_then_solved. -# When you examine polymer five, you can see that the two conformations are extended random coils, and is quite different to the conformations of polymers one to four. - -# The shuffled conformation has some steric clashes. For example, monomers 47 and 54 are overlapping. -# The shuffled_then_solved conformation has been adjusted to remove these clashes. - -# TODO make branched tutorial \ No newline at end of file diff --git a/polyconf_examples/01a_build_PEI_with_extend.py b/polyconf_examples/01a_build_PEI_with_extend.py new file mode 100644 index 0000000..885936e --- /dev/null +++ b/polyconf_examples/01a_build_PEI_with_extend.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + + +############################################################################### +# # +# Polymer one: Building a polymer with extend() # +# # +############################################################################### + +# this polymer will be 128 monomers long. +# there will be an initial monomer ( PEI_start.pdb ), 126 middle monomers (PEI_monomer.pdb) , and a final monomer ('PEI_end.pdb') + +# first, we initialise a new polymer with the starting monomer + +polymer1=Polymer(Monomer('PEI_start.pdb')) + +# then, we will extend our polymer by adding new monomers +# we will repeat this 126 times + +adds=126 + +for _ in range (0,adds): # repeat 126 times + + n=polymer1.maxresid() + new_n=polymer1.newresid() + + polymer1.extend( # extend with one monomer + Monomer('PEI_monomer.pdb'), # we are adding a new monomer using the coordinates from the file 'PEI_monomer.pdb' + n=n, + nn=new_n, + names=dict(P='C1',Q='CX',R='NX',S='N1'), # this defines the mapping used to connect our monomers. The new polymer will be created so the PR bond in the new monomer is colinear with the QS bond in the existing monomer + joins=[('N1','C1')], # new connection between N1_n and C1_nn + ) + + # after adding each monomer, we will flag dummy atoms + # in this case, the dummy atoms are atoms CX in monomer n, and atom NX in monomer nn + + polymer1.renamer(n,'CX') # convert atom CX in monomer n to a dummy atom + polymer1.renamer(new_n,'NX') # convert atom NX in monomer new_n to a dummy atom + + + +# we now have a chain of 127 monomers +# now we add the final monomer +# the process is the same, but we are using a different pdb file + +n = polymer1.maxresid() +new_n = polymer1.newresid() + +polymer1.extend( + Monomer('PEI_end.pdb'), + n=polymer1.maxresid(), + nn=polymer1.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + +polymer1.renamer(n,'CX') +polymer1.renamer(new_n,'NX') + +# Now we can save our polymer as a pdb + +# we start by making a PDB object, which will save our file +print('saving') + +Saver = PDB(polymer1) + +Saver.cleanup() # this puts our polymer in the center of the box +Saver.save(fname='PEI_linear_01_vanilla-extend') + +# When you examine polymer one , you can see that the resulting strucure is a tightly coiled helix. +# This structure is highly ordered, and the turns of the helix are very close. +# It's not a very likely conformation + +# In the next tutorial, we will use some other methods to make different conformations + diff --git a/polyconf_examples/01b_build_PEI_with_linear_extend.py b/polyconf_examples/01b_build_PEI_with_linear_extend.py new file mode 100644 index 0000000..a0d0a24 --- /dev/null +++ b/polyconf_examples/01b_build_PEI_with_linear_extend.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + + +############################################################################### +# # +# Polymer two: Building a linear polymer with linear extend() # +# # +############################################################################### + +# this polymer will also be 128 monomers long, and we will use the same input files as before +# however, this time we will extend the polymer by along a vector [1,0,0]. +# this will result in a pseudolinear conformation, which should greatly reduce the number of overlapping atoms + +polymer2=Polymer(Monomer('PEI_start.pdb')) # initialise the polymer the same way + +adds=126 +aligner=[1,0,0] # linearization vector + +for _ in range (0,adds): + + polymer2.extend( + Monomer('PEI_monomer.pdb'), + n=polymer2.maxresid(), + nn=polymer2.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1', V1='C1',V2='N1'), # we are aligning each monomer so the vector A1 A2 lies along aligner + linearise=True, # linearize incoming monomer so monomer vector lies along a specified vector, named ortho + ortho=aligner, + joins=[('N1','C1')], + ) + +# now we add the final monomer +# the process is the same, but we are using a different pdb file + +polymer2.extend( + Monomer('PEI_end.pdb'), + n=polymer2.maxresid(), + nn=polymer2.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1', V1='C1',V2='N1'), # we are aligning each monomer so the vector A1 A2 lies along aligner + linearise=True, # linearize incoming monomer so monomer vector lies along along ortho + ortho=aligner, + joins=[('N1','C1')], + ) + + +# you may have noticed some differences + +# first, we did not define n and new_n separately +# we just called maxresid() and newresid() during extend() + +# we also didn't use renamer to flag dummy atoms +# this is because our monomer coordinate files were designed so every monomer has the same atoms names, +# so that our dummy atoms would be be every atom named NX or CX + +# when we save our polymer, we can explicitly name our dummy atoms + +Saver = PDB(polymer2) +Saver.cleanup() +Saver.save(dummies='CX NX',fname='PEI_linear_02_simple-linear-extend') # save, excluding dummy atoms + +# When you examine polymer two , you can see that the resulting structure is extended in a long linear chain. +# The linearization process has introduced some unrealistic bond angles, the angle [ N1_n , C1_n+1 , C2_n+1] is completely linear. +# Additionally, the structure is still highly ordered. + +# These pseudolinear structures can be useful for inspecting your polymer to ensure you have built it correctly. +# They may also be useful starting points for generating valid conformations with shuffle() and dihedral_solver() +# If your simulation parameters are correct, the unrealistic angles will be corrected during energy minimzation and equilibration. + +############################################################################### +# # +# Polymer three: Building a linear polymer with alternating extend() steps # +# # +############################################################################### + +# this polymer will also be 128 monomers long, and we will use the same input files as before +# however, this time we will align to alternating vectors +# even monomers will be aligned to the vector [1,0,0] +# odd monomers will be aligned to the vector [0,1,0] + +# they will continue to switch between the two until they reach the end + +polymer3=Polymer(Monomer('PEI_start.pdb')) # initialise +adds=126 # lay 24 additional monomers +alternator=True # this is a boolean we will use to alternate between our two alignment vectors + +vector1=[1,0,0] # linearization vector +vector2=[0,1,0] # linearization vector + +for i in range (0,adds): + # first, choose which alignment vector to use + + if alternator: + aligner=vector1 + else: + aligner=vector2 + + alternator=not alternator #flip the alternator + + polymer3.extend( + Monomer('PEI_monomer.pdb'), + n=polymer3.maxresid(), + nn=polymer3.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1', V1='C1',V2='N1'), + linearise=True, + ortho=aligner, + joins=[('N1','C1')], + ) + +# then add the final monomer + + if alternator: + aligner=vector1 + else: + aligner=vector2 + +polymer3.extend( + Monomer('PEI_end.pdb'), + n=polymer3.maxresid(), + nn=polymer3.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1', V1='C1',V2='N1'), + linearise=True, + ortho=aligner, + joins=[('N1','C1')], + ) + +# save this polymer + +Saver = PDB(polymer3) +Saver.cleanup() # center in box +Saver.save(dummies='CX NX',fname='PEI_linear_03_alternating-linear-extend') # save, excluding dummy atoms + +# When you examine polymer three , you can see that the resulting structure is again extended in a long linear chain. +# The conformation of this chain is different, and the bond angles are less unrealistic. If your simulation parameters are correct, energy minimzation would ajust the bond angles +# however there are still some problems. +# This conformation has ammonium groups positioned close to each other, despite the repulsive electrostatic interactions betweeen them. +# Additionally, the structure is still highly ordered. + +# For our next polymer, we will address the problem of highly ordered starting structures. diff --git a/polyconf_examples/01c_build_PEI_then_solve.py b/polyconf_examples/01c_build_PEI_then_solve.py new file mode 100644 index 0000000..6fd2875 --- /dev/null +++ b/polyconf_examples/01c_build_PEI_then_solve.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + + +############################################################################### +# # +# Polymer four: Building a reasonable conformation with dihedral_solver() # +# # +############################################################################### + + + +# in this example, we will make a copy through the same process as in polymer 1, then use dihedral_solver() to generate a more reasonable starting conformation + +polymer4=Polymer(Monomer('PEI_start.pdb')) + +adds=126 + +for i in range (0,adds): + polymer4.extend( + Monomer('PEI_monomer.pdb'), + n=polymer4.maxresid(), + nn=polymer4.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + +polymer4.extend( + Monomer('PEI_end.pdb'), + n=polymer4.maxresid(), + nn=polymer4.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + +# we need to choose a set of dihedrals to shuffle + +# I am going to choose every dihedral centered on the torsion between atoms N1 and C1. This is the connection between monomers. + +# we use gen pairlist to generate a list of all of these dihedrals +# note that we only specify the two atoms in the center of the dihedral + +dihedrals=polymer4.gen_pairlist(J='N1',K='C1',first_resid=1,last_resid=127,mult=3,same_res=False) + +# next, we use dihedral solver to try to solve this conformation +# solver will rotate each dihedral to attempt to find a conformation where no atoms on one side of the dihedral are within a specified distance of any atom on the other side of the dihedral + +cutoff = 0.9 + +# we want the generated conformation to have no two atoms within 0.9 A of each other + +# this means that when you attempt to solve dihedral between monomers 4 and 5, it splits the polymer into two halves, separated by the dihedral +# then rotates the torsion tries to find a conformation where no atom on one side of the dihedral is within 0.9 A of any atom on the other side of the dihedral, considering only atoms in monomers 1 to 5 + +# it does not look beyond residue 5. +# the reason for this is that we want to solve the polymer starting from one end, and working along the chain + +# In my experience, this sort of distance is typically sufficient to produce a structure suitable for energy minimization, followed by equilibration. + +polymer4.dihedral_solver(dihedrals,dummies='CX NX',cutoff=cutoff) +# this will iterate through every N1-C1 torsion, rotating them so that atoms one either side of the torsion are always more than 0.7 A apart. +# when it checks interatom distances, it ignores dummy atoms provided we have specified the names of dummy atoms + +Saver = PDB(polymer4) + +Saver.cleanup() # center in box +Saver.save(dummies='CX NX',fname='PEI_linear_04_solved') # save, excluding dummy atoms + +# When you examine polymer four, you can see that the resulting structure is much more reasonable than polymer 1. +# The polymer is arranged in a twisting helix. The loops are much less tightly packed, and the distances between adjacent loops are much more reasonable. +# Furthermore, because we did not use a linear extend, the bonds and angles in this structure are consistent with those in the original monomers +# However, dihedal_solver() is algorithmic, and so it has still resulted in a highly ordered structure. +# This polymer may be reasonable for illustrative purposes, but we may wish to generate an ensemble of less ordered structure for production simulations. diff --git a/polyconf_examples/01d_build_PEI_conformation.py b/polyconf_examples/01d_build_PEI_conformation.py new file mode 100644 index 0000000..7a55450 --- /dev/null +++ b/polyconf_examples/01d_build_PEI_conformation.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + + +############################################################################### +# # +# Polymer five: Generating a randomised starting conformation shuffler() # +# # +############################################################################### + +# in this example, we will create a PEI polymer as in example 01a, then use shuffler() to generate a random conformation +# we will then use dihedral_solver() to convert the random conformation into something reasonable + +#first, make a copy of polymer one + +polymer5=Polymer(Monomer('PEI_start.pdb')) + +# then, we will add the middle monomers +# we will do this by adding one middle monomer, and repeating 126 times + +adds=126 + +for i in range (0,adds): + polymer5.extend( # extend with one monomer, aligned along this step's linearization vector + Monomer('PEI_monomer.pdb'), # extend with this monomer + n=polymer5.maxresid(), # we will allways add onto the existing monomer with the highest resid + nn=polymer5.newresid(), # the incoming monomer needs a new resid + names=dict(P='C1',Q='CX',R='NX',S='N1'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + joins=[('N1','C1')],# new connection between N1_i and C1_i+1 + ) + +# now we add the final monomer +# the process is the same, but we are using a different pdb file + +polymer5.extend( + Monomer('PEI_end.pdb'), + n=polymer5.maxresid(), + nn=polymer5.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + + + +import random + +random.seed(8) # by setting the random seed we can ensure the random conformations are replicable. + +# we are going to generate a conformation in a two step process +# first, we will are going to randomly rotate the dihedrals centered on the bond between atoms N1_n and C1_n+1 +# then, we are going to solve the polymer by rotating each dihedral centered on a C1-C2 bond + +# in principle, you can shuffle as many dihedrals as you like +# however, the process of randomly rotating dihedrals does not consider interatomic distances +# in my experience, randomising more dihedrals results in very tight compacted conformations with may overlapping atoms +# it is best to choose only one ot two bonds to shuffle +# if you are not having success shuffling one bond, it may be useful to try shuffling a different bond instead +# or shuffling a second time to find a more amenable starting conformation + +NC_dihedrals=polymer5.gen_pairlist(J='N1',K='C1',first_resid=1,last_resid=127,mult=3,same_res=False) +# these dihedrals are centered on the bond that links two monomers, so we need to set same_res to False + +alkane_dihedrals=polymer5.gen_pairlist(J='C1',K='C2',first_resid=1,last_resid=128,mult=3) + +polymer5.shuffler(NC_dihedrals,) # this generates a random conformation + +# lets look at the intermediate conformation we got by random shuffling + +Saver = PDB(polymer5) +Saver.cleanup() # center in box +Saver.save(dummies='CX NX',fname='PEI_linear_05_shuffled') # save, excluding dummy atoms + +# now solve the conformation +polymer5.dihedral_solver(alkane_dihedrals,dummies='CX NX',cutoff=0.9) # this converts the shuffled conformation into one without overlapping atoms, by rotating the C1-C2 torsions + +Saver = PDB(polymer5) +Saver.cleanup() # center in box +Saver.save(dummies='CX NX',fname='PEI_linear_05_shuffled_then_solved') # save, excluding dummy atoms + +# Polymer five has two output files, shuffled, and shuffled_then_solved. +# When you examine polymer five, you can see that the two conformations are extended random coils, and is quite different to the conformations of polymers one to four. + +# The shuffled conformation has some steric clashes. For example, monomers 47 and 54 are overlapping. +# The shuffled_then_solved conformation has been adjusted to remove these clashes. + diff --git a/polyconf_examples/01e_build_PEI_conformation_ensemble.py b/polyconf_examples/01e_build_PEI_conformation_ensemble.py new file mode 100644 index 0000000..080fd06 --- /dev/null +++ b/polyconf_examples/01e_build_PEI_conformation_ensemble.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + +import random + +random.seed(12) # by setting the random seed we can ensure the random conformations are replicable. + +############################################################################### +# # +# Polymer six: Generating an ensemble of conformations # +# # +############################################################################### + +# In this example, we will follow the process used in polymer 5 to generate a conformation, +# but we will repeat it five times to generate an ensemble of starting conformations. + +# Using different starting conformations reduces sampling bias arising from the starting conformation. +# This is very useful when you wish to study polymer dynamics + +# First, make a PEI polymer through the same process as polymer 1 + +polymer6=Polymer(Monomer('PEI_start.pdb')) + +adds=126 + +for i in range (0,adds): + polymer6.extend( # extend with one monomer, aligned along this step's linearization vector + Monomer('PEI_monomer.pdb'), # extend with this monomer + n=polymer6.maxresid(), # we will allways add onto the existing monomer with the highest resid + nn=polymer6.newresid(), # the incoming monomer needs a new resid + names=dict(P='C1',Q='CX',R='NX',S='N1'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + joins=[('N1','C1')],# new connection between N1_i and C1_i+1 + ) +polymer6.extend( + Monomer('PEI_end.pdb'), + n=polymer6.maxresid(), + nn=polymer6.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + + +# Generate lists of dihedrals as for polymer 5 + +NC_dihedrals=polymer6.gen_pairlist(J='N1',K='C1',first_resid=1,last_resid=127,mult=3,same_res=False) +alkane_dihedrals=polymer6.gen_pairlist(J='C1',K='C2',first_resid=1,last_resid=128,mult=3) + +# Then, we will generate five starting conformations + +# For each conformation: +# start by making a copy of our initial conformation +# randomise the conformation by shuffling the NC dihedrals +# solve the conformation by rotating the alkane dihedrals + +for conf in range(1,6): + + print(f'generating conformation {conf}') + newconf=polymer6.copy() # make a duplicate of the original polymer + + newconf.shuffler(NC_dihedrals) # as before, randomise the NC dihedrals to generate a random conformation + newconf.dihedral_solver(alkane_dihedrals,dummies='CX NX',cutoff=0.9) # solve by rotating the alkane dihedrals + + Saver = PDB(newconf) + Saver.cleanup() + Saver.save(dummies='CX NX',fname=f'PEI_linear_06_conformation_{conf}') + +# This should produce an ensemble of five starting conformations diff --git a/polyconf_examples/02a_build_PMMA_isotactic.py b/polyconf_examples/02a_build_PMMA_isotactic.py index a5598ae..f78b634 100644 --- a/polyconf_examples/02a_build_PMMA_isotactic.py +++ b/polyconf_examples/02a_build_PMMA_isotactic.py @@ -25,31 +25,33 @@ Monomer('MMAD_bonds.pdb'), # extend with this monomer n=PMMA_isotactic.maxresid(), # extend existing residue i nn=PMMA_isotactic.newresid(), # incoming monomer will have resid i+1 - names=dict(Q='CA',P='CMA',S='C',R='CN',), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i - joins=[('C','CA')],# new connection between N1_i and C1_i+1 - #linearise=True, - #ortho=[1,0,0] + names=dict(P='CMA',Q='CA',R='CN',S='C'), # CA_i+1 fit to CMA_i, then rotate so C_i+1 fit to CN_i + joins=[('C','CA')],# new connection between C_i and CA_i+1 ) -alkanes=PMMA_isotactic.gen_pairlist(J='C',K='CA',first_resid=1,same_res=False,last_resid=49,mult=3) -sidechains=PMMA_isotactic.gen_pairlist(J='CA',K='CB',first_resid=1,same_res=True,last_resid=50,mult=6) - Saver = PDB(PMMA_isotactic) Saver.cleanup() # center in box -Saver.save(dummies=dummies,fname='PMMA_isotactic_vanilla-extend') # save, excluding dummies atoms +Saver.save(dummies=dummies,fname='PMMA_isotactic_vanilla-extend') # save, excluding dummy atoms + +# next, generate conformations by adjusting the dihedrals -PMMA_isotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=0.8) # this converts the shuffled conformation into one without overlapping atoms -PMMA_isotactic.dihedral_solver(sidechains,dummies=dummies,cutoff=0.8) # this converts the shuffled conformation into one without overlapping atoms +alkanes=PMMA_isotactic.gen_pairlist(J='C',K='CA',first_resid=1,same_res=False,last_resid=49,mult=3) +sidechains=PMMA_isotactic.gen_pairlist(J='CA',K='CB',first_resid=1,same_res=True,last_resid=50,mult=6) + +PMMA_isotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=0.8) # this converts the shuffled conformation into one without overlapping atoms by rotating the alkanes Saver = PDB(PMMA_isotactic) Saver.cleanup() # center in box -Saver.save(dummies=dummies,fname='PMMA_isotactic_vanilla-solved') # save, excluding dummies atoms +Saver.save(dummies=dummies,fname='PMMA_isotactic_vanilla-solved') -PMMA_isotactic.shuffler(alkanes) +PMMA_isotactic.shuffler(alkanes) # randomise the rotation of C-CA torsions PMMA_isotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms +PMMA_isotactic.shuffler(sidechains) # randomise the position of the sidechains +PMMA_isotactic.dihedral_solver(sidechains,dummies=dummies,cutoff=0.8) + Saver = PDB(PMMA_isotactic) Saver.cleanup() # center in box -Saver.save(dummies=dummies,fname='PMMA_isotactic_shuffled_and_solved') # save, excluding dummies atoms +Saver.save(dummies=dummies,fname='PMMA_isotactic_shuffled_and_solved') diff --git a/polyconf_examples/03_build_branched_PEI.py b/polyconf_examples/03_build_branched_PEI.py new file mode 100644 index 0000000..ae02fe6 --- /dev/null +++ b/polyconf_examples/03_build_branched_PEI.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB + +import random + +random.seed(10) # by setting the random seed we can ensure the random conformations are replicable. + +############################################################################### +# # +# Branched PEI # +# # +############################################################################### + +# In this example, we will follow the process used to generate a linear PEI, 128 monomers long + +# but then we will add branches every 5 monomers + +PEI_branched=Polymer(Monomer('PEI_start.pdb')) + +adds=126 + +for i in range (0,adds): + PEI_branched.extend( + Monomer('PEI_monomer.pdb'), + n=PEI_branched.maxresid(), + nn=PEI_branched.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + +# add the final monomer + +PEI_branched.extend( + Monomer('PEI_end.pdb'), + n=PEI_branched.maxresid(), + nn=PEI_branched.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + ) + + +# now we have a linear chain, we can add our branches +# for this example each branch will be three units long + +branch_bonds=[] # make a list to keep track of all of our branches +branch=1 + +for n in range(5,128,5): + + branch_resid=PEI_branched.newresid() + branch_bonds+=[(n,branch_resid)] + # add the first monomer of the branch + PEI_branched.extend( + Monomer('PEI_monomer.pdb'), + n=n, + nn=branch_resid, + names=dict(P='C1',Q='H1',R='NX',S='N1'), + joins=[('N1','C1')], + beta=branch # we are using beta factors to label our branches + ) + + PEI_branched.renamer(n,'H1') + + # add the second monomer of the branch + + PEI_branched.extend( + Monomer('PEI_monomer.pdb'), + n=PEI_branched.maxresid(), + nn=PEI_branched.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + beta=branch + + ) + + # add the third monomer of the branch + + PEI_branched.extend( + Monomer('PEI_end.pdb'), + n=PEI_branched.maxresid(), + nn=PEI_branched.newresid(), + names=dict(P='C1',Q='CX',R='NX',S='N1'), + joins=[('N1','C1')], + beta=branch + + ) + + branch+=1 # change the branch number for the next branch + +# Generate lists of dihedrals + +mainchain_NC= PEI_branched.gen_pairlist(J='N1',K='C1',first_resid=1,last_resid=127,mult=3,same_res=False) +mainchain_alkanes= PEI_branched.gen_pairlist(J='C1',K='C2',first_resid=1,last_resid=128,mult=3) + +# we are going to manually generate a dihedral list for the connections between the main chain and the branches +branch_dihedrals=[{'J': 'N1', 'J_resid': i, 'K': 'C1', 'K_resid': j, 'mult': 3} for i,j in branch_bonds] + +branch_alkanes=PEI_branched.gen_pairlist(J='C1',K='C2',first_resid=129,last_resid=PEI_branched.maxresid(),mult=3) + +for conf in range(1,6): + + print(f'generating conformation {conf}') + newconf=PEI_branched.copy() # make a duplicate of the original polymer + + print(f'shuffling backbone') + + newconf.shuffler(mainchain_NC) # as before, randomise the NC dihedrals to generate a random conformation + + print(f'solving backbone') + + newconf.dihedral_solver(mainchain_alkanes,dummies='CX NX X*',cutoff=0.9) # solve by rotating the alkane dihedrals + + print(f'shuffling branchpoints') + + newconf.shuffler(branch_dihedrals) + + print(f'solving branchpoints') + + newconf.dihedral_solver(branch_dihedrals,dummies='CX NX X*',cutoff=0.9,backwards_only=False) + + print(f'solving branches') + + newconf.dihedral_solver(branch_alkanes,dummies='CX NX X*',cutoff=0.9,backwards_only=False) + + + + Saver = PDB(newconf) + Saver.cleanup() + Saver.save(dummies='CX NX X1',fname=f'PEI_branched_conformation_{conf}') + +# This should produce an ensemble of five starting conformations + +# when you visualize these, try colouring by beta factor +# you can see each branch chain is labelled with different beta factor +# this can be useful for manipulating and keeping track of branches \ No newline at end of file diff --git a/polyconf_examples/README.md b/polyconf_examples/README.md index ef8e161..60c0f3b 100644 --- a/polyconf_examples/README.md +++ b/polyconf_examples/README.md @@ -21,7 +21,7 @@ This tutorial is all contained in a single script. There are detailed comments ## Tutorial 02: Strategies for building conformations of polymers with different types of tacticity, or copolymers with different monomer distributions -The second tutorial shows strategies for making polymers with different types of [tacticity](https://en.wikipedia.org/wiki/Tacticity). Using poly methyl methacrylate, this lesson showcases a strategy for making an isotactic polymer, a syndiotactic polymer, and an atactic polymer. These strategies use two monomers, `MMAD_bonds.pdb` and `MMAL_bonds.pdb`. +The second tutorial shows strategies for making polymers with different types of [tacticity](https://en.wikipedia.org/wiki/Tacticity). Using [poly(methyl methacrylate)](https://en.wikipedia.org/wiki/Poly(methyl_methacrylate)) as an example, this lesson showcases a strategy for making an isotactic polymer, a syndiotactic polymer, and an atactic polymer. These strategies use two monomers, `MMAD_bonds.pdb` and `MMAL_bonds.pdb`. In these examples, the two monomers are enantiomers. However, these strategies are also applicable to the creation of [copolymers](https://en.wikipedia.org/wiki/Copolymer), where the monomers are chemically distinct. @@ -29,19 +29,19 @@ This tutorial is split into four scripts. Each script includes detailed comment ### Building an isotactic PMMA polymer -`02a_build_PMMA_isotactic.py` shows a strategy for building an isotactic PMMA polymer, in which all monomers have the same tacticity. The method is very similar to the approach used in tutorial 01. +`02a_build_PMMA_isotactic.py` shows a strategy for building an [isotactic](https://en.wikipedia.org/wiki/Tacticity#Isotactic_polymers) PMMA polymer, in which all monomers have the same tacticity. The method is very similar to the approach used in tutorial 01. ### Building a sydniotactic PMMA polymer -`02b_build_PMMA_syndiotactic.py` shows a strategy for building a syndiotactic PMMA polymer, where monomers alternate between two different enantiomers. +`02b_build_PMMA_syndiotactic.py` shows a strategy for building a [syndiotactic](https://en.wikipedia.org/wiki/Tacticity#Syndiotactic_polymers) PMMA polymer, where monomers alternate between two different enantiomers. -This strategy could be adapted to make an [alternating copolymer](https://en.wikipedia.org/wiki/Copolymer). +This strategy could be adapted to make an [alternating copolymer](https://en.wikipedia.org/wiki/Copolymer#Alternating_copolymers). ### Failing to build an atactic PMMA polymer -`02c_build_PMMA_atactic.py` shows how to build an atactic PMMA polymer, with tacticity of the monomers distributed randomly along the chain. +`02c_build_PMMA_atactic.py` shows how to build an [atactic](https://en.wikipedia.org/wiki/Tacticity#Atactic_polymers) PMMA polymer, with tacticity of the monomers distributed randomly along the chain. -This approach could also be adapted to generate a [statistical copolymer](https://en.wikipedia.org/wiki/Copolymer). +This approach could also be adapted to generate a [statistical copolymer](https://en.wikipedia.org/wiki/Copolymer#Statistical_copolymers). The script for this part of the tutorial **will fail!** This is intentional. We have chosen a specific random seed that will reliably fail to generate a valid starting conformations, in order to showcase some of the challenges in generating polymers where monomers are arranged in a random or stochastic order. @@ -53,7 +53,7 @@ There is a discussion at the end of this script, which contain several suggestio This example is one solution to resolve the difficulties showcased in `02c_build_PMMA_atactic.py`. This strategy uses `extend()` with `linearise=True` to produce a pseudo-linear conformation. It then applies `shuffler()` and `dihedral_solver()` to generate an ensemble of five starting conformations for use in replicate molecular dynamics simulations. -This type of approach could also be used to generate a [statistical copolymer](https://en.wikipedia.org/wiki/Copolymer). +Again, this type of approach could also be used to generate a [statistical copolymer](https://en.wikipedia.org/wiki/Copolymer#Statistical_copolymers). ## Tutorial 03: Building a branched polyethyleneimine polymer conformation using PolyConf From ba0e8776f232ded7d29c740a5dfd881b52eefe60 Mon Sep 17 00:00:00 2001 From: recombinatrix Date: Mon, 17 Feb 2025 22:51:53 +1100 Subject: [PATCH 2/2] extended polyconf examples --- polyconf_examples/README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/polyconf_examples/README.md b/polyconf_examples/README.md index 60c0f3b..25ca1c4 100644 --- a/polyconf_examples/README.md +++ b/polyconf_examples/README.md @@ -8,6 +8,8 @@ These tutorials are not intended as a primer on polymer design or stereochemistr ## Tutorial 01: Building a linear polyethyleneimine polymer conformation using PolyConf +**TODO:** extend based on rewrite + The first tutorial is an example of how to make a linear 128 unit polyethyleneimine chain using *PolyConf*. It is contained in the file `01_build_PEI_linear.py` This tutorial showcases several approaches, including: @@ -57,8 +59,8 @@ Again, this type of approach could also be used to generate a [statistical copol ## Tutorial 03: Building a branched polyethyleneimine polymer conformation using PolyConf -To be written +**TODO:** Ada to write description -## Tutorial 04: Building a hyperbranched polyethyleneimine polymer conformation using PolyConf +## Tutorial 04: Building a polyethyleneimine dendrimer conformation using PolyConf To be written