Andrew J. Schultz

Sixth, seventh and eighth virial coefficients of the Lennard-Jones model

We report values of the virial coefficients B n of the Lennard-Jones (LJ) model, as computed by the Mayer Sampling Monte Carlo method. For n = 4 and 5, values are reported for 103 temperatures T = 0.62 to 40.0 (in LJ units); for n = 6, 31 values are reported for T = 0.625 to 20.0; for n = 7, 15 values are reported from T = 0.625 to 10; and for n = 8, four values are reported from T = 0.75 to 10. Data are used to estimate the location of the LJ critical point, and the critical temperature estimated this way is given to within 0.8% of the established value, while the critical density is too low by 10%. Data derived from the virial equation of state (VEOS) are compared to pressures and internal energies calculated by Monte Carlo simulation. Simulations of systems ranging from 125 to 30,000 particles are extrapolated to infinite system size, and it is shown that the VEOS–when applied at densities where the series has reached convergence–provides results closer to the infinite-system values than obtained by any of the finite-system simulations. For n = 6, convergence of VEOS (within a 1% tolerance) is obtained for densities up to the spinodal for subcritical temperatures and up to ρ = 0.4 (in LJ units) in the vicinity of the critical temperature; the range of applicability of VEOS increases with temperature, reaching for example densities of 0.65 for T = 5.0 and 0.8 for T = 8.0 when truncated at n = 6.

Computer simulation of copolymer phase behavior

Discontinuous molecular dynamics simulation is used to study the phase behavior of diblock copolymers modeled as chains of tangent hard spheres with square shoulder repulsions between unlike species as a function of chain length, volume fraction and interaction strength (χ). The location of the order–disorder transition for a symmetric copolymer is close to the predictions of Fredrickson and Helfand. Our simulation results for packing fractions of 0.35, 0.40, and 0.45 and chain lengths 10 and 20 are summarized in phase diagrams which display disordered, lamellae, perforated lamellae (PL), cylindrical, and BCC spherical (S) phases in the χN versus f plane. These phase diagrams are consistent with phase diagrams from other simulation studies. Contrary to theoretical predictions we observe the PL phase near regions of predicted gyroid stability, and the S phase only in the systems with high packing fraction and long chain length. These discrepancies may be due to the short chain lengths considered, as they a...

Fourth and fifth virial coefficients of polarizable water.

We report values of the viral coefficients B4 and B5 for the Gaussian charge polarizable model (GCPM) of water using the overlap-sampling implementation of Mayer sampling molecular simulation. These results supplement values for the lower-order coefficients B2 and B3 reported previously, and in the present work, we provide more precise values of these coefficients as well. The precision of all viral coefficients is such that the standard error in the calculated pressure is significantly less than 1% for most temperatures, with the exception of temperatures near the critical, where the error approaches 100% at the critical density, and supercritical, where the uncertainty in B5 introduces an error of about 5% in the pressure at the critical density. We examine these coefficients in the context of the equation of state and molecular clustering. Comparisons are made to established molecular simulation data, quantum chemical calculations, and experimental data for real water. Over both sub- and supercritical temperatures, the viral series to B5 is accurate for densities only up to about half the critical density. In this regime, deviation is observed from experimental data for real water, and it is suggested that further development of the model might do well to further improve the relatively good agreement it has with the experimental second viral coefficient of water. The viral coefficients are used to characterize molecular clusters (dimers through pentamers) in GCPM water under supercritical and saturated vapor conditions between 210 and 673 K.

Virial coefficients of model alkanes

We report the results from Mayer-sampling Monte Carlo calculations of the virial coefficients of the united-atom TraPPE-UA model of normal alkanes. For alkane chain lengths from n=2 to 20 (where n is the number of carbon atoms), results are given for the virial coefficients B(2), B(3), and B(4); results for B(5) are given for chains up to length n=12; and results for B(6) are given for chains of length n=2, 3, and 4. In all cases, values are given for temperatures ranging from 200 K to 2000 K in 20-50 K increments. The values are used to calculate the equation of state for butane and the pressure-density behavior is compared to experimental data at 350 and 550 K. Critical points are calculated for all systems and compared to simulation data previously taken for the same molecular model, and to experiment. The comparison with temperature is very good (within 1.5% for all chain lengths up to n=12), while the critical density is underestimated by about 5%-15% and the critical pressure is given within about 10%. The convergence behavior of the virial equation of state as applied across the n-alkane series is well characterized by corresponding states, meaning that the accuracy at a given density relative to the critical density does not deteriorate with increasing chain length.

Mayer Sampling Monte Carlo calculation of virial coefficients on graphics processors

Virial coefficients for various molecular models are calculated up to B 8 using the Mayer Sampling Monte Carlo method and implemented on a graphics processing unit (GPU). The execution time and performance of these calculations is compared with equivalent computations done on a CPU. The speedup between virial coefficient computations on a CPU (w/optimized C code) and a GPU (w/CUDA) is roughly two orders of magnitude. We report values of B 6, B 7, and B 8 of the Lennard-Jones (LJ) model, as computed on the GPU, for temperatures T = 0.6 to 40 (in LJ units).

Lattice Strain Due to an Atomic Vacancy

Volumetric strain can be divided into two parts: strain due to bond distance change and strain due to vacancy sources and sinks. In this paper, efforts are focused on studying the atomic lattice strain due to a vacancy in an FCC metal lattice with molecular dynamics simulation (MDS). The result has been compared with that from a continuum mechanics method. It is shown that using a continuum mechanics approach yields constitutive results similar to the ones obtained based purely on molecular dynamics considerations.

Path-integral Mayer-sampling calculations of the quantum Boltzmann contribution to virial coefficients of helium-4

We present Mayer-sampling Monte Carlo calculations of the quantum Boltzmann contribution to the virial coefficients B(n), as defined by path integrals, for n = 2 to 4 and for temperatures from 2.6 K to 1000 K, using state-of-the-art ab initio potentials for interactions within pairs and triplets of helium-4 atoms. Effects of exchange are not included. The vapor-liquid critical temperature of the resulting fourth-order virial equation of state is 5.033(16) K, a value only 3% less than the critical temperature of helium-4: 5.19 K. We describe an approach for parsing the Boltzmann contribution into components that reduce the number of Mayer-sampling Monte Carlo steps required for components with large per-step time requirements. We estimate that in this manner the calculation of the Boltzmann contribution to B(3) at 2.6 K is completed at least 100 times faster than the previously reported approach.

A comparative study of methods to compute the free energy of an ordered assembly by molecular simulation

We present a comparative study of methods to compute the absolute free energy of a crystalline assembly of hard particles by molecular simulation. We consider all combinations of three choices defining the methodology: (1) the reference system: Einstein crystal (EC), interacting harmonic (IH), or r(-12) soft spheres (SS); (2) the integration path: Frenkel-Ladd (FL) or penetrable ramp (PR); and (3) the free-energy method: overlap-sampling free-energy perturbation (OS) or thermodynamic integration (TI). We apply the methods to FCC hard spheres at the melting state. The study shows that, in the best cases, OS and TI are roughly equivalent in efficiency, with a slight advantage to TI. We also examine the multistate Bennett acceptance ratio method, and find that it offers no advantage for this particular application. The PR path shows advantage in general over FL, providing results of the same precision with 2-9 times less computation, depending on the choice of a common reference. The best combination for the FL path is TI+EC, which is how the FL method is usually implemented. For the PR path, the SS system (with either TI or OS) proves to be most effective; it gives equivalent precision to TI+FL+EC with about 6 times less computation (or 12 times less, if discounting the computational effort required to establish the SS reference free energy). Both the SS and IH references show great advantage in capturing finite-size effects, providing a variation in free-energy difference with system size that is about 10 times less than EC. This result further confirms previous work for soft-particle crystals, and suggests that free-energy calculations for a structured assembly be performed using a hybrid method, in which the finite-system free-energy difference is added to the extrapolated (1/N→0) absolute free energy of the reference system, to obtain a result that is nearly independent of system size.

Effect of Cu and Ag solute segregation on βSn grain boundary diffusivity

We investigate the effect of various amounts of Ag and Cu solute atoms on the self-diffusivity of Sn in the (101) symmetric tilt βSn grain boundary. Using molecular dynamics simulations over a temperature range of 300 K to 450 K, we show that both Ag and Cu decrease the grain boundary self-diffusivity of Sn as the amount of solute in the interface increases. Additionally, the presence of Ag at the grain boundary interface causes a greater reduction in the self-diffusivity of Sn when compared to Cu. We also analyze the solute effect on the diffusive width of the interface and find that low concentrations of both Ag and Cu shrink the width relative to the pure βSn interface. However, adding Cu in excess density greater than 5 × 10−3 atoms/A2 causes the interface to expand to values almost twice the original width, possibly caused by the larger cohesive energy of Cu-Sn versus Sn-Sn.

Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation

We examine a method for computing the change in free energy with temperature of a crystalline solid. In the method, the free-energy difference between nearby temperatures is calculated via overlap-sampling free-energy perturbation with Bennetts optimization. Coupled to this is a harmonically targeted perturbation that displaces the atoms in a manner consistent with the temperature change, such that for a harmonic system, the free-energy difference would be recovered with no error. A series of such perturbations can be assembled to bridge larger gaps in temperature. We test this harmonically targeted temperature perturbation (HTTP) method through the application to the inverse-power soft potential, u(r)=ε(σ/r)(n), over a range of temperatures up to the melting condition. Three exponent values (n=12, 9, and 6) for the potential are studied with different crystal structures, specifically face-centered cubic (fcc), body-centered cubic (bcc), and hexagonal close packing. Absolute free energies (classical only) for each system are obtained by implementing the series to near-zero temperature, where the harmonic model becomes very accurate. The HTTP method is shown to provide very precise results, with errors in the free energy smaller than two parts in 10(5). An analysis of the thermodynamic stability of the various structures in the infinite-system limit confirms previous findings. In particular, for n=12 and 9, the fcc structure is stable for all temperatures up to melting, and for n=6, the bcc crystal becomes stable relative to fcc for temperatures above kT/ε=0.802±0.001. The effects of vacancies and other defects are not considered in the analysis.