Supporting information for: "State-resolved infrared spectrum of the
protonated water dimer: Revisiting the characteristic proton transfer
doublet peak"
Henrik R. Larsson,
1, 2
Markus Schröder,
3
Richard Beckmann,
4
Fabien Brieuc,
4,
a)
Christoph Schran,
4,
b)
Dominik
Marx,
4
and Oriol Vendrell
3
1)
Department of Chemistry and Biochemistry, University of California, Merced, CA 95343,
USA
2)
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125,
USA
3)
Theoretische Chemie, Physikalisch-Chemisches Institut, Universität Heidelberg, Im Neuenheimer Feld 229, D - 69120 Heidelberg,
Germany
4)
Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, 44780 Bochum, Germany
S1. SETUP
We use the same basis and polyspherical coordinates as
Vendrell et al.
1–3
. The names and basis size of the used co-
ordinates are shown in Table S1.
TABLE S1. Coordinates of the Zundel cation and their basis size.
Symbol basis size
meaning
푅
20 water-water stretch
푧
19
proton transfer
푥
7
푦
7
훼
11
water rotation
훾
푎
19
water wagging
훾
푏
19
훽
푎
9
water rocking
훽
푏
9
푟
1
푎
9
water
푎
푟
1
푏
9
휃
푎
9
푟
2
푎
9
water
푏
푟
2
푏
9
휃
푏
9
A. Tensor Tree
Both the tensor tree network state TTNS
4
and the PES repre-
sentations (see Section S4) uses the same mode combination,
that is, several coordinates are combined to “logical” coordi-
nates to decrease the dimensionality of the problem from
15
to
6
. Mode combination is used because the previous multicon-
figuration time-dependent Hartree (MCTDH) simulations
1–3
benefited from it. This is not necessarily the case for multi-
layer (ML) MCTDH and TTNSs but cumbersome attempts to
refit the PESs without mode combination have not been done.
a)
Current Address:Laboratoire Matière en Conditions Extrêmes, Université
Paris-Saclay, CEA, DAM, DIF, 91297 Arpajon, France
b)
Current Address:Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK
Due to the limit overall dimensionality of the system, we did
not use a sophisticated tree with many layers but rather use a
so-called Tucker matrix product state (MPS), as displayed in
Fig. S1. A Tucker MPS is a MPS with additional transforma-
tion matrices added to the legs of the MPS, shown in green
in Fig. S1. These transformation matrices are added to reduce
the overall size of each tensor, which is required because the
basis size of the logical coordinates is very large (up to
729
),
compared to the used bond dimension
퐷
(size here between
50
and
150
). Transformation matrices have not been added to the
ends of the MPS because there the tensors are two dimensional
and overall have small size (see Ref. [5] for a detailed discus-
sion). We optimized the order of the six logical coordinates
using the procedure from Ref. [4].
FIG. S1. Topology of the tensor tree used in the simulations. Dis-
played are the coordinates and the combined basis size.
퐷
푢
(
퐷
푙
)
is
the bond dimension of the upper (lower) tree.
B. Eigenstate optimization
We use a density matrix renormalization group-like proce-
dure to optimize the TTNS for each eigenstate. See Ref. [4] for
a detailed discussion of the used procedure. For each eigen-
state the bond dimensions
퐷
are adjusted based on a singular
value threshold
휖
, which we set to
10
−6
. In addition, we fixed
the maximal used bond dimension
퐷
max
to be either
50
,
70
,
100
, or
150
. For Fig. 1 in the main text, on the BBSM surface,
we used
퐷
max
= 100
for the lowest 383 states up to
1510cm
−1
,
followed by
퐷
max
=70
for the lowest 540 states up to
1665cm
−1
.
On the HBB surface, we used
퐷
max
= 150
for the lowest 150
Electronic
Supplementary
Material
(ESI)
for
Chemical
Science.
This
journal
is
©
The
Royal
Society
of
Chemistry
2022
2
states up to
1107cm
−1
, followed by
퐷
max
=100
for the lowest
249 states up to
1305cm
−1
and by
퐷
max
=50
for the lowest 892
states up to
1920cm
−1
.
To test the convergence, we compared the energies from the
optimizations with larger bond dimensionsions to those with
smaller bond dimensions. A comparison of the energy terms
up to
∼1030cm
−1
computed with a max. bond dimension of
either
100
or
150
leads to an error of less than
1
.
5cm
−1
, much
less than the errors of the PES representation (see Section S4).
A comparison of the energy terms up to
∼1300cm
−1
com-
puted with a max. bond dimension of either
70
or
100
leads to
an error of less than
8
.
0cm
−1
(
5cm
−1
) on the HBB (BBSM)
surface.
C. Time evolution
We generated the infrared spectrum obtained by time evo-
lution by Fourier transforming the convoluted autocorrela-
tion function of the dipole-operated ground state, as shown in
Ref. [1]. As damping function we chose
푔
(
푡
) = cos
2
(
휋푡
∕2
푇
)
.
The final propagation times
푇
were
2742fs
for the BBSM PES
and
2362fs
for the HBB PES.
The TTNSs were propagated using the time-depended
DMRG-based algorithm described in Refs. [6,7]. We gener-
ated the initial wavefunction
|
Ψ(
푡
= 0)
⟩
=
̂휇
푖
|
Ψ
0
⟩
by fitting
the action of one of the three dipole operators
̂휇
푥
, ̂휇
푦
, ̂휇
푧
to the
ground state
|
Ψ
0
⟩
using a DMRG-like procedure.
8
Unlike the
eigenstate optimizations, to simplify the simulation, we kept
bond dimensions for the MPS (
퐷
푢
) and for the transformation
matrices (
퐷
푙
) fixed at
퐷
푢
= 50
and
퐷
푙
= 80
, respectively.
S2. WAVEFUNCTIONS CUT
Plotting the two-dimensional wavefunction cuts requires de-
ciding which positions of the other coordinates need to be
fixed. As positions we chose the maxima of the diagonal of
the one-dimensional reduced density matrices (
Ψ
2
integrated
over all but one coordinate). In the cases where there were sev-
eral maxima, we analysed all possible wavefunction cuts and
chose the most representative one, typically that of the global
maximum of the reduced density matrix diagonal. The partic-
ular positions are listed in Table S2.
To reveal the excitations along the water-water stretch mo-
tion (
푅
), Fig. S2 shows additional cuts of the three states anal-
ysed in the main text.
S3. OVERLAP WITH TESTSTATES
Following the procedure introduced in Ref. [9] we confirm
our analysis of the decomposition of the eigenstates by com-
puting the overlap of the wavefunctions with zero-order states.
While there is some arbitrariness in defining these zero-order
states, they nevertheless provide a semi-quantitative confirma-
tion of the assignment. The teststates have been prepared sim-
ilar to Ref. [1].
|
1
푧
⟩
and
|
1
푅
⟩
(
|
04 − 40
⟩
) have been generated
TABLE S2. Positions of the cuts of the wavefunctions displayed
in this work. For all wavefunctions, the following positions are the
same:
푥,푦
= 0
,
훼
= 1
.
57
,
푟
1
푎
= 1
.
15
,
푟
2
푎
= 3
,
푟
1
푏
= 1
.
15
,
푟
2
푏
= 3
,
휃
푎
=
휃
푏
= 0
, and
훽
푎
=
훽
푏
= 0
. The angles are shown in radian.
Otherwise we use atomic units.
figure state cuts
2
Ψ
BBSM
푎
푅
= 4
.
78
, ,
훾
푏
= 0
Ψ
BBSM
푏
푅
= 4
.
67
,
훾
푏
= 0
Ψ
BBSM
푐
푅
= 4
.
56
,
훾
푏
= 0
Ψ
HBB
푎
푅
= 4
.
78
,
훾
푏
= 0
Ψ
HBB
푏
푅
= 4
.
56
,
훾
푏
= 0
Ψ
HBB
푐
푅
= 4
.
67
,
훾
푏
= 0
3
Ψ
BBSM
푏
푧
= 0
.
20
,
푅
= 4
.
67
Ψ
BBSM
푐
푧
= 0
,
푅
= 4
.
56
Ψ
HBB
푏
푧
= −0
.
10
,
푅
= 4
.
56
Ψ
HBB
푐
푧
= 0
.
15
,
푅
= 4
.
67
S2
Ψ
BBSM
푎
푧
= 0
.
15
,
훾
푏
= 0
Ψ
BBSM
푏
푧
= −0
.
20
,
훾
푏
= 0
Ψ
BBSM
푐
푧
= 0
,
훾
푏
= 0
Ψ
HBB
푎
푧
= 0
,
훾
푏
= 0
Ψ
HBB
푏
푧
= 0
.
10
,
훾
푏
= 0
Ψ
HBB
푐
푧
= −0
.
15
,
훾
푏
= 0
TABLE S3. Overlap of zero-order states with the three wavefunctions
Ψ
푖
around the doublet.
name & PES
̃
퐸
∕cm
−1
⟨
1
푅,
02 − 20
|
Ψ
⟩
2
⟨
1
푧
|
Ψ
⟩
2
⟨
04 − 40
|
Ψ
⟩
2
Ψ
푎
HBB
897
0.34
0.09
0.00
Ψ
푎
BBSM
920
0.46
0.19
0.00
Ψ
푏
HBB
1043
0.04
0.20
0.30
Ψ
푏
BBSM
1041
0.20
0.36
0.00
Ψ
푐
HBB
1060
0.11
0.21
0.13
Ψ
푐
BBSM
1095
0.02
0.00
0.44
from the 15-dimensional Hamiltonian with all but the
푧
and
푅
(
훾
푎
and
훾
푏
coordinates fixed at the maxima of the ground-
state).
|
02 − 20;1
푅
⟩
was generated by multiplying the cut of
|
1
푅
⟩
in
푅
onto the
|
02 − 20
⟩
state of the full Hamiltonian. For
all states up to a transition of
1310cm
−1
, the zero-order states
only had significant overlap with the ones mentioned in the
main text.
Overlaps of the zero-order states with the wavefunctions
considered here are shown in Table S3. The table confirms
the previous, wavefunction-based analysis. While Fig. S2 dis-
plays a small contribution from
|
1
푧
⟩
in
Ψ
HBB
푐
, there is no over-
lap between the
|
1
푧
⟩
zero-order state and
Ψ
HBB
푐
. Although this
could be attributed to the quality and arbitrariness of the zero-
order state, a closer inspection of
Ψ
HBB
푐
does not reveal a sig-
nificant excitations along the
1
푧
modes, except for the small
contribution displayed in Fig. S2. Interestingly, compared to
Ψ
HBB
푎
and
Ψ
HBB
푏
, the
|
1
푧
⟩
state is more dominant in
Ψ
BBSM
푎
and
Ψ
BBSM
푏
. This can be attributed to the missing contribu-
tions from
|
04 − 40
⟩
for these states.
3
FIG. S2. Representative cuts of the wavefunctions corresponding to either the doublet and a satellite peak (BBSM, upper panels), or the triplet
(HBB PES, lower panels) in the IR spectrum. Panel
(
푖
)
corresponds to wavefunction
Ψ
푖
, see text for details. The abscissa shows the wagging
(pyramidalization) motion of one of the water molecules. The ordinate shows the water-water stretch motion. The red lines denote the zero
contours.
S4. FITTING OF THE POTENTIAL ENERGY SURFACE
AND DIPOLE MOMENTS IN SUM-OF-PRODUCTS FORM
The multi-dimensional potential and dipole moment sur-
faces are expressed in the form of a canonical polyadic
decomposition
10
as
푉
(
푞
1
,
...
,푞
푓
) =
푅
∑
푟
=1
푐
푟
∏
휅
휗
푟,휅
(
푞
휅
)
,
S1
where
푞
휅
are (logical) coordinates,
푐
푟
are expansion coeffi-
cients and
휗
are basis functions. The basis functions are not
required to be orthogonal, which allows for a very compact
representation of the surfaces.
The generation of the decomposition is performed using a
Monte-Carlo variant (MCCPD) of the alternating least squares
(ALS) algorithm as outlined in Ref. [11].
A. HBB surface
The PES from Ref. [12] has been fitted with
푅
= 2048
terms. All symmetries have been implemented as in Ref. [11].
The PES differs from Ref. [11] by the use of more terms and
thus a slightly higher accuracy. TTNS-based reoptimization
of the first 239 states on the PES from Ref. [11] leads to dif-
ferences in the frequencies of max.
3
.
5cm
−1
.
The sampling points for the Monte-Carlo integration within
the ALS algorithm were created with a Metropolis algorithm
such they resemble a Boltzmann distribution
푊
(
푞
1
,
...
,푞
푓
) = exp
(
훽푉
(
푞
1
,
...
,푞
푓
)
)
S2
where
훽
= 1∕
푘
B
푇
may be interpreted as an inverse tempera-
ture
푇
with
푘
B
being the Boltzmann constant.
In the present case three different sets of sampling points,
each containing 10
6
points, have been created with a Metropo-
lis algorithm at inverse temperatures of
1∕
훽
=
1000 cm
−1
,
1∕
훽
=
2000 cm
−1
and
1∕
훽
=
4000 cm
−1
. These sets have been
subsequently combined into one single set for creating the fit.
The validation of the fit was performed using independent
path integral molecular dynamics (PIMD) configurations gen-
erated with the same settings as in Ref. [13]; see the next Sec-
tion S4 B for details. The mean and root-mean-square (RMS)
errors obtained with the validation sets for different tempera-
tures are outlined in Table S4.
All three dipole moment surfaces,
휇
푥
,
휇
푦
, and
휇
푧
, have been
fitted with
푅
= 1024
terms using the same sampling points as
for generating the PES fit. The obtained errors using the same
validation sets as for the PES fit are given in Table S5.
B. BBSM surface
The PES
14
and dipole moments
15
have been fitted using
PIMD configurations from Ref. [14] as sampling points. The
sets of PIMD configurations have been created at different tem-
peratures from 1.6K to 300K as outlined in Tables S6 and S7.