Search

Using the generalized Langevin equation (GLE) is a promising approach to build coarse-grained (CG) models of molecular systems since the GLE model often leads to more accurate thermodynamic and kinetic predictions than Brownian dynamics or Langevin models by including a more sophisticated friction with memory. The GLE approach has been used for CG coordinates such as the center of mass of a group of atoms with pairwise decomposition and for a single CG coordinate. We present a GLE approach when CG coordinates are multiple generalized coordinates, defined, in general, as nonlinear functions of microscopic atomic coordinates. The CG model for multiple generalized coordinates is described by the multidimensional GLE from the Mori-Zwanzig formalism, which includes an exact memory matrix. We first present a method to compute the memory matrix in a multidimensional GLE using trajectories of a full system. Then, in order to reduce the computational cost of computing the multidimensional friction with memory, we introduce a method that maps the GLE to an extended Markovian system. In addition, we study the effect of using a nonconstant mass matrix in the CG model. In particular, we include mass-dependent terms in the mean force. We used the proposed CG model to describe the conformational motion of a solvated alanine dipeptide system, with two dihedral angles as the CG coordinates. We showed that the CG model can accurately reproduce two important kinetic quantities: the velocity autocorrelation function and the distribution of first passage times. [ABSTRACT FROM AUTHOR]

Grate and co-workers at Pacific Northwest National Laboratory recently developed high information content triazine-based sequence-defined polymers that are robust by not having hydrolyzable bonds and can encode structure and functionality by having various side chains. Through molecular dynamics (MD) simulations, the triazine polymers have been shown to form particular sequential stacks, have stable backbone-backbone interactions through hydrogen bonding and π–π interactions, and conserve their cis/trans conformations throughout the simulation. However, we do not know the effects of having different side chains and backbone structures on the entire conformation and whether the cis or trans conformation is more stable for the triazine polymers. For this reason, we investigate the role of non-covalent interactions for different side chains and backbone structures on the conformation and assembly of triazine polymers in MD simulations. Since there is a high energy barrier associated with the cis-trans isomerization, we use replica exchange molecular dynamics (REMD) to sample various conformations of triazine hexamers. To obtain rates and intermediate conformations, we use the recently developed concurrent adaptive sampling (CAS) algorithm for dimers of triazine trimers. We found that the hydrogen bonding ability of the backbone structure is critical for the triazine polymers to self-assemble into nanorod-like structures, rather than that of the side chains, which can help researchers design more robust materials. [ABSTRACT FROM AUTHOR]

Molecular dynamics simulations are useful in obtaining thermodynamic and kinetic properties of bio-molecules, but they are limited by the time scale barrier. That is, we may not obtain properties' efficiently because we need to run microseconds or longer simulations using femtosecond time steps. To overcome this time scale barrier, we can use the weighted ensemble (WE) method, a powerful enhanced sampling method that efficiently samples thermodynamic and kinetic properties. However, the WE method requires an appropriate partitioning of phase space into discrete macrostates, which can be problematic when we have a high-dimensional collective space or when little is known a priori about the molecular system. Hence, we developed a new WE-based method, called the "Concurrent Adaptive Sampling (CAS) algorithm," to tackle these issues. The CAS algorithm is not constrained to use only one or two collective variables, unlike most reaction coordinate-dependent methods. Instead, it can use a large number of collective variables and adaptive macrostates to enhance the sampling in the high-dimensional space. This is especially useful for systems in which we do not know what the right reaction coordinates are, in which case we can use many collective variables to sample conformations and pathways. In addition, a clustering technique based on the committor function is used to accelerate sampling the slowest process in the molecular system. In this paper, we introduce the new method and show results from two-dimensional models and bio-molecules, specifically penta-alanine and a triazine trimer. [ABSTRACT FROM AUTHOR]

Computation of reaction rates and elucidation of reaction mechanisms are two of the main goals of molecular dynamics (MD) and related simulation methods. Since it is time consuming to study reaction mechanisms over long time scales using brute force MD simulations, two ensemble methods, Markov State Models (MSMs) and Weighted Ensemble (WE), have been proposed to accelerate the procedure. Both approaches require clustering of microscopic configurations into networks of "macro-states" for different purposes. MSMs model a discretization of the original dynamics on the macro-states. Accuracy of the model significantly relies on the boundaries of macro-states. On the other hand, WE uses macro-states to formulate a resampling procedure that kills and splits MD simulations for achieving better efficiency of sampling. Comparing to MSMs, accuracy of WE rate predictions is less sensitive to the definition of macro-states. Rigorous numerical experiments using alanine dipeptide and penta-alanine support our analyses. It is shown that MSMs introduce significant biases in the computation of reaction rates, which depend on the boundaries of macro-states, and Accelerated Weighted Ensemble (AWE), a formulation of weighted ensemble that uses the notion of colors to compute fluxes, has reliable flux estimation on varying definitions of macro-states. Our results suggest that whereas MSMs provide a good idea of the metastable sets and visualization of overall dynamics, AWE provides reliable rate estimations requiring less efforts on defining macro-states on the high dimensional conformational space. [ABSTRACT FROM AUTHOR]

In free energy calculations based on thermodynamic integration, it is necessary to compute the derivatives of the free energy as a function of one (scalar case) or several (vector case) order parameters. We derive in a compact way a general formulation for evaluating these derivatives as the average of a mean force acting on the order parameters, which involves first derivatives with respect to both Cartesian coordinates and time. This is in contrast with the previously derived formulas, which require first and second derivatives of the order parameter with respect to Cartesian coordinates. As illustrated in a concrete example, the main advantage of this new formulation is the simplicity of its use, especially for complicated order parameters. It is also straightforward to implement in a molecular dynamics code, as can be seen from the pseudocode given at the end. We further discuss how the approach based on time derivatives can be combined with the adaptive biasing force method, an enhanced sampling technique that rapidly yields uniform sampling of the order parameters, and by doing so greatly improves the efficiency of free energy calculations. Using the backbone dihedral angles [uppercase_phi_synonym] and Ψ in N-acetylalanyl-N′-methylamide as a numerical example, we present a technique to reconstruct the free energy from its derivatives, a calculation that presents some difficulties in the vector case because of the statistical errors affecting the derivatives. [ABSTRACT FROM AUTHOR]

The efficiencies of two recently developed methods for calculating free energy changes along a generalized coordinate in a system are discussed in the context of other, related approaches. One method is based on Jarzynski’s identity [Phys. Rev. Lett. 78, 2690 (1997)]. The second method relies on thermodynamic integration of the average force and is called the adaptive biasing force method [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)]. Both methods are designed such that the system evolves along the chosen coordinate(s) without experiencing free energy barriers and they require calculating the instantaneous, unconstrained force acting on this coordinate using the formula derived by Darve and Pohorille. Efficiencies are analyzed by comparing analytical estimates of statistical errors and by considering two numerical examples—internal rotation of hydrated 1,2-dichloroethane and transfer of fluoromethane across a water-hexane interface. The efficiencies of both methods are approximately equal in the first but not in the second case. During transfer of fluoromethane the system is easily driven away from equilibrium and, therefore, the performance of the method based on Jarzynski’s identity is poor. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A new, general formula that connects the derivatives of the free energy along the selected, generalized coordinates of the system with the instantaneous force acting on these coordinates is derived. The instantaneous force is defined as the force acting on the coordinate of interest so that when it is subtracted from the equations of motion the acceleration along this coordinate is zero. The formula applies to simulations in which the selected coordinates are either unconstrained or constrained to fixed values. It is shown that in the latter case the formula reduces to the expression previously derived by den Otter and Briels [Mol. Phys. 98, 773 (2000)]. If simulations are carried out without constraining the coordinates of interest, the formula leads to a new method for calculating the free energy changes along these coordinates. This method is tested in two examples — rotation around the C–C bond of 1,2-dichloroethane immersed in water and transfer of fluoromethane across the water-hexane interface. The calculated free energies are compared with those obtained by two commonly used methods. One of them relies on determining the probability density function of finding the system at different values of the selected coordinate and the other requires calculating the average force at discrete locations along this coordinate in a series of constrained simulations. The free energies calculated by these three methods are in excellent agreement. The relative advantages of each method are discussed. © 2001 American Institute of Physics. [ABSTRACT FROM AUTHOR]

It is well known that a dilute suspension of spheroids sedimenting under gravity at low Reynolds number is unstable to density fluctuations as a result of hydrodynamic interactions [D. L. Koch and E. S. G. Shaqfeh, J. Fluid Mech. 209, 521 (1989)]. Using a linear stability analysis, it is shown that a vertical density gradient in such a suspension can lead to a wave number selection by damping fluctuations at long wavelengths. A scaling for the most unstable wavelength, or characteristic size of the density fluctuations, is obtained in terms of the background stratification and volume fraction, and is compared to results from numerical simulations in stratified particulate suspensions using methods that we have developed previously. In initially homogeneous suspensions, simulations show a continuous decay of the size of the density fluctuations over time, which we demonstrate can be attributed to the development of stratification inside the suspension. [ABSTRACT FROM AUTHOR]

We use numerical simulations to investigate the dynamics in suspensions of ideally polarizable rods sedimenting under gravity in a vertical electric field. While such suspensions are unstable to concentration fluctuations when no field is applied, we show that the induced-charge electrophoresis that results from the application of the field provides control over the concentration instability by causing particle alignment in the field direction. A phase diagram is obtained for the occurrence of the instability in terms of field strength and volume fraction. In stable suspensions velocity hindrance is shown to occur, and results for the hindered settling function are presented. [ABSTRACT FROM AUTHOR]

Large-scale simulations of non-Brownian rigid fibers sedimenting under gravity at zero Reynolds number have been performed using a fast algorithm. The mathematical formulation follows the previous simulations by Butler and Shaqfeh [“Dynamic simulations of the inhomogeneous sedimentation of rigid fibres,” J. Fluid Mech. 468, 205 (2002)]. The motion of the fibers is described using slender-body theory, and the line distribution of point forces along their lengths is approximated by a Legendre polynomial in which only the total force, torque, and particle stresslet are retained. Periodic boundary conditions are used to simulate an infinite suspension, and both far-field hydrodynamic interactions and short-range lubrication forces are considered in all simulations. The calculation of the hydrodynamic interactions, which is typically the bottleneck for large systems with periodic boundary conditions, is accelerated using a smooth particle-mesh Ewald (SPME) algorithm previously used in molecular dynamics simulations. In SPME the slowly decaying Green’s function is split into two fast-converging sums: the first involves the distribution of point forces and accounts for the singular short-range part of the interactions, while the second is expressed in terms of the Fourier transform of the force distribution and accounts for the smooth and long-range part. Because of its smoothness, the second sum can be computed efficiently on an underlying grid using the fast Fourier transform algorithm, resulting in a significant speed-up of the calculations. Systems of up to 512 fibers were simulated on a single-processor workstation, providing a different insight into the formation, structure, and dynamics of the inhomogeneities that occur in sedimenting fiber suspensions. [ABSTRACT FROM AUTHOR]

Memory effects are often introduced during coarse-graining of a complex dynamical system. In particular, a generalized Langevin equation (GLE) for the coarse-grained (CG) system arises in the context of Mori-Zwanzig formalism. Upon a pairwise decomposition, GLE can be reformulated into its pairwise version, i.e., non-Markovian dissipative particle dynamics (DPD). GLE models the dynamics of a single coarse particle, while DPD considers the dynamics of many interacting CG particles, with both CG systems governed by non-Markovian interactions. We compare two different methods for the practical implementation of the non-Markovian interactions in GLE and DPD systems. More specifically, a direct evaluation of the non-Markovian (NM) terms is performed in LE-NM and DPD-NM models, which requires the storage of historical information that significantly increases computational complexity. Alternatively, we use a few auxiliary variables in LE-AUX and DPD-AUX models to replace the non-Markovian dynamics with a Markovian dynamics in a higher dimensional space, leading to a much reduced memory footprint and computational cost. In our numerical benchmarks, the GLE and non-Markovian DPD models are constructed from molecular dynamics (MD) simulations of star-polymer melts. Results show that a Markovian dynamics with auxiliary variables successfully generates equivalent non-Markovian dynamics consistent with the reference MD system, while maintaining a tractable computational cost. Also, transient subdiffusion of the star-polymers observed in the MD system can be reproduced by the coarse-grained models. The non-interacting particle models, LE-NM/AUX, are computationally much cheaper than the interacting particle models, DPD-NM/AUX. However, the pairwise models with momentum conservation are more appropriate for correctly reproducing the long-time hydrodynamics characterised by an algebraic decay in the velocity autocorrelation function.