Simulating equilibrium surface forces in polymer solutions using a canonical grid method

A new simulation method for nonuniform polymer solutions between planar surfaces at full chemical equilibrium is described. The technique uses a grid of points in a two-dimensional thermodynamic space, labeled by surface area and surface separations. Free energy differences between these points are determined via Bennett's optimized rates method in the canonical ensemble. Subsequently, loci of constant chemical potential are determined within the grid via simple numerical interpolation. In this way, a series of free energy versus separation curves are determined for a number of different chemical potentials. The method is applied to the case of hard sphere polymers between attractive surfaces, and its veracity is confirmed via comparisons with established alternative simulation techniques, namely, the grand canonical ensemble and isotension ensemble methods. The former method is shown to fail when the degree of polymerization is too large. An interesting interplay between repulsive steric interactions and attractive bridging forces occurs as the surface attraction and bulk monomer density are varied. This behavior is further explored using polymer density functional theory, which is shown to be in good agreement with the simulations. Our results are also discussed in light of recent self-consistent field calculations which correct the original deGennes results for infinitely long polymers. In particular, we look at the role of chain ends by investigating the behavior of ring polymers.     

Force-momentum-based self-guided Langevin dynamics: a rapid sampling method that approaches the canonical ensemble

The self-guided Langevin dynamics (SGLD) is a method to accelerate conformational searching. This method is unique in the way that it selectively enhances and suppresses molecular motions based on their frequency to accelerate conformational searching without modifying energy surfaces or raising temperatures. It has been applied to studies of many long time scale events, such as protein folding. Recent progress in the understanding of the conformational distribution in SGLD simulations makes SGLD also an accurate method for quantitative studies. The SGLD partition function provides a way to convert the SGLD conformational distribution to the canonical ensemble distribution and to calculate ensemble average properties through reweighting. Based on the SGLD partition function, this work presents a force-momentum-based self-guided Langevin dynamics (SGLDfp) simulation method to directly sample the canonical ensemble. This method includes interaction forces in its guiding force to compensate the perturbation caused by the momentum-based guiding force so that it can approximately sample the canonical ensemble. Using several example systems, we demonstrate that SGLDfp simulations can approximately maintain the canonical ensemble distribution and significantly accelerate conformational searching. With optimal parameters, SGLDfp and SGLD simulations can cross energy barriers of more than 15 kT and 20 kT, respectively, at similar rates for LD simulations to cross energy barriers of 10 kT. The SGLDfp method is size extensive and works well for large systems. For studies where preserving accessible conformational space is critical, such as free energy calculations and protein folding studies, SGLDfp is an efficient approach to search and sample the conformational space.     

Application of the constrained fluid lambda-integration path to the calculation of high temperature Au(110) surface free energies

Recently a method termed constrained fluid lambda-integration was proposed for calculating the free energy difference between bulk solid and liquid reference states via the construction of a reversible thermodynamic integration path; coupling the two states in question. The present work shows how the application of the constrained fluid lambda-integration concept to solid/liquid slab simulation cells makes possible a generally applicable computer simulation methodology for calculating the free energy of any surface and/or surface defect structure, including surfaces requiring variations in surface atom or density number, such as the (1 x 5) Au(100) or (1 x 2) missing row Au(110) reconstructed surfaces or excess adatom/vacancy/step populated surfaces. We evaluate the methodology by calculating the free energy of various disordered high temperature Au(110) embedded atom method surfaces constrained to differing excess surface atom numbers [including those corresponding to the (1 x 2) missing row reconstructed surface] and obtained the interesting result that at 1000 K (as distinct from lower temperatures) the free energy difference between these surfaces is reduced to zero; a result which is consistent with an expected order-disorder phase transition for the Au(110) surface at such high temperatures.     

On interfacial tension calculation from the test-area methodology in the grand canonical ensemble

We propose the extension of the test-area methodology, originally proposed to evaluate the surface tension of planar fluid-fluid interfaces along a computer simulation in the canonical ensemble, to deal with the solid-fluid interfacial tension of systems adsorbed on slitlike pores using the grand canonical ensemble. In order to check the adequacy of the proposed extension, we apply the method for determining the density profiles and interfacial tension of spherical molecules adsorbed in slitlike pore with different pore sizes and solid-fluid dispersive energy parameters along the same simulation. We also calculate the solid-fluid interfacial tension using the original test-area method in the canonical ensemble. Agreement between the results obtained from both methods indicate that both methods are fully equivalent. The advantage of the new methodology is that allows to calculate simultaneously the density profiles and the amount of molecules adsorbed onto a slitlike pore, as well as the solid-fluid interfacial tension. This ensures that the chemical potential at which all properties are evaluated during the simulation is exactly the same since simulations can be performed in the grand canonical ensemble, mimicking the conditions at which the adsorption experiments are most usually carried out in the laboratory.     

Communication: Hybrid ensembles for improved force matching

Force matching is a method for parameterizing empirical potentials in which the empirical parameters are fitted to a reference potential energy surface (PES). Typically, training data are sampled from a canonical ensemble generated with either the empirical potential or the reference PES. In this Communication, we show that sampling from either ensemble risks excluding critical regions of configuration space, leading to fitted potentials that deviate significantly from the reference PES. We present a hybrid ensemble which combines the Boltzmann probabilities of both potential surfaces into the fitting procedure, and we demonstrate that this technique improves the quality and stability of empirical potentials.     

Effect of surface grafted polymers on the adsorption of different model proteins

Adsorption of a model protein to a surface with end-grafted polymers was studied by Monte Carlo simulations. In the model the effect on protein adsorption in the presence of end-grafted polymers was evaluated by calculating the change in free energy between an end-grafted surface and a surface without polymers. The change in free energy was calculated using statistical mechanical perturbation theory. Apart from ordinary athermal polymer-polymer and protein-polymer interactions we also study a broad selection of systems by varying the interaction between proteins and polymers and effective polymer-solvent interactions. The interactions between the molecules span an interval from -0.5 to +0.5 kT. Consequently, general features of protein adsorption to end-grafted surfaces is investigated by systematically changing properties like hydrophilicity/hydrophobicity of the polymer, protein and surface as well as grafting density, degree of polymerization and protein size. Increasing grafting density as well as degree of polymerization decreases the adsorption of protein except in systems with attractive polymer-protein interactions, where adsorption increases with increasing chain length and higher grafting density. At a critical polymer-protein interaction neither chain length nor grafting density affects the free energy of adsorption. Hydrophilic polymers were found to prevent adsorption better than hydrophobic polymers. Very small particles with radii comparable to the size of a polymer segment were, however, better excluded from the surface when using hydrophobic than hydrophilic polymers. For systems with attractive polymer-protein interaction not only the volume of the protein was shown to be of importance but also the size of the exposed surface.     

Modeling of adsorption and nucleation in infinite cylindrical pores by two-dimensional density functional theory

In this paper, we present an analysis of argon adsorption in cylindrical pores having amorphous silica structure by means of a nonlocal density functional theory (NLDFT). In the modeling, we account for the radial and longitudinal density distributions, which allow us to consider the interface between the liquidlike and vaporlike fluids separated by a hemispherical meniscus in the canonical ensemble. The Helmholtz free energy of the meniscus was determined as a function of pore diameter. The canonical NLDFT simulations show the details of density rearrangement at the vaporlike and liquidlike spinodal points. The limits of stability of the smallest bridge and the smallest bubble were also determined with the canonical NLDFT. The energy of nucleation as a function of the bulk pressure and the pore diameter was determined with the grand canonical NLDFT using an additional external potential field. It was shown that the experimentally observed reversibility of argon adsorption isotherms at its boiling point up to the pore diameter of 4 nm is possible if the potential barrier of 22kT is overcome due to density fluctuations.     

Problems of contact angle and solid surface free energy determination

The current general problems of formulation and determination of surface free energy are discussed. So far several theories and approaches have been proposed, but formulation of surface and interfacial free energy, as regards its components, is still a very debatable issue. However, as long as no method for determination of real surface free energy quantities is known, even relative values charged with many simplified assumptions are useful for better understanding of the wetting processes. In this paper special focus is concentrated on powdered solids for which direct measurement of the contact angles is not possible. For such solids the porous layer imbibition techniques are most frequently applied. Then, using the wicking results the contact angle is calculated from Washburn's equation. However, such a procedure leads to overestimated contact angle values in comparison to those measured directly on smooth surfaces of the same solid, if such surface can be obtained at all. As a consequence, the solid surface free energy components calculated via such overestimated contact angles are significantly lower than those obtained from contact angles measured directly. Methodologies to avoid this problem are also described.     

Simulations of surface forces in polyelectrolyte solutions

We have simulated interactions between charged surfaces in the presence of oppositely charged polyelectrolytes by coupling perturbations in the isotension ensemble to a free energy variance minimization scheme. For polymeric systems, this method completely outperforms configurationally biased versions of grand canonical simulations. Proper diffusive equilibrium between bulk and slit has been established for polyelectrolytes with up to 60 monomers per chain. A consequence of imposing diffusive equilibrium conditions, in contrast to previous more restricted models, is the possibility of surface charge inversion; ion-ion correlation and the cooperativity of monomer adsorption drive the formation of a polyion layer close to the surface, that overcompensates the nominal surface charge. This is observed even at modest surface charge densities, and leads to a build up of a long ranged electrostatic barrier. In addition, the onset of charge inversion requires very low bulk polymer densities. Due to screening effects, this leads to a higher and more long-ranged free energy barrier at low, compared to high, bulk densities. Oscillatory forces, reminiscent of those found in simple hard sphere systems, are resolved in the high concentration regime. As a consequence of a second surface charge inversion, the system "stratifies" to form a stable polyelectrolyte layer in the central part of the slit, stabilized by the adsorbed surface layers.     

Adsorption of proteins onto charged surfaces: A Monte Carlo approach with explicit ions

A computer model has been developed to simulate the adsorption of proteins onto charged surfaces displaying an electric double layer. Coadsorption of ions onto the surface is included by means of explicit ions. Only electrostatic interactions are considered. Monte Carlo simulations in the canonical ensemble of the enzyme cutinase and 15 variants (modeled from the X-ray tertiary structure of the wild-type) were performed. Adsorption free energies for all variants were calculated by the thermodynamic integration method. Distributions of the electric moment and the vector pointing toward the protein active site and parallel to its central β-sheet were determined to elucidate the mean orientation of the protein with respect to the surface as a function of its distance from the surface. It was found that the free energy of adsorption varied linearly with the total charge of the protein, while the electric moment (dipole moment) had a second-order but significant effect. Though an increase of the electric moment generally resulted in a slightly increased affinity of the protein for the surface, close to the surface the mean force acting on the protein clearly varied linearly with the strength of the electric moment, such that a clear correlation between the latter and the protein orientation with respect to the surface could be established. Wild-type cutinase displayed the highest affinity for the charged surface amongst all proteins having the same total charge, even though it did not have the largest electric moment. © 1996 by John Wiley & Sons, Inc.     

Monte Carlo molecular simulation of the hydration of K-montmorillonite at 353 K and 625 bar

Monte Carlo molecular simulations of the hydration of K-saturated Wyoming-type montmorillonite at constant stress in the NPzzT ensemble and at constant chemical potential in the grand canonical muVT ensemble, under basin-like conditions of 353 K and 625 bar, show a strong tendency of the K+ ions to adhere to the siloxane surface, forming predominant inner-sphere complexes with tetrahedral oxygen atoms and adsorbed water molecules. Simulations in the grand canonical ensemble predict that none of the K-montmorillonite hydrates, the one-, two-, and three-layer hydrates, are stable in this environment of high depth, temperature, and pressure. The most nearly stable configuration corresponds to the one-layer hydrate, characterized by a d001 spacing of 12.75 A, the adsorbed water being 60 molecules/layer or 180.83 mg of H2O/g of clay, an internal energy of -22.73 kcal/mol, an interlayer density of 0.365 g/mL, and a pressure tensor, Pzz, of 1999.9 bar. The interlayer structure consists of two close layers of water molecules 0.50 A from the midplane, with broad shoulders on the sides, the protons oriented toward the midplane and the siloxane surfaces, and the K+ ions close to the clay surfaces and on the interlayer midplane.     

Atomistic simulations to compute surface properties of poly(N-vinyl-2-pyrrolidone) (PVP) and blends of PVP/chitosan

Atomistic simulations were performed on poly(N-vinyl-2-pyrrolidone) (PVP) and its blends with chitosan (CS) in different ratios using molecular mechanics (MM) and molecular dynamics (MD) simulations in three-dimensionally periodic and effective two-dimensionally periodic condensed phases. Four independent microstructures were generated to analyze their surface properties. The calculated surface-energy values for PVP compared quite well with the experimental data reported in the literature. The density profile was analyzed, and the structure of the films showed an interior region of the bulk density. Various components of the energetic interactions (torsional, van der Waals, etc.) were examined to gain deeper insight into the nature of regular and anomalous interactions between the bulk and the surface films. Surface energies of PVP/CS blends were computed by MD simulations using the bulk pressure-volume-temperature (PVT) parameters. Bulk properties such as the cohesive energy density (CED) and solubility parameter (delta) were calculated using MM and MD simulations in the NVT ensemble under periodic boundary conditions. The Flory equation of state was used to compute the thermal expansion coefficient as well as PVT parameters. These surface-energy values agreed well with the surface-energy data calculated using the Zisman equation, which were also in accordance with the experimental observations. The results from this study suggest that computer simulations would provide valuable information on polymers and polymer-blend surfaces.     

Numerical approaches to determine the interface tension of curved interfaces from free energy calculations

A recently proposed method to obtain the surface free energy σ(R) of spherical droplets and bubbles of fluids, using a thermodynamic analysis of two-phase coexistence in finite boxes at fixed total density, is reconsidered and extended. Building on a comprehensive review of the basic thermodynamic theory, it is shown that from this analysis one can extract both the equimolar radius R(e) as well as the radius R(s) of the surface of tension. Hence the free energy barrier that needs to be overcome in nucleation events where critical droplets and bubbles are formed can be reliably estimated for the range of radii that is of physical interest. It is found that the conventional theory of nucleation, where the interface tension of planar liquid-vapor interfaces is used to predict nucleation barriers, leads to a significant overestimation, and this failure is particularly large for bubbles. Furthermore, different routes to estimate the effective radius-dependent Tolman length δ(R(s)) from simulations in the canonical ensemble are discussed. Thus we obtain an instructive exemplification of the basic quantities and relations of the thermodynamic theory of metastable droplets/bubbles using simulations. However, the simulation results for δ(R(s)) employing a truncated Lennard-Jones system suffer to some extent from unexplained finite size effects, while no such finite size effects are found in corresponding density functional calculations. The numerical results are compatible with the expectation that δ(R(s) → ∞) is slightly negative and of the order of one tenth of a Lennard-Jones diameter, but much larger systems need to be simulated to allow more precise estimates of δ(R(s) → ∞).     

Solid surface mapping by inverse gas chromatography

Inverse gas chromatography (IGC) at infinite dilution, is a technique for characterising solid surfaces. Current practice is the injection of n-alkane homologous series to obtain the free energy of adsorption of the CH2 group, from which the London component of the solid surface free energy, gamma(d)s, is calculated. A value around 40 mJ/m2 is obtained for poly(ethylene), and 30 mJ/m2 for a clean glass fibre, while the potential surface interactivity of a glass fibre is far greater than that of poly(ethylene). A specific component of the surface, in mJ/m2, should be calculated in order to obtain significant parameters. As applied up to date, when calculating the specific component of the surface energy, the fact that W(sp)a energy values are in a totally different scale than AN or DN values is a major drawback. Consequently, Ka and Kb values obtained are in arbitrary energy units, different from those of the London component measured by injecting the n-alkane series. This paper proposes a method to obtain Ka and Kb values of the surface in the same energetic scale than the London component. The method enables us to correct the traditional London component of a solid, obtaining a new value, where the amount of WaCH2 accounting for Debye interactions with polar sites, is excluded. As a result, an approach to surface mapping is performed in several different substrate materials. We show results obtained on different solid surfaces: poly(ethylene), clean glass fibre, glass beads, chemically modified glass beads and carbon fibre.     

Effective Debye length in closed nanoscopic systems: a competition between two length scales

The Poisson-Boltzmann equation (PBE) is widely employed in fields where the thermal motion of free ions is relevant, in particular in situations involving electrolytes in the vicinity of charged surfaces. The applications of this non-linear differential equation usually concern open systems (in osmotic equilibrium with an electrolyte reservoir, a semi-grand canonical ensemble), while solutions for closed systems (where the number of ions is fixed, a canonical ensemble) are either not appropriately distinguished from the former or are dismissed as a numerical calculation exercise. We consider herein the PBE for a confined, symmetric, univalent electrolyte and quantify how, in addition to the Debye length, its solution also depends on a second length scale, which embodies the contribution of ions by the surface (which may be significant in high surface-to-volume ratio micro- or nanofluidic capillaries). We thus establish that there are four distinct regimes for such systems, corresponding to the limits of the two parameters. We also show how the PBE in this case can be formulated in a familiar way by simply replacing the traditional Debye length by an effective Debye length, the value of which is obtained numerically from conservation conditions. But we also show that a simple expression for the value of the effective Debye length, obtained within a crude approximation, remains accurate even as the system size is reduced to nanoscopic dimensions, and well beyond the validity range typically associated with the solution of the PBE.     

Effect of surface functionalities on gas adsorption in microporous carbons: a grand canonical Monte Carlo study

In an attempt to offer a more realistic picture of adsorption in highly heterogeneous porous systems, such as oxygen functionalized porous carbons, we consider a series of carbon surfaces baring different amounts of oxygen functionalities (hydroxyl and epoxy). These surfaces are used to construct “oxidized” slit pores of varying width and functionality. With the aid of such inhomogeneous structures we study the interaction of Ar (87 K) inside “functionalized” pores and report grand canonical Monte Carlo adsorption simulations results. Based on our simulation data, we discuss the role of chemical heterogeneity on adsorbed/gas phase equilibrium properties such as density, heat of adsorption, and molecular packing within the pores. Comparisons are made with the case of the oxygen–free (completely homogeneous) slit pore models and conclusions on the suitability of Ar based pore size distributions for functionalized porous carbons are drawn.