NVE Molecular Dynamics Research Articles

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions.

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, rc. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, Es, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 1010 are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h4) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

Shadow Hamiltonian in classical NVE molecular dynamics simulations: A path to long time stability.

The shadow energy, Es, is the conserved quantity in microcanonical ensemble (NVE) molecular dynamics simulations carried out with the position Verlet central-difference algorithm. A new methodology for calculating precise and accurate values of Es is presented. It is shown for the first time that Es rather than E is constant during structural changes occurring within a supercooled liquid. It is also explained how to prepare and conduct microsecond range bulk-phase NVE simulations with essentially zero energy drift without the need for thermostating. The drift is analyzed with block averaging and new drift functions of the shadow energy. With such minimal drift, extremely small and accurate standard errors in the mean for quantities like Es, E, and temperature, T, can be obtained. Values of the standard error for Es of ≈10-10 in molecule-based reduced units can be routinely achieved for simulations of 108 time steps. This corresponds to a simulation temperature drift of ≈10-6 K/μs, six orders of magnitude smaller than generally considered to be acceptable for protein simulations. We also show for the first time how these treatments can be extended with no loss of accuracy to polyatomic systems with both flexible degrees of freedom and arbitrary geometric constraints imposed via the SHAKE algorithm. As a bonus, estimates of simulation-average kinetic and total energies from high order velocity expressions can be obtained to a good approximation from 2nd order velocities and the average mean square force (for polyatomics, this refers to per site, including any constraint forces).

Molecular dynamics simulation investigation of the microwave heating supercritical water

Nonequilibrium NVE molecular dynamics (MD) simulations were performed to investigate the heating process of supercritical water (SCW) by the microwaves (the field intensities are from 3 × 105V/m to 3 × 109V/m with the frequencies of 915 MHz and 2.45 GHz, respectively). The energy distribution, microscopic molecular structure, diffusivity and hydrogen bond dynamics were analyzed. Microwaves field dependence on the rotational and translational coupling mechanisms of water molecules have been detected. The microwave energy has been converted to the kinetic and intermolecular energies of water. There are about 40% absorbed microwave energy is stored as the intermolecular energy, and the proportion decline at higher field intensity. The diffusion of supercritical water is non-linearly fluctuation due to the coupling of rotation and translation motions. And this work may provide a theoretical basis to improve the efficiency of industrial application with microwave and SCW.

SQDFT: Spectral Quadrature method for large-scale parallel [formula omitted] Kohn–Sham calculations at high temperature

We present SQDFT: a large-scale parallel implementation of the Spectral Quadrature (SQ) method for O(N) Kohn–Sham Density Functional Theory (DFT) calculations at high temperature. Specifically, we develop an efficient and scalable finite-difference implementation of the infinite-cell Clenshaw–Curtis SQ approach, in which results for the infinite crystal are obtained by expressing quantities of interest as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw–Curtis quadrature rules. We demonstrate the accuracy of SQDFT by showing systematic convergence of energies and atomic forces with respect to SQ parameters to reference diagonalization results, and convergence with discretization to established planewave results, for both metallic and insulating systems. We further demonstrate that SQDFT achieves excellent strong and weak parallel scaling on computer systems consisting of tens of thousands of processors, with near perfect O(N) scaling with system size and wall times as low as a few seconds per self-consistent field iteration. Finally, we verify the accuracy of SQDFT in large-scale quantum molecular dynamics simulations of aluminum at high temperature. Program summaryProgram Title: SQDFTProgram Files doi:http://dx.doi.org/10.17632/h4jnrmh9v3.1Licensing provisions: GNU General Public License 3 (GPL)Programming language: C/C++External routines/libraries: MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/) or later.Nature of problem: Quantum mechanical calculation of static and dynamic properties of condensed matter at high temperature in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method: High-order finite-difference discretization. Calculation of the electronic ground state using the self-consistent field (SCF) iteration in conjunction with Periodic Pulay extrapolation/mixing. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Solution of the Poisson equation using the Alternating Anderson–Richardson (AAR) method. Calculation of the electron density, band structure energy, electronic entropy energy, and atomic forces using the infinite-cell Clenshaw–Curtis Spectral Quadrature (SQ) approach, in which results for the infinite crystal are obtained by expressing quantities as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw–Curtis quadrature rules. NVE molecular dynamics using the leapfrog method and NVT molecular dynamics using the Verlet algorithm with the Nose–Hoover thermostat. Extrapolation of electron density between molecular dynamics steps. Parallelization via domain decomposition.Restrictions: Cubical domain. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections.

SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Extended systems

As the second component of SPARC (Simulation Package for Ab-initio Real-space Calculations), we present an accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory (DFT) for extended systems. Specifically, employing a local formulation of the electrostatics, the Chebyshev polynomial filtered self-consistent field iteration, and a reformulation of the non-local force component, we develop a finite-difference framework wherein both the energy and atomic forces can be efficiently calculated to within desired accuracies in DFT. We demonstrate using a wide variety of materials systems that SPARC achieves high convergence rates in energy and forces with respect to spatial discretization to reference plane-wave result; exponential convergence in energies and forces with respect to vacuum size for slabs and wires; energies and forces that are consistent and display negligible ‘egg-box’ effect; accurate properties of crystals, slabs, and wires; and negligible drift in molecular dynamics simulations. We also demonstrate that the weak and strong scaling behavior of SPARC is similar to well-established and optimized plane-wave implementations for systems consisting up to thousands of electrons, but with a significantly reduced prefactor. Overall, SPARC represents an attractive alternative to plane-wave codes for performing DFT simulations of extended systems. Program summaryProgram title: SPARCCatalogue identifier: AFBR_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFBR_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU GPL v3No. of lines in distributed program, including test data, etc.: 93822No. of bytes in distributed program, including test data, etc.: 1386659Distribution format: tar.gzProgramming language: C/C++.Computer: Any system with C/C++ compiler.Operating system: Linux.RAM: Problem dependent. Ranges from 80 GB to 800 GB for a system with 2500 electrons.Classification: 7.3.External routines: PETSc 3.5.3 (http://www.mcs.anl.gov/petsc), MKL 11.2 (https://software.intel.com/en-us/intel-mkl), and MVAPICH2 2.1 (http://mvapich.cse.ohio-state.edu/).Does the new version supersede the previous version?: YesNature of problem: Calculation of the static and dynamic properties of isolated and extended systems in the framework of Kohn–Sham Density Functional Theory (DFT).Solution method:High-order finite-difference discretization. Local reformulation of the electrostatics in terms of the electrostatic potential and pseudocharge densities. Application of Bloch-periodic and zero-Dirichlet boundary conditions on the orbitals in the direction of periodicity and vacuum, respectively. Application of periodic and Dirichlet boundary conditions on the electrostatic potential in the direction of periodicity and vacuum, respectively. Integration over the Brillouin zone for extended systems using the Monkhorst–Pack grid. Calculation of the electronic ground-state using the Chebyshev polynomial filtered self-consistent field iteration in conjunction with Anderson based extrapolation/mixing schemes. Reformulation of the non-local component of the force. Geometry optimization using the Polak–Ribiere variant of non-linear conjugate gradients with secant line search. NVE molecular dynamics using the leapfrog method. Parallelization via domain decomposition and over Brillouin zone integration.Reasons for new version:To enable the study of extended systems like crystals, slabs, and wires using SPARC.Summary of revisions:Incorporated the ability to study the static and dynamic properties of crystals, slabs, and wires.Restrictions:System size less than ∼4000 electrons. Local Density Approximation (LDA). Troullier–Martins pseudopotentials without relativistic or non-linear core corrections. Domain has to be cuboidal.Running time:Problem dependent. Timing results for selected examples provided in the paper.

All-atom molecular dynamics simulations of multiphase segregated polyurea under quasistatic, adiabatic, uniaxial compression

We present an approach to model the mechanical response of multiphase segregated polyurea under large quasistatic uniaxial compression. The approach is based on a two-step procedure. We first conduct all-atom NVE molecular dynamics simulations of the fully-mixed and hard phases of polyurea. We then put forth a model of hard-phase activation based on Taylor, or meanfield, averaging and compute the composite response. The predictions of the present approach, with or without Taylor averaging, show remarkable agreement with data from plate-impact experiments.

Enhanced Gas Absorption in the Ionic Liquid 1-n-Hexyl-3-methylimidazolium Bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]) Confined in Silica Slit Pores: A Molecular Simulation Study

Two-dimensional NPxyT and isostress-osmotic (N2PxyTf1) Monte Carlo simulations were used to compute the density and gas absorption properties of the ionic liquid (IL) 1-n-hexyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]) confined in silica slit pores (25-45 Å). Self-diffusivity values for both gas and IL were calculated from NVE molecular dynamics simulations using both smooth and atomistic potential models for silica. The simulations showed that the molar volume of [hmim][Tf2N] confined in 25-45-Å silica slit pores is 12-31% larger than that of the bulk IL at 313-573 K and 1 bar. The amounts of CO2, H2, and N2 absorbed in the confined IL are 1.1-3 times larger than those in the bulk IL because of the larger molar volume of the confined IL compared to the bulk IL. The CO2, N2, and H2 molecules are generally absorbed close to the silica wall where the IL density is very low. This arrangement causes the self-diffusivities of these gases in the confined IL to be 2-8 times larger than those in the bulk IL at 298-573 K. The solubilities of water in the confined and bulk ILs are similar, which is likely due to strong water interactions with [hmim][Tf2N] through hydrogen bonding, so that the molar volume of the confined IL plays a less important role in determining the H2O solubility. Water molecules are largely absorbed in the IL-rich region rather than close to the silica wall. The self-diffusivities of water correlate with those of the confined IL. The confined IL exhibits self-diffusivities larger than those of the bulk IL at lower temperatures, but smaller than those of the bulk IL at higher temperatures. The findings from our simulations are consistent with available experimental data for similar confined IL systems.

MOLECULAR DYNAMICS SIMULATIONS OF CARBON NANOTUBE-BASED OSCILLATORS HAVING TOPOLOGICAL DEFECTS

Effect of vacancy and Stone–Wales defects on the oscillatory behavior of (5,5)/(10,10) carbon nanotube-based oscillator are studied using NVE molecular dynamics simulations. Results show that defects reduce stability of the oscillators. Effect of single vacancy defect on stability is very small, whereas Stone–Wales defect considerably reduces the stability thereby damping the oscillations quickly. Further increase in density of vacancy defects causes a monotonic decrease of stability of oscillator. In all cases the initial temperature (1 and 300 K) had almost no effect on the oscillation stability.

Molecular dynamics simulation of diffusion of O 2 and CO 2 in blends of amorphous poly(ethylene terephthalate) and related polyesters

The diffusion of small molecules through polymers is important in many areas of polymer science, such as gas barrier and separation membrane materials, polymeric foams, and in the processing and properties of polymers. Molecular dynamics simulation techniques have been applied to study the diffusion of oxygen and carbon dioxide as small molecule penetrants in models polyester blends of bulk amorphous poly(ethylene terephthalate) and related aromatic polyesters. A bulk amorphous configuration with periodic boundary conditions was generated into a unit cell whose dimensions were determined for each of the simulated polyester blends in the cell having the experimental density. The diffusion coefficients for O 2 and CO 2 were determined via NVE molecular dynamics simulations using the Dreiding 2.21 molecular mechanics force field over a range of temperatures (300, 500 and 600 K) using up to 40 ns simulation time. We have focussed on the influence of the temperature, polymer dynamics, density and free volume distribution on the diffusion properties. Correlation of diffusion coefficients with free volume distribution was found.

Molecular dynamics simulation of diffusion of O 2 and CO 2 in amorphous poly(ethylene terephthalate) and related aromatic polyesters

The diffusion of small molecules through polymers is important in many areas of polymer science, such as gas barrier and separation membrane materials, polymeric foams, and in the processing and properties of polymers. Molecular simulation techniques have been applied to study the diffusion of oxygen and dioxide of carbon as small molecule penetrants in models of bulk amorphous poly(ethylene terephthalate) and related aromatic polyesters. A bulk amorphous configuration with periodic boundary conditions is generated into a unit cell whose dimensions are determined for each of the simulated aromatic polyesters in the cell to have the experimental density. The aim for this research is to explore and investigate the diffusion of gases through bulk amorphous poly(ethylene terephthalate) and related aromatic polyesters. The diffusion coefficients for O 2 and CO 2 were determined via NVE molecular dynamics simulations using the Dreiding 2.21 molecular mechanics force field over a range of temperatures (300, 500 and 600 K) using up to 30 ns simulation time. We have focussed on the influence of the temperature, polymer dynamics, number of aromatic rings, ortho-, meta-, para-isomers, density and free volume distribution on the diffusion properties. Correlation of diffusion coefficients with free volume, temperature, number of aromatic rings, ortho-, meta- and para-isomers was found.

Simulation of diffusion of O 2 and CO 2 in amorphous poly(ethylene terephthalate) and related alkylene and isomeric polyesters

The diffusion of small molecules through polymers is important in many areas of polymer science, such as gas barrier and separation membrane materials, polymeric foams, and in the processing and properties of polymers. Molecular simulation techniques have been applied to study the diffusion of oxygen and carbon dioxide as small molecule penetrants in models of bulk amorphous poly(ethylene terephthalate) (PET) and related alkylene and isomeric polyesters. A bulk amorphous configuration with periodic boundary conditions made into a unit cell whose dimensions were determined for each of the simulated polyesters in the cell having the experimental density. The diffusion coefficients for O 2 and CO 2 were determined via NVE molecular dynamics simulations using the Dreiding 2.21 molecular mechanics force field over a range of temperatures (300, 500 and 600 K) using up to 3 ns simulation time. We have focussed on the influence of the temperature, polymer dynamics, number of CH 2 groups, density and free volume distribution on the diffusion properties. Correlation of diffusion coefficients with free volume and number of CH 2 groups was found.

Vibrational energy relaxation in classical fluids. I. High-frequency spectra in gases

We study high-frequency spectra of classical gases interacting via Lennard-Jones and similar potentials with applications to the problem of vibrational energy relaxation. We derive expressions for the spectra which are asymptotically exact at high frequencies and compare these with classical trajectory and molecular dynamics simulations. Trajectory calculations are in quantitative agreement with the theory, but conventional NVE molecular dynamics simulations are qualitatively in error at high frequencies due to insufficient simulation time.

Application of the hypervolume Monte Carlo methods to a molten ionic system

The Hyper Volume Monte Carlo Method (HVMC) proposed by Fernandes and Ramalho initially developed for the NVE ensemble, has recently been extended to the NpH ensemble. The methods have been applied to simple molecular and magnetic systems. The purpose of this work is to apply the HVMC method to an ionic system in the NVE and NpH ensembles and to compare the results with a NVE molecular dynamics (MD) simulation. The Born–Mayer–Huggins pair potential with parameters for KCl has been used. The long range electrostatic interactions were calculated using the modified Ewald sum of Adams and Dubey (D. Adams, G. Dubey, CCP5 Newsletter 22 (1986) 39; D. Adams, G. Dubey, J. Comp. Phys. 72 (1987) 156). A molten KCl system at 1045, 1306, 2090 and 2874 K and at several densities has been studied using 216 ions. At 1306 K, a number-dependence analysis using 64, 216, 512 and 1000 ions has been performed. The thermodynamics averages obtained by NVE MD simulations were used as inputs for the HVMC method. The thermal, structural properties and speed distribution functions have been calculated and an excellent agreement has been found between the MD and the HVMC methods at different temperatures and numbers of ions.

A molecular dynamics study of a bi-disperse liquid crystal mixture using a generalized Gay–Berne potential

We report results from a constant NVE molecular dynamics simulation of a 50:50 mixture of Gay–Berne rods with axial ratios of 3·5 to 1 and 3 to 1. The simulation makes use of a recently developed extension of the Gay–Berne potential which allows simulation of non-equivalent anisotropic molecules with attractive interactions. There is clear evidence of pre-smectic ordering in the nematic phase. The system undergoes a well defined transition from the nematic phase to a randomly mixed smectic phase. There is some local compositional ordering within the smectic layers but no evidence of demixing. In both the nematic and smectic phases it is found that the longer molecules have a higher ⟨S⟩ order parameter than the shorter molecules. This difference in order parameter is at its greatest a little below the isotropic-nematic transition temperature and decreases as the smectic phase is approached.

A molecular dynamics study of a bi-disperse liquid crystal mixture using a generalized Gay-Berne potential

We report results from a constant NVE molecular dynamics simulation of a 50:50 mixture of Gay-Berne rods with axial ratios of 35 to 1 and 3 to 1. The simulation makes use of a recently developed extension of the Gay-Berne potential which allows simulation of non-equivalent anisotropic molecules with attractive interactions. There is clear evidence of pre-smectic ordering in the nematic phase. The system undergoes a well defined transition from the nematic phase to a randomly mixed smectic phase. There is some local compositional ordering within the smectic layers but no evidence of demixing. In both the nematic and smectic phases it is found that the longer molecules have a higher [S] order parameter than the shorter molecules. This difference in order parameter is at its greatest a little below the isotropic-nematic transition temperature and decreases as the smectic phase is approached.

A Hybrid Density Functional-Classical Molecular Dynamics Simulation of a Water Molecule in Liquid Water

A hybrid NVE Molecular Dynamics simulation of liquid water is presented using a coupled Density Functional/Molecular Mechanics hamiltonian. The quantum subsystem is a single water molecule described by means of a triple-zeta quality basis set with polarization orbitals on oxygen and hydrogen atoms. Non-local exchange-correlation corrections are included self-consistently. The classical system is constituted by 128 classical TIP3P water molecules. Results are in reasonable agreement with experimental data and in particular a good description of the solute polarization is obtained. Large fluctuations of the instantaneous value of the dipole moment of the quantum molecule are predicted.

Hypervolumes in microcanonical Monte Carlo

A Monte Carlo method to perform microcanonical simulations by sampling the configurational and momenta spaces is presented. The technique was inspired by the method of hypervolumes for calculating the entropy in a microcanonical ensemble. Although this method is strictly proven in the thermodynamic limit, the hypervolume Monte Carlo (HVMC) algorithm, presented in this article, works well with a relatively small number of particles. In contrast to other algorithms for microcanonical Monte Carlo simulations, the HVMC method does not involve previous integrations over the momenta space or demons. Therefore, it can be used with any form of Hamiltonian. Thermal and structural properties for the Lennard-Jones system obtained by NVE molecular dynamics are compared with results from the HVMC method. The agreement is excellent. Additionally, the method provides the speed distribution functions of the system which are, also, in excellent agreement with the results from molecular dynamics. A discussion of the HVMC method in the context of the statistical mechanical theory of the microcanonical ensemble is presented.

Comparison of the Structure of Liquid Amides as Determined by Diffraction Experiments and Molecular Dynamics Simulations

The intermodular structures of liquid formamide, N-methylformamide and N,N-dimethylformamide at room temperature are studied by means of NVE molecular dynamics computer simulations. Newly developed flexible models are used. X-ray and neutron weighted structure and total radial pair distribution functions are computed from the simulated site-site pair distribution functions. They are compared with experimental results. The agreement is usually satisfactory as far as heavy atom pairs are concerned while the lengths of the hydrogen bonds are found to be systematically too long in the simulations.

Dynamics of an ammonium ion in water: Molecular dynamics simulation

The translational and rotational dynamics of an aqueous ammonium ion are examined using an NVE molecular dynamics computer simulation and a rigid model of water. The linear momentum, angular momentum and reorientational autocorrelation functions of the ammonium ion have been measured. From these, the memory functions for the linear and angular momentum correlation functions are calculated. Reorientation of the ammonium ion about its center of mass shows unambiguous departure from Debye rotational behavior. The translational and rotational diffusion constants have been calculated.

Extensions of the molecular dynamics simulation method. II. Isothermal systems

In this paper we systematically develop two general classes of equations of motion for performing molecular dynamics simulations at constant temperature. One member of each class is found to be identical to one of the two commonly used methods for performing such simulations. Of those two, one, namely the momentum scaling procedure, is proven to have no effect on the static properties of the simulated system. In addition, we have tested the two methods via equilibrium molecular dynamics on the Lennard-Jones fluid and find the resulting velocity autocorrelation functions and self-diffusion coefficients to be in good agreement with those determined by the standard NVE molecular dynamics method.