Good + Bad = Ugly, and other pitfalls to avoid when calculating free energies

Molecular simulation is an experiment performed on a system defined by a molecular model. In recent years much of the effort to advance molecular simulation has been put toward improving the quality and quantity of information that it provides. These efforts have been very successful, and perhaps have improved the utility of molecular simulation even more than advances in raw computing power occurring over the same period. One important application of molecular simulation is the calculation of free-energy differences, which are required for analyses of phase and reaction equilibria, solvation, binding affinity, stability, kinetics, and so on. Some of the most popular approaches to calculating the free energy are highly prone to systematic errors, and simple countermeasures designed to remedy these inaccuracies often do not improve the outcome. We show that the key consideration influencing the accuracy is the overlap of the important regions of phase space for the systems of interest. We are developing measures to quantify this overlap and we examine the connection between them and the performance of the calculations. We use these ideas to formulate simple variants of the basic technique that can be applied to increase the likelihood of obtaining a good result

Evaluation of bridge-function diagrams via mayer-sampling Monte Carlo simulation

We report coefficients of the h-bond expansion of the bridge function of the hard-sphere system up to order rho(4) (where rho is the density in units of the hard-sphere diameter), which in the highest-order term includes 88 cluster diagrams with bonds representing the total correlation function h(r). Calculations are performed using the recently introduced Mayer-sampling method for evaluation of cluster integrals, and an iterative scheme is applied in which the h(r) used in the cluster integrals is determined by solution of the Ornstein-Zernike equation with a closure given by the calculated clusters. Calculations are performed for reduced densities from 0.1 to 0.9 in increments of 0.1. Comparison with molecular simulation data shows that the convergence is very slow for the density expansion of the bridge function calculated this way.open9

Calculation of high-order virial coefficients for the square-well potential

Accurate virial coefficients BN(λ,ε) (where ε is the well depth) for the three-dimensional square-well and square-step potentials are calculated for orders N = 5–9 and well widths λ = 1.1−2.0 using a very fast recursive method. The efficiency of the algorithm is enhanced significantly by exploiting permutation symmetry and by storing integrands for reuse during the calculation. For N = 9 the storage requirements become sufficiently large that a parallel algorithm is developed. The methodology is general and is applicable to other discrete potentials. The computed coefficients are precise even near the critical temperature, and thus open up possibilities for analysis of criticality of the system, which is currently not accessible by any other means

Cluster integrals and virial coefficients for realistic molecular models

We present a concise, general, and efficient procedure for calculating the cluster integrals that relate thermodynamic virial coefficients to molecular interactions. The approach encompasses nonpairwise intermolecular potentials generated from quantum chemistry or other sources; a simple extension permits efficient evaluation of temperature and other derivatives of the virial coefficients. We demonstrate with a polarizable model of water. We argue that cluster-integral methods are a potent yet underutilized instrument for the development and application of first-principles molecular models and methods

Elastic constants and the effect of strain on monovacancy concentration in fcc hard-sphere crystals

We investigate the free energy and the concentration of monovacancies in strained face-centered-cubic (fee) hard-sphere crystals for several densities at and above melting. We use the conventional molecular dynamics method for simulations and employ a bias insertion method to extract properties of a monovacancy. We study two distinct constant-volume strains, considering a simple shear and an orthogonal expansion and contraction. Strains are examined across the linear elastic region and include also some nonlinear elastic deformations. Second-order elastic constants are reported as a function of density. The concentration of monovacancies decreases as density increases for both strained and unstrained crystals. The effect of strain is to cause the monovacancy concentration to increase by up to 72% for the expansion-contraction strain at the largest deformation studied. The effect of the shear strain is considerably less, and produces an increase in monovacancy concentration of at most 9% for the conditions studied here.open5

Efficiency of spinal anesthesia versus general anesthesia for lumbar spinal surgery: a retrospective analysis of 544 patients.

BACKGROUND: Previous studies have shown varying results in selected outcomes when directly comparing spinal anesthesia to general in lumbar surgery. Some studies have shown reduced surgical time, postoperative pain, time in the postanesthesia care unit (PACU), incidence of urinary retention, postoperative nausea, and more favorable cost-effectiveness with spinal anesthesia. Despite these results, the current literature has also shown contradictory results in between-group comparisons. MATERIALS AND METHODS: A retrospective analysis was performed by querying the electronic medical record database for surgeries performed by a single surgeon between 2007 and 2011 using procedural codes 63030 for diskectomy and 63047 for laminectomy: 544 lumbar laminectomy and diskectomy surgeries were identified, with 183 undergoing general anesthesia and 361 undergoing spinal anesthesia (SA). Linear and multivariate regression analyses were performed to identify differences in blood loss, operative time, time from entering the operating room (OR) until incision, time from bandage placement to exiting the OR, total anesthesia time, PACU time, and total hospital stay. Secondary outcomes of interest included incidence of postoperative spinal hematoma and death, incidence of paraparesis, plegia, post-dural puncture headache, and paresthesia, among the SA patients. RESULTS: SA was associated with significantly lower operative time, blood loss, total anesthesia time, time from entering the OR until incision, time from bandage placement until exiting the OR, and total duration of hospital stay, but a longer stay in the PACU. The SA group experienced one spinal hematoma, which was evacuated without any long-term neurological deficits, and neither group experienced a death. The SA group had no episodes of paraparesis or plegia, post-dural puncture headaches, or episodes of persistent postoperative paresthesia or weakness. CONCLUSION: SA is effective for use in patients undergoing elective lumbar laminectomy and/or diskectomy spinal surgery, and was shown to be the more expedient anesthetic choice in the perioperative setting

Calculation of surface tension via area sampling

We examine the performance of several molecular simulation techniques aimed at evaluation of the surface tension through its thermodynamic definition. For all methods explored, the surface tension is calculated by approximating the change in Helmholtz free energy associated with a change in interfacial area through simulation of a liquid slab at constant particle number, volume, and temperature. The methods explored fall within three general classes: free-energy perturbation, the Bennett acceptance-ratio scheme, and the expanded ensemble technique. Calculations are performed for both the truncated Lennard-Jones and square-well fluids at select temperatures spaced along their respective liquid-vapor saturation lines. Overall, we find that Bennett and expanded ensemble approaches provide the best combination of accuracy and precision. All of the methods, when applied using sufficiently small area perturbation, generate equivalent results for the Lennard-Jones fluid. However, single-stage free-energy-perturbation methods and the closely related test-area technique recently introduced by Gloor et al. [J. Chem. Phys. 123, 134703 (2005)], generate surface tension values for the square-well fluid that are not consistent with those obtained from the more robust expanded ensemble and Bennett approaches, regardless of the size of the area perturbation. Single-stage perturbation methods fail also for the Lennard-Jones system when applied using large area perturbations. Here an analysis of phase-space overlap produces a quantitative explanation of the observed inaccuracy, and shows that the satisfactory results obtained in these cases from the test-area method arise from a cancellation of errors that cannot be expected in general.Comment: 28 pages, 7 figures, to appear in J. Chem. Phys. (07 Nov 2007 issue

Communication: Analytic continuation of the virial series through the critical point using parametric approximants

The mathematical structure imposed by the thermodynamic critical point motivates an approximant that synthesizes two theoretically sound equations of state: the parametric and the virial. The former is constructed to describe the critical region, incorporating all scaling laws; the latter is an expansion about zero density, developed from molecular considerations. The approximant is shown to yield an equation of state capable of accurately describing properties over a large portion of the thermodynamic parameter space, far greater than that covered by each treatment alone

Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line

Thermodynamic integration along a path that coincides with the saturation line is proposed as an efficient means for evaluation of phase equilibria by molecular simulation. The technique allows coexistence to be determined by just one simulation, without ever attempting or performing particle insertions. Prior knowledge of one coexistence point is required to start the procedure. Integration then advances from this state according to the Clapeyron formula-a first-order ordinary differential equation that prescribes how the pressure must change with temperature to maintain coexistence. The method is unusual in the context of thermodynamic integration in that the path is not known at the outset of the process; results from each simulation determine the course that the integration subsequently takes. Predictor-corrector methods among standard numerical techniques are shown to be particularly well suited for this type of integration. A typical integration step along the saturation line proceeds as follows: An increment in the temperature is chosen, and the saturation pressure at the new temperature is "predicted" from previous data (the initial coexistence datum and/or previous simulations). Simultaneous but independent NPT simulations of the coexisting phases are initiated at the said conditions. Averages taken throughout the simulations are repeatedly used to "correct" the estimate of the pressure to convergence. Thus strictly the pressure is not fixed during the simulation. Vapor-liquid coexistence of the van der Waals model is first used to study the numerical integration method without the complications of molecular simulation. In a second application the phase envelope of the Lennard-Jones model fluid is computed, and many variations of the technique are examined. Overall, the results are remarkably consistent and in agreement with previous simulation studies. Difficulty is encountered upon approach of the critical point, but, by artificially coupling the simulation volumes, the method remains effective in this regime so long as a suitably small integration step is employed. Many extensions and improvements of the technique are discussed