Stability analysis of multiple nonequilibrium fixed points in self-consistent electron transport calculations

We present a method to perform stability analysis of nonequilibrium fixed points appearing in self-consistent electron transport calculations. The nonequilibrium fixed points are given by the self-consistent solution of stationary, nonlinear kinetic equation for single-particle density matrix. We obtain the stability matrix by linearizing the kinetic equation around the fixed points and analyze the real part of its spectrum to assess the asymptotic time behavior of the fixed points. We derive expressions for the stability matrices within Hartree-Fock and linear response adiabatic time-dependent density functional theory. The stability analysis of multiple fixed points is performed within the nonequilibrium Hartree-Fock approximation for the electron transport through a molecule with a spin-degenerate single level with local Coulomb interaction.     

A gentler approach to RNEMD: nonisotropic velocity scaling for computing thermal conductivity and shear viscosity

We present a new method for introducing stable nonequilibrium velocity and temperature gradients in molecular dynamics simulations of heterogeneous systems. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods which use momentum exchange swapping moves. The standard swapping moves can create nonthermal velocity distributions and are difficult to use for interfacial calculations. By using nonisotropic velocity scaling (NIVS) on the molecules in specific regions of a system, it is possible to impose momentum or thermal flux between regions of a simulation while conserving the linear momentum and total energy of the system. To test the method, we have computed the thermal conductivity of model liquid and solid systems as well as the interfacial thermal conductivity of a metal-water interface. We find that the NIVS-RNEMD improves the problematic velocity distributions that develop in other RNEMD methods.     

The influence of thermostats and manostats on reverse nonequilibrium molecular dynamics calculations of fluid viscosities

Reverse nonequilibrium molecular dynamics to calculate the shear viscosity of Lennard-Jones liquids was extended to simulations at constant number of particles, constant volume, and constant pressure using a Berendsen thermostat and a Berendsen manostat. Using additional systems such as water and hexane, we also report on the performance of shear viscosity calculations of systems with electrostatic and nontrivial intramolecular interactions when a manostat is applied. We compare the shear viscosities of simulations using no coupling, only temperature coupling, and temperature and pressure coupling and characterize discrepancies, where observed. From this, we deduce guidelines for when and how manostats can be usefully applied in reverse nonequilibrium simulations.     

Universality of charge-carrier transport in molecularly doped polymers

The charge-carrier transport model based on the multiple-trapping quasi-band theory with the Gaussian or exponential energy distributions of traps for the two-layer structure of the molecularly doped polymer sample is used to explain the experimentally observed constancy of the flat plateau on time-of-flight curves in a wide electric-field range. The constant shape of the time-of-flight curve under the condition of nonequilibrium transport is observed for the exponential, rather than Gaussian, energy trap distribution. This observation may be used to distinguish between these trap distributions. Finding the true Poole–Frenkel constant requires that the nonequilibrium transport in molecularly doped polymers be taken into account during treatment of the data on the field dependence of mobility.     

Multiple Markov transition matrix method: Obtaining the stationary probability distribution from multiple simulations

We herein propose the multiple Markov transition matrix method (MMMM), an algorithm by which to estimate the stationary probability distribution from independent multiple molecular dynamics simulations with different Hamiltonians. Applications to the potential of mean force calculation in combination with the umbrella sampling method are presented. First, the performance of the MMMM is examined in the case of butane. Compared with the weighted histogram analysis method (WHAM), the MMMM has an advantage with respect to the reasonable evaluation of the stationary probability distribution even from nonequilibrium trajectories. This method is then applied to Met‐enkephalin nonequilibrium simulation. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Dynamic of Exchange Sorption of Transition Metal Ions on Phosphate Cationite KFP-12

The dynamics of ion-exchange sorption of double-charged metal ions on a phosphate cationite KFP-12 was studied. Theoretical calculations of stationary fronts were carried out, the motion velocity of the stationary fronts was calculated, and the rate constant of sorption was evaluated.     

Equilibrium/nonequilibrium initial configurations in forward/reverse electron transfer within mixed-metal hemoglobin hybrids

In a protein-protein electron transfer (ET) photocycle, the "forward" ET reaction is initiated with the excited complex, [3DA], in an equilibrium ensemble of configurations, the majority of which exhibit less than the maximal ET matrix element. In contrast, the charge-separated intermediate complex is formed in a nonequilibrium set of configurations with maximal ET matrix elements and would be expected to return to the ground state with the largest rate constant possible unless conformational interconversion first "breaks the connection" and the complex converts to less-reactive substates. According to this analysis, the forward and back ET reactions should show a differential response to viscosity, and the latter could even show an increased rate constant under conditions which suppress departure from the reactive configuration(s). We now report that the viscosity dependences of forward and back ET rate constants for the photocycle within the [alpha2(Zn),beta2(Fe3+N3-)] mixed-metal hemoglobin hybrid at pH 7 show the anticipated behavior: kf decreases as viscosity increases, but, in sharp contrast, kb increases strongly.     

On the chaperon mechanism: application to ClO + ClO (+N2) --> ClOOCl (+N2)

The dynamics of the ClO + ClO (+N(2)) radical complex (or chaperon) mechanism is studied by electronic structure methods and quasi-classical trajectory calculations. The geometries and frequencies of the stationary points on the potential energy surface (PES) are optimized at the B3LYP/6-311+G(3df) level of theory, and the energies are refined at the CCSD(T)/6-311+G(3df) (single-point) level of theory. Basis set superposition error (BSSE) corrections are applied to obtain 1.5 kcal mol(-1) for the binding energy of the ClO.N(2) van der Waals (VDW) complex. A model PES is developed and used in quasi-classical trajectory calculations to obtain the capture rate constant and nascent energy distributions of ClOOCl* produced via the chaperon mechanism. A range of VDW binding energies from 1.5 to 9.0 kcal mol(-1) are investigated. The anisotropic PES for the ClO.N(2) complex and a separable anharmonic oscillator approximation are used to estimate the equilibrium constant for formation of the VDW complex. Rate constants, branching ratios to produce ClOOCl, and nascent energy distributions of excited ClOOCl* are discussed with respect to the VDW binding energy and temperature. Interestingly, even for weak VDW binding energies, the N(2) usually carries away enough energy to stabilize the nascent ClOOCl*, although the VDW equilibrium constant is small. For stronger binding energies, the stabilization efficiency is reduced, but the capture rate constant is increased commensurately. The resulting rate constants for forming ClOOCl* from the title reaction are only weakly dependent on the VDW binding energy and temperature. As a result, the relative importance of the chaperon mechanism is mostly dependent on the VDW equilibrium constant. For the calculated ClO.N(2) binding energy of 1.5 kcal mol(-1), the VDW equilibrium constant is small, and the chaperon mechanism is only important at very high pressures.     

Experimental designs for finding stationary points and for computing force constant matrix

We use an experimental design method for computing a local quadratic form which is adequate to find stationary points and to evaluate the force constant matrix on a k-dimensional potential energy hypersurface. Two types of designs are particularly discussed: the composite and the Doehiert's plane. We give examples for illustrating the methodology. Some thermodynamical properties are deduced from our results. Such a method may be used with every type of theoretical calculations (SCF or CI ) and with any program without modification.     

Exact solutions for the entropy production rate of several irreversible processes

We investigate thermal conduction described by Newton's law of cooling and by Fourier's transport equation and chemical reactions based on mass action kinetics where we detail a simple example of a reaction mechanism with one intermediate. In these cases we derive exact expressions for the entropy production rate and its differential. We show that at a stationary state the entropy production rate is an extremum if and only if the stationary state is a state of thermodynamic equilibrium. These results are exact and independent of any expansions of the entropy production rate. In the case of thermal conduction we compare our exact approach with the conventional approach based on the expansion of the entropy production rate near equilibrium. If we expand the entropy production rate in a series and keep terms up to the third order in the deviation variables and then differentiate, we find out that the entropy production rate is not an extremum at a nonequilibrium steady state. If there is a strict proportionality between fluxes and forces, then the entropy production rate is an extremum at the stationary state even if the stationary state is far away from equilibrium.     

Nonequilibrium capillary electrophoresis of equilibrium mixtures--a single experiment reveals equilibrium and kinetic parameters of protein-DNA interactions

We introduce a novel electrophoretic method, nonequilibrium capillary electrophoresis of equilibrium mixtures (NECEEM), and demonstrate its use for studying protein-DNA interactions. The equilibrium mixture of protein and DNA contains three components: free protein, free DNA, and the protein-DNA complex. A short plug of such a mixture is injected into the capillary, and the three components are separated under nonequilibrium conditions. The resulting electropherograms are composed of characteristic peaks and exponential curves. An easy nonnumerical analysis of a single electropherogram reveals two parameters: the equilibrium binding constant and the monomolecular rate constant of complex decay. The bimolecular rate constant of complex formation can then be calculated as the product of the two experimentally determined constants. NECEEM was applied to study the interaction between single-stranded DNA binding protein and a fluorescently labeled 15-mer oligonucleotide. It allowed us to measure for the first time the rate constant of complex decay for this important protein-DNA pair, k-1 = 0.03 s-1. The value of the equilibrium binding constant, Kb = 3.6 x 10-6 M-1, was in good agreement with those measured by other methods. As low as 10-18 mol of the protein was sufficient for the measurements. Thus, the new method is simple, informative, and highly sensitive. Moreover, it can be equally applied to other noncovalent protein-ligand complexes. These features of NECEEM make this method an indispensable tool in studies of macromolecular interactions. They also emphasize the potential role of NECEEM in the development of extremely sensitive protein assays using nucleotide aptamers.     

Accounting for memory effects in calculating the rate constant of a chemical reaction

The first-escape-time method has been used in solving the problem of calculating chemical reaction rate constants on the basis of a model in which a classical particle passes from one potential well to the other, this model describing non-Markovian Brownian motion by means of a generalized Langevin equation with the correlation function of the Ornstein-Uhlenbeck process. Singular perturbation theory was used in the solution. An analytical formula was obtained for the rate constant; calculations using this formula require the solution of a cubic equation. Analysis of the results shows that within a certain interval of parameters, non-Markov character is not manifested; i.e., the rate constant of the process is very little different from the result of the Kramers model. It has been established that the result obtained by this method is the same as that obtained by calculating the rate constant on the basis of the stationary flux method.Translated from Teoreticheskaya i Éksperimental'naya Khimiya, Vol. 24, No. 2, pp. 138–143, March–April, 1988.     

Equilibrium calculations of viscosity and thermal conductivity across a solid-liquid interface using boundary fluctuations

We calculate viscosity and thermal conductivity in systems of Lennard-Jones particles consisting of coexisting solid and liquid with different interface wetting properties using the recently developed equilibrium boundary fluctuation theory. We compare the slip length and equivalent liquid length obtained from these calculations with those obtained from nonequilibrium molecular dynamics. The equilibrium and nonequilibrium calculations of the slip length and the sum of the thermal equivalent lengths are in good agreement. We conclude that for both interfacial properties, the nonequilibrium simulations were probing the linear response. The significant dependence of the intrinsic equivalence length on the interfacial temperature difference used to generate the thermal gradient is explained as a consequence of the different thermodynamic states of the two interfaces.     

Comparison of Accuracy and Convergence Rate between Equilibrium and Nonequilibrium Alchemical Transformations for Calculation of Relative Binding Free Energy

Estimation of protein-ligand binding affinity within chemical accuracy is one of the grand challenges in structure-based rational drug design. With the efforts over three decades, free energy methods based on equilibrium molecular dynamics (MD) simulations have become mature and are nowadays routinely applied in the community of computational chemistry. On the contrary, nonequilibrium MD simulation methods have attracted less attention, despite their underlying rigor in mathematics and potential advantage in efficiency. In this work, the equilibrium and nonequilibrium simulation methods are compared in terms of accuracy and convergence rate in the calculations of relative binding free energies. The proteins studied are T4-lysozyme mutant L99A and COX-2. For each protein, two ligands are studied. The results show that the nonequilibrium simulation method can be competitively as accurate as the equilibrium method, and the former is more efficient than the latter by considering the convergence rate with respect to the cost of wall clock time. In addition, Bennett acceptance ratio, which is a bidirectional post-processing method, converges faster than the unidirectional Jarzynski equality for the nonequilibrium simulations.     

The role of temperature in nucleation processes

Heat and mass transfers are coupled processes, also in nucleation. In principle, a nucleating cluster would have a different temperature compared to the surrounding supersaturated old phase because of the heat release involved with attaching molecules to the cluster. In turn a difference in temperature across the cluster surface is a driving force for the mass transfer to and from the cluster. This coupling of forces in nonisothermal nucleation is described using mesoscopic nonequilibrium thermodynamics, emphasizing measurable heat effects. An expression was obtained for the nonisothermal nucleation rate in a one-component system, in the case where a temperature difference exists between a cluster distribution and the condensed phase. The temperature is chosen as a function of the cluster size only, while the temperature of the condensed phase is held constant by a bath. The generally accepted expression for isothermal stationary nucleation is contained as a limiting case of the nonisothermal stationary nucleation rate. The equations for heat and mass transport were solved for stationary nucleation with the given cluster distribution and with the temperature controlled at the boundaries. A factor was defined for these conditions, determined by the heat conductivity of the surrounding phase and the phase transition enthalpy, to predict the deviation between isothermal and nonisothermal nucleation. For the stationary state described, the factor appears to give small deviations, even for primary nucleation of droplets in vapor, making the nonisothermal rate smaller than the isothermal one. The set of equations may lead to larger and different thermal effects under different boundary conditions, however.     

Extension of solution‐reaction‐surface description to examination of nonequilibrium solvation effect for microsolvated reaction

A microscopic method to examine a nonequilibrium solvation effect is reported. The solution reaction is simplified as a barrier‐crossing reaction within a solution reaction surface that corresponds to a two‐dimensional space determined by solute and solvent reactive coordinates. For this simplification, the motions within the space spanned by nonreactive coordinates are frozen. We derive three rate constant expressions: (1) in the nonadiabatic solvation limit, (2) in the equilibrium solvation limit, and (3) of the transition‐state theory. This method was applied to the examination of the contact‐ion‐pair formation of t‐BuCl in four waters. We found that the nonadiabatic solvation picture overestimates the nonequilibrium solvation effect. © 2000 John Wiley & Sons, Inc. Int J Quant Chem 77: 791–796, 2000     

Reversibility and diffusion in mandelythiamin decarboxylation. Searching dynamical effects in decarboxylation reactions

Decarboxylation of mandelylthiamin in aqueous solution is analyzed by means of quantum mechanics/molecular mechanics simulations including solvent effects. The free energy profile for the decarboxylation reaction was traced, assuming equilibrium solvation, while reaction trajectories allowed us to incorporate nonequilibrium effects due to the solvent degrees of freedom as well as to evaluate the rate of the diffusion process in competition with the backward reaction. Our calculations that reproduce the experimental rate constant show that decarboxylation takes place with a non-negligible free energy barrier for the backward reaction and that diffusion of carbon dioxide is very fast compared to the chemical step. According to these findings catalysts would not act by preventing the backward reaction.     

The inelastic hard dimer gas: a nonspherical model for granular matter

We study a two-dimensional gas of inelastic smooth hard dimers. Since the collisions between dimers are dissipative, being characterized by a coefficient of restitution alpha<1, and no external driving force is present, the energy of the system decreases in time and no stationary state is achieved. However, the resulting nonequilibrium state of the system displays several interesting properties in close analogy with systems of inelastic hard spheres, whose relaxational dynamics has been thoroughly explored. We generalize to inelastic systems a recently method introduced [G. Ciccotti and G. Kalibaeva, J. Stat. Phys. 115, 701 (2004)] to study the dynamics of rigid elastic bodies made up of different spheres held together by rigid bonds. Each dimer consists of two hard disks of diameter d, whose centers are separated by a fixed distance a. By describing the rigid bonds by means of holonomic constraints and deriving the appropriate collision rules between dimers, we reduce the dynamics to a set of equations which can be solved by means of event-driven simulation. After deriving the algorithm we study the decay of the total kinetic energy, and of the ratio between the rotational and the translational kinetic energy of inelastic dimers. We show numerically that the celebrated Haff's homogeneous cooling law t(-2), describing how the kinetic energy of an inelastic hard-sphere system with a constant coefficient of restitution decreases in time, holds even in the case of these nonspherical particles. We fully characterize this homogeneous decay process in terms of appropriate decay constants and confirm numerically the scaling behavior of the velocity distributions.     

Forward flux sampling-type schemes for simulating rare events: efficiency analysis

We analyze the efficiency of several simulation methods which we have recently proposed for calculating rate constants for rare events in stochastic dynamical systems in or out of equilibrium. We derive analytical expressions for the computational cost of using these methods and for the statistical error in the final estimate of the rate constant for a given computational cost. These expressions can be used to determine which method to use for a given problem, to optimize the choice of parameters, and to evaluate the significance of the results obtained. We apply the expressions to the two-dimensional nonequilibrium rare event problem proposed by Maier and Stein [Phys. Rev. E 48, 931 (1993)]. For this problem, our analysis gives accurate quantitative predictions for the computational efficiency of the three methods.