WangâLandau Reaction Ensemble Method - ACS Publications

Article pubs.acs.org/JCTC

Wang−Landau Reaction Ensemble Method: Simulation of Weak Polyelectrolytes and General Acid−Base Reactions Jonas Landsgesell, Christian Holm, and Jens Smiatek* Institute for Computational Physics, University of Stuttgart, D-70569 Stuttgart, Germany S Supporting Information *

ABSTRACT: We present a novel method for the study of weak polyelectrolytes and general acid−base reactions in molecular dynamics and Monte Carlo simulations. The approach combines the advantages of the reaction ensemble and the Wang−Landau sampling method. Deprotonation and protonation reactions are simulated explicitly with the help of the reaction ensemble method, while the accurate sampling of the corresponding phase space is achieved by the Wang−Landau approach. The combination of both techniques provides a sufficient statistical accuracy such that meaningful estimates for the density of states and the partition sum can be obtained. With regard to these estimates, several thermodynamic observables like the heat capacity or reaction free energies can be calculated. We demonstrate that the computation times for the calculation of titration curves with a high statistical accuracy can be significantly decreased when compared to the original reaction ensemble method. The applicability of our approach is validated by the study of weak polyelectrolytes and their thermodynamic properties.

1. INTRODUCTION Polyelectrolytes with acidic side groups can be distinguished by their affinity to dissociate protons. Strong polyelectrolytes like DNA are fully deprotonated or ionized for nearly all pH values. In contrast, weak polyelectrolytes like poly(acrylic acid) or most proteins1 bear titrable groups whose protonation state explicitly depends on the pH value of the solution.2 This strong dependency on the pH value gives rise to some interesting phenomena for weak polyelectrolytes like charge regulation3,4 or conformation-protonation state effects.5,6 Hence, the detailed simulation of chemical reactions is of fundamental importance to achieve deeper insights into configurational properties of weak polyelectrolytes as well as the influence of solution effects.7−14 A protonation/deprotonation reaction of a titrable group can be regarded as a Brønsted acid−base reaction

n̅ =

where HA denotes the titrable group, A− denotes the deprotonated group, and H+ denotes the dissociated proton. It has to be noted that the presence of water is omitted for the sake of simplicity.15 A simple definition for the pH value of the solution reads pH = −log10[H+]. The corresponding apparent reaction constant for strong acids is given by [A−][H+] [HA]

(1)

where the brackets [·] denote the concentration of the species,15 whereas the apparent logarithmic reaction constant can be written as pKa = −log10 Ka with values between pKa = 2−10 for standard weak polyelectrolytes.16 The degree of association © XXXX American Chemical Society

(2)

represents the number of associated titrable groups NHA when compared to the total number of titrable groups N0 = NHA + NA− with the number of deprotonated units NA−.17,18 For higher pH values, the titrable groups are more dissociated (n̅ → 0), whereas for lower pH values, the titrable groups remain in the associated state (n̅ → 1). Standard molecular dynamics (MD) simulations are not well suited to study acid−base reactions. This can be mostly rationalized by the impracticability of atomistic force fields to model the cleavage of covalent bonds. More sophisticated approaches like ReaxFF19 or empirical valence bond (EVB) molecular dynamics20 were specifically designed to consider chemical reactions in classical simulations but are often computationally more expensive when compared to standard force fields. Moreover, these approaches are intended to be used for small systems as they often rely on ab initio reference data. Thus, for larger solutes, for example weak polyelectrolytes, reduced and computationally cheap models with an effective and thermodynamically consistent description of the occurring chemical reactions have to be used. These models are commonly called coarse-grained models due to their inherent reduction in complexity and the underlying ability of using larger time steps. Hence, a full resolution of the chemical reaction is not required, while the concentration of the protonated and deprotonated configurations has to obey eq 1 and therefore the principles of thermodynamics. Nevertheless, also for coarse-grained models,8−10,21−30 the influence of

HA ⇌ A− + H+

Ka =

NHA N0

Received: August 10, 2016 Published: December 28, 2016 A

DOI: 10.1021/acs.jctc.6b00791 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

2. THE WANG−LANDAU REACTION ENSEMBLE METHOD In the following section, the basic principles of our new WLRE method are described. The WLRE approach is a combination of the Wang−Landau and the reaction ensemble method. Therefore, we first give a brief introduction to the two original methods before we discuss the properties of the new combined algorithm. For a more detailed description of the reaction ensemble and the Wang−Landau approach, we refer the reader to the original publications.31,32,38,39 2.1. The Reaction Ensemble Method. The reaction ensemble (RE) method provides the possibility to simulate acid−base reactions in chemical equilibrium by proposing changes in the particle numbers due to forward (deprotonation) and backward (protonation) reactions.31 As outlined by Turner et al.,32 the underlying reaction ensemble with fluctuating particle numbers can be derived from the grand canonical ensemble via the separation of the kinetic and the configurational canonical partition sum. A chemical reaction can be written as

protonation/deprotonation reactions and varying pH values can be only studied with the help of more elaborate algorithms. Most of these approaches consider a constant pH value which results in fluctuations of the corresponding protonation states. In contrast, other approaches rely on the explicit value of the chemical reaction constant Ka according to eq 1. As a representative example, the reaction ensemble method31,32 provides the possibility to consider protonated and deprotonated configurations for weak polyelectrolytes with the correct thermodynamic behavior.33−36 A common feature of all these methods is the dominating occurrence of the most probable dissociation state. Less probable protonation states, e.g. a fully protonated or a fully deprotonated polyelectrolyte, respectively, are rare events which do not necessarily have to be realized in simulation time due to their low occurrence probability. Despite this fact, several rare event sampling techniques which explicity focus on less probable simulation states were already introduced over the last decades.37 Specifically the Wang−Landau (WL) method38,39 has been shown to be widely applicable in the context of polymer studies using Monte Carlo simulations.40−43 The WL approach can be regarded as a flat histogram method which provides an accurate sampling of the accessible phase space. The method introduces a bias potential energy which enables the equal occurrence of less probable states such that the density of states or the partition sum becomes accessible. In this article, we present a new algorithm which can be regarded as a combination of the reaction ensemble method and the Wang−Landau approach. The algorithm is applied to study the influence of weak acidic groups in macromolecules or at surfaces, but it can be also used for the thermodynamic study of arbitrary acid−base reactions without loss of accuracy. As an extension, our method provides an accurate sampling of the density of states in contrast to the original reaction ensemble method. In addition, a simple approach for the calculation of titration curves is described which significantly decreases computation times when compared to the original reaction ensemble method for identical model systems. The advantages of the new method are demonstrated by the calculation of thermodynamic properties for weakly charged polyelectrolytes. Thus, our method can be regarded as a significant extension of the reaction ensemble method in order to be used in combination with coarse-grained models. The stochastic nature of these models prohibits the explicit description of chemical reactions such that the resulting protonated or deprotonated configurations can be considered as Markovian states. In fact, our extension provides access to the partition sum as the most important quantity for the study of thermodynamic and statistical mechanical properties of the systems. We will demonstrate the benefits of the proposed method by the discussion of the density of potential energy states, the heat capacity, and the free energy of specific model systems. Moreover, our method provides an efficient and reliable sampling of the phase space and significantly accelerates the calculation of titration curves with high resolution. The article is organized as follows. In the next section we discuss the general properties of the Wang−Landau Reaction Ensemble (WLRE) method. In section three, we present the simulation details of three test systems: the ideal reacting gas, a fixed rodlike and a flexible weak polyelectrolyte with titrable groups. The simulation results are discussed in the fourth section. We briefly conclude and summarize in the last section.

z

∑ νisi = 0 i=1

(3)

for z chemical species of type si with stoichiometric coefficients νi.15 Particle insertion or deletion can be brought into agreement with the law of mass conservation Ni = Ni0 + νξ i

(4)

where Ni is the number of particles after a reaction, N0i is the number of particles prior to a reaction, and ξ is the “extent” of the reaction. The value for ξ is selected randomly with ξ ± 1, where a deprotonation (forward) reaction is given by ξ = +1 and a protonation (backward) reaction by ξ = −1. The transition probability in the reaction ensemble for a forward reaction from state k to l in terms of an individual reaction step ξ = +1 with regard to detailed balance conditions32 reads z ⎡ ⎫ ⎧ ⎤ N 0! ξ ⎥ exp(−β ΔEpot, k → l)⎬ = min⎨1, (βP 0V )νξ̅ K ξ ∏ ⎢ 0 i pkRE, →l (Ni + ξνi)! ⎦ ⎭ ⎩ i=1 ⎣ ⎪

⎪

⎪

⎪

(5)

with the dimensionless reaction constant K, which is proportional to the apparent reaction constant in eq 1, the standard pressure P0, the potential energy difference with ΔEpot,k→l = Epot, k − Epot,l, the volume of the system V, the total change in the number of molecules ν̅ = ∑iνi, and the inverse temperature β = 1/kBT with the Boltzmann constant kB and the temperature T. The dimensionless reaction constant is an input parameter which can be calculated via ⎛ ∑z νμ 0 ⎞ i K = exp⎜⎜ i = 1 i ⎟⎟ k T ⎝ ⎠ B

(6)

with the chemical standard potential μ0 for each species. The corresponding protonation and deprotonation reactions are usually performed after a fixed number of MD simulation steps with constant particle numbers. 2.2. The Wang−Landau Method. The Wang−Landau (WL) approach38,39 is an importance sampling algorithm which was developed for Monte Carlo simulations. The algorithm forces system configurations rN⃗ to occur in accordance with a probability weight B


Article

Journal of Chemical Theory and Computation p( r ⃗ N ) =

1 Ω(E( r ⃗ N ))

with the partition sum47 Z(n̅) = Zkin(n̅)∑EpotΩ(n̅, Epot) exp(−βEpot). The kinetic partition sum Zkin(n), ̅ which is related ⎤ N0 ! z ⎡ νξ̅ ξ 0 to the reaction ensemble factor (βP V ) K ∏i = 1 ⎣⎢ (N 0 +i ξν ) ! ⎦⎥, i i can be calculated as outlined in the Supporting Information. Optionally, it is also possible to study a two-dimensional set of reaction coordinates Γ(n,̅ Epot). The acceptance probability between two reaction states k = (n̅k, Epot,k) and l = (nl̅ , Epot,l) reads

(7)

which is inversely proportional to the corresponding density of states Ω(E), more accurately coined the degeneracy of states with energy E. Since the density of states, and therefore also the probabilities, are a priori unknown, they have to be constructed iteratively during the simulation39 by counting the number of occurrences. The algorithm typically evaluates an actual WL potential

Sk(Ek ) = ln Ω(Ek )

z ⎡ ⎧ ⎫ ⎤ Ni0! 0 ξ νξ̅ ξ ⎬ ⎨ ⎢ ⎥ pkWLRE, min 1, ( P V ) K exp( S S ) = β − ∏ k l 0 →l (Ni + ξνi)! ⎦ ⎭ ⎩ i=1 ⎣ (12)

(8)

for each configuration k which allows the transition from one state k to another state l with energies Ek and El via pkWL = min{1, exp(Sk − Sl)} →l

⟨(⟩ =

∑ N ... ∑ N ∑E (({Ni}, E)Ω({Ni}, E) exp(−β(E + ∑i Niμi ) z

1

∑ N ... ∑ N ∑E Ω({Ni}, E) exp(−β(E + ∑i Niμi )) 1

z

(14)

where ∑N1...∑Nz denotes the sum over the particles for all species {Ni} with the chemical potentials μi and the joint degeneracy of the total energy Ω({Ni}, E). Thus, we considered the individual spatial configurations of the system by including them in the corresponding associated total energies. The connection between the chemical potentials μi and the number of particles Ni for an acid−base reaction in the reaction ensemble yields z

∑ Niμi = μHA NHA + μH+(N0 − NHA) + μ A−(N0 − NHA) = μHA N0 i=1

(15)

z

⎫ exp( −β ΔEpot, k → l) exp(Sk − Sl)⎬ ⎭

with regard to the definition of the chemical equilibrium ∑zi νiμi = 0 such that μA− + μH+ = μHA. Due to eq 15, the exponential factor exp(−β∑iNiμi) in eq 14 disappears and ensemble averages can be calculated via

⎪

⎪

(13)

with the potential energy and the degree of association as reaction coordinates. It has to be noted that discretization errors in the Wang−Landau potential, i.e. the discrete range covered by a single bin, can be minimized as described in ref 48. In order to avoid errors due to wrongly chosen values for the potential energy boundaries,49 the strategy of preliminary simulation runs as described in ref 47 is advised. 2.4. Efficient Calculation of Ensemble Averages and Titration Curves. It was previously discussed that the reaction ensemble is closely connected to the grand canonical ensemble.32 Therefore, ensemble averages for a quantity ( can be calculated via

⎪

(10)

⟨(⟩ =

which is the actual acceptance probability for the WLRE approach. In fact, the consideration of the single reaction coordinate n̅ is achieved by the occurrence of the factor exp(−βΔEpot, k →1) in eq 10 which therefore slightly differs from the original WL approach represented by eq 9. The resulting Wang−Landau potential is given by S(n ̅ ) = ln(Z(n ̅ ))

⎪

S(n ̅ , Epot) = ln(Z kin(n ̅ )Ω(n ̅ , Epot))

⎧ ⎤ ⎡ N 0! ⎥× = min⎨1, (βP 0V )νξ̅ K ξ ∏ ⎢ 0 i ⎩ i = 1 ⎣ (Ni + ξνi)! ⎦ ⎪

⎪

⎪

in accordance to the original WL approach in terms of eq 9. Hence, the WL potential can be written as

(9)

which represents the actual acceptance transition probability between two states.38 For the accepted state, which can be either j = k or j = l, the Wang−Landau potential S(Ej) is increased by the Wang−Landau parameter f WL in accordance to S(Ej) → S(Ej) + f WL. Furthermore, the histogram which tracks the number of occurrences is updated according to H(Ej) → H(Ej) + 1. If the shape of the histogram satisfies a global flatness criterion, the Wang−Landau parameter is reduced,44 and the histogram entries are reset to zero. This iterative procedure continues until the Wang−Landau parameter f WL converges to zero within an upper threshold44 which is usually between f cWL = 10−4−10−8. An alternative updating rule was recently proposed in ref 45. In this publication, the usage of Gaussian kernels was advised in order to minimize finite size effects for the histogram bins. The evaluation of continuous functions for discrete values of the reaction coordinate was also addressed in other free energy methods.46 2.3. Description of the Wang−Landau Reaction Ensemble Algorithm. In order to combine the Wang− Landau method with the reaction ensemble algorithm, we perform reaction steps in agreement with the original RE method.31,32 We randomly select between a forward and a backward reaction and perform trial moves in accordance to the suggested reaction. In fact, the most obvious reaction coordinate for a protonation/deprotonation reaction is represented by the degree of association n̅ defined in eq 2. Herewith, the corresponding one-dimensional reaction coordinate space can be transformed to Γ(n)̅ in contrast to Γ(E) for the original WL approach. The acceptance probability for the transition into the new state is then given by the combination of eq 5 and eq 9 and reads ξ pkWLRE, →l

⎪

∑n ∑E ((n ̅ , E)Ω(n ̅ , E) exp( −βE) ∑n ∑E Ω(n ̅ , E) exp( −βE)

(16)

as a sum over all occurring total energies E and over all degrees of association n̅, where Ω(n̅, E) can be interpreted as the degeneracy of the total energy at a given degree of association. For a static observable ((n ̅ , Epot), which only depends on n̅ and the potential energy Epot, eq 16 can be rewritten as

(11) C


Article

Journal of Chemical Theory and Computation ⟨(⟩ =

∑n Z kin(n ̅ ) ∑ E ((n ̅ , Epot)Ω(n ̅ , Epot) exp(− βEpot) pot

⟨((n ̅ , Epot)⟩ for arbitrarily chosen input reaction constants. In fact, this property is beneficial for the calculation of titration curves (section 4) where a large number of sampling points can be obtained without additional computational cost. More 0 z ⎡ N ! ⎤ details on the factor (βP 0V )νξ̅ K ξ ∏i = 1 ⎢⎣ (N 0 +i ν ) ! ⎥⎦ with the i i input reaction constant K, which is necessary for the calculation of the kinetic partition sum Zkin(n), ̅ are outlined in the Supporting Information. In addition, Wang−Landau potentials for varying input reaction constants can be also obtained by using the analytic expression for the kinetic partition sum Zkin(n,̅ K). In order to follow this approach, the logarithm of ln(ZKkin1 (n̅)) with the reaction constant K1 is subtracted from the Wang−Landau potential S, and the logarithm of the kinetic K partition sum ln(Zkin2 (n)) ̅ with the new value K2 has then to be added. 2.5. Implementation of the WLRE Method. The general flowchart of the Wang−Landau reaction ensemble algorithm can be described as follows. First, the histogram H(Γj) has to be initialized with 0 for all entries Γj in the reaction coordinate space. Furthermore, also the corresponding Wang−Landau potential S(Γj) has to be set to 0 for all values of the reaction coordinates, and the Wang−Landau parameter is initialized with f WL = 1. 1. The chemical composition of the system is changed with a reaction trial move in accordance to the forward or backward reaction for a specific titrable group which is selected uniformly at random. This step is identical to the original reaction ensemble method, and more details can be found in refs 31 and 32. The trial move is accepted or rejected with regard to the Wang−Landau reaction ensemble transition probability in accordance to eq 10 or eq 12. The histogram entry H(Γj) of the final state j is then increased by one, and the Wang−Landau potential S(Γj) is updated with S(Γj) → S(Γj) + f WL. 2. A series of simulation steps to change the particle positions are performed before step 1 is repeated (see below for additional information and ref 47 for more details). 3. If the histogram converges with deviations within a flatness criterion, and the system shows a random walk behavior,38,44 all entries of the histogram are reset to zero, and the Wang− Landau parameter is reduced in accordance to ref 44 with

∑n Z(n ̅ )

(17)

where we used the joint degeneracy of the potential energy Ω(n,̅ Epot) and the identity ∑n̅∑EΩ(n,̅ E) exp(−βE) = ∑nZ(n ). ̅ ̅ Thus, the canonical partition sum can be separated into kinetic and configurational contributions in accordance to Z(n̅) = Zkin(n̅)Zpot(n̅). Moreover, the static observable ((n ̅ , Epot) can be interpreted as the microcanonical ensemble average50,51 which can be directly obtained when using the sampling scheme as defined in eq 12 to yield the sampling average M ((n ̅ , Epot) ≈ 1/M ∑i = 1 ( i(n ̅ , Epot) for M samples. Alternatively, when using the sampling scheme as denoted by eq 10, the canonical ensemble average can be directly written as M ⟨((n ̅ , Epot)⟩n̅ , V , T ≈ 1/M ∑i = 1 ( i(n ̅ , Epot) which yields ⟨(⟩ =

∑n ⟨((n ̅ , E)⟩n̅ , V , T Z(n ̅ ) ∑ n Z( n ̅ )

(18)

as an exact definition for the observable ⟨A⟩ in the WLRE approach. It has to be noted that the above-discussed ensemble average directly depends on the reaction constant K which is 32 included in the kinetic canonical partition sum Zkin(n). ̅ Therefore, an analytic expression for Zkin(n̅) would allow us to calculate ensemble averages ⟨A(n,̅ Epot)⟩ for any arbitrarily chosen value of K. In the Supporting Information, we thus derive an expression for the kinetic partition sum for each discrete degree of association n̅ which reads p

ln(Z kin(n ̅ = pΔn ̅ )) =

∑ j=1

⎛ K (j Δ n ) ⎞ ̅ ⎟ ln⎜⎜ a ⎟ K ⎝ a,input ⎠

(19)

for integer values j = 1,.., N0. The apparent input dissociation constant in the denominator reads K a,input = K (βP 0)

(20)

and the apparent dissociation constant according to the law of mass action (eq 1) is given by K a(n ̅ ) = c0(1 − n ̅ )2 / n ̅

(21)

which has a predefined value in chemical equilibrium. The value for c0 denotes the concentration of titrable groups in the solution according to c0 = N0/V with the volume V. A closed expression for the kinetic partition sum of a simple acid−base reaction can be derived (as shown in the Supporting Information) in order to read

fWL, i + 1

⎧1 f ⎪ if WL, i ≤ 1 ⎪t 2 t =⎨ ⎪ fWL, i else ⎪ ⎩ 2

(23)

where t is the number of trial moves divided by the number of bins in the reaction coordinate space. If the upper option is reached, the simulation stays in this branch. As it was mentioned in step 2, configuration changing moves have to be performed in order to provide an accurate sampling of the configuration space. In more detail, for a twodimensional reaction coordinate space Γ(n,̅ Epot) associated with eq 12, the Monte Carlo moves for new spatial configurations rN1,2 are accepted in accordance to

⎛ ⎛ c ⎞⎞ ⎛ K a,inputx ⎞ ⎟ + ln⎜⎜ 0 ⎟⎟⎟⎟ ln(Z kin(x)) = − N0⎜⎜x + ln x + (x − 1) ln⎜ 2 ⎝ c0(x − 1) ⎠ ⎝ K a,input ⎠⎠ ⎝

(22)

for x = pΔn,̅ which is valid in the thermodynamic limit. Thus, instead of performing multiple independent Wang− Landau simulations for different input reaction constants as in the original reaction ensemble method, we exploit that kinetic partition sums for different reaction constants can be calculated in accordance to eq 19 or eq 22. After the configurational partition sum, or alternatively, the density of states is sufficiently sampled by the Wang−Landau procedure in one simulation run, the kinetic partition sum can be calculated and combined with eq 17 or eq 18 to estimate ensemble averages

⎧ Ω(n1 , Epot,1) ⎫ ̅ ⎬ pr N → r N = min⎨1, 1 2 ⎩ Ω(n2̅ , Epot,2) ⎭ ⎪

⎪

⎪

⎪

(24)

which results in a modification of the WL potential and also in an update of the corresponding histogram entries after each D


Article

Journal of Chemical Theory and Computation ⎛ ⎛ r − r0 ⎞2 ⎞ 1 2 UFENE(r ) = − krmax log⎜⎜1 − ⎜ ⎟⎟ 2 ⎝ rmax ⎠ ⎟⎠ ⎝

trial move.47 In contrast, the Wang−Landau potential remains 47 unchanged for a one-dimensional reaction coordinate Γ(n). ̅ Hence, Γ(n)̅ does not depend on the potential energy Epot, and therefore the configurational phase space can be sampled by the acceptance probability ⎧ exp( −βEpot,2) ⎫ ⎬ pr N → r N = min⎨1, 1 2 ⎩ exp( −βEpot,1) ⎭ ⎪

⎪

⎪

⎪

(26)

with the spring constant k = 10 ϵ/σ , the maximal elongation rmax = 1. 5σ, and an equilibrium length r0 = 21/6σ ≈ 1.12σ. For the rodlike weak polyelectrolyte, the number of titrable groups was N0 = 50 corresponding to a monomer concentration of c0 = 0.00028σ−3. An apparent reaction constant Ka = KβP0 with K = 0.0088, a standard pressure of P0 = 1 bar = 0. 00108 ϵ/σ3, and β = 1/kBT = 1ϵ−1 were used as input parameters. The proton concentration of the solution can be calculated by [H+] = c0(1 − ⟨ n⟩) ̅ and pH = − log10(c0(1 − ⟨n ̅ ⟩)) (27) 2

(25)

which is the original Metropolis criterion. Moreover, new configurations can be also realized by ordinary molecular dynamics moves or hybrid molecular dynamics/Monte Carlo schemes as an alternative. Finally, it is worth mentioning that the proposed algorithm can be trivially parallelized by dividing the reaction coordinate space into multiple windows which can be sampled simultaneously.52 In addition, also any other arbitrarily chosen reaction coordinate instead of n̅ can be used53 to circumvent the problem of hidden complexities in low-dimensional reaction coordinate spaces.54

which has to be combined with pKa = −log10 Ka in order to obtain the actual pKa − pH value. We used a Langevin thermostat (according to step 2 in section 2.5) when using the one-dimensional reaction coordinate space mi ri ⃗ ̈ = −γ ri ⃗ ̇ + ηi⃗ + Fi ⃗

3. SIMULATION DETAILS As a test case for our method, we studied the properties of a weak rodlike polyelectrolyte by the WLRE method in terms of a coarse-grained model with one-dimensional periodic boundary conditions. A rod-shaped weak polyelectrolyte was simulated by fixing discrete monomers into a rod shape conformation, and the cubic simulation box was constrained to a cylinder of height b = 56.3124σ and radius R = b/√π in agreement with the cell model approach55 (see Figure 1). In

(28)

with the conservative forces F⃗i, the friction forces −γ ri ⃗ ̇, and the random forces η⃗i. All masses were set to mi = 1 due to the fact that the evaluation of the kinetic partition sum effectively relies on the consideration of the input reaction constant, as it is shown in the Supporting Information. The random force acts on each particle independently and obeys the fluctuation− dissipation theorem with ⟨ηiα⟩ = 0 and ⟨ηiα(t)ηjβ(t′)⟩ = 2γkBTδijδαβδ(t − t′) which ensures the presence of Gaussian white noise for particles i and j in the spatial directions α and β. The friction coefficient had a value of γ = 1.0σ−1(mϵ)1/2. The temperature was set to T = 1ϵ/kB, and the Langevin equation was integrated by a Velocity Verlet algorithm with a time step of δt = 0.01 σ(m/ϵ)1/2. Electrostatics for the rodlike polyelectrolyte with periodic boundary conditions in one direction were calculated by the MMM1D method58 and for the flexible polyelectrolyte by the P3 M method.59 The Bjerrum length was set to λB = e2/(4πε0εrkBT) = 2σ if not otherwise mentioned, with the dielectric constant ϵr and the elementary charge e. The Wang−Landau reaction ensemble method has been implemented in the freely available software package ESPResSo.60,61

4. NUMERICAL RESULTS 4.1. Ideal Reacting Gas. We start this section by the discussion of the titration curve for the ideal reacting gas. The ideal reacting gas offers the opportunity to test the new algorithm due to the existence of an analytical solution for the titration curve 1 n̅ = 1 − (29) 1 + 10 pK a − pH

Figure 1. Schematic drawing of the rodlike weak polyelectrolyte with titrable groups for the cell model simulations.

and further allows us to compare our results with the original reaction ensemble method. All conservative interactions are absent, which significantly simplifies the calculations, and the particles can only perform protonation or deprotonation steps. The analytical ideal titration curve in combination with the RE and WLRE simulation results for different numbers N0 = 80 and 200 of titrable units is shown in Figure 2. One can observe that the simulated titration data shows a perfect agreement with the analytical curve presented in eq 29. Furthermore, a good

contrast, the ideal reacting gas and the flexible polyelectrolyte were simulated in a cubic simulation box with full periodic boundary conditions and identical parameters and simulation protocols as used for the rodlike weak polyelectrolyte. All titrable groups repel each other by a WCA potential56 with amplitude ϵ = 1kBT, diameter σ, and a cutoff radius rc = 21/6σ. Bonds between adjacent titrable groups (beads) were modeled by a FENE potential57 according to E


Article


show a good agreement with the results obtained by the WLRE method. In the presence of electrostatic interactions, and for pKa − pH ≤ 1, the associated state is more favored compared to the ideal analytical results (ns̅ im ≫ n̅id) due to strong repulsive electrostatic interactions between like-charged species. Hence, the resulting curve in the presence of electrostatic interactions is less steep when compared to the ideal titration curve. In order to benefit from the combination with the Wang− Landau sampling approach, we now focus on the calculation of thermodynamic properties. Thus, we considered a twodimensional reaction coordinate space Γ(n̅, Epot) in accordance to eq 12 with a Wang−Landau potential which depends on both reaction coordinates (eq 13). From the corresponding two-dimensional Wang−Landau potential S(n̅, Epot), the logarithm of the density of potential energy states ln Ω(n̅, Epot) is obtained by subtracting the logarithm of the kinetic partition sum in accordance to eq 19. The density of potential energy states for a rodlike weak polyelectrolyte with N0 = 50 titrable groups and λB = 2σ is shown in Figure 4. When

Figure 2. Titration curves for values pH − pKa with degrees of association n̅ for N0 = 80 (red squares) and N0 = 200 (blue circles) as calculated by the WLRE method and for N0 = 200 by the RE method (black filled squares). The black solid line represents the analytical solution for the ideal titration curve according to eq 29.

agreement between the results obtained by the WLRE and the standard reaction ensemble method can be recognized. In addition, identical outcomes for different numbers of titrable groups indicate that finite-size effects are sufficiently small. As a rough time estimate, the reaction ensemble took around two CPU hours for the computation of one data point in the ideal system, whereas the Wang−Landau reaction ensemble needs one CPU day until histogram convergence is reached. Hence, for this simple system, it becomes obvious that the WLRE method is superior to the RE method if the analytic expressions of eq 19 or eq 22 are used and 12 or more data points are needed. 4.2. Rodlike Weak Polyelectrolyte. As another test system, we simulated an infinite linear rodlike weak polyelectrolyte with fixed equidistant positions for the titrable groups in a confining cylinder. This model is the equivalent of the often studied quenched charged rod cell model,55,62,63 that has been solved using the Poisson−Boltzmann approach62,63 and by simulations.55 Electrostatic interactions between charged species are omnipresent, and the Bjerrum length is λB = 2σ. The degrees of association for different pKa − pH values (titration curves) are shown in Figure 3. The results of the original reaction ensemble method for different pKa values

Figure 4. Logarithm of the density of potential energy states ln Ω(n̅, Epot) for the rodlike weak polyelectrolyte with the reaction coordinates Epot and degree of association n̅. The number of titrable groups was N0 = 50, the Bjerrum length was λB = 2σ, and the values for ln Ω(n,̅ Epot) include an arbitrary shift value. Regions with a high and a low value for ln Ω(n,̅ Epot) are denoted in the top of the figure.

calculating the sum over the degeneracy at a given degree of association n̅, it becomes evident that for low values of n̅, a larger number of states exist than for a degree of association n̅ ≈ 1. In addition, the extreme value n̅ = 1 is only realized by a single state with a well-defined value for the potential energy. The corresponding free energy landscape of the rodlike weak polyelectrolyte for the one-dimensional reaction coordinate space n̅ is calculated by F(n ̅ ) = −kBT ln(Z(n ̅ )) (30) and compared to the results for an ideal rod. The outcomes are shown in Figure 5. It can be clearly seen that a global minimum for the free energy of a charged polyelectrolyte at n̅ ≈ 0.88 exists. The minimum in the free energy F({Ni}, V, T) for a fixed volume V and a temperature T is directly induced by the chemical equilibrium condition dF = ∑iνiNi = 0. Thus, for an ideal reacting gas without interactions, the degree of association n̅min which is associated with the minimum free energy directly reflects the chosen value for the apparent reaction constant via

Figure 3. Degree of association n̅ for a weak rodlike polyelectrolyte for different pKa − pH values (red squares) with N0 = 50 and λB = 2σ. The black filled circles denote the results for the original reaction ensemble method, whereas the red squares represent the values obtained by the Wang−Landau reaction ensemble method. The black solid line represents the solution for the ideal gas in accordance to eq 29. F


Article


Figure 6. Heat capacity divided by the number of free protons as a function of pKa − pH calculated by the Wang−Landau reaction ensemble method for the fixed rodlike polyelectrolyte (N0 = 50 and λB = 2σ) with electrostatic interactions (red solid line) and a fixed ideal polyelectrolyte without electrostatic interactions (blue dashed line). The black square is the solution for the heat capacity obtained by the Poisson−Boltzmann cell model theory (more details in the text and the Supporting Information). The black circle is the heat capacity obtained by a canonical simulation (without reactions) for a fully dissociated charged rod in a cylinder.

Figure 5. Free energy of the rodlike weak polyelectrolyte (red squares) as a function of the degree of association for N0 = 50, K = 0.0088, P0 = 0.00108ϵ/σ3, β = 1/ϵ, and λB = 2σ. The free energy of an ideal reacting gas is denoted by the blue circles.

the ideal titration curve (eq 29). The corresponding deviations between the ideal reacting gas and the rodlike weak polyelectrolyte can be attributed to the presence of electrostatic interactions. As an additional standard thermodynamic quantity, the heat capacity at constant volume can be calculated from statistical mechanics by the knowledge of the partition function or the density of states. A previous publication64 demonstrated the benefits of this quantity for polyelectrolytes in order to detect transitions in the counterion condensation behavior. Moreover, the heat capacity is also used frequently in numerical studies on phase transitions in polymeric systems (see for example in refs 51 and 65−67). Thus, the introduction of an accurate method to study the behavior of this quantity is highly beneficial. The heat capacity is defined by

CV =

∂⟨E⟩ ∂T

V

exists which is in agreement to the properties found for strong polyelectrolytes. The differences between the heat capacity for the rodlike weak polyelectrolyte and the ideal rod for values pKa − pH ≤ − 5 can be fully attributed to the presence of electrostatic interactions. Exact analytical solutions cannot be easily derived due to the electrostatic coupling between like- and differently charged species. However, for all values pKa − pH ≤ −5, one can assume Ka ≫ 1 which means that the rod is fully charged and can be regarded as a strong polyelectrolyte without any protonation reactions. An exact solution for the Poisson− Boltzmann potential can be obtained by employing the Poisson−Boltzmann cell model (PB)55,62,63 with an inner cylinder radius r0 = 1σ and an outer cylinder radius R as depicted in Figure 1. The internal mean field potential energy ⟨Epot⟩PB = ∫ Vd3r(ϵr/2)[∇Ψ(r, T)]2 with the dielectric constant ϵ can be calculated from the electrostatic potential Ψ(r) as derived from the Poisson−Boltzmann cell model55 via integration over the cylindrical volume V which allows us to

(31)

which gives CV =

⎡ 2 ⎤ 3⟨N ⟩ 3 3 + kBβ 2⎢⟨Epot ⟩ − ⟨Epot⟩2 + ⟨NEpot⟩ − ⟨N ⟩⟨Epot⟩⎥ ⎣ ⎦ 2βT 2β 2β

(32)

∂⟨Epot⟩PB

compute the excess heat capacity CV ,pot,PB = ∂T as it was described in ref 64. Finally, the total heat capacity of the rodlike polyelectrolyte can be obtained by adding the kinetic ideal gas contribution 3/2⟨N⟩kB. The resulting value CV/⟨N⟩ ≈ 174.26 kB/50 is shown in Figure 6 as a black square. The differences between the results of the WLRE method and the PB approach can be mainly attributed to correlation effects between the protons that arise for electrostatic coupling parameters Ξ = (q2λB/κ) > 0 with the valency q and the Gouy−Chapman length κ.64 In fact, the PB approach is only strictly valid for Ξ → 064 which explains the deviations from the simulation results due to correlation effects between the protons. In order to further assess the validity of the WLRE simulations, a canonical simulation of a strong rodlike polyelectrolyte without reactions was performed to estimate the heat capacity in the limit of strong dissociation. As it can be seen in Figure 6, the values for the heat capacity obtained by the WLRE and the NVT simulations (strong polyelectrolyte) reveal a perfect coincidence. More details, for example the proton distribution around the rods and exact expressions for the heat capacity in the PB approach, can be found in the Supporting Information.

by using eq 17 for the evaluation of the ensemble averages in the brackets where N denotes the number of freely moving particles (protons). The first term in eq 32 is associated with the kinetic energy contributions. For the studied system, the titrable groups N0 = 50 are fixed in their position, such that only the dissociated protons with an actual number of N = N0(1 − ⟨ n⟩) ̅ contribute to the kinetic energy. Thus, we assume that the kinetic partition sum is not altered although it was derived for free particles.32 The heat capacity for the rodlike weak polyelectrolyte is shown in Figure 6 in comparison to the heat capacity of a rodlike ideal polyelectrolyte without any electrostatic interactions. A maximum value of TCV/⟨N⟩ at pKa − pH ≈ −2.57 can be observed. In contrast, the results for the ideal fixed weak polyelectrolyte reveal a constant value of 3/2 kB for all pKa − pH values. In fact, the maximum peak in the heat capacity for the rodlike polyelectrolyte between values pKa − pH ≈ − 4.5 and pKa − pH ≈ −2 was discussed as an indicator for the onset of counterion condensation behavior around strong polyelectrolytes.64 However, although this effect will be studied in more detail in the future, it might indicate that a second order phase transition64 for weak polyelectrolytes G


Article

Journal of Chemical Theory and Computation 4.3. Flexible Weak Polyelectrolytes. In this subsection, we study the influence of different Bjerrum lengths on the titration curves of a flexible weak polyelectrolyte in order to test the applicability of our method for more realistic systems. The corresponding results for a flexible weak polyelectrolyte with N0 = 50 titrable groups and different Bjerrum lengths are shown in Figure 7. In comparison, we also depict the ideal titration curve

Figure 8. Canonical free energy difference ΔFn̅,V,T between the individual points and the free energy minimum as a function of the degree of association n̅ for different values of the dissociation constant Ka at Bjerrum length λB = 2σ. The gray line denotes the corresponding titration curve as depicted in Figure 7.

coordinates in order to study the free energy landscape around the resulting titration curves (Figure 7). It is clearly visible that the corresponding titration curve for the flexible weak polyelectrolyte follows the shape of a pronounced and steep free energy basin. This finding can be rationalized by the fact that the degree of association mostly fluctuates around the free energy minimum in accordance with the definition of the chemical equilibrium dF = ∑iνiNi = 0. The occurrence of the deep free energy basin further rationalizes the applicability of the standard reaction ensemble method in order to determine titration curves. Hence, regions with a large distance to the titration curve are inaccessible for brute-force methods as demonstrated by pronounced free energy differences of ΔFn̅,V,T ≥ 60 kBT. The results of the average end-to-end distance for flexible weak polyelectrolytes

Figure 7. Degree of association n̅ (titration curves) for a flexible weak polyelectrolyte (N0 = 50 titrable groups) at different Bjerrum lengths λB = 0.5, 2.0, 4.0, 8.0, and 10.0σ for varying pKa − pH values. The black line indicates the ideal titration curve in accordance to eq 29.

(black solid line) for a polyelectrolyte without any conservative interactions. As one would expect, the titrable groups at all Bjerrum lengths are more dissociated (n̅ → 0) for smaller pKa − pH values, which correspond to higher pH values, whereas for lower pH values, the titrable groups remain in the associated state (n̅ → 1). Furthermore, it can be clearly seen that stronger electrostatic interactions in terms of increasing Bjerrum lengths induce nonmonotonic deviations from the ideal titration curve. Hence, the degree of association is strongly affected by the magnitude of the electrostatic interactions. In fact, for different values of λB, the corresponding value for n̅ = 0.5 is shifted from pKa − pH = 0 (ideal curve) to pKa − pH = −0.75 (λB = 0.5σ), pKa − pH = −1.75 (λB = 2σ), pKa − pH = −1.7 (λB = 4σ), pKa − pH = 0 (λB = 8σ), and pKa − pH = 0.75 for λB = 10σ, respectively. These results indicate that stronger electrostatic interactions favor the association of titrable groups at equal pKa − pH values for low and moderate Bjerrum lengths λB ≤ 4σ. In this regime, stronger repulsive interactions between likecharged species for increasing Bjerrum lengths become energetically more unfavorable such that the associated state (n̅ = 1) is preferred. As a possible interpretation, one has to notice the reciprocal dependence of the Bjerrum length on the dielectric constant according to λB ∼ 1/ϵr. Thus, it can be concluded that apolar solvents with lower dielectric constants favor the presence of the associated state in agreement with recent results obtained by all-atom molecular dynamics simulations.68 For even stronger stronger electrostatic interactions (λB ≥ 8σ), the titration curves reveal a “counterintuitive” shift in the opposite direction of the ideal titration curve. It was already discussed that a deprotonation reaction in this regime becomes more favorable due to charge correlation effects, which are accompanied by high electrostatic energy gains.69 As a further powerful application of our method, we depict the free energy landscape for a flexible weak polyelectrolyte at λB = 2σ in Figure 8. The degree of association n̅ and the corresponding pKa − pH values are chosen as reaction

⟨re⟩ = ⟨| r1⃗ − rN⃗ |⟩

(33)

with the position of the first r1⃗ and the last titrable group rN⃗ as a function of the pKa − pH value at different Bjerrum lengths are shown in Figure 9. Our findings indicate that for highly associated weak polyelectrolytes with n̅ → 1 and pKa − pH ≥ 1.5 according to Figure 7, the end-to-end distances are comparable for nearly all Bjerrum lengths. This outcome can be attributed to the occurrence of a low net charge along the

Figure 9. Average end-to-end distances ⟨re⟩ for a flexible weak polyelectrolyte at different Bjerrum lengths. H


Article


explicit water molecules. Although the realization of this approach would require more sophisticated water models, recent advances in this direction are promising.80

polyelectrolyte backbone due to the favored presence of the associated state. For values pKa − pH ≤ 1.5, the presence of two regimes can be observed. For low and moderate Bjerrum lengths λB ≤ 4σ, pronounced differences in the end-to-end distances can be observed for values pKa − pH ≤ − 1. Thus, the value of the end-to-end distance decreases in a nonmonotonic trend according to re(λB = 2σ) > re(λB = 0.5σ) > re(λB = 4σ). In fact, the increase of the end-to-end distance for highly charged polyelectrolytes (n̅ → 0) between λB = 0.5σ and λB = 2σ can be explained by a stronger electrostatic repulsion between the charged titrable groups of the backbone. In contrast, the decrease of the end-to-end distance for highly charged polyelectrolytes and Bjerrum lengths λB ≥ 4σ was also reported in ref 70. It was discussed,69,71 that stronger electrostatic interactions induce counterion condensation/ion correlation effects. At a certain electrostatic interaction strength, even a first order phase transition occurs due to a polymer collapse in terms of correlation effects.72 Thus, the decrease of the end-to-end distance at higher Bjerrum lengths can be attributed to exactly these correlation effects.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00791. Information about the calculation of the kinetic partition sum, the calculation of the Poisson−Boltzmann heat capacity and proton probability densities around the rodlike weak polyelectrolyte (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jens Smiatek: 0000-0002-3821-0690 Notes

5. SUMMARY AND CONCLUSION We presented a new Wang−Landau reaction ensemble (WLRE) approach which can be used for the simulation of weak polyelectrolytes and general Brønstedt acid−base reactions. With this method, we were able to consider the influence of explicit deprotonation and protonation reactions while obtaining at the same time information about the density of states. Our approach reproduces the titration curves of weak polyelectrolytes in agreement with the original reaction ensemble method. We propose a simple expression for the kinetic partition sum which allows us to calculate the full titration curve or further thermodynamic variables by a single simulation run. Therefore, our method significantly decreases the computation time for the calculation of ensemble averages when compared to the original RE approach. Another advantage of our method is the knowledge of the density of states after sufficient sampling of the phase space. Thus, thermodynamic properties like free energies and ensemble averages are accessible, and provide deeper insights into the system behavior. We calculated the heat capacity and the reaction free energy of a weak polyelectrolyte to highlight the benefits of our approach. Moreover, also unfavorable protonation states can be accessed such that an accurate sampling of the phase space is achieved. Hence, the proposed method is fully consistent with detailed balance conditions when histogram convergence is reached. It has to be mentioned that our method is applicable for a realistic system like flexible weak polyelectrolytes. The huge advantage of our algorithm is the access to thermodynamic information. For instance, entropic and enthalpic contributions to the free energy can be separated, which is of specific importance for the understanding of weak polyelectrolyte gels73 and brushes.74−76 Moreover, free energy differences between polymer conformations can be estimated, which is relevant for the mapping between coarse-grained and real proteins.77 With regard to the broad applicability of our algorithm, also modified coarse-grained models can be designed in order to study specific pH dependent unfolding mechanisms.78,79 In principle, the method can be also used for the study of pH dependent protonation states for proteins and other biomacromolecules in atomistic simulations. As a necessary extension, the deprotonation reactions must be complemented by the presence of

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Tobias Richter and Peter Košovan for helpful discussions. Financial funding by the Deutsche Forschungsgemeinschaft through the cluster of excellence initiative ”Simulation Technology” (EXC 310), the SFB 716 and the grant AR 593/7-1 is gratefully acknowledged.

■

REFERENCES

(1) Elaissari, A. Thermally sensitive colloidal particles: From preparation to biomedical application. Prog. Colloid Polym. Sci. 2006, 133, 9−14. (2) Biochemistry, 8th ed.; Berg, J. M., Ed.; Freeman: New York, NY, USA, 2015. (3) Lund, M.; Jönsson, B. On the charge regulation of proteins. Biochemistry 2005, 44, 5722−5727. (4) Lund, M.; Jönsson, B. Charge regulation in biomolecular solution. Q. Rev. Biophys. 2013, 46, 265−281. (5) Castelnovo, M.; Sens, P.; Joanny, J.-F. Charge distribution on annealed polyelectrolytes. Eur. Phys. J. E: Soft Matter Biol. Phys. 2000, 1, 115−125. (6) Shi, C.; Wallace, J. A.; Shen, J. K. Thermodynamic coupling of protonation and conformational equilibria in proteins: theory and simulation. Biophys. J. 2012, 102, 1590−1597. (7) Siegel, R. A.; Firestone, B. A. pH-dependent equilibrium swelling properties of hydrophobic polyelectrolyte copolymer gels. Macromolecules 1988, 21, 3254−3259. (8) Reed, C. E.; Reed, W. F. Monte Carlo study of titration of linear polyelectrolytes. J. Chem. Phys. 1992, 96, 1609−1620. (9) Ullner, M.; Jönsson, B.; Söderberg, B.; Peterson, C. A Monte Carlo study of titrating polyelectrolytes. J. Chem. Phys. 1996, 104, 3048−3057. (10) Ullner, M.; Jönsson, B. A Monte Carlo study of titrating polyelectrolytes in the presence of salt. Macromolecules 1996, 29, 6645−6655. (11) Berghold, G.; van der Schoot, P.; Seidel, C. Equilibrium charge distribution on weak polyelectrolytes. J. Chem. Phys. 1997, 107, 8083− 8088. (12) Zito, T.; Seidel, C. Equilibrium charge distribution on annealed polyelectrolytes. Eur. Phys. J. E: Soft Matter Biol. Phys. 2002, 8, 339− 346. (13) Encyclopedia of smart materials; Schwartz, M. M., Ed.; A Wiley Interscience publication; Wiley: New York, NY, USA, 2002. I


Article


(36) Uhlík, F.; Kosovan, P.; Limpouchová, Z.; Procházka, K.; Borisov, O. V.; Leermakers, F. A. Modeling of ionization and conformations of starlike weak polyelectrolytes. Macromolecules 2014, 47, 4004−4016. (37) Smiatek, J.; Hansen, N.; Kästner, J. In Simulating Enzyme Reactivity: Computational Methods in Enzye Catalysis; Tunon, I., Moliner, V., Eds.; RSC Publishing: Cambridge, UK, 2016. (38) Wang, F.; Landau, D. P. Efficient, multiple-range random walk algorithm to calculate the density of states. Phys. Rev. Lett. 2001, 86, 2050. (39) Wang, F.; Landau, D. Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64, 056101. (40) Bachmann, M.; Janke, W. Thermodynamics of lattice heteropolymers. J. Chem. Phys. 2004, 120, 6779−6791. (41) Taylor, M. P.; Paul, W.; Binder, K. Phase transitions of a single polymer chain: A Wang-Landau simulation study. J. Chem. Phys. 2009, 131, 114907. (42) Schnabel, S.; Vogel, T.; Bachmann, M.; Janke, W. Surface effects in the crystallization process of elastic flexible polymers. Chem. Phys. Lett. 2009, 476, 201−204. (43) Swetnam, A. D.; Allen, M. P. Improving the Wang-Landau algorithm for polymers and proteins. J. Comput. Chem. 2011, 32, 816− 821. (44) Belardinelli, R.; Pereyra, V. Fast algorithm to calculate density of states. Phys. Rev. E 2007, 75, 046701. (45) Junghans, C.; Perez, D.; Vogel, T. Molecular dynamics in the multicanonical ensemble: Equivalence of Wang-Landau sampling, statistical temperature molecular dynamics, and metadynamics. J. Chem. Theory Comput. 2014, 10, 1843−1847. (46) Smiatek, J.; Heuer, A. Calculation of free energy landscapes: A histogram reweighted metadynamics approach. J. Comput. Chem. 2011, 32, 2084−2096. (47) Yan, Q.; Faller, R.; de Pablo, J. J. Density-of-states Monte Carlo method for simulation of fluids. J. Chem. Phys. 2002, 116, 8745−8749. (48) Shell, M. S.; Debenedetti, P. G.; Panagiotopoulos, A. Z. Generalization of the Wang-Landau method for off-lattice simulations. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2002, 66, 056703. (49) Seaton, D.; Mitchell, S.; Landau, D. Developments in WangLandau simulations of a simple continuous homopolymer. Braz. J. Phys. 2008, 38, 48−53. (50) Cunha Netto, A.; Silva, C.; Caparica, A.; Dickman, R. WangLandau sampling in three-dimensional polymers. Braz. J. Phys. 2006, 36, 619−622. (51) Seaton, D.; Wüst, T.; Landau, D. Collapse transitions in a flexible homopolymer chain: Application of the Wang-Landau algorithm. Phys. Rev. E 2010, 81, 011802. (52) Rathore, N.; Knotts, T. A., IV; de Pablo, J. J. Density of states simulations of proteins. J. Chem. Phys. 2003, 118, 4285−4290. (53) Calvo, F. Sampling along reaction coordinates with the WangLandau method. Mol. Phys. 2002, 100, 3421−3427. (54) Smiatek, J.; Janssen-Müller, D.; Friedrich, R.; Heuer, A. Systematic detection of hidden complexities in the unfolding mechanism of a cytosine-rich DNA strand. Phys. A 2014, 394, 136− 144. (55) Deserno, M.; Holm, C.; May, S. Fraction of condensed counterions around a charged rod: Comparison of Poisson-Boltzmann theory and computer simulations. Macromolecules 2000, 33, 199−206. (56) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237−5247. (57) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. J. Chem. Phys. 1990, 92, 5057−5086. (58) Arnold, A.; Holm, C. MMM1D: A method for calculating electrostatic interactions in one-dimensional periodic geometries. J. Chem. Phys. 2005, 123, 144103.

(14) Uyaver, S.; Seidel, C. Effect of varying salt concentration on the behavior of weak polyelectrolytes in a poor solvent. Macromolecules 2009, 42, 1352−1361. (15) Atkins, P. W.; de Paula, J. Physical Chemistry; Oxford Univ. Press: Oxford, UK, 2010. (16) Polyelectrolyte Complexes in the Dispersed and Solid State II: Application Aspects; Müller, M., Ed.; Adv. Polym. Sci.; Springer: Berlin, 2013; Vol. 256. (17) Harnung, S. E.; Johnson, M. S. Chemistry and the environment; Cambridge Univ. Press: Cambridge, UK, 2012. (18) Butler, J. N. Ionic equilibrium: solubility and pH calculations; John Wiley & Sons: New York, NY, USA, 1998. (19) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: a reactive force field for hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (20) Åqvist, J.; Warshel, A. Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches. Chem. Rev. 1993, 93, 2523−2544. (21) Baptista, A. M.; Martel, P. J.; Petersen, S. B. Simulation of protein conformational freedom as a function of pH: constant-pH molecular dynamics using implicit titration. Proteins: Struct., Funct., Genet. 1997, 27, 523−544. (22) Börjesson, U.; Hünenberger, P. H. Explicit-solvent molecular dynamics simulation at constant pH: methodology and application to small amines. J. Chem. Phys. 2001, 114, 9706−9719. (23) Bürgi, R.; Kollman, P. A.; van Gunsteren, W. F. Simulating proteins at constant pH: an approach combining molecular dynamics and Monte Carlo simulation. Proteins: Struct., Funct., Genet. 2002, 47, 469−480. (24) Baptista, A. M.; Teixeira, V. H.; Soares, C. M. Constant-pH molecular dynamics using stochastic titration. J. Chem. Phys. 2002, 117, 4184−4200. (25) Walczak, A. M.; Antosiewicz, J. M. Langevin dynamics of proteins at constant pH. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2002, 66, 051911. (26) Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem. 2004, 25, 2038−2048. (27) Lee, M. S.; Salsbury, F. R.; Brooks, C. L. Constant-pH molecular dynamics using continuous titration coordinates. Proteins: Struct., Funct., Genet. 2004, 56, 738−752. (28) Machuqueiro, M.; Baptista, A. M. Constant-pH molecular dynamics with ionic strength effects: protonation-conformation coupling in decalysine. J. Phys. Chem. B 2006, 110, 2927−2933. (29) Stern, H. A. Molecular simulation with variable protonation states at constant pH. J. Chem. Phys. 2007, 126, 164112. (30) Goh, G. B.; Hulbert, B. S.; Zhou, H.; Brooks, C. L. Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism. Proteins: Struct., Funct., Genet. 2014, 82, 1319−1331. (31) Smith, W.; Triska, B. The reaction ensemble method for the computer simulation of chemical and phase equilibria. I. Theory and basic examples. J. Chem. Phys. 1994, 100, 3019−3027. (32) Heath Turner, C.; Brennan, J. K.; Lisal, M.; Smith, W. R.; Karl Johnson, J.; Gubbins, K. E. Simulation of chemical reaction equilibria by the reaction ensemble Monte Carlo method: a review. Mol. Simul. 2008, 34, 119−146. (33) Lísal, M.; Brennan, J. K.; Smith, W. R. Mesoscale simulation of polymer reaction equilibrium: Combining dissipative particle dynamics with reaction ensemble Monte Carlo. I. Polydispersed polymer systems. J. Chem. Phys. 2006, 125, 164905. (34) Košovan, P.; Limpouchová, Z.; Procházka, K. Charge Distribution and Conformations of Weak Polyelectrolyte Chains in Poor Solvents. Collect. Czech. Chem. Commun. 2008, 73, 439−458. (35) Lísal, M.; Brennan, J. K.; Smith, W. R. Mesoscale simulation of polymer reaction equilibrium: Combining dissipative particle dynamics with reaction ensemble Monte Carlo. II. Supramolecular diblock copolymers. J. Chem. Phys. 2009, 130, 104902. J


Article

Journal of Chemical Theory and Computation (59) Hockney, R.; Eastwood, J. Computer Simulation Using Particles; McGraw-Hill: New York, USA, 1981. (60) Limbach, H.-J.; Arnold, A.; Mann, B. A.; Holm, C. ESPResSo an extensible simulation package for research on soft matter systems. Comput. Phys. Commun. 2006, 174, 704−727. (61) Arnold, A.; Lenz, O.; Kesselheim, S.; Weeber, R.; Fahrenberger, F.; Roehm, D.; Košovan, P.; Holm, C. In Meshfree methods for partial differential equations VI; Schweitzer, M. A., Ed.; Lecture Notes in Computational Science and Engineering; Springer: Berlin, 2013; Vol. 89; pp 1−23. (62) Fuoss, R. M.; Katchalsky, A.; Lifson, S. The potential of an infinite rod-like molecule and the distribution of the counter ions. Proc. Natl. Acad. Sci. U. S. A. 1951, 37, 579−589. (63) Alfrey, T.; Berg, P. W.; Morawetz, H. The counterion distribution in solutions of rod-shaped polyelectrolytes. J. Polym. Sci. 1951, 7, 543−547. (64) Naji, A.; Netz, R. R. Scaling and universality in the counterioncondensation transition at charged cylinders. Phys. Rev. E 2006, 73, 056105. (65) Parsons, D. F.; Williams, D. R. An off-lattice Wang-Landau study of the coil-globule and melting transitions of a flexible homopolymer. J. Chem. Phys. 2006, 124, 221103. (66) Seaton, D.; Wüst, T.; Landau, D. A Wang−Landau study of the phase transitions in a flexible homopolymer. Comput. Phys. Commun. 2009, 180, 587−589. (67) Wüst, T.; Li, Y. W.; Landau, D. P. Unraveling the beautiful complexity of simple lattice model polymers and proteins using WangLandau sampling. J. Stat. Phys. 2011, 144, 638−651. (68) Smiatek, J.; Wohlfarth, A.; Holm, C. The solvation and ion condensation properties for sulfonated polyelectrolytes in different solvents-a computational study. New J. Phys. 2014, 16, 025001. (69) Panagiotopoulos, A. Charge correlation effects on ionization of weak polyelectrolytes. J. Phys.: Condens. Matter 2009, 21, 424113. (70) Limbach, H. J.; Holm, C. Single-chain properties of polyelectrolytes in poor solvent. J. Phys. Chem. B 2003, 107, 8041− 8055. (71) Kundagrami, A.; Muthukumar, M. Effective Charge and CoilGlobule Transition of a Polyelectrolyte Chain. Macromolecules 2010, 43, 2574−2581. (72) Brilliantov, N.; Kuznetsov, D.; Klein, R. Chain collapse and counterion condensation in dilute polyelectrolyte solutions. Phys. Rev. Lett. 1998, 81, 1433. (73) Longo, G. S.; Olvera de La Cruz, M.; Szleifer, I. Molecular theory of weak polyelectrolyte gels: the role of pH and salt concentration. Macromolecules 2011, 44, 147−158. (74) Leonforte, F.; Müller, M. Functional Poly (N-isopropylacrylamide)/Poly (acrylic acid) Mixed Brushes for Controlled Manipulation of Nanoparticles. Macromolecules 2016, 49, 5256−5265. (75) Smiatek, J.; Heuer, A.; Wagner, H.; Studer, A.; Hentschel, C.; Chi, L. Coat thickness dependent adsorption of hydrophobic molecules at polymer brushes. J. Chem. Phys. 2013, 138, 044904. (76) Hentschel, C.; Wagner, H.; Smiatek, J.; Heuer, A.; Fuchs, H.; Zhang, X.; Studer, A.; Chi, L. AFM-based force spectroscopy on polystyrene brushes: effect of brush thickness on protein adsorption. Langmuir 2013, 29, 1850−1856. (77) Bereau, T.; Deserno, M. Generic coarse-grained model for protein folding and aggregation. J. Chem. Phys. 2009, 130, 235106. (78) Smiatek, J.; Chen, C.; Liu, D.; Heuer, A. Stable conformations of a single stranded deprotonated DNA i-motif. J. Phys. Chem. B 2011, 115, 13788−13795. (79) Smiatek, J.; Heuer, A. Deprotonation mechanism of a singlestranded DNA i-motif. RSC Adv. 2014, 4, 17110−17113. (80) Raju, M.; Kim, S.-Y.; van Duin, A. C. T.; Fichthorn, K. A. ReaxFF Reactive Force Field Study of the Dissociation of Water on Titania Surfaces. J. Phys. Chem. C 2013, 117, 10558−10572.

K


WangâLandau Reaction Ensemble Method - ACS Publications

Recommend Documents

WangâLandau Reaction Ensemble Method - ACS Publications