Grcarma: A fully automated task‐oriented interface for the analysis of molecular dynamics trajectories

We report the availability of grcarma, a program encoding for a fully automated set of tasks aiming to simplify the analysis of molecular dynamics trajectories of biological macromolecules. It is a cross‐platform, Perl/Tk‐based front‐end to the program carma and is designed to facilitate the needs of the novice as well as those of the expert user, while at the same time maintaining a user‐friendly and intuitive design. Particular emphasis was given to the automation of several tedious tasks, such as extraction of clusters of structures based on dihedral and Cartesian principal component analysis, secondary structure analysis, calculation and display of root‐mean‐square deviation (RMSD) matrices, calculation of entropy, calculation and analysis of variance–covariance matrices, calculation of the fraction of native contacts, etc. The program is free‐open source software available immediately for download. © 2013 Wiley Periodicals, Inc.     

Essential dynamics/factor analysis for the interpretation of molecular dynamics trajectories

Subject of this work is the analysis of molecular dynamics (MD) trajectories of neurophysins I (NPI) and II (NPII) and their complexes with the neurophyseal nonapeptide hormones oxytocin (OT) and vasopresssin (VP), respectively, simulated in water. NPs serve in the neurosecretory granules as carrier proteins for the hormones before their release to the blood. The starting data consisted of two pairs of different trajectories for each of the (NPII/VP)2 and (NPI/OT)2 heterotetramers and two more trajectories for the NPII2 and NPI2 homodimers (six trajectories in total). Using essential dynamics which, to our judgement, is equivalent to factor analysis, we found that only about 10 degrees of freedom per trajectory are necessary and sufficient to describe in full the motions relevant for the function of the protein. This is consistent with these motions to explain about 90% of the total variance of the system. These principal degrees of freedom represent slow anharmonic motional modes, clearly pointing at distinguished mobility of the atoms involved in the protein's functionality.     

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     

Mechanisms of double ionization in strong laser field from simulation with Coupled Coherent States: Beyond reduced dimensionality models

The recently developed Coupled Coherent States method for solution of the time-dependent Schrödinger equation is used to simulate the strong laser field double ionization of a helium atom in six spatial dimensions. The calculated double ionization yield, as a function of laser intensity, reproduces the “shoulder”, attributable to electron recollisions. The method, which employs a Coherent State representation guided by classical trajectories also provides physical insight into the ionization mechanism. The appearance of the ‘shoulder’, is seen by analysis of the guiding trajectories to arise predominantly from sequential excitation of one electron into a highly excited orbit, followed by an electron–electron collision leading to correlated escape of the two electrons. A second little known mechanism involves recollisions resulting in ionization of the secondary electron, leaving the first in a highly excited state, from which it is later ionized by the laser field.     

New insights into the meaning and usefulness of principal component analysis of concatenated trajectories

A comparison between different conformations of a given protein, relating both structure and dynamics, can be performed in terms of combined principal component analysis (combined‐PCA). To that end, a trajectory is obtained by concatenating molecular dynamics trajectories of the individual conformations under comparison. Then, the principal components are calculated by diagonalizing the correlation matrix of the concatenated trajectory. Since the introduction of this approach in 1995 it has had a large number of applications. However, the interpretation of the eigenvectors and eigenvalues so obtained is based on intuitive foundations, because analytical expressions relating the concatenated correlation matrix with those of the individual trajectories under consideration have not been provided yet. In this article, we present such expressions for the cases of two, three, and an arbitrary number of concatenated trajectories. The formulas are simple and show what is to be expected and what is not to be expected from a combined‐PCA. Their correctness and usefulness is demonstrated by discussing some representative examples. The results can be summarized in a simple sentence: the correlation matrix of a concatenated trajectory is given by the average of the individual correlation matrices plus the correlation matrix of the individual averages. From this it follows that the combined‐PCA of trajectories belonging to different free energy basins provides information that could also be obtained by alternative and more straightforward means. © 2014 Wiley Periodicals, Inc.     

An application of error reduction and harmonic inversion schemes to the semiclassical calculation of molecular vibrational energy levels

A singular value decomposition based harmonic inversion signal processing scheme is applied to the semiclassical initial value representation (IVR) calculation of molecular vibrational states. Relative to usual IVR procedure of Fourier analysis of a signal made from the Monte Carlo evaluation of the phase space integral in which many trajectories are needed, the new procedure obtains acceptable results with many fewer trajectories. Calculations are carried out for vibrational energy levels of H2O to illustrate the overall procedure.     

Active-site mobility inhibits reductive dehalogenation of 1,1,1-trichloroethane by cytochrome P450cam

Summary Recent studies by Wackett and co-workers have shown that cytochrome P450cam is capable of reductively dehalogenating hexachloroethane at a significant rate, but that no appreciable dehalogenation of 1,1,1-trichloroethane is observed. A growing body of evidence indicates that differences in intrinsic reactivity can not completely explain this observation. We therefore explored the possible role of differences in preferred binding orientation and in active-site mobility. A detailed analysis of molecular dynamics trajectories with each of these substrates bound at the active site of P450cam is presented. While the dynamics and overall time-average structure calculated for the protein are similar in the two trajectories, the two substrates behave quite differently. The smaller substrate, 1,1,1-trichloroethane, is significantly more mobile than hexachloroethane and has a preferred orientation in which the substituted carbon is generally far from the heme iron. In contrast, for hexachloroethane, one of the chlorine atoms is nearly always in van der Waals contact with the heme iron, which should favor the initial electron transfer step.     

Energy exchange network of inter‐residue interactions within a thermally fluctuating protein molecule: A computational study

Protein function is regulated not only by the structure but also by physical dynamics and thermal fluctuations. We have developed the computer program, CURrent calculation for proteins (CURP), for the flow analysis of physical quantities within thermally fluctuating protein media. The CURP program was used to calculate the energy flow within the third PDZ domain of the neuronal protein PSD‐95, and the results were used to illustrate the energy exchange network of inter‐residue interactions based on atomistic molecular dynamics simulations. The removal of the α3 helix is known to decrease ligand affinity by 21‐fold without changing the overall protein structure; nevertheless, we demonstrated that the helix constitutes an essential part of the network graph. © 2015 Wiley Periodicals, Inc.     

A_2B模型分子经典轨迹的辛算法计算

采用辛算法计算了Ａ２Ｂ模型分子的经典轨迹并与传统Ｒｕｎｇｅ－Ｋｕｔｔａ　（Ｒ－Ｋ）算法进行了比较．结果表明，在微观反应动力学研究所应考虑的时间范围内，辛算法的结果与理论分析一致，Ｒ－Ｋ法的结果则面目全非．因此，用辛算法替代传统数值方法有可能克服目前经典轨迹计算存在的困难，从根本上改进微观反应动力学研究的经典轨迹方法．     

Internal dynamics of flexible molecules: Cyclohexane

Summary The international motions of a single cyclohexane molecule are studied by molecular dynamics calculations. Classical trajectories are calculated by integrating Newton's equation of motion. The potential functions used are essentially the same as in Allinger's MM2 program which is widely applied for calculations on conformational energies of organic molecules.Geometries and relative energies are reported for all stable low energy conformers and some transition states of cyclohexane. Vibrational frequencies of classical oscillations of individual bonds — computed for ethane as reference system — are close to the experimental values.Two trajectories of the molecular dynamics of cyclohexane were simulated. In the first we were able to follow the process of ring inversion: chair twisted forms inverted chair. The reaction path is analysed in detail and compared with static approaches. The second trajectory shows the correlated reorientation of the possible twisted forms alone. This process is known as pseudorotation.Dedicated to Prof. Dr. Karl Schlögl     

CHARMM-GUI: a web-based graphical user interface for CHARMM

CHARMM is an academic research program used widely for macromolecular mechanics and dynamics with versatile analysis and manipulation tools of atomic coordinates and dynamics trajectories. CHARMM-GUI, http://www.charmm-gui.org, has been developed to provide a web-based graphical user interface to generate various input files and molecular systems to facilitate and standardize the usage of common and advanced simulation techniques in CHARMM. The web environment provides an ideal platform to build and validate a molecular model system in an interactive fashion such that, if a problem is found through visual inspection, one can go back to the previous setup and regenerate the whole system again. In this article, we describe the currently available functional modules of CHARMM-GUI Input Generator that form a basis for the advanced simulation techniques. Future directions of the CHARMM-GUI development project are also discussed briefly together with other features in the CHARMM-GUI website, such as Archive and Movie Gallery.     

Computing ensembles of transitions from stable states: Dynamic importance sampling

There is an increasing dataset of solved biomolecular structures in more than one conformation and increasing evidence that large-scale conformational change is critical for biomolecular function. In this article, we present our implementation of a dynamic importance sampling (DIMS) algorithm that is directed toward improving our understanding of important intermediate states between experimentally defined starting and ending points. This complements traditional molecular dynamics methods where most of the sampling time is spent in the stable free energy wells defined by these initial and final points. As such, the algorithm creates a candidate set of transitions that provide insights for the much slower and probably most important, functionally relevant degrees of freedom. The method is implemented in the program CHARMM and is tested on six systems of growing size and complexity. These systems, the folding of Protein A and of Protein G, the conformational changes in the calcium sensor S100A6, the glucose-galactose-binding protein, maltodextrin, and lactoferrin, are also compared against other approaches that have been suggested in the literature. The results suggest good sampling on a diverse set of intermediates for all six systems with an ability to control the bias and thus to sample distributions of trajectories for the analysis of intermediate states.     

Interference and quantization in semiclassical response functions

Application of the Herman-Kluk semiclassical propagator to the calculation of spectroscopic response functions for anharmonic oscillators has demonstrated the quantitative accuracy of these approximate dynamics. In this approach, spectroscopic response functions are expressed as multiple phase-space integrals over pairs of classical trajectories and their associated stability matrices. Here we analyze the Herman-Kluk semiclassical approximation to a linear response function and determine the origin of the capacity of this method to reproduce quantum effects in a response function from classical dynamical information. Our analysis identifies those classical trajectories that contribute most significantly to the response function on different time scales. This finding motivates a procedure for computing the linear response function in which the interference between pairs of classical trajectories is treated approximately, resulting in an integral over a single average trajectory, as in a purely classical calculation.     

Exploring potential energy surfaces for chemical reactions: an overview of some practical methods

Potential energy surfaces form a central concept in the application of electronic structure methods to the study of molecular structures, properties, and reactivities. Recent advances in tools for exploring potential energy surfaces are surveyed. Methods for geometry optimization of equilibrium structures, searching for transition states, following reaction paths and ab initio molecular dynamics are discussed. For geometry optimization, topics include methods for large molecules, QM/MM calculations, and simultaneous optimization of the wave function and the geometry. Path optimization methods and dynamics based techniques for transition state searching and reaction path following are outlined. Developments in the calculation of ab initio classical trajectories in the Born-Oppenheimer and Car-Parrinello approaches are described.     

On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: A review

Methods to search for low-energy conformations, to generate a Boltzmann-weighted ensemble of configurations, or to generate classical-dynamical trajectories for molecular systems in the condensed liquid phase are briefly reviewed with an eye to application to biomolecular systems. After having chosen the degrees of freedom and method to generate molecular configurations, the efficiency of the search or sampling can be enhanced in various ways: (i) efficient calculation of the energy function and forces, (ii) application of a plethora of search enhancement techniques, (iii) use of a biasing potential energy term, and (iv) guiding the sampling using a reaction or transition pathway. The overview of the available methods should help the reader to choose the combination that is most suitable for the biomolecular system, degrees of freedom, interaction function, and molecular or thermodynamic properties of interest.     

Dcdftbmd: Divide-and-Conquer Density Functional Tight-Binding Program for Huge-System Quantum Mechanical Molecular Dynamics Simulations

Dcdftbmd is a Fortran 90/95 program that enables efficient quantum mechanical molecular dynamics (MD) simulations using divide-and-conquer density functional tight-binding (DC-DFTB) method. Based on the remarkable performance of previous massively parallel DC-DFTB energy and gradient calculations for huge systems, the code has been specialized to MD simulations. Recent implementations and modifications including DFTB extensions, improved computational speed in the DC-DFTB computational steps, algorithms for efficient initial guess charge prediction, and free energy calculations via metadynamics technique have enhanced the capability to obtain atomistic insights in novel applications to nanomaterials and biomolecules. The energy, structure, and other molecular properties are also accessible through the single-point calculation, geometry optimization, and vibrational frequency analysis. The available functionalities are outlined together with efficiency tests and simulation examples. © 2019 Wiley Periodicals, Inc.     

E−Z Isomerization in Guanidine: Second-order Saddle Dynamics,Non-statisticality,and Time-frequency Analysis

Our recent work on the E−Z isomerization reaction of guanidine using ab initio chemical dynamics simulations [Rashmi et al., Regul. Chaotic Dyn. 2021 , 26, 119] emphasized the role of second-order saddle ( SOS ) in the isomerization reaction; however, we could not unequivocally establish the non-statistical nature of the dynamics followed in the reaction. In the present study, we performed thousands of on-the-fly trajectories using forces computed at the MNDO level to investigate the influence of second-order saddle in the E−Z isomerization reaction of guanidine and the role of intramolecular vibrational energy redistribution (IVR) on the reaction dynamics. The simulations reveal that while majority of the trajectories follow the traditional transition state pathways, 15 % of the trajectories follow the SOS path. The dynamics was found to be highly non-statistical with the survival probabilities of the reactants showing large deviations from those obtained within the RRKM assumptions. In addition, a detailed analysis of the dynamics using time-dependent frequencies and the frequency ratio spaces reveal the existence of multiple resonance junctions that indicate the existence of regular dynamics and long-lived quasi-periodic trajectories in the phase space associated with non-RRKM behavior.