{'fraction_non_alphanumeric': 0.02596565066865793, 'fraction_numerical': 0.01035231824044532, 'mean_word_length': 4.736565420560748, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient analytic second derivative of electrostatic embedding QM/MM energy: normal mode analysis of plant cryptochrome

Analytic second derivatives of electrostatic embedding (EE) quantum mechanics/molecular mechanics (QM/MM) energy are important for performing vibrational analysis and simulating vibrational spectra of quantum systems interacting with an environment represented as a classical electrostatic potential. The main bottleneck of EE-QM/MM second derivatives is the solution of coupled perturbed equations for each MM atom perturbation. Here, we exploit the Q-vector method {[J. Chem. Phys., 151, 041102 (2019)]} to workaround this bottleneck. We derive the full analytic second derivative of the EE-QM/MM energy, which allows to compute QM, MM and QM-MM Hessian blocks in an efficient and easy to implement manner. To show the capabilities of our method, we compute the normal modes for the full \textit{arabidopsis thaliana} plant cryptochrome. We show that the flavin adenine dinucleotide vibrations (QM subsystem) strongly mix with protein modes. We compute approximate vibronic couplings for the lowest bright transition, from which we extract spectral densities and the homogeneous broadening of FAD absorption spectrum in protein using vibrationally resolved electronic spectrum simulations.


Introduction
Analytic derivatives of the energy revolutionized quantum chemistry, 1 giving access to gradients and Hessians (respectively first and second derivative with respect to nuclear positions) to explore potential energy surfaces, 2 in addition to derivatives with respect to electromagnetic fields, from which properties are nowadays routinely extracted.
The development of mixed quantum mechanics (QM) and molecular mechanics (MM) models by Warshel, Levitt and Karplus, allowed a more realistic simulation of quantum systems embedded in complex environments. [3][4][5][6] Analytic derivatives of QM/MM energy have essentially focused on analytic gradient, which have been implemented for the main self-consistent field (SCF) and post-SCF QM methods. [7][8][9][10][11][12][13][14] Analytic second derivatives of QM/MM energy are less commonly available. As a matter of fact, most of the QM/MM Hessian constructions are based on the so-called Partial Hessian approach (see 15 and references therein) in which only the QM degrees of freedom are considered. A first derivation of an analytic Hessian was introduced by Morokuma and coworkers in the framework of the substractive ONIOM QM/MM scheme. 16 One year later, Cui and Karplus proposed an analytic Hessian evaluation in their additive QM/MM approach. 17,18 The latter method is expensive, mainly due to the fact that the Hessian calculation involves solving a coupled perturbed (CP) SCF equation not only for the QM atoms but also for the MM atoms. Thus, the method of Cui and Karplus resorted to approximations of the CPSCF equations. More recently, Ghysels, Woodcock III and coworkers developed a coarse-graining method for the MM environment to reduce the number of CPSCF equations to solve. 19 A Hessian with fragment molecular orbital method has been derived by Nakata and coworkers. 20 Other approximate second derivative methods have been proposed by Sakai and Morita, 21 and Gadre and coworkers. 22 The most up-to-date formulation of the Hessian is certainly the one by Capelli and coworkers, especially derived for their very promising polarizable QM/MM method which can be used to calculate infra-red spectra of very large systems in the so-called Partial Hessian Vibrational Approach, hence ignoring the major part of the Hessian that only depends on the MM atoms. 23 Herein, we develop a new method for calculating the complete second derivatives of the QM/MM energy with additive QM-MM interaction terms. We base our derivation on singledeterminant SCF wavefunctions characterizing a QM charge density interacting with a classical point-like atomic environment via the Electrostatic Potential Fitting (ESPF) interaction Hamiltonian developed by one of us. 9,10 For that purpose, we apply the CPSCF Q-vector method for MM atoms, 24 avoiding one of the main bottlenecks of second derivatives of QM-MM energies. The method works even when the QM and MM regions are covalently linked via the use of the link atom method of Morokuma and coworkers. 16 For testing our implementation, we apply this method to perform a normal mode analysis of arabidopsis thaliana cryptochrome, 25 and we compute the vibronic couplings, which allows us to efficiently construct spectral densities and account of homogeneous broadening of the spectrum via vibrationally resolved spectrum simulations.
2 Analytic Derivatives of the ESPF QM/MM energy We use the notation R = (x, y, z) to symbolize the coordinates of QM atom, andR = (x,ỹ,z) to symbolize the coordinates of a MM atom. The partial derivatives with respect to a coordinate will be written in the short-hand notation The indexes i, j, · · · will refer to occupied molecular orbitals (MO), a, b, · · · will refer to virtual MOs, p, q, · · · to generic MOs (both occupied and virtual), and µ, ν, · · · to atomic orbital (AO) index. Capital roman letters will be used for QM atom indexes.

ESPF QM/MM energy
In the ESPF method, the QM/MM energy can be written as in which E QM is the internal energy of the QM subsystem whose electron density is polarized by the MM subsystem (thus containing the energy contribution due to the electrostatic interaction with the MM subsystem), while E M M is the internal energy of the MM part complemented with classical interaction terms between QM and MM atoms. The electrostatic interaction energy in the ESPF method is evaluated as ∆E = Tr [Ph], in which P is the atomic density and h is the QM/MM interaction hamiltonian. 9,10 This interaction energy is also expressed as the classical electrostatic interaction between quantum-derived atomic multipoles q A , µ A , . . . and the MM electrostatic potential φ A , the electric field vector ∇φ A , . . .

∆E =
in which the index A runs over the number of QM atoms N QM . For the sake of simplicity, we will restrict the multipolar expansion (Eq. 3) to the order 0 in the following, even though our method can be easily generalized to higher orders. 24 The expressions for the derivatives of the QM-only and MM-only energy terms will be omitted here, since they have been lengthy discussed in the literature. 9,26,27 Here, we will focus on the electrostatic interaction term and its derivatives. At the order 0 of the QM density of charges, each atomic charge q A is expressed as the sum of the nuclear charge Z A and the mean value of an atomic one-electron operator, whose matrix elements are denoted Q A,µν . The interaction energy in atomic orbital basis reads, in which we have introduced a new one-electron operator whose elements are:, In a nutshell, the ESPF method consists in fitting the Q A,µν matrix elements to reproduce the purely electronic QM electrostatic integrals calculated on a grid surrounding the QM subsystem. 9

First derivative
The analytic first derivatives of the ESPF interaction term were first derived in Ref. 9. For completeness, we reexpress them with the current notation. The first derivative of Eq. 4 can be expressed as a vector gradient, in which the first block corresponds to the derivative with respect to QM atom positions and the second block to MM atoms. The expression of the energy gradient with respect to QM coordinates is given by As usual, the first term is not calculated thanks to the introduction of the so-called energyweighted density matrix −PhPS x . 9,10,28 The derivative of the ESPF operator element h µν is given in Appendix 5.1. Regarding the derivative with respect to the MM coordinates, the energy gradient can be expressed as noting that Sx = 0 for MM perturbations.

Second derivative
The analytic second derivative of Eq. 4 can be expressed as a Hessian matrix in which the ∆E xy is the block of second derivatives with respect to QM atom coordinates, ∆Exỹ is the second derivative block with respect to MM atom coordinates, and ∆E xỹ = ∆Ex y is the second derivative block with respect to one MM and one QM atom coordinates. The analytic formula for the ∆E xy second derivative is given by The second derivative of the Fock operator requires the second derivative of the ESPF operator, which is given in 5.1. This second derivative block also requires the density matrix derivative P y with respect to QM atom coordinates, which is obtained by solving a CPSCF equation for each QM coordinate. 29 The ∆E xỹ block over a QM and an MM atom coordinates is given by (11) The explicit form of the second derivative of the ESPF operator is given in Appendix 5.1. This energy derivative requires the explicit construction of Px. These would require the solution of a CPSCF equation for all MM atom coordinates, which is prohibitive. A more convenient form to implement, is to compute the equivalent ∆Ex y derivative block, which is given by (12) in which the charge derivative is given by In this case, only derivatives of the density matrix with respect to QM coordinates are required, which is already constructed for the QM energy second derivative. 26 Finally, the ∆Exỹ derivatives over MM atom coordinates is given by requires the derivatives of the ESPF atomic charges, which would in principle require the density derivatives The construction of Pỹ would imply solving a coupled-perturbed type equation for each MM perturbation. Usually, the number of MM atoms is hundreds to thousands, becoming rapidly intractable. In Ref. 24, we showed how to avoid this by defining an auxiliary set of coupled-perturbed equations (Q-vector method) that only scales with the number of QM atoms, and from which we can efficiently reconstruct the charge derivatives, whereQ B,ia is the solution of the Q-vector coupled-perturbed equation, 24 which allows to rewrite the MM contribution to second derivative as This equation can be interpreted in terms of classical atomistic mechanics, Then, Q AB can be interpreted as a tensor that, given an electric charge Q A on QM atom A interacting with a gradient of the external potential, will give the "effective" charge on QM center B that interacts with another gradient of the external field.

Link atoms
For molecular systems in which there exist covalent bonds at the interface between the MM and QM subsystems, unphysical free valences are created in the QM side. Morokuma and coworkers proposed to use a fictitious link atom in order to saturate these dangling bonds, most often using a hydrogen atom. The link atom is constrained to remain positioned between the frontier MM and QM atoms, 16 in which g is a parameter that fixes the coordinates of the link atom between the covalently bonded QM and MM atoms. From the QM perspective, the "auxiliary" link atom is considered as a QM atom, and therefore, N QM includes all real QM atoms plus the link atoms. Consequently, there will be an extra contribution for the QM and MM frontier atoms that comes from the link atom. In the case of first energy derivatives, this is given by In other words, the link atom gradient is projected onto the frontier QM and MM atoms. Similarly, the contribution of the link atom to the MM and QM second derivatives can be expressed as

Computational Details
All electronic structure computations have been performed with a modified version of Gaus-sian16 including the Q-vector method and the ESPF energy, gradient and Hessian. The MM information is obtained via an interface with a modified version of Tinker 6.3.3. 9,30,31 The vibrationally resolved spectrum for FAD has been obtained by Fourier transforming the time-dependent auto-correlation function for the multi-dimensional harmonic oscillator, as implemented in Gaussian16. 30,32 The Franck-Condon approximation of the auto-correlation function has been integrated on 2 20 time steps, over an interval of about 10.5 ps. A Gaussian broadening of 235 cm −1 has been applied to the vibrational peaks.
All QM calculations have been performed at the DFT and full TDDFT levels, 33

Results and Discussion
As a first example of our method, we have computed the frequencies of flavin adenine dinucleotide at the QM/MM level. In the QM region, we have placed the isoalloxazine ring while the rest of the molecule is considered MM. This selection requires the introduction of a link atom to separate the QM and MM regions, which has been placed 3 bonds far away from the isoalloxazine ring to avoid overpolarization, as shown in Fig. 1. In this figure, we show the frequencies computed at the QM/MM level compared to the QM frequencies, which we take as reference. For the sake of comparison, we also show the pure MM frequencies. Overall, the quality of the QM/MM and MM frequencies is in very good accordance with the QM frequencies. This is expected since the MM parameters have been obtained from purely ab initio data. Some appreciable errors are observed for the highest frequency modes, as well as the mid-range. The QM/MM frequencies follow a similar trend than the MM frequencies, indicating that these modes have simply a small mixing and remain rather purely MM or purely QM. Indeed, the most delocalized frequencies are found only in the low frequency range (see further below). This results indicate that the choice of force field is indeed critical for obtaining a good spectrum of frequencies that is comparable to the ab initio spectrum.
To further illustrate the potential of our method, we have performed a normal mode analysis of plant cryptochrome, including vi- bronic coupling and homogeneous broadening of absorption spectrum using vibrationally resolved spectra simulations. For this, we have selected chain A of arabidopsis thaliana cryptochrome (PDB code 2J4D, see Fig. 2). Crystallographic water molecules have been kept in the model. Hydrogens have been added for standard aminoacid protonation states at physiological pH. The entire FAD coenzyme was defined as the QM subsystem. FAD is not covalently linked to the rest of the protein, thus not requiring any special care of the QM/MM boundary. The initial PDB structure was pre-optimized at the Amber99 level using microiterations. Subsequent geometry optimization, hessian calculation and normal mode analysis were carried out at B3LYP/TZVP/TZVPFit/Amber99 level for the entire protein model. One of the main advantages of the method presented here is that the CPSCF equations for the MM perturbations are not scaling with the MM subsystem size. As we showed, we can decouple this dependence by solving the Q-vector equations for N QM atoms and reconstruct the atomic charge derivatives with respect to MM perturbations by contracting with the external field derivatives. This induces an obvious saving in computational time. For the present case (plant cryptochrome, 9432 atoms at the B3LYP/TZVP/TZVPFit/Amber99 level using 32 cores of 2.6 GHz), the solution of conventional CPSCF equations for the full system would take around 86.5 hours, while the Qvector method (84 atom) just takes 17 minutes. This allows the construction of the full Hessian including the charge polarization terms at the cost of a mechanical embedding QM/MM full Hessian.

Normal mode analysis
We performed the normal mode analysis of the plant cryptochrome shown in Fig. 2 at the minimum energy structure of the QM/MM internal energy. This geometry was obtained from a quadratically convergent optimization of the QM subsystem, with a steepest descent optimization of the MM susbsystem at every step of optimization (microiterations). At the minimum structure, the FAD remains in a Ufolded conformation due to the electrostatic interactions with the protein environment. The protein model has 8121 atoms, from which 84 atoms are in the QM subsystem. This leads to a total of 24357 normal modes. In Fig. 3, we analyse the degree of delocalization as well as the QM character of these modes. First, we computed the participation ratio in order to measure the degree of delocalization. For this, we defined the participation ratio as in which 0 L i is the mass-weighted normal mode displacement of the ground state. For a fully delocalized mode, all atoms contribute equally to the mass-weighted displacement, and the participation ratio is equal to one, whereas if only few atoms contribute to the displacement, the participation ratio has the value of (3N ) −1 (which is close to 0 in the case of plant cryptochrome.) Following this analysis, we can observe that the low frequency modes (<200 cm −1 ) show on average the largest degree of delocalization, involving thus a large number of atoms in each normal mode frequency. These modes are usually involved in describing protein phase transitions. The rest of the normal modes involve only a few atoms. The participation ratio analysis is not sufficient to determine the degree of mixing between QM and MM atom motions. If one would perform a QM subsystem normal mode analysis, one would obtain 246 full QM normal modes. For the full protein diagonalization, these modes are strongly mixed. In Fig. 3, we report the QM character of each normal mode, defined as As we can observe, the low frequency modes below 200 cm −1 have a negligible QM character. This implies a strong mixing between the QM and MM atoms in this region. Since the number of MM atoms is so large, the QM character is "diluted" within the MM motion. Only two normal modes between 1200-1400 cm −1 keep a clear QM character (>90%), and another 23 normal modes have a QM character larger than 60%, comprised in the range of 400-1600 cm −1 . Therefore, only 10% of the FAD normal modes keep in essence its original character. For the rest, QM and MM atom motions are essentially mixed, even if the modes are not fully delocalized over the whole protein. Importantly, this mixing is present in the whole range of frequencies, and it is not only appearing for the fully delocalized low-frequency modes.

Vibronic couplings, Huang-Rhys factors and spectral densities
The Huang-Rhys couplings are a fundamental measure for the electron-nuclear (vibronic) coupling. The HR factor is directly related to the reorganization energy of one mode in the excited state. They can be defined as in which ω i is the normal mode frequency and ∆q i, is the projection in the ground state normal mode direction 0 L i of the Cartesian geometry difference between excited (R I ) and ground state (R 0 ) minima, mass weighted by the mass metric M = δ ij m i , in which m i is the atomic mass of atom i. The ∆q I i can be approximately obtained by the so-called vertical gradient approximation, in which the mass-weighted energy gradient of the excited state I (∆q I i ) is projected to the normal mode direction, leading to In order to define the effect of the bath modes, the notion of spectral density J(ω) is useful. It can be defined as Indeed, the spectral density contains information about each normal mode frequency, and its coupling to the system, which suffices to characterize the internal dissipation of an excited state due to nuclear motion. Using the information of the QM/MM analytic second derivatives together with the excited-state gradient at the ground-state minimum, we can construct the spectral densities characterizing the same cryptochrome model. In Figure 4, we have plotted such spectral density for the lowest energy bright excited state. Two spectral densities are shown, one corresponding to the normal modes and gradient of the full protein, and the other corresponding to the normal modes and gradient of the QM subsystem.
The overall contribution of MM atoms to the broadened spectral density is only minor. This is expected because the transition density is mainly located on the chromophore, which is entirely located in the QM subsystem, and thus only modes with substantial QM motion should be contributing to the spectral density. Still, an indirect contribution of the MM atoms appear in the spectral density due to the mode mixing, as shown for example in the most intense peak around 1620-1625 cm −1 . While the QM spectral density has only one peak of 40 meV, the full protein spectral density shows a spectral density contribution of 4 modes. This is due to the mixing of QM and MM atom motions, alike to the intensity borrowing mechanism of electronic spectra, in which dark states couple to bright states borrowing dipole intensity and thus becoming observable. This can be observed for almost all peaks in the spectral density, which reflects once again the strong normal mode mixed character that we showed in Fig. 3. It is important to notice though that even though there is a strong mixing, the spec-tral densities of the QM and the full protein are essentially equivalent. This is because the normal mode mixing occurs only between QM and MM modes at almost the same frequency. Eventually, the effect of mixing has only a minor impact on the frequency of the modes (less than 1 cm −1 for the most intense peak). This results justify the extraction of vibronic couplings from the normal mode analysis of the QM subsystem.

Vibrationally resolved absorption spectrum
Flavoproteins like cryptochrome act as bluelight sensors in biological systems. 43 These proteins have a well-resolved vibrational progression of the absorption spectrum, characteristic of the isoalloxazine chromophore. Unlike other types of proteins in which inhomogeneous broadening dominates, the homogeneous broadening has to be explicitly taken into account to reproduce the color of absorption of the flavin chromophore. 44,45 The approximate vibronic couplings shown in Eq. 26 can be used to compute the absorption spectrum including the vibrational structure due to homogeneous broadening.
Here, we compute the vibrational resolved spectrum of FAD in plant cryptochrome, by using the QM subsystem normal modes in conjunction with the vibrationally-resolved timedependent approach of Bloino and coworkers. 32 In Fig. 5, we compare the experimental and theoretical absorption spectra for arabidopsis thaliana plant cryptochrome with FAD in its fully oxidized form. 25 The experimental absorption spectrum is characterized by two main absorption peaks of π → π * character. The lowest bright transition has a peak maximum around 450 nm and the second bright state at around 360 nm. 25 These peaks have a marked vibrational structure, especially due to in-plane C-N stretchings. 46 Indeed, the first electronic transition is distributed in resolved peaks, approximately found at 470, 450 and 420 nm. Our simulated spectrum perfectly matches the relative position and intensities of the experimental peaks. From our simulations, we observe that these peaks are in fact formed from two electronic transitions, a HOMO→LUMO transition which give the basic vibrational structure, and a HOMO-1→LUMO transition which is centered at 420 nm. For the second absorption peak, we have a reasonable agreement between simulated and experimental spectra. Indeed, the relative peak position and the relative intensities between the two vibronic peaks are well reproduced by the theory. However, our theoretical simulations underestimate by around 40% the relative intensity with respect to the lowest energy absorption peak. We attribute this either to the approximations of the level of theory or some missing contributions of low-intensity high-energy peak. For the sake of comparison, both spectra have been normalized for the highest intensity peak between the range of 350 nm and 500 nm. The simulated spectrum has been red-shifted by 12 nm to match the experimental maximum. Experimental spectrum has been taken from Ref. 25

Conclusion
We have derived a fully analytic expression for the second derivative of the electrostatic embedding ESPF QM/MM energy. In this approach, the atomic charges representing the QM charge distribution are not only polarizable by the usual QM/MM electrostatic interactions, but also by the motion of the MM atoms. QM/MM electrostatic energy second derivative contains MM-block contributions too expensive to compute, and usually simply ignored or approximated. Using the ESPF method, we have shown that these terms can be equivalently computed with a set of auxiliary coupledperturbed self-consistent field equations that only scales with the QM subsystem size, 24 reducing thus the total computational cost of the full electrostatic embedding QM/MM hessian to a similar cost than mechanical embedding QM/MM hessian.
To show the capabilities of our method, we have computed the full hessian of arabidopsis thaliana plant cryptochrome. We have evidenced that most of the FAD chromophore (QM) vibrational modes are strongly mixed with the modes of the protein, and only a few high-frequency modes around 1200-1500 cm −1 keep a full QM character. In addition, we have computed the vibronic couplings and spectral densities for the bright excited state, together with the vibrationally resolved spectra of FAD in the protein. These applications open up the use of QM/MM hessians construct vibronic models for photoactive proteins and QM systems in complex environments.
Several limitations are currently imposed in our current model for a practical use. First, one of the main limitations of our QM/MM method is the lack of periodic boundary conditions to include solvent effects. In the future, we plan to implement our method in an extension of the ESPF QM/MM method to include Ewald summations. 47 Meantime, the only strategy is to freeze a state of the solvent, and relax the protein within a solvation shell. Since the QM equations do not scale with the MM subsystem size, the calculations considering a large number of solvent molecules should be tractable at a reasonable computational cost. Second, the current simulation of vibrationally resolved absorption spectrum requires a positive semi-definite Hessian matrix. This implies that any initial geometry must be optimized to a minimum energy structure so that all frequencies are positive. The geometry optimization step is time-consuming, hence severely limiting the usual efficiency associated with the inhomogenous broadening as it is frequently reported in the literature, i.e., by selection of snapshots from a MD simulation and performing the QM/MM analysis at each of them. This limitation could be overcome by using for example instantaneous normal modes, 48 although further investigation in this direction is required.
Finally, we plan to extend our model to address computed vibrational spectroscopies. Such simulations face frequently the problem of comparing with experimental spectra, which could ultimately be attributed to the lack of anharmonic effects, both from the QM and the MM frequency spectrum. To better match the QM frequencies to experimental data, it is common to scale the Hessian via a scaling factor that depends on the method and the basis set. Morokuma and coworkers applied a similar strategy to their ONIOM QM/MM Hessian, in which three different scaling factors for the QM, MM and QM-MM Hessian blocks are applied. 16 Such scheme would be readily applied to our analytic QM/MM Hessian too. In addition, some recent developments at the level of force fields for better describing infra-red spectra can also be trivially integrated in our Hessian scheme. 49,50

Appendix: Derivatives of the ESPF operator
The ESPF atomic charge operator elements Q A,µν are obtained by fitting the corresponding electrostatic potential integrals V k,µν , V k,µν = drχ µ (r)|r − R k | −1 χ ν (r) calculated over a grid of N grid points k, at position R k , around the QM atoms: The kernel T is a rectangular matrix of dimensions N grid ·N QM . The kernel definition depends on the multipolar expansion order. For example, for atomic charges it contains the inverse distance between QM atom and the grid point T A,k = |R A − R k | −1 . Thus, we can give an expression for the Q A,µν = χ µ | Q A |χ ν atomic charge operator elements by solving the previous set of linear equations, leading to The derivative of the ESPF hamiltonian involves the derivative of the atomic charge operator. A recursive formula for the n-derivative of the charge operator can be obtained by directly deriving n-times Eq. 29, in which Here we give the explicit first and second derivatives of the h operator defined in Eq. 5. The first derivative of the operator is given by The second term involves the derivative of the classical external field, which have simple expression given the form of the force field. The derivative of the atomic charge operator can be in term expanded (in matrix notation) as we can easily see that the computation of firstorder derivatives of h operators requires the computation of φ x A , V x k,µν and T x in addition to the zero-order field and charges.
The second derivatives of the h operator is The second term involves the second derivative of the classical external field, which again has a simple expression given the form of the force field, while the third and fourth term are constructed with elements of the first derivative of the operator. Hence, one needs only to develop the second derivative of the atomic charge operator. This is given by Therefore, in order to obtain the second derivative of the ESPF operator, we need to compute φ xy A , V xy k,µν and T xy in addition to the first derivatives and the zero-order quantities.