We present a numerically exact approach for evaluating vibrationally resolved electronic spectra at finite temperatures using the coherence thermofield dynamics. In this method, which avoids implementing an algorithm for solving the von Neumann equation for coherence, the thermal vibrational ensemble is first mapped to a pure-state wavepacket in an augmented space, and this wavepacket is then propagated by solving the standard, zero-temperature Schrödinger equation with the split-operator Fourier method. We show that the finite-temperature spectra obtained with the coherence thermofield dynamics in a Morse potential agree exactly with those computed by Boltzmann-averaging the spectra of individual vibrational levels. Because the split-operator thermofield dynamics on a full tensor-product grid is restricted to low-dimensional systems, we briefly discuss how the accessible dimensionality can be increased by various techniques developed for the zero-temperature split-operator Fourier method.

Vibrationally resolved electronic spectroscopy has proven to be a powerful tool for probing molecular systems and understanding photochemical processes. These experiments are typically conducted at ambient or higher temperatures, which are relevant to understanding biochemical and atmospheric processes as well as developing new materials.1–4 Other than the appearance of hot bands and the thermal Doppler broadening effects, non-Condon effects can also noticeably alter spectral features at finite temperatures.5–8 These thermal effects may become more pronounced in two-dimensional spectroscopy experiments.9 Including temperature effects in theoretical simulations is thus essential to accurately interpret and predict the results of vibronic spectroscopy.3,4,9–16

Given the large number of vibrational levels that may be occupied, time-independent sum-over-state expressions become impractical at higher temperatures for more than a few vibrational degrees of freedom. Time-dependent approaches offer a promising alternative, in which the spectrum is evaluated as the Fourier transform of an appropriate autocorrelation function without the need to select relevant transitions and compute their Franck–Condon factors.11,17 However, the most widely used quantum mechanical18–22 and semiclassical17,23–28 methods were developed only for pure states and cannot directly handle a thermal ensemble of states.

Standard approaches to include temperature effects either directly evolve the thermal density matrix or employ statistical averaging over pure states. An alternative way to represent quantum dynamics at finite temperatures, called thermofield dynamics,29,30 maps the thermal ensemble to a pure state in an augmented space, which can be propagated by conventional wavepacket methods for solving the time-dependent Schrödinger equation. Thermofield dynamics has been applied to problems in electronic structure31 and molecular quantum dynamics,32–38 including the simulation of vibronic spectra.9,12,39–43 In a similar spirit, Grossmann performed exact calculations of the thermal dipole time correlation function for vibrational spectra by propagating two thermal wavepackets and evaluating the matrix element of the dipole moment between them.44 

In vibronic spectroscopy applications, the evolution of the thermofield wavepacket replaces the dynamics of the coherence between two electronic states. Due to the doubling of dimensions, the coherence thermofield dynamics obviously benefits from a combination with semiclassical or other approximate propagation methods that scale favorably with the number of degrees of freedom. Combining the thermofield dynamics with the extended thawed Gaussian approximation45–47 has enabled the simultaneous description of finite-temperature, non-Condon,40 and even anharmonicity effects.9,12,41 Yet, one would sometimes like to avoid dynamical approximations and obtain the exact quantum solution.

For small systems with a known potential energy surface, the dynamical Fourier method may be used to solve the time-dependent Schrödinger equation on a grid with the split-operator algorithm.18,22,26,48,49 Although limited to low dimensions, the split-operator Fourier method offers a numerically exact solution, which can validate semiclassical and other more approximate approaches.

Here, we combine the coherence thermofield dynamics with the split-operator Fourier method in order to simulate vibronic spectra at finite temperatures. This approach is validated on a one-dimensional anharmonic Morse model by comparing the spectra computed with the split-operator coherence thermofield dynamics to the results of the Boltzmann averaging of spectral contributions from different initial vibrational levels. We also examine the thermofield wavepacket’s evolution in the augmented configuration space consisting of both the physical and fictitious degrees of freedom.

Let us consider a molecular system with D nuclear degrees of freedom and two uncoupled adiabatic electronic states: the ground state |g⟩ and the excited state |e⟩. The absorption spectrum for vibronic transitions from |g⟩ to |e⟩ is obtained as the Fourier transform,
σ(ω)=constC(t)eiωtdt,
(1)
of the dipole time correlation function,
C(t)=Tr(μ̂eiĤet/μ̂ρ̂eiĤgt/),
(2)
where μ̂ represents the transition dipole moment operator, Ĥg and Ĥe are the vibrational Hamiltonians associated with the ground and excited electronic states, ρ̂ is the vibrational density operator of the ground electronic state, and const = 2πω/(ℏc).17,22 As we are interested only in the dependence of the spectrum on frequency and temperature, in the following, we set const = 1/(2π).
At a finite temperature T, the density of occupied vibrational levels in the initial population follows the Boltzmann distribution,
ρ̂=eβĤg/Tr(eβĤg),
(3)
with the inverse temperature β = 1/(kBT). We can then rewrite C(t) as
C(t)=npneiωn,gtn,g|μ̂eiĤet/μ̂|n,g,
(4)
where pn=eβωn,g/keβωk,g is the Boltzmann probability, |n, g⟩ denotes the nth vibrational eigenfunction of Ĥg, and ωn,g is the corresponding vibrational energy.11 The thermal spectrum can be obtained from
C(t)=npneiωn,gtφn(0)|φn(t)
(5)
by propagating the initial vibrational wavepackets |φn(0)=μ̂|n,g individually as |φn(t)=eiĤet/|φn(0). In high-dimensional systems, propagating wavepackets associated with all occupied vibrational levels becomes impractical, and a statistical sampling approach is needed.15,16,50–52
Using the coherence thermofield dynamics, the thermal ensemble is, instead, mapped to a pure state in a doubled configuration space. The thermal dipole time correlation function (2) can be rewritten exactly12 (see the  Appendix) as a wavepacket autocorrelation function,
C(t)=ψ̄|eiH̄̂t/|ψ̄,
(6)
of an initial thermofield wavepacket, represented in the augmented thermofield space (denoted by a bar) as
ψ̄(q̄)=q|μ̂ρ̂1/2|q,
(7)
where q̄=(q,q)T is a 2D-dimensional coordinate vector in the augmented configuration space and |q′⟩ represents the position states in the “fictitious” Hilbert space (denoted by a prime). The wavepacket is propagated according to the time-dependent Schrödinger equation with the “augmented” Hamiltonian H̄̂ expressed in the position representation as
H̄(q̄)=He(q)Hg(q),
(8)
with a 2D × 2D mass matrix,
m̄=m00m,
(9)
replacing the D × D mass matrix m. In contrast to standard thermofield dynamics, where the two Hamiltonians in the right-hand side of Eq. (8) are identical, the thermofield Hamiltonian (8) in vibronic spectroscopy depends on the two Hamiltonians associated with the two electronic states involved in the transition. In essence, the thermofield wavepacket encodes the coherence between the dynamics on the two potential energy surfaces.
If the ground surface
Vg(q)=v0+(qqref)Tκ(qqref)/2
(10)
is harmonic, the initial thermofield wavepacket12,
ψ(q̄)=e(i/)[(q̄q̄0)TĀ0(q̄q̄0)/2+p̄0T(q̄q̄0)+γ̄0]
(11)
is a Gaussian centered at
q̄0=qrefqref,p̄0=00
(12)
in the phase space. This Gaussian has a temperature-dependent width matrix,12 
Ā0=iABBA,
(13)
where A and B are the D × D block matrices,
A=m1/2Ωcoth(βΩ/2)m1/2,
(14)
B=m1/2Ωsinh(βΩ/2)1m1/2,
(15)
with Ω=m1/2κm1/2. The phase factor,
γ̄0=(i/2)ln[det(mΩ/π)],
(16)
serves to ensure normalization.

The spectrum at a finite temperature can now be obtained by applying standard techniques for solving the time-dependent Schrödinger equation to the thermofield wavepacket with a doubled number of coordinates. Although the increase in dimensionality implies that an exact thermofield solution does not save computational time compared to the direct propagation of the density matrix, the wavepacket approach permits the use of zero-temperature numerical methods and avoids implementing methods for solving the von Neumann equation.

Among various numerical methods for propagating wavepackets, we choose to combine the thermofield dynamics with the second-order split-operator algorithm,18,22,26,48,49 a popular choice for separable Hamiltonians, because it is explicit, numerically stable, and easy to implement and because it preserves several geometric properties of the exact evolution operator, such as unitarity, symplecticity, and time reversibility.

In the following, we use natural units with m==kB=μ̂=1. As in Ref. 12, the excited-state surface is a one-dimensional Morse potential,
Ve(q)=Ve(qref)+ωe4χ1e2mωeχ(qqref)2,
(17)
with the equilibrium position qref = 1.5, the minimum energy Ve(qref) = 10, the fundamental frequency ωe = 0.9 at qref, and the anharmonicity measure χ = 0.02. The initial, ground-state potential energy surface is harmonic and given by Eq. (10) with qref=0 and κ = 1.

The spectra are computed either by Boltzmann averaging [Eq. (4)] or with thermofield dynamics [Eq. (6)] at scaled temperatures Tω = kBT/(ℏωg) = 0, 0.5, 1, and 2, where ωg=κ/m is the ground-state vibrational frequency. Here, as a demonstration, we apply the standard second-order TVT algorithm on a grid with 256 points for each degree of freedom. The wavepacket is propagated for a total time of 1000 with a time step of 0.1. A Gaussian damping function with a half-width at half-maximum of 15 is applied to the autocorrelation functions before the spectra are evaluated with Eq. (1).

As shown in Fig. 1, the thermofield split-operator Fourier method gives the same results as the Boltzmann-weighted average. As the temperature increases, the spectral range broadens and additional peaks (the “hot bands”) appear. In Ref. 12, a comparison was made between the thermal spectra obtained using the thermofield thawed Gaussian approximation and the “exact” spectra obtained by computing individual Franck–Condon factors by numerical integration. With the time-dependent split-operator approach, we avoid the need to preselect the relevant transitions and can give numerically exact results for arbitrary global potentials, even if their vibrational eigenfunctions do not have analytical expressions. This approach thus has the potential to validate semiclassical methods in a wider range of situations, including those with a fitted global potential energy surface.

FIG. 1.

Comparison of spectra obtained via coherence thermofield dynamics (black solid line) and Boltzmann averaging (red dashed line) at four different scaled temperatures Tω.

FIG. 1.

Comparison of spectra obtained via coherence thermofield dynamics (black solid line) and Boltzmann averaging (red dashed line) at four different scaled temperatures Tω.

Close modal

In Fig. 2, we show the decomposition of the thermal spectrum at Tω = 1 into the contributions from different initial vibrational levels. Only spectra from levels n = 0, 1, 2, and 3 are shown since the contributions from n > 3 are negligible here (the Boltzmann weight is 0.0002 for n = 4). As expected, the higher vibrational levels are the source of hot bands and expand the spectral range. Whereas the Boltzmann averaging requires propagating each initial vibrational level separately, a single thermofield wavepacket trajectory takes into account all thermal contributions.

FIG. 2.

Contributions from different vibrational levels, scaled by their Boltzmann probability of occupation, to the Boltzmann-averaged spectrum at Tω = 1.

FIG. 2.

Contributions from different vibrational levels, scaled by their Boltzmann probability of occupation, to the Boltzmann-averaged spectrum at Tω = 1.

Close modal

With a grid-based method, we can easily visualize the wavepacket and its evolution. Figure 3 shows the initial thermofield wavepacket at different temperatures in the augmented coordinate space. As the temperature increases, the off-diagonal block B of the initial width matrix (13) increases and the initial thermofield wavepacket rotates and stretches as a result.

FIG. 3.

Initial thermofield wavepackets at different scaled temperatures Tω; q is the physical degree of freedom, whereas q′ is fictitious.

FIG. 3.

Initial thermofield wavepackets at different scaled temperatures Tω; q is the physical degree of freedom, whereas q′ is fictitious.

Close modal

In Fig. 4, we show the evolution of the wavepacket at Tω = 1. Whereas the extent of the wavepacket over the fictitious degree of freedom q′ remains approximately constant, the initial excitation of the coordinate q on the anharmonic excited-state surface leads to a significant evolution of the wavepacket in the physical dimension.

FIG. 4.

Evolution of the thermofield wavepacket (for Tω = 1) in the augmented configuration space at different times t; q is the physical degree of freedom, whereas q′ is fictitious.

FIG. 4.

Evolution of the thermofield wavepacket (for Tω = 1) in the augmented configuration space at different times t; q is the physical degree of freedom, whereas q′ is fictitious.

Close modal

We have described a very simple method for computing vibronic spectra at finite temperatures. This method applies the Fourier split-operator propagation algorithm to the thermofield wavepacket representation of the electronic coherence, avoiding both Boltzmann averaging and solving the von Neumann equation. Instead, we solve the zero-temperature Schrödinger equation in an augmented space. The method yields quantum mechanically exact spectra of arbitrary global potentials, making it a valuable tool for validating semiclassical and other approximate techniques for finite-temperature spectra calculations. In contrast to Boltzmann averaging spectral contributions from different initial vibrational levels, the thermofield dynamics approach obtains the exact result for a given temperature from a single propagation. At high temperatures, as the number of contributing states increases, the thermofield approach becomes more favorable. If highly accurate results are desired, the simple approach described here can be made both more accurate and efficient by implementing high-order geometric integrators obtained by symmetric composition of the elementary second-order algorithm, in the same way as for the zero-temperature split-operator calculations.53–56 

The method presented here is based on the full tensor-product grid and, therefore, is limited to systems with very few degrees of freedom. To extend its applicability to higher-dimensional systems, one can combine the split-operator thermofield dynamics with adaptive moving grids,57–59 matching-pursuit coherent-state basis,60 tensor-train representation,61 or any other technique designed for the zero-temperature case. In fact, tensor-train representations have already been used extensively32,36 for the thermofield dynamics, although not in conjunction with the split-operator Fourier method. Here, we have intentionally restrained ourselves to the simplest, full-grid method, because we did not want to obscure the main idea by combining it with various strategies for increasing the accessible dimensionality. We believe that even the full-grid version can be interesting for researchers who would like to quickly check the thermal effects on their zero-temperature spectra of low-dimensional systems without implementing new codes.

The method is not restricted to linear spectra. In the same way as the Gaussian,9 multiconfigurational Ehrenfest,42 or optimized mean-trajectory43 semiclassical thermofield dynamics, the method can be applied to more sophisticated experimental techniques, such as pump–probe and two-dimensional spectroscopies. Beyond vibronic spectra, this approach may also find applications in other finite-temperature dynamics problems, such as the calculation of the rates of internal conversion.27,37,62,63

The authors thank Tomislav Begušić for discussions and acknowledge the financial support from the EPFL.

The authors have no conflicts to disclose.

Zhan Tong Zhang: Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Jiří J. L. Vaníček: Conceptualization (lead); Formal analysis (equal); Methodology (equal); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

The correlation function (2) can be re-expressed as12 
C(t)=Tr(ρ̂1/2μ̂eiĤet/μ̂ρ̂1/2eiĤgt/)=dqdqq|ρ̂1/2μ̂|qq|eiĤet/μ̂ρ̂1/2eiĤgt/|q=dq̄ψ̄0(q̄)*ψ̄t(q̄),
(A1)
where ψ̄t(q̄) is the propagated thermofield wavefunction,
ψ̄t(q̄)=q|eiĤet/μ̂ρ̂1/2eiĤgt/|q=eiĤe(q)t/eiĤg(q)t/q|μ̂ρ̂1/2|q=eiH̄(q̄)t/ψ̄0(q̄),
(A2)
q̄=(q,q)T is a 2D-dimensional coordinate vector, and H̄(q̄) is the augmented Hamiltonian (8).
1.
E. J. G.
Peterman
,
T.
Pullerits
,
R.
Van Grondelle
, and
H.
Van Amerongen
,
J. Phys. Chem. B
101
,
4448
(
1997
).
2.
A.
Gaiduk
,
M.
Yorulmaz
,
P. V.
Ruijgrok
, and
M.
Orrit
,
Science
330
,
353
(
2010
).
3.
Y.
Qian
,
X.
Li
,
A. R.
Harutyunyan
,
G.
Chen
,
Y.
Rao
, and
H.
Chen
,
J. Phys. Chem. A
124
,
9156
(
2020
).
4.
S.
Ndengué
,
E.
Quintas-Sánchez
,
R.
Dawes
,
C. C.
Blackstone
, and
D. L.
Osborn
,
J. Phys. Chem. A
127
,
6051
(
2023
).
6.
Y.
Tanimura
and
S.
Mukamel
,
J. Opt. Soc. Am. B
10
,
2263
(
1993
).
7.
S.
Kundu
,
P. P.
Roy
,
G. R.
Fleming
, and
N.
Makri
,
J. Phys. Chem. B
126
,
2899
(
2022
).
8.
P. P.
Roy
,
S.
Kundu
,
N.
Makri
, and
G. R.
Fleming
,
J. Phys. Chem. Lett.
13
,
7413
(
2022
).
9.
T.
Begušić
and
J.
Vaníček
,
J. Phys. Chem. Lett.
12
,
2997
(
2021
).
10.
G.
Zucchelli
,
R. C.
Jennings
,
F. M.
Garlaschi
,
G.
Cinque
,
R.
Bassi
, and
O.
Cremonesi
,
Biophys. J.
82
,
378
(
2002
).
11.
A.
Baiardi
,
J.
Bloino
, and
V.
Barone
,
J. Chem. Theory Comput.
9
,
4097
(
2013
).
12.
T.
Begušić
and
J.
Vaníček
,
J. Chem. Phys.
153
,
024105
(
2020
).
13.
J.
Xun
,
J.
Deng
, and
R.
He
,
Chem. Phys. Lett.
759
,
138045
(
2020
).
14.
Š.
Sršeň
,
J.
Sita
,
P.
Slavíček
,
V.
Ladányi
, and
D.
Heger
,
J. Chem. Theory Comput.
16
,
6428
(
2020
).
15.
A. M.
Wernbacher
and
L.
González
,
Phys. Chem. Chem. Phys.
23
,
17724
(
2021
).
16.
A.
Prlj
,
E.
Marsili
,
L.
Hutton
,
D.
Hollas
,
D.
Shchepanovska
,
D. R.
Glowacki
,
P.
Slavíček
, and
B. F. E.
Curchod
,
ACS Earth Space Chem.
6
,
207
(
2022
).
17.
E. J.
Heller
,
Acc. Chem. Res.
14
,
368
(
1981
).
18.
M. D.
Feit
,
J. A.
Fleck
, Jr.
, and
A.
Steiger
,
J. Comput. Phys.
47
,
412
(
1982
).
19.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
,
Chem. Phys. Lett.
165
,
73
(
1990
).
20.
M.
Ben-Nun
and
T. J.
Martínez
,
J. Phys. Chem. A
103
,
10517
(
1999
).
21.
G. A.
Worth
,
M. A.
Robb
, and
I.
Burghardt
,
Faraday Discuss.
127
,
307
(
2004
).
22.
D. J.
Tannor
,
Introduction to Quantum Mechanics: A Time-dependent Perspective
(
University Science Books
,
Sausalito
,
2007
).
23.
E. J.
Heller
,
J. Chem. Phys.
75
,
2923
(
1981
).
24.
M. F.
Herman
and
E.
Kluk
,
Chem. Phys.
91
,
27
(
1984
).
25.
26.
C.
Lubich
,
From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis
, 12th ed. (
European Mathematical Society
,
Zürich
,
2008
).
27.
E. J.
Heller
,
The Semiclassical Way to Dynamics and Spectroscopy
(
Princeton University Press
,
Princeton, NJ
,
2018
).
28.
M.
Bonfanti
,
J.
Petersen
,
P.
Eisenbrandt
,
I.
Burghardt
, and
E.
Pollak
,
J. Chem. Theory Comput.
14
,
5310
(
2018
).
29.
M.
Suzuki
,
J. Phys. Soc. Jpn.
54
,
4483
(
1985
).
30.
Y.
Takahashi
and
H.
Umezawa
,
Int. J. Mod. Phys. B
10
,
1755
(
1996
).
31.
G.
Harsha
,
T. M.
Henderson
, and
G. E.
Scuseria
,
J. Chem. Theory Comput.
15
,
6127
(
2019
).
32.
R.
Borrelli
and
M. F.
Gelin
,
J. Chem. Phys.
145
,
224101
(
2016
).
33.
R.
Borrelli
and
M. F.
Gelin
,
Sci. Rep.
7
,
9127
(
2017
).
34.
M. F.
Gelin
and
R.
Borrelli
,
Ann. Phys.
529
,
1700200
(
2017
).
35.
L.
Chen
and
Y.
Zhao
,
J. Chem. Phys.
147
,
214102
(
2017
).
36.
R.
Borrelli
and
M. F.
Gelin
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
11
,
e1539
(
2021
).
37.
E. W.
Fischer
and
P.
Saalfrank
,
J. Chem. Phys.
155
,
134109
(
2021
).
38.
M. F.
Gelin
and
R.
Borrelli
,
J. Chem. Theory Comput.
19
,
6402
(
2023
).
39.
C. S.
Reddy
and
M. D.
Prasad
,
Mol. Phys.
113
,
3023
(
2015
).
40.
C. S.
Reddy
and
M. D.
Prasad
,
J. Phys. Chem. A
120
,
2583
(
2016
).
41.
T.
Begušić
and
J.
Vaníček
,
Chimia
75
,
261
(
2021
).
42.
L.
Chen
,
R.
Borrelli
,
D. V.
Shalashilin
,
Y.
Zhao
, and
M. F.
Gelin
,
J. Chem. Theory Comput.
17
,
4359
(
2021
).
43.
K.
Polley
and
R. F.
Loring
,
J. Chem. Phys.
156
,
124108
(
2022
).
44.
F.
Grossmann
,
J. Chem. Phys.
141
,
144305
(
2014
).
45.
S.-Y.
Lee
and
E. J.
Heller
,
J. Chem. Phys.
76
,
3035
(
1982
).
46.
A.
Patoz
,
T.
Begušić
, and
J.
Vaníček
,
J. Phys. Chem. Lett.
9
,
2367
(
2018
).
47.
A.
Prlj
,
T.
Begušić
,
Z. T.
Zhang
,
G. C.
Fish
,
M.
Wehrle
,
T.
Zimmermann
,
S.
Choi
,
J.
Roulet
,
J.-E.
Moser
, and
J.
Vaníček
,
J. Chem. Theory Comput.
16
,
2617
(
2020
).
48.
D.
Kosloff
and
R.
Kosloff
,
J. Comput. Phys.
52
,
35
(
1983
).
49.
R.
Kosloff
and
D.
Kosloff
,
J. Chem. Phys.
79
,
1823
(
1983
).
50.
U.
Manthe
and
F.
Huarte-Larrañaga
,
Chem. Phys. Lett.
349
,
321
(
2001
).
51.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, 1st ed., edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
Wiley VCH
,
Weinheim
,
2009
), https://onlinelibrary.wiley.com/doi/book/10.1002/9783527627400
52.
M.
Barbatti
and
K.
Sen
,
Int. J. Quantum Chem.
116
,
762
(
2016
).
53.
M.
Wehrle
,
M.
Šulc
, and
J.
Vaníček
,
Chimia
65
,
334
(
2011
).
54.
S.
Choi
and
J.
Vaníček
,
J. Chem. Phys.
150
,
204112
(
2019
).
55.
J.
Roulet
,
S.
Choi
, and
J.
Vaníček
,
J. Chem. Phys.
150
,
204113
(
2019
).
56.
J.
Roulet
and
J.
Vaníček
,
J. Chem. Phys.
155
,
204109
(
2021
).
57.
J. F.
Thompson
,
Z. U.
Warsi
, and
C. W.
Mastin
,
Numerical Grid Generation: Foundations and Applications
(
North-Holland
,
Amsterdam
,
1997
), Vol.
45
, Chap. 11.
58.
L. R.
Pettey
and
R. E.
Wyatt
,
Chem. Phys. Lett.
424
,
443
(
2006
).
59.
S.
Choi
and
J.
Vaníček
,
J. Chem. Phys.
151
,
234102
(
2019
).
60.
X.
Chen
and
V. S.
Batista
,
J. Chem. Phys.
125
,
124313
(
2006
).
61.
S. M.
Greene
and
V. S.
Batista
,
J. Chem. Theory Comput.
13
,
4034
(
2017
).
62.
Y.
Wang
,
J.
Ren
, and
Z.
Shuai
,
J. Chem. Phys.
154
,
214109
(
2021
).
63.
M.
Wenzel
and
R.
Mitric
,
J. Chem. Phys.
158
,
034105
(
2023
).
Published open access through an agreement with EPFL