Thermostat Scheme Research Articles

Extensive numerical tests of leapfrog integrator in middle thermostat scheme in molecular simulations

Accurate and efficient integration of the equations of motion is indispensable for molecular dynamics (MD) simulations. Despite the massive use of the conventional leapfrog (LF) integrator in modern computational tools within the framework of MD propagation, further development for better performance is still possible. The alternative version of LF in the middle thermostat scheme (LF-middle) achieves a higher order of accuracy and efficiency and maintains stable dynamics even with the integration time stepsize extended by several folds. In this work, we perform a benchmark test of the two integrators (LF and LF-middle) in extensive conventional and enhanced sampling simulations, aiming at quantifying the time-stepsize-induced variations of global properties (e.g., detailed potential energy terms) as well as of local observables (e.g., free energy changes or bondlengths) in practical simulations of complex systems. The test set is composed of six chemically and biologically relevant systems, including the conformational change of dihedral flipping in the N-methylacetamide and an AT (Adenine-Thymine) tract, the intra-molecular proton transfer inside malonaldehyde, the binding free energy calculations of benzene and phenol targeting T4 lysozyme L99A, the hydroxyl bond variations in ethaline deep eutectic solvent, and the potential energy of the blue-light using flavin photoreceptor. It is observed that the time-step-induced error is smaller for the LF-middle scheme. The outperformance of LF-middle over the conventional LF integrator is much more significant for global properties than local observables. Overall, the current work demonstrates that the LF-middle scheme should be preferably applied to obtain accurate thermodynamics in the simulation of practical chemical and biological systems.

Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat.

An explicit treatment of electronic polarization is critically important to accurate simulations of highly charged or interfacial systems. Compared to the iterative self-consistent field (SCF) scheme, extended Lagrangian approaches are computationally more efficient for simulations that employ a polarizable force field. However, an appropriate thermostat must be chosen to minimize heat flow and ensure an equipartition of kinetic energy among all unconstrained system degrees of freedom. Here we investigate the effects of different thermostats on the simulation of condensed phase systems with the Drude polarizable force field using several examples that include water, NaCl/water, acetone, and an ionic liquid (IL) BMIM+/BF4-. We show that conventional dual-temperature thermostat schemes often suffer from violations of equipartitioning and adiabatic electronic state, leading to considerable errors in both static and dynamic properties. Heat flow from the real degrees of freedom to the Drude degrees of freedom leads to a steady temperature gradient and puts the system at an incorrect effective temperature. Systems with high-frequency internal degrees of freedom such as planar improper dihedrals or C-H bond stretches are most vulnerable; this issue has been largely overlooked in the literature because of the primary focus on simulations of rigid water molecules. We present a new temperature-grouped dual-Nosé-Hoover thermostat, where the molecular center of mass translations are assigned to a temperature group separated from the rest degrees of freedom. We demonstrate that this scheme predicts correct static and dynamic properties for all the systems tested here, regardless of the thermostat coupling strength. This new thermostat has been implemented into the GPU-accelerated OpenMM simulation package and maintains a significant speedup relative to the SCF scheme.

Unified Efficient Thermostat Scheme for the Canonical Ensemble with Holonomic or Isokinetic Constraints via Molecular Dynamics.

We have recently proposed a new unified theoretical scheme (the "middle" scheme) for thermostat algorithms for efficient and accurate configurational sampling of the canonical ensemble. In this paper, we extend the "middle" scheme to molecular dynamics algorithms for configurational sampling in systems subject to constraints. Holonomic constraints and isokinetic constraints are used for demonstration. Numerical examples indicate that the "middle" scheme presents a promising approach to calculate configuration-dependent thermodynamic properties and their thermal fluctuations.

A leap-frog algorithm-based efficient unified thermostat scheme for molecular dynamics

An important challenge for modern molecular simulations is to develop robust algorithms for molecular dynamics (MD), which employ a time interval (i.e., step size) of evolution as large as possible to achieve high efficiency while maintaining accuracy for sampling the distribution of configuration and momentum, e.g., for the canonical ensemble where the number of particles ( N ), the volume ( V ), and the temperature ( T ) are constant. Starting from the velocity-Verlet algorithm, in the recent paper we have proposed a unified scheme—the middle scheme for constructing efficient thermostatting integrators for configurational sampling of the canonical ensemble. In the present paper we extend the middle scheme to derive new thermostatting integrators based on the leap-frog algorithm. The leap-frog middle scheme is able to recover the exact marginal distribution of configuration and that of momentum in the harmonic limit, regardless of the time interval as long as the evolution of the MD trajectory is stable. It is also proved that the leap-frog middle scheme leads to the same marginal distribution of configuration as the velocity- Verlet middle scheme does even for general anharmonic systems. As illustrated in numerical examples (such as a 1-dimensional quartic potential, liquid water at room temperature with or without bond-length constraints, etc.), while both middle schemes perform better than conventional approaches (e.g., the side scheme and the default Langevin dynamics method in AMBER) and can calculate configuration-dependent properties with the same accuracy by using a larger time interval for from at least 4 times to even an order of magnitude, the leap-frog middle scheme produces a more accurate marginal distribution of momentum than that the velocity-Verlet middle scheme does. The sampling efficiency analysis by employing the characteristic correlation time step number (of the potential or of the Hamiltonian) suggests that both middle schemes share the same efficiency in sampling the configuration space and lead to similar efficiency in sampling the whole phase space (i.e., both configuration and momentum). When both only-configuration- dependent and only-momentum-dependent physical observables are considered, we suggest that the leap-frog middle scheme should in principle offer a better approach. When we study properties involving cross-terms of configuration and momentum, e.g., the fluctuation of the Hamiltonian (related to the heat capacity) and the energy current and its fluctuation (related to the thermal conductivity), analytic results in the harmonic limit and numerical data for anharmonic systems indicate that which of the two middle schemes performs better depends on the system and on the specific property. Note that the leap-frog middle scheme is of first-order accuracy for sampling the whole phase space while the velocity-Verlet middle scheme is of second-order accuracy. It is reasonable to choose the velocity-Verlet middle scheme to perform simulations for such properties for general molecular systems. Because either the leap-frog or the velocity-Verlet algorithm is the default integrator in molecular dynamics simulation software packages, the two unified middle schemes can offer novel, simple, accurate, and efficient thermostatting integrators for performing MD (or path integral MD (PIMD)) for the classical (or quantum) canonical ensemble.

Dynamic principle for ensemble control tools.

Dynamical equations describing physical systems in contact with a thermal bath are commonly extended by mathematical tools called "thermostats." These tools are designed for sampling ensembles in statistical mechanics. Here we propose a dynamic principle underlying a range of thermostats which is derived using fundamental laws of statistical physics and ensures invariance of the canonical measure. The principle covers both stochastic and deterministic thermostat schemes. Our method has a clear advantage over a range of proposed and widely used thermostat schemes that are based on formal mathematical reasoning. Following the derivation of the proposed principle, we show its generality and illustrate its applications including design of temperature control tools that differ from the Nosé-Hoover-Langevin scheme.

A unified thermostat scheme for efficient configurational sampling for classical/quantum canonical ensembles via molecular dynamics.

We show a unified second-order scheme for constructing simple, robust, and accurate algorithms for typical thermostats for configurational sampling for the canonical ensemble. When Langevin dynamics is used, the scheme leads to the BAOAB algorithm that has been recently investigated. We show that the scheme is also useful for other types of thermostats, such as the Andersen thermostat and Nosé-Hoover chain, regardless of whether the thermostat is deterministic or stochastic. In addition to analytical analysis, two 1-dimensional models and three typical real molecular systems that range from the gas phase, clusters, to the condensed phase are used in numerical examples for demonstration. Accuracy may be increased by an order of magnitude for estimating coordinate-dependent properties in molecular dynamics (when the same time interval is used), irrespective of which type of thermostat is applied. The scheme is especially useful for path integral molecular dynamics because it consistently improves the efficiency for evaluating all thermodynamic properties for any type of thermostat.

Water flow in carbon nanotubes: The effect of tube flexibility and thermostat.

Although the importance of temperature control in nonequilibrium molecular dynamics simulations is widely accepted, the consequences of the thermostatting approach in the case of strongly confined fluids are underappreciated. We show the strong influence of the thermostatting method on the water transport in carbon nanotubes (CNTs) by considering simulations in which the system temperature is controlled via the walls or via the fluid. Streaming velocities and mass flow rates are found to depend on the tube flexibility and on the thermostatting algorithm, with flow rates up to 20% larger when the walls are flexible. The larger flow rates in flexible CNTs are explained by a lower friction coefficient between water and the wall. Despite the lower friction, a larger solid-fluid interaction energy is found for flexible CNTs than for rigid ones. Furthermore, a comparison of thermostat schemes has shown that the Berendsen and Nosé-Hoover thermostats result in very similar transport rates, while lower flow rates are found under the influence of the Langevin thermostat. These findings illustrate the significant influence of the thermostatting methods on the simulated confined fluid transport.

An ergodic configurational thermostat using selective control of higher order temperatures.

The conventional Nosé-Hoover type deterministic thermostat scheme for controlling temperature by configurational variables (Braga-Travis (BT) thermostat) is non-ergodic for systems with a few degrees of freedom. While for the original Nosé-Hoover kinetic thermostat ergodicity has been achieved by controlling the higher order moments of kinetic energy, the issues of nonergodicity of BT thermostat persists. In this paper, we introduce two new measures of configurational temperature (second and third order) based on the generalized temperature-curvature relationship and obtain a family of deterministic thermostatting schemes by selectively (and simultaneously) controlling the different orders of temperatures through pseudo-friction terms. The ergodic characteristics of the proposed thermostats are tested using a single harmonic oscillator through statistical (normality of joint distributions at different Poincare sections) as well as dynamical tests (difference of the minimum and maximum largest Lyapunov exponent). Our results indicate that simultaneously controlling the first and the second order configurational temperatures (C(1,2) thermostat) is sufficient to make the dynamics ergodic. A 2000 particle Lennard-Jones system is subjected to (i) equilibrium and (ii) sudden temperature change under BT and C(1,2) thermostatting schemes. The C(1,2) thermostat is found to be more robust than the BT thermostat without increasing computational costs.

On the numerical treatment of dissipative particle dynamics and related systems

We review and compare numerical methods that simultaneously control temperature while preserving the momentum, a family of particle simulation methods commonly used for the modelling of complex fluids and polymers. The class of methods considered includes dissipative particle dynamics (DPD) as well as extended stochastic-dynamics models incorporating a generalized pairwise thermostat scheme in which stochastic forces are eliminated and the coefficient of dissipation is treated as an additional auxiliary variable subject to a feedback (kinetic energy) control mechanism. In the latter case, we consider the addition of a coupling of the auxiliary variable, as in the Nosé–Hoover–Langevin (NHL) method, with stochastic dynamics to ensure ergodicity, and find that the convergence of ensemble averages is substantially improved. To this end, splitting methods are developed and studied in terms of their thermodynamic accuracy, two-point correlation functions, and convergence. In terms of computational efficiency as measured by the ratio of thermodynamic accuracy to CPU time, we report significant advantages in simulation for the pairwise NHL method compared to popular alternative schemes (up to an 80% improvement), without degradation of convergence rate. The momentum-conserving thermostat technique described here provides a consistent hydrodynamic model in the low-friction regime, but it will also be of use in both equilibrium and nonequilibrium molecular simulation applications owing to its efficiency and simple numerical implementation.

Invited: Ab Initio Molecular Dynamics Study of Electrochemical Reaction at EC-Li2O2 Interfaces

Development of next-generation secondary batteries beyond Li-ion ones has been progressing all over the world. Li-air batteries are one of attracted candidates because of its large energy density. Organic liquids are normally used as a solvent of the electrolyte because of the lithium metal as an anode in the battery. However, it has been well known that some kinds of organic liquids such as ethylene carbonate (EC) are not stable on Li2O2 that grows up on the cathode in the discharging process. Previous ab initio simulation based on density functional theory (DFT) showed that a ring in the EC molecule is open by electron reduction [1]. Nevertheless, large dipole moment fluctuation by movement of EC molecule must change the electrostatic potential at the Li2O2-EC interface and will affect the chemical reaction process. In this presentation, we demonstrate EC adsorption and ring-opening process of the EC molecules on Li2O2 substrate by ab initio molecular dynamics. We find out that ring opening of EC molecule leads O-O bond breaking in an O2 2- unit in the Li2O2. We use openMX code [2] in order to carry out ab initiomolecular dynamics (AIMD). The exchange-correlation functional is treated within the generalized gradient approximation proposed by Perdew, Burke, and Ernzerhof (PBE). The norm-conserving pseudopotential scheme combined with pseudo-atomic basis sets imposing cutoff energy of 300 eV. The Brillouin-zone summation is evaluated only on Γ point. In our simulation model, there are thirty-two EC molecules on Li2O2 slab of six atomic layers. We first run classical molecular dynamics with force field named “pcff” [3] during 100 ps with a time step of 1 fs to relax the interface structure. Bottom layer of the Li2O2slab is fixed through the simulation. To eliminate electrostatic interaction between slab images in the repeated slab model, effective screening medium (ESM) method is adopted. Vacuum(ESM)/unit cell/metal(ESM) boundary conditions are imposed. No extra charge is introduced into the system. Simulation temperature keeps 350 K controlled by Nose-Hoover thermostat scheme. Time step is 1.0 fs. After the initial large fluctuation in the total energy and system temperature were converged (in 100 fs), we analyzed the trajectory. Figure 1 (a) shows the adsorbed structure of EC molecule onto the Li2O2. In the adsorbed EC molecule, sp2 orbital of carbon changes to sp3 like one with a surface oxygen of Li2O2. Figure 2 shows time evolution of the molecular charge estimated by Bader analysis.. The adsorption is reduction process. There are three other adsorbed EC molecules in the system. They are the same characteristics with each other. After adsorption, ring opening of one adsorbed EC molecule occurs around 600 fs in simulation time depicted in Fig. 1 (b). The reactant contains a single EC molecules and a O atom that is an adsorption site. 0.8 electron flows into the reactant as is shown in Fig.3. Therefore, the ring opening is one-electron reduction process. Interestingly, O2 2- bond in Li2O2 is broken and carboxylate structure with the EC molecule is made from O-. Finally, O-(CH2)2OCOO- is produced. The left O- exchange another O2 2-pair circled by white broken line in Fig.1 (b). It indicates that EC decomposition makes charging process impossible in Li-Air batteries. Reference [1] T. Laino andA. Curioni, Chem. Eur. J. 18, (2012) 3510.[2] See http://www.openmx-square.org/for the code, OpenMX, pseudo-atomic basis functions and pseudopotentials. [3] H. Sun, S. J. Mumby, J. R. Maple, A. T. J. Hagler, Am. Chem. Soc. 116 , (1994), 2978.

Effect of the thermostat in the molecular dynamics simulation on the folding of the model protein chignolin

Molecular dynamics simulations of the model protein chignolin with explicit solvent were carried out, in order to analyze the influence of the Berendsen thermostat on the evolution and folding of the peptide. The dependence of the peptide behavior on temperature was tested with the commonly employed thermostat scheme consisting of one thermostat for the protein and another for the solvent. The thermostat coupling time of the protein was increased to infinity, when the protein is not in direct contact with the thermal bath, a situation known as minimally invasive thermostat. In agreement with other works, it was observed that only in the last situation the instantaneous temperature of the model protein obeys a canonical distribution. As for the folding studies, it was shown that, in the applications of the commonly utilized thermostat schemes, the systems are trapped in local minima regions from which it has difficulty escaping. With the minimally invasive thermostat the time that the protein needs to fold was reduced by two to three times. These results show that the obstacles to the evolution of the extended peptide to the folded structure can be overcome when the temperature of the peptide is not directly controlled.

High-Performance Scalable Molecular Dynamics Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD

Incorporating the influence of induced polarization in large-scale atomistic molecular dynamics (MD) simulations is a critical challenge in the progress toward computations of increased accuracy. One computationally efficient treatment is based on the classical Drude oscillator in which an auxiliary charged particle is attached by a spring to each nucleus. Here, we report the first implementation of this model in the program NAMD. An extended Lagrangian dynamics with a dual-Langevin thermostat scheme applied to the Drude-nucleus pairs is employed to efficiently generate classical dynamic propagation near the self-consistent field limit. Large-scale MD simulations based on the Drude polarizable force field scale very well on massively distributed supercomputing platforms, the computational demand being only about 50-100% higher than for nonpolarizable models. As an illustration, a large-scale 150 mM NaCl aqueous salt solution is simulated, and the calculated ionic conductivity is shown to be in excellent agreement with experiment.

Molecular Dynamics Simulation of a Core−Shell Structured Metallic Nanoparticle

This study investigates the heating and cooling of a shell−core structured Ni−Al functional particle through a molecular dynamics simulation. The functional nanoparitcle is constructed from their perfect crystal measuring 5.6 nm in diameter with a shell thickness of 0.5 nm. The simulation is carried out within a Nose−Hoover thermostat scheme of the NVT canonical ensemble and performed by an Embedded-Atom-Method (EAM) force field. The thermodynamic properties and structure evolution during continuous heating and cooling processes are investigated through the characterization of the gyration radius, radial distribution function (RDF), atom number distribution, mean-square-distance (MSD), and layered potential energy distribution. Some unique behaviors related to nanometer scale functional particles are identified that include two-way diffusion of aluminum and nickel atoms, enhanced melting of the aluminum core through the shell structure, glass phase formation for the fast cooling rate, and the Ll2 Ni3Al ph...

Thermal stability of single and multi-walled carbon nanotubes

This paper investigates the thermal stability of single and multi-walled carbon nanotubes (CNTs). The modeling and simulation of the problem are carried out using molecular dynamics simulation incorporating the Berendsen thermostat scheme for both open-ended and close-capped CNTs. It is found that single-walled CNTs are thermally more stable than multi-walled CNTs. The simulations also indicate that shorter CNTs are able to withstand higher thermal loads. CNTs with a larger diameter are also more resistant to thermal loads.