key: cord-0788290-4j5pzjjz authors: Tao, Yunwen; Zou, Wenli; Nanayakkara, Sadisha; Freindorf, Marek; Kraka, Elfi title: A revised formulation of the generalized subsystem vibrational analysis (GSVA) date: 2021-03-09 journal: Theor Chem Acc DOI: 10.1007/s00214-021-02727-y sha: ca6720f5965ca801fa74a0959c58ad1538afe693 doc_id: 788290 cord_uid: 4j5pzjjz In this work, a simplified formulation of our recently developed generalized subsystem vibrational analysis (GSVA) for obtaining intrinsic fragmental vibrations (J Chem Theory Comput 14:2558, 2018) is presented. In contrast to the earlier implementation, which requires the explicit definition of a non-redundant set of internal coordinate parameters to be constructed for the subsystem, the new implementation circumvents this process by employing massless Eckart conditions to the subsystem fragment paired with a Gram–Schmidt orthogonalization to span the same internal vibration space indirectly. This revised version of GSVA (rev-GSVA) can be applied to equilibrium structure as well as transition state structure, and it has been incorporated into the open-source package UniMoVib (https://github.com/zorkzou/UniMoVib). SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s00214-021-02727-y. Molecular vibrations described by normal modes are in general delocalized over the molecular structure of interest [1] . This delocalization hampers the direct comparison between the vibrations of a molecule in different environments (e.g., gas phase and solution) or when this molecule binds to a host compound. Over the past two decades, enormous effort has been put into finding "localized" normal vibrational modes (which can also be called "localized modals" in analogue to "localized orbitals" [2] ), belonging to a molecule within a noncovalently bonded complex or a fragment of one molecule, including partial Hessian diagonalization [3] , partial Hessian vibrational analysis (PHVA) [4, 5] , mobile block Hessian (MBH) [6] [7] [8] [9] [10] , vibrational subsystem analysis (VSA) [11] [12] [13] , local Hessian transformation [14] , just to name a few [15] . However, all these approaches share a common deficiency-the unphysical partitioning of the full Hessian matrix, which causes information loss about the interaction between the subsystem and its environment. Recently, we proposed the generalized subsystem vibrational analysis (GSVA) as a new solution to obtain intrinsic fragmental vibrations [16, 17] . The key feature of GSVA compared to its predecessors lies in avoiding the partitioning of the full Hessian matrix. Instead, GSVA extracts for the subsystem a unique effective Hessian matrix x sub where x is the full Hessian matrix expressed in Cartesian coordinates of dimension ( 3N × 3N ) for the whole molecular system being composed of N atoms including n atoms of the target subsystem and N − n atoms of the environment. The Wilson -matrices [1] ′ and ′ sub define a non-redundant set of ( 3n − k sub ) internal coordinates for the subsystem fragment in rows with full 3N columns and truncated 3n columns (excluding the environment atoms), respectively. k sub is the total number of rotations and translations for the subsystem being 5 or 6 depending on whether the subsystem geometry is linear or nonlinear. ( x ) + is the Moore-Penrose inverse [18] of x , which is singular. x sub on the left-hand side is a symmetric matrix of dimension (3n × 3n), and it has exactly k sub zero eigenvalues. The † superscript denotes matrix transpose. With the effective Hessian matrix x sub expressed in Cartesian coordinates, the conventional normal mode analysis (NMA) machinery [19, 20] which is widely implemented in quantum chemical packages can be employed to calculate for the subsystem a new type of localized normal modes, which we coined intrinsic fragmental vibrations. The reason why these normal vibrations are called "intrinsic" is due to the fact that the effective Hessian matrix x sub retains the curvature of the potential energy surface (PES) in the direction defined by any internal coordinate within the subsystem [16, 17] . In other words, the subsystem fragment "feels" exactly the same curvature of the PES as the whole system being described with the full Hessian matrix x . This property of x sub endows our GSVA method with a solid physical basis [16] . As discussed in our earlier work on the original GSVA implementation [16] , GSVA requires a complete and nonredundant set of internal coordinate parameters, and its Wilson-matrices (see Eq. 1) to span the internal vibration space of the subsystem. However, the construction of the non-redundant parameter set is nontrivial and it needs either judicious selection of parameters manually with expert knowledge, or a dedicated algorithm which automatically selects the non-redundant parameter set from a series of redundant set of parameters in a trial-and-error manner. In this work, we propose an alternative formulation of GSVA, which can save the effort of constructing the nonredundant parameter set for the subsystem. The new implementation replaces the Wilson-matrices in Eq. 1 with a different matrix, which also spans the internal vibrational space of the subsystem via the following procedure. First, we apply to the subsystem fragment with its Cartesian coordinates collected in a 3n × 1 column vector cart the massless (assuming all atomic masses are identical) Eckart conditions [21, 22] to generate a set of five or six translational and rotational vectors which are orthonormal to each other; see Eq. 2: where k equals 5 or 6 depending on whether the subsystem is linear or not and i is a column vector of length 3n. Next, a Gram-Schmidt orthonormalization is conducted on tr.+ro. to generate n vib = 3n − k remaining vectors collected in where j is a column vector of length 3n. It has to be noted that matrix is equivalent to matrix � † sub in spanning the internal coordinate/vibration space of the subsystem. In order to obtain the equivalent matrix of � † , we pad each column vector j with 3(N − n) zeros associated with environmental atoms, resulting in the matrix full with the dimension (3N × n vib ) . In this way, the effective Hessian matrix x sub for the subsystem can be written as Then, the conventional NMA machinery is applied to obtain the intrinsic fragmental vibrations for the subsystem as in our earlier formulation [16] . The new formulation has two major advantages for practical implementation. (1) It avoids the complicated process of finding the complete and nonredundant internal coordinate parameter set for the subsystem; (2) the code for finding translation/rotation vectors from the Eckart conditions and conducting Gram-Schmidt orthonormalization in a modern quantum chemical package can be reused, which facilitates implementing GSVA. In order to distinguish from the original implementation of GSVA, this revised implementation was named rev-GSVA. As a showcase example, we have employed rev-GSVA to calculate the intrinsic fragmental vibrations of the methane (CH 4 ) molecule in (1) methane-intercalated B 36 N 36 complex (Fig. 1a) , (2) methane-intercalated C 60 structure (Fig. 1b ) [23] and (3) gas phase as reference. Unlike the methaneintercalated C 60 complex, the methane-intercalated B 36 N 36 system [24] has not been synthesized experimentally so far, and it is interesting to compare the intrinsic fragmental vibrations of the methane molecule in B 36 N 36 and C 60 in order to explore the different encapsulation effect. These three molecular systems were modeled at the M06-2X/6-31G(d,p) level [25] [26] [27] with Grimme's D3(0) dispersion correction [28] using the Gaussian 16 package [29] . The results in Table 1 show that the methane molecules encapsulated inside the two cages retain T d symmetry, as the reference methane molecule in gas phase. The non-degenerate A 1 vibration describes the symmetric stretching of four C-H bonds. The doubly degenerate E modes specify the relative turnstile twisting motions of two H-C-H fragments. The triply degenerate T 2 modes (1-3) with lower frequencies denote the bending of methane, while the other triply degenerate T 2 modes (1'-3') with higher frequencies are antisymmetric stretching motions of four C-H bonds. We found the (4) x sub = ( † full ( x ) + full ) −1 † . respectively. This means the fullerene cage could strengthen the C-H bonds of methane, while the B 36 N 36 weakens the C-H bonds of the contained methane molecule. One might argue that comparing the normal mode frequencies of the whole system (the way spectroscopists usually adopt) could lead to similar conclusion because the methane molecule is well separated from the cage structure; however, one needs to note that only the (rev-)GSVA method could provide the localized normal modes and frequencies which can be legitimately comparable across different systems containing the same target subsystem. In addition to local minima on the PES, we have also tested a first-order saddle point (i.e., transition state, TS) structure. As it has been proven in our earlier work [16] that GSVA retains the curvature of the PES, it is of interest to explore whether (rev-)GSVA can retain the imaginary vibrational mode specifying the bond breaking/forming in the subsystem. In this pilot study, we investigated the TS of the chemical reaction involving a potential -ketoamide inhibitor of the SARS-CoV-2 main protease (Mpro), which is assumed to inhibit the activity of SARS-CoV-2 virus by blocking viral replication [30] . According to a recent X-ray structure of the -ketoamide SARS-CoV-2 main protease (Mpro) complex [30], ketoamide and enzyme are linked via a cysteine side chain of the enzyme. According to the suggested catalytic mechanism, the chemical reaction starts with a nucleophilic attack of the sulfur atom onto a C=O carbon atom of -ketoamide moiety, which is followed by proton transfer from the -SH group to a nearby oxygen atom of the inhibitor, as shown in Fig. 2 . Work is in progress to model the reaction in the enzyme. We started with a minimal subsystem of 3 atoms containing the proton and its donor/acceptor atoms, and rev-GSVA was applied to calculate the corresponding intrinsic fragmental vibrations. Surprisingly, no imaginary frequency exists for this small 3-atom subsystem. This is probably due to the fact that the reaction center should also include the carbon atom to draw a more complete picture about bond forming and breaking. However, one imaginary frequency starts to emerge when more surrounding atoms are included into the subsystem (see Table 2 ) and the imaginary frequency value quickly converges to that of the full system when the subsystem contains 20 atoms. This result indicates that if the normal vibration is localized in a particular part of the molecule (e.g., bond breaking/forming or C=O bond stretching), (rev-)GSVA is expected to reproduce this vibration using a subsystem containing this segment and a few surrounding atoms. This valuable feature of (rev-)GSVA can lead to important first insights into the role of surrounding atoms for the reaction mechanism by analyzing the TS, before starting a more complex reaction path following procedure. We have implemented the new formulation of GSVA (rev-GSVA) introduced in this work into the open-source package UniMoVib [20] for interested readers to use in their own research. The online version contains supplementary material available at https ://doi.org/10.1007/s0021 4-021-02727 -y. Molecular vibrations: the theory of infrared and raman vibrational spectra SCF CI calculations for vibrational eigenvalues and wavefunctions of systems exhibiting fermi resonance Computation of vibrational frequencies for adsorbates on surfaces Partial Hessian vibrational analysis: the localization of the molecular vibrational energy and entropy Partial Hessian vibrational analysis of organic molecules adsorbed on Si(100) Vibrational modes in partially optimized molecular systems Cartesian formulation of the mobile block Hessian approach to vibrational analysis in partially optimized systems Normal modes for large molecules with arbitrary link constraints in the mobile block Hessian approach Mobile block Hessian approach with adjoined blocks: an efficient approach for the calculation of frequencies in macromolecules Calculating reaction rates with partial Hessians: validation of the mobile block Hessian approach Vibrational subsystem analysis: a method for probing free energies and correlations in the harmonic limit Comparative study of various normal mode analysis techniques based on partial Hessians Comparative study of various normal mode analysis techniques based on partial Hessians An effective procedure for analyzing molecular vibrations in terms of local fragment modes Localizing normal modes in large molecules Recovering intrinsic fragmental vibrations using the generalized subsystem vibrational analysis Decoding chemical information from vibrational spectroscopy data: local vibrational mode theory A generalized inverse for matrices White paper: vibrational analysis in Gaussian UniMoVib: a unified interface for molecular harmonic vibrational frequency calculations Some studies concerning rotating axes and polyatomic molecules The kinetic energy of polyatomic molecules Methane-intercalated C 60 : preparation, orientational ordering, and structure Theoretical infrared, Raman, and optical spectra of the B 36 N 36 cage The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four m06-class functionals and 12 other functionals Self-consistent molecular orbital methods. XV. Extended Gaussian-type basis sets for lithium, beryllium, and boron Self-consistent molecular orbital methods. XII. Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved -ketoamide inhibitors Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations Acknowledgements This work was financially supported by the National Science Foundation Grants CHE 1464906. W. Z. was supported by the National Natural Science Foundation of China (Grant No. 21673175) and the Double First-Class University Construction Project of Northwest University. We thank SMU for providing supercomputing resources.