Rapid anharmonic vibrational corrections derived from partial Hessian analysis

Vibrational analysis within a partial Hessian framework can successfully describe the vibrational properties of a variety of systems where the vibrational modes of interest are localized within a specific region of the system. We have developed a new approach to calculating anharmonic frequencies based on vibrational frequencies and normal modes obtained from a partial Hessian analysis using second-order vibrational perturbation theory and the transition optimized shifted Hermite method. This allows anharmonic frequencies for vibrational modes that are spatially localized to be determined at a significantly reduced computational cost. Several molecular systems are examined in order to demonstrate the effectiveness of this method including organic molecules adsorbed on the Si(100)-2×1 surface, model peptides in solution, and the C-H stretching region of polycyclic aromatic hydrocarbons. Overall, for a range of systems, anharmonic frequencies calculated using the partial Hessian approach are found to be in close agreement with the results obtained using full anharmonic calculations while providing a significant reduction in computational cost.     

Optimizing the performance of the multiconfiguration molecular mechanics method

Multiconfiguration molecular mechanics (MCMM) is a general algorithm for constructing potential energy surfaces for reactive systems (Kim, Y.; Corchado, J. C.; Villà, J.; Xing, J.; Truhlar, D. G. J. Chem. Phys. 2000, 112, 2718). This paper illustrates how the performance of the MCMM method can be improved by adopting improved molecular mechanics parameters. We carry out calculations of reaction rate constants using variational transition state theory with optimized multidimensional tunneling on the MCMM PESs for three hydrogen transfer reactions, and we compare the results to direct dynamics. We find that the MCMM method with as little as one electronic structure Hessian can describe the dynamically important regions of the ground-electronic state PES, including the corner-cutting-tunneling region of the reaction swath, with practical accuracy.     

Another way to implement the Powell formula for updating Hessian matrices related to transition structures

A way to update the Hessian matrix according to the Powell formula is given. With this formula one does not need to store the full Hessian matrix at any iteration. A method to find transition structures, which is a combination of the quasiNewton–Raphson augmented Hessian algorithm with the proposed Powell update scheme, is also given. The diagonalization of the augmented Hessian matrix is carried out by Lanczoslike methods. In this way, during all the optimization process, one avoids to store full matrices.     

Comparison of methods for finding saddle points without knowledge of the final states

Within the harmonic approximation to transition state theory, the biggest challenge involved in finding the mechanism or rate of transitions is the location of the relevant saddle points on the multidimensional potential energy surface. The saddle point search is particularly challenging when the final state of the transition is not specified. In this article we report on a comparison of several methods for locating saddle points under these conditions and compare, in particular, the well-established rational function optimization (RFO) methods using either exact or approximate Hessians with the more recently proposed minimum mode following methods where only the minimum eigenvalue mode is found, either by the dimer or the Lanczos method. A test problem involving transitions in a seven-atom Pt island on a Pt(111) surface using a simple Morse pairwise potential function is used and the number of degrees of freedom varied by varying the number of movable atoms. In the full system, 175 atoms can move so 525 degrees of freedom need to be optimized to find the saddle points. For testing purposes, we have also restricted the number of movable atoms to 7 and 1. Our results indicate that if attempting to make a map of all relevant saddle points for a large system (as would be necessary when simulating the long time scale evolution of a thermal system) the minimum mode following methods are preferred. The minimum mode following methods are also more efficient when searching for the lowest saddle points in a large system, and if the force can be obtained cheaply. However, if only the lowest saddle points are sought and the calculation of the force is expensive but a good approximation for the Hessian at the starting position of the search can be obtained at low cost, then the RFO approaches employing an approximate Hessian represent the preferred choice. For small and medium sized systems where the force is expensive to calculate, the RFO approaches employing an approximate Hessian is also the more efficient, but when the force and Hessian can be obtained cheaply and only the lowest saddle points are sought the RFO approach using an exact Hessian is the better choice. These conclusions have been reached based on a comparison of the total computational effort needed to find the saddle points and the number of saddle points found for each of the methods. The RFO methods do not perform very well with respect to the latter aspect, but starting the searches further away from the initial minimum or using the hybrid RFO version presented here improves this behavior considerably in most cases.     

Comparative study of various normal mode analysis techniques based on partial Hessians

Standard normal mode analysis becomes problematic for complex molecular systems, as a result of both the high computational cost and the excessive amount of information when the full Hessian matrix is used. Several partial Hessian methods have been proposed in the literature, yielding approximate normal modes. These methods aim at reducing the computational load and/or calculating only the relevant normal modes of interest in a specific application. Each method has its own (dis)advantages and application field but guidelines for the most suitable choice are lacking. We have investigated several partial Hessian methods, including the Partial Hessian Vibrational Analysis (PHVA), the Mobile Block Hessian (MBH), and the Vibrational Subsystem Analysis (VSA). In this article, we focus on the benefits and drawbacks of these methods, in terms of the reproduction of localized modes, collective modes, and the performance in partially optimized structures. We find that the PHVA is suitable for describing localized modes, that the MBH not only reproduces localized and global modes but also serves as an analysis tool of the spectrum, and that the VSA is mostly useful for the reproduction of the low frequency spectrum. These guidelines are illustrated with the reproduction of the localized amine‐stretch, the spectrum of quinine and a bis‐cinchona derivative, and the low frequency modes of the LAO binding protein. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Analysis of the updated Hessian matrices for locating transition structures

We present an analysis of the behavior of different updating Hessian formulas when they are used for the location and optimization of transition structures. The analysis is based on the number of iterations, the minimum of the weighted Euclidean matrix norm, and first-order perturbation theory applied to each type of Hessian correction. Finally, we give a derivation of a family of updated Hessians from the variational method proposed by Greenstadt. We conclude that the proposed family of updated Hessians is useful for the optimization of transition structures. © 1995 John Wiley & Sons, Inc.     

Searching for saddle points of potential energy surfaces by following a reduced gradient

The old coordinate driving procedure to find transition structures in chemical systems is revisited. The well-known gradient criterion, ∇E( x )= 0 , which defines the stationary points of the potential energy surface (PES), is reduced by one equation corresponding to one search direction. In this manner, abstract curves can be defined connecting stationary points of the PES. Starting at a given minimum, one follows a well-selected coordinate to reach the saddle of interest. Usually, but not necessarily, this coordinate will be related to the reaction progress. The method, called reduced gradient following (RGF), locally has an explicit analytical definition. We present a predictor–corrector method for tracing such curves. RGF uses the gradient and the Hessian matrix or updates of the latter at every curve point. For the purpose of testing a whole surface, the six-dimensional PES of formaldehyde, H₂CO, was explored by RGF using the restricted Hartree–Fock (RHF) method and the STO-3G basis set. Forty-nine minima and saddle points of different indices were found. At least seven stationary points representing bonded structures were detected in addition to those located using another search algorithm on the same level of theory. Further examples are the localization of the saddle for the HCN⇌CNH isomerization (used for steplength tests) and for the ring closure of azidoazomethine to 1H-tetrazole. The results show that following the reduced gradient may represent a serious alternative to other methods used to locate saddle points in quantum chemistry. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1087–1100, 1998     

Calculating rate constants with updated Hessians using variational transition state theory with multidimensional tunneling

Variational transition state theory with multidimensional tunneling (VTST/MT) has been used for calculating the rate constants of reactions. The updated Hessians have been used to reduce the computational costs for both geometry optimization and trajectory following procedures. In this paper, updated Hessians are used to reduce the computational costs while calculating the rate constants applying VTST/MT. Although we found that directly applying the updated Hessians will not generate good vibrational frequencies along the minimum energy path (MEP), however, we can either re-compute the full Hessian matrices at fixed intervals or calculate the Block Hessians, which is constructed by numerical one-side difference for the Hessian elements in the "critical" region and Bofill updating scheme for the rest of the Hessian elements. Due to the numerical instability of the Bofill update method near the saddle point region, we have suggested a simple strategy in which we follow the MEP until certain percentage of the classical barrier height from the barrier top with full Hessians computed and then performing rate constant calculation with the extended MEP using Block Hessians. This strategy results a mean unsigned percentage deviation (MUPD) around 10% with full Hessians computed till the point with 80% classical barrier height for four studied reactions. This proposed strategy is attractive not only it can be implemented as an automatic procedure but also speeds up the VTST/MT calculation via embarrassingly parallelization to a personal computer cluster.     

Multipath variational transition state theory: rate constant of the 1,4-hydrogen shift isomerization of the 2-cyclohexylethyl radical

We propose a new formulation of variational transition state theory called multipath variational transition state theory (MP-VTST). We employ this new formulation to calculate the forward and reverse thermal rate constant of the 1,4-hydrogen shift isomerization of the 2-cyclohexylethyl radical in the gas phase. First, we find and optimize all the local-minimum-energy structures of the reaction, product, and transition state. Then, for the lowest-energy transition state structures, we calculate the reaction path by using multiconfiguration Shepard interpolation (MSCI) method to represent the potential energy surface, and, from this representation, we also calculate the ground-state vibrationally adiabatic potential energy curve, the reaction-path curvature vector, and the generalized free energy of activation profile. With this information, the path-averaged generalized transmission coefficients <γ> are evaluated. Then, thermal rate constant containing the multiple-structure anharmonicity and torsional anharmonicity effects is calculated using multistructural transition state theory (MS-TST). The final MP-VTST thermal rate constant is obtained by multiplying k(MS-T)(MS-TST) by <γ>. In these calculations, the M06 density functional is utilized to compute the energy, gradient, and Hessian at the Shepard points, and the M06-2X density functional is used to obtain the structures (conformers) of the reactant, product, and the saddle point for computing the multistructural anharmonicity factors.     

Using bonding to guide transition state optimization

Optimization of a transition state typically requires both a good initial guess of the molecular structure and one or more computationally demanding Hessian calculations to converge reliably. Often, the transition state being optimized corresponds to the barrier in a chemical reaction where bonds are being broken and formed. Utilizing the geometries and bonding information for reactants and products, an algorithm is outlined to reliably interpolate an initial guess for the transition state geometry. Additionally, the change in bonding is also used to increase the reliability of transition state optimizations that utilize approximate and updated Hessian information. These methods are described and compared against standard transition state optimization methods. © 2015 Wiley Periodicals, Inc.     

Improved configuration space sampling: Langevin dynamics with alternative mobility

We present a new and efficient method for determining optimal configurations of a large number (N) of interacting particles. We use a coarse-grained stochastic Langevin equation in the overdamped limit to describe the dynamics of this system and replace the standard mobility by an effective space dependent inverse Hessian correlation matrix. Due to the analogy of the drift term in the Langevin equation and the update scheme in Newton's method, we expect accelerated dynamics or improved convergence in the convex part of the potential energy surface Phi. The stochastic noise term, however, is not only essential for proper thermodynamic sampling but also allows the system to access transition states in the concave parts of Phi. We employ a Broyden-Fletcher-Goldfarb-Shannon method for updating the local mobility matrix. Quantitative analysis for one and two dimensional systems shows that the new method is indeed more efficient than standard methods with constant effective friction. Due to the construction, our effective mobility adapts high values/low friction in configurations which are less optimal and low values/high friction in configurations that are more optimal.     

C library for topological study of the electronic charge density

The topological study of the electronic charge density is useful to obtain information about the kinds of bonds (ionic or covalent) and the atom charges on a molecule or crystal. For this study, it is necessary to calculate, at every space point, the electronic density and its electronic density derivatives values up to second order. In this work, a grid‐based method for these calculations is described. The library, implemented for three dimensions, is based on a multidimensional Lagrange interpolation in a regular grid; by differentiating the resulting polynomial, the gradient vector, the Hessian matrix and the Laplacian formulas were obtained for every space point. More complex functions such as the Newton–Raphson method (to find the critical points, where the gradient is null) and the Cash–Karp Runge–Kutta method (used to make the gradient paths) were programmed. As in some crystals, the unit cell has angles different from 90°, the described library includes linear transformations to correct the gradient and Hessian when the grid is distorted (inclined). Functions were also developed to handle grid containing files (grd from DMol® program, CUBE from Gaussian® program and CHGCAR from VASP® program). Each one of these files contains the data for a molecular or crystal electronic property (such as charge density, spin density, electrostatic potential, and others) in a three‐dimensional (3D) grid. The library can be adapted to make the topological study in any regular 3D grid by modifying the code of these functions. © 2012 Wiley Periodicals, Inc.     

Comparison of double-ended transition state search methods

While a variety of double-ended transition state search methods have been developed, their relative performance in characterizing complex multistep pathways between structurally disparate molecular conformations remains unclear. Three such methods (doubly-nudged elastic band, a string method, and a growing string method) are compared for a series of benchmarks ranging from permutational isomerizations of the seven-atom Lennard-Jones cluster (LJ(7)) to highly cooperative LJ(38) and LJ(75) rearrangements, and the folding pathways of two peptides. A database of short paths between LJ(13) local minima is used to explore the effects of parameters and suggest reasonable default values. Each double-ended method was employed within the framework of a missing connection network flow algorithm to construct more complicated multistep pathways. We find that in our implementation none of the three methods definitively outperforms the others, and that their relative effectiveness is strongly system and parameter dependent.     

A detailed theoretical study on the excited‐state hydrogen‐bonding dynamics and the proton transfer mechanism for a novel white‐light fluorophore

In this work, density functional theory (DFT) and time‐dependent DFT (TDDFT) methods were used to investigate the excited‐state dynamics of the excited‐state hydrogen‐bonding variations and proton transfer mechanism for a novel white‐light fluorophore 2‐(4‐[dimethylamino]phenyl)‐7‐hyroxy‐6‐(3‐phenylpropanoyl)‐4H‐chromen‐4‐one ( 1 ). The methods we adopted could successfully reproduce the experimental electronic spectra, which shows the appropriateness of the theoretical level in this work. Using molecular electrostatic potential (MEP) as well as the reduced density gradient (RDG) versus the product of the sign of the second largest eigenvalue of the electron density Hessian matrix and electron density (sign[λ₂]ρ), we demonstrate that an intramolecular hydrogen bond O1–H2···O3 should be formed spontaneously in the S₀ state. By analyzing the chemical structures, infrared vibrational spectra, and hydrogen‐bonding energies, we confirm that O1–H2·O3 should be strengthened in the S₁ state, which reveals the possibility of an excited‐state intramolecular proton transfer (ESIPT) process. On investigating the excitation process, we find the S₀ → S₁ transition corresponding to the charge transfer, which provides the driving force for ESIPT. By constructing the potential energy curves, we show that the ESIPT reaction results in a dynamic equilibrium in the S₁ state between the forward and backward processes, which facilitates the emission of white light.     

How good is a Broyden–Fletcher–Goldfarb–Shanno-like update Hessian formula to locate transition structures? Specific reformulation of Broyden–Fletcher–Goldfarb–Shanno for optimizing saddle points

Based on a study of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update Hessian formula to locate minima on a hypersurface potential energy, we present an updated Hessian formula to locate and optimize saddle points of any order that in some sense preserves the initial structure of the BFGS update formula. The performance and efficiency of this new formula is compared with the previous updated Hessian formulae such as the Powell and MSP ones. We conclude that the proposed update is quite competitive but no more efficient than the normal updates normally used in any optimization of saddle points. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 349–362, 1998     

Parallel algorithm for efficient calculation of second derivatives of conformational energy function in internal coordinates

A parallel algorithm for efficient calculation of the second derivatives (Hessian) of the conformational energy in internal coordinates is proposed. This parallel algorithm is based on the master/slave model. A master processor distributes the calculations of components of the Hessian to one or more slave processors that, after finishing their calculations, send the results to the master processor that assembles all the components of the Hessian. Our previously developed molecular analysis system for conformational energy optimization, normal mode analysis, and Monte Carlo simulation for internal coordinates is extended to use this parallel algorithm for Hessian calculation on a massively parallel computer. The implementation of our algorithm uses the message passing interface and works effectively on both distributed-memory parallel computers and shared-memory parallel computers. We applied this system to the Newton–Raphson energy optimization of the structures of glutaminyl transfer RNA (Gln-tRNA) with 74 nucleotides and glutaminyl-tRNA synthetase (GlnRS) with 540 residues to analyze the performance of our system. The parallel speedups for the Hessian calculation were 6.8 for Gln-tRNA with 24 processors and 11.2 for GlnRS with 54 processors. The parallel speedups for the Newton–Raphson optimization were 6.3 for Gln-tRNA with 30 processors and 12.0 for GlnRS with 62 processors. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1716–1723, 1998     

Identification and thermodynamic treatment of several types of large-amplitude motions

We present a partially automated method for the thermodynamic treatment of large-amplitude motions. Starting from the molecular geometry and the Hessian matrix, we evaluate anharmonic partition functions for selected vibrational degrees of freedom. Supported anharmonic vibration types are internal rotation and inversion (oscillation in a double-well potential). By heuristic algorithms, we identify internal rotations in most cases automatically from the Hessian eigenvectors, and we also estimate the parameters of anharmonic partition functions (e.g., potential barrier, periodicity, and symmetry number) with thermodynamically sufficient precision. We demonstrate the validity of our schemes by comparison to pointwise calculated ab initio potential curves.     

On a partial differential equation method for determining the free energies and coexisting phase compositions of ternary mixtures from light scattering data

In this paper we present a method for determining the free energies of ternary mixtures from light scattering data. We use an approximation that is appropriate for liquid mixtures, which we formulate as a second-order nonlinear partial differential equation. This partial differential equation (PDE) relates the Hessian of the intensive free energy to the efficiency of light scattering in the forward direction. This basic equation applies in regions of the phase diagram in which the mixtures are thermodynamically stable. In regions in which the mixtures are unstable or metastable, the appropriate PDE is the nonlinear equation for the convex hull. We formulate this equation along with continuity conditions for the transition between the two equations at cloud point loci. We show how to discretize this problem to obtain a finite-difference approximation to it, and we present an iterative method for solving the discretized problem. We present the results of calculations that were done with a computer program that implements our method. These calculations show that our method is capable of reconstructing test free energy functions from simulated light scattering data. If the cloud point loci are known, the method also finds the tie lines and tie triangles that describe thermodynamic equilibrium between two or among three liquid phases. A robust method for solving this PDE problem, such as the one presented here, can be a basis for optical, noninvasive means of characterizing the thermodynamics of multicomponent mixtures.     

Convergence optimization of restricted open-shell self-consistent field calculations

Different self-consistent field (SCF) iteration schemes for open-shell systems are discussed. After a brief summary of the well-known level shifting and damping procedure, we describe the quadratically convergent SCF (QCSCF) approach based on the gradient and the Hessian matrix in a space of orbital rotation parameters. An analytical expression for the latter is derived for the general many-shell case. Starting from the expression for the energy change obtained by the QCSCF method, we then present a simplified direct procedure avoiding matrix diagonalization but also the difficulties of the QCSCF method in handling the Hessian matrix. Numerical calculations on some open-shell systems involving transition-metal complexes show that this method leads to rapid and reliable convergence of the iteration process in cases where the usual SCF procedure of iterative diagonalization tends to diverge. © 1997 John Wiley & Sons, Inc. Int J Quant Chem 62: 617–637, 1997     

The cholesteric-smectic A transition

Instead of studying the cholesteric-smectic A transition from the analogy with superconductors, we present a direct investigation of the properties of this transition. The transition is of first order and a very simple relationship between the four basic lengths (correlation length, penetration length, layer spacing of the smectic structure, and pitch of the cholesteric structure) exists at the transition temperature. We calculate the wall energy separating a cholesteric region from a smectic region in two different configurations. For one of them, we recover the definition of the type I (positive wall energy) and type II (negative wall energy) smectics A. We do not find a specific value of the Ginzburg-Landau parameter kc which gives the distinction between the two types of smectics; kc is a function of the product of the perpendicular correlation length and the wave vector of the smectic density modulation. We recover the well known value from the superconductor analogy, kc 1/tick 2, only if this product is large and the rotation of the director inside the wall is small. For the second configuration, the wall energy is always positive and there is no analogy with superconductors. Finally, we discuss briefly two different possibilities of mixed phases.