From 812d1eb07780927e1709ed0dc834a58716f45558 Mon Sep 17 00:00:00 2001 From: recombinatrix Date: Mon, 17 Feb 2025 15:18:04 +1100 Subject: [PATCH 1/3] reworked atom names to align with manuscript --- polyconf/polyconf/PDB.py | 22 +- polyconf/polyconf/Polymer.py | 222 +++++++------- polyconf_examples/build_PEI_linear.py | 290 +++++++++++------- polyconf_examples/build_PMMA_01-isotactic.py | 55 ++++ .../build_PMMA_02-syndiotactic.py | 64 ++++ polyconf_examples/build_PMMA_03-atactic.py | 122 ++++++++ .../build_PMMA_04-atactic-linearized.py | 106 +++++++ 7 files changed, 654 insertions(+), 227 deletions(-) create mode 100644 polyconf_examples/build_PMMA_01-isotactic.py create mode 100644 polyconf_examples/build_PMMA_02-syndiotactic.py create mode 100644 polyconf_examples/build_PMMA_03-atactic.py create mode 100644 polyconf_examples/build_PMMA_04-atactic-linearized.py diff --git a/polyconf/polyconf/PDB.py b/polyconf/polyconf/PDB.py index 18a7e99..ea3b099 100644 --- a/polyconf/polyconf/PDB.py +++ b/polyconf/polyconf/PDB.py @@ -33,15 +33,15 @@ def _write(self, selection, name): atoms = self.polymer.select_atoms(selection) atoms.write(name) - def save(self, dummyAtoms="X*", fname="polymer", selectionString = None, gmx = False): + def save(self, dummies="X*", fname="polymer", selectionString = None, gmx = False): """ Save polymer as a PDB or GROMACS file with dummy atoms excluded. Optionally, select a subset of the polymer to save. - :param dummyAtoms: names of all the dummy atoms, for use in the + :param dummies: names of all the dummy atoms, for use in the selection string (e.g. 'CN CMA CP CQ' to exclude these four dummy atom types), defaults to "X*" - :type dummyAtoms: str, optional + :type dummies: str, optional :param fname: name of the output file, defaults to "polymer" :type fname: str, optional :param selectionString: MDAnalysis atom selection string, for more @@ -55,14 +55,26 @@ def save(self, dummyAtoms="X*", fname="polymer", selectionString = None, gmx = F """ if selectionString: if gmx: +<<<<<<< HEAD + self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.gro") + else: + self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.pdb") +======= +<<<<<<< HEAD self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.gro") else: self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.pdb") +======= + self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.gro") + else: + self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.pdb") +>>>>>>> d62405c (rewrote all polyconf functions to use updated atom names) +>>>>>>> 9e689d2 (Testing rebase against main worked) else: if gmx: - self._write(f"not name {dummyAtoms}", f"{fname}.gro") + self._write(f"not name {dummies}", f"{fname}.gro") else: - self._write(f"not name {dummyAtoms}", f"{fname}.pdb") + self._write(f"not name {dummies}", f"{fname}.pdb") def crudesave(self,fname="polymer_crude"): """ diff --git a/polyconf/polyconf/Polymer.py b/polyconf/polyconf/Polymer.py index dd9fbd7..1f07340 100644 --- a/polyconf/polyconf/Polymer.py +++ b/polyconf/polyconf/Polymer.py @@ -128,20 +128,26 @@ def extend(self, monomer, n, nn, names, joins, ortho=[1,1,1], :param nn: residue (identified by its resid) to extend to (i.e. the new residue) :type nn: int - :param names: four single 'key:value' pairs with keys - P1, Q1, P2 and Q2 for the four dummy atoms that are - 'overlaid' during polymer extension. P1 is the real atom in - the polymer while P2 is its corresponding dummy atom in the - incoming monomer. Similarly, Q1 is the real atom in the - polymer while Q2 is the corresponding dummy atom in the - incoming polymer. Such that P1 and Q1 and bonded, and P2 - and Q2 have an equivalent bond. The result will be a new - bond between Q1 and the next atom beside Q2 that is NOT P2. + :param names: a dictionary of single 'key:value' pairs with keys + P, Q, R and S for the atoms that are used to define the mapping + during polymer extension. + + There are two additional optional values, V1 and V2, which are only required in linear extend. + + The definition is as follows: + + * Atoms P and R are a pair of bonded atoms located in the incoming monomer. + * Atoms Q and S are a pair of bonded atoms located in residue n of the existing polymer. + * Atoms P and Q are equivalent in the final polymer, and either P or Q is a dummy atom. + * Atoms R and S are equivalent in the final polymer, and either R or S is a dummy atom. + * Bonds PR and QS are equivalent in the final polymer. + `extend()` is agnostic as to which atom in each pair is real, and which atom is a dummy + after extension, you can use renamer() to explicilty set dummy atoms + * Atoms V1 and V2 a pair of atoms located in the incoming monomer. :type names: dict - :param joins: defines which two atoms of residue - 'n' (first value) and residue 'nn' (second value) will - be bonded - :type joins: list of a 2 item tuple + :param joins: defines a pair of atoms to connected after extension + the first atom in a pair belongs to residue 'n', and the second belongs to residue 'nn' + :type joins: list of 2 item tuples :param ortho: used for linearise and defines which axis to align along (e.g. align along x with ortho=[1,0,0]), defaults to [1,1,1] @@ -189,47 +195,47 @@ def extend(self, monomer, n, nn, names, joins, ortho=[1,1,1], # cool. # love that for me. - P1 = self.polymer.select_atoms(f'resid {n} and name {names["P1"]}').positions[-1] - Q1 = self.polymer.select_atoms(f'resid {n} and name {names["Q1"]}').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) u_.residues.resids = nn - P2 = u_.select_atoms(f'resid {nn} and name {names["P2"]}').positions[0] + P = u_.select_atoms(f'resid {nn} and name {names["P"]}').positions[0] # first, translate u_ by vector P2->P1 - T = P1 - P2 + PQ = Q - P - u_.atoms.translate(T) + u_.atoms.translate(PQ) # next, rotate around cross product of backbone vectors to align C_n to CN_n+1 - P2 = u_.select_atoms(f'resid {nn} and name {names["P2"]}').positions[0] - Q2 = u_.select_atoms(f'resid {nn} and name {names["Q2"]}').positions[0] + P = u_.select_atoms(f'resid {nn} and name {names["P"]}').positions[0] + R = u_.select_atoms(f'resid {nn} and name {names["R"]}').positions[0] - v1 = Q2 - P2 - v1_n = np.linalg.norm(v1) + PR = R - P + PR_n = np.linalg.norm(PR) - v2 = Q1 - P1 - v2_n = np.linalg.norm(v2) + QS = S - Q + QS_n = np.linalg.norm(QS) - theta = degrees(np.arccos(np.dot(v1,v2)/(v1_n * v2_n))) # TODO there are edge cases where this can fail, I think if v1 and v2 are exactly antiparallel. I need to add a checker to resolve that. + theta = degrees(np.arccos(np.dot(PR,QS)/(PR_n * QS_n))) # TODO there are edge cases where this can fail, I think if v1 and v2 are exactly antiparallel. I need to add a checker to resolve that. - k = np.cross(v1,v2) - u_r1 = u_.atoms.rotateby(theta,axis=k,point=P1) + k = np.cross(PR,QS) + u_r1 = u_.atoms.rotateby(theta,axis=k,point=P) if linearise: - R= u_.select_atoms('resid '+str(nn)+' and name '+names['R']).positions[0] - S= u_.select_atoms('resid '+str(nn)+' and name '+names['S']).positions[0] - RS=S-R - RS_n=np.linalg.norm(RS) + V1= u_.select_atoms('resid '+str(nn)+' and name '+names['V1']).positions[0] + V2= u_.select_atoms('resid '+str(nn)+' and name '+names['V2']).positions[0] + V=V2-V1 + V_n=np.linalg.norm(V) ortho_n=np.linalg.norm(ortho) - theta = degrees(np.arccos(np.dot(RS,ortho)/(RS_n * ortho_n))) + theta = degrees(np.arccos(np.dot(V,ortho)/(V_n * ortho_n))) - k=np.cross(RS,ortho) + k=np.cross(V,ortho) - u_r2 = u_r1.atoms.rotateby(theta,axis=k,point=R) + u_r2 = u_r1.atoms.rotateby(theta,axis=k,point=V1) # R= u_.select_atoms('resid '+str(nn)+' and name '+names['R']).positions[0] # S= u_.select_atoms('resid '+str(nn)+' and name '+names['S']).positions[0] @@ -253,24 +259,24 @@ def extend(self, monomer, n, nn, names, joins, ortho=[1,1,1], new.dimensions = list(new.atoms.positions.max(axis=0) + [0.5,0.5,0.5]) + [90]*3 self.polymer = new.copy() - def _split_pol(self,a1,a1_resid,a2,a2_resid): + def _split_pol(self,J,J_resid,K,K_resid): """ Given a pair of bonded atoms, uses a graph representation to identify groups of connected atoms on either side of the bond returns atomgroups corresponding to all atoms on either side of the bond. Use with caution if your polymer contains rings or closed loops, the atoms on either side of the bond must not be connected through other parts of the molecule Args: - a1 (atom name): the name of the first atom in the bond - a1_resid: the resid of the first atom in the bond - a2 (atom name): the name of the second atom in the bond - a2_resid: the resid of the second atom in the bond + J (atom name): the name of the first atom in the bond + J_resid: the resid of the first atom in the bond + K (atom name): the name of the second atom in the bond + K_resid: the resid of the second atom in the bond Returns two atomgroups, fore and aft - fore is all atoms connected to a1 - aft is all atoms connected to a2 + fore is all atoms connected to J + aft is all atoms connected to K """ - pair = self.polymer.select_atoms(f'(resid {a1_resid} and name {a1}) or (resid {a2_resid} and name {a2})' ) + pair = self.polymer.select_atoms(f'(resid {J_resid} and name {J}) or (resid {K_resid} and name {K})' ) bond = self.polymer.atoms.bonds.atomgroup_intersection(pair,strict=True)[0] g = nx.Graph() g.add_edges_from(self.polymer.atoms.bonds.to_indices()) @@ -281,21 +287,21 @@ def _split_pol(self,a1,a1_resid,a2,a2_resid): return(fore,aft) - def _rotate(self,a1,a1_resid,a2,a2_resid,mult=3,step=1): + def rotate(self,J,J_resid,K,K_resid,mult=3,step=1): """ - Given a pair of bonded atoms a1 and a2, uses _split_pol() to identify all atoms connected to a1, then rotates them around the bond by (step * int(360/mult)) degrees, chanding the state of the dihedral centered over a1-b1. + Given a pair of bonded atoms J and K within some torsion, uses _split_pol() to identify all atoms connected to J, then rotates them around vector JK by (step * int(360/mult)) degrees, rotating the dihedral centered over J-K by one step. Args: - a1 (atom name): the name of the first atom in the bond - a1_resid: the resid of the first atom in the bond - a2 (atom name): the name of the second atom in the bond - a2_resid: the resid of the second atom in the bond - mult: the multiplicity of the dihedral centered over a1-b1 + J (atom name): the name of the first atom in the bond + J_resid: the resid of the first atom in the bond + K (atom name): the name of the second atom in the bond + K_resid: the resid of the second atom in the bond + mult: the multiplicity of the dihedral centered over J-K step: how many steps around the dihedral to rotate """ - fore,_=self._split_pol(a1,a1_resid,a2,a2_resid) - pair = self.polymer.select_atoms(f'(resid {a1_resid} and name {a1}) or (resid {a2_resid} and name {a2})' ) + fore,_=self._split_pol(J,J_resid,K,K_resid) + pair = self.polymer.select_atoms(f'(resid {J_resid} and name {J}) or (resid {K_resid} and name {K})' ) bond = self.polymer.atoms.bonds.atomgroup_intersection(pair,strict=True)[0] v = bond[1].position - bond[0].position o = (fore & bond.atoms)[0] @@ -303,25 +309,25 @@ def _rotate(self,a1,a1_resid,a2,a2_resid,mult=3,step=1): fore.rotateby(rot, v, point=o.position) - def dist(self,a1,a1_resid,a2,a2_resid,dummy='X*',backwards_only=True): + def dist(self,J,J_resid,K,K_resid,dummies='X*',backwards_only=True): """ - Given a pair of bonded atoms a1 and a2, get minimum distance between + Given a pair of bonded atoms J and K, get minimum distance between atoms on one side of a bond, and atoms on the other side of the bond. This is useful for detecting overlapping atoms. - :param a1: the name of the first atom in the bond - :type a1: str - :param a1_resid: the resid of the first atom in the bond - :type a1_resid: int - :param a2: the name of the second atom in the bond - :type a2: str - :param a2_resid: the resid of the second atom in the bond - :type a2_resid: int - :param dummy: the names of dummy atoms, to be discluded from the + :param J: the name of the first atom in the bond + :type J: str + :param J_resid: the resid of the first atom in the bond + :type J_resid: int + :param K: the name of the second atom in the bond + :type K: str + :param K_resid: the resid of the second atom in the bond + :type K_resid: int + :param dummies: the names of dummy atoms, to be excluded from the distance calculation, defaults to 'X*' - :type dummy: str, optional + :type dummies: str, optional :param backwards_only: only consider atoms from residues 1 to - max([a1_resid,a2_resid]), and do not consider atoms further + max([J_resid,K_resid]), and do not consider atoms further along the chain. This is useful for solving dihedrals algorithmically, as only clashes in the solved region of the polymer are considered, defaults to True @@ -331,34 +337,34 @@ def dist(self,a1,a1_resid,a2,a2_resid,dummy='X*',backwards_only=True): between atoms on both halves of the bond :rtype: mda.analysis.distance_array """ - fore,aft=self._split_pol(a1,a1_resid,a2,a2_resid) + fore,aft=self._split_pol(J,J_resid,K,K_resid) trim='' if backwards_only: - maxres=max([a1_resid,a2_resid]) + maxres=max([J_resid,K_resid]) trim=f' and (resid 0 to {maxres})' - fore_trim=fore.select_atoms(f'(not name {dummy}) {trim}') - aft_trim=aft.select_atoms(f'(not name {dummy}) {trim}') + fore_trim=fore.select_atoms(f'(not name {dummies}) {trim}') + aft_trim=aft.select_atoms(f'(not name {dummies}) {trim}') return(distances.distance_array(fore_trim.atoms.positions, aft_trim.atoms.positions).min()) # this is the clash detection - def shuffle(self,a1,a1_resid,a2,a2_resid,dummy='X*',mult=3,cutoff=0.5,clashcheck=False,backwards_only=False): + def shuffle(self,J,J_resid,K,K_resid,dummies='X*',mult=3,cutoff=0.5,clashcheck=False,backwards_only=False): """ - Given a pair of bonded atoms a1 and a2, randomly rotates the dihedral + Given a pair of bonded atoms J and K, randomly rotates the dihedral between them a random amount with clashcheck=True, will check if the resulting structure has overlapping atoms, and will undo the rotation if so. - :param a1: the name of the first atom in the bond - :type a1: str - :param a1_resid: the resid of the first atom in the bond - :type a1_resid: int - :param a2: the name of the second atom in the bond - :type a2: str - :param a2_resid: the resid of the second atom in the bond - :type a2_resid: int - :param dummy: the names of dummy atoms, to be discluded from the + :param J: the name of the first atom in the bond + :type J: str + :param J_resid: the resid of the first atom in the bond + :type J_resid: int + :param K: the name of the second atom in the bond + :type K: str + :param K_resid: the resid of the second atom in the bond + :type K_resid: int + :param dummies: the names of dummy atoms, to be excluded from the distance calculation. Passed to :func:`dist`, defaults to 'X*' - :type dummy: str, optional - :param mult: the multiplicity of the dihedral centered over a1-b1, + :type dummies: str, optional + :param mult: the multiplicity of the dihedral centered over J-K, defaults to 3 :type mult: int, optional :param cutoff: maximum interatomic distance for atoms to be considered @@ -368,8 +374,8 @@ def shuffle(self,a1,a1_resid,a2,a2_resid,dummy='X*',mult=3,cutoff=0.5,clashcheck so reverses the shuffle, defaults to False :type clashcheck: bool, optional :param backwards_only: passed to :func:`dist` to decide if clash - checking is done along the entire polymer (if True) or just - from residue 1 up to max([a1_resid,a2_resid]) (if False), + checking is done from residue 1 up to max([J_resid,K_resid]) (if True), + or along the entire polymer (if False) defaults to False :type backwards_only: bool, optional @@ -377,26 +383,31 @@ def shuffle(self,a1,a1_resid,a2,a2_resid,dummy='X*',mult=3,cutoff=0.5,clashcheck :rtype: bool """ step = random.randrange(1,mult) # rotate fore by a random multiplicity - self._rotate(a1,a1_resid,a2,a2_resid,mult,step) - clash= (self.dist(a1,a1_resid,a2,a2_resid,dummy,backwards_only=backwards_only) <= cutoff ) + self.rotate(J,J_resid,K,K_resid,mult,step) + clash= (self.dist(J,J_resid,K,K_resid,dummies,backwards_only=backwards_only) <= cutoff ) if clash and clashcheck: - self._rotate(a1,a1_resid,a2,a2_resid,mult,-1 * step,) + self.rotate(J,J_resid,K,K_resid,mult,-1 * step,) return(clash) - def dihedral_solver(self,pairlist,dummy='X*',cutoff=0.7): + def dihedral_solver(self,pairlist,dummies='X*',cutoff=0.7,backwards_only=True): """ Converts the shuffled conformation (i.e. after :func:`shuffle`) into one without overlapping atoms by resolving dihedrals :param pairlist: list of dicts of atom pairs, created with - :func:`gen_pairlist`, e.g. [{a1,a1_resid,a2,a2_resid,mult}] + :func:`gen_pairlist`, e.g. [{J,J_resid,K,K_resid,mult}] :type pairlist: list of dicts of atom pairs - :param dummy: the names of dummy atoms, to be discluded from the + :param dummies: the names of dummy atoms, to be discluded from the distance calculation. Passed to :func:`dist`, defaults to 'X*' - :type dummy: str, optional + :type dummies: str, optional :param cutoff: maximum interatomic distance for atoms to be considered overlapping, defaults to 0.7 :type cutoff: float, optional + :param backwards_only: passed to :func:`dist` to decide if clash + checking is done from residue 1 up to max([J_resid,K_resid]) for the current dihedral (if True), + or along the entire polymer (if False) + defaults to True + :type backwards_only: bool, optional :return: True if unable to resolve dihedrals or False if dihedrals all resolved and no clashes detected @@ -416,7 +427,7 @@ def dihedral_solver(self,pairlist,dummy='X*',cutoff=0.7): with tqdm(total=steps) as pbar: while i < steps and i >=0: dh=pairlist[i] - check=self.dist(a1=dh['a1'],a1_resid=dh['a1_resid'],a2=dh['a2'],a2_resid=dh['a2_resid'],dummy=dummy) + check=self.dist(J=dh['J'],J_resid=dh['J_resid'],K=dh['K'],K_resid=dh['K_resid'],dummies=dummies,backwards_only=backwards_only) if check > cutoff and not retry: i+=1 pbar.update(1) @@ -437,7 +448,7 @@ def dihedral_solver(self,pairlist,dummy='X*',cutoff=0.7): repeats += 1 else: # no, there are more tries tries[i] += 1 - self._rotate(a1=dh['a1'],a1_resid=dh['a1_resid'],a2=dh['a2'],a2_resid=dh['a2_resid'],mult=dh['mult']) + self.rotate(J=dh['J'],J_resid=dh['J_resid'],K=dh['K'],K_resid=dh['K_resid'],mult=dh['mult']) done=True if failed or i<0: # hard coded to detect failure if you stop at i<=0 because detecting this automatically wasn't working print('Could not reach a valid conformation') @@ -446,16 +457,16 @@ def dihedral_solver(self,pairlist,dummy='X*',cutoff=0.7): else: return False - def shuffler(self,pairlist,dummy='X*',cutoff=0.5,clashcheck=False): + def shuffler(self,pairlist,dummies='X*',cutoff=0.5,clashcheck=False): """ Shuffles each pair of bonded atoms in the provided pairlist with :func:`shuffle` :param pairlist: list of dicts of atom pairs, created with - :func:`gen_pairlist`, e.g. [{a1,a1_resid,a2,a2_resid,mult}] + :func:`gen_pairlist`, e.g. [{J,J_resid,K,K_resid,mult}] :type pairlist: list of dicts of atom pairs - :param dummy: _the names of dummy atoms, to be discluded from the + :param dummies: _the names of dummy atoms, to be discluded from the distance calculation. Passed to :func:`shuffle`, defaults to 'X*' - :type dummy: str, optional + :type dummies: str, optional :param cutoff: maximum interatomic distance for atoms to be considered overlapping, defaults to 0.5 :type cutoff: float, optional @@ -464,17 +475,17 @@ def shuffler(self,pairlist,dummy='X*',cutoff=0.5,clashcheck=False): :type clashcheck: bool, optional """ for dh in tqdm(pairlist): - self.shuffle(a1=dh['a1'],a1_resid=dh['a1_resid'],a2=dh['a2'],a2_resid=dh['a2_resid'],mult=dh['mult'],dummy=dummy,clashcheck=clashcheck,cutoff=cutoff) + self.shuffle(J=dh['J'],J_resid=dh['J_resid'],K=dh['K'],K_resid=dh['K_resid'],mult=dh['mult'],dummies=dummies,clashcheck=clashcheck,cutoff=cutoff) - def gen_pairlist(self,a1,a2,first_resid=1,last_resid=999,resid_step=1,same_res=True,mult=3): + def gen_pairlist(self,J,K,first_resid=1,last_resid=999,resid_step=1,same_res=True,mult=3): """ Generates a list of dicts of atom pairs for use in :func:`shuffler` or :func:`dihedral_solver`. - :param a1: the name of the first atom in the bond - :type a1: str - :param a2: the name of the second atom in the bond - :type a2: str + :param J: the name of the first atom in the bond + :type J: str + :param K: the name of the second atom in the bond + :type K: str :param first_resid: resid of first residue to include in pairlist, defaults to 1 :type first_resid: int, optional @@ -486,23 +497,24 @@ def gen_pairlist(self,a1,a2,first_resid=1,last_resid=999,resid_step=1,same_res=T defaults to 1 :type resid_step: int, optional :param same_res: if True, paired atoms are bonded within the same - residue, if False assumes that paired atoms a1 and a2 form a + residue, if False assumes that paired atoms J and K form a bond between residue n and residue n+1, defaults to True :type same_res: bool, optional - :param mult: the multiplicity of the dihedral centered over a1-b1, + :param mult: the multiplicity of the dihedral centered over J-K, defaults to 3 :type mult: int, optional :return: a list of dicts of atom pairs :rtype: list of dicts """ - pairlist=[{'a1':a1,'a1_resid':i,'a2':a2,'a2_resid':i+int(not same_res),'mult':mult} for i in range(first_resid,last_resid+1,resid_step)] + pairlist=[{'J':J,'J_resid':i,'K':K,'K_resid':i+int(not same_res),'mult':mult} for i in range(first_resid,last_resid+1,resid_step)] # crude, doesn't actually check if the pairs exist # want to extend this to cover more than one dihedral return(pairlist) def gencomp(mdict, length, fill, middle, count = True, frac = False): """ + Deprecated and unsupported, use with caution. Generate the linear conformation of a polymer from a given length, specified monomer composition and a dictionary of monomers. Only used by 'polyconf_automatic' and should not be used otherwise. diff --git a/polyconf_examples/build_PEI_linear.py b/polyconf_examples/build_PEI_linear.py index 82cf2ba..3a18aa5 100644 --- a/polyconf_examples/build_PEI_linear.py +++ b/polyconf_examples/build_PEI_linear.py @@ -11,32 +11,58 @@ # # ############################################################################### -polymer1=Polymer(Monomer('PEI_start.pdb')) # initialise -imax=127 # we will lay 127 additional monomers +# 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') -for i in range (0,imax): - if not i==imax: # do this for everything except the first one - polymer1.extend( # extend with one monomer, aligned along this step's linearization vector +# 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(), # extend existing residue i - nn=polymer1.newresid(), # incoming monomer will have resid i+1 - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) - else: - # add final monomer - polymer1.extend(Monomer('PEI_end.pdb'),i,i+1, - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), - joins=[('N1','C1')],) + + +# 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() # center in box -Saver.save(dummyAtoms='CX NX',fname='polymer_01_vanilla-extend') # save, excluding dummy atoms + +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 ############################################################################### # # @@ -44,37 +70,47 @@ # # ############################################################################### +# 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 -imax=127 # lay 127 additional monomers -ortho=[1,0,0] # linearization vector +polymer2=Polymer(Monomer('PEI_start.pdb')) # initialise the polymer the same way -for i in range (1,imax+1): - if not i==imax: # do this for everything except the first one - polymer2.extend( - Monomer('PEI_monomer.pdb'), # extend with this monomer - i, # extend existing residue i - i+1, # incoming monomer will have resid i+1 - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX', # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i - R='C1',S='N1'), # monomer vector is vector C1_i+1 N1_i+1, +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=ortho,# specify linearization vector - joins=[('N1','C1')],# new connection between N1_i and C1_i+1 + 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')], ) - else: - # add final monomer - polymer2.extend(Monomer('PEI_end.pdb'),i,i+1, - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), - joins=[('N1','C1')],) -# polymer.genconf([('C2','N1')],length=10, cutoff = 0.5) - 'C1' 'NX' 'CX' 'N1' + + +# 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 linear. +# 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. @@ -88,52 +124,70 @@ # # ############################################################################### -# In this example, we will use a linear extend, but we will alternate between two vectors -# the first linear extend will use vector A -# the second linear extend will use vector B +# 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 -imax=127 # lay 24 additional monomers -alternator=True # lay monomers along alternating vectors to get a clean linear conformation, because shuffle isn't working - - - -for i in range (1,imax+1): - if not i==imax: # do this for everything except the first one - if alternator: - ortho=[1,0,0] # linearization vector A - else: - ortho=[0,1,0] # linearization vector B - alternator=not alternator # flip the alternator +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')], + ) - polymer3.extend(Monomer( # extend with one monomer, aligned along this step's linearization vector - 'PEI_monomer.pdb'), # extend with this monomer - i, # extend existing residue i - i+1, # incoming monomer will have resid i+1 - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX', # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i - R='C1',S='N1'), # monomer vector is vector C1_i+1 N1_i+1, - linearise=True, # linearize incoming monomer so monomer vector lies along along ortho - ortho=ortho,# specify linearization vector - joins=[('N1','C1')],# new connection between N1_i and C1_i+1 +# 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')], ) - else: - # add final monomer - polymer3.extend(Monomer('PEI_end.pdb'),i,i+1, - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), - joins=[('N1','C1')],) -# polymer.genconf([('C2','N1')],length=10, cutoff = 0.5) - 'C1' 'NX' 'CX' 'N1' + +# 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 more realistic. +# 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. ############################################################################### # # @@ -143,46 +197,48 @@ # in this example, we will make a copy of polymer 1, then use dihedral_solver() to generate a more reasonable starting conformation -polymer4=Polymer(Monomer('PEI_start.pdb')) # initialise -imax=127 # lay 127 additional monomers +polymer4=polymer1.copy() # make a copy of our first polymer, which was a tightly packed helix -# TODO replace with deepcopy +# 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 -for i in range (1,imax+1): - if not i==imax: # do this for everything except the first one - polymer4.extend( # extend with one monomer, aligned along this step's linearization vector - Monomer('PEI_monomer.pdb'), # extend with this monomer - i, # extend existing residue i - i+1, # incoming monomer will have resid i+1 - names=dict(P1='CX',P2='C1',Q1='N1',Q2='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 - ) - else: - # add final monomer - polymer4.extend(Monomer('PEI_end.pdb'),i,i+1, - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), - joins=[('N1','C1')],) +# 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 -CN=polymer4.gen_pairlist(a1='C2',a2='N1',first_resid=1,last_resid=128,mult=3) +# 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 -polymer4.dihedral_solver(CN,dummy='CX,NX',cutoff=0.7) +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, the algorithmic solver has still resulted in a highly ordered structure. +# 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() # @@ -194,42 +250,40 @@ #first, make a copy of polymer one -polymer5=Polymer(Monomer('PEI_start.pdb')) # initialise -imax=127 # lay 127 additional monomers +polymer5=polymer1.copy() -for i in range (1,imax+1): - if not i==imax: # do this for everything except the first one - polymer5.extend( # extend with one monomer, aligned along this step's linearization vector - Monomer('PEI_monomer.pdb'), # extend with this monomer - i, # extend existing residue i - i+1, # incoming monomer will have resid i+1 - names=dict(P1='CX',P2='C1',Q1='N1',Q2='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 - ) - else: - # add final monomer - polymer5.extend(Monomer('PEI_end.pdb'),i,i+1, - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), - joins=[('N1','C1')],) +import random +random.seed(8) # by setting the random seed we can ensure the random conformations are replicable. -import random +# 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 -random.seed(8) # by setting the random seed we can ensure the random conformations are replicable +alkane_dhedrals=polymer5.gen_pairlist(J='C1',K='C2',first_resid=1,last_resid=128,mult=3) -alkanes=polymer5.gen_pairlist(a1='C1',a2='C2',first_resid=1,last_resid=128,mult=3) -NC=polymer5.gen_pairlist(a1='N1',a2='C1',first_resid=1,last_resid=127,mult=6,same_res=False) +polymer5.shuffler(NC_dihedrals,) # this generates a random conformation -polymer5.shuffler(NC,) # 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 - -polymer5.dihedral_solver(alkanes,dummy='CX NX',cutoff=1) # this converts the shuffled conformation into one without overlapping 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 @@ -239,4 +293,6 @@ # 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. \ No newline at end of file +# 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/build_PMMA_01-isotactic.py b/polyconf_examples/build_PMMA_01-isotactic.py new file mode 100644 index 0000000..a5598ae --- /dev/null +++ b/polyconf_examples/build_PMMA_01-isotactic.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB +import random +random.seed(1) + +############################################################################### +# # +# Polymer one: Building an isotactic PMMA polymer # +# # +############################################################################### + +print ('creating an isotactic PMMA polymer') + +dummies="CMA CN CP CQ" + +PMMA_isotactic=Polymer(Monomer('MMAD_bonds.pdb')) # initialise + +adds=49 # we will lay 49 additional monomers + +for i in range (0,adds): + PMMA_isotactic.extend( # extend with one monomer, aligned along this step's linearization vector + 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] + ) + +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 + +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 + +Saver = PDB(PMMA_isotactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_isotactic_vanilla-solved') # save, excluding dummies atoms + +PMMA_isotactic.shuffler(alkanes) +PMMA_isotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms + +Saver = PDB(PMMA_isotactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_isotactic_shuffled_and_solved') # save, excluding dummies atoms + + diff --git a/polyconf_examples/build_PMMA_02-syndiotactic.py b/polyconf_examples/build_PMMA_02-syndiotactic.py new file mode 100644 index 0000000..249aab2 --- /dev/null +++ b/polyconf_examples/build_PMMA_02-syndiotactic.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB +import random +random.seed(1) + + +############################################################################### +# # +# Polymer two: Building a syndiotactic PMMA polymer # +# # +############################################################################### + +print ('creating a syndiotactic PMMA polymer') + +dummies="CMA CN CP CQ" + +Monomer_D='MMAD_bonds.pdb' +Monomer_L='MMAL_bonds.pdb' +alternator=True + +PMMA_syndiotactic=Polymer(Monomer(Monomer_D)) # initialise + +adds=49 # we will lay 49 additional monomers + +for i in range (0,adds): + if alternator: + monomer=Monomer_L + else: + monomer=Monomer_D + alternator=not alternator + PMMA_syndiotactic.extend( # extend with one monomer, aligned along this step's linearization vector + Monomer(monomer), # extend with this monomer + n=PMMA_syndiotactic.maxresid(), # extend existing residue i + nn=PMMA_syndiotactic.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 + joins=[('C','CA')],# new connection between N1_i and C1_i+1 + ) + +Saver = PDB(PMMA_syndiotactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_syndiotactic_vanilla-extend') # save, excluding dummies atoms + +alkanes=PMMA_syndiotactic.gen_pairlist(J='C',K='CA',first_resid=1,same_res=False,last_resid=49,mult=3) +sidechains=PMMA_syndiotactic.gen_pairlist(J='CA',K='CB',first_resid=1,same_res=True,last_resid=50,mult=6) + +PMMA_syndiotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms +PMMA_syndiotactic.dihedral_solver(sidechains,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms + +Saver = PDB(PMMA_syndiotactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_syndiotactic_vanilla-solved') # save, excluding dummies atoms + + +PMMA_syndiotactic.shuffler(alkanes) +PMMA_syndiotactic.dihedral_solver(alkanes,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms +PMMA_syndiotactic.shuffler(sidechains) +PMMA_syndiotactic.dihedral_solver(sidechains,dummies=dummies,cutoff=1.1) # this converts the shuffled conformation into one without overlapping atoms + +Saver = PDB(PMMA_syndiotactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_syndiotactic_shuffled_and_solved') # save, excluding dummies atoms diff --git a/polyconf_examples/build_PMMA_03-atactic.py b/polyconf_examples/build_PMMA_03-atactic.py new file mode 100644 index 0000000..87b5928 --- /dev/null +++ b/polyconf_examples/build_PMMA_03-atactic.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB +import random +random.seed(10) + +############################################################################### +# # +# Polymer three: Building an atactic PMMA polymer # +# # +############################################################################### + +print ('creating an atactic PMMA polymer') + +dummies="CMA CN CP CQ" +Monomer_D='MMAD_bonds.pdb' +Monomer_L='MMAL_bonds.pdb' + +composition=25*[Monomer_D] + 25*[Monomer_L] # a list containing 50 monomers, 25 with each taticity +composition = random.sample(composition,len(composition)) # randomise the order of the monomers + +PMMA_atactic=Polymer(Monomer(composition[0])) # initialise with first monomer + +for monomer in composition[1:]: # extend with the remain 49 monomers + PMMA_atactic.extend( # extend with one monomer, aligned along this step's linearization vector + Monomer(monomer), # extend with this monomer + n=PMMA_atactic.maxresid(), # extend existing residue i + nn=PMMA_atactic.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, + joins=[('C','CA')],# new connection between N1_i and C1_i+1 + ) + +Saver = PDB(PMMA_atactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_atactic_vanilla-extend') # save, excluding dummies atoms + +# have a look at this initial polymer, compared to the vanilla isotactic and syndiotactic polymers +# because of the particular conformations of the two monomers, and the random order of the monomer composition +# the backbone of this polymer is much more coiled.absyou might think this is desifreable, as it reflects a less ordered starting conformation +# however, many monomers have overlapping atoms, and it is unlikely that we could minimise this starting confromation + +alkanes=PMMA_atactic.gen_pairlist(J='CA',K='C',first_resid=1,same_res=True,last_resid=50,mult=6) +sidechains=PMMA_atactic.gen_pairlist(J='CA',K='CB',first_resid=1,same_res=True,last_resid=50,mult=6) + +dh=[] #combine alkanes and sidechains into one list of dihedrals, so we can shuffle the alkane and dihedral in each monomer in order + +for i in range(0,50): + dh += [alkanes[i]] + dh += [sidechains[i]] + +print('attempting to solve initial conformation') + +PMMA_atactic.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +# this failed due to overlapping atoms + +# we can attempt to rotate and shuffle dihedrals to find a valid conformation +# however the brute force search of dihedral space used by dihedral_solver() will fail for this particular random +# it is likely that repeated applications of shuffle and dihedral_solver will eventually find a reasonable conformation +# but we cannot know how many attempts this will take + +print('shuffle attempt 01') + +PMMA_atactic.shuffler(dh) + +PMMA_atactic.dihedral_solver(dh,dummies=dummies,cutoff=1.1 ) + +Saver = PDB(PMMA_atactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_atactic_shuffle_01') # save, excluding dummies atoms + +print('shuffle attempt 02') + +PMMA_atactic.shuffler(dh) +PMMA_atactic.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_atactic_shuffle_02') # save, excluding dummies atoms + +print('shuffle attempt 03') + +PMMA_atactic.shuffler(dh) + +PMMA_atactic.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_atactic_shuffle_03') # save, excluding dummies atoms + + +print('shuffle attempt 04') + +PMMA_atactic.shuffler(dh) + +PMMA_atactic.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic) +Saver.cleanup() # center in box +Saver.save(dummies=dummies,fname='PMMA_atactic_shuffle_04') # save, excluding dummies atoms + +# even after four steps, this is still unable to solve the conformation + +# If your initial conformation contains many overlapping monomers, I have found you are much less likely to find valid starting conformations +# There are several steps we could take to resolve this +# One approach would be to reduce the cutoff, which allows conformations with smaller interatomic distances + +# In this example, using a cutoff of 0.8 A will allow dihedral_solver() to generate a starting starting conformation after every step + +# Reducing the cutoff means some atoms are likely to be unrealistically close, but an energy minimzation may address this and result in structures that can be used in molecular dynamics + +# Another approach is to vary the multiplicity of your dihedrals, to alter the search parameters through dihedral space +# In this case the alkane dihedrals have been permitted to explore a multiplicity of 6, but changing it to 3 (restricting the search space) or to 12 (increasing the search space) may facilitate a solution + +# You could also try solving the alkanes first, and then the sidechains, or doing the reverse + +# By varying your approach, you will eventually stumble upon a method that identifies valid starting conformations + +# However, I have found it is quite reliable to start by generating a psuedolinear conformation, and then shuffling and solving that conformation +# An example of this process is given in build_PMMA_04-atactic-linearized \ No newline at end of file diff --git a/polyconf_examples/build_PMMA_04-atactic-linearized.py b/polyconf_examples/build_PMMA_04-atactic-linearized.py new file mode 100644 index 0000000..6e51fbf --- /dev/null +++ b/polyconf_examples/build_PMMA_04-atactic-linearized.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python + +from polyconf.Monomer import Monomer +from polyconf.Polymer import Polymer +from polyconf.PDB import PDB +import random +random.seed(0) + +print ('creating an atactic PMMA polymer using a linear extend()') + +dummies="CMA CN CP CQ" + +Monomer_D='MMAD_bonds.pdb' +Monomer_L='MMAL_bonds.pdb' + +composition=25*[Monomer_D] + 25*[Monomer_L] # a list containing 50 monomers, 25 with each taticity +composition = random.sample(composition,len(composition)) # randomise the order of the monomers + +PMMA_atactic_linear=Polymer(Monomer(composition[0])) # initialise with first monomer + +for monomer in composition[1:]: # extend with the remain 49 monomers + PMMA_atactic_linear.extend( # extend with one monomer, aligned along this step's linearization vector + Monomer(monomer), # extend with this monomer + n=PMMA_atactic_linear.maxresid(), # extend existing residue i + nn=PMMA_atactic_linear.newresid(), # incoming monomer will have resid i+1 + names=dict(Q='CA',P='CMA',S='C',R='CN',V1='CN',V2='C'), # C1_i+1 fit to CX_i, + joins=[('C','CA')],# new connection between N1_i and C1_i+1 + linearise=True, + ortho=[1,0,0] + ) + +Saver = PDB(PMMA_atactic_linear) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-extend') # save, excluding dummies atoms + +# when you look at this polymer, you can see the alkane backbone is extremely regular +# while this is not realistic, we can use it to generate reasonable conformations + +alkanes=PMMA_atactic_linear.gen_pairlist(J='CA',K='C',first_resid=1,same_res=True,last_resid=50,mult=6) +sidechains=PMMA_atactic_linear.gen_pairlist(J='CA',K='CB',first_resid=1,same_res=True,last_resid=50,mult=6) + +dh=[] #combine alkanes and sidechains into one list of dihedrals, so we can shuffle the alkane and dihedral in each monomer in order + +for i in range(0,50): + dh += [alkanes[i]] + dh += [sidechains[i]] + +print ('solving pseudolinear atactic PMMA ') + +PMMA_atactic_linear.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +# unsurprisingly it was very easy to solve this pseudolinear conformation + +Saver = PDB(PMMA_atactic_linear) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-solved') # save, excluding dummies atoms + +print ('shuffle attempt 01') + +PMMA_atactic_linear_01 = PMMA_atactic_linear.copy() + +PMMA_atactic_linear_01.shuffler(dh) +PMMA_atactic_linear_01.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic_linear_01) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-shuffled_and_solved_01') # save, excluding dummies atoms + +# this generated a valid conformation +# lets repeat the process three more times + +print ('shuffle attempt 02') + +PMMA_atactic_linear_02 = PMMA_atactic_linear.copy() + +PMMA_atactic_linear_02.shuffler(dh) +PMMA_atactic_linear_02.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic_linear_02) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-shuffled_and_solved_02') # save, excluding dummies atoms + + +print ('shuffle attempt 03') + +PMMA_atactic_linear_03 = PMMA_atactic_linear.copy() + +PMMA_atactic_linear_03.shuffler(dh) +PMMA_atactic_linear_03.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic_linear_03) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-shuffled_and_solved_03') # save, excluding dummies atoms + +print ('shuffle attempt 04') + +PMMA_atactic_linear_04 = PMMA_atactic_linear.copy() + +PMMA_atactic_linear_04.shuffler(dh) +PMMA_atactic_linear_04.dihedral_solver(dh,dummies=dummies,cutoff=1.1) + +Saver = PDB(PMMA_atactic_linear_04) +Saver.cleanup() +Saver.save(dummies=dummies,fname='PMMA_atactic_linear-shuffled_and_solved_04') # save, excluding dummies atoms + +# all four attempts were able to be solved, generating four possible starting conformations for energy minimzation From 38390ee23d2a281d35fb66c5832a2bb0e5781cd4 Mon Sep 17 00:00:00 2001 From: recombinatrix Date: Mon, 17 Feb 2025 15:23:43 +1100 Subject: [PATCH 2/3] resolve conflicts, i hate github --- polyconf/polyconf/PDB.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/polyconf/polyconf/PDB.py b/polyconf/polyconf/PDB.py index ea3b099..2a00895 100644 --- a/polyconf/polyconf/PDB.py +++ b/polyconf/polyconf/PDB.py @@ -55,21 +55,9 @@ def save(self, dummies="X*", fname="polymer", selectionString = None, gmx = Fals """ if selectionString: if gmx: -<<<<<<< HEAD self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.gro") else: self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.pdb") -======= -<<<<<<< HEAD - self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.gro") - else: - self.select_atoms(f"{selectionString} and not name {dummyAtoms}")._write(f"{fname}.pdb") -======= - self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.gro") - else: - self.select_atoms(f"{selectionString} and not name {dummies}").atoms._write(f"{fname}.pdb") ->>>>>>> d62405c (rewrote all polyconf functions to use updated atom names) ->>>>>>> 9e689d2 (Testing rebase against main worked) else: if gmx: self._write(f"not name {dummies}", f"{fname}.gro") From 120144bdd5d5e13874c119b86d4b8863862a02b7 Mon Sep 17 00:00:00 2001 From: recombinatrix Date: Mon, 17 Feb 2025 15:28:13 +1100 Subject: [PATCH 3/3] update tests for new nomenclature --- tests/polyconf/test_polymer.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/tests/polyconf/test_polymer.py b/tests/polyconf/test_polymer.py index b54e31d..32f2520 100644 --- a/tests/polyconf/test_polymer.py +++ b/tests/polyconf/test_polymer.py @@ -99,7 +99,7 @@ def test_polymer_extend(data_dir): Monomer(data_dir / 'PEI_monomer.pdb'), # extend with this monomer n=polymer.maxresid(), # we will allways add onto the existing monomer with the highest resid nn=polymer.newresid(), # the incoming monomer needs a new resid - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) @@ -107,7 +107,7 @@ def test_polymer_extend(data_dir): Monomer(data_dir / 'PEI_end.pdb'), n=polymer.maxresid(), nn=polymer.newresid(), - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), + names=dict(Q='CX',P='C1',S='N1',R='NX'), joins=[('N1','C1')], ) @@ -126,7 +126,7 @@ def test_polymer_maxresid(data_dir): Monomer(data_dir / 'PEI_monomer.pdb'), # extend with this monomer n=polymer.maxresid(), # we will allways add onto the existing monomer with the highest resid nn=polymer.newresid(), # the incoming monomer needs a new resid - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) @@ -134,7 +134,7 @@ def test_polymer_maxresid(data_dir): Monomer(data_dir / 'PEI_end.pdb'), n=polymer.maxresid(), nn=polymer.newresid(), - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), + names=dict(Q='CX',P='C1',S='N1',R='NX'), joins=[('N1','C1')], ) @@ -156,7 +156,7 @@ def test_polymer_newresid(data_dir): Monomer(data_dir / 'PEI_monomer.pdb'), # extend with this monomer n=polymer.maxresid(), # we will allways add onto the existing monomer with the highest resid nn=polymer.newresid(), # the incoming monomer needs a new resid - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) assert polymer.newresid() == i+3 @@ -165,7 +165,7 @@ def test_polymer_newresid(data_dir): Monomer(data_dir / 'PEI_end.pdb'), n=polymer.maxresid(), nn=polymer.newresid(), - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), + names=dict(Q='CX',P='C1',S='N1',R='NX'), joins=[('N1','C1')], ) @@ -183,7 +183,7 @@ def test_polymer_dihedralsolver_and_genpairlist(data_dir): Monomer(data_dir / 'PEI_monomer.pdb'), # extend with this monomer n=polymer.maxresid(), # we will allways add onto the existing monomer with the highest resid nn=polymer.newresid(), # the incoming monomer needs a new resid - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) @@ -191,19 +191,19 @@ def test_polymer_dihedralsolver_and_genpairlist(data_dir): Monomer(data_dir / 'PEI_end.pdb'), n=polymer.maxresid(), nn=polymer.newresid(), - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), + names=dict(Q='CX',P='C1',S='N1',R='NX'), joins=[('N1','C1')], ) - CN_dihedrals=polymer.gen_pairlist(a1='C2',a2='N1',first_resid=1,last_resid=130,mult=3) + CN_dihedrals=polymer.gen_pairlist(J='C2',K='N1',first_resid=1,last_resid=130,mult=3) cutoff = 0.7 - polymer.dihedral_solver(CN_dihedrals,dummy='CX,NX',cutoff=cutoff) + polymer.dihedral_solver(CN_dihedrals,dummies='CX,NX',cutoff=cutoff) for dh in CN_dihedrals: # check that the dihedral_solver has removed all clashes - assert polymer.dist('C2',dh['a1_resid'],'N1',dh['a2_resid'],'CX,NX',backwards_only=False) >= cutoff + assert polymer.dist('C2',dh['J_resid'],'N1',dh['K_resid'],'CX,NX',backwards_only=False) >= cutoff def test_polymer_shuffler(data_dir): """ @@ -216,7 +216,7 @@ def test_polymer_shuffler(data_dir): Monomer(data_dir / 'PEI_monomer.pdb'), # extend with this monomer n=polymer.maxresid(), # we will allways add onto the existing monomer with the highest resid nn=polymer.newresid(), # the incoming monomer needs a new resid - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), # C1_i+1 fit to CX_i, then rotate so NX_i+1 fit to N1_i + 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 ) @@ -224,21 +224,21 @@ def test_polymer_shuffler(data_dir): Monomer(data_dir / 'PEI_end.pdb'), n=polymer.maxresid(), nn=polymer.newresid(), - names=dict(P1='CX',P2='C1',Q1='N1',Q2='NX'), + names=dict(Q='CX',P='C1',S='N1',R='NX'), joins=[('N1','C1')], ) - CN_dihedrals=polymer.gen_pairlist(a1='C2',a2='N1',first_resid=1,last_resid=130,mult=3) + CN_dihedrals=polymer.gen_pairlist(J='C2',K='N1',first_resid=1,last_resid=130,mult=3) distances_before = [] for dh in CN_dihedrals: - distances_before.append(polymer.dist('C2',dh['a1_resid'],'N1',dh['a2_resid'],'CX,NX',backwards_only=False)) + distances_before.append(polymer.dist('C2',dh['J_resid'],'N1',dh['K_resid'],'CX,NX',backwards_only=False)) polymer.shuffler(CN_dihedrals) for i, dh in enumerate(CN_dihedrals): # check that the shuffler has chaged all the dihedrals - assert polymer.dist('C2',dh['a1_resid'],'N1',dh['a2_resid'],'CX,NX',backwards_only=False) != distances_before[i] + assert polymer.dist('C2',dh['J_resid'],'N1',dh['K_resid'],'CX,NX',backwards_only=False) != distances_before[i] # TODO: Add more tests for other methods and functionalities of the Polymer class as needed \ No newline at end of file