Chemically Realistic Tetrahedral Lattice Models for ... - ACS Publications

Article pubs.acs.org/JCTC

Chemically Realistic Tetrahedral Lattice Models for Polymer Chains: Application to Polyethylene Oxide Johannes C. B. Dietschreit,† Dennis J. Diestler,†,‡ and Ernst W. Knapp*,† †

Department of Biology, Chemistry and Pharmacy, Institute of Chemistry and Biochemistry, Freie Universität Berlin, Fabeckstrasse 36A, D-14195 Berlin, Germany ‡ University of Nebraska-Lincoln, Lincoln, Nebraska 68583, United States S Supporting Information *

ABSTRACT: To speed up the generation of an ensemble of poly(ethylene oxide) (PEO) polymer chains in solution, a tetrahedral lattice model possessing the appropriate bond angles is used. The distance between noncovalently bonded atoms is maintained at realistic values by generating chains with an enhanced degree of self-avoidance by a very efficient Monte Carlo (MC) algorithm. Potential energy parameters characterizing this lattice model are adjusted so as to mimic realistic PEO polymer chains in water simulated by molecular dynamics (MD), which serves as a benchmark. The MD data show that PEO chains have a fractal dimension of about two, in contrast to self-avoiding walk lattice models, which exhibit the fractal dimension of 1.7. The potential energy accounts for a mild hydrophobic effect (HYEF) of PEO and for a proper setting of the distribution between trans and gauche conformers. The potential energy parameters are determined by matching the Flory radius, the radius of gyration, and the fraction of trans torsion angles in the chain. A gratifying result is the excellent agreement of the pair distribution function and the angular correlation for the lattice model with the benchmark distribution. The lattice model allows for the precise computation of the torsional entropy of the chain. The generation of polymer conformations of the adjusted lattice model is at least 2 orders of magnitude more efficient than MD simulations of the PEO chain in explicit water. This method of generating chain conformations on a tetrahedral lattice can also be applied to other types of polymers with appropriate adjustment of the potential energy function. The efficient MC algorithm for generating chain conformations on a tetrahedral lattice is available for download at https://github.com/Roulattice/Roulattice.

■

INTRODUCTION Polymer chains are frequently represented formally as selfavoiding walks (SAW) on a lattice. In the simplest version, the chain is generated stochastically by placing atoms sequentially on the sites of a regular, usually simple cubic, lattice under the constraints that the distance between successive atoms is fixed (the covalent bond length) and that no two atoms can occupy the same site. Such polymer models can lead to valuable insights into phenomena such as universality class behavior1,2 and characterization of the glass transition,3 which are fundamental to the physics of polymers. It is unfortunate that such models are often inappropriate for the description of specific polymers of interest to the chemical community. To overcome this problem, rotational isomeric state (RIS) polymer models were used on high-coordination lattices.4 These coarse-grained lattice polymer models can be converted to atomistic models in continuous space by reverse mapping, thus allowing efficient Monte Carlo (MC) simulation of large polymer systems.5 In the following, the RIS polymer model was further developed and adjusted to specific polymers as, for instance, polyethylene oxide (PEO).6−8 An alternative approach uses special localized MC moves in (ϕ,ψ)-torsion angle space to model the flexibility of a polypeptide backbone to simulate protein folding.9−14 © 2016 American Chemical Society

Here, we propose a modified chain model on a diamond lattice that accounts for all atoms of the chain backbone explicitly. Its energy terms are therefore more easily interpretable on an atomistic level than are those of a coarsegrained model. We demonstrate that the properties of the lattice model can reproduce those of realistic polymer chains of chemical relevance. To guarantee realistic values for the angles between successive covalent bonds of the backbone of a polymer chain, we employ a tetrahedral lattice. This lattice accounts for the sp3-hybridization of atomic orbitals belonging to the backbone atoms of polymer models. Examples from the broad class of chemically significant polymers to which the tetrahedral lattice model applies are polyethylene, PEO, and polyamides, whose “backbone” atoms are a combination carbon, nitrogen, or oxygen. In these polymers, the angle between successive covalent bonds is generally very near the tetrahedral angle of 109.47°. A critical deficiency of the conventional self-avoiding random walk (henceforward referred to as SAW1) is that it permits conformations in which two noncovalently bonded (to which Received: February 8, 2016 Published: April 5, 2016 2388

DOI: 10.1021/acs.jctc.6b00144 J. Chem. Theory Comput. 2016, 12, 2388−2400

Article

Journal of Chemical Theory and Computation

atom pair correlation function, angular correlation, and SASA are presented.

we refer simply as nonbonded for brevity’s sake) atoms can occupy neighboring lattice sites (i.e., sites separated just by a single bond length). However, in a realistic atomistic description of polymers, the mean distance between nonbonded atoms should be considerably greater than the bond length. However, if the polymer model is a coarse-grained description, where occupied lattice sites represent polyatomic segments, i.e. monomers,15 rather than single atoms, the SAW1 model may be useful, provided that the length of the monomer (i.e., the distance between nearest-neighbor lattice sites) is close to its mean cross-sectional dimension.16,17 Indeed, such coarsegrained polymer models have been applied to both chemical18−23 and biological24−27 processes. They can also provide useful descriptions of nanoscopic systems28−30 and multivalent interactions in biology and medicine.31−33 To correct the deficiency of the conventional SAW1 model for polymers, where distances between nonbonded atoms can be as short as covalent bond lengths, we recently introduced a hierarchy of SAWn models.34 These models are characterized by the degree of self-avoidance, which increases systematically with the index n, which governs the distance of closest approach permitted between nonbonded pairs of atoms on the tetrahedral lattice. In the SAW3 models, for example, the distance of closest approach of two nonbonded atoms is (11/ 3)1/2 lbond = 1.915 lbond, where lbond is the covalent bond length, i.e., the lattice constant.34 Our previous work focused primarily on the development of efficient MC algorithms for the generation of SAWn chains with different degrees (n) of self-avoidance.34 A comparison with the results of preliminary molecular dynamics (MD) simulations of PEO polymer models in explicit water with the results of systematic studies of SAWn models for n ≤ 3 suggests that the lattice models overemphasize the repulsive character of nonbonded interactions. As a consequence, characteristic quantities of the SAWn models, such as radius of gyration and Flory radius, are too large, whereas the fraction of trans torsions is too small. In the present article, extensive MD simulations of PEO chains in explicit water up to the microsecond time scale were performed to provide a benchmark for modifying the SAW3 model by the addition of appropriate potential energy terms. The SAW3 model was chosen because it possesses an appropriate minimum distance between nonbonded pairs of atoms. The added potential energy conduces to a more compact structure through attractive interactions and to a correct ratio between gauche and trans conformers through torsional interactions. The added potential is characterized by two parameters, which are determined by fitting selected properties (radius of gyration, Flory radius, and fraction of trans conformers) of the lattice model to those of the benchmark, i.e., the realistic PEO chain model simulated by MD. Besides the properties used for the fitting of the potential energy parameters, the entropy loss of the polymer chain, the atom pair distribution function, the angular correlation of bond vectors, and the solvent accessible surface area (SASA) due to the added potential energy are also considered. The remainder of the article is organized as follows. The realistic PEO chain and its lattice counterpart are described. The properties used for determining the parameters of the added potential energy are discussed. Details of the numerical methods are presented. The efficiencies of the lattice model and the MD benchmark are presented. Results on chain entropy,

■

MATERIALS AND METHODS Polymer Chain Models. The focus of the present study is on a single polymer chain in a good solvent. The ends of the polymer chain may be functionalized, serving, for instance, as linkers that connect ligands to form multivalent binding units.31 Properties of the polymer chain, represented by an appropriate microscopic model, are computed by means of classical statistical thermodynamics. Ensembles of chain conformations are generated, and the required statistical mechanical averages are performed. Model of Polyethylene Oxide for MD Simulation. For a realistic polymer model, which can serve as a benchmark, we consider a single chain of polyethylene oxide (PEO) immersed in a large box of water molecules with periodic boundary conditions. The PEO chain alone is represented by H − (C(1)H 2 − O(2) − C(3)H 2) − (C(4)H 2 − O(5) − C(6)H 2) ···(C(3n − 2)H 2 − O(3n − 1) − C(3n)H 2) − H

(1)

where the superscripts label the backbone atoms, and n is the number of − CH2−O−CH2− monomer units. The number of backbone atoms is therefore Na = 3n. We generate the ensemble of chain conformations by MD simulation, using a potential energy that has been optimized for PEO (see Computer Simulations Section ).35,36 MD simulations were carried out for four PEO chains comprising 9, 17, 25, and 33 monomers and corresponding chain lengths of 27, 51, 75, and 99 backbone atoms, respectively. Details of the preparation and conditions for MD simulation are given below. Lattice Model of PEO Chain. Only the polymer chain itself is simulated (i.e., PEO is represented on the lattice only by its heavy backbone atoms). The chains are grown by means of an MC algorithm on a tetrahedral (diamond) lattice, which enforces the tetrahedral bond angle of 109.47° for sp3hybridized atoms. This value is very close to the bond angles of PEO. For instance, the C−O−C and the C−C−O bond angles are 109.7° and 111.5°, respectively.35 Chains constructed on this lattice have a single bond length, lbond, unlike the chains simulated by MD, where the C−O and the C−C bond lengths are 1.415 and 1.528 Å, respectively.35,37 Therefore, we take lbond = 1.453 Å, the average of the backbone bond lengths of PEO. We generate chain conformations having Na identical atoms, where Na corresponds to the number of backbone atoms in PEO. We focus on the SAW3 model, which is detailed in our previous article.34 In essence, a pair of nonbonded atoms that are separated by at least four (backbone) bonds must be separated from each other by at least two unoccupied sites (i.e., sites on the tetrahedral lattice). The resulting minimum distance between nonbonded pairs of atoms is (11/3)1/2 lbond = 2.78 Å, which corresponds, for instance, to the distance of nonbonded atoms that are involved in an H-bond. With the SAW3 model, more realistic chain conformations can be generated since it implicitly accounts also for the steric effects of hydrogen atoms attached to the carbon atoms of the PEO chain, which are explicitly considered in the MD simulations. However, comparing results for the pure SAW3 model with results from the benchmark PEO in more detail, we find that the SAW3 model places too much weight on the repulsive 2389


Article


obvious that not all possible pairs of atoms (e.g., directly bonded atoms) can form nonbonded contacts. However, there are additional pairs that we do not consider. For example, every 1−4 pair (i.e., atoms which are connected through three covalent bonds) would satisfy the condition of nearest-neighbor contact, as would gauche−gauche 1−5 pairs. However, to include these nonbonded interactions would distort the covalent torsion−angle interactions. The 1−4 atom pair interactions belong to bonded energies. These are excluded in most standard molecular force fields. Also the 1−5 interactions are still practically of bonded character. It turns out that 1−6, 1−7, and 1−9 atom pairs of a SAW3 chain model cannot form contacts for the SAW3 model on account of geometric constraints imposed by the diamond lattice. The 1−8 interaction is a more special case. Figure 1 shows a plot of the fraction of contacts between a pair of atoms versus

character of the nonbonded interactions. Therefore, we define a new model, SAW3(pot), by adding a conformation (conf) dependent potential energy to the pure SAW3 model consisting of two terms, a contact term and a torsional term. The contact potential reads Na − 9

Ucontact(conf ) = εnn

Na

∑ ∑

u[rconf (i , j)], εnn < 0

i=1 j=i+9

(2a)

where ⎧ 19/3 lbond ⎪1, rconf (i , j) ≤ u[rconf (i , j)] = ⎨ ⎪ otherwise ⎩ 0,

(2b)

It is an attractive interaction, which is constant over the range of the “contact” distance [(19/3)1/2 − (11/3)1/2 ] lbond = 0.60 lbond = 0.87 Å for pairs of nonbonded atoms that are separated by nine or more covalent bonds, as defined above. The contact potential, eq 2a, is intended to account for the hydrophobic effect (HYEF). It describes the tendency of water to minimize its contact surface area with a solute molecule such as a polymer chain. As a result, the entropy of water remains high, maximizing the entropy of the whole system, although the chain entropy is reduced. The net effect is a reduction of the chain’s effective surface area, which can be described by attractions between nonbonded atoms of the polymer chain according to eq (2). The HYEF generates short-range attractive interactions between molecular components of the solute. Since the evaluation of the chain surface is computationally more expensive, it can best be described by a contact energy term. Pairs of nonbonded atoms that are at a distance of (19/3)1/2 lbond or closer (eq 2b) are “nearest neighbors” and are said to be in contact. For the SAW3 model, there are three types of such contacts corresponding to different allowed nearest neighbor distances. In two cases, contacting atom pairs are separated by two empty lattice sites. Together with the two empty lattice sites this pair of atoms defines a torsion angle, which can be either in gauche or trans conformation corresponding to distances of (11/3)1/2 lbond and (19/3)1/2 lbond, respectively. In the third type, the contacting atoms are separated by three empty lattice sites with a distance of (16/3)1/2 lbond. Every contacting pair of atoms in a given chain conformation (conf) lowers the potential energy by εnn. Hydrogen atoms on the backbone of the PEO chain tend to favor trans conformations by steric hindrance.36 Since hydrogen atoms are neglected, a second contribution to the potential energy is introduced that promotes a proper distribution between trans and gauche conformations Utorsion(conf ) = εtorPG(conf )Ntor , εtor > 0

Figure 1. Fraction of atoms that are in contact with other atoms in the SAW3 chain versus the number of bonds between the atoms in contact. The histogram results from an ensemble average of chain conformations of Na = 99 atoms: black bars, SAW3; gray bars, SAW3(pot). For Nb = Na − 1 > 7 the heights of the bars shown in the histogram are approximately independent of the chain length for chains with Na = 50 atoms or more. The alternating height of the bars is caused by the two interpenetrating fcc lattices forming the tetrahedral lattice.

the number of bonds that separate the pair. The contacts between 1 and 8 atom pairs occur far more frequently than most other types of contacts, except 1−5 pairs. Contact between 1 and 8 pairs is possible only with a unique sequence of torsion angles (Figure 2) forming a hairpin-like structure. The central part of this loop involving the atoms 2−7 forms a gauche−trans−gauche (GTG) conformation. Such conformations occur undetectably low according to Raman spectroscopy data.38 An attractive contact potential between 1 and 8 atom pairs would enhance the occurrence of this GTG conformer, which is in conflict with the measurements. Therefore, the value of the 1−8 pair contact potential is set to zero in eq 2a. Obviously a polymer model on the diamond lattice may foster for short chain segments neighborhoods that differ from a corresponding polymer model in continuum space. Hence, interactions of atom pairs that are separated only by few covalent bonds may differ for lattice and continuum space polymer models.

(3)

The two potential energy terms (eqs 2 and 3) are combined to give the total potential energy Upot(conf ) = Ucontact(conf ) + Utorsion(conf )

(4)

The quantity PG(conf) in eq 3 is the fraction of local gauche conformations (G) for a given global chain conformation conf. Here, Ntor is the number of torsion angles along the chain backbone (Ntor = Na − 3). Since the potential energy of the chain increases with the fraction of gauche torsions, conformations with trans torsions are preferred. The contact potential defined in eq 2 applies only to pairs of atoms that are separated by nine or more covalent bonds. It is 2390


Article


where ⟨F⟩n refers to the average based on the nth ensemble, which is defined by eqs 5 and 6 for the MD and MC simulations, respectively. The measure of the statistical error is then taken to be ⎡ 1 δF = ⎢ 2 ⎢⎣ nens

⎤1/2

nens

2⎥

∑ (⟨F ⟩n − ⟨F ⟩all ) ⎥ ⎦

n=1

This formula assumes that the ⟨F⟩n values are normally distributed around the overall average ⟨F⟩all. The square of nens appearing in the denominator in eq 9 accounts for the fact that the statistical error is estimated with ensembles whose size is 1/nens of the size if the total ensemble used to compute ⟨F⟩all. Statistical Independence of Conformations from MD Simulation. In order to estimate the number of statistically independent conformations in the MD trajectories, we examine the time autocorrelation function (TAF) of several properties. These are defined by

Figure 2. Unique sequence of torsion angles (blue portion of the chain) between 1 and 8 contact atom pairs at a distance of (19/3)1/2 lbond = 2.516 lbond. The 1−8 contact is indicated by a dashed line. The central part of this hairpin-like conformation involves a gauche−trans− gauche (GTG) conformation.

■

⟨F(k Δt )F(0)⟩MD =

PROPERTIES OF THE PEO CHAIN The potential energy parameters εnn and εtor of the SAW3(pot) model are determined by fitting the following properties to those of the PEO benchmark from MD simulation: Flory radius (root-mean-square end-to-end distance), RF; radius of gyration, Rg; and fraction of gauche and trans conformations, PG(conf) and PT(conf). We also compute the absolute torsional entropy, the SASA, the angular correlation of bond vectors, and the pair distribution function, as outlined below. Ensemble Averages for MD and MC Simulations. The ensemble of NMD MD-generated chain conformations is denoted by the set {conf MD(k), k = 1, 2, . . . NMD}, where conf MD(k) signifies the conformation of the Na backbone atoms at time tk. The ensemble average of a property F[conf MD(k)] is given by ⟨F ⟩MD =

1 NMD

k=1

×

k=1

(10)

CF(k Δt ) = ⟨F(k Δt )F(0)⟩MD /⟨F 2⟩MD

(MC,0) Stor = kB ln(3Ntor ) = NtorkB ln(3)

(5)

N

(MC,3) Stor = kB ln(2.39 Ntor) = NtorkB0.87

1 nens

(MC,3,pot) (MC,3,pot) (MC,3) Stor = δStor + Stor = NtorkBm(pot)

n=1

(14)

(pot)

where m < 0.87 is the presumed slope for the SAW3(pot) model. Here, the superscripts (MC, 3) and (MC, 3, pot) refer to the SAW3 and the SAW3(pot) model, respectively, where the latter includes the influence of the potential energy Upot. The additional entropy decrease δS(MC,3,pot) < 0, due to Upot, is tor computed according to

(7)

(MC,3,pot) (MC,3,pot) (MC,3) δStor = Stor − Stor

= kB ln⟨e−Upot(conf )/ kBT ⟩MC,3

nens

∑ ⟨F ⟩n

(13)

Including also the potential energy Upot (eq 4), we expect the behavior

(6)

is the probability of occurrence of the MC conformation conf MC(k), where Upot is the total potential energy (eq 4). Estimation of Statistical Errors. In general, we construct several ensembles for each system, as described in the next section. We estimate the statistical errors associated with properties computed from nens = 10 ensembles as follows. The overall average value of the property F is given by

⟨F ⟩all =

(12)

With the self-avoidance of the SAW3 model, the entropy increases more slowly with Ntor, as was demonstrated in ref 34, according to

exp(−Upot[confMC (k)]/kBT ) ∑l =MC exp(−Upot[confMC (l)]/kBT ) 1

(11)

to yield CF(0) = 1. Torsional Entropy of a Chain Molecule. According to our previous work,34 the torsional entropy of a free (nonselfavoiding) linear chain on the tetrahedral lattice with Ntor torsional angles generated by a forward random walk (FRW) is

where the index MC indicates that the expression applies to the lattice model, and p[confMC (k)] =

[F(t j + k Δt ) − ⟨F ⟩MD ][F(t j) − ⟨F ⟩MD ]

where the summation runs over time origins t = 0, and kΔt is the displacement from the origin. The TAF in eq 10 is normalized according to

NMC

∑ p[confMC (k)]F[confMC (k)]

∑ j=0

where the index MD indicates that the expression applies to the MD simulation. The corresponding formula for the lattice model is ⟨F ⟩MC =

NMD

1 −k−1

NMD − k − 1

NMD

∑ F[confMD (k)]

(9)

+ (8) 2391

−Upot(conf )/ kBT ⟩MC,3 1 ⟨Upot(conf )e − U conf k T ( )/ B T ⟨e pot ⟩MC,3

(15)


Article

Journal of Chemical Theory and Computation where ⟨. . . . .⟩MC,3 denotes the ensemble average for the SAW3 model. Characteristic Radii of Chain Conformations. The square of the Flory radius RF2 = ⟨Ree2⟩ is defined as the ensemble average of the squared end-to-end distance Ree between atoms 1 and Na of the chain backbone R ee2 = ( rN⃗ a − r1⃗)2

atom as a reference. The atom pair distribution function is then defined by g (r ) = n(r )/V (r )

where r is the distance of a given atom of the chain from the center (reference) atom, n(r) is the number of atoms in a range Δr around r, and V(r) is the volume of the spherical shell of radius r and thickness Δr centered on the reference atom. The distribution function g(r) is computed by assigning the atoms of a given conformation to the proper bin (r − Δr/2, r + Δr/2) and averaging over conformations according to eqs 5 and 6. Angular Correlation of Bond Vectors. To complete the analysis of the three-dimensional polymer structure, we define the correlation of the angles between bonds along the backbone. This property is closely related to the persistence length, which can be physically interpreted in terms of local angular correlation for short “chemical” distances along the backbone.40 Let bi⃗ ≡ ( ri ⃗ − ri −⃗ 1)/| ri ⃗ − ri −⃗ 1| be the unit vector parallel to the bond between atoms i and i − 1 and s = |i − j| be defined as the “chemical” distance between the ith and jth bonds. Analogously to ref 41, we define

(16)

We define the center of geometry (cog) of the chain with Na backbone atoms, whose coordinates are ri ⃗ by R⃗ cog =

1 Na

Na

∑ ri ⃗

(17)

i=1

as well as the following moments of order pow by R gpow = [

1 Na

Na

∑ ( ri ⃗ − R⃗cog)2 ]pow/2 i=1

(18)

where i runs over all atoms of the chain. Then the radius of gyration is

R̅ g = ⟨R g2⟩1/2

⟨cos θ(s)⟩ = ⟨bi⃗ ·bj⃗ ⟩, s = |i − j|

(19)

1 Ntor

PT(conf ) =

1 Ntor

Na − 3

ΔSASA ≈ −0.19608

Na

∑ ∑

Aij

i=1 j=i+3

(23)

where Aij is the surface area lost by the overlap of the spheres centered on atoms i and j. The coefficient in front of the double sum in eq 23 accounts for the effects due to the bonding of tetrahedral carbon atoms to two other non-hydrogen atoms.23 The sphere diameter is chosen to be (19/3)1/2 lbond, which is the largest contact distance.

Ntor

∑ g(φi(conf )) i=1

⎧ −120o ≤ φi(conf ) < 120o ⎪ 1, g (φi(conf )) = ⎨ ⎪ otherwise ⎩ 0,

(22)

The ensemble average of the cosine is examined for short chemical distances s = 1,...,5. Computation of Solvent Accessible Surface Area. We estimate the difference in SASA between a given conformation of the chain and the conformation where all torsion angles are trans, for which the SASA is maximal. We estimate the change in the SASA using the following approximate version of an expression derived in ref 42

Note that the above definition of R̅ g referring to the center of geometry differs only marginally from one referring to the center of gravity since the effective masses of the backbone atoms in the PEO chain differ only slightly [CH2 (14 amu) and O (16 amu)]. Fraction of Gauche and Trans Torsion Angles. We take gauche and antigauche conformations to correspond respectively to torsion angles in the ranges −120° ≤ φ < 0° and 0° ≤ φ < 120°, respectively. Torsion angles that correspond to trans conformations lie in the complementary ranges −180° ≤ φ < − 120° and 120° ≤ φ < 180°. We define the fraction of gauche and antigauche, PG(conf), as well as the fraction of trans torsion angles, PT(conf), of a single conformation conf by PG(conf ) =

(21)

■

COMPUTER SIMULATIONS Molecular Dynamics Simulations. The molecular system, comprising a single PEO chain and a large number of water molecules, is assembled by means of the program CHARMM36b1.43 We use the CHARMM ether force field C35r developed by Vorobyov et al.35 and improved by Lee et al.36 The edge lengths of the cubes used for imposition of periodic boundary conditions are set at 35, 50, 65, and 75 Å for chains containing 27, 51, 75, and 99 backbone atoms, respectively. MD simulations were performed using the NAMD44 software on graphics processing units (GeForce GTX 670, GeForce GTX Titan, and GeForce 780 Ti). The temperature of the system is held constant at 300 K by Langevin dynamics with a friction constant of β = 1 ps−1 and the pressure at 1 atm by the Nosé−Hoover Langevin piston method.45,46 The electrostatic interactions are evaluated explicitly with dielectric constant of unity and cutoff distance of 12 Å. The van der Waals interactions are cut off by means of a switching function in the interval from 10 to13.5 Å. The time step of the MD simulation is 2 fs. The nonbonded interactions

(20a)

Ntor

∑ t(φi(conf )) i=1

⎧ −180o ≤ φi(conf ) < −120oor ⎪ ⎪1, t(φi(conf )) = ⎨ 120o ≤ φi(conf ) < 180o ⎪ ⎪ 0, otherwise ⎩ (20b)

where Ntor = Na − 3 is the number of torsion angles in the polymer chain. Here, φi(conf) denotes the value of the ith torsion angle of conformation conf. Atom Pair Distribution Function. To obtain a more detailed picture of the three-dimensional structure of the polymer chain, we define an atom pair distribution function that shows the degree of local order.39 Since the backbone atoms of a polymer chain are not equivalent, we take the center 2392


Article


7). The averages are computed “on the fly” (i.e., no conformation is stored). For chain lengths of 27, 51, 75, 99, 120, 150, 193, and 257 backbone atoms, we generated 10 ensembles. Each ensemble consists of 108 independent conformations of the SAW3 model. For the longest chain (321 atoms), 30 ensembles, each one consisting of 108 conformations, were used in order to mitigate the statistical error, which increases with the chain length. Depending on the chain length, different algorithms, or adjustments, were employed to maximize the efficiency. A summary of the algorithms can be found in Table SI of the Supporting Information. Optimization of Lattice Model. The SAW3(pot) model is optimized by adjusting the energy parameters εnn and εtor so that RF, Rg, and the mean fraction of torsion angles, ⟨PT(conf)⟩, agree with the corresponding properties of the PEO benchmark. We define a deviation function by

are evaluated every time step, and a complete electrostatic energy evaluation is performed every second time step. The Verlet neighbor list is updated every 10 time steps. All bond lengths between hydrogen and other atoms are fixed by the SHAKE algorithm. The atomic coordinates of the molecular system are saved every 100 ps. The initial PEO chain conformation is generated with CHARMM43as follows: (1) The PEO chain is assembled in the all-trans conformation with the monomers and end groups described by the ether topology files of the C35r ether force field.35,36 (2) To obtain a more densely packed initial conformation, gauche torsion angles are introduced (i.e., for chains with 27, 51, 75, and 99 atoms 3, 4, 8, and 9 trans torsion angles were changed to gauche, respectively). (3) The energy of the chain conformation is minimized by application of 20,000 steps of the steepest descent algorithm, followed by 10,000 steps of the ABNR algorithm. (4) The chain is immersed in the center of a box with periodic boundary conditions containing pre-equilibrated water. Water molecules that overlap with polymer atoms, as determined from the atomic van der Waals radii, are removed. (5) The configurational energy of the resulting molecular system of the PEO polymer chain plus water is minimized by 3000 steps in NAMD.44 (6) The following MD simulation is performed with NAMD, where the initial velocities are randomly assigned to the atoms according to the Boltzmann distribution at absolute temperature of 20 K. The temperature of the system is subsequently raised to 300 K in steps of 20 K by velocity rescaling over a time period of 50 ps. Simultaneously, Langevin dynamics with a small friction constant of 1 ps−1 at 300 K is applied. The resulting equilibrated molecular system is used for “production” MD simulations. The MD simulation of each of the 9-, 17-, and 25-monomer PEO chains was is performed as a single, continuous run starting from the initial state, prepared as described above. Because of the large number of atoms (40,988) in the 33monomer PEO chain, the MD simulation is too slow to be performed as a single trajectory. Therefore, parallel simulations are carried out as follows. At 300 ns the conformation of the initial trajectory is taken as the starting point for two new trajectories, while the original trajectory is continued to 824 ns. For the two new trajectories, the initial velocities are randomized at a temperature of 20 K. The subsequent heating and equilibration is performed in the same way as for the initial trajectory. These two additional trajectories are propagated for 100 and 228 ns. Ensemble averages of the PEO chain of 33 monomers are based on all three trajectories. The first nanosecond of each production run is discarded in order to avoid a bias from the initial conformation. The trajectories are partitioned into 10 subsets of roughly equal size (i.e., ⟨F⟩ is the average of the nens = 10 “subensemble” averages; eq 8). Thereby the three trajectories of the largest chain are concatenated before they are partitioned. Monte Carlo Lattice Simulations. Basic Method. The polymer chains are grown efficiently by means of fragmentfusion and brick-fusion algorithms, as detailed in our previous article.34 We do not use such techniques as “pivoting”47 or “sample enrichment”,48 which manipulate pre-existing conformations. The energy of a given conformation is computed according to the potential energy U pot (eq 4). The thermodynamic average ⟨F⟩ of a property F is evaluated with the corresponding Boltzmann factors, according to eqs 6 and

Δ(εnn , εtor) =

∑ |1 − ⟨F ⟩MC,L (εnn , εtor)/⟨F ⟩MD,L | F ,L

(24)

where the summation on F runs over the three considered properties [RF, Rg, ⟨PT(conf)⟩] and that on L runs over the three longer of the four considered PEO chains involving 51, 75, and 99 backbone atoms for which MD simulations were performed. We exclude the shortest PEO chain (27 backbone atoms) in order to avoid anomalous influences from the end atoms, which are dominant for too short chains. Note that ⟨F⟩MC,L is a function of the energy parameters εnn and εtor. We first determine the minimum of Δ(εnn, εtor) on a relatively crude 10 × 10 grid in the (εnn, εtor)-parameter space. Next, we center a refined 10 × 10 grid on the minimum and find a new minimum. This procedure, which is repeated until Δ(εnn, εtor) is less than the sum of relative errors of the considered properties (eq 9), yields εnn = −0.44 kBT = −1097 J/mol and εtor = 0.49 kBT = 1222 J/mol. If we consider also the shortest 27 atom chain in the optimization technique based on eq 24, we obtain the parameters εnn = −0.42 kBT and εtor = 0.40 kBT.

■

RESULTS AND DISCUSSION Comparison of Efficiencies of MD and MC Simulations. In order to determine the optimum time lag between chain conformations of the MD trajectory (i.e., to estimate the number of statistically independent conformations), we examine the time autocorrelation function (TAF) (eq 10) of the end-to-end distance Ree (eq 16) and the square of the radius of gyration, Rg2 (eq 19). Because the TAF of Ree is found to decay faster than that of Rg2, the former is used to estimate the TAFs for the four PEO chains considered (i.e., 9, 17, 25, and 33 monomers). Figure 3 displays plots of the TAFs. The correlation time, τcorr, is defined as the time required for the TAF to drop to e−1. The number of statistically independent conformations, Neff MD, is then taken to be the ratio of the trajectory time length, ttraj, to the correlation time, τcorr, i.e., eff NMD = t traj/τcorr

(25)

The numbers of independent conformations of PEO chains of length 27, 51, 75, and 99 are found to be 10,000, 2758, 862, and 968 for MD simulation times (ttraj) of 1, 0.8, 0.63, and 1.17 μs, respectively (Figure 3). The efficiency of the MD simulation (Ef f MD) is defined as the ratio of number of independent conformations (Neff MD) per unit wall time, 2393


Article


Figure 3. Time autocorrelation function of end-to-end distance (Ree, eq 16) versus time based on MD simulations of the PEO chains comprising 27, 51, 75, and 99 backbone atoms. The horizontal line indicates 1/e level. Vertical lines indicate estimates of correlation times, τcorr: 0.10 ns (27), 0.29 ns (51), 0.73 ns (75), 1.21 ns (99). eff EffMD = NMD /twall(MD)

Figure 4. Comparison of efficiency of MD (⊙) (eqs 25 and 26) and MC (×) (eq 28) techniques versus chain length. The MD data are all scaled to the fastest GPU (GeForce 780 Ti) used for the MD simulations. The efficiency does not necessarily decrease monotonously with system size, as explained in the text.

may therefore appear in the MC efficiency for longer chain lengths (Figure 4). Torsional Entropy. Figure 5 shows the torsional entropy of different chain models as a function of the number of torsion

(26)

The wall times of the MD trajectories generated on different GPUs were scaled to the graphics card GeForce 780 Ti, the fastest GPU used in the present study. The conformations of the SAW3 chains on the tetrahedral lattice are generated independently and are therefore uncorrelated. The chain conformations obtained by MC on the lattice and by MD in continuous space are assumed to represent the same system with the same underlying distribution. Since the conformations of the SAW3(pot) model are weighted differently than those of the pure SAW3 model (in order to account for the potential energy in the former), the effective number of conformations in the ensemble for the SAW3(pot) model (Neff MC) differs from the number of conformations in the ensemble for the SAW3 model (NMC). To estimate the effective number of independent conformations for the SAW3(pot) ensemble, we invoke the connection between entropy and ensemble size given by the Boltzmann-equation S = kB ln N

Figure 5. Torsional entropy versus number of torsion angles Ntor = Na − 3 in the chain for the lattice model. FRW (dashed line) nonselfavoiding chain generated as forward random walk, slope m = 1.10, eq 12. SAW3 (solid black line) chain model generated by self-avoiding random walk, slope m = 0.87, eq 13. SAW3(pot), (solid red line): SAW3, chain model that includes the potential energy, eq 4, slope m(pot) = 0.792 ± 0.008, eq 14.

(27)

Thus, the effective number Neff MC of SAW3(pot) conforma(MC,3,pot) tions is estimated from the entropy difference δStor according to eff (MC,3,pot) NMC = NMC exp(δStor /kB)

(28)

angles (Ntor). The total torsional entropy S(MC,3,pot) (eq 14) tor consists of two terms, S(MC,3) (eq 13) for the pure SAW3 model tor and δS(MC,3,pot) = S(MC,3,pot) − S(MC,3) (eq 15), the correction due tor tor tor to the potential energy Upot (eq 4). The nonself-avoiding chain, which is equivalent to the forward random walk (FRW) model, exhibits the greatest slope, m = ln (3) = 1.0986 as a function of the number of torsion angles, Ntor. It describes the increase in entropy as the chain grows by one atom in the absence of interactions between nonbonded atoms, whose presence would diminish the increase of entropy. For the SAW3 model, the slope (m = 0.87) is smaller.34 Inclusion of the potential energy (eq 4) reduces the slope further to m(pot) = 0.79. Hence, the entropy is greatest for the FRW and drops in the order FRW > SAW3 > SAW3(pot). To connect the torsional entropy of the PEO chain model on the tetrahedral lattice with the corresponding chain model in continuous space used in MD simulations, one must add a correction, which according to ref 34 gives

The values of the parameters used in eq 28 are listed in Table SII of the Supporting Information. On the basis of NMC = 109 independent chain conformations of the SAW3 ensemble, eq 28 yields about 1.07 × 108, 2.17 × 108, and 3.61 × 105, 9491 independent conformations for the chains of 27, 55, 75, and 99 atoms, respectively. In analogy to eq 26, the efficiency of the MC method is defined as Eff MC = Neff MC/twall(MC). The efficiency of both simulation techniques (MD and MC) is compared in Figure 4. The comparison favors the MD method since the MC method used just one CPU core (Intel Xeon L5420), whereas the MD simulations were run on a fast GPU (GeForce 780 Ti). Nevertheless, it is obvious that MC is much more efficient than MD. Equation 28 provides a measure for the effective number of SAW3(pot) conformations but may involve significant errors, since δS(MC,3,pot) appears in the exponential function. Jumps tor 2394


Article


where Nb = Na − 1 is the number backbone bonds. The shortest chain (Nb = 26) is not considered in this fit since it is anomalously stiff because it is so short. Note the excellent agreement of the optimized SAW3(pot) model with the PEO benchmark. It is particularly noteworthy that the lattice model results lie very close to the extension of the line determined by fitting only the MD results, which pertain to relatively short chains. Hence, for the time being, the lattice model can be employed in ranges that are well beyond the reach of practical MD simulations. According to Flory, the power law is directly connected to the fractal dimension Df of the polymer chain, which describes the scaling behavior of the effective volume of the chain.49 For nonself-avoiding chains Df = 2, whereas for self-avoiding chains, like all SAWn models presented in ref 34, Df = 1.7006.1 Previous explicit solvent MD simulations for PEO find Df = 236 based on a trajectory much shorter than in the present study. Note that this finding (i.e., Df = 2) conflicts with the experimentally determined self-avoiding scaling behavior found for long PEO chains with several thousand backbone atoms.50 Even though much longer simulation times are used in the present study to determine Df for the PEO chain, the result is in accordance with that of Lee et al.36 The main difference between the MD simulations and the experimental results is the chain length. Devan and Selser50 considered PEO polymer chains that are 2 orders of magnitude longer than those employed in the present study. Polymer chains of this length are simply beyond the scope of MD simulations with explicit water. The different behavior of short and long PEO chains may be explained by a consideration of the nature of the HYEF. In aqueous solution, the presence of solute molecules causes water molecules to lose entropy. Water molecules in contact with solute molecules are forced into a fixed orientation and have thus a reduced ability to assume alternative conformations with different hydrogen-bonding patterns involving other water molecules. In order to minimize the entropy loss of the water molecules caused by the presence of solute molecules, the system reacts by reducing the SASA of solute molecules. As a consequence, polymer chains form loops and more selfcontacts in order to reduce their SASA and thus adopt structures that are more compact (and effectively of higher dimension) than those generated by self-avoiding-walk polymer models. This scenario apparently applies to short polymer chains. For the lattice model of the polymer chain, this HYEF is accounted for by the attractive contact energy (eq 2). Much longer polymer chains occupy a volume much larger than that of a densely packed polymer with formal dimension of Df = 3. Hence, the long polymer chains necessarily form a sponge-like structure and contain a large number of small cavities. To increase the overall entropy of the polymer chain water system, these cavities need to be filled with water molecules, even though these molecules lose entropy. The expansion of these cavities that results from their taking on more water means that the additional water molecules are no longer in contact with the polymer and therefore behave bulklike. This process facilitates a further swelling extension of the polymer chain structure up to the lower dimension typical of a self-avoiding chain.51 In this way, large polymer chains assume a globular structure and thus create their own environment in their interior, which vitiates the HYEF normally created by bulk water. The phenomenon relates to the swelling of polymers when they are solvated.49,52,53 The polymer chains considered in the MD simulations are too short to form such cavities and

(MD) (MC,3,pot) Stor = Stor + kBNtor[ln(C /Λ) + 1/2 − ln 3]

(29)

where C is the circumference of the circle of radius lbond sin(γbond) on which the atom moves as the torsion angle varies, Λ ≡ h/(2πmkBT)1/2 is the thermal de Broglie wavelength, m is the mass of the atom, h is Planck’s constant, kB is Boltzmann’s constant, and T is the absolute temperature. Flory Radius and Radius of Gyration. Figures 6 and 7 show log−log plots of RF and Rg as functions of chain length,

Figure 6. Flory radius RF versus number of bonds, Nb = Na − 1, in chain as a double logarithmic plot. Circles refer to MD simulation of PEO and crosses to the SAW3(pot) lattice model optimized for PEO. The power law (eq 30) (solid red line, Df = 2.08 ± 0.15) is fitted to the MD results for the three longest PEO chains with 51, 75, and 99 backbone atoms. For comparison, results of the forward random walk (FRW, dashed line, Df = 2) and SAW3 (solid black line, Df = 1.7) are also shown.

Figure 7. Radius of gyration Rg versus number of bonds, Nb = Na − 1, in chain in a double logarithmic plot. Circles refer to MD simulation of PEO and crosses to the optimized SAW3(pot) model. The power law (eq 30) (solid red line, Df = 1.95 ± 0.15) is fitted to the results of MD simulations for three longest PEO chains with 51, 75, and 99 backbone atoms. For comparison, results of the forward random walk (FRW, dashed line, Df = 2) and SAW3 (solid black line, Df = 1.7) are also shown.

Na. A simple computation shows that for freely jointed chains the two radii scale as RF2 = 6 Rg2. This is evidently also true of the chain models considered in Figures 6 and 7. The MD results for PEO chains with 51, 75, and 99 backbone atoms can be fitted by the following power law

RF ∝ R̅ g ∝ Nb1/ Df

(30) 2395


Article

Journal of Chemical Theory and Computation are therefore under the influence of the HYEF of bulk water. Hence, they are more compact and have a larger fractal dimension. Distribution of Radius of Gyration. Figure 8 displays the quantity Mg4 = ⟨Rg4⟩/⟨Rg2⟩2, which is the reduced fourth-order

Figure 9. Mean fraction of trans torsion angles, PT, versus number of torsion angles. Circles refer to MD for PEO chains with 27, 51, 75, and 99 atoms and crosses to MC for optimized SAW3(pot) model. Horizontal (red) line at 0.58 indicates mean PT from MD. The statistical error for all data points is less than 0.0035. Figure 8. The reduced moment Mg4 = ⟨Rg4⟩/⟨Rg2⟩2 versus the inverse of number of bonds, 1/Nb. The averages are defined by eqs 17 and 18. The black solid line refers to PEO benchmark and the red line to the optimized SAW3(pot) model. For both models, the error bars are displayed. For the two shorter chains, the error bars for SAW3(pot) model are smaller than the line width. This is also the case for all other data shown in the figure, which refer to both self-avoiding (SAW1 and SAW3) and nonself-avoiding models [random walk (RW) and forward random walk (FRW)].

approximately independent of the length of the chain. Again, the optimized SAW3(pot) model agrees well with the MD benchmark as well as with spectroscopic measurements on dimethoxyethane in solution, where PT ≈ 0.563.38 Atom Pair Distribution Function. In Figure 10, the atom pair distribution g(r) relative to the center atom 50 in a

moment as a function of the inverse chain length 1/Nb. It characterizes the shape of the distribution of Rg. For a nonselfavoiding chain on a simple cubic lattice, this ratio is 19/15 = 1.267 in the limit of an infinitely long chain (see eq 8.44 of ref 54). The same value is obtained for a chain on the tetrahedral lattice in the present study. Hence, the value of Mg4 is independent of the underlying grid. For the Flory radius, the ratio ⟨RF4⟩/⟨RF2⟩2 is 5/3 = 1.667, which is equal to the ratio for the random Gaussian chain model.55 For a self-avoiding chain, the ratio ⟨RF4⟩/⟨RF2⟩2 = 1.513,56 slightly smaller than for the random Gaussian chain. For the self-avoiding chain, SAW1, on a tetrahedral lattice, we estimated by extrapolation using a maximum chain length of Nb = 300 the value 1.5, which is close to the theoretical limit. For short chains, Mg4 moment is systematically underestimated. It seems that the chains need to be sufficiently long to represent conformations at the boundaries appropriately. In the absence of self-avoidance the plot whose abscissa is 1/Nb seems most suitable for extrapolation to long chains. With an increasing degree of self-avoidance, Mg4 becomes smaller. This can be explained as follows. Self-avoidance broadens the distribution of radii primarily around the center of the chain, whereas it has little effect at large distances from the chain center. It is clear from Figure 8 that this trend is inverted if attractive interactions are added to the SAW3 model to give the SAW3(pot) model. The ⟨Rg4⟩/⟨Rg2⟩2 from the MD simulation exhibits a dependence on chain length similar to that of the SAW3(pot) model, but it is uniformly lower (Figure 8). This may be due to subtle influences of weak attractive forces at longer distances, which are absent from the simple SAW3(pot) model. Fraction of Trans Torsion Angles. Figure 9 displays a plot of the mean fraction of trans torsion angles, PT, for the SAW3(pot) model and the MD benchmark versus the number of torsion angles in the chain, Ntor = Na − 3. Here, PT is

Figure 10. Comparison of the pair correlation function g(r) defined in eq 21 for different lattice models and the MD benchmark on a semilogarithmic plot. The atom pair distribution function is shown as a function of atom pair distances relative to the center atom of the 99 atom chain. The figure was generated with bin thickness of Δr = 1.0 Å. The errors of FRW, SAW1, and SAW3 are smaller than the thickness of the lines. For the benchmark (solid black line) and the SAW3(pot) model (solid red line), the error bars are indicated.

polymer chain of 99 backbone atoms is shown for the various models in a semi-logarithmic plot. The atom pair distribution functions for chains with 27, 51, and 75 backbone atoms can be found in the Supporting Information (Figure SI−SIII). All atom pair distributions decay exponentially initially and then more rapidly as r becomes large. The chains with the most compact conformations, smallest Flory radius, and radius of gyration are, of course, manifested by FRW model, for which more than one atom of the chain is allowed to occupy the same lattice site. These chains show fractal dimension Df = 2. Thus, g(r) for the FRW decays most rapidly (Figure 10). The least compact chains belong to the SAW3 ensemble with the fractal dimension Df = 1.7.1 The g(r) for SAW3 decays least rapidly of all models. The atom pair distribution function g(r) of the SAW3(pot) model decays only slightly faster than that of the SAW3 model 2396


Article


MD simulations, we compute the difference in SASA between a given conformation and the all-trans conformer according to eq 23. The expected loss of SASA is related to the decrease in the contact potential energy, eq 2, by

but still much slower than that of the FRW model. This is surprising since the fractal dimensions for both the SAW3(pot) and FRW models are practically identical and are larger than that of the SAW3 model. However, it is even more surprising that g(r) of the SAW1 model decays much faster than that of the SAW3(pot) model, although the fractal dimension of the former is smaller than that of the latter. Hence, the fractal dimensions suggest that SAW1 chains are more stretched than SAW3(pot) chains, while the opposite would be concluded from the atom pair distribution function. The apparent contradiction can be resolved by observing that the fractal dimension of a chain emphasizes the behavior of chains at large distances, while the atom pair distribution function emphasizes the behavior at short distances. The atom pair distribution function of the SAW3(pot) model is virtually identical with the g(r) obtained for the benchmark, although we did not account for the g(r) in optimizing the potential parameters of the SAW3(pot) model. The deviation between the g(r) of the SAW3(pot) model from the MD simulation data at the distance of 5 Å is due to the discrete character of the lattice, which appears at small r (Figure 10). Angular Correlation. In Figure 11, the angular correlation defined in eq 22 as a function of the chemical distance s in a

Na − 9

εs = εnn⟨ ∑

Na

∑

u[rij(conf )]⟩⟨ΔSASA (conf )⟩−1 (31)

i=1 j=i+9

The idea is that the contact potential of the SAW3(pot) model could, in principle, be replaced by a surface potential with the energy parameter εs. The values of both ensemble averages in eq 31 depend linearly on the number of atoms. Thus, the surface energy parameter εs is approximated by the contact parameter εnn = −1097 J/mol times the mean number of contacts (0.0151 Na) divided by the mean loss of solvent accessible surface area (−1.1195Na Å2) εs =

−1097 J mol−1 (0.0151 Na) −1.120 Na Å2

≈ 3.5cal/(mol Å2)

which correlates with other estimated values of this parameter, as explained below. In a recent adjustment of force fields for implicit solvent MD simulations, the value of εs for aliphatic carbon atoms is taken to be about 9.5 cal/(mol Å2).57 However, the interquartile range of the εs values derived from MD simulations indicates that individual values of εs can vary considerably. In much earlier work, measured enthalpies of solvation of amino acids were corrected for entropy contributions and transformed on an atomic basis to yield 4 or 12 cal/(mol Å2), depending on the applied entropy correction.58,59 Thus, the uncertainty in εs is quite large. The optimized surface energy parameter εs = 3.5 cal/(mol Å2) for the SAW3(pot) model applies equally to all backbone atoms. For the ether oxygen present in PEO εs is unavailable. In ref 57, the values of εs are estimated to be negative for oxygen atoms. Although no value is provided for ether oxygen, one can expect that εs is also negative in this case. Hence, the low value of εs for the SAW3(pot) model seems reasonable. Thus, PEO seems to be only slightly hydrophobic so that it is still quite soluble in water. Figure 12 shows the loss of SASA relative to the all trans conformer computed according to eq 23 as a function of the chain length. The SASA loss increases in the following order:

Figure 11. Comparison of angular correlation ⟨cos θ(s)⟩ versus chemical distance s (on a semilogarithmic plot), according to eq 22, for different lattice models with the PEO benchmark. The angular correlation is shown as a function of the chemical distance along the backbone of the 99 atom chain. Errors for all polymer models are smaller than the line thickness.

polymer chain of 99 backbone is shown for various polymer models. Plots for the three shorter chains lengths can be found in the Supporting Information (Figure SIV−SVI). The nonselfavoiding polymer model, FRW, decays exponentially with the chemical distance. The rate of decay of ⟨cos θ(s)⟩ for the selfavoiding polymer models decreases with an increasing degree of self-avoidance. The MD benchmark shows a distinct increase for s = 2, which is only matched by the SAW3(pot) model. This indicates that the improved polymer model represents the benchmark also for short distances along the backbone. Thus, SAW3(pot) behaves like the MD benchmark for both short and long chains. Relation of Contact Potential to SASA Model Potential. The contact potential energy Ucontact (eq 2) is introduced to account for the HYEF in the SAW3 model. In MD simulations with implicit solvent, the HYEF is normally modeled by a surface potential, which tends to reduce the effective solvent accessible surface area (SASA). In order to compare the contact potential with the surface potential used in

Figure 12. Comparison of the loss of solvent accessible surface area (SASA) according to eq 23 relative to the all-trans conformation for different lattice models and for the MD benchmark (black). The SAW3(pot) model (red) is the closest to the MD. Errors are smaller than the line thickness. 2397


Article


agrees very well with g(r) for the MD benchmark of PEO, although it was not among the properties used in fitting the SAW3(pot) model potential parameters to the MD data. The comparison of the CPU times for the MD simulation with those for the SAW3(pot) model favor the latter method. The generation of conformations for the longer chains with 75 and 99 backbone atoms is at least 2 orders of magnitude faster for the SAW3(pot) model than for the MD benchmark. The present approach (i.e., optimization of a lattice model using MD data) can be extended to other types of polymers. We expect that the tetrahedral lattice models introduced in the present study will be very useful in the investigation of computationally demanding problems, for example, the study of multivalent ligand binding with polymer linkers. The present SAW3(pot) model may also be useful for describing PEO polymer solutions of low concentration. The SAW3(pot) model is flexible enough to be applicable to different types of chemical polymers and different solvents by varying contact energies and torsion potentials and adjusting them to known distributions of chain conformers. However, in solvents of low dielectric constant, explicit electrostatic interactions among polymer atoms would need to be considered as well. The Monte Carlo program to generate SAW3 model ensembles is freely available for download at https://github. com/Roulattice/Roulattice.

MD benchmark, SAW3(pot), SAW3, SAW1, and FRW chain models. The SASA loss is proportional to the magnitude of the HYEF, which evidently diminishes with increasing selfavoidance from FRW to SAW3. Hence, the SASA loss suggests that the MD benchmark is the least hydrophobic model. It seems counterintuitive that SAW3(pot) is more hydrophilic than SAW3, even though a hydrophobic potential term is used. This contradiction can be resolved by looking at the interplay between contact and torsion energy terms. The latter forces the chain to take on more trans conformations, which increase SASA. Neglecting the torsion term in the SAW3(pot) model increases the SASA loss compared with the SAW3 model. The discrepancy between MD and SAW3(pot) is small and may be due to the approximation implicit in eq 23.

■

SUMMARY The present study demonstrates that lattice polymer models are useful not only for investigations of fundamental phenomena in physics, such as universality class behavior and glass transition, but also for characterization of the physicochemical behavior of specific, chemically relevant polymer chains in atomic detail. This requires, however, a special lattice-supporting realistic bond geometry. For this purpose, the diamond lattice is most appropriate for carbon-based polymers since polymer models whose backbone atoms are on these lattice sites possess the tetrahedral bond angle of 109°, a value which is very close to bond angles of polymers with sp3 hybridization of backbone atoms. Use of the conventional self-avoiding walk (SAW1), however, where each lattice site can be occupied only once, leads to nonbonded atom pair distances that are too small. Thus, the SAW3 polymer model was designed to keep nonbonded atom pairs at larger distances.34 The chains of SAWn models exhibit the fractal dimension Df = 1.7.34 In contrast PEO chains in explicit water, simulated by MD in the present study, are more compact, having a fractal dimension around 2 and a higher fraction of trans conformers of the torsion angles than the SAW3 chains. To modify the SAW3 model on the tetrahedral lattice so that its behavior conforms to that of PEO simulated by MD, we introduce an additional potential energy comprising two contributions: an attractive contact potential to account for the HYEF and a torsion potential, designed to increase the fraction of trans conformers. The potential parameters are determined by adjusting the Flory radius, the radius of gyration and the fraction of trans torsion angles to the MD simulation data. In essence, the PEO is a benchmark for the design of the SAW lattice model that behaves like the PEO. The attractive contact potential corresponds to a weak HYEF of solvation. The strength of contact potential that we find for the PEO benchmark correlates well with the literature values. It is interesting to note that the contact potential, albeit rather weak, changes the universality class behavior of short PEO chains, which now exhibit a fractal dimension of about 2. This is in conflict with measurements, which are however performed for much longer chains.50 We may resolve the discrepancy by considering that large polymers can generate their own environment with cavities in their interior large enough to host water that behaves as if it were bulk. In contrast to the Flory radius and the radius of gyration, whose values decrease significantly for the SAW3(pot) model, the torsional entropy decreases only slowly, and the atom pair distribution function g(r) decays only moderately faster with distance r. It is also notable that g(r) of the SAW3(pot) model

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00144. List of algorithms used for the Monte Carlo chain generation. Table with detailed values of the efficiency calculation. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

* E-mail: [email protected]. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS J.C.B.D. acknowledges his student position financed by the graduate college 1772 “Computational Systems Biology” of the Deutsche Forschungsgemeinschaft (DFG). The initial phase of this work was financially supported by a project in the Sfb765 of the DFG.

■

ABBREVIATIONS FRW, forward random walk; HYEF, hydrophobic effect; MC, Monte Carlo; MD, molecular dynamic; PEO, polyethylene oxide; RIS, rotational isomeric state; RW, random walk; SASA, solvent accessible surface area; SAWn, self-avoiding walk of degree n

■

REFERENCES

(1) LeGuillou, J. C.; Zinn-Justin, J. Critical Exponents for the nVector Model in Three Dimensions from Field Theory. Phys. Rev. Lett. 1977, 39 (2), 95−98.

2398


Article


(23) Karimi-Varzaneh, H. A.; van der Vegt, N. F. A.; Mueller-Plathe, F.; Carbone, P. How Good Are Coarse-Grained Polymer Models? A Comparison for Atactic Polystyrene. ChemPhysChem 2012, 13 (15), 3428−3439. (24) Coarse-Graining of Condensed Phase and Biomolecular Systems; Voth, G. A., Ed.; CRC Press: Boca Raton, FL, 2008. (25) Mogilner, A.; Oster, G. Polymer Motors: Pushing out the Front and Pulling up the Back. Curr. Biol. 2003, 13 (18), R721−R733. (26) Tark-Dame, M.; van Driel, R.; Heermann, D. W. Chromatin folding - from biology to polymer models and back. J. Cell Sci. 2011, 124 (6), 839−845. (27) Binder, K.; Milchev, A. Off-lattice Monte Carlo methods for coarse-grained models of polymeric materials and selected applications. J. Comput.-Aided Mater. Des. 2002, 9, 33−74. (28) Odegard, G. M.; Gates, T. S.; Wise, K. E.; Park, C.; Siochi, E. J. Constitutive modeling of nanotube-reinforced polymer composites. Compos. Sci. Technol. 2003, 63 (11), 1671−1687. (29) Martinez-Veracoechea, F. J.; Frenkel, D. Designing super selectivity in multivalent nano-particle binding. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (27), 10963−10968. (30) Martinez-Veracoechea, F. J.; Leunissen, M. E. The entropic impact of tethering multivalency and dynamic recruitment in systems with specific binding groups. Soft Matter 2013, 9 (12), 3213−3219. (31) Fasting, C.; Schalley, C. A.; Weber, M.; Seitz, O.; Hecht, S.; Koksch, B.; Dernedde, J.; Graf, C.; Knapp, E. W.; Haag, R. Multivalency as a Chemical Organization and Action Principle. Angew. Chem., Int. Ed. 2012, 51 (42), 10472−10498. (32) Numata, J.; Juneja, A.; Diestler, D. J.; Knapp, E. W. Influence of spacer-receptor interactions on the stability of bivalent ligand-receptor complexes. J. Phys. Chem. B 2012, 116 (8), 2595−604. (33) Turrin, C.-O.; Caminade, A.-M. Unexpected Biological Applications of Dendrimers and Specific Multivalency Activities. In Dendrimers: Towards Catalytic, Material and Biomedical Uses; Caminade, A.-M., Turrin, C.-O., Laurent, R., Ouali, A., DelavauxNicot, B., Eds.; John Wiley & Sons, Ltd.: Chichester, U.K., 2011. (34) Dietschreit, J.; Diestler, D. J.; Knapp, E. W. Models for SelfAvoiding Polymer Chains on the Tetrahedral Lattice. Macromol. Theory Simul. 2014, 23 (7), 452−463. (35) Vorobyov, I.; Anisimov, V. M.; Greene, S.; Venable, R. M.; Moser, A.; Pastor, R. W.; MacKerell, A. D. Additive and classical drude polarizable force fields for linear and cyclic ethers. J. Chem. Theory Comput. 2007, 3 (3), 1120−1133. (36) Lee, H.; Venable, R. M.; MacKerell, A. D.; Pastor, R. W. Molecular dynamics studies of polyethylene oxide and polyethylene glycol: Hydrodynamic radius and shape anisotropy. Biophys. J. 2008, 95 (4), 1590−1599. (37) Mackerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Nguyen, D. T.; Ngo, T.; Prodhom, B.; Roux, B.; Schlenkrich, M.; Smith, J.; Stote, R.; Straub, J.; Wiorkiewiczkuczera, J.; Karplus, M. Self-Consistent Parameterization of Biomolecules for Molecular Modeling and Condensed Phase Simulations. Faseb J. 1992, 6 (1), A143. (38) Goutev, N.; Ohno, K.; Matsuura, H. Raman spectroscopic study on the conformation of 1,2-dimethoxyethane in the liquid phase and in aqueous solutions. J. Phys. Chem. A 2000, 104, 9226−9232. (39) Stillinger, F. H.; Weber, T. A. Inherent structure theory of liquids in the hard-sphere limit. J. Chem. Phys. 1985, 83, 4767−4775. (40) Huang, A.; Hsu, H.-P.; Bhattacharya, A.; Binder, K. Semiflexible macromolecules in quasi-one-dimensional confinement: Discrete versus continuous bond angles. J. Chem. Phys. 2015, 143, 243102. (41) Hsu, H.-P.; Paul, W.; Binder, K. Standard Definitions of Persistence Length Do Not Describe the Local ″Intrinsic″ Stiffness of Real Polymer Chains. Macromolecules 2010, 43, 3094−3102. (42) Weiser, J.; Shenkin, P. S.; Still, W. C. Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 1999, 20 (2), 217−230. (43) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch,

(2) Eizenberg, N.; Klafter, J. Critical exponents of self-avoiding walks in three dimensions. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53 (9), 5078−5081. (3) Binder, K.; Baschnagel, J.; Paul, W. Glass transition of polymer melts: test of theoretical concepts, by computer simulation. Prog. Polym. Sci. 2003, 28 (1), 115−172. (4) Rapold, R. F.; Mattice, W. L. New high-coordination lattice model for rotational isomeric state polymer chains. J. Chem. Soc., Faraday Trans. 1995, 91 (16), 2435−2441. (5) Doruker, P.; Mattice, W. L. Reverse mapping of coarse-grained polyethylene chains from the second nearest neighbor diamond lattice to an atomistic model in continuous space. Macromolecules 1997, 30 (18), 5520−5526. (6) Helfer, C. A.; Xu, G.; Mattice, W. L.; Pugh, C. Monte Carlo simulations investigating the threading of cyclic polyethylene oxide by linear chains in the melt. Macromolecules 2003, 36, 10071−10078. (7) Prasitnok, K.; Wilson, M. R. A coarse-grained model for polyethylene glycol in bulk water and at a water/air interface. Phys. Chem. Chem. Phys. 2013, 15, 17093−17100. (8) Vao-Soongnern, V. A multiscale simulation model for polyethylene oxide. Polym. Sci., Ser. A 2014, 56 (6), 928−935. (9) Knapp, E. W.; Irgens-Defregger, A. Off-lattice Monte-Carlo method with constraints - long time dynamics of a protein model without nonbonded interactions. J. Comput. Chem. 1993, 14, 19−29. (10) Hoffmann, D.; Knapp, E. W. Polypeptide folding with off-lattice Monte Carlo dynamics: the method. Eur. Biophys. J. 1996, 24, 387− 403. (11) Hoffmann, D.; Knapp, E. W. Folding pathways of a helix-turnhelix model protein. J. Phys. Chem. B 1997, 101, 6734−6740. (12) Schofield, J.; Ratner, M. A. Monte Carlo methods for short polypeptides. J. Chem. Phys. 1998, 109, 9177−9191. (13) Ulmschneider, J. P.; Jorgensen, W. L. Monte Carlo backbone sampling for polypeptides with variable bond angles and dihedral angles using concerted rotations and a Gaussian bias. J. Chem. Phys. 2003, 118, 4261−4271. (14) Ulmschneider, J. P.; Ulmschneider, M. B.; Di Nola, A. Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations. J. Phys. Chem. B 2006, 110, 16733−16742. (15) Lee, H.; de Vries, A. H.; Marrink, S.-J.; Pastor, R. W. A CoarseGrained Model for Polyethylene Oxide and Polyethylene Glycol: Conformation and Hydrodynamics. J. Phys. Chem. B 2009, 113, 13186−13194. (16) Schluter, A. D. The tenth anniversary of Suzuki polycondensation (SPC). J. Polym. Sci., Part A: Polym. Chem. 2001, 39 (10), 1533− 1556. (17) Hu, X.-D.; Jenkins, S. E.; Min, B. G.; Polk, M. B.; Kumar, S. Rigid-Rod Polymers: Synthesis, Processing, Simulation, Structure, and Properties. Macromol. Mater. Eng. 2003, 288 (11), 823−843. (18) Baschnagel, J.; Binder, K.; Doruker, P.; Gusev, A. A.; Hahn, O.; Kremer, K.; Mattice, W. L.; Mueller-Plathe, F.; Murat, M.; Paul, W.; Santos, S.; Suter, U. W.; Tries, V. Bridging the Gap Between Atomistic and Coarse-Grained Models of Polymers: Status and Perspectives. Viscoelasticity, Atomistic Models, Statistical Chemistry 2000, 152, 41− 156. (19) Mueller-Plathe, F. Coarse-Graining in Polymer Simulation: From the Atomistic to the Mesoscopic Scale and Back. ChemPhysChem 2002, 3 (9), 754−769. (20) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K. Hierarchical Modeling of Polystyrene: From Atomistic to Coarse-Grained Simulations. Macromolecules 2006, 39 (19), 6708− 6719. (21) Harmandaris, V. A.; Reith, D.; van der Vegt, N. F. A.; Kremer, K. Comparison Between Coarse-Graining Models for Polymer Systems: Two Mapping Schemes for Polystyrene. Macromol. Chem. Phys. 2007, 208 (19−20), 2109−2120. (22) Wang, Q.; Keffer, D. J.; Nicholson, D. M. A coarse-grained model for polyethylene glycol polymer. J. Chem. Phys. 2011, 135 (21), 214903. 2399


Article

Journal of Chemical Theory and Computation S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30 (10), 1545−1614. (44) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781−1802. (45) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613. (46) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177. (47) Lal, M. Monte Carlo Computer Simulation of Chain Molecules.I. Mol. Phys. 1969, 17 (1), 57. (48) Wall, F. T.; Erpenbeck, J. J. New Method for the Statistical Computation of Polymer Dimensions. J. Chem. Phys. 1959, 30 (3), 634−637. (49) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: New York, 1953; p 672. (50) Devanand, K.; Selser, J. C. Asymptotic-Behavior and LongRange Interactions in Aqueous-Solutions of Polyethylene Oxide. Macromolecules 1991, 24 (22), 5943−5947. (51) Netz, R. R.; Andelman, D. Neutral and charged polymers at interfaces. Phys. Rep. 2003, 380, 1−95. (52) Hermans, J. J. Statistical Thermodynamics of Swollen Polymer Networks. J. Polym. Sci. 1962, 59 (167), 191. (53) Dusek, K.; Patterson, D. Transition in Swollen Polymer Networks Induced by Intramolecular Condensation. J. Polym. Sci. A2 1968, 6 (7), 1209. (54) Yamakawa, H. Modern Theory of Polymer Solutions; Harper & Row Publishers, Inc.: New York, 1971. (55) Flory, P. J. Statistical Mechanics of Chain Molecules; Wiley Interscience: New York, 1969. (56) Caracciolo, S.; Causo, M. S.; Pelissetto, A. End-to-ebd Distribution Function for Dilute Polymers. J. Chem. Phys. 2000, 112, 7693−7710. (57) Kleinjung, J.; Scott, W. R. P.; Allison, J. R.; van Gunsteren, W. F.; Fraternali, F. Implicit Solvation Parameters Derived from Explicit Water Forces in Large-Scale Molecular Dynamics Simulations. J. Chem. Theory Comput. 2012, 8, 2391−2403. (58) Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B. Afiinities of amino acids side chains for solvent water. Biochemistry 1981, 20, 849−855. (59) Wesson, L.; Eisenberg, D. Atomic solvation parameters applied to molecular dynamics of proteins in solution. Protein Sci. 1992, 1, 227−235.

2400


Chemically Realistic Tetrahedral Lattice Models for ... - ACS Publications

Recommend Documents