Simple and accurate scheme to compute electrostatic interaction: Zero-dipole summation technique for molecular system and application to bulk water

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 A for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.     

Application of the Wolf damped Coulomb method to simulations of SiC

A multibody interatomic potential is developed for bulk SiC using a modification of the Wolf et al. summation technique [D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999)] for the electrostatic interaction. The technique is modified to account for the short-range nonpoint charge effect. The nonelectrostatic interaction is modeled by a simple Morse-stretch term. This potential is then applied to beta-SiC to calculate various bulk properties using molecular dynamics simulations. The simulated x-ray diffraction pattern, radial distribution functions, lattice constant, elastic constants, and defect energy agree well with experimental data.     

Computer simulation of charged hard spherocylinders

In this work we present a computer simulation study of charged hard spherocylinders of aspect ratio L/sigma=5, using NVT and NPT Monte Carlo methods. Coulombic interactions are handled using the Wolf method [D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999)]. Thermodynamic and structural properties are in excellent agreement with the results obtained with the standard Ewald summation method. A partial prediction of the corresponding phase diagram is obtained by studying two isotherms of this system. The stability of the liquid crystalline phases is examined and compared with the phase diagrams of neutral hard spherocylinders and dipolar hard spherocylinders.     

Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics

We investigate pairwise electrostatic interaction methods and show that there are viable computationally efficient (O(N)) alternatives to the Ewald summation for typical modern molecular simulations. These methods are extended from the damped and cutoff-neutralized Coulombic sum originally proposed by Wolf et al. [J. Chem. Phys. 110, 8255 (1999)]. One of these, the damped shifted force method, shows a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. Comparisons were performed with this and other pairwise methods against the smooth particle-mesh Ewald summation to see how well they reproduce the energetics and dynamics of a variety of molecular simulations.     

Semiclassical calculation of nonadiabatic thermal rate constants: application to condensed phase reactions

The nonadiabatic transition state theory proposed recently by Zhao et al. [J. Chem. Phys. 121, 8854 (2004)] is extended to calculate rate constants of complex systems by using the Monte Carlo and umbrella sampling methods. Surface hopping molecular dynamics technique is incorporated to take into account the dynamic recrossing effect. A nontrivial benchmark model of the nonadiabatic reaction in the condensed phase is used for the numerical test. It is found that our semiclassical results agree well with those produced by the rigorous quantum mechanical method. Comparing with available analytical approaches, we find that the simple statistical theory proposed by Straub and Berne [J. Chem. Phys. 87, 6111 (1987)] is applicable for a wide friction region although their formula is obtained using Landau-Zener [Phys. Z. Sowjetunion 2, 46 (1932); Proc. R. Soc. London, Ser. A 137, 696 (1932)] nonadiabatic transition probability along a one-dimensional diffusive coordinate. We also investigate how the nuclear tunneling events affect the dependence of the rate constant on the friction.     

Contribution of inter- and intramolecular energy transfers to heat conduction in liquids

The molecular dynamics expression of heat flux, originally derived by Irving and Kirkwood [J. Chem. Phys. 18, 817 (1950)] for pairwise potentials, is generalized in this paper for systems with many-body potentials. The original formula consists of a kinetic part and a potential part, and the latter term is found in the present study to be expressible as a summation of contributions from all the many-body potentials defined in the system. The energy transfer among a set of sites for which a many-body potential is defined is discussed and evaluated by the rate of increase in the kinetic energy of each site due to the potential, and its accumulation over all the potentials in the system is shown to make up the potential part of the generalized expression. A molecular dynamics simulation for liquid n-octane was performed to demonstrate the applicability of the new expression obtained in this study to measure the heat flux and to elucidate the contributions of inter- and intramolecular potentials to heat conduction.     

Integral equation theory for correcting truncation errors in molecular simulations

Various strategies for correcting structural and energetic artefacts of molecular simulations with truncated potentials based on integral equation theory are described and applied to liquid water. The performance of the methods is examined for a range of cutoff distances and different shifted-force potentials. With the recently enhanced damped Coulomb potential (D. Zahn, B. Schilling, S.M. Kast, J. Phys. Chem. B, 106 (2002) 10725), parameterised and corrected by integral equation theory, radial distribution functions and excess internal energy very close to the Ewald simulation limit are obtained from a simulation with a cutoff distance of only 6 Å.     

Revisiting the Frenkel-Ladd method to compute the free energy of solids: the Einstein molecule approach

In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. [J. Chem. Phys. 112, 5339 (2000)]. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p*=11.54(4), which is slightly lower than the classical value of Hoover and Ree (p*=11.70) [J. Chem. Phys. 49, 3609 (1968)]. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.     

Rovibrational spectroscopy calculations of neon dimer using a phase space truncated Weyl-Heisenberg wavelet basis

In a series of earlier articles [B. Poirier J. Theor. Comput. Chem. 2, 65 (2003); B. Poirier and A. Salam J. Chem. Phys. 121, 1690 (2004); B. Poirier and A. Salam J. Chem. Phys. 121, 1740 (2004)], a new method was introduced for performing exact quantum dynamics calculations in a manner that formally defeats exponential scaling with system dimensionality. The method combines an optimally localized, orthogonal Weyl-Heisenberg wavelet basis set with a simple phase space truncation scheme, and has already been applied to model systems up to 17 degrees of freedom (DOF's). In this paper, the approach is applied for the first time to a real molecular system (neon dimer), necessitating the development of an efficient numerical scheme for representing arbitrary potential energy functions in the wavelet representation. All bound rovibrational energy levels of neon dimer are computed, using both one DOF radial coordinate calculations and a three DOF Cartesian coordinate calculation. Even at such low dimensionalities, the approach is found to be competitive with another state-of-the-art method applied to the same system [J. Montgomery and B. Poirier J. Chem. Phys. 119, 6609 (2003)].     

Soft core thermodynamics from self-consistent hard core fluids

In an effort to generalize the self-consistent Ornstein-Zernike approximation (SCOZA)-an accurate liquid state theory that has been restricted so far to hard core systems-to arbitrary soft core systems we study a combination of SCOZA with a recently developed perturbation theory. The latter was constructed by Ben-Amotz and Stell [J. Phys. Chem. B 108, 6877 (2004)] as a reformulation of the Weeks-Chandler-Andersen [J. Chem. Phys. 54, 5237 (1971)] perturbation theory directly in terms of an arbitrary hard sphere reference system. We investigate the accuracy of the combined approach for the Lennard-Jones fluid in comparison with simulation data and pure perturbation theory predictions and determine the dependence of the thermodynamic properties and the phase behavior on the choice of the effective hard core diameter of the reference system.     

The nature of the calculation of the pressure in molecular simulations of continuous models from volume perturbations

We consider some fundamental aspects of the calculation of the pressure from simulations by performing volume perturbations. The method, initially proposed for hard-core potentials by Eppenga and Frenkel [Mol. Phys.52, 1303 (1984)] and then extended to continuous potentials by Harismiadis et al. [J. Chem. Phys. 105, 8469 (1996)], is based on the numerical estimate of the change in Helmholtz free energy associated with the perturbation which, in turn, can be expressed as an ensemble average of the corresponding Boltzmann factor. The approach can be easily generalized to the calculation of components of the pressure tensor and also to ensembles other than the canonical ensemble. The accuracy of the method is assessed by comparing simulation results obtained from the volume-perturbation route with those obtained from the usual virial expression for several prototype fluid models. Monte Carlo simulation data are reported for bulk fluids and for inhomogeneous systems containing a vapor-liquid interface.     

Ultrafast vibrational spectroscopy of water and aqueous N-methylacetamide: Comparison of different electronic structure/molecular dynamics approaches

Kwac and Cho [J. Chem. Phys. 119, 2247 (2003)] have recently developed a combined electronic structure/molecular dynamics approach to vibrational spectroscopy in liquids. The method involves fitting ab initio vibrational frequencies for a solute in a cluster of solvent molecules to a linear combination of the electrostatic potentials on the solute atoms due to the charges on the solvent molecules. These authors applied their method to the N-methylacetamide-D/D(2)O system. We (S. A. Corcelli, C. P. Lawrence, and J. L. Skinner, [J. Chem. Phys. 120, 8107 (2004)]) have recently explored a closely related method, where instead of the electrostatic potential, the solute vibrational frequencies are fit to the components of the electric fields on the solute atoms due to the solvent molecules. We applied our method to the HOD/D(2)O and HOD/H(2)O systems. In order to make a direct comparison of these two approaches, in this paper we apply their method to the water system, and our method to the N-methylacetamide system. For the water system we find that the electric field method is superior to the potential approach, as judged by comparison with experiments for the absorption line shape. For the N-methylacetamide system the two methods are comparable.     

Theoretical scheme for the shear viscosity of Lennard-Jones fluids

We use the shear viscosity expression from the Enskog theory of dense gases in a perturbative scheme for the Lennard-Jones (LJ) fluid. This perturbative scheme is formulated by combining the analytic rational function approximation method of Bravo Yuste and Santos [Phys. Rev. A 43, 5418 (1991)] for the radial distribution function of hard-sphere fluids and the well known Mansoori-Canfield/Rasaiah-Stell perturbation theory to determine an effective diameter for the LJ fluid. The scheme is reliable on a wide range of temperatures and densities, and is very accurate around the critical point. Using this information, we build an accurate empirical formula for the shear viscosity in the liquid phase, which fits the recent data [K. Meier et al., J. Chem. Phys. 121, 3671 (2004)] in the whole simulation range.     

A semiclassical hybrid approach to many particle quantum dynamics

We analytically derive a correlated approach for a mixed semiclassical many particle dynamics, treating a fraction of the degrees of freedom by the multitrajectory semiclassical initial value method of Herman and Kluk [Chem. Phys. 91, 27 (1984)] while approximately treating the dynamics of the remaining degrees of freedom with fixed initial phase space variables, analogously to the thawed Gaussian wave packet dynamics of Heller [J. Chem. Phys. 62, 1544 (1975)]. A first application of this hybrid approach to the well studied Secrest-Johnson [J. Chem. Phys. 45, 4556 (1966)] model of atom-diatomic collisions is promising. Results close to the quantum ones for correlation functions as well as scattering probabilities could be gained with considerably reduced numerical effort as compared to the full semiclassical Herman-Kluk approach. Furthermore, the harmonic nature of the different degrees of freedom can be determined a posteriori by comparing results with and without the additional approximation.     

Analytical density-dependent representation of Hartree—Fock atomic potentials

A procedure to represent atomic electron charge densities [L. Fernandez Pacios, J. Phys. Chem., 95 , 10653 (1991); J. Phys. Chem., 96 , 7294 (1992)] is here generalized to obtain simple analytical functions for potential energy contributions. Based upon suitable functions to describe atomic electron densities in a physically meaningful form, the procedure is developed to define density-dependent analytical expressions for the electrostatic (classical) and exchange (quantum) potentials by means of proper approximate functionals. Calculations of correlation energies by using various density-functional approaches are also performed. The whole scheme is used to represent Hartree–Fock limit atomic wave functions by Clementi–Roetti. This way, a set of analytically simple, nonbasis set-dependent functions are defined with the aim to be further implemented in energy decomposition schemes for molecular interactions studies using atomic instead of electronic building blocks. © 1993 John Wiley & Sons, Inc.     

Accurate calculations of free-energy differences by the distribution method

We employ the strategy used in the successive umbrella sampling method [P. Virnau and M. Muller, J. Chem. Phys. 120, 10925 (2004)] to obtain the energy-difference distribution over its desired range. This is very helpful in calculating free-energy differences, where the source of the error is well recognized as the insufficient sampling over the relevant tail region in the energy-difference distribution. The distribution method proposed here employs the idea of restricting the sampling within an appropriate energy range, as was presented by Shing and Gubbins in their restricted umbrella sampling method [Mol. Phys. 46, 1109 (1982)]. We demonstrate the efficiency of the distribution method by calculating the free-energy difference of a model of harmonic oscillators where the systems exhibit nonoverlap features in their important phase spaces through the original Metropolis sampling. For this particular case, we show that the distribution method outperforms the free-energy perturbation method and even the Bennett's acceptance ratio method [J. Comput. Phys. 22, 245 (1976)] with the fastest convergence and the smallest relative errors. We further demonstrate the application of the distribution method with a simple point charge water model.     

Reconsidering an analytical gradient expression within a divide-and-conquer self-consistent field approach: exact formula and its approximate treatment

An analytical energy gradient formula for the density-matrix-based linear-scaling divide-and-conquer (DC) self-consistent field (SCF) method was proposed in a previous paper by Yang and Lee (YL) [J. Chem. Phys. 103, 5674 (1995)]. Since the formula by YL does not correspond to the exact gradient of the DC-SCF energy, we derive the exact formula by direct differentiation, which requires solving the coupled-perturbed equations while including the inter-subsystem coupling terms. Next, we present an alternative formula for approximately evaluating the DC-SCF energy gradient, assuming the variational condition for the subsystem density matrices. Numerical assessments confirmed that the DC-SCF energy gradient values obtained by the present formula are in reasonable agreement with the conventional SCF values when adopting a reliable buffer region. Furthermore, the performance of the present method was found to be better than that of the YL method.     

Ionic solvation studied by image-charge reaction field method

In a preceding paper [J. Chem. Phys. 131, 154103 (2009)], we introduced a new, hybrid explicit/implicit method to treat electrostatic interactions in computer simulations, and tested its performance for liquid water. In this paper, we report further tests of this method, termed the image-charge solvation model (ICSM), in simulations of ions solvated in water. We find that our model can faithfully reproduce known solvation properties of sodium and chloride ions. The charging free energy of a single sodium ion is in excellent agreement with the estimates by other electrostatics methods, while offering much lower finite-size errors. Similarly, the potentials of mean force computed for Na-Cl, Na-Na, and Cl-Cl pairs closely reproduce those reported previously. Collectively, our results demonstrate the superior accuracy of the proposed ICSM method for simulations of mixed media.     

Generalized principle of corresponding states and the scale invariant mean-field approach

In this paper we apply the relations between the critical points of the Lennard-Jones fluids and lattice gas model found in [V. L. Kulinskii, J. Phys. Chem. B 114, 2852 (2010)] to other short-ranged potentials like Buckingham and the Mie-potentials. The estimates for the corresponding critical point loci correlate quite satisfactory with the available numerical data for these potentials. The explanation for the correlation between the value of the second virial coefficient at the critical temperature and the particle volume found in [G. A. Vliegenthart and H. N. W. Lekkerkerker, J. Chem. Phys. 112, 5364 (2000)] is proposed. The connection of the stability of the liquid phase with the short range character of the potentials is discussed on the basis of the global isomorphism approach.