Titration of hydrophobic polyelectrolytes using Monte Carlo simulations

The conformation and titration curves of weak (or annealed) hydrophobic polyelectrolytes have been examined using Monte Carlo simulations with screened Coulomb potentials in the grand canonical ensemble. The influence of the ionic concentration pH and presence of hydrophobic interactions has been systematically investigated. A large number of conformations such as extended, pearl-necklace, cigar-shape, and collapsed structures resulting from the subtle balance of short-range hydrophobic attractive interactions and long-range electrostatic repulsive interactions between the monomers have been observed. Titration curves were calculated by adjusting the pH-pK(0) values (pK(0) represents the intrinsic dissociation constant of an isolated monomer) and then calculating the ionization degree alpha of the polyelectrolyte. Important transitions related to cascades of conformational changes were observed in the titration curves, mainly at low ionic concentration and with the presence of strong hydrophobic interactions. We demonstrated that the presence of hydrophobic interactions plays an important role in the acid-base properties of a polyelectrolyte in promoting the formation of compact conformations and hence decreasing the polyelectrolyte degree of ionization for a given pH-pK(0) value.     

Molecular organization in liquids consisting of anisometric molecules: Monte Carlo simulations and comparison with theories

Monte Carlo computer simulations of systems of rigid anisometric molecules with axial ratio x 3, 4 and 5 have been performed at densities comparable to those of the mesogenic groups in liquid crystalline materials. The agreement between computer experiments and Flory-Ronca-Irvine theory is quite good. In particular, the nematic-isotropic transition temperatures found in the simulations are within 20 K of those predicted by the theory in all cases examined with x 3 and 4, and only slightly higher when x 5. Furthermore, the orientational distributions of the molecular long axes with respect to the nematic director are always nearly coincident with those evaluated according to the FRI theory, provided that the order parameter is the same. However, the order parameters at the transition points are smaller than expected on the basis of the theory, and closer to experimental values. Calculations performed by changing the size of the simulated systems prove that the reported results are size independent, while size dependent results are obtained with smaller systems.     

Microscopic solvation: spectroscopic results vs. Monte Carlo simulations

Van-der-Waals clusters of carbazole (representing the ‘solute’) with up to 40 nitrogen or methane solvent molecules were characterized using two-color resonant two-photon ionization spectroscopy. Features in these spectra (redshift, homogeneous and heterogeneous broadening, etc.) are interpreted as being caused by various static and dynamic effects of the solvent shell surrounding the aromatic substrate. For a better understanding of such effects, Monte Carlo simulations of the clusters were performed: Statics: Using a Monte Carlo simulated annealing minimization procedure, minimum energy structures (local, probably global) were found for the various cluster species. Using a simple empirical additivity rule, spectral shifts are rationalized from these structures.Dynamics: Starting from these minimum configurations, canonical ensemble simulations were carried out in a temperature range from 0 to 35 K. Severalorder-disorder transitions were identified including:

1. orientational isomerization or ‘melting’
2. surface isomerization or decoupling
3. rigid → fluxional transitions
4. full cluster ‘melting’

We present some of our experimental results on the systems carbazole · (N₂) _n and carbazole · (CH₄) _n together with the corresponding simulation data.     

Monte Carlo simulations of liquid crystals near rough walls

The effect of surface roughness on the structure of liquid crystalline fluids near solid substrates is studied by Monte Carlo simulations. The liquid crystal is modeled as a fluid of soft ellipsoidal molecules and the substrate is modeled as a hard wall that excludes the centers of mass of the fluid molecules. Surface roughness is introduced by embedding a number of molecules with random positions and orientations within the wall. It is found that the density and order near the wall are reduced as the wall becomes rougher, i.e., the number of embedded molecules is increased). Anchoring coefficients are determined from fluctuations in the reciprocal space order tensor. It is found that the anchoring strength decreases with increasing surface roughness.     

Monte Carlo simulations of the hydrophobic effect in aqueous electrolyte solutions

The hydrophobic interaction between two methane molecules in salt-free and high salt-containing aqueous solutions and the structure in such solutions have been investigated using an atomistic model solved by Monte Carlo simulations. Monovalent salt representing NaCl and divalent salt with the same nonelectrostatic properties as the monovalent salt have been used to examine the influence of the valence of the salt species. In salt-free solution the effective interaction between the two methane molecules displayed a global minimum at close contact of the two methane molecules and a solvent-separated secondary minimum. In 3 and 5 M monovalent salt solution the potential of mean force became slightly more attractive, and in a 3 M divalent salt solution the attraction became considerably stronger. The structure of the aqueous solutions was determined by radial distribution functions and angular probability functions. The distortion of the native water structure increased with ion valence. The increase of the hydrophobic attraction was associated with (i) a breakdown of the tetrahedral structure formed by neighboring water molecules and of the hydrogen bonds between them and (i) the concomitant increase of the solution density.     

Ion binding to polyelectrolytes: Monte Carlo simulations versus classical mean field theories

The influence of ion size and surface charge model in titrations of ionizable polyelectrolytes is studied by means of the Semi Grand Canonical Monte Carlo simulation method in the context of the primitive model. Three models describing a discrete distribution of charged functional groups on the polyelectrolyte and different values for the radius of the background electrolyte spanning from ionic to hydrated radii values were analyzed. The polyelectrolyte titrations were simulated by calculating the degree of ionization versus pH curves at two ionic strengths. The results allow us to quantify the impact of the sizes of the background salt ions and surface functional groups of the polyelectrolyte on the dissociation degree. This influence is explained in terms of the effectiveness of the screening of the charged surface sites. Finally, by comparison with the Non-Linear Poisson–Boltzmann model, the influence of ionic correlations and finite size of the solution ions is assessed.     

Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations

An efficient Monte Carlo (MC) algorithm including concerted rotations is directly compared to molecular dynamics (MD) in all-atom statistical mechanics folding simulations of small polypeptides. The previously reported algorithm "concerted rotations with flexible bond angles" (CRA) has been shown to successfully locate the native state of small polypeptides. In this study, the folding of three small polypeptides (trpzip2/H1/Trp-cage) is investigated using MC and MD, for a combined sampling time of approximately 10(11) MC configurations and 8 micros, respectively. Both methods successfully locate the experimentally determined native states of the three systems, but they do so at different speed, with 2-2.5 times faster folding of the MC runs. The comparison reveals that thermodynamic and dynamic properties can reliably be obtained by both and that results from folding simulations do not depend on the algorithm used. Similar to previous comparisons of MC and MD, it is found that one MD integration step of 2 fs corresponds to one MC scan, revealing the good sampling of MC. The simplicity and efficiency of the MC method will enable its future use in folding studies involving larger systems and the combination with replica exchange algorithms.     

Thermodynamic and structural properties of repulsive hard-core Yukawa fluid: integral equation theory, perturbation theory and Monte Carlo simulations

The thermodynamic and structural properties of purely repulsive hard-core Yukawa particles in the fluid state are determined through Monte Carlo simulation and modeled using perturbation theory and integral equation theory in the mean spherical approximation (MSA). Systems of particles with Yukawa screening lengths of 1.8, 3.0, and 5.0 are examined with results compared to variations of MSA and perturbation theory. Thermodynamic properties were predicted well by both theories in the fluid region up to the fluid-solid phase boundary. Further, we found that a simplified exponential version of the MSA is the most accurate at predicting radial distribution function at contact. Radial distribution function of repulsive hard-core Yukawa particles are also reported. The results show that methods based on MSA and perturbation theory that are typically applied to the attractive hard-core Yukawa potential can also be extended to the purely repulsive hard-core Yukawa potential.     

Quantum mechanical single molecule partition function from path integral Monte Carlo simulations

An algorithm for calculating the partition function of a molecule with the path integral Monte Carlo method is presented. Staged thermodynamic perturbation with respect to a reference harmonic potential is utilized to evaluate the ratio of partition functions. Parallel tempering and a new Monte Carlo estimator for the ratio of partition functions are implemented here to achieve well converged simulations that give an accuracy of 0.04 kcal/mol in the reported free energies. The method is applied to various test systems, including a catalytic system composed of 18 atoms. Absolute free energies calculated by this method lead to corrections as large as 2.6 kcal/mol at 300 K for some of the examples presented.     

Grand canonical Monte Carlo simulations of water in protein environments

The grand canonical simulation algorithm is considered as a general methodology to sample the configuration of water molecules confined within protein environments. First, the probability distribution of the number of water molecules and their configuration in a region of interest for biochemical simulations, such as the active site of a protein, is derived by considering a finite subvolume in open equilibrium with a large system serving as a bulk reservoir. It is shown that the influence of the bulk reservoir can be represented as a many-body potential of mean force acting on the atoms located inside the subvolume. The grand canonical Monte Carlo (GCMC) algorithm, augmented by a number of technical advances to increase the acceptance of insertion attempts, is implemented, and tested for simple systems. In particular, the method is illustrated in the case of a pure water box with periodic boundary conditions. In addition, finite spherical systems of pure water and containing a dialanine peptide, are simulated with GCMC while the influence of the surrounding infinite bulk is incorporated using the generalized solvent boundary potential [W. Im, S. Berneche, and B. Roux, J. Chem. Phys. 114, 2924 (2001)]. As a last illustration of water confined in the interior of a protein, the hydration of the central cavity of the KcsA potassium channel is simulated.     

Monte Carlo simulations of conformations of chain molecules in a cylindrical pore

Off-lattice Monte Carlo simulations are employed to study the behavior of chain molecules confined in a long cylindrical pore under the condition of hard-body interaction. The emphasis is placed on the chain and bead distributions as well as the location dependence of the chain conformation and anisotropy. The simulation results show that the chains very near the pore surface tend to wrap around the surface in various configurations. This behavior is qualitatively similar to that of the chains near but outside a cylindrical rod. Moreover, the bead concentration near the pore surface increases with increasing surface curvature.     

Spatially adaptive lattice coarse-grained Monte Carlo simulations for diffusion of interacting molecules

While lattice kinetic Monte Carlo (KMC) methods provide insight into numerous complex physical systems governed by interatomic interactions, they are limited to relatively short length and time scales. Recently introduced coarse-grained Monte Carlo (CGMC) simulations can reach much larger length and time scales at considerably lower computational cost. In this paper we extend the CGMC methods to spatially adaptive meshes for the case of surface diffusion (canonical ensemble). We introduce a systematic methodology to derive the transition probabilities for the coarse-grained diffusion process that ensure the correct dynamics and noise, give the correct continuum mesoscopic equations, and satisfy detailed balance. Substantial savings in CPU time are demonstrated compared to microscopic KMC while retaining high accuracy.     

Monte Carlo simulations of systems of rodlike molecules with semiflexible terminal groups

The results of Monte Carlo simulations of systems of rodlike molecules with semiflexible terminal groups are compared with those of systems with terminal groups at a fixed angle τ with respect to the rigid core and with those of a system containing the rigid cores alone in the same total volume. The nematic/isotropic transition temperature decreases with increasing flexibility or with increasing τ. For semiflexible groups, the transition is coupled with a conformational change favouring more extended conformations in the nematic phase. The results of the simulations are discussed in connection with those obtained for similar models in the recent literature and with the predictions of the Maier‐Saupe and Flory‐Ronca‐Irvine theories.     

Surface tension of water and acid gases from Monte Carlo simulations

We report direct Monte Carlo (MC) simulations on the liquid-vapor interfaces of pure water, carbon dioxide, and hydrogen sulfide. In the case of water, the recent TIP4P/2005 potential model used with the MC method is shown to reproduce the experimental surface tension and to accurately describe the coexistence curves. The agreement with experiments is also excellent for CO(2) and H(2)S with standard nonpolarizable models. The surface tensions are calculated by using the mechanical and the thermodynamic definitions via profiles along the direction normal to the surface. We also discuss the different contributions to the surface tension due to the repulsion-dispersion and electrostatic interactions. The different profiles of these contributions are proposed in the case of water.     

Structure of penetrable-rod fluids: exact properties and comparison between Monte Carlo simulations and two analytic theories

Bounded potentials are good models to represent the effective two-body interaction in some colloidal systems, such as the dilute solutions of polymer chains in good solvents. The simplest bounded potential is that of penetrable spheres, which takes a positive finite value if the two spheres are overlapped, being 0 otherwise. Even in the one-dimensional case, the penetrable-rod model is far from trivial, since interactions are not restricted to nearest neighbors and so its exact solution is not known. In this paper the structural properties of one-dimensional penetrable rods are studied. We first derive the exact correlation functions of the penetrable-rod fluids to second order in density at any temperature, as well as in the high-temperature and zero-temperature limits at any density. It is seen that, in contrast to what is generally believed, the Percus-Yevick equation does not yield the exact cavity function in the hard-rod limit. Next, two simple analytic theories are constructed: a high-temperature approximation based on the exact asymptotic behavior in the limit T--> infinity and a low-temperature approximation inspired by the exact result in the opposite limit T--> 0. Finally, we perform Monte Carlo simulations for a wide range of temperatures and densities to assess the validity of both theories. It is found that they complement each other quite well, exhibiting a good agreement with the simulation data within their respective domains of applicability and becoming practically equivalent on the borderline of those domains. A comparison with numerical solutions of the Percus-Yevick and the hypernetted-chain approximations is also carried out. Finally, a perspective on the extension of our two heuristic theories to the more realistic three-dimensional case is provided.     

Monte Carlo simulations of the polymer glass transition: From the test of theories to material modeling

We present results on the glass transition in polymer melts using Monte Carlo simulations of the bond fluctuation lattice model. There are two questions we address in this work. What is the temperature dependence of the entropy density in such a model polymer melt and how well is it described by theories like the Gibbs-DiMarzio theory of the glass transition? And to what degree is one able to map the Hamiltonian of such an abstract lattice model onto a specific polymer material and use it to model the large scale and long time properties of a realistic polymer melt?     

Thermodynamic properties of small rare gas clusters by path integral Monte Carlo simulations: a preliminary study

A strategy for reducing the risk of non-ergodic simulations in Monte Carlo calculations of the thermodynamic properties of clusters is discussed with the support of some examples. The results obtained attest the significance of the approach for the low-temperature regime, as non-ergodic sampling of potential energy surfaces is a particularly insidious occurrence. Fourier path integral Monte Carlo techniques for taking into account quantum effects are adopted, in conjunction with suitable tricks for improving the procedure reliability. Applications are restricted to Lennard-Jones clusters of rare-gas systems.     

Monte Carlo simulations of critical cluster sizes and nucleation rates of water

We have calculated the critical cluster sizes and homogeneous nucleation rates of water at temperatures and vapor densities corresponding to experiments by Wolk and Strey [J. Phys. Chem B 105, 11683 (2001)]. The calculations have been done with an expanded version of a Monte Carlo method originally developed by Vehkamaki and Ford [J. Chem. Phys. 112, 4193 (2000)]. Their method calculates the statistical growth and decay probabilities of molecular clusters. We have derived a connection between these probabilities and kinetic condensation and evaporation rates, and introduce a new way for the calculation of the work of formation of clusters. Three different interaction potential models of water have been used in the simulations. These include the unpolarizable SPC/E [J. Phys. Chem. 91, 6269 (1987)] and TIP4P [J. Chem. Phys. 79, 926 (1983)] models and a polarizable model by Guillot and Guissani [J. Chem. Phys. 114, 6720 (2001)]. We show that TIP4P produces critical cluster sizes and a temperature and vapor density dependence for the nucleation rate that agree well with the experimental data, although the magnitude of nucleation rate is constantly overestimated by a factor of 2 x 10(4). Guissani and Guillot's model is somewhat less successful, but both the TIP4P and Guillot and Guissani models are able to reproduce a much better experimental temperature dependency of the nucleation rate than the classical nucleation theory. Using SPC/E results in dramatically too small critical clusters and high nucleation rates. The water models give different average binding energies for clusters. We show that stronger binding between cluster molecules suppresses the decay probability of a cluster, while the growth probability is not affected. This explains the differences in results from different water models.