Equation of State of Charged Rod Dispersions

We study the accuracy of the theory of Stroobants, Lekkerkerker, and Odijk [Macromolecules1986, 19, 2232–2238], called SLO theory, to describe the thermodynamic properties of an isotropic fluid of charged rods. By incorporation of the effective diameter of the rods according to SLO theory into scaled particle theory (SPT), we obtain an expression for the rod concentration-dependent free volume fraction and the osmotic pressure of a collection of charged hard spherocylinders. The results are compared to Monte Carlo simulations. We find close agreement between the simulation results and the SLO-SPT predictions for not too large values of the Debye length and for high rod charge densities. The deviations increase with rod density, particularly at concentrations above which isotropic–nematic phase transitions are expected.


■ INTRODUCTION
Dispersions of rodlike particles are of great interest because of the relatively large excluded volume between rods.−19 Understanding solutions of rodlike macromolecules is also of biological interest: think of the dense packing of rodlike DNA molecules in virus heads 20,21 or in the organization of cells. 22,23Micro-organisms sometimes assume a spherocylinder shape, such as Escherichia coli, 24 Campylobacter fetus, Salmonella typhimurium, and Mycobacterium tuberculosis bacteria.Assuming rodlike morphologies is considered to be a tool to gain a competitive advantage. 25,26odlike colloids often assume charges when they are dispersed in a polar solvent.CNCs are rodlike and are charged when dispersed in aqueous solution. 27,28Several types of viruses such as tobacco mosaic virus (TMV) 29−31 or the filamentous bacteriophages fd virus, 32,33 and its variants 34,35 like Pf1 and Pf4, are long and thin rods, assuming charged surfaces in water.In aqueous salt solutions charged rods are surrounded by double layers with an inhomogeneous distribution of co-and counterions. 36The screening length of the extent of double layers is often described by using the Debye length, while the magnitude of the inhomogeneity is determined by the surface charge density at the rod surface.The resulting double-layer forces between charged colloidal rods in a polar solvent are specified by the range and the strength of the repulsive interaction, in a similar fashion as for colloidal spheres. 37Because the ionic strength affects the range of the double-layer repulsion between the rods, it influences the rod concentration at which the dispersion undergoes phase transitions.
Fraden and others 32,38−41 showed that the salt concentration affects the isotropic−nematic phase transition of tobacco mosaic 30 and fd 42 viruses.It has also been shown that the phase transition concentration of CNC dispersions depends on the ionic strength. 43o theoretically predict such phase transitions, Onsager 2 proposed to describe charged rods as hard rods with an ef fective diameter D eff > D, to mimic the hard core plus soft repulsion of charged rods.The Stroobants−Lekkerkerker−Odijk (SLO) approach 44 nowadays is a standard method to estimate D eff .From this theory, it follows that D eff becomes larger when the Debye length and/or the surface charge are larger.
Although the SLO theory is applied frequently, its accuracy has, as far as we are aware, not been evaluated in detail.Vologodskii and Cozzarelli 45 tested different approaches for the electrostatic interactions in closed DNA chains modeled as chains of connected straight segments.They found good agreement between the Debye−Huckel theory and the effective diameter model for the conformational properties of interest.In this paper, we verify the accuracy of SLO theory against computer simulations.We first show that determining the free volume fraction available directly provides the osmotic pressure and thereby the equation of state.Subsequently, we analyze the accuracy of this approach.

Theory. Connection between α and the Equation of State.
The volume available in a system of interest is an important thermodynamic quantity.From Widom's particle insertion theorem 46 it follows that the relative volume available for a particle upon insertion into a dispersion of particles, called the free volume fraction α, is directly related to the immersion free energy W. The excess chemical potential of the particle μ ex = k B Tμẽ x = μ − μ 0 − k B T ln ϕ, where μ 0 is the (reference) chemical potential and ϕ is the particle volume fraction, is related to W through Applying the Gibbs−Duhem relation to eq 1 yields 3)   which provides a direct relation 47 between the free volume fraction for a particle in a system and the osmotic pressure Π ̃, hence the equation of state.We apply this to the cases of hard rods and charged hard rods.Rods are described as spherocylinders, which are cylinders (with hard core volume v r ) with length L capped by hemispheres with diameter D. It is assumed the spherocylinders only adopt isotropic configurations for the range of rod volume fractions studied here. 1,13ressure of Hard Spherocylinders: SPT Prediction.We use scaled particle theory (SPT) to quantify the volume fraction ϕ r dependence of the osmotic pressure Π of a fluid of hard spherocylinders (HSCs):  = and y = ϕ r /(1 − ϕ r ).The dimensionless parameter γ is related to the aspect ratio of the rods L/D as γ = 1 + L/D.This result for the osmotic pressure accurately describes computer simulations of hard spherocylinders 1,49 (see ref 13).The free volume fraction α available to a HSC can be determined with SPT and reads 50 Ä where the last expression results from the insertion of eq 4 into eq 5a.
Pressure of Charged Spherocylinders.We now consider charged rods dispersed in a polar solvent that interact through a double-layer repulsion next to the hard-core excluded-volume interaction.This double-layer interaction gives rise to an orientation-dependent soft repulsive interaction between the rods.The density of surface charge groups, which is directly related to the electrostatic surface potential Ψ at the rod surface, determines the strength of the repulsion.The ionic strength of the medium dictates the Debye length, which mediates the range of the double-layer repulsion. 36Stroobants, Lekkerkerker, and Odijk 44 used the following expression for the pair interaction U(x,θ) between two charged rods: where κ −1 is the ionic-strength-dependent Debye screening length, x is the closest distance between the center lines of the rods, and θ is the angle between the rods (cf. Figure 1).
Incorporation of eq 6 into the Onsager free energy of infinitely long rods reveals that an effective rod diameter is given by 44 where k E is Euler's constant ≈ 0.5772156649 and A′ = Ae −κD follows from )   with λ B the Bjerrum length.The parameter ζ is the proportionality constant of the outer part of the double-layer electrostatic potential profile near a charged rod: 51 The Journal of Physical Chemistry B where e is the elementary charge, r is the distance from the center line of the rod, and K 0 is the modified second kind of order 0 Bessel function.For a weakly charged rod the Debye−Huckel approximation provides 52 where Z is the linear charge density (per unit length) of the rod and K 1 is the modified Bessel function of the second kind of order 1.Comparison of eqs 9 and 10 provides an expression for ζ which, after insertion into eq 8, gives For thick and thin double layers the following asymptotic analytical expressions have been derived: 53 We found that eq 11 can be approximated accurately for arbitrary double-layer thickness as The accuracy of this approximation is visualized in Figure 2. The relative difference between the approximation and the exact expression is below 5% for the whole range of κD (cf. Figure S1).
The expression for the osmotic pressure in a system of charged rods is found by insertion of D eff for D into eq 4: Computer Simulations.In computer simulations, one cannot easily deal with infinitely ranged interactions and, thus, has to define a maximum range, up to which pair interactions are taken into account.To avoid further spurious effects, the interaction potential is also shifted so that it is zero when the distance between particles equals the defined range.This results in the following adaptation of the interaction potential where d c is the cutoff length.
To estimate the effective free volume fraction α eff available for the rods and the immersion free energy W, we applied test particle insertions.To test the relation between α eff and W given in eq 1, we compute , where the brackets denote an average over n t trial insertions and n c independent configurations.U c is the system's energy difference due to the particle insertion and an insertion is done by randomly choosing a position and an orientation for the inserted rod (which is removed again after the trial).For the effective free volume fraction α eff , the relative number of successful trial insertions, n s /n t , is computed, where successful means that there is no overlap between the inserted rod and any other rod in the system assuming that the rods have diameter D eff .The quantity α eff is then approximated by the average ⟨n s /n t ⟩ nd c .
The systems contain N r = 1352 rods and have periodic boundary conditions in all three dimensions.For each set (A, κ, d c ), we initialize the system at a density high enough to have a very small effective free volume fraction and with all rods pointing in the same direction.The system then relaxes to its equilibrium state, in which subsequently the measurements are done.The simulation boxes of equilibrated states are expanded to reduce the rod volume fraction and, thus, obtain the density dependence of the observables (after re-equilibrating at the new density).The simulation method is Metropolis Monte Carlo, where single particle translations and rotations are used as trial moves, which are accepted with probability min(1, e ) U c with the energy difference before and after the move, U c .The acceptance rate for the trial moves is adjusted to be between 0.3 and 0.4 by changing the maximum displacement and rotation (unless the density is too low, which resulted in overall high acceptance rates).

■ RESULTS AND DISCUSSION
In this section we first show the accuracy of SPT to predict the free volume fraction and osmotic pressure using eqs 3, 4, and 5 for hard spherocylinders by comparing the results with computer simulations.Because the simulations can only be performed by defining a cutoff length of the range of the soft interaction of eq 7, we next study its proper value.Subsequently, we compare Monte Carlo simulation results for the free volume fractions of a dispersion of rods with an effective diameter against that of rods with a hard core but soft repulsion, following eq 7. Finally, we compare the theoretical predictions of the equation of state from theory (eq 14) with computer simulations.

The Journal of Physical Chemistry B
Benchmarking the Approach for Hard Spherocylinders.First we evaluated the accuracy of our method for hard spherocylinders for various values of γ.In a previous study, 54 we already confirmed the accuracy of eq 5a for a specific aspect ratio of γ = 6.In Figure 3a, the data points are Monte Carlo computer simulation results for the free volume fraction, and the curves are predictions of eq 5b.These MC simulation data points were used to compute the osmotic pressure using the exact relation eq 3, and those results are the data in Figure 3b.It follows that SPT correctly predicts the free volume fraction (eq 5) and osmotic pressure (eq 4) of hard spherocylinders in the isotropic phase state for a wide range of aspect ratios (cf.Figures 3a and 3b  )   where the exponential integral is defined through Ei(z) = ∫ −∞ z e t /t dt.For d c → ∞, eq 17 becomes In this case, the relative difference between the approximation, eq 7, and the exact value of the effective diameter D eff , eq 16, is less than 2.2% for A′ ≥ 2. For a finite cutoff length, one can determine the value of d c , for which the relative difference between the effective diameter using the approximation of an infinite interaction range and the exact effective diameter is also less than 2.2%.This threshold cutoff length, d c t , decreases with κ and increases with A. We provide some exemplary values for d c t / D in Table 1.Thus, if d c > d c t and A′ ≥ 2, eq 7 gives a very close approximation for the effective diameter; if d c > d c t and A′ < 2, D eff can be approximated by using eq 18 for evaluating eq 16; and if d c ≤ d c t , the exact value from using eq 17 should be applied.Comparison of ⟨e −U ̃c⟩ and the Effective Free Volume Fraction.The effective free volume fraction is less expensive to compute compared with the true immersion free energy.Thus, it would be a great advantage, if the D eff approach gives a good approximation of the immersion free energy.Figure 4 (and Figures S2 and S3) shows the relative difference between α eff and α ref for different sets (A, κD) and different cutoff lengths versus the effective volume fraction of rods, The difference varies greatly depending on the specific set of parameters.It increases with the rod density and with the cutoff length, which can be explained as follows: The rods in the system can have a smaller shortest distance of their axes than the effective diameter, whereas the inserted rod only fits when its shortest distance to all other rods is at least D eff .That means that if the rods truly had a hard diameter of D eff , the respective free volume would be smaller.This effect increases with the density because at higher density more rods have a shortest distance between D and D eff .And the effect increases when the effective diameter is larger, which is the case for increasing d c .Because the estimation of the effective diameter stems from a second virial approximation, it is not too surprising that the accuracy breaks down for denser systems.
We have assumed that the effective rod diameter is independent of the charged rod density.When it comes to The Journal of Physical Chemistry B studying the isotropic−nematic transition and other phase equilibria at higher rod concentrations, this may lead to deviations.In practice, the assumption that eq 6 is independent of the rod concentration will become inaccurate as the charged rods themselves will also contribute to the effective value of the screening length κ −1 .In addition, also the twisting effect 44,55,56 needs to be taken into account as soon as the rods undergo phase transitions.To improve the D eff approach, one could, for example, follow the avoidance model for soft particles that results in a reduction of the effective diameter with increasing concentration of the particles. 57,58n order to further quantify the accuracy of the D eff approach for different combinations of A, κD, and d c , we do the following.(i) We interpolate the data for α eff (φ r eff ) and α ref (φ r eff ).(ii) We find the effective rod volume fraction (φ r eff,0.3 ) for which the reference free volume fraction value equals α ref = 0.3.(iii) We compute the relative difference between the effective free volume fraction at this rod volume fraction and the reference The decreasing difference seen for some of the largest rod volume fractions is due to both free volume fraction estimates getting close to zero.)The Journal of Physical Chemistry B free volume fraction value of 0.3. 59The results are shown in Figure 5 (and in Figures S4 and S5).
Figures 5a and 5b provide an overview of the deviations of the D eff approach from the results obtained using soft repulsive rods.It also identifies (Figures 5c and 5d) the range of D eff /D values over which the D eff can be accurately applied for some illustrative (A, κD) values.

Comparison of Theoretical Predictions and Simulations.
In Figure 6 we directly compare the theoretical predictions (SLO + SPT) to the simulation results for the effective free volume fraction and the resulting osmotic pressure.The agreement is almost perfect for all studied aspect ratios, cutoff lengths, and interaction parameters.Only when the effective free volume fraction goes to very small values (i.e., high effective rod concentrations) does the relative difference  The effective diameter ratio D eff /D is indicated in the legend.Points are shown in gray when the relative difference between the effective free volume fraction and the reference from the immersion free energy is larger than 2.2%, i.e., when there already is a systematic error in using the effective free volume to study thermodynamic properties.

The Journal of Physical Chemistry B
between simulation and theory start to increase (cf. Figure 7).The conclusion from Figure 7 is that when the system volume fraction is small enough so that the effective free volume fraction gives a good estimate of the immersion free energy, then the theoretically predicted value for the effective free volume fraction is very close to the value obtained from simulations.The deviation between theory and simulation for larger volume fractions is again related to overlaps of the effective shells; the theory assumes no overlap, whereas there can be overlaps in the simulated system so that the actual volume fractions differ.

■ CONCLUSIONS
In this paper, we have studied the equation of state of charged rods.A frequently used theory for the thermodynamic properties of charged rods is that by Odijk, Stroobants, and Lekkerkerker. 44his SLO theory showed that the free energy for infinitely long rods interacting through a hard-core double-layer potential can be described via the effective diameter of the rods.Here we verified the accuracy of this approach for finite sized rods.
We first showed that scaled particle theory (SPT) provides accurate results for the concentration and aspect ratio dependence of rods in the isotropic phase by a comparison of the free volume fraction and osmotic pressure of hard rods.Next, we extended this approach by incorporating SLO theory into the SPT to predict the osmotic pressure for charged rods and compare the results with computer simulations.
It is found that the accuracy of using the effective diameter following SLO theory depends on the Debye length, charge density of the rods, and the rod concentration.Deviations decrease when the rod charge density is higher and the Debye length is smaller.For denser rods, the deviations increase.It is realized, however, that for dense rods phase transitions are expected to nematic, smectic, or crystalline phases where rodmediated screening will reduce the effective diameter.We have the ambition to explore these phase transitions of charged rods in a subsequent work.
Relative difference between the approximation and the exact expression for A′ and some computer simulation results for other rod aspect ratios as those presented in the main text (PDF) ■

Figure 1 .
Figure 1.Sketch to indicate the closest distance x and the angle θ between two rods.The blue shape shows the hard spherocylinders (no overlap allowed), while the black solid lines mark their central axes which are decisive for the pair interaction (eq 6).
with γ* = 1 + L/D eff and y eff = φ r eff /(1 − φ r eff ), where φ r eff is the effective volume fraction of rods with diameter D eff .

Figure 2 .
Figure 2. Comparison of approximations for A DH ′ : exact expression (eq 11), proposed approximation (eq 13), and often used approximations for small and large κD.
). Effective Diameter for a Finite Ranged Interaction.Equation 7 is valid for A′ ≳ 2 and d c → ∞ (see ref 44).The exact result for the effective diameter and an arbitrary cutoff length d c can be found by (numerically) evaluating the following integral: distribution function in the isotropic state, f(Ω) = 1/(4π), and the scaled electrostatic excluded volume

Figure 3 .
Figure 3. Theoretical SPT predictions (curves) of (a) eq 5 and (b) eq 4 compared to simulation results (symbols) of the free volume fraction (a) and the osmotic pressure (b) of (uncharged) hard spherocylinders vs rod volume fraction for different rod aspect ratios γ = 1 + L/D as indicated in the legend.

Table 1 .
Examples for the Threshold of the Cutoff Length, d c T /D, beyond Which the Relative Difference between the Exact Result and the Infinite Range Approximation for the Effective Diameter Is Smaller than 2.2% (the Corresponding Effective Diameter, D eff /D, Is Given in Square Brackets)

Figure 4 .
Figure 4. Relative difference between the effective free volume fraction and the reference from the immersion free energy vs effective volume fraction for different cutoff lengths as indicated in the legend and (A, κD) = (1.0,0.5) (a), (16.0, 0.5) (b), (1.0, 2.0) (c), and (16.0, 2.0) (d).The hard rod aspect ratio is L/D + 1 = 11, and the effective diameter ratio D eff /D is indicated in the legend.The inset in (a) shows the whole range of relative differences, while the main figures show the same range for (a−d).(Note that we only show data with at least 10 successful insertions during the simulation run.The decreasing difference seen for some of the largest rod volume fractions is due to both free volume fraction estimates getting close to zero.)

Figure 5 .
Figure 5. (a, b) Relative difference between the effective free volume fraction and the reference from the immersion free energy at the effective volume fraction for which the reference equals 0.3 vs effective diameter for different A and κD.Identical symbols refer to the same (A, κD) but different cutoff lengths (resulting in different effective diameter ratios).The hard rod aspect ratio is L/D + 1 = 11.(a) Full range of effective diameters.(b) Small effective diameters.(c, d) Accessible range of effective diameter ratios for given (A, κD) and varying d c ≥ 2.

Figure 6 .
Figure 6.Theoretical prediction and simulation result of the effective free volume fraction (a, b) and the osmotic pressure (c, d) vs effective rod volume fraction for different cutoff lengths as indicated in the legend and (A, κD) = (1.0,2.0) (a, c) and (16.0, 0.5) (b, d).The hard rod aspect ratio is L/D + 1 = 11, and the effective diameter ratio D eff /D is indicated in the legend.Points in(a, b) are shown in gray when the relative difference between the effective free volume fraction and the reference from the immersion free energy is larger than 2.2%, i.e., when there already is a systematic error in using the effective free volume to study thermodynamic properties.

Figure 7 .
Figure 7. Relative difference between the theoretical and simulation value of the effective free volume fraction vs effective rod volume fraction for different cutoff lengths as indicated in the legend and (A, κD, L/D) = (1.0,2.0, 5) (a), (16.0, 0.5, 5) (b), (1.0, 2.0, 10) (c), and (16.0, 0.5, 10) (d).The effective diameter ratio D eff /D is indicated in the legend.Points are shown in gray when the relative difference between the effective free volume fraction and the reference from the immersion free energy is larger than 2.2%, i.e., when there already is a systematic error in using the effective free volume to study thermodynamic properties.