Efficient Test Molecule Sampling in Molecular Simulation

Abstract

Excluded volume map sampling (EVMS) is a particularly efficient means of performing test molecule sampling to estimate or impose chemical potential in molecular simulations. This paper discusses the motivation and applications of excluded volume map sampling, presents computer code demonstrating its implementation, and gives an example of its application.     

Computing free energies using nested sampling-based approaches

Abstract

Nested sampling (NS) has emerged as a powerful statistical mechanical sampling technique to compute the partition function of atomic and molecular systems. From the partition function all thermodynamic quantities can be calculated in absolute terms, including absolute free energies and entropies. In this article, we provide a brief overview of NS within a Bayesian context, as well as overviews of how NS is used to compute the partition functions and thermodynamic quantities in the canonical and isothermal-isobaric ensembles. Then we introduce a new scheme, Coupling Parameter Path Nested Sampling, to estimate the free energy difference between two systems with different potential energy functions. The method uses a NS simulation to traverse the same path through phase space as would be covered in traditional coupling parameter-based methods such as thermodynamic integration and perturbation approaches. We demonstrate the new method with two case studies and confirm its accuracy by comparison to conventional methods, including Widom test particle insertion and thermodynamic integration. The proposed method provides a powerful alternative to traditional coupling parameter-based free energy simulation methods.     

A New Potential Model for Carbon Dioxide from AB Initio Calculations

Abstract

Ab initio quantum chemical calculations have been carried out for carbon dioxide dimer and the results have been used to establish potential functions usable in molecular simulations. Since the intermolecular interaction in carbon dioxide is fairly weak, careful treatment is required: this study uses 6–31G* basis set and takes electron correlations by the 2nd order Møller-Plesset theory into account. The potential energy surface is elucidated using the four representative relative configurations of the dimer. A new potential function model has been proposed on the basis of these ab initio data. In the super-critical region, this model is used to calculate the PVT relation of carbon dioxide fluid by the Monte Carlo simulations and confirmed to reproduce reasonably well the experimental isotherms.     

Chemical Potential,Partial Enthalpy and Partial Volume of Mixtures by NPT Molecular Dynamics

Abstract

Equilibrium NPT molecular dynamics computer simulations have been used to determine the chemical potential, partial enthalpy and partial volume of model Ar-Kr mixtures using newly devised non-intrusive particle insertion and particle swap techniques [P. Sindzingre et al. Chemical Physics, 129 (1989) 213]. In this report we examine, for the first time, in some detail the relative convergence statistics of the particle swap and particle insertion methods for these properties for binary Lennard-Jones (LJ) mixtures. Both species are represented by single-site Lennard-Jones pair potentials with Lorentz-Berthelot rules for the cross-species interactions. We show that, over the whole phase diagram and especially in the vicinity of the fluid-solid coexistence line, the particle swap method gives significantly better statistics than the particle insertion method for the difference in chemical potential of the two species, partial enthalpy and partial volume of each species. Also, we find that, using the particle swap method, the difference in the chemical potential converges more rapidly than the differences in the partial enthalpy and volume.     

A New Monte Carlo Method for Direct Calculation of the Critical Size and the Formation Work of a Microdrop

New Monte Carlo procedures in open ensembles are proposed. Non-stationary Markov chain procedure in the μl;pT - ensemble provides a direct estimation for the critical size of a condensation nucleus at given p and T. A stationary procedure in the μlpT ensemble with two allowed particle numbers n and n + 1 provides the direct way to calculate the chemical potential and Gibbs free energy of a cluster; in the grand canonical (μlVT) ensemble the same approach gives μl and the Helmholtz free energy. The same procedures are readily applicable to periodic systems representing bulk phases.     

Predicting Liquid–Vapour Equilibria for Water Using an ab-initio Potential from Histogram Reweighting Monte Carlo Simulations

Abstract

The coexisting densities for an ab-initio model for water have been calculated using grand canonical Monte Carlo simulations with the histogram reweighting technique. Although good agreement with experimental data is found for the radial distribution function at room temperature, the predicted critical density and temperature are well below both the experimental value as well as predictions from semi-empirical potentials. Improvement in the repulsive part of the ab-initio potential is suggested as a way to obtain better agreement with experiment.     

On the Application of Widom's Test Particle Method to Homogeneous and Inhomogeneous Fluids

Abstract

Chemical potentials of a homogeneous and an inhomogeneous Lennard-Jones fluid have been determined by molecular dynamics simulations on the vector computer CYBER 205 by applying essentially the fictitious test particle method of Widom. For the homogeneous fluid we find, contrary to the previous result of Guillot and Guissani, that the simulated chemical potential is independent of the particle number. The crucial point, however, is a sufficiently large cut-off radius in the evaluation of the Boltzmann factor. Comparing with our WCA-type perturbation theory, we get agreement in the chemical potentials within 0.1 kT up to the density n[sgrave]³ = 0.80 and a difference of 0.2 kT at n[sgrave]³ = 0.85. For the inhomogeneous case we consider a fluid in a cylindrical pore and integrate Widom's equation over a certain probe volume as suggested earlier by us. Chemical potentials are then calculated independently in five different probe volumes, which are cylindrical shells. The results agree well from the second to the fourth shell. Inaccuracies in the innermost cylinder can be easily explained by bad statistics. In the shell close to the wall the extremely high local density is responsible for the inaccuracies. Extending the probe volume over all cylindrical shells besides the one closest to the wall is thought to yield rather reliable results for the chemical potential. As a by-product of the simulations we also obtained diffusion coefficients, which are given in an appendix.     

Monte Carlo Simulations of the Chemical Potential and Free Energy for Trimer and Hexamer Rings

Abstract

The chemical potential of a trimer and hexamer model ring system was determined by computer simulation over a range of temperatures and densities. Such ring molecules are important as model aromatic and naphthenic hydrocarbons. Thermodynamic integration of the pressure along a reversible path, Widom's ghost particle insertion method and Kirkwood's charging parameter method were used over a molecular density range of 0.05 to 0.30. Data were obtained by Monte Carlo simulation of a 96 molecule system that was modelled with a Lennard-Jones 6-12 truncated potential. The original insertion method, which does not take into account the orientation of the molecule when it is inserted, gives results for the chemical potential which deviate from that obtained using the thermodynamic pressure integration. At high density or temperature the deviation is significant. We have modified the Widom insertion technique to account for this short range orientation and find good agreement between this technique and the thermodynamic integration method for the chemical potential. We also calculated the free energy difference between our model ring molecules and ring molecules made up of hard spheres.     

The Lennard-Jones Fluid Revisited: Computer Simulation Results

Abstract

Pseudoexperimental data of high accuracy on the pressure and the internal energy of the Lennard-Jones fluid have been generated both by the Monte Carlo and molecular dynamics methods for five subcritical and three supercritical isotherms. Values of the chemical potential of the Lennard-Jones fluid computed by a new version of the gradual insertion particle method for two isotherms up to very high densities are also reported and discussed, and compared with existing data.     

The Chemical Potential of Dense Liquids by the Coupled Test Particle Method

Abstract

A new method for chemical potential estimation is proposed which is based on the coupled particle approach. The coupled particle method defines an attractive solution to the weighting function problem in umbrella sampling, bridging the gap between f and g distributions at high density. A way of eliminating the origin singularity is suggested, which is similar in spirit to the restricted umbrella sampling of Shing and Gubbins, but which is based on geometric rather than energetic criteria.

The method is illustrated on the Lennard-Jones system up to a reduced density ρ? = 1.1 along the isotherm T? = 1.2 and results are compared with the test particle insertion method and empirical equations of state. The new method is particularly useful at high liquid densities where it is superior to the other methods relying on the degree of overlap of f and g distributions. It gives reliable estimates of the chemical potential in the whole range of liquid densities.     

Grand Molecular Dynamics: A Method for Open Systems

We present a new molecular dynamics method for studying the dynamics of open systems. The method couples a classical system to a chemical potential reservior. In the formulation, following the extended system dynamics approach, we introduce a variable, v to represent the coupling to the chemical potential reservoir. The new variable governs the dynamics of the variation of number of particles in the system. The number of particles is determined by taking the integer part of v. The fractional part of the new variable is used to scale the potential energy and the kinetic energy of an additional particle: i.e., we introduce a fractional particle. We give the ansatz Lagrangians and equations of motion for both the isothermal and the adiabatic forms of grand molecular dynamics. The averages calculated over the trajectories generated by these equations of motion represent the classical grand canonical ensemble (μVT) and the constant chemical potential adiabatic ensemble (μVL) averages, respectively. The microcanonical phase space densities of the adiabatic and isothermal forms the molecular dynamics method are shown to be equivalent to adiabatic constant chemical potential ensemble, and grand canonical ensemble partition functions. We also discuss the extension to multi-component systems, molecular fluids, ionic solutions and the problems and solutions associated with the implementation of the method. The statistical expressions for thermodynamic functions such as specific heat; adiabatic bulk modulus, Grüneissen parameter and number fluctuations are derived. These expressions are used to analyse trajectories of constant chemical potential systems.     

Monte Carlo Study of the Buckingham Exponential-six Fluid

Abstract

In this paper we report the results of extensive Monte Carlo simulations of a pure fluid of Buckingham modified exponential-six molecules. Data are presented for the configurational energy and pressure covering a wide range of temperatures and densities. These data are interpreted using the generalized van der Waals partition function with a novel separation into free volume and mean potential terms. We find, surprisingly, that the Buckingham fluid is described by a simple van der Waals-like equation of state provided that the b parameter is temperature dependent and chosen in a theoretically correct manner.     

Direct Evaluation of Vapour-Liquid Equilibria of Mixtures by Molecular Dynamics Using Gibbs-Duhem Integration

Abstract

We present an extension of the Gibbs-Duhem integration method that permits direct evaluation of vapour-liquid equilibria of mixtures by molecular dynamics. The Gibbs-Duhem integration combines the best elements of the Gibbs ensemble Monte Carlo technique and thermodynamic integration. Given conditions of coexistence of pure substances, simultaneous but independent molecular dynamics simulations of each phase at constant number of particles, constant pressure, constant temperature and constant fugacity fraction of species 2 are carried out in succession along coexistence lines. In each simulation, the coexistence pressure is adjusted to satisfy the Clapeyron-type equation. The Clapeyron-type equation is a first-order nonlinear differential equation that prescribes how the pressure must change with the fugacity fraction of species 2 to maintain coexistence at constant temperature. The Clapeyron-type equation is solved by the predictor-corrector method. Running averages of mole fraction and compressibility factor for the two phases are used to evaluate the right-hand side of the Clapeyron-type equation. The Gibbs-Duhem integration method is applied to three prototypes of binary mixtures of the two-centre Lennard-Jones fluid having various elongations. The starting points on the coexistence curve were taken from published data.     

A New Version of the Insertion Particle Method for Determining the Chemical Potential by Monte Carlo Simulation

A new version of the test particle method for determining the chemical potential by Monte Carlo simulations is proposed. The method, applicable to any fluid at any density, combines the Widom's test particle insertion method with the ideas of the scaled particle theory, gradual insertion method and multistage sampling. Its applicability is exemplified by evaluating the chemical potential of the hard sphere fluid at a very high density in semi-grand-canonical and grand-canonical ensembles. A theory estimating the efficiency (i.e. statistical errors) of the method is proposed and the results are compared with the Widom's and gradual insertion methods, and the analytic results.     

Monte-Carlo Calculations for Methane and Argon over a Wide Range of Density and Temperature,Including the Two-Phase Vapor-Liquid Region

Abstract

The Gibbs-ensemble simulation technique provides a powerful method to calculate vapor-liquid phase behavior [1]. To evaluate the configurational energy of a system of molecules, commonly used experessions describe the interaction between two molecules. Contributions from higher-body forces are usually implicitly taken into account by adjusting two-body potential parameters to give agreement with experimental data. Explicit expressions for higher-body potentials are not commonly used in simulations [8]. The work by Smit et al. [9] gives the appropriate expressions to evaluate the pressure as well as the chemical potential from a density-dependent two-body potential in an NVT ensemble.

In the present work, contributions to the potential from two-body interactions are separated from those due to higher-body interactions; to take higher-body forces into account, a mean-field term, proportional to (density)^0.9, is added to the two-body potential. NPT-simulations over a wide range of temperature and density, as well as Gibbs-ensemble simulations, are used to evaluate phase behavior of argon and of methane. The results indicate that a simple mean-field correction to the “true” two-body Kihara potential provides good agreement between experiment and simulation.     

Computational approaches to study protein unfolding: Hen egg white lysozyme as a case study

Four methods are compared to drive the unfolding of a protein: (1) high temperature (T-run), (2) high pressure (P-run), (3) by imposing a gradual increase in the mean radius of the protein using a penalty function added to the physical interaction function (F-run, radial force driven unfolding), and (4) by weak coupling of the difference between the temperature of the radially outward moving atoms and the radially inward moving atoms to an external temperature bath (K-run, kinetic energy driven unfolding). The characteristic features of the four unfolding pathways are analyzed in order to detect distortions due to the size or the type of the applied perturbation, as well as the features that are common to all of them. Hen egg white lysozyme is used as a test system. The simulations are analyzed and compared to experimental data like ¹H-NMR amide proton exchange-folding competition, heat capacity, and compressibility measurements. © 1995 Wiley-Liss, Inc.     

Hysteresis and Statistical Errors in Free Energy Perturbation L to D Amino Acid Conversion

Abstract

A theory based on a Langevin equation along the reaction coordinate is developed to explain and calculate systematic and statistical errors in free energy perturbation simulations. The errors are calculated exactly when both the perturbation potential and the mean potential from the surrounding degrees of freedom are harmonic in the reaction coordinate. The effect of the mean potential is small as long as the force constant is small compared to the force constant of the perturbation potential. This indicates that the results obtained with zero mean force may still be valid as long as the second derivate of the mean potential is small compared to that of the perturbation potential. The theory is applied to conversion between L and D amino acids by changing the position of the minimum of the harmonic improper dihedral potential between ±35.264 degrees. For phenylalanine bound in the active site of a protein (thermolysin) we find from 20 psec. simulations statistical errors and hysteresis that both are about 2.5 kJ/mol in agreement with what is obtained from the theoretical predictions. The statistical errors are proportional to the square root of the coupling to the heat bath and inversely proportional to the square root of integration time while the (positive) hysteresis due to that the reaction coordinate lags behind is linear in the same quantities. This shows that the systematic errors will dominate in short simulations while the statistical ones will dominate for long simulations. The treatment is based on that the systematic influence of the surroundings can be represented by a mean force upon the reaction coordinate. If the relaxation processes of the environment are slow this may not be true. Then additional errors have to be considered.     

Mass and Energy Transport Through Slit Pores: Application to Planar Poiseuille Flow

Abstract

We develop a simple, efficient and general statistical mechanical technique for calculating the pressure tensor and the heat flux vector in atomic fluids. The method is applied to the case of planar Poiseuille flow through a narrow slit pore and the results indicate that our technique is accurate and relatively efficient. A second method to calculate shear stress is derived from the momentum continuity equation. This mesoscopic method again is seen to be accurate with good computational efficiency.

We also find that the commonly used approximation to the Irving-Kirkwood expression for the heat flux and the pressure tensor (where the Irving-Kirkwood O_ij operator is set equal to unity-the so-called IK1 approximation), leads to incorrect results for highly inhomogeneous fluids. In such cases the pressure tensor and heat flux vector display spurious oscillations.

We calculate the spatially dependent viscosity across a narrow pore and find that it exhibits real but weak oscillations, a consequence of oscillations in the number density. Finally we point out that if the heat flux vector is coupled to the gradient of the square of the strain rate tensor such an effect will only affect the shape of the temperature profile. For planar Poiseuille flow, the temperature profile should deviate from the classical quartic form and include an additional quadratic component. The actual magnitude and shape of the heat flux vector remain exactly as they would if such a coupling did not exist.     

On the Implementation of Friedman Boundary Conditions in Liquid Water Simulations

Abstract

We have applied the image approximation to the reaction field as suggested by H.L. Friedman [Mol. Phys., 29, 1533 (1975)] by investigating appropriate cavity sizes and system parameters for use in molecular simulations. The energy of and the structure around a central simple point charge (SPC) water molecule in a dielectric cavity was found to be in good agreement with the properties of a liquid sample. To confine the water molecules within the cavity, we introduced a short-range repulsion between a real charge and its image as the Lennard-Jones repulsive potential between oxygen atoms of the SPC potential. For a system of 65 water molecules a cavity radius of 10.45 Å is appropriate; this radius is altered to 12.00 Å for a cavity surrounding 113 molecules. The effect of the boundary is restricted to the outer-most water layer which is in contact with the dielectric continuum.     

Combined Diffusive and Viscous Transport of Methane in a Carbon Slit Pore

Abstract

The transport of mass through porous materials can occur by essentially two different mechanisms: (1) diffusion and (2) viscous flow. The former occurs when there is a gradient in chemical potential of the pore fluid, while the latter occurs in the presence of a pressure gradient. In general, fluid transport occurs by both of these mechanisms and their respective contributions to the total intra-pore flux are approximately additive. Experimentally, there is no unambiguous way of determining the individual contributions to the total flux of these two modes of transport. Fortunately, molecular simulations does provide a solution.

We present a novel simulation method in which the separate contributions to the total flux are determined. The method involves the use of two non-equilibrium molecular dynamics techniques: dual control volume grand canonical molecular dynamics (DCV GCMD) and an algorithm for simulating planar Poiseuille flow. We apply this technique to study the combined (viscous and diffusive) transport of methane through single slit-shaped graphite pores of width 2.5, 5.0 and 10.0 methane diameters. We find that the viscous contribution to the total intrapore flux through each of these pores is 10%, 15% and 34%, respectively.