Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. II. Simulation results

The results of a series of constant pressure and temperature molecular-dynamics (MD) simulation studies based on the rigorous shell particle formulation of the isothermal-isobaric (NpT) ensemble are presented. These MD simulations validate the newly proposed constant pressure equations of motion in which a "shell" particle is used to define uniquely the volume of the system [M. J. Uline and D. S. Corti, J. Chem. Phys. (to be published), preceding paper]. Ensemble averages obtained with the new MD NpT algorithm match the ensemble averages obtained using the previously derived shell particle Monte Carlo NpT method [D. S. Corti, Mol. Phys. 100, 1887 (2002)]. In addition, we also verify that the Hoover NpT MD algorithm [W. G. Hoover, Phys. Rev. A 31, 1695 (1985); 34, 2499 (1986)] generates the correct ensemble averages, though only when periodic boundary conditions are employed. The extension of the shell particle MD algorithm to multicomponent systems is also discussed, in which we show for equilibrium properties that the identity of the shell particle is completely arbitrary when periodic boundary conditions are applied. Self-diffusion coefficients determined with the shell particle equations of motion are also identical to those obtained in other ensembles. Finally, since the mass of the shell particle is known, the system itself, and not a piston of arbitrary mass, controls the time scales for internal pressure and volume fluctuations. We therefore consider the effects of the shell particle on the dynamics of the system. Overall, the shell particle MD algorithm is an effective simulation method for studying systems exposed to a constant external pressure and may provide an advantage over other existing constant pressure approaches when developing nonequilibrium MD methods.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. III. Discontinuous potentials

Based on the approach of Gruhn and Monson [Phys. Rev. E 63, 061106 (2001)], we present a new method for deriving the collisions dynamics for particles that interact via discontinuous potentials. By invoking the conservation of the extended Hamiltonian, we generate molecular dynamics (MD) algorithms for simulating the hard-sphere and square-well fluids within the isothermal-isobaric (NpT) ensemble. Consistent with the recent rigorous reformulation of the NpT ensemble partition function, the equations of motion impose a constant external pressure via the introduction of a shell particle of known mass [M. J. Uline and D. S. Corti, J. Chem. Phys. 123, 164101 (2005); 123, 164102 (2005)], which serves to define uniquely the volume of the system. The particles are also connected to a temperature reservoir through the use of a chain of Nose-Hoover thermostats, the properties of which are not affected by a hard-sphere or square-well collision. By using the Liouville operator formalism and the Trotter expansion theorem to integrate the equations of motion, the update of the thermostat variables can be decoupled from the update of the positions of the particles and the momentum changes upon a collision. Hence, once the appropriate collision dynamics for the isobaric-isenthalpic (NpH) equations of motion is known, the adaptation of the algorithm to the NpT ensemble is straightforward. Results of MD simulations for the pure component square-well fluid are presented and serve to validate our algorithm. Finally, since the mass of the shell particle is known, the system itself, and not a piston of arbitrary mass, controls the time scales for internal pressure and volume fluctuations. We therefore consider the influence of the shell particle algorithm on the dynamics of the square-well fluid.     

Density functional theory of fluids in the isothermal-isobaric ensemble

We present a density functional theory for inhomogeneous fluids at constant external pressure. The theory is formulated for a volume-dependent density, n(r,V), defined as the conjugate variable of a generalized external potential, nu(r,V), that conveys the information on the pressure. An exact expression for the isothermal-isobaric free-energy density functional is obtained in terms of the corresponding canonical ensemble functional. As an application we consider a hard-sphere system in a spherical pore with fluctuating radius. In general we obtain very good agreement with simulation. However, in some situations a peak develops in the center of the cavity and the agreement between theory and simulation becomes worse. This happens for systems where the number of particles is close to the magic numbers N=13, 55, and 147.     

Non-Hamiltonian molecular dynamics implementation of the Gibbs ensemble method. I. Algorithm

In this paper (paper I) and the following paper (paper II) [C. Bratschi, H. Huber, and D. J. Searles, J. Chem. Phys. 126, 164105 (2007)], a new molecular dynamics algorithm implementing the Gibbs ensemble will be presented and then on applied to the liquid-vapor coexistence curve for two ab initio CO2 potentials. In paper I, the Gibbs ensemble molecular dynamics algorithm using non-Hamiltonian molecular dynamics techniques is introduced. It is shown that states of the correct probability density function are sampled and the correct exchange probability is generated. The extended system Nose-Hoover formalism is used to generate a constant temperature ensemble with equal pressures in the subsystems, combined with single particle transfers between the subsystems, over several time steps, to get equal chemical potentials.     

Recovering the Crooks equation for dynamical systems in the isothermal-isobaric ensemble: a strategy based on the equations of motion

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)] relates the work done on a system during a nonequilibrium transformation to the free energy difference between the final and the initial state of the transformation. Recently, the authors have derived the Crooks equation for systems in the canonical ensemble thermostatted by the Nose-Hoover or Nose-Hoover chain method [P. Procacci et al., J. Chem. Phys. 125, 164101 (2006)]. That proof is essentially based on the fluctuation theorem by Evans and Searles [Adv. Phys. 51, 1529 (2002)] and on the equations of motion. Following an analogous approach, the authors derive here the Crooks equation in the context of molecular dynamics simulations of systems in the isothermal-isobaric (NPT) ensemble, whose dynamics is regulated by the Martyna-Tobias-Klein algorithm [J. Chem. Phys. 101, 4177 (1994)]. Their present derivation of the Crooks equation correlates to the demonstration of the Jarzynski identity for NPT systems recently proposed by Cuendet [J. Chem. Phys. 125, 144109 (2006)].     

Molecular simulation of dendrimers and their mixtures under shear: comparison of isothermal-isobaric (NpT) and isothermal-isochoric (NVT) ensemble systems

Flow properties of dendrimers are studied with the aid of nonequilibrium molecular dynamics techniques. Simulations are performed in the NpT ensemble using the NpT-SLLOD algorithm [P. J. Davis and D. J. Evans, J. Chem. Phys. 100, 541 (1994)] and are compared to the results from simulations performed in the NVT ensemble reported earlier [J. T. Bosko, B. D. Todd, and R. J. Sadus, Chem. Phys. 121, 12050 (2004)]. Shear thickening observed at high strain rates vanishes in systems kept under constant pressure. Also the exponents in the power-law dependencies of the viscosity and the normal stress coefficients change. The variations are significant only at high strain rates and do not affect largely microscopic properties such as shape, alignment, or rotation of molecules. The NpT-SLLOD algorithm has been applied to study various systems including dendrimers in solution and their blends with linear chain molecules of the same molecular mass, and some results for these systems are presented.     

Theory for the atomic shell structure of the cluster magnetic moment and magnetoresistance of a cluster ensemble

We present a simple theory for the cluster size dependence of the average cluster magnetic moment of transition metal clusters. Assuming a local environmental dependence of the atomic magnetic moments, the cluster magnetization exhibits a magnetic shell structure, reflecting the atomic structure of the cluster. Thus, the observed oscillations of the average cluster magnet moment may serve as a fingerprint of the cluster geometry. We also discuss the giant magnetoresistance (GMR) exhibited by an ensemble of magnetic clusters embedded in a metallic matrix. It is shown that the magnetic anisotropy affects strongly the magnetization of the cluster ensemble under certain conditions. Since the GMR depends on the cluster ensemble magnetization, it can be used to determine the cluster magnetic anisotropy energy.     

Relaxation dynamics of a linear molecule in a random static medium: a scaling analysis

We present extensive molecular dynamics simulations of the motion of a single linear rigid molecule in a two-dimensional random array of fixed overlapping disklike obstacles. The diffusion constants for the center of mass translation, D(CM), and for rotation, D(R), are calculated for a wide range of the molecular length, L, and the density of obstacles, rho. The obtained results follow a master curve Drho(micro) approximately (L(2)rho)(-nu) with an exponent micro=-3/4 and 1/4 for D(R) and D(CM), respectively, that can be deduced from simple scaling and kinematic arguments. The nontrivial positive exponent nu shows an abrupt crossover at L(2)rho=zeta(1). For D(CM) we find a second crossover at L(2)rho=zeta(2). The values of zeta(1) and zeta(2) correspond to the average minor and major axis of the elliptic holes that characterize the random configuration of the obstacles. A violation of the Stokes-Einstein-Debye relation is observed for L(2)rho>zeta(1), in analogy with the phenomenon of enhanced translational diffusion observed in supercooled liquids close to the glass transition temperature.     

Molecular dynamics of apo-adenylate kinase: a principal component analysis

Adenylate kinase from E. coli (AKE) is studied with molecular dynamics. AKE undergoes large-scale motions of its Lid and AMP-binding domains when its open form closes over its substrates, AMP and Mg2+-ATP. The third domain, the Core, is relatively stable during closing. The resulting trajectory is analyzed with a principal component analysis method that decomposes the atom motions into modes ordered by their decreasing contributions to the total protein fluctuation. Simulations at 303 K (normal T) and 500 K (high T) reveal that at both temperatures the first three modes account for 70% of the total fluctuation. The residues that contribute the most to these three modes are concentrated in the Lid and AMP-binding domains. Analysis of the normal T modes indicates that the Lid and AMP-binding domains sample a broad distribution of conformations indicating that AKE is designed to provide its substrates with a large set of conformations. The high T results show that the Lid initially closes toward the Core. Subsequently, the Lid rotates to a new stable conformation that is different from what is observed in the substrate-bound AKE. These results are discussed in the context of experimental data that indicate that adenylate kinases do sample more than one conformational state in solution and that each of these conformational states undergoes substantial fluctuations. A pair of residues is suggested for labeling that would be useful for monitoring distance fluctuations by energy transfer experiments.     

Non-Hamiltonian molecular dynamics implementation of the Gibbs ensemble method. II. Molecular liquid-vapor results for carbon dioxide

The Gibbs ensemble molecular dynamics algorithm introduced in the preceding paper (paper I) [C. Bratschi and H. Huber, J. Chem. Phys. v126, 164104 (2007)] is applied to two recently published CO2 ab initio pair potentials, the Bock-Bich-Vogel and symmetry-adapted perturbation theory site-site potentials. The critical properties of these potentials are calculated for the first time. Critical values and points in the single and two-phase zones are compared with Monte Carlo results to demonstrate the accuracy of the molecular dynamics algorithm, and are compared with experiment to test the accuracy of the potentials. Pressure calculations in the liquid, gas, and supercritical states are carried out and are used to explain potential-related effects and systematic discrepancies. The best ab initio potential yields results in good agreement with experiment.     

Molecular simulation of protein dynamics in nanopores. I. Stability and folding

Discontinuous molecular dynamics simulations, together with the protein intermediate resolution model, an intermediate-resolution model of proteins, are used to carry out several microsecond-long simulations and study folding transition and stability of alpha-de novo-designed proteins in slit nanopores. Both attractive and repulsive interaction potentials between the proteins and the pore walls are considered. Near the folding temperature T(f) and in the presence of the attractive potential, the proteins undergo a repeating sequence of folding/partially folding/unfolding transitions, with T(f) decreasing with decreasing pore sizes. The unfolded states may even be completely adsorbed on the pore's walls with a negative potential energy. In such pores the energetic effects dominate the entropic effects. As a result, the unfolded state is stabilized, with a folding temperature T(f) which is lower than its value in the bulk and that, compared with the bulk, the folding rate decreases. The opposite is true in the presence of a repulsive interaction potential between the proteins and the walls. Moreover, for short proteins in very tight pores with attractive walls, there exists an unfolded state with only one alpha-helical hydrogen bond and an energy nearly equal to that of the folded state. The proteins have, however, high entropies, implying that they cannot fold onto their native structure, whereas in the presence of repulsive walls the proteins do attain their native structure. There is a pronounced asymmetry between the two termini of the protein with respect to their interaction with the pore walls. The effect of a variety of factors, including the pore size and the proteins' length, as well as the temperature, is studied in detail.     

Molecular dynamics simulations of a liquid crystalline molecule in the smectic A and E phases and in vacuum

Molecular dynamics simulations are performed in this work at 393 and 323 K for a mesogenic molecule ( R )-1-methylheptyl 4[4-(2-allyloxyethoxy)biphenyl-4-carbonyloxy]benzoate in the simulated smectics A and E, respectively, and in a vacuum at 300 K, for a period of 1.0ns. The trajectories obtained from molecular dynamics simulations allow us to investigate the dynamical behaviour of this mesogenic molecule in the simulated smectic phases. This dynamical behaviour of a single molecule is presented using the distributions of dihedral angles and rotational diffusion around the C-axis defined by the simulated cells. Simulation results indicate that, except for the bonds near the end of the spacer segment, the dihedral angles all exhibit a single Gaussian-like distribution in the smectic A and E phases. Fluctuations of a dihedral angle about its mean value are more restricted in the smectics A and E than in those simulated in a vacuum. The average value of the fluctuations of the dihedral angles at the bonds in the spacer is found to be about 2 fold larger than that of fluctuations in the tail of the same molecule in the smectic A and E phases. In the smectic A phase, the distribution of orientations of a molecule about its long axis in a 36 molecule cell in which the outer molecules are fixed is found to have three distinct peaks. This result shows that the orientational fluctuations of single molecules are limited by confinement due to neighbouring molecules, i.e. that the layers have short-range structural correlations. The orientational distributions show larger fluctuations at the ends of the molecules.     

Molecular dynamics simulation study of nonconcatenated ring polymers in a melt. I. Statics

Molecular dynamics simulations were conducted to investigate the structural properties of melts of nonconcatenated ring polymers and compared to melts of linear polymers. The longest rings were composed of N = 1600 monomers per chain which corresponds to roughly 57 entanglement lengths for comparable linear polymers. For the rings, the radius of gyration squared, [linear span]R(g)(2)[linear span], was found to scale as N(4/5) for an intermediate regime and N(2/3) for the larger rings indicating an overall conformation of a crumpled globule. However, almost all beads of the rings are "surface beads" interacting with beads of other rings, a result also in agreement with a primitive path analysis performed in the next paper [J. D. Halverson, W. Lee, G. S. Grest, A. Y. Grosberg, and K. Kremer, J. Chem. Phys. 134, 204905 (2011)]. Details of the internal conformational properties of the ring and linear polymers as well as their packing are analyzed and compared to current theoretical models.     

Single molecule kinetics. I. Theoretical analysis of indicators

Single molecule experiments reveal intriguing phenomenon in chemical and biological systems. Several indicators of complex dynamics, including "intensity" correlations, "event" correlations, and characteristic functions have been proposed, but extraction of information from these indicators can be difficult since these indicators only observe certain characteristics of the system. Generally, for systems that follow Poisson kinetics, all of these indicators contain similar information about the relaxation times of the system and the connections between different relaxation times, but the information is convoluted in different ways so the strength of various indicators is system specific. The paper discusses the theoretical implications and information content of various data analysis methods for single molecule experiments and demonstrates the relationships between indicators. Under certain conditions, common indicators contain all available information about systems with Poisson kinetics between degenerate states, but extraction of this information is generally not numerically feasible. The paper also discusses practical issues associated with these analyses, which motivates a numerical study based on Bayes' formula in the companion paper [J. Witkoskie and J. S. Cao, J. Chem. Phys. 121, 6373 (2004), following paper].     

Short-time dynamics through conical intersections in macrosystems. I. Theory: effective-mode formulation

The short-time dynamics through a conical intersection of a macrosystem comprising a large number of nuclear degrees of freedom (modes) is investigated. The macrosystem is decomposed into a "system" part carrying a limited number of modes, and an "environment" part. An orthogonal transformation in the environment's space is introduced, as a result of which a subset of three effective modes can be identified which couple directly to the electronic subsystem. Together with the system's modes, these govern the short-time dynamics of the overall macrosystem. The remaining environmental modes couple, in turn, to the effective modes and become relevant at longer times. In this paper, we present the derivation of the effective Hamiltonian, first introduced by Cederbaum et al. [Phys. Rev. Lett. 94, 113003 (2005)], and analyze its properties in some detail. Several special cases and topological aspects are discussed.     

Spontaneous and induced dynamic fluctuations in glass formers. I. General results and dependence on ensemble and dynamics

We study theoretically and numerically a family of multipoint dynamic susceptibilities that quantify the strength and characteristic length scales of dynamic heterogeneities in glass-forming materials. We use general theoretical arguments (fluctuation-dissipation relations and symmetries of relevant dynamical field theories) to relate the sensitivity of averaged two-time correlators to temperature and density to spontaneous fluctuations of the local dynamics. Our theoretical results are then compared to molecular dynamics simulations of the Newtonian, Brownian, and Monte Carlo dynamics of two representative glass-forming liquids, a fragile binary Lennard-Jones mixture, and a model for the strong glass-former silica. We justify in detail the claim made by Berthier et al. [Science 310, 1797 (2005)] that the temperature dependence of correlation functions allows one to extract useful information on dynamic length scales in glassy systems. We also discuss some subtle issues associated with the choice of microscopic dynamics and of statistical ensemble through conserved quantities, which are found to play an important role in determining dynamic correlations.     

Molecular diffusion in liquid crystals and chiral discrimination. I. Theory

The possibility of using cholesteric phases for discriminating enantiomers of a chiral solute on the basis of their different transport properties, motivates the investigation of the translational diffusion by taking fully into account the roto-translational coupling. In this article a detailed theoretical analysis is presented for the transport properties evaluated according to the asymptotic limit of the mean-squared displacement. A general relation is derived for the transport coefficients, having as main ingredients the mean-field potential due to the mesophase, and the diffusion tensor with its purely translational and rotational components, and with the blocks describing the roto-translational coupling. The application of the theory to nematic phases shows that the roto-translational coupling generates a dynamical contribution reducing the transport coefficients evaluated by taking into account only the translational diffusion components in the center of diffusion. The theory is also specialized to a cholesteric phase with a given helical pitch for the director arrangement, in a form which is suitable for calculations of model systems of chiral solutes to be presented in a forthcoming paper.     

The dynamics of a bubble ensemble in a cavitation field

A system of equations was obtained to describe the dynamics of bubbles in a cavitation cloud taking into account the interaction of pulsating bubbles involved in translational motion. The kinetics of cavitation bubble concentration changes, changes in the compressibility of the liquid, and phase transitions within a cavitation bubble and in the neighboring volume of the liquid were taken into account. The role played by bubble deformation in a cavitation cloud was considered. The Bernoulli pressure effect was shown to be negligible. The interaction of cavitation bubbles was a substantial factor that strongly influenced the dynamics of bubbles. It was suggested that there was at least one more mechanism that reduced sonoluminescence intensity from the multiple-bubble cavitation field, namely, a fairly high efficiency of sonoluminescence quenching could additionally be related to the arrival of a cumulative liquid stream at the central cavitation bubble region, where the concentration of active species was high. The dynamics of bubbles in the cavitation field is not only related to the expansion and compression of cavitation bubbles in the acoustic field, but also governed to a great extent by their interaction, translational motion, deformation, and the influence of cumulative streams penetrating the bubbles.