Coordinate reduction for exploring chemical reaction paths

The potential energy surface for the reaction of a typical molecular system composed of N atoms is defined uniquely by 3N-6 coordinates. These coordinates can be defined by the Cartesian coordinates of the atomic centers (minus overall translation and rotation), or a set of internally defined coordinates such as bond stretches, angle bends, and torsions. By applying principal component analysis to the geometries along a reaction path, a reduced set of coordinates, d ≪ 3N-6, can be obtained. This reduced set of coordinates can reproduce the changes in geometry along the reaction path with chemical accuracy and may help improve the efficiency of reaction path optimization algorithms.     

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual‐system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free‐energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.     

Model hysteresis dimer molecule II: deductions from probability profile due to system coordinates

The hysteresis dimer reaction of the first sequel is applied to test the Gibbs density-in-phase hypothesis for a canonical distribution at equilibrium. The probability distribution of variously defined internal and external variables is probed using the algorithms described, in particular the novel probing of the energy states of a labeled particle where it is found that there is compliance with the Gibbs’ hypothesis for the stated equilibrium condition and where the probability data strongly suggests that an extended equipartition principle may be formulated for some specific molecular coordinates, whose equipartition temperature does not equal the mean system temperature and a conjecture concerning which coordinates may be suitable is provided. Evidence of violations to the mesoscopic nonequilibrium thermodynamics (MNET) assumptions used without clear qualifications for a canonical distribution for internal variables are described, and possible reasons outlined, where it is found that the free dimer and atom particle kinetic energy distributions agree fully with Maxwell–Boltzmann statistics but the distribution for the relative kinetic energy of bonded atoms does not. The principle of local equilibrium (PLE) commonly used in nonequilibrium theories to model irreversible systems is investigated through NEMD simulation at extreme conditions of bond formation and breakup at the reservoir ends in the presence of a temperature gradient, where for this study a simple and novel difference equation algorithm to test the divergence theorem for mass conservation is utilized, where mass is found to be conserved from the algorithm in the presence of flux currents, in contradiction to at least one aspect of PLE in the linear domain. It is concluded therefore that this principle can be a good approximation at best, corroborating previous purely theoretical results derived from the generalized Clausius Inequality, which proved that the PLE cannot be an exact principle for nonequilibrium systems.      

On the Use of Accelerated Molecular Dynamics to Enhance Configurational Sampling in Ab Initio Simulations

We have implemented the accelerated molecular dynamics approach (Hamelberg, D.; Mongan, J.; McCammon, J. A. J. Chem. Phys. 2004, 120 (24), 11919) in the framework of ab initio MD (AIMD). Using three simple examples, we demonstrate that accelerated AIMD (A-AIMD) can be used to accelerate solvent relaxation in AIMD simulations and facilitate the detection of reaction coordinates: (i) We show, for one cyclohexane molecule in the gas phase, that the method can be used to accelerate the rate of the chair-to-chair interconversion by a factor of ～1 × 10(5), while allowing for the reconstruction of the correct canonical distribution of low-energy states; (ii) We then show, for a water box of 64 H(2)O molecules, that A-AIMD can also be used in the condensed phase to accelerate the sampling of water conformations, without affecting the structural properties of the solvent; and (iii) The method is then used to compute the potential of mean force (PMF) for the dissociation of Na-Cl in water, accelerating the convergence by a factor of ～3-4 compared to conventional AIMD simulations.(2) These results suggest that A-AIMD is a useful addition to existing methods for enhanced conformational and phase-space sampling in solution. While the method does not make the use of collective variables superfluous, it also does not require the user to define a set of collective variables that can capture all the low-energy minima on the potential energy surface. This property may prove very useful when dealing with highly complex multidimensional systems that require a quantum mechanical treatment.     

Correlated motions analysis from molecular dynamics trajectories: Statistical accuracy on the determination of canonical correlation coefficients

A numerical study of the accuracy on the determination of a unique global canonical correlation coefficient C between two groups of random variables is presented as a function of the number of variables of both groups (n and m, respectively), of the sampling size N and of the actual level of correlation C between the groups. The method used to estimate C has been already described (Briki, F.; Genest, D. Biophys Chem 1994, 52, 35–43; and J Biomol Struct Dynam 1995, 12, 1063–1082), and is implemented in the home made program TECOR. To check the accuracy on the estimation of C for given values of n, m, N, and C, samples of the random variables are synthesized (with known C), then TECOR is used to get an estimated value M of the global canonical coefficient, which is compared to the actual value C. Special attention is paid to the application of the method to the analysis from molecular dynamics simulation trajectories of concerted motions of two groups of atoms (not larger than about 20 atoms) in the course of internal deformation of biopolymers. It is found that there is a good agreement between M and C for moderate and high correlation (C≥0.35), provided that at least about 2000 configurations are stored during the molecular dynamics simulation. If C is smaller than 0.35, the method can overestimate its value if the number of configurations is not increased, especially for larger groups. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1571–1576, 1999     

Localization measure and maximum delocalization in molecular systems

A mathematically well-defined measure of localization is presented based on Mulliken's orbital populations. It is shown that this quantity equals 1 for core- and lone-pair orbitals, 2 for two-atomic bonds, 6 for benzene rings, etc., and it is applicable for delocalized canonical HF orbitals as well. The definition of this quantity is general in the sense that ab initio MOS with overlapping AO expansion, and semiempirical wave functions using the ZDO approximation as well, can be treated. The localization quantity is essentially “intrinsic,” i.e., no subdivision of the molecule is required. For N-electron wave functions, mean delocalization can be defined. This measure is not invariant to unitary transformations of the one-electron orbitals, characterizing in this way the localized or extended representation of the N-electron wave function. It can be proven, however, that for unitary transformed wave functions a maximum delocalization exists which depends only on the physical (N-electron) properties of the molecule. It is shown that inhomogeneous charge distribution can cause strong electron localization in molecular systems. The delocalization of the canonical Hartree–Fock orbitals, the Parr–Chen circulant orbitals, and the optimum delocalized orbitals is studied by numerical calculations in extended systems.     

Sampling free energy surfaces as slices by combining umbrella sampling and metadynamics

Metadynamics (MTD) is a very powerful technique to sample high‐dimensional free energy landscapes, and due to its self‐guiding property, the method has been successful in studying complex reactions and conformational changes. MTD sampling is based on filling the free energy basins by biasing potentials and thus for cases with flat, broad, and unbound free energy wells, the computational time to sample them becomes very large. To alleviate this problem, we combine the standard Umbrella Sampling (US) technique with MTD to sample orthogonal collective variables (CVs) in a simultaneous way. Within this scheme, we construct the equilibrium distribution of CVs from biased distributions obtained from independent MTD simulations with umbrella potentials. Reweighting is carried out by a procedure that combines US reweighting and Tiwary–Parrinello MTD reweighting within the Weighted Histogram Analysis Method (WHAM). The approach is ideal for a controlled sampling of a CV in a MTD simulation, making it computationally efficient in sampling flat, broad, and unbound free energy surfaces. This technique also allows for a distributed sampling of a high‐dimensional free energy surface, further increasing the computational efficiency in sampling. We demonstrate the application of this technique in sampling high‐dimensional surface for various chemical reactions using ab initio and QM/MM hybrid molecular dynamics simulations. Further, to carry out MTD bias reweighting for computing forward reaction barriers in ab initio or QM/MM simulations, we propose a computationally affordable approach that does not require recrossing trajectories. © 2016 Wiley Periodicals, Inc.     

Monte Carlo simulation of chain length distribution in radical polymerization with transfer reaction

In this paper, the basic principle and a Monte Carlo method are described for numerically simulating the chain-length distribution in radical polymerization with transfer reaction to monomer. The agreement between the simulated and analytical results shows that our algorithm is suitable for systems with transfer reaction. With the simulation algorithm, we confirm that transfer reaction has a similar effect as disproportionation on the molecular weight distribution in radical polymerization with continuous initiation. In the pulsed laser (PL) initiated radical polymerization with transfer reaction, the ‘waves’ on the chain-length distribution profile become weaker as the ratio of transfer reaction rate constant, k_tr, to the propagation rate constant, k_p, is increased in the case with either combination-type or disproportionation-type termination. Moreover, it seems that the combination termination has a broadening effect on the waves. Therefore, k_p can also be determined by precisely locating the inflection point L_o on the chain-length distribution profile for radical polymerization with transfer reaction, unless k_tr is large enough to smear out the waves on the chain-length distribution.     

Repeated‐annealing sampling combined with multicanonical algorithm for conformational sampling of bio‐molecules

A novel conformational sampling method (repeated‐annealing sampling method) is proposed to execute an efficient conformational sampling at a reasonable computational cost. In the method, a molecular dynamics simulation is done with repeating an elemental process. An elemental process consists of four subprocesses: high‐temperature run, annealing, room‐temperature run, and fast heating. The sampling is done automatically according to a temperature‐control schedule. The room‐temperature run is treated with the multicanonical algorithm, and the other subprocesses are done with the conventional molecular dynamics algorithm. The method, differing from the generalized ensemble methods recently developed, is not warrantable to give the canonical ensemble because of the nonphysical process in the annealing. However, we observed that the slower the annealing and the longer the high‐temperature run, the closer the sampled conformations to those of the canonical ensemble. A test was performed with tri‐N‐acetyl‐D ‐glucosamine in vacuo, and the results were compared with those from the conventional multicanonical simulation. Not only the reweighted canonical distribution function but also the energy landscape were in good agreement with those from the conventional multicanonical simulation. The potential of mean force also showed a fairly good agreement with that from the conventional multicanonical simulation in the room‐temperature region. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1098–1106, 2001     

Enhanced sampling method in molecular simulations using genetic algorithm for biomolecular systems

We propose a molecular simulation method using genetic algorithm (GA) for biomolecular systems to obtain ensemble averages efficiently. In this method, we incorporate the genetic crossover, which is one of the operations of GA, to any simulation method such as conventional molecular dynamics (MD), Monte Carlo, and other simulation methods. The genetic crossover proposes candidate conformations by exchanging parts of conformations of a target molecule between a pair of conformations during the simulation. If the candidate conformations are accepted, the simulation resumes from the accepted ones. While conventional simulations are based on local update of conformations, the genetic crossover introduces global update of conformations. As an example of the present approach, we incorporated genetic crossover to MD simulations. We tested the validity of the method by calculating ensemble averages and the sampling efficiency by using two kinds of peptides, ALA3 and (AAQAA)₃. The results show that for ALA3 system, the distribution probabilities of backbone dihedral angles are in good agreement with those of the conventional MD and replica-exchange MD simulations. In the case of (AAQAA)₃ system, our method showed lower structural correlation of α-helix structures than the other two methods and more flexibility in the backbone ψ angles than the conventional MD simulation. These results suggest that our method gives more efficient conformational sampling than conventional simulation methods based on local update of conformations. © 2018 Wiley Periodicals, Inc.     

Efficient calculation of two-dimensional adiabatic and free energy maps: Application to the isomerization of the C13C14 and C15N16 bonds in the retinal of bacteriorhodopsin

Accurate calculation of potential energy and free-energy profiles along reaction coordinates of biological processes such as enzymatic reactions or conformational changes is fundamental to the obtention of theoretical insight into protein function. We describe here the practical implementation of the Automatic Map Refinement Procedure (AMRP) and two-dimensional Weighted Histogram Analysis Method (WHAM) for efficient computation of adiabatic potential energy and free-energy maps, respectively. Methods for efficiently sampling configuration space with high-energy barriers and for removing hysteresis in the case of periodic reaction coordinates are presented. The application of these techniques to the isomerization of the C13C14 and C15N16 bonds in the retinal of bacteriorhodopsin is described. In dark-adapted bacteriorhodopsin (bR), the retinal moiety exists in two conformers, all-trans and (13,15)cis, with the latter making ≃67% of the population. This experimental free energy difference is reproduced here to within k_BT. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1644–1658, 1999     

Narrow molecular weight distribution of condensation polymers achieved by adding monomers in many batches

The Z transform method has been used to calculate the molecular weight distribution (MWD) of condensation polymers. The MWD obtained by using Z transform is explicitly discrete. The method is illustrated for two cases: (1) further polycondensation of AB prepolymers with certain initial MWD, and (2) polycondensation of AB and Ar (r is the number of A type functional groups) monomers where AB monomers are added in several batches. In the latter case, it is found that the resulting MWD is much narrower than that of one-batch polycondensation. The trick of producing narrow MWDs of condensation polymers is merely a consequence of keeping AB monomer concentration as low as possible during the reaction in order to suppress the condensation reaction between monomeric AB molecules. The theoretical prediction has been confirmed by Monte Carlo simulation. Therefore, it provides a new possible technique for obtaining narrow MWD polymers through polycondensation reactions.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

New out-of-plane angle and bond angle internal coordinates and related potential energy functions for molecular mechanics and dynamics simulations

With currently used definitions of out-of-plane angle and bond angle internal coordinates, Cartesian derivatives have singularities, at ±π/2 in the former case and π in the latter. If either of these occur during molecular mechanics or dynamics simulations, the forces are not well defined. To avoid such difficulties, we provide new out-of-plane and bond angle coordinates and associated potential energy functions that inherently avoid these singularities. The application of these coordinates is illustrated by ab initio calculations on ammonia, water, and carbon dioxide. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1067–1084, 1999     

A method for the calculation of normal-mode vibrational frequencies using symmetry coordinates. Application to the calculation of secondary deuterium isotope effects on carbocations

A method is introduced for the calculation of normal-mode vibrational frequencies of polyatomic molecules based on numerical differencing of analytical gradients in symmetry coordinates. This procedure requires a number of gradient evaluations equal to the largest number of symmetry coordinates belonging to any single irreducible representation of the molecular point group (plus a single gradient evaluation at the equilibrium configuration), which is fewer than the 3N-6 (N atoms) gradient evaluations needed for schemes based on Cartesian or internal coordinates. While the proposed method will not generally be as efficient as procedures which involve the direct calculation of energy second derivatives analytically (as are now available for single-determinant wavefunctions) it appears to be equally accurate, and it should be the method of choice for frequency calculations involving multideterminant wavefunctions for which analytical second-derivative algorithms have yet to be developed. The method is illustrated by the calculation of equilibrium secondary deuterium-isotope effects on a number of reactions involving simple carbocations.     

A radial probability density function for analysis of canonical molecular orbitals

A one‐dimensional probability density function, analogous to the atomic radial density for the hydrogen atom, r²R_nl(r), is defined for an arbitrary three‐dimensional density. It is obtained numerically by taking the derivative of a cumulative probability distribution with respect to the cubic root of the volume enclosed by each in a series of isosurfaces. Each point in the function is associated with a unique isosurface, and the isosurface associated with the maximum of the defined function represents the most probable isosurface with respect to the putative radius. This function therefore provides an objective selection criterion for a single isosurface to represent a three‐dimensional density. This technique is applied to set of canonical molecular orbitals. The selected threshold value varies from orbital to orbital, but the enclosed probability falls in the range of 20% to 55% for the reported orbitals. In all cases, the enclosed probability is much smaller than the common choices found in the literature. The concomitant smaller volume often makes possible a more localized interpretation and helps to clarify the conventional delocalized interpretation of molecular orbitals. For example, the isosurface plots selected by this method distinguish the formally bonding orbital in He₂ from the true bonding orbital in H₂. Examples from N₂, F₂, HF, H₂O, C₂H₆, and Ni(CO)₄ are also presented. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 310–321, 2000     

On the Quantitative Click Conjugation of Molecular Weight Distributions: What Can Theoretically Be Expected?

The assessment of the success of a click conjugation for block copolymer formation often occurs via ill‐defined criteria based on the qualitative shape of the size exclusion chromatography (SEC) distribution as well as unsuitable characteristics of these molecular weight distribution such as the peak molecular weight, M_p. The data presented herein illustrate that from the shape of the SEC distribution, w(log M) versus log M, of the click product—under the assumption of a 100% efficient conjugation reaction—a conclusion about the effectiveness of the click conjugation cannot be derived. It is demonstrated that under certain conditions multi‐modal molecular weight distributions are obtained from two well‐defined monomodal initial distributions. Similarly, the comparison of peak maxima between the conjugate and the initial SEC distributions is a poorly defined criterion to assess the success of a conjugation. The only reliable assessment for the success of the click coupling reaction is the true number average molecular weight, , of the SEC distribution of the click product, as well as a plot of the concentration or number distributions, c versus M, of the precursors and product polymer.

     

Graphical approach for defining natural internal coordinates

A simple, customizable connectivity scheme is rigorously defined in which pairs of atoms are classified into three categories. The tools of graph theory are used to analyze the molecular graph and to efficiently find rings and ring assemblies through a combination of pruning and homeomorphic reduction. The definition of natural internal coordinates is extended in a nonredundant fashion for the various cases of weakly interacting components and for fused ring systems. The ring system coordinates were tested and found to be superior to Z-matrix coordinates. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 504–511, 1998