Volume 33 Issue 5
Nov.  2020
Turn off MathJax
Article Contents

Zhi-jun Zhang, Zi-fei Chen, Jian Liu. Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide†[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 613-622. doi: 10.1063/1674-0068/cjcp2006099
Citation: Zhi-jun Zhang, Zi-fei Chen, Jian Liu. Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 613-622. doi: 10.1063/1674-0068/cjcp2006099

Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide

doi: 10.1063/1674-0068/cjcp2006099
More Information
  • Corresponding author: Jian Liu, E-mail: jianliupku@pku.edu.cn
  • Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
  • Received Date: 2020-06-16
  • Accepted Date: 2020-07-13
  • Publish Date: 2020-10-27
  • Formaldehyde and hydrogen peroxide are two important realistic molecules in atmospheric chemistry. We implement path integral Liouville dynamics (PILD) to calculate the dipole-derivative autocorrelation function for obtaining the infrared spectrum. In comparison to exact vibrational frequencies, PILD faithfully captures most nuclear quantum effects in vibrational dynamics as temperature changes and as the isotopic substitution occurs.
  • Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
  • 加载中
  • [1] J. Liu and W. H. Miller, J. Chem. Phys. 126, 234110 (2007). doi:  10.1063/1.2743023
    [2] J. Liu and W. H. Miller, J. Chem. Phys. 128, 144511 (2008). doi:  10.1063/1.2889945
    [3] J. Liu and W. H. Miller, J. Chem. Phys. 131, 074113 (2009). doi:  10.1063/1.3202438
    [4] J. Liu and W. H. Miller, J. Chem. Phys. 134, 104102 (2011). doi:  10.1063/1.3555274
    [5] J. Liu and W. H. Miller, J. Chem. Phys. 134, 104101 (2011). doi:  10.1063/1.3555273
    [6] J. Liu, J. Chem. Phys. 134, 194110 (2011). doi:  10.1063/1.3589406
    [7] Q. Shi and E. Geva, J. Chem. Phys. 119, 9030 (2003). doi:  10.1063/1.1613636
    [8] A. Horikoshi and K. Kinugawa, J. Chem. Phys. 122, 174104 (2005). doi:  10.1063/1.1888576
    [9] J. Liu and Z. J. Zhang, J. Chem. Phys. 144, 034307 (2016). doi:  10.1063/1.4939953
    [10] J. P. Pinto, G. R. Gladstone, and Y. L. Yung, Science 210, 183 (1980). doi:  10.1126/science.210.4466.183
    [11] D. J. Luecken, S. L. Napelenok, M. Strum, R. Scheffe, and S. Phillips, Environ. Sci. Technol. 52, 4668 (2018). doi:  10.1021/acs.est.7b05509
    [12] Z. Li, A. N. Schwier, N. Sareen, and V. F. McNeill, Atmos. Chem. Phys. 11, 11617 (2011). doi:  10.5194/acp-11-11617-2011
    [13] T. Salthammer, S. Mentese, and R. Marutzky, Chem. Rev. 110, 2536 (2010). doi:  10.1021/cr800399g
    [14] T. Salthammer, Angew. Chem. Int. Ed. 52, 3320 (2013). doi:  10.1002/anie.201205984
    [15] L. Zhu, D. J. Jacob, F. N. Keutsch, L. J. Mickley, R. Scheffe, M. Strum, G. G. Abad, K. Chance, K. Yang, B. Rappenglück, D. B. Millet, M. Baasandorj, L. Jaeglé, and V. Shah, Environ. Sci. Technol. 51, 5650 (2017). doi:  10.1021/acs.est.7b01356
    [16] R. H. Hunt, R. A. Leacock, C. W. Peters, and K. T. Hecht, J. Chem. Phys. 42, 1931 (1965). doi:  10.1063/1.1696228
    [17] J. E. Carpenter and F. Weinhold, J. Phys. Chem. 90, 6405 (1986). doi:  10.1021/j100282a001
    [18] M. Praprotnik and D. Janežič, J. Chem. Phys. 122, 174102 (2005). doi:  10.1063/1.1884608
    [19] H. Romanowski, J M. Bowman, and L. B. Harding, J. Chem. Phys. 82, 4155 (1985). doi:  10.1063/1.448858
    [20] M. J. Bramley and N. C. Handy, J. Chem. Phys. 98, 1378 (1993). doi:  10.1063/1.464305
    [21] J. D. Rogers and J. J. Hillman, J. Chem. Phys. 75, 1085 (1981). doi:  10.1063/1.442181
    [22] R. Chen, G. Ma, and H. Guo, Chem. Phys. Lett. 320, 567 (2000). doi:  10.1016/S0009-2614(00)00254-2
    [23] R. Chen, G. Ma, and H. Guo, J. Chem. Phys. 114, 4763 (2001). doi:  10.1063/1.1348274
    [24] E. L. Klinting, D. Lauvergnat, and O. Christiansen, J. Chem. Theo. Comp. 16, 4505 (2020). doi:  10.1021/acs.jctc.0c00261
    [25] P. A. Giguére and T. K. K. Srinivasan, J. Raman Spectrosc. 2, 125 (1974). doi:  10.1002/jrs.1250020203
    [26] C. Camy-Peyret, J. M. Flaud, J. W. C. Johns, and M. Noël, J. Mol. Spectrosc. 155, 84 (1992). doi:  10.1016/0022-2852(92)90550-8
    [27] J. M. Flaud, C. Camy-Peyret, J. W. C. Johns, and B. Carli, J. Chem. Phys. 91, 1504 (1989). doi:  10.1063/1.457110
    [28] A. Perrin, A. Valentin, J. M. Flaud, C. Camypeyret, L. Schriver, A. Schriver, and P. Arcas, J. Mol. Spectrosc. 171, 358 (1995). doi:  10.1006/jmsp.1995.1125
    [29] A. Perrin, J. M. Flaud, C. Camy-Peyret, A. Goldman, F. J. Murcray, and R. D. Blatherwick, J. Mol. Spectrosc. 142, 129 (1990). doi:  10.1016/0022-2852(90)90297-4
    [30] J. J. Hillman, D. E. Jennings, W. B. Olson, and A. Goldman, J. Mol. Spectrosc. 117, 46 (1986). doi:  10.1016/0022-2852(86)90091-3
    [31] W. B. Olson, R. H. Hunt, B. W. Young, A. G. Maki and J. W. Brault, J. Mol. Spectrosc. 127, 12 (1988). doi:  10.1016/0022-2852(88)90004-5
    [32] W. B. Cook, R. H. Hunt, W. N. Shelton, and F. A. Flaherty, J. Mol. Spectrosc. 171, 91 (1995). doi:  10.1006/jmsp.1995.1104
    [33] D. Vione, V. Maurino, C. Minero, and E. Pelizzetti, Ann. Chim. 93, 477 (2003).
    [34] H. Sakugawa, I. R. Kaplan, W. Tsai, and Y. Cohen, Environ. Sci. Technol. 24, 1452 (1990). doi:  10.1021/es00080a002
    [35] G. W. Robinson and V. E. DiGiorgio, Can. J. Chem. 36, 31 (1958). doi:  10.1139/v58-004
    [36] J. T. Massey and R. W. Hart, J. Chem. Phys. 23, 942 (1955). doi:  10.1063/1.1742152
    [37] J. M. L. Martin, T. J. Lee, and P. R. Taylor, J. Mol. Spectrosc. 160, 105 (1993). doi:  10.1006/jmsp.1993.1161
    [38] J. Koput, S. Carter, and N. C. Handy, J. Phys. Chem. A 102, 6325 (1998). doi:  10.1021/jp9812583
    [39] M. Bonfanti, J. Petersen, P. Eisenbrandt, I. Burghardt, and E. Pollak, J. Chem. Theory Comput. 14, 5310 (2018). doi:  10.1021/acs.jctc.8b00355
    [40] J. Tatchen and E. Pollak, J. Chem. Phys. 130, 041103 (2009). doi:  10.1063/1.3074100
    [41] D. P. Chong and Y. Takahata, Chem. Phys. Lett. 418, 286 (2006). doi:  10.1016/j.cplett.2005.10.091
    [42] M. Mladenović, Spectrochim. Acta A Mol. Biomol. Spectrosc. 58, 809 (2002). doi:  10.1016/S1386-1425(01)00670-9
    [43] R. M. Levy, D. B. Kitchen, J. T. Blair, and K. Krogh-Jespersen, J. Phys. Chem. 94, 4470 (1990). doi:  10.1021/j100374a023
    [44] T. H. Dunning and V. McKoy, J. Chem. Phys. 48, 5263 (1968). doi:  10.1063/1.1668203
    [45] J. Koput, S. Carter, and N. C. Handy, J. Chem. Phys. 115, 8345 (2001). doi:  10.1063/1.1410976
    [46] D. V. Ilyin, W. A. Goddard Ⅲ, J. J. Oppenheim, and T. Cheng, Proc. Natl. Acad. Sci. USA 116, 18202 (2019). doi:  10.1073/pnas.1701383115
    [47] M. T. C. Martins-Costa and M. F. Ruiz-López, Chem. Phys. 332, 341 (2007). doi:  10.1016/j.chemphys.2006.12.018
    [48] J. Liu, W. H. Miller, G. S. Fanourgakis, S. S. Xantheas, S. Imoto, and S. Saito, J. Chem. Phys. 135, 244503 (2011). doi:  10.1063/1.3670960
    [49] J. Liu, J. Chem. Phys. 140, 224107 (2014). doi:  10.1063/1.4881518
    [50] E. L. Pollock and D. M. Ceperley, Phys. Rev. B 30, 2555 (1984). doi:  10.1103/PhysRevB.30.2555
    [51] M. E. Tuckerman, B. J. Berne, G. J. Martyna, and M. L. Klein, J. Chem. Phys. 99, 2796 (1993). doi:  10.1063/1.465188
    [52] J. Liu, D. Z. Li, and X. Liu, J. Chem. Phys. 145, 024103 (2016). doi:  10.1063/1.4954990
    [53] D. Z. Li, X. Han, Y. C. Chai, C. Wang, Z. J. Zhang, Z. F. Chen, J. Liu, and J. S. Shao, J. Chem. Phys. 147, 184104 (2017). doi:  10.1063/1.4996204
    [54] D. Z. Li, Z. F. Chen, Z. J. Zhang, and J. Liu, Chin. J. Chem. Phys. 30, 735 (2017). doi:  10.1063/1674-0068/30/cjcp1711223
    [55] Z. J. Zhang, X. Liu, Z. F. Chen, H. F. Zheng, K. Y. Yan, and J. Liu, J. Chem. Phys. 147, 034109 (2017). doi:  10.1063/1.4991621
    [56] Z. J. Zhang, K. Y. Yan, X. Liu, and J. Liu, Chin. Sci. Bull. 63, 3467 (2018). doi:  10.1360/N972018-00908
    [57] J. Liu, D. Z. Li, and X. Liu, Sci. Sin. Chim. 46, 27 (2016). doi:  10.1360/N032015-00143
    [58] X. Liu and J. Liu, J. Chem. Phys. 148, 102319 (2018). doi:  10.1063/1.5005059
    [59] H. B. Wang, X. Liu, and J. Liu, Chin. J. Chem. Phys. 31, 446 (2018). doi:  10.1063/1674-0068/31/cjcp1805122
    [60] L. Hascoet and V. Pascual, ACM Trans. Mathem. Softw. 39, 20 (2013).
    [61] J. W. Ochterski, http://gaussian.com/vib/ (1999).
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Figures(5)  / Tables(6)

Article Metrics

Article views(10) PDF downloads(1) Cited by()

Proportional views
Related

Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide

doi: 10.1063/1674-0068/cjcp2006099

Abstract: Formaldehyde and hydrogen peroxide are two important realistic molecules in atmospheric chemistry. We implement path integral Liouville dynamics (PILD) to calculate the dipole-derivative autocorrelation function for obtaining the infrared spectrum. In comparison to exact vibrational frequencies, PILD faithfully captures most nuclear quantum effects in vibrational dynamics as temperature changes and as the isotopic substitution occurs.

Part of the special issue for "the Chinese Chemical Society's 16th National Chemical Dynamics Symposium".
Zhi-jun Zhang, Zi-fei Chen, Jian Liu. Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide†[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 613-622. doi: 10.1063/1674-0068/cjcp2006099
Citation: Zhi-jun Zhang, Zi-fei Chen, Jian Liu. Path Integral Liouville Dynamics Simulations of Vibrational Spectra of Formaldehyde and Hydrogen Peroxide[J]. Chinese Journal of Chemical Physics , 2020, 33(5): 613-622. doi: 10.1063/1674-0068/cjcp2006099
  • It is challenging to include nuclear quantum effects (NQEs) for chemical dynamics of molecular systems. Most trajectory-based quantum dynamics methods fail to satisfy the two important properties simultaneously [1-8]: (ⅰ) conserve the quantum canonical distribution for a thermal equilibrium system, and (ⅱ) recover exact quantum correlation functions in the harmonic limit for both linear and nonlinear operators (i.e. linear or nonlinear functions of position or momentum operators). Recently developed path integral Liouville dynamics (PILD) is able to have these two properties, as well as yields exact correlation functions in the classical ($ \hbar \rightarrow $0) and high-temperature ($ \beta \rightarrow $0) limits. PILD has been shown to be a useful approach to compute vibrational spectra of realistic molecular systems, which leads to a reasonably accurate peak position with a relatively small full width at half maximum (FWHM) for an excited mode [9].

    In this work, we investigate the numerical performance of PILD for its applications to vibrational spectra for formaldehyde (HCHO) and hydrogen peroxide (H2O2), two realistic small but representative molecules that play important roles in atmospheric chemistry. Formaldehyde is the smallest organic molecule containing a carbonyl group. The photolysis of HCHO generates new radicals (OH and HO2) that drive ozone (O3) production [10, 11]. Formaldehyde has a strong influence on secondary organic aerosols (SOAs) by yielding radicals that increase gas-phase oxidation of hydrocarbons as well as by promoting formation of surface-active organic products that depress the surface tension [12]. Formaldehyde is a hazardous air pollutant and a kind of carcinogen [13-15]. In addition, formaldehyde is a chemical feedstock used in numerous industrial productions, e.g. methanol, wood-based materials, coatings, etc. [13].

    Hydrogen peroxide is a nonlinear and nonplanar molecule, a prototype for studying internal rotation [16-32]. H2O2 is one of the main final products of photochemical reactions in atmosphere [33]. It is often generated by the photolysis of formaldehyde as well as produced from the dismutation of the hydroperoxyl radical (HO2) formed from reactions between hydroxyl radicals and hydrocarbons [33]. Hydrogen peroxide is an efficient oxidant for SO2, which is responsible for the phenomenon of acid fog and acid rain [33, 34].

    Computational study on the two molecules has been performed since 1950s [16, 17, 19-24, 35-47]. Vibrational spectra produced by either conventional normal-mode analysis (NMA) or molecular dynamics (MD) are often significantly deviated from exact results (especially in the high frequency domain), because anharmonicity is strong in these two systems. As discussed in Ref.[48], MD explores only a small region around the equilibrium configuration due to the lack of the zero-point energy. This is why MD and NMA lead to similar results when the vibrational frequency is relatively high.

    The outline of the paper is as follows. Section II begins by reviewing the PILD approach for calculating the quantum correlation function. It then presents the PILD formulation for the dipole-derivative correlation function for computing the infrared (IR) spectrum. Section III describes the simulation details for the HCHO and H2O2 molecules. Section IV investigates the PILD spectrum as a function of temperature as well as isotopic substitution. Conclusion remarks and outlook follow in Section V.

  • Most dynamic quantities, such as spectra, rates and transport coefficients, can be described by quantum correlation functions, of which the form is

    where $ \hat A $ and $ \hat B $ are relevant operators, $ Z $ = $ {{\rm Tr}}\left({{ {\rm e}^{ - \beta \hat H}}} \right) $ is the partition function, $ \beta $ = $ (1/{k_ {\rm B}}T) $ with $ k_ {\rm B} $ being the Boltzmann constant and $ T $ being the temperature of the system, $ \hat A_{{{\rm Kubo}}}^\beta $ = $ {\frac{1}{\beta }}\int_0^\beta { {\rm d}\lambda {{\rm }}{ {\rm e}^{(\lambda - \beta)\hat H}}\hat A{ {\rm e}^{ - \lambda \hat H}}} $ for the Kubo-transformed version, or $ \hat A_{{{\rm std}}}^\beta $ = $ { {\rm e}^{ - \beta \hat H}}\hat A $ for the standard version of the correlation function. The time-independent Hamiltonian $ \hat H $ is assumed to be of standard Cartesian form for a general system with the total number of degrees of freedom $ N $,

    where $ {\bf{M}} $ is the diagonal "mass matrix" with elements $ \left\{ {{m_j}} \right\} $, $ V\left({{\bf{\hat x}}} \right) $ is the potential energy surface, and $ {\bf{\hat p}} $ and $ {\bf{\hat x}} $ are the momentum and coordinate operators, respectively.

    For instance, the Fourier transform of the (Kubo-transformed) collective dipole-derivative autocorrelation function

    with $ \hat {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}} $ the first-order derivative of the total dipole moment operator with respect to time, is directly connected to the experimental IR spectrum with the relation

    where $ \alpha (\omega) $ is the Beer-Lambert absorption constant, $ n(\omega) $ is the refractive index of the medium, $ \omega $ is the angular frequency, $ {\varepsilon _0} $ the vacuum permittivity, $ V $ the volume of the system, and $ c $ the speed of light. Calculation of $ {\langle {{\mathit{\boldsymbol{\hat {{\dot{{{\rm \mu }}} }}}}}(0){\hat {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}}}(t)} \rangle _{{{\rm Kubo}}}} $ is then the key to the atomic level simulation of the IR spectrum.

    PILD was derived from the phase space formulation of quantum mechanics [4-6] by employing imaginary time path integral to obtain the effective force [49]. The PILD expression of the quantum correlation function in Eq.(1) is

    where $ \rho _W^{{{\rm eq}}}\left({{\bf{x}}, {\bf{p}}} \right) $, $ {B_W}\left({{\bf{x}}, {\bf{p}}} \right) $ and $ f_{{A^\beta }}^W\left({{\bf{x}}, {\bf{p}}} \right) $ are defined as

    In Eq.(6) $ \rho _W^{{{\rm eq}}}\left({{\bf{x}}, {\bf{p}}} \right) $ is the density distribution function in thermal equilibrium in the Wigner space, which can be expressed with the Gaussian approximation of the momentum distribution [3, 6, 9, 49] as

    For the Kubo-transformed collective dipole-derivative auto-correlation function $ {\langle {\hat {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}}(0)\hat {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}}(t)} \rangle _{{{\rm Kubo}}}} $, the product $ f_{{A^\beta }}^W\left({{{\bf{x}}_0}, {{\bf{p}}_0}} \right){B_W}\left({{{\bf{x}}_t}, {{\bf{p}}_t}} \right) $ in Eq.(5) can be then obtained as

    Here the thermal mass matrix $ {{\bf{M}}_{{{\rm therm}}}} $ is defined in Appendix B, which is reduced to the (diagonal) physical mass matrix $ {\bf{M}} $ in the high-temperature or classical limit [49].

    In Eq.(5) the equations of motion of the trajectory $ \left({{{\bf{x}}_t}, {{\bf{p}}_t}} \right) $ are

    with the effective force $ -\frac{ {\partial V_{{{\rm eff}}}^{{{\rm PILD}}}\left({{\bf{x}}, {\bf{p}}} \right)}}{{\partial {\bf{x}}}} $ given by

    Path integral representation of $ \langle {\bf{x}}\left| {{{\rm{e}}^{ - \beta \hat H}}} \right.|{\bf{x}}\rangle $ is

    where $ {\bf{x}} \equiv {{\bf{x}}_1} \equiv {{\bf{x}}_{P + 1}} $ with $ P $ the number of path integral beads. After employing the staging transformation [50, 51] of the beads,

    Eq.(13) becomes

    where

    $ \omega _{{{\rm ad}}} $ is the adiabatic frequency, $ {\bf{\tilde M}}_j $ represents the fictitious masses corresponding to the staging variables $ \left({{{\mathit{\boldsymbol{\xi}}}_2}, \cdots, {{\mathit{\boldsymbol{\xi}}}_P}} \right) $, and the adiabatic parameter $ {\gamma _{{\rm{ad}}}} \in (0,\left. 1 \right]$ is chosen to separate the time scales for all the staging variables $ \left({{{\mathit{\boldsymbol{\xi}}}_2}, \cdots, {{\mathit{\boldsymbol{\xi}}}_P}} \right) $ from the one for $ {{\mathit{\boldsymbol{\xi}}}_1} \equiv {\bf{x}} $. The forces on the staging variables can be expressed recursively as

    It is easy to show that the effective force of Eq.(12) can be obtained from

    Inserting fictitious momenta $ \left({{{\bf{p}}_2}, \cdots, {{\bf{p}}_P}} \right) $ into Eq.(18) leads to

    Eq.(22) or the effective force Eq.(11) can be produced by sampling $ \left({{{\mathit{\boldsymbol{\xi}}}_2}, \cdots, {{\mathit{\boldsymbol{\xi}}}_P}, {{\bf{p}}_2}, \cdots, {{\bf{p}}_P}} \right) $ in a MD scheme. One can choose the adiabatic parameter $ {\gamma _{ {\rm ad}}} $ such that all the staging variables $ \left({{{\mathit{\boldsymbol{\xi}}}_2}, \cdots, {{\mathit{\boldsymbol{\xi}}}_P}} \right) $ share the same time scale that is well separated from the time scale of $ {{\mathit{\boldsymbol{\xi}}}_1} \equiv {\bf{x}} $. Then the PILD equations of motion (Eq.(11)) become

    While Eq.(22) suggests that the staging variables $ \left({{{\mathit{\boldsymbol{\xi}}}_2}, \cdots, {{\mathit{\boldsymbol{\xi}}}_P}, {{\bf{p}}_2}, \cdots, {{\bf{p}}_P}} \right) $ should be thermostatted for obtaining the effective force, the quantum phase space variables $ \left({{{\mathit{\boldsymbol{\xi}}}_1}, {{\bf{p}}_1}} \right) \equiv \left({{\bf{x}}, {\bf{p}}} \right) $ should not, because they account for real time dynamics. When Langevin dynamics is used as the thermostatting method, the PILD equations of motion are

    where $ {\gamma _{{{\rm Lang}}}} $ is the friction coefficient for Langevin dynamics.

    Because the unified "middle" thermostat scheme is efficient in sampling the configuration space [52-56], it is employed for the propagation of the PILD trajectory [9, 57]. When Langevin dynamics is used as the thermostat, the "BAOAB-num" algorithm for path integral molecular dynamics (PIMD) proposed in Refs. [52, 58, 59] is recommended. The integrator for Eq.(24) through one time step $ \Delta t $ can be expressed as

    Here, $ {{\mathit{\boldsymbol{\eta}}}_j} $ is a vector of random numbers with independent Gaussian distribution, and the coefficients $ {c_1} $ and $ {c_2} $ are

  • PIMD is used for equilibrating the molecular system and for evaluating the thermal mass matrix $ {{\bf{M}}_{{{\rm therm}}}} $ before PILD is employed for real time dynamics. The dipole moment $ {\mathit{\boldsymbol{{{\rm \mu }}}}} $ and its time derivative $ {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}} $ are obtained by the point charge model ($ {\mathit{\boldsymbol{{{\rm \mu }}}}} = {\sum\limits_{i = 1}^{{N_{ {\rm atom}}}} {{q_i}{{\mathit{\boldsymbol{x}}}_i}} ~~ {\rm{and}}~~ {\mathit{\boldsymbol{\dot{{{\rm \mu }}}}}}} = {\sum\limits_{i = 1}^{{N_{ {\rm atom}}}} {{q_i}{{{\mathit{\boldsymbol{\dot x}}}}_i}}} $, where $ {q_i} $, $ {{\mathit{\boldsymbol{x}}}_i} $, and $ {{\mathit{\boldsymbol{\dot x}}}_i} $ represent the point charge, position, and velocity of the $ i $-th atom, respectively). The Fourier transform of the Kubo-transformed collective dipole-derivative autocorrelation function is used to calculate the IR spectrum (i.e., Eq.(3) and Eq.(4)).

    We first apply PILD to simulate HCHO (formaldehyde) and its various isotope molecules with the accurate potential energy surface (PES) developed by Martin et al. [37]. As the explicit form of the PES is complicated, the code for the force and that for the Hessian of the PES are obtained by Tapenade [60] with automatic differentiation (AD) technique. $ P $ = 720 path integral beads are used in PIMD and PILD for $ T $ = 100 K, while $ P $ = 240 beads for $ T $ = 300 K. Converged results are obtained with the adiabatic parameter $ {\gamma _{ {\rm ad}}} $ = $ {10^{ - 4}} $. While the time interval for PIMD is $ ~ $0.024 fs, that for PILD is decreased to $ ~ $0.0012 fs in order to achieve $ ~ $1 cm-1 accuracy for the peak positions of the spectrum. The quantum correction factors are evaluated every 100 steps ($ ~ $2.4 fs) in the PIMD trajectory. 48 PIMD trajectories with each propagated up to $ ~ $480 ps are used for obtaining the averaged thermal mass matrix. Each PILD trajectory is propagated up to $ ~ $14.4 ps and 48 such trajectories (from different initial conditions) are collected for evaluating the dipole-derivative correlation function with the time averaging technique. The IR spectrum of HCHO is simulated for both $ T $ = 300 K and $ T $ = 100 K to test temperature effects. We then study the isotope molecules (HCDO and DCDO) for $ T $ = 300 K. "Exact" vibrational frequencies of these systems are available in Ref.[37] for comparison. The "exact" results were obtained using the vibrational self-consistent field-configuration interaction (VSCF-CI) method described in Ref.[19].

    We then investigate H216O2 and its isotope molecules. It is one of the simplest molecules involving a large-amplitude motion: internal rotation (torsion) around O-O bond. The molecule system presents a challenging test for both experimental and theoretical study. The accurate six dimensional PES was developed by Koput, Carter and Handy [38]. We use Tapenade [60] to generate the code for the force as well as that for the Hessian of the PES. The parameters used for converged results for the system are the same as those used in HCHO molecule. The evaluation of thermal mass matrix for H2O2 is different from that for HCHO, which is described in the Appendix B. We first study the spectrum of H2O2 at $ T $ = 300 K and that at $ T $ = 100 K. H218O2, HDO2, and D2O2 are then simulated for $ T $ = 300 K for studying isotope effects in the spectrum. The simulation results are compared to the "exact" vibrational frequencies in Refs.[38, 45].

    NMA data are obtained by diagonalization of the mass-weighted Hessian matrix at the equilibrium configuration (i.e., the minimum of the PES). Because vibrational frequencies produced by MD are similar to NMA results, we only compare PILD results to "exact" or NMA values in the work. The "exact" results were obtained using the fully symmetry-adapted finite-basis variational method described in Ref.[20].

  • FIG. 1 presents the PILD infrared spectra for the HCHO molecule at $ T $ = 300 K and $ T $ = 100 K. Table Ⅰ then depicts comparison of the exact vibrational frequencies to the spectra obtained by PILD at $ T $ = 100 K and 300 K and by the normal mode analysis (NMA). In comparison to the exact values of this molecular system, the NMA results are reasonable for low-frequency vibrational modes, but significantly blue-shifted for high-frequency modes. It indicates that nuclear quantum effects are remarkably strong for C-H bond stretching motions in the high-frequency region. The peak position of the symmetric stretching mode for the two C-H bonds is only about 50 cm-1 from that of the asymmetric stretching mode. The two peaks are readily distinguished in the PILD spectrum because the full width at half maximum (FWHM) of the PILD peak is about 10-15 cm-1. The peak positions of PILD spectrum are relatively insensitive to the temperature---no more than 5 cm-1 difference between the results for 300 K and for 100 K. This is in good agreement with the fact that exact vibrational frequencies of an isolated small molecule are independent of the temperature change.

    Figure 1.  The IR spectra for HCHO at both 300 K and 100 K. Purple dotted lines represent exact frequencies.

    Table Ⅰ.  Peak positions of the HCHO molecule at different temperatures. PES and exact results from Ref.[37] (unit: cm−1).

    FIG. 2 and Table Ⅱ then demonstrate the spectra and peak positions of the isotope molecules (HCDO and DCDO) for 300 K. All peaks are effectively sharp and their positions agree well with the exact values (the difference is no more than 11 cm-1). When the hydrogen atoms are substituted with the deuterium atoms in the formaldehyde molecule, the peak positions in HCDO and in DCDO are red-shifted from the corresponding ones in HCHO. The red-shift data are shown in Table Ⅲ. For instance, the red shifts for O-H stretching motions for the exact vibrational frequencies are 689 cm-1 in HCDO, and 725.9 cm-1, 684.2 cm-1 in DCDO, respectively. By comparison, the red shifts for the PILD results are 713 cm-1 in HCDO, and 737 cm-1, 693 cm-1 in DCDO. The differences between PILD and exact red shifts are small, no more than 24 cm-1. As a contrast, the NMA red shifts are 780.5 cm-1 in HCDO, and 795.8 cm-1, 758.2 cm-1 in DCDO, respectively. That is, NMA considerably overestimates the red shift, of which the deviation from the exact value is always greater than 70 cm-1. Comparison of the PILD peak positions of HCHO, HCDO, and DCDO to the exact vibrational frequencies in Table Ⅱ indicates that PILD robustly captures isotope effects for this system.

    Figure 2.  The IR spectra for isotope molecules of HCHO for 300 K. Purple dotted lines represent exact frequencies.

    Table Ⅱ.  Peak positions of isotope molecules of the HCHO molecule at 300 K. PES and exact results from Ref.[37] (unit: cm−1).

    Table Ⅲ.  Corresponding peak position red-shifts of the isotope molecules to the HCHO molecule at 300 K. Calculated based on Table Ⅰ and Table Ⅱ (unit: cm−1):

  • FIG. 3 shows the IR spectra obtained by PILD at $T = 300 $K and $T = 100$ K and Table Ⅳ lists the peak positions in comparison to exact and NMA results. In addition, experimental values are listed. Both PILD and NMA yield reasonable results that are close to exact vibrational frequencies in the low-frequency region. In the high-frequency region, PILD is much superior to NMA. While the error of the PILD peak position for the highest frequency is no more than ~25 cm-1 at 300 K and ~ 11 cm-1 at 100 K, that of the NMA result is as great as ~ 187 cm-1. Any PILD peak position for 100 K is no more than 16 cm-1 different from that for 300 K, which demonstrates that PILD is reasonably stable as the temperature changes.

    Figure 3.  The IR spectra for H2O2 at both 300 K and 100 K. Purple dotted lines represent exact frequencies.

    Table Ⅳ.  Peak positions of the H2O2 molecule at different temperatures. PES and exact results from Refs.[38, 45] (unit: cm−1). Tunneling splitting values for the PES are available in Ref.[45]. NMA captures no tunneling splitting effects. Experimental results are listed for comparison.

    FIG. 4 and Table Ⅴ demonstrate the results of the isotope molecules of H2O2 (HDO2, D2O2, and H218O2). All PILD peak positions show reasonable agreement with the exact results for these isotope molecules. The peak of the O-O stretching mode (the second fundamental, 864.5/884.9 cm-1 for H2O2, 825.7 cm-1 for H218O2, or 877.4 cm-1 for D2O2) is extremely weak in the infrared spectrum of H2O2, H218O2, or D2O2 (Refs.[21, 25, 26]. Since HDO2 breaks the symmetry of H2O2, one more distinct peak appears around 880 cm-1, as well as the two O-H/O-D bond stretching peaks become noticeably distinguishable. The isotopic substitution yields a significant red shift from the O-H stretching peak to the O-D one. Table Ⅵ lists the red-shift values of the peak positions of D2O2 in comparison to H2O2. PILD faithfully describes the red shift with an error no larger than 27 cm-1, while the difference between the NMA red shift and the exact value can be as large as 89 cm-1. The characteristic features in the spectrum of D2O2 are similar to those in the spectrum of H2O2. When the hydrogen atoms are replaced by the heavier deuterium atoms, the symmetry and asymmetry stretching modes (with highest frequencies) are considerably red-shifted from those of H2O2. The deviations of the PILD results from the exact red shift values are smaller than 3 cm-1, and in contrast, the error of the red shift caused by NMA is larger than 84 cm-1. Isotope effects of the substitution for the oxygen element are less significant for the high-frequency stretching vibrational modes.

    Figure 4.  The IR spectra for isotope molecules of H2O2 for 300 K. Purple dotted lines represent exact frequencies.

    Table Ⅴ.  Peak positions of isotope molecules of the H2O2 molecule at 300 K. PES and exact results from Refs.[38, 45] (unit: cm−1). The tunneling splitting value is available only for the first fundamental of H218O2, HDO2, or D2O2 in Ref.[45]. NMA captures no tunneling splitting effects.

    Table Ⅵ.  Corresponding peak position red-shifts of the D2O2 molecule from the H2O2 molecule at 300 K. Calculated based on Table Ⅳ and Table Ⅴ (unit: cm−1). Except for the first fundamental corresponding to the torsion mode, the exact red-shift value is obtained by using the average of the two tunneling splitting peak positions of the fundamental of the H2O2 molecule as the reference.

    It is worth pointing out that PILD fails to capture true long-time quantum coherence effects [9] in the correlation function, so it is not able to faithfully describe the tunneling splitting when its value is smaller than 25 cm-1. PILD semi-quantitatively captures the tunneling splitting of the torsion motion (during 200-400 cm-1) in the spectrum, which is relatively large and ranges from 39.3 to 120.3 cm-1 for the PES as shown in Tables IV and V. But PILD is not able to semi-quantitatively depict the tunneling splitting of any other fundamental as shown in Table Ⅳ. Neither is PILD capable of describing the difference between the exact vibrational frequencies around 3614-3632 cm-1 for H2O2 or those around 2676 cm-1 for D2O2 corresponding to the symmetric and asymmetric O-H (or O-D) stretching modes that are near-degenerate, e.g., the difference is less than 3 cm-1 for D2O2. This drawback of PILD can hardly be improved in the present framework, because no unique definition of trajectory in the phase space exists in quantum mechanics.

    Finally, the overtone bands of the torsion mode appear in the region of about 500-800 cm-1 in the PILD spectra. For example, the exact first overtone 776.2 cm-1 is reasonably accurately captured in the IR spectra of H2O2 in FIG. 3. As the intensity is relatively low, we do not explicitly discuss it in the paper.

  • As a trajectory-based approximate quantum dynamics approach (in the Wigner phase space), PILD has been demonstrated to be reasonably accurate for computing vibrational spectra for small molecule systems. The error of the PILD peak position is no more than ~ 42 cm-1 for the low-frequency region (500 cm-1), and no more than ~ 30 cm-1 for the relatively high frequency regions (< 500 cm-1) for all applications to formaldehyde, hydrogen peroxide and their isotope molecules. The FWHM of the PILD peak is relatively small (no larger than ~ 15 cm-1) such that the two highest-frequency modes with only ~ 50 cm-1 difference in the HCHO spectrum can be distinguished in the PILD simulation. The PILD vibrational frequencies for the two benchmark molecular systems are rather robust as the temperature changes. In summary, as long as true quantum coherence effects are not important, PILD performs well in studying the vibrational spectrum as a function of temperature and of isotopic substitution. In future work, it will be interesting to implement PILD for studying molecular vibrational spectra in different environments such as at interfaces and in solution.

  • This work was supported by the National Natural Science Foundation of China (No.21961142017), the Ministry of Science and Technology of China (No.2017YFA0204901 and No.2016YFC0202803), and the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase)(No.U1501501). We acknowledge the High-performance Computing Platform of Peking University and Beijng PARATERA Tech Co., Ltd. for providing computational resources.

  • In the standard normal-mode analysis, mass-weighted Hessian matrix elements are given by

    Where ${m_k}$ represents the mass of the $k$ -th degree of freedom with $N$ the total number of degrees of freedom ($N$ = $3{N_{ \rm atom }}$, $N_{ \rm atom }$ is the number of atoms). The mass-weighted coordinate vector is defined as ${{\bf{q}}}$ = ${{{\bf{M}}}^{1/2}}{{\bf{x}}}$.

    The Hamiltonian around ${{\bf{x}}}$ can be expanded to the 2nd order as

    The eigenvalues of the mass-weighted Hessian matrix produce normal-mode frequencies {${\omega _k}$}, i.e.,

    with ${{\bf{\Lambda}}}$ a diagonal matrix with the elements ${\lambda _k}$ = ${({\omega _k})^2}$ and ${{\bf{T}}}$ an orthogonal matrix. If ${{{\bf{b}}}_k}$ is the $k$ -th column of the matrix ${{\bf{T}}}$, then

    Rotational modes are not always easy to be distinguished from low-frequency vibrational modes (especially when imaginary frequencies are encountered) according to the eigenvalues of $\left\{ {{\lambda _k}} \right\}$ from directly diagonalization of the mass-weighted Hessian matrix. We adopt the vibrational analysis approach [61] to separate low-frequency vibrational modes from rotational ones.

    The essential idea is to construct a matrix ${{\bf{D}}}$ to block-diagonalize the mass-weighted Hessian matrix according to their modes. The three vectors $({{{\bf{D}}}_{{\bf{1}}}}, {{{\bf{D}}}_{{\bf{2}}}}, {{{\bf{D}}}_{{\bf{3}}}})$ of length $N$ corresponding to the translation of the center of mass are simply $\sqrt {{m_i}} $ times the corresponding coordinate axis. For example, for the water molecule (using ${m_{\rm H} }$ =1 and ${m_{\rm O}}$ =16) the translational vectors are

    Generating vectors corresponding to the rotational motion of the whole molecule is more complicated. The moment of inertia tensor ${\bf{I}}$ is defined as

    where ${{{\bf{r}}}_i}$ is the three-dimensional vector of the Cartesian coordinates shifted from the center of mass, for the $i$ -th atom. Diagonalize the moment of inertia tensor ${\bf{I}}$, and use the eigenvectors to form a 3 $\times$ 3 matrix ${\bf{W}}$, i.e.,

    where ${\bf{\Phi}}$ is a 3 $\times$ 3 diagonal matrix with its diagonal elements being the eigenvalues of ${\bf{I}}$.

    For the $i$ -th atom and the $\alpha$ -th dimension, (${P_{\alpha }}$)$_i$ is given by the dot product of ${{{\bf{r}}}_i}$ and the $\alpha$ -th column of ${\bf{W}}$,

    where $\alpha$ = $x, y, z$ represents the three spatial dimensions.

    The vectors for these rotational modes are then defined as

    where $\alpha = x, y, z$; $i$ = $1, \cdots, {N_{ \rm atom }}$.

    A Schmidt orthogonalization procedure is used to generate ${N_{ \rm vib }}$ = $N - 6$ (or $N - 5$ for linear molecule) remaining vectors (in matrix ${{\bf{D}}}$) for the vibrational modes of the system, which are orthogonal to the five or six rotational/translational vectors. Use the transformation matrix ${\bf{D}}$ to obtain the internal coordinates ${{\bf{S}}}$ = ${{\bf{Dq}}}$ from the mass weighted Cartesian coordinates ${{\bf{q}}} $ = ${{{\bf{M}}}^{1/2}}{{\bf{x}}}$, where the rotational and translational modes have been separated out.

    We then diagonalize submatrix ${\left({{{{\bf{D}}}^T}{{\bf{HD}}}} \right)_{ \rm vib }}$ corresponding to the vibrational modes, which is a ${N_{ \rm vib }}$ $ \times $ ${N_{ \rm vib }}$ matrix. i.e., the eigenvalues {${\lambda _k} $ = ${({\omega _k})^2}$} for vibrational modes are obtained from

  • Following the appendix of Ref.[9], below we describe the construction procedure of the thermal mass matrix ${{{\bf{M}}}_{{ \rm therm }}}$ defined in Eq.(9).

    The quantum correction factor $Q(u)$ with the local Gaussian approximation (LGA) ansatz [3] for both real and imaginary frequencies is as follows:

    for real $u$ ($u$ = $\beta \hbar \omega$)

    for imaginary $u$ ($u$ = $i{u_i}$)

    The frequency $\omega$ can be obtained using the standard normal mode analysis or the vibrational analysis (Eqs.(35)-(40)) in Appendix A.

    Consider the equilibrium configuration ${{{\bf{x}}}_0}$ of the molecule. ${{{\bf{b}}}_{0k}}$ is the $k$ -th column of the orthogonal matrix ${{{\bf{T}}}_0}$ for diagonalizing its Hessian matrix $H({{{\bf{x}}}_0})$. It is natural to use ${{{\bf{x}}}_0}$ as the reference configuration. That is, project the vector $\left\{ {{{{\bf{b}}}_j}({{\bf{x}}})} \right\}$ onto the vector $\left\{ {{{{\bf{b}}}_{0k}}({{{\bf{x}}}_0})} \right\}$ to define the average quantum correction factor matrix $\langle {{\bf{Q}}} \rangle $ as

    The bracket represents the canonical ensemble average over ${{\bf{x}}}$. The thermal mass matrix then becomes

    which is independent of the position ${{\bf{x}}}$. It is necessary to align the molecular configuration ${{\bf{x}}}$ and the equilibrium one ${{{\bf{x}}}_0}$ before evaluating the term ${\left| {{{\bf{b}}}_j^T({{\bf{x}}}) \cdot {{{\bf{b}}}_{0k}}({{{\bf{x}}}_0})} \right|^2}$ for the average quantum correction factor matrix $\langle {{\bf{Q}}} \rangle $. The Kabsch algorithm is employed as suggested in Ref.[9].

    As discussed in Eqs.(43) and (44), the thermal mass matrix depends on the evaluation of ${Q_j}({{\bf{x}}})$ and ${\left| {{{\bf{b}}}_j^T({{\bf{x}}}) \cdot {{{\bf{b}}}_{0k}}({{{\bf{x}}}_0})} \right|^2}$. The standard normal mode analysis and the vibrational analysis procedure (Eqs.(33)-(38)) yield nearly the same values of ${Q_j}({{\bf{x}}})$ for HCHO and its isotope molecules, because the vibrational modes are well separated from the rotational ones. In the paper, the standard normal mode analysis is used for the formaldehyde system.

    In the H2O2 molecule the torsional vibrational mode is entangled with the rotational modes, which are difficult to be distinguished by direct diagonalization of the mass-weighted Hessian matrix, i.e., the standard normal mode analysis. The vibrational analysis approach is then used to separate vibrational and rotational modes in the diagonalization of the mass-weighted Hessian matrix, though the separation is still not complete as it only works perfectly at the equilibrium configuration. In the evaluation of the thermal mass matrix for the hydrogen peroxide system, we use vibrational analysis to calculate frequencies lower than 1200 cm-1 and normal mode analysis for the others.

    The potential of the H2O2 molecule in the torsional dimension has a double well. That is, two equilibrium configurations with different H-O-O-H dihedral angles exist in the system. Both are used as the reference configurations in evaluation of the thermal mass matrix. When we use PIMD to sample the configuration, each sample is aligned to the two reference configurations, and we select the reference configuration that is more similar to the aligned sample (i.e., the one with the lower RMSD) for estimating the quantum correction factor. Since the low-frequency vibrational modes are not easy to be completely separated from the rotational modes, the overlapped element ${\left| {{{\bf{b}}}_j^T({{\bf{x}}}) \cdot {{{\bf{b}}}_{0k}}({{{\bf{x}}}_0})} \right|^2}$ is renormalized for the vibrational modes in order to cancel the error due to the coupling between the rotational modes and the vibrational ones.

Reference (61)

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return