Tests of Monte Carlo perturbation theory for the free energy of liquid copper

Monte Carlo perturbation theory, in which terms in the thermodynamic perturbation series are evaluated by Monte Carlo averaging, has potentially large advantages in efficiency for calculating free energies of liquids from ab initio potential surfaces. In order to test the accuracy of perturbation theory for liquid metals, a series of calculations has been done on liquid copper, modeled by an embedded atom potential. A simple 1/r(12) pair potential is used as the reference system. The free energy is calculated to third order in perturbation theory, and the results are compared to an exact formula. It is found that for optimal reference potential parameters, second order perturbation theory is essentially exact. Second and third order theories give accurate results for significantly nonoptimal reference parameters. The relation between perturbation theory and reweighting is discussed, and an approximate formula is derived that shows an exponential dependence of the efficiency of reweighting on the second order free energy correction. Finally, techniques for application to ab initio potentials are discussed. It is shown that with samples of 100 configurations, it is possible to obtain accuracy and precision at the level of approximately 1 meV/atom.     

Solvation thermodynamics: theory and applications

Potential distribution and coupling parameter theories are combined to interrelate previous solvation thermodynamic results and derive several new expressions for the solvent reorganization energy at both constant volume and constant pressure. We further demonstrate that the usual decomposition of the chemical potential into noncompensating energetic and entropic contributions may be extended to obtain a Gaussian fluctuation approximation for the chemical potential plus an exact cumulant expansion for the remainder. These exact expressions are further related to approximate first-order thermodynamic perturbation theory predictions and used to obtain a coupling-parameter integral expression for the sum of all higher-order terms in the perturbation series. The results are compared with the experimental global solvation thermodynamic functions for xenon dissolved in n-hexane and water (under ambient conditions). These comparisons imply that the constant-volume solvent reorganization energy has a magnitude of at most approximately kT in both experimental solutions. The results are used to extract numerical values of the solute-solvent mean interaction energy and associated fluctuation entropy directly from experimental solvation thermodynamic measurements.     

Thermodynamic perturbation theory of simple liquids: Monte Carlo simulation of a hard sphere system and the Helmholtz free energy of SW fluids

The Monte Carlo method is used to simulate similar sized hard sphere systems in a wide range of densities (from η 0.005 to 0.530 with a step of η = 0.005). The models are used to calculate the coefficients of the thermodynamic perturbation theory series for SW fluids up to the fourth order. The width of the attraction zone of the SW potential λ is varied from 1 to 2.5 sphere diameters. The analytical expressions approximating the obtained coefficients by polynomials with respect to the variables η and λ are determined. The absolute accuracy of the approximation is estimated to be better than ±0.001. All the necessary data for the calculation of the Helmholtz free energy of SW fluids up to the fourth-order perturbation theory are given.     

Non-hard sphere thermodynamic perturbation theory

A non-hard sphere (HS) perturbation scheme, recently advanced by the present author, is elaborated for several technical matters, which are key mathematical details for implementation of the non-HS perturbation scheme in a coupling parameter expansion (CPE) thermodynamic perturbation framework. NVT-Monte Carlo simulation is carried out for a generalized Lennard-Jones (LJ) 2n-n potential to obtain routine thermodynamic quantities such as excess internal energy, pressure, excess chemical potential, excess Helmholtz free energy, and excess constant volume heat capacity. Then, these new simulation data, and available simulation data in literatures about a hard core attractive Yukawa fluid and a Sutherland fluid, are used to test the non-HS CPE 3rd-order thermodynamic perturbation theory (TPT) and give a comparison between the non-HS CPE 3rd-order TPT and other theoretical approaches. It is indicated that the non-HS CPE 3rd-order TPT is superior to other traditional TPT such as van der Waals/HS (vdW/HS), perturbation theory 2 (PT2)/HS, and vdW/Yukawa (vdW/Y) theory or analytical equation of state such as mean spherical approximation (MSA)-equation of state and is at least comparable to several currently the most accurate Ornstein-Zernike integral equation theories. It is discovered that three technical issues, i.e., opening up new bridge function approximation for the reference potential, choosing proper reference potential, and/or using proper thermodynamic route for calculation of f(ex-ref), chiefly decide the quality of the non-HS CPE TPT. Considering that the non-HS perturbation scheme applies for a wide variety of model fluids, and its implementation in the CPE thermodynamic perturbation framework is amenable to high-order truncation, the non-HS CPE 3rd-order or higher order TPT will be more promising once the above-mentioned three technological advances are established.     

Thermodynamic perturbation theory of simple liquids. A hierarchy of equations for expansion coefficients

The properties of an expansion of the statistical sum of a simple liquid with respect to the potential in thermodynamic perturbation theory are analyzed. The coefficients of this expansion are determined by the unperturbed potential, depend on temperature and density, and can be calculated by means of mathematical modeling. It is shown here that the derivatives of these coefficients with respect to temperature and density are expressed through the higher expansion coefficient (these relations are usually called a hierarchy of equations). These coefficients determine the expansion of the Helmholtz free energy and RDF with respect to the perturbation potential. The thermodynamic characteristics of the system (entropy, internal energy, pressure) are expressed through both the differential relations for the Helmholtz free energy and the integral expressions containing RDF. It is found that the hierarchy of equations obtained in this work makes these different methods equivalent. This is important for the application of thermodynamic perturbation theory because it becomes unnecessary to model any other equilibrium properties of the system apart from the expansion coefficients.     

g-Hartree ab-initio calculations of the total energies of the four-electron isoelectronic series

The atomic total energies of the four-electron isoelectronic series are calculated by theg-Hartree 2nd order perturbation theory and the Dirac-Hartree-Fock-Rayleigh-Schrödinger 2nd order perturbation theory. The Coulomb correlation energy is calculated by these theories. The Breit interaction, vacuum polarization, self-energy and Q.E.D. corrections are calculated by the lowest order approximation. The results show that theg-Hartree approach overestimates the Coulomb correlation energy. However, with an increase of the nuclear charge, it overestimates much less. In the case of the Hartree-Fock 2nd order calculation, it underestimates the Coulomb correlation energy. With an increase of the nuclear charge, it underestimates much more.     

Can we predict the structure and stability of molecular crystals under increased pressure? First‐principles study of glycine phase transitions

The aim of this study was to determine whether the periodic density functional theory (DFT) calculations can be used for accurate prediction of the influence of the increased pressure on crystal structure and stability of molecular solids. To achieve this goal a series of geometry optimization and thermodynamic parameters calculations were performed for γ‐glycine and δ‐glycine structures at different pressure values using CASTEP program. In order to perform most accurate geometry optimization various exchange‐correlation functionals defined within generalized gradient approximation (GGA): PBE, PW91, RPBE, WC, PBESOL as well as defined within local density approximation (LDA), i.e. CAPZ, were tested. Geometry optimization was carried out using different dispersion correction methods (i.e. Grimme, TS, OBS) or without them. The linear response density functional perturbation theory (DFPT) was used to obtain the phonon dispersion curves and phonon density of states from which thermodynamic parameters, such as: free energy (ΔG), enthalpy (ΔH) and entropy (ΔS) were evaluated. The results of the geometry optimization depend strongly on the choice of the DFT functional. Calculated differences between the free energy of the studied polymorphic forms at the studied pressure values were consistent with experimental observations on their stability. The computations of thermodynamic properties not only confirmed the order of stability of two studied forms, but also enabled to predict the pressure at which this order is reversed. The results obtained in this study have proven that the plane‐wave basis set first principles calculations under periodic conditions is suitable for accurate prediction of crystal structure and stability. © 2018 Wiley Periodicals, Inc.     

Structure and thermodynamics of hard-core Yukawa fluids: thermodynamic perturbation approaches

The thermodynamic perturbation theories, which are based on the power series of a coupling constant (λ-expansion), have been proposed for studying the structural and thermodynamic properties of a hard-core Yukawa (HCY) fluid: one (A1-approximation) is the perturbation theory based on the hard-sphere repulsion as a reference system. The other (A2-approximation) is the perturbation theory based on the reference system which incorporates both the repulsive and short-range attractive interactions. The first-order mean-spherical approximation (FMSA) provided by Tang and Lu [J. Chem. Phys. 99, 9828 (1993)] has been employed for investigating the thermodynamic properties of a HCY fluid using the alternative method via the direct correlation function. The calculated results show that (i) the A1 and A2 approximations are in excellent agreements with previous computer simulation results in the literature and compare with the semi-empirical works of Shukla including the higher-order free energy terms, (ii) the A1 and A2 approximations are better than the FMSA and the mean-spherical approximation, (iii) the A2-approximation compares with the A1-approximation, even though the perturbation effect of an A2-approximation is much smaller than that of an A1-approximation, and that (iv) the FMSA study is particularly of advantage in providing the structure and thermodynamics in a simple and analytic manner.     

Second order of perturbation theory correction to the radial distribution function of a liquid with the interaction potential of hard spheres plus a square-well

A liquid with the interaction potential of hard spheres plus a square-well is analyzed using the Monte-Carlo technique. Numerical results for the perturbation theory series over a square-well potential are obtained in the form of the Barker and Henderson discrete representation. Approximating expressions for the correction to a liquid radial distribution function in the second order of perturbation theory are presented. The obtained results allow us to define this correction with a root-mean-square deviation of about 0.007. It is shown that the given approach provides a complete calculation in the second order of perturbation theory, and also the determination of the third order correction to the free energy for a liquid interacting with the potential of the Lennard-Jones type.     

Thermodynamics of the polarizable stockmayer fluid

Pressure and internal energy of a fluid composed of polarizable Stockmayer molecules have been calculated by a molecular dynamics computation as well as by thermodynamic perturbation theory. It is found that the effect of molecular polarizability is underestimated by perturbation theory.     

Evaluating perturbation contributions in SAFT models by comparing to molecular simulation of n-alkanes

SAFT models are generally written as a perturbation series of the Helmholtz energy with reciprocal temperature as the argument. The perturbation coefficients are then functions of density and molecular size. The variation of the perturbation coefficients with molecular size is given primarily by Wertheim's theory [6], [7], [8] and [9], but there may be additional variations as in the PC-SAFT model. In the present work, we compare the characterization of perturbation coefficients inferred from PC-SAFT to those derived from molecular simulations.The molecular simulations are based on Discontinuous Molecular Dynamics (DMD) and second order Thermodynamic Perturbation Theory (TPT). DMD simulation is applied to the repulsive part of the potential model with molecular details like fused hard spheres for the interaction sites and 110° bond angles. The thermodynamic effects of disperse attractions are treated by rigorous application of TPT. The present work re-examines the related work of Elliott and Gray [35] in the low density and critical regions, focusing on n-alkanes with carbon numbers ranging from 3 to 80.We find that SAFT theory overestimates the repulsive contribution (A₀) and underestimates the first order contribution (A₁) of Helmholtz energy relative to simulation. Nevertheless, the correlations are qualitatively reasonable. Significant inconsistencies arise when considering the second order contribution (A₂). For example, the PC-SAFT characterization of A₂ becomes larger than A₁ in the low density, long chain limit, raising concerns about the convergence of the series. Furthermore, fluctuations are underestimated in the critical region and overestimated in the liquid region. In each case, we can suggest improved characterizations. Altogether, these results suggest ways to modify the SAFT formalism to achieve greater consistency between atomistic and coarse-grained models.     

Fifth-order constant denominator perturbation theory within a localized bond model: Contributions from triple and quadruple excitations

The contributions of the triple and quadruple excitations to the fifth-order perturbation energy for the perturbation configuration interaction using localized orbitals (PCILO ) method are derived. This completes the development of a fifth-order constant denominator perturbation theory initiated in a previous paper [5] with the single and double excitations. This theory is tested on molecules containing strained ring geometries, stretched bonds, strongly polarized bonds, and delocalized pi systems: cases where the starting zero order reference wave function poorly describes the system. Although the perturbation expansions turn out to be slowly convergent, the Padé approximant taken from an energy series which itself is constructed from Padé approximants provides results accurate to within a few kilocalories/mole of benchmark calculations. Computational times as in the original PCILO procedure remain proportional to N³, where N is the number of bonds.     

Modeling intercalation through the sandwich‐type interactions between benzene and 14 polyaromatic molecules: DFT and ab initio results

We use second order Moller Plesset perturbation theory and several density functional theory methods to calculate the counterpoise corrected electronic interaction energies between benzene and a series of polyaromatic molecules. These systems serve as a simple model for DNA intercalation. We show that addition of nitrogen atoms to the polyaromatic molecules always increases sandwich‐type interactions, and that, of the density functional theory methods studied, only SVWN can mimic the interaction energies and optimal separations obtained with perturbation theory. SVWN reproduces the optimal molecular distances obtained with perturbation theory very well, and often comes within less than 10% of the interaction energy. © 2007 Wiley Periodicals, Inc. Int J Quantum Chem, 2008     

Phase Transitions of the Liquid–Liquid Type and a Change in the Particle Charge in Colloidal Solutions

The thermodynamic properties are studied for the solutions of charged colloidal particles with ionizable surface groups. The microscopic mechanism of microion binding at surface groups is considered. The free energy of the system in the parameter range where the usual theory of such solutions is inadequate (a range of practical interest) is calculated using the method of the thermodynamic perturbation theory. The first-order phase transition of the liquid–liquid type is shown to be possible; in this phase transition, a phase with a high concentration of colloidal particles that have a higher charge coexists with a phase with a lower concentration of particles that have a lower charge.     

A change of variable for the perturbation parameter in rayleigh–Schrödinger perturbation theory

The Feenberg–Goldhammer change of scale, whereby the H₀ in perturbation theory is replaced by (1/μ)H₀ with μ a scaling parameter, is shown to be equivalent to a change of variable for the perturbation parameter. A more general change of variable is shown to lead to a perturbation series with perturbation energies E₃, …, E_2n+1 equal to zero. The resulting energy through (2n + 1)th order has the same form as that found from the Brillouin–Wigner series by different methods.     

A completely analytic equation of state for the square-well chain fluid of variable well width

A completely analytic perturbation theory equation of state for the freely-jointed square-well chain fluid of variable well width (1 ≤ λ ≤ 2) is developed and tested against Monte Carlo simulation data. The equation of state is based on second-order Barker and Henderson perturbation theory to calculate the thermodynamic properties of the reference monomer fluid, and on first-order Wertheim thermodynamic perturbation theory to account for the connectivity of monomers to form chains. By using a recently developed real function expression for the radial distribution function of hard spheres in perturbation theory, we obtain analytic, closed form expressions for the Helmholtz free energy and the radial distribution function of square-well monomers of any well width. This information is used as the reference fluid in the perturbation theory of Wertheim to obtain an analytic equation of state, without adjustable parameters, that leads to good predictions of the compressibility factors and residual internal energies for 4-mer, 8-mer and 16-mer square-well fluids when compared with the simulation results. Further, very good results are obtained when this equation of state with temperature-independent parameters is used to correlate the vapor pressures and critical points of the linear alkanes from methane to n-decane.     

Thermodynamic properties of van der Waals fluids from Monte Carlo simulations and perturbative Monte Carlo theory

Computer simulations have been performed for fluids with van der Waals potential, that is, hard spheres with attractive inverse power tails, to determine the equation of state and the excess energy. On the other hand, the first- and second-order perturbative contributions to the energy and the zero- and first-order perturbative contributions to the compressibility factor have been determined too from Monte Carlo simulations performed on the reference hard-sphere system. The aim was to test the reliability of this "exact" perturbation theory. It has been found that the results obtained from the Monte Carlo perturbation theory for these two thermodynamic properties agree well with the direct Monte Carlo simulations. Moreover, it has been found that results from the Barker-Henderson [J. Chem. Phys. 47, 2856 (1967)] perturbation theory are in good agreement with those from the exact perturbation theory.     

Perturbation theory for liquids with a hard sphere plus square well potential

Numerical data have been obtained for a series of perturbation theory in the form of Barker-Henderson discrete representation. A hard sphere liquid (864 particles) was modeled by the Monte Carlo method for 106 values of occupancy from the range η = 0.005–0.530 (step 0.005); the well width was varied from 0 to 1.5 hard sphere diameters. For each occupancy, averaging was done over 60·10⁶ ensembles. Approximating analytical expressions are given as functions of well density and width for the first three terms of the free energy decomposition. For reduced temperature T* = 1.0, the standard error in the free energy calculation (three orders of perturbation theory if the above expressions are used) is of the order of ±0.001. Analysis of four-particle correlations revealed peaks at distances of 1.35 Å and 1.65 Å (η = 0.53) in the high-density region of hard sphere liquids, which are absent in crystals and dense liquids. It is analyzed whether the results are useful for realization of WCA theory of simple liquids in second order perturbation theory.