RESEARCH ARTICLE | APRIL 22 2020

TRAVIS—A free analyzer for trajectories from molecular

simulation
Special Collection: Classical Molecular Dynamics (MD) Simulations: Codes, Algorithms, Force fields, and Applications ,
Chemical Physics Software Collection

M. Brehm  ; M. Thomas; S. Gehrke ; B. Kirchner

J. Chem. Phys. 152, 164105 (2020)

https://doi.org/10.1063/5.0005078

CrossMark

 
View Export
Online Citation

Articles You May Be Interested In

Observation of out-of-plane ambient noise on two vector sensor moorings in Lake Travis
J. Acoust. Soc. Am. (October 2019)

A Beamline Matching application based on open source software

AIP Conference Proceedings (August 2001)

Computer simulation of the role of torsional flexibility on mass and momentum transport for a series of

01 November 2023 13:58:33

linear alkanes
J. Chem. Phys. (August 2012)
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

TRAVIS—A free analyzer for trajectories from

molecular simulation
Cite as: J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078
Submitted: 17 February 2020 • Accepted: 23 March 2020 •
Published Online: 22 April 2020

M. Brehm,1,a) M. Thomas,1 S. Gehrke,2 and B. Kirchner2

AFFILIATIONS
1
Institut für Chemie, Martin-Luther-Universität Halle–Wittenberg, von-Danckelmann-Platz 4, D-06120 Halle (Saale), Germany
2
Mulliken Center for Theoretical Chemistry, Rheinische Friedrich-Wilhelms-Universität Bonn, Beringstr. 4+6,
D-53115 Bonn, Germany

Note: This paper is part of the JCP Special Topic on Classical Molecular Dynamics (MD) Simulations: Codes, Algorithms,
Force fields, and Applications.
a)
Author to whom correspondence should be addressed: [email protected]. URL: https://brehm-research.de/

ABSTRACT
TRAVIS (“Trajectory Analyzer and Visualizer”) is a program package for post-processing and analyzing trajectories from molecular dynamics

01 November 2023 13:58:33

and Monte Carlo simulations, mostly focused on molecular condensed phase systems. It is an open source free software licensed under the
GNU GPL, is platform independent, and does not require any external libraries. Nine years after the original publication of TRAVIS, we
highlight some of the recent new functions and features in this article. At the same time, we shortly present some of the underlying algorithms
in TRAVIS, which contribute to make trajectory analysis more efficient. Some modern visualization techniques such as Sankey diagrams are
also demonstrated. Many analysis functions are implemented, covering structural analyses, dynamical analyses, and functions for predicting
vibrational spectra from molecular dynamics simulations. While some of the analyses are known since several decades, others are very recent.
For example, TRAVIS has been used to compute the first ab initio predictions in the literature of bulk phase vibrational circular dichroism
spectra, bulk phase Raman optical activity spectra, and bulk phase resonance Raman spectra within the last few years.
© 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license
(http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0005078., s

I. INTRODUCTION It is a commonly known phenomenon that many working

groups from the field of computational science use their own in-
Within the last few decades, the methods of molecular dynam- house scripts and codes to evaluate and analyze simulation trajec-
ics (MD) and Monte Carlo simulations have become a central part tories. These codes all have their specific scopes, strengths, and lim-
of computational chemistry and physics. In contrast to static calcu- itations, and despite being developed in independent groups, many
lations, they implicitly include the effect of entropy at finite temper- methods and algorithms are for sure almost equivalent in all these
atures. Most experimentally accessible observables are only defined tools. This comes with several disadvantages. On the one hand, it
as ensemble averages, rendering it indispensable to rely on ensem- wastes a lot of valuable scientific manpower by “reinventing the
ble averaging also in computations if experimental results shall be wheel” hundreds of times. On the other hand, it increases the risk
reproduced or predicted. The direct result of a simulation is a tra- of subtle scientific or coding problems not being detected because
jectory that contains the positions and velocities of the particles the user base of the code is often only one working group and there-
in equidistant discrete time steps. Mathematically speaking, such a fore rather small. Hence, it is highly desirable to have common tools
trajectory is a path through a 6n-dimensional space if n particles for analyzing trajectories, which are utilized by a large group of
have been simulated. Such a high-dimensional object is not suitable scientists.
for direct evaluation—methods for reducing the dimensionality are In 2011, we published the TRAVIS (“Trajectory Analyzer and
required. These methods will be termed trajectory analyses here. Visualizer”) program package1 as an open source free software,

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-1

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

which specifically aims at providing as many trajectory analyses as image files, but rather input files for visualization programs such
possible within a unified program so that many scientists are hope- as VMD,2 xmgrace,8 Gnuplot,9 Povray,10 and Mathematica.11 By
fully no longer forced to develop their own in-house solutions to supplying these files to the corresponding programs, high-quality
analyze trajectories. There already exist several other software pack- images can be easily obtained without requiring much knowledge
ages for analyzing trajectories, e.g., the Visual Molecular Dynamics of the visualization program.
(VMD) program package2 or Python frameworks such as MDAnal-
ysis.3 However, by its huge number of available analyses and by A. Fast molecule recognition via domain
its interactive command-line user interface, TRAVIS is in a certain decomposition
sense unique. Since the release in 2011, the original TRAVIS article1
TRAVIS only reads atom labels, positions, and optionally
has already been cited more than 400 times by more than 800 differ-
velocities from the input trajectory; any topological information
ent co-authors (according to the journal homepage), and around 95%
(molecules, bonds, etc.) is ignored. Instead, a distance-based bond
of these publications actually used the software to produce results
detection is executed in one trajectory frame, and molecules are
in the manuscript. This gives us the impression that our code is
defined on the basis of the detected covalent bonds. A pair of atoms
somehow useful for the community.
is defined to be covalently bonded if its distance is smaller than
In this article, we want to highlight some of the recently added
the sum of the covalent radii of the two atom types (taken from
methods and features in TRAVIS. At the same time, we are going
the literature)12 multiplied with a constant to account for thermal
to shortly introduce some of the underlying algorithms that are
vibrations (by default 1.15, but can be chosen differently). To detect
typically invisible for users but are, of course, very important to
bonds, a straight-forward (“canonical”) approach of checking the
achieve the desired efficiency and accuracy of trajectory analyses.
distance of all pairs of atoms in the system is very inefficient, as its
This manuscript is structured as follows: After a discussion of gen-
computational time scales with O(n2 ). TRAVIS uses a more sophis-
eral features of TRAVIS in Sec. II, we will discuss specific features
ticated approach, which is based on a spatial domain decomposition.
from the fields of structural, dynamical, and spectroscopic analy-
Domain decompositions and similar approaches such as the linked
ses in Secs. III–V. In this context, structural analyses are defined as
list method13 are commonly used in computational chemistry (in
analyses that are invariant under any re-shuffling of the trajectory
particular, in parallelized force field molecular dynamics codes such as
frames, while dynamical analyses somehow depend on the ordering
LAMMPS5 and Gromacs14 ) but have rarely been applied to trajectory
of the frames. Of course, spectroscopic analyses are also dynamical
post-processing. First, the longest possible bond distance is deter-

01 November 2023 13:58:33

analyses, but they are discussed separately. We end our article with
mined by finding the element in the system with the largest covalent
some conclusions.
radius, and taking this radius twice, multiplied by the constant men-
tioned above. Then, the simulation cell is tessellated into cubes of
equal size and shape. The number of cubes along each cell vector is
II. GENERAL FEATURES chosen as large as possible while still satisfying the constraint that
the edge length of the cube is larger than the longest possible bond
Being licensed under the GNU GPL v3 license, TRAVIS distance. All atoms in the system are then binned into these cubes.
is an open source free software and can be downloaded from To find all covalent bond partners of a certain atom inside of such
http://www.travis-analyzer.de/. a cube, only atoms within the same cube and within the 26 neigh-
Apart from the source code, Windows executables can also boring cubes need to be considered. This yields a very fast bond
be obtained from the homepage. New versions, including bug fixes detection algorithm that essentially scales with O(n).
and additional features, are regularly supplied. TRAVIS is written in The efficiency of this approach is easily visible in Table I,
C++, and the developers’ version currently contains around 270 000 where the molecule recognition times for similar systems (con-
lines of source code. It does not require any external libraries or taining ion pairs of the ionic liquid [EMIm][OAc], i.e., 1-ethyl-3-
dependencies (except from a C++ compiler) and is platform inde- methylimidazolium acetate) of increasing size are shown. While the
pendent. TRAVIS is currently not parallelized, i.e., it runs on one fourth column gives the molecule recognition time with the canon-
single central processing unit (CPU) core. However, we are work- ical algorithm, the fifth column depicts the timing of the domain
ing on introducing optional OpenMP-based parallelization for some decomposition approach. For relatively small systems, both methods
suitable modules. are acceptable. However, for a system with around 500 000 atoms in
TRAVIS reads many different input trajectory formats,
currently including XYZ, PDB, mol2, DL_POLY,4 LAMMPS,5
Amber,6 Gaussian Cube, and compressed BQB trajectories.7 Non-
TABLE I. Execution times of the molecule recognition for different system sizes with
orthorhombic cell geometries and variable cell vectors (NpT ensem- canonical implementation (t Can ) and domain decomposition implementation (t Dom ).
ble) are supported. The cell geometry can directly be read from The snapshots contain ion pairs of the ionic liquid [EMIm][OAc].
the trajectory if present. Very large systems with several hundred
thousand of atoms can be easily handled. Ion pairs Atoms Cell (nm) t Can (s) t Dom (s)
TRAVIS can be used as an interactive program, asking ques-
tions to the user in a text-mode console interface. Most of the ques- 36 936 2.121 0.003 0.001
tions are more or less self-explanatory, and reasonable default values 288 7 488 4.242 1.036 0.001
are provided in most cases. On the other hand, TRAVIS can also 2 304 59 904 8.485 83.009 0.011
be used for batch processing with the help of input files. When 18 432 479 232 16.969 3240.000 1.002
creating plots and visualizations, TRAVIS does not directly output

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-2

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

a 17 nm cubic cell, the canonical approach takes almost 1 h for the which is, of course, not the case with multi-purpose compression
molecule recognition, while the implementation based on domain formats. The high compression efficiency is achieved by applying
decomposition only runs for 1 s. a sequence of methods, including multi-dimensional polynomial
After the molecule recognition, topological atom codes are extrapolation, Hilbert curve re-ordering, and canonical multi-table
assigned to the atoms by using Shelley’s extension15 of the Mor- Huffman encoding.7
gan algorithm.16 The algorithm guarantees to give identical codes to While our original aim was to compress volumetric data trajec-
topologically equivalent atoms and different codes otherwise. Atom tories, we noted that the same methodology is also able to losslessly
numbers are assigned by descending atom codes. As a by-product, compress standard position trajectories. When starting from a XYZ
TRAVIS maintains a list of topologically identical atoms, which is trajectory file with a frame interval of 0.5 fs, the BQB format reaches
helpful for many analyses. Furthermore, a full ring system detection a compression ratio of around 20:1 here. This is much more effi-
and refinement17 is carried out, as some of the analyses are well- cient than other existing trajectory formats such as DCD or Gromacs
suited for cyclic structures. Finally, all individual molecules with XTC14 while not sacrificing accuracy. More details can be found in
identical sets of atom codes are considered to be topologically equiv- Ref. 7.
alent and grouped together as molecule types. We would like to point The BQB format and the methods described above are imple-
out again that this full protocol only takes 1 s on one CPU core for a mented in TRAVIS, which can be used to compress or decompress
system with around 500 000 atoms and 40 000 molecules, as shown volumetric or position trajectories, or directly read and analyze such
in Table I. trajectories. Apart from that, there is also a stand-alone tool for
compressing and decompressing trajectories, called the “bqbtool.”
B. The BQB format—Lossless compression A library (“libbqb”) to add support for the BQB format to other pro-
of trajectories and volumetric data gram packages will be released soon. For more information, refer to
the website of the BQB format.18
In 2018, we developed and published an approach for a highly
efficient lossless compression of volumetric data trajectories7 and
introduced the BQB file format to store the compressed data (“loss-
less” indeed means that the input Gaussian Cube file is bit-wise repro- III. STRUCTURAL ANALYSES
duced). The methods in TRAVIS for computing vibrational spectra In this section, some structural analyses implemented in
(see below) require the total electron density on a grid along the tra- TRAVIS will be presented. As mentioned above, structural analyses

01 November 2023 13:58:33

jectory, easily amounting to several terabytes of raw data. With our are defined as analyses that are invariant under any re-shuffling of
novel compression scheme, these data can be compressed signifi- the trajectory frames in this context. Due to this property, structural
cantly without losing any information. Figure 1 gives some details analyses can also be applied to Monte Carlo simulations because they
on the performance of the method. While multi-purpose compres- do not require that a physical time elapses along the trajectory.
sion tools such as bzip2 and xz fail to exploit the particular structure
of the data (smoothness in both space and time) and only reach a
compression ratio of around 5:1, the BQB format yields a lossless A. Particle density histograms
compression ratio of around 40:1 while still being faster than the
former two methods in compression and decompression time. Fur- Probably, the most common kinds of trajectory analyses are
thermore, the BQB format allows for random access to the frames, particle density histograms. In particular, the radial distribution
function (RDF), also known as g(r), is found in countless publica-
tions on the subject of liquid phase simulation. This was also one of
the first analyses that has been implemented in TRAVIS. Apart from
the RDF itself, TRAVIS also writes out the corresponding number
integral curve so that, e.g., the number of molecules within the first
solvation shell can be easily determined. Recently, we added a func-
tion that automatically recognizes the first maximum and the first
minimum of the RDF and prints out the distance as well as the inten-
sity of these two points in the log file together with the particle num-
ber found within the first shell. This is useful for batch processing if
many different trajectories are evaluated. Another recent addition is
an algorithm for computing RDFs up to larger distances than half of
the cell diameter,19 which can save some computer time. RDFs com-
puted with TRAVIS have appeared many times in the literature,20–32
including number integrals and coordination numbers.33–38
Apart from RDFs, TRAVIS can also create other one-
dimensional particle density histograms, e.g., involving angles and
dihedral angles. The corresponding functions are then called angu-
lar distribution functions (ADFs) and dihedral distribution func-
FIG. 1. Compression times (orange), decompression times (green), and com- tions (DDFs), respectively. Furthermore, the distance between a
pressed sizes (blue) of a volumetric electron density frame with different lossless
point and a line can be monitored, where the point can be any atom
compression algorithms. The BQB format7 discussed here is shown on the right.
in some molecule and the line can be defined by any other two atoms

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-3

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

in the same or a different molecule. The corresponding analysis is

called point–line distribution function (LiDF). Similarly, the dis-
tance between a point and a plane can be observed, where the plane
can be defined either by three points in some molecules or by spec-
ifying a base point and a normal vector. This analysis is termed
point–plane distribution function (PlDF). For all these functions,
the temporal development of individual pairs can also be observed.
In the case of RDFs, this would yield a plot with selected (or all) pair
distances against the simulation time to give an example.
If the simulation is not in thermal equilibrium, some of the
histograms might change over simulation time. TRAVIS offers to
create a 2D contour plot with the histogram quantity on one axis and
the simulation time on the other axis. This analysis, which we call
time-dependent distribution function (TDDF), investigates if some
processes (phase separation or crystallization) take place during a
simulation. Currently, only RDFs are supported for TDDFs, but the
other quantities will be implemented soon.
Another frequently used class of particle density histograms are
spatial distribution functions (SDFs). Here, a reference molecule
is fixed by defining a local coordinate system based on three points
from the molecule. Afterward, the average particle density of some
other molecules around the fixed reference molecule is computed. FIG. 2. Combined distribution function (CDF) depicting the O⋯H distance and the
As the resulting particle density is a volumetric dataset, only iso- O⋯H–C angle in the ionic liquid [EMIm][OAc].23,24 A possible geometric hydrogen
surfaces of these data can be visualized. TRAVIS offers to output bond criterion is shown by the yellow rectangle.
the resulting volumetric data as a particle density (nm−3 ) or rel-
ative to the uniform density, similar to RDFs. Furthermore, there

01 November 2023 13:58:33

are some advanced functions implemented such as smoothing SDFs, indicating a hydrogen bond, at low distances and angles of around
inverting SDFs, creating cut planes, and forcing SDFs to be symmet- 180○ . The plot gives the information that at small distances, the
ric to certain planes to improve sampling. There are many exam- angle is preferentially found to be close to 180○ . This is a correlation
ples where SDFs computed by TRAVIS have been published in the between both quantities, and the standard 1D histograms would not
literature.20,21,24,25,32,39–48 be able to give this information.
If the simulation cell is non-isotropic (e.g., systems with a sur- TRAVIS combines a multitude of scalar quantities into multi-
face or lipid bi-layer), it can be interesting to compute density pro- dimensional histograms, which we call combined distribution
files (DProf) along a fixed vector in the simulation cell. The result is functions (CDFs) here. By using combinations of distances, angles,
a histogram that gives the particle density of a selected particle type dihedral angles, point–plane distances, and point–line distances,
(either in nm−3 or relative to uniform density) in thin slices of the sys- thousands of interesting CDFs can be computed from a single
tem perpendicular to the chosen vector. Density profiles computed trajectory. There are many examples in the literature of CDFs
by TRAVIS have already been published.29,49,50 computed with TRAVIS.21–26,30,33,37,40–42,46,51–54 It is even supported
Recently, we implemented a function called geometric den- to combine three scalar quantities into a three-dimensional his-
sity analysis (GeoDens) to observe the temporal development of the togram in TRAVIS. The resulting histogram can be visualized by
average particle density within a geometric region along the trajec- means of iso-surfaces, as it has been shown in the literature a few
tory. In the specific application, we defined a cylinder between the times.24,33,55
two amino acid molecules involved in a salt bridge.31 The cylinder An even more direct way for understanding correlation
has a fixed radius, but the starting and ending points are tied to the between quantities is the so-called correlation plot in TRAVIS.
molecules so that the cylinder moves together with the molecules. By Consider Fig. 3 where a CDF with two distances is shown in the left
monitoring the particle density of anions within this cylinder relative panel. While the horizontal axis depicts the C–O bond length in an
to the uniform density, we were able to conclude that the anions are acetate anion within a bulk phase simulation,23,24 the vertical axis
repelled from the salt bridge (see Fig. S-18 of Ref. 31). shows the length of the other C–O bond in the same acetate anion. If
both quantities would be uncorrelated, the 2D histogram would sim-
ply be the Cartesian product of the two 1D bond length histograms.
B. Combined distribution functions To reveal the correlation, TRAVIS internally creates this Cartesian
If two observables are correlated, the individual histograms of product of the two histograms and subtracts it from the actual CDF.
the observables do not give a complete picture of the situation, and The resulting correlation plot is shown in the right-hand panel of
a multi-dimensional histogram is required to understand the corre- Fig. 3. Positive values indicate that the probability of finding this
lation effects (for example, see Fig. 2). Within a liquid phase tra- configuration is larger (i.e., positive correlation) than if the two quan-
jectory,23,24 the length of a certain hydrogen bond is depicted on tities would be uncorrelated, while negative values depict situations
the horizontal axis, while the hydrogen bond angle is shown on the less probable than in the uncorrelated case (negative correlation). We
vertical axis. The 2D contour plot possesses a very strong signal, find that configurations with one elongated and one shortened bond

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-4

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

FIG. 3. Two-dimensional histogram of

the two C–O bond lengths in an acetate
anion in [EMIm][OAc] (left panel);23,24
correlation plot for both quantities
reveals significant correlation (right
panel).

are more likely than expected from the individual bond length statis- a few of them, which are suitable to characterize lipid bi-layers that
tics, while cases where both bonds are either long or short are less we have been simulating.29 The first one is the tail vector tilting
likely—a clear correlation effect. analysis. If a lipid bi-layer is in the sol phase, the fatty acid tail
vectors are oriented in random directions and constantly fluctuat-
C. Structure factors ing. However, when the bi-layer undergoes a transition into the gel

01 November 2023 13:58:33

phase, the tail vectors are all tilted with the same angle into the
A very important experimental technique for condensed phase
same direction. By projecting all tail vectors into the X–Y plane
systems is measuring structure factors, also called S(q). The
(where the Z direction is the normal vector of the bi-layer), it can
structure factor is related to g(r) of the system by a Fourier
be easily seen if the bi-layer is in the sol or gel state (see Fig. 3 of
transform,
Ref. 29).
The second order parameter that is implemented in TRAVIS
S(q) = 1 + ρ ∫ e−iqr (g(r) − 1) dr. (1) is SCD , also called the deuterium order parameter, as it can be
directly measured by nuclear magnetic resonance (NMR) experi-
ments of deuterated samples.63,64 It is computed by the following
Despite this easy relation, it is slightly more involved to compute
equation:
S(q) from a simulation trajectory, as all pairs of atoms contribute
to it, and therefore, all possible pairwise RDFs of the system have to
be computed, correctly weighted by the scattering cross sections of 1 N T 3 2 1
SkCD = C–H
∑∑ cos ((rik (t), n)) − , (2)
the corresponding element types, and finally summed up. TRAVIS NT i=1 t=1 2 2
does this automatically and directly outputs the total S(q), together
with all partial structure factors (i.e., contribution to the total struc-
ture factor from a certain pair of elements).27 Several commonly where k depicts the position in the alkyl chain, N is the number
used normalization factors are available. Furthermore, it is possible of lipid molecules, T is the total number of time steps, rikC–H (t)
to compute individual intramolecular and intermolecular contribu- is one of the C–H (or C–D) bond vectors of the kth carbon
tions to the structure factor. Neutron scattering cross sections for atom of the chain in time step t, and n is the membrane normal
the most important isotopes56 as well as x-ray scattering cross sec- vector. SCD can take values between −0.5 and 1, where 0 corre-
tions (given by Cromer–Mann coefficients) are included in TRAVIS, sponds to the statistical distribution (i.e., no ordering) and −0.5
so both neutron and x-ray structure factors can be computed.27 indicates that the bond vectors are strictly within the membrane
Recently, the option to apply a Lorch-type window function57 to plane. An application of this analysis can be found in Fig. 4 of
the structure factor calculation was implemented. Structure factors Ref. 29.
computed by TRAVIS have been published in the literature sev- Finally, we implemented a function to compute the gauche ratio
eral times.27,52,55,58–62 TRAVIS can also compute dynamic structure along alkyl chains, i.e., the average fraction of C–C–C–C dihedral
factors S(q, ω), which is discussed in Sec. IV E. angles at a given chain position being close to 60○ , which is called the
gauche configuration. This configuration is energetically less favor-
able than the anti configuration with dihedral angles of around 180○
D. Order parameters
and therefore also statistically less populated. For an application of
There exists a large variety of so-called order parameters to this analysis and interpretation of the results, see Fig. 5 and Table 2
determine perturbed or ordered structures. TRAVIS only contains of Ref. 29.

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-5

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

E. Solvent orientation F. Hydrogen bond network topology

In condensed phase simulations of solutions, it is often very Many molecular liquids form hydrogen bonds, which often
interesting to understand how the solvent molecules are arranged have a significant influence on the properties (just consider water).
and oriented around the solute. Information on the solvent shells If there are several different hydrogen bond donors and acceptors in
can be gained from RDFs and SDFs, but no details on the orien- a complex system, a subtle balance of interactions is found, and even
tation of the individual solvent molecules are gained from these small changes in hydrogen bond donor or acceptor strengths can
analyses. CDFs involving distance and angle can give information lead to large effects. Therefore, it is very important to understand the
on the solvent orientation, but it is still not apparent what the pre- hydrogen bond topology and the competition for donor/acceptor
ferred solvent orientation at a specific position around the solute positions in such systems. A simple estimate for the strength of a
molecule is. To overcome this, we implemented the plane projec- hydrogen bond is the distance and height of the first maximum in the
tion analysis (PlProj). Consider Fig. 4 where an [EMIm]+ cation X⋯H RDF, together with the number integral up to the first mini-
from a simulation of one ion pair [EMIm][OAc] in water is fixed. mum, which gives the average number of hydrogen bonds of atom X.
The color scale of the plot depicts the average particle density of Let us consider a system where a cellulose oligomer is dissolved in a
water at the given position relative to uniform density (considering mixture of the ionic liquid [EMIm][OAc] and water. There are three
only a thin slice of the simulation cell defined by the [EMIm]+ ring hydrogen bond donors (cellulose alcohol H, water H, and [EMIm]+
plane). A value larger than 1 means that water is more frequently ring protons) and three acceptors (cellulose O, water O, and [OAc]−
found here than at some arbitrary position, while a value below 1 O). The hydrogen bond estimates from the pairwise RDFs could
indicates a depletion of water molecules at that point. It can be seen be written as a matrix of numbers, but this is not very illustrative.
that there exists a well-defined first solvation layer with particle den- Recently, we found a much more appealing way to represent such
sities above 1, which is separated from the water bulk phase by a data, which is based on Sankey diagrams.65 The concept is more than
clear region of depletion with particle densities below 1. The vec- 120 years old now, but it is rarely found in computational chem-
tors in the plot depict the average orientation of water dipole vectors istry and never has been used to visualize hydrogen bond topology
at each point. As expected, it can be seen that water forms hydro- to the best of our knowledge. Therefore, despite its early origin, we
gen bonds to the three protic ring hydrogen atoms, and therefore, consider this a “modern” visualization technique. Some common
the water dipole vector points away from these atoms on average. plotting libraries such as MatPlotLib66 are able to produce Sankey
Next to the non-polar regions of the molecule, there is no clear pref- diagrams, but not the circular variant described below. For a visu-

01 November 2023 13:58:33

erential orientation of the water molecules, resulting in very short alization of the hydrogen bond topology of the mixture described
average vectors. above, see Fig. 5. In the left-hand panel, a linear Sankey diagram is
shown, which contains the hydrogen bond donors on the left-hand
side and the acceptors on the right-hand side. The numbers corre-
spond to the average hydrogen bond count of the donor/acceptor.
The width of the connection bars is proportional to the hydrogen
bond number. Note that the change in bar width from left to right
arises from the different molecule counts. In the right-hand panel of
Fig. 5, the same data are visualized as a circular Sankey diagram.
Again, the numbers represent the total hydrogen bond count of the
donors/acceptors, and the connection width is proportional to the
hydrogen bond number for the two connected groups. We believe
that Sankey diagrams are a very illustrative and clear tool to rep-
resent such topologies, and we hope to see them in many future
computational chemistry papers from now on. As no common plot-
ting program is able to produce both linear and circular Sankey dia-
grams as shown, TRAVIS directly outputs the images as SVG vector
graphics.
For cases where an even more fine-grained analysis is required,
we developed another analysis, which reveals the hydrogen bond
pattern between several dozens of hydrogen bond donors and accep-
tors. For example, if it would be of interest which of the 16 dif-
ferent oxygen atoms in the cellulose pentamer forms hydrogen
bonds to which acceptors, this analysis will give a shortcut to the
desired results. Here, we consider a system with just cellulose and
[EMIm][OAc], but no water (see Fig. 6). We call this analysis the
connection matrix analysis (CMat). The rows of the matrix on the
left-hand side correspond to all different hydrogen bond acceptors
FIG. 4. Average particle density of water molecules around an [EMIm]+ cation in the system, while the columns represent the different hydrogen
(color scale) together with the average orientation of the water dipole vectors bond donors. For each of the matrix elements, TRAVIS internally
(arrows) in a simulation of one ion pair of [EMIm][OAc] in water.
computes an RDF and extracts the distance and height of the first

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-6

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

FIG. 5. Linear Sankey diagram representing the hydrogen bonding topology of cellulose dissolved in a mixture of [EMIm][OAc] and water (left). Circular Sankey diagram
depicting the same dataset (right), both were created by TRAVIS. Numbers correspond to the average hydrogen bond count per donor/acceptor.

01 November 2023 13:58:33

maximum. If the RDF does not indicate the presence of a hydro- distance, which is therefore not considered very strong. From look-
gen bond, the matrix element is simply filled with a black cross. If a ing at the matrix, it is immediately visible that the acetate anion
hydrogen bond exists, the square is filled with a color according to forms strong hydrogen bonds to many donors, while many cellu-
a two-dimensional color scale that is shown on the right-hand side lose oxygen atoms are not involved in strong hydrogen bonds. There
of Fig. 6. Both the distance and the height of the first RDF maxi- are two strong intramolecular hydrogen bonds clearly visible in cel-
mum are encoded in color. For example, red indicates a hydrogen lulose, namely, O3⋯H25 and O8⋯H31. Obtaining this information
bond with a very small distance and a very large maximum height just with standard RDF analyses would have required to compute
(i.e., a very strong hydrogen bond). Yellow, on the other hand, corre- and manually look at many hundred RDFs. We are convinced that
sponds to a hydrogen bond with a large maximum height but large the connection matrix analysis greatly simplifies that process.

FIG. 6. Connection matrix analysis of cellulose dissolved in [EMIm][OAc].32 Rows represent hydrogen bond acceptors, and columns stand for hydrogen bond donors. The
color in each square represents both the intensity and distance of the first maximum in the corresponding RDF (for color scale, see the right-hand side).

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-7

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

G. Voronoi analysis oxygen atoms and the protic ring hydrogen atoms can be directly
seen, while some interaction between the hydrogen atoms of the
Despite being more than 110 years old now, the Voronoi tessel- aliphatic side chains is also visible. Such a neighborhood matrix
lation67 is still a fascinating mathematical formalism. When applied captures the full neighborhood relationship of the system in one
to a simulation cell with the atoms as Voronoi sites, in simple words, image and can therefore be considered a fingerprint of the liquid
it partitions the box into Voronoi cells belonging to each atom in structure.
a way that all points in space that are closer to a particular atom
than to any other atom will be assigned to that atom’s Voronoi H. Domain analysis
cell. TRAVIS uses the free Voro++ library68,69 (integrated in the
To study microheterogeneity and microphase separation in
TRAVIS source code so that it does not have to be provided exter-
complex liquid systems, we developed the Voronoi-based domain
nally) to perform a periodic radical Voronoi tessellation.70,71 From
analysis (DomA) in 2015.28 First, the atoms of the system are cat-
such a tessellation, many interesting properties can be derived. On
egorized into different classes. This can be, e.g., polar, non-polar,
the one hand, Voronoi cells can be directly visualized by a 3D ren-
side chain, aromatic, fluorinated, etc. Then, a radical Voronoi tes-
dering via Povray.10 On the other hand, one can create statistics
sellation of the simulation cell is performed using all atoms of the
on the Voronoi cells’ volumes, surface areas, isoperimetric quo-
system as Voronoi sites. The Voronoi cells of atoms from the same
tients (measures for the “sphericity” of a cell), and so on. It is also
class, which are neighbors, i.e., share a common Voronoi face, are
possible to monitor Voronoi surface coverage, i.e., which percent-
merged. This leads to domains of different sizes and shapes. A
age of the faces of a molecule’s Voronoi cell is shared with other
two-dimensional illustration of this concept is presented in Fig. 8,
molecules of a certain type. Both the visualization and the statis-
where the different fill colors represent the different atom classes
tics can be performed with TRAVIS, as already demonstrated in the
and the bold black lines correspond to the domain boundaries.
literature.37,42,72,73
While in a perfectly homogeneous system, there will be many
Another interesting result from a Voronoi tessellation is a sim-
small domains, a certain degree of microheterogeneity will lead to
ple and unbiased neighborhood criterion: If the Voronoi cells of
a smaller number of larger domains because the atoms from the
two atoms share a common face, these two atoms are said to be
same class tend to form clusters. From performing this analysis in
neighbors. In contrast to other neighborhood criteria, e.g., based on
all trajectory frames and looking at the domain statistics (average
the first minimum of the RDF, this approach convinces through its
domain count, volume, surface area, isoperimetric quotient, etc.), one

01 November 2023 13:58:33

great simplicity, as it does not possess any empirical parameters. It
can gain much insight into the microphase behavior of the sys-
is therefore an unbiased criterion and transferable to all kinds of
tem. This has already been successfully applied many times in the
systems. In Fig. 7, the neighborhood probabilities for all pairs of
literature.27,60,73–78
atoms in a bulk phase simulation of [EMIm][OAc]23,24 are visual-
ized as a matrix. White corresponds to a neighborhood probability
I. Void analysis
of zero, i.e., these two atoms never share a common Voronoi face,
while purple represents a neighborhood probability of one, i.e., the In a solution, the exclusion volume of a solute perturbs the
one atom type always has at least one shared Voronoi face with an structure of the pure liquid. Obviously, the required amount of
atom of the other type. The strong hydrogen bonds between acetate energy for this perturbation is decreased if there already exist void

FIG. 7. Neighborhood probability matrix from the Voronoi analysis. The color in
each square depicts the probability that the two atom types from the row/column FIG. 8. Schematic representation of the Voronoi-based domain analysis.28 Colors
share a common Voronoi face at a time in a simulation of [EMIm][OAc]. represent different classes of atoms. Domains are highlighted by bold black lines.

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-8

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

cavities in the pure liquid, which are suited to accommodate the is slow to compute so many unnecessary distances. A more opti-
solute. Therefore, it can be expected that the size, shape, and mobil- mized approach would first define a maximum observation range.
ity of void cavities play a crucial role in many fundamental processes All distances beyond this range are neglected. Then, it would uti-
of solvation and mass transport in fluid systems. lize a spatial domain decomposition, as explained for the molecule
With the aim to quantify this effect, we have implemented a recognition in Sec. II A, using the maximum observation range as
Voronoi-based void analysis in TRAVIS, where a set of spheres is minimal cube diameter. Now, for a particular water molecule, only
grown into the voids between the molecules in a simulation cell the ions within the central cube and the 26 surrounding cubes need
(see Fig. 9).79 In contrast to other algorithms,80,81 we enhance this to be considered for forming three-tuples with that water molecule,
by the definition of a void domain: The spheres are handled as excluding all the ions that would not fulfill the conditions anyway.
pseudo-atoms and the corresponding Voronoi cells are combined This leads to an implementation that essentially scales linearly with
(cf. the domain analysis). With the aid of the isoperimetric quo- system size.
tient and autocorrelation functions, information on the shape and We are currently working on an unified formulation for all such
mobility of the void space can be extracted. analyses, which we call the generalized tuple analysis (GenTup).
Based on a central particle, as many particles as required can be
J. Generalized tuple analysis involved in the analysis at the same time, i.e., iterating over n-tuples
of particles for any number n, which all have to be located within a
Most of the particle density histograms discussed above only
maximum observation range around the central particle (the maxi-
involve two molecules at a time. Any RDF is based on the distance
mum observation range can be chosen by the user). It will be possible
between two atoms. In CDFs, there can be more atoms involved (e.g.,
to use many different condition relations between the particles (dis-
when observing a distance and an angle, as shown in Fig. 2), but there
tances, angles, dihedral angles, triangle surface area, tetrahedron vol-
are still only two molecules involved at a time. One can easily think
ume, etc.). To the best of our knowledge, such a generalized analysis
of analyses that involve more than two molecules.
has not yet been published in the literature. Our implementation in
For example, let us consider a simulation of a concentrated
TRAVIS will be available in some months.
aqueous sodium chloride solution. We want to investigate how the
O–H bond length histogram depends on the distance from the H
atom to Cl− ions, but only under the condition that the oxygen
IV. DYNAMICAL ANALYSES
atom is not further than 250 pm away from a Na+ ion. This anal- In this section, some dynamical analyses contained in TRAVIS

01 November 2023 13:58:33

ysis involves three particles at a time (one water molecule, one Na+ , will be shortly discussed. As already described above, dynamical
and one Cl− ). A straightforward approach would loop over all pos- analyses are defined as analyses that are not invariant under re-
sible 3-tuples of these particles, skipping all tuples that do not ful- shuffling of the trajectory frames, i.e., they depend on the ordering
fill the conditions and using the remaining tuples to populate the of the frames in some way. Most dynamical analyses should not be
multi-dimensional histogram. However, this is an inefficient algo- applied to Monte Carlo simulations, as there exists no physical time
rithm, as most of these tuples will not fulfill the conditions, and it in standard Monte Carlo approaches.

A. Diffusion coefficients
A common quantity when studying dynamics of condensed
phase systems are diffusion coefficients. They can be determined
both experimentally and via simulation, and are therefore a valu-
able property for validating simulations. The diffusion coefficient D
of some molecule can be computed from the mean square displace-
ment (MSD) via the Einstein relation83

1 ⟨∥p(t) − p(t + τ)∥2 ⟩t

D= lim , (3)
2d τ→∞ τ

where p(t) is the position vector of a particle at time t and d rep-

resents the dimensionality of the system (usually d = 3). TRAVIS
can compute MSDs and automatically performs a linear fit to the
second half of the obtained function to determine the diffusion coef-
ficient. Some years ago, a fast method for computing MSDs via
Fourier transform has been proposed in the literature;84 we are cur-
rently implementing this in TRAVIS, which will be available in some
months. MSDs and diffusion coefficients computed by TRAVIS have
appeared in the literature several times.22,29,33,40,49,58,77,85–87
FIG. 9. Example of the detected void spheres in a simulation box containing Another approach to diffusion coefficients of simulations is
N-methylimidazole and acetic acid (see Ref. 82 for details). Several single the Green–Kubo diffusion relation,88 which relates the integral of
spherical voids are combined into a much larger void space domain. the velocity autocorrelation function to the diffusion coefficient.

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-9

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

It reads criterion, it has shown to be effective to define a criterion for both the
∞ distance and angle (see the yellow rectangle in Fig. 2 for a reasonable
1
D= ⟨v(t) ⋅ v(t + τ)⟩t dτ, (4) hydrogen bond criterion). Then, an auxiliary function
d∫
0
1 if criteria fulfilled
where v(t) is the velocity vector of the considered particle at time βij (t) ∶= { (7)
0 otherwise
t and d again stands for the dimensionality of the system. In addi-
tion, the computation of velocity autocorrelation functions is imple- is defined, where βij (t) = 1 exactly when an aggregate between
mented in TRAVIS so that results from both approaches can be com- molecules i and j at time t is intact according to the chosen geometric
pared (they are typically not identical for finite-length trajectories, criterion. From this, the autocorrelation
and it is unclear which one is “better”).
∞
1 N N
B. Reorientation dynamics c(τ) = 2 ∑ ∑ ∫ βij (t) ⋅ βij (t + τ) dt (8)
N i=1 j=1
0
Another very interesting dynamical property is the reorien-
tation dynamics of certain vectors in a simulation, also known as is computed. The lifetime T of the aggregate can then be obtained as
rotational relaxation time. It can directly be determined by certain twice the total integral of this autocorrelation function
NMR and electron paramagnetic resonance (EPR) techniques and ∞
therefore offers a good opportunity for comparison between exper-
iment and simulation. Based on a trajectory, this quantity can easily T = 2 ⋅ ∫ c(τ) dτ. (9)
be expressed as an integral over the corresponding vector autocor- 0

relation function C(τ). For comparison with some experiments, it is The factor of 2 can be understood from the following simple con-
desirable to apply a Legendre polynomial Pn of order n to the vec- sideration: Imagine an aggregate that existed only once for a time
tor dot product in the correlation, then termed nth order correlation interval a. The corresponding correlation function c(τ) is then a
function Cn (τ) and nth order reorientation time T n . This leads to piecewise-linear function that runs through the two points c(0) = 1
the equations and c(a) = 0. The integral of this function will be a2 . A factor of 2 is
∞ required to retain the original lifetime of a.
u(t) ⋅ u(t + τ)

01 November 2023 13:58:33

Cn (τ) = ∫ Pn ( ) dt, (5) In practice, the correlation function c(τ) will often not decay
∥u(t)∥ ∥u(t + τ)∥ to zero within the simulation time. As explained in Sec. IV B,
0
∞
TRAVIS automatically fits a poly-exponential function to the cor-
relation function via a least-squares procedure (using a Levenberg–
Tn = ∫ Cn (τ) dτ, (6)
Marquardt minimizer),89 which is integrated analytically to obtain
0
the lifetime of the aggregate. Furthermore, in finite-size simulation
where u(t) is the chosen vector in the system at time t. As the cells, the autocorrelation function never decays to zero because any
first-order Legendre polynomial P1 is the identity, the reorientation two molecules that formed an aggregate will eventually meet again.
dynamics for n = 1 equals the case without the Legendre polynomial. To obtain correct lifetimes, it is important to subtract the ensemble
For typical simulation lengths, the vector autocorrelation function average value from the autocorrelation function before perform-
C(τ) does not decay to zero so that the integration to obtain T can- ing the exponential fit and the integration. TRAVIS automatically
not be directly performed. To overcome this, TRAVIS automatically performs this correction.
fits a poly-exponential function to C(τ) via a least-squares pro- The above-description corresponds to so-called intermittent
cedure (using a Levenberg–Marquardt minimizer),89 which is then autocorrelation functions. TRAVIS can also compute continuous
analytically integrated to determine T as a result. Vector reorienta- autocorrelation functions, which are not discussed here. Aggregate
tion times computed with TRAVIS have already been published in lifetimes and autocorrelation functions computed by TRAVIS have
the literature.20,31,58,90–92 been published in the literature several times.20,22,25,37,39,40,43,47,77,78

C. Aggregation dynamics 1. Efficient implementation

In computer simulation, it is often desirable to determine the Computing aggregate lifetimes as described above would be
average lifetimes of certain aggregates (e.g., ion pairs or hydrogen highly inefficient. Imagine a liquid water force field trajectory with
bond complexes). It is for sure not a valid approach to simply observe 1000 molecules and 106 trajectory steps. Storing the value of the
the lifetime of such a complex once in a trajectory—a molecular βij (t) function for all pairs and all time steps would require 1012
dynamics simulation is a chaotic system, and by chance, any arbi- memory cells, i.e., 1 TB of RAM if every function value is stored
trary lifetime can be obtained from such an “experiment” with some in 1 byte. A much better approach is to keep a list of aggregates.
non-zero probability. A valid approach is to run a simulation until For any aggregate between two molecules, it is only stored which
the investigated aggregate has formed and broken sufficiently many two molecules are involved in that aggregate and from which to
times and then use the autocorrelation formalism to determine its which time step the aggregate was intact. Such a list requires only
lifetime.93,94 To do so, first a geometric criterion for aggregation is a tiny fraction of the memory described above, and for every pair i
chosen. This can be a simple distance criterion or a more complex and j, the function βij (t) can be reconstructed from the list. By this
variant. For example, when defining hydrogen bonds via a geometric approach, only the data for one pair of i and j have to be stored in

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-10

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

memory at a time—all the resulting autocorrelation functions are

summed up anyway.
This is still not an optimal approach. It requires to compute
autocorrelations of i ⋅ j different time series. Autocorrelations can be
efficiently computed by using fast Fourier transform (FFT) due to
the Wiener–Khintchine theorem,95 but it still takes some time. An
even more efficient method can exploit the fact that βij (t) can only
take the values 0 and 1 and can therefore be written as the sum of
rectangular functions. Let
FIG. 10. Criteria for the reactive flux analysis. (Left) Hydrogen bond analysis:
1 a≤t<b aggregation if distance r AH and angle α are below a cutoff, vicinity defined via the
Rab (t) ∶= { (10) distance r AD . (Right) Ion pair analysis: aggregation defined via the next neighbor
0 otherwise
condition, vicinity defined via the distance r IC of the first solvation shell.
be a rectangular function. The cross-correlation of two such func-
tions Rab (t) and Rcd (t) has a simple piecewise-linear form (a trape- The inverse of the rate constant of decay kd is defined as the average
zoid) lifetime τlt = k−1
d of the interaction.
⎧ τ−E
⎪
⎪
⎪ F−E
, E≤τ<F The current implementation in TRAVIS30,99 enables the analy-
⎪
⎪
⎪1,
⎪ F≤τ<G sis for hydrogen bonding—based on a geometrical criterion as visu-
⟨Rab (t) ⋅ Rcd (t + τ)⟩t = ⎨ τ−H (11) alized in the left part of Fig. 10—as well as for ion pair dynamics
⎪
⎪
⎪ G−H
, G≤τ<H
⎪
⎪ based on a next neighbor criterion and a distance, as illustrated in
⎩0, otherwise
⎪
⎪
the right part of Fig. 10. The analysis of the interaction dynamics
with the substitutions allows a deeper insight into the molecular cosmos of chemical sys-
E ∶= c − b, tems,100 especially in organocatalysis,101 electrolytes,102 and metal
extractants.77
F ∶= min(c − a, d − b),
G ∶= max(c − a, d − b), E. Van Hove correlation functions

01 November 2023 13:58:33

H ∶= d − a. (12) As described above, an RDF is computed by taking the posi-
tions of two particles in the same trajectory frame, measuring the
As βij (t) can be written as a sum of rectangular functions (by know-
ing the aggregate lifetime intervals from the list), the resulting prod-
uct in the autocorrelation can be expanded (“multiplied out”) to a
sum of trapezoid functions as shown above. Then, the autocorre-
lation function for a pair of i and j can be directly constructed as
a sum of piecewise-linear functions without performing any actual
correlation at all. This highly efficient approach is implemented in
TRAVIS. To the best of our knowledge, it has not been described in
the literature before.

D. Reactive flux analysis

The lifetimes obtained from the approach described in
Sec. IV C are in some cases suffering from a diffusion bias. The
reactive flux approach of Luzar and Chandler96–98 eliminates this
perturbation with the aid of a second operator H(t). By defini-
tion, it becomes unity if a pair of interactants is close enough to
form the interaction in question and zero otherwise. This allows the
calculation of the function
τ
⟨H(t) [1 − h(t)] ḣ(0)⟩
n(τ) = ∫ (− )dt, (13)
⟨h⟩
0

which gives the probability of finding a pair that was initially inter-
acting at time t = 0 as not interacting anymore, but still in a distance
at which an interaction would be possible at time t = τ. Finally,
the two functions c(t) and n(t) can be fitted on the simple kinetic
equation FIG. 11. Van Hove correlation function103 (i.e., RDF with lag time) for the O⋯H
distance in liquid water, computed by TRAVIS. The 2D Fourier transform of this
− ċ(t) = kd c(t) − kf n(t). (14) function corresponds to the dynamic structure factor S(q, ω).

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-11

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

distance between the two positions, and adding the result to a his- transforming the time series, taking the absolute square, and trans-
togram. It is interesting to take the positions of the two particles forming back,
from different trajectory frames with a “lag time” of τ between them. ∞ 2
If this is performed for different lag times and represented as a C(τ) = ∫ f (t) ⋅ f (t + τ) dt = F −1 ∣F( f (τ))∣ .
⎛ ⎞
(15)
contour plot with the histogram value on the horizontal axis and ⎝ ⎠
−∞
the lag time on the vertical axis, the result is called a van Hove
correlation function (VHCF).103 In Fig. 11, a van Hove correla- As the Fourier transform of a discrete dataset of length n can be effi-
tion function for the O⋯H distance in a simulation of liquid ciently computed by the fast Fourier transform (FFT) method that
water is presented. It can be seen that the standard RDF is repro- only scales with O(n ⋅ log n), this saves a lot of computational time.
duced at τ = 0, but the features in the RDF quickly are blurred Equation (15) is only a special case of the cross-correlation theorem
out for an increase in lag time. Van Hove correlation functions so that a similar approach can also be used to speed up the cal-
are particularly interesting because their two-dimensional Fourier culation of cross-correlation functions via FFT. TRAVIS computes
transform is the dynamic structure factor S(q, ω), which can be all autocorrelation and cross-correlation functions based on these
measured by neutron scattering experiments. This offers another approaches using the KISS FFT library108 (integrated in the TRAVIS
opportunity to directly compare experimental and computational source code).
results. 2. Window function and zero padding
When considering Fourier transforms of data, it is very impor-
V. SPECTROSCOPIC ANALYSES tant to choose a suitable window function. There exists no Fourier
It is known since a long time that vibrational spectra can be transform without a window function—using no window function
computed from molecular dynamics simulations by time corre- means using a rectangular window, which is a very poor choice.
lation functions.104 This comes with several advantages over the Applying a window function simply means element-wise multipli-
widely used static–harmonic approach to spectra. Automatic con- cation of the discrete input data with the window function before
former sampling takes place during the MD simulation, the full Fourier transforming. By default, TRAVIS uses the Hann window
solvent effect can be included, realistic band shapes are directly function109 for computing spectra, which is given by
obtained, and even some anharmonic effects (such as temperature- πn
Wn = cos2 ( (16)

01 November 2023 13:58:33

dependent frequency shifts as well as approximate overtones and com- ),
2(N − 1)
bination bands) are covered.105 During the last few years, many
functions for predicting vibrational spectra have been implemented but other window functions (exponential, Gaussian) are also imple-
in TRAVIS. This section will give an overview on these meth- mented. The total integral of the Fourier transform of a discrete
ods and techniques. A detailed step-by-step tutorial on computing dataset ai depends only on the first value a0 . Therefore, any win-
vibrational spectra with CP2k106 and TRAVIS can be found on the dow function that fulfills the condition W 0 = 1 does not modify the
Internet.107 total integral at all. In other words, window functions alter the peak
shapes of the spectrum but maintain the integral intensities of the
A. General processing techniques bands.
This section describes some general techniques used in By adding zeroes to the tail of the correlation function before
TRAVIS, which are relevant for the calculation of all types of vibra- Fourier transforming, the resolution of the resulting spectrum can
tional spectra. While most of these techniques are basic knowledge be increased, which is called zero padding. No new spectral infor-
of signal processing, they are rarely discussed in articles on computer mation is gained by this technique, it is just an interpolation between
simulations. the data points of the original spectrum, i.e., a cosmetic measure. By
default, TRAVIS uses a zero padding factor of 4.
1. Fast correlation 3. Finite difference correction
All prediction methods for vibrational spectra from molec- For most spectra, the time derivative instead of the quan-
ular dynamics simulations are based on either autocorrelation or tity itself enters the correlation functions. For a cosine function
cross-correlation functions of some time series along the trajectory. cos(ωt + φ), it can be shown105 that the second-order central finite
Typically, these correlations are computed for molecular proper- difference derivative is given by
ties, not for the properties of the whole simulation box. The rea-
cos(ω(t+Δt)+φ) − cos(ω(t−Δt)+φ)
son is that correlating properties of the whole system is equiva- D cos(ωt+φ) =
lent to taking into account all cross-correlations between different 2Δt
molecules. However, the motion of very distant molecules is not sin(ωΔt)
= −ω sin(ωt + φ) , (17)
correlated at all so that these cross-correlations mostly introduce ωΔt
noise into the spectrum. Computing all the molecular correlation i.e., the exact value of the derivative is multiplied by a sinus cardinalis
functions requires much computational power because each such
(sinc) function sin(ωΔt)
ωΔt
. Since the product of two such time deriva-
calculation scales with O(m ⋅ n), where n is the length of the time
tives is formed in the correlation functions, the final spectra need to
series and m is the chosen correlation depth. Fortunately, there exists 2
the Wiener–Khintchine theorem,95 which states that computing the be divided by a factor of ( sin(ωΔt)
ωΔt
) to obtain the correct intensities.
autocorrelation C(τ) of a time series f (t) is equivalent to Fourier This correction is implemented in TRAVIS and active by default.

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-12

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

Due to the correction, high-quality spectra with accurate intensities time step Δt, it can be derived105 that the observed frequency ω̃ is
can be obtained even if the dipole moment only was computed every given by
4.0 fs.
1 1
4. Improved sampling via time reversibility ω̃ = arccos(1 − ω20 Δt 2 ). (21)
Δt 2
Similar to most other properties, vibrational spectra should be
Fortunately, the effect is rather small. When simulating an oscillator
invariant under reversal of the time direction: The trajectory with
with ω0 = 3000.0 cm−1 and Δt = 0.5 fs, the observed frequency is
reversed time samples the same thermodynamic ensemble as the for-
ω̃ = 2990.0 cm−1 , which is a deviation of around 0.3%. When com-
ward trajectory. As vibrational spectra are ensemble averages, they
puting spectra from a trajectory, the observed frequencies ω̃ are
cannot differ for these two cases. Consider the cross-correlations
given and the true frequencies ω0 are sought after. The inverse of
of two time series a(t) and b(t) along the trajectory. Let ã(t)
Eq. (21) can be easily written as
∶= a(T − t) and b̃(t) ∶= b(T − t) be the two time series of the tra-
√
jectory with reversed time, where T is the total trajectory length. We 2 − 2 cos(ω̃Δt)
find that ω0 = . (22)
Δt
⟨ã(t + τ) ⋅ b̃(t)⟩t = ∫ ã(t + τ) ⋅ b̃(t)dt By using Eq. (22), the frequency axis of the spectra can be corrected
for the frequency shift of the Verlet integrator. The correction is
= ∫ a(T − t − τ) ⋅ b(T − t)dt implemented in TRAVIS.

= ∫ a(t̃) ⋅ b(t̃ + τ)dt̃ B. Voronoi integration

= ⟨b(t̃ + τ) ⋅ a(t̃)⟩t̃ (18) When vibrational spectra shall be computed from molecu-
with the substitution t̃ ∶= T − t − τ. Thus, the reversal of the time lar dynamics simulations, molecular electric properties (e.g., dipole
direction is equivalent to swapping the first and second argument moments) are required. The standard approach to obtain those is
of the cross-correlations. If the spectra are invariant under time the so-called Wannier localization.111,112 A unitary transformation
reversal, the underlying cross-correlation functions should also be that minimizes the spatial spread is applied to the set of molecular
invariant. Therefore, swapping the two arguments will yield identi- orbitals, yielding localized Wannier orbitals. For each such orbital,

01 November 2023 13:58:33

cal correlation functions for an infinite trajectory length. However, the center of mass (called “Wannier center”) can be considered an
in practice, sampling is not perfect, and the two correlation functions integer point charge. By assigning the Wannier centers to individ-
differ due to statistical noise. By simply taking the arithmetic average ual molecules, the molecular dipole moment can be computed. In
of both correlation functions (with forward time and reversed time) the case of one isolated molecule, this is not even an approxima-
[see Eq. (19)], the sampling is effectively increased by a factor of 2. tion. However, this approach comes with several limitations. First, it
We call this approach the “commutator trick,” requires a lot of computational time (for systems in the order of 1000
atoms, easily more than the actual SCF cycle). Even worse is the fact
1 that the iterative approach to determine the Wannier centers some-
CCF(a, b) ∶= (⟨a(t + τ) ⋅ b(t)⟩t + ⟨b(t + τ) ⋅ a(t)⟩t ). (19)
2 times does not converge at all. Apart from that, there arise prob-
It only improves sampling if a(t) and b(t) are different time series, lems with aromatic systems. The Wannier localization enforces an
i.e., in the case of true cross-correlations. If a and b are identical (i.e., alternating single bond–double bond pattern in aromatic systems,
autocorrelation, as with correlating the dipole moment for infrared which can introduce artificial spectral bands, as we have recently
spectra), nothing is gained. reported.113 When computing molecular polarizabilities for Raman
Another subtlety arises if one of both properties a and b is the spectra of aromatic systems, the numerical noise introduced by the
magnetic dipole moment [which is required for vibrational circu- Wannier localization is several orders of magnitude larger than
lar dichroism (VCD) and Raman optical activity (ROA) spectra; see the spectrum itself. Finally, to the best of our knowledge, it is not
below]. This quantity depends linearly on the atomic velocities. A possible to derive molecular quadrupole moments from Wannier
time reversal inverts all atomic velocities and therefore swaps the centers; however, these quadrupole moments are required for the
sign of the magnetic dipole moment. In this case, taking the average computation of ROA spectra.
of both correlation functions would only yield noise, as both contri- To overcome these issues when computing molecular proper-
butions would cancel out. The sign of the correlation function with ties, we developed the Voronoi integration in 2015.113 After a peri-
reversed time has to be changed to obtain the desired effect, which is odic radical Voronoi tessellation67,70,71 of the system (with the atoms
given in the following equation: as Voronoi sites and atomic van der Waals radii114 as Voronoi radii)
has been computed by using the Voro++ library,68,69 the total elec-
CCF∗(a, b) ∶= (⟨a(t + τ) ⋅ b(t)⟩t − ⟨b(t + τ) ⋅ a(t)⟩t ).
1
(20) tron density ρ(r), which is supplied on a grid, is directly integrated
2 within the Voronoi cell of a molecule to yield the desired molecular
5. Correcting the Verlet frequency shift moments (e.g., electric dipole moment pMol and electric quadrupole
moment QMol ij ),
Describing the dynamics of a system with finite time steps
does unfortunately not leave the vibrational frequencies of har- NMol
monic oscillators invariant. If a harmonic oscillator with “true” fre- pMol = ∑ qn rn − ∫ ρ(r)r d3 r, (23)
quency ω0 is solved by a Verlet integration scheme110 with integrator n=1
Mol

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-13

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

NMol
different diatomic molecules with harmonic bonds (black curve) is
QMol
ij = ∑ qn (3rn,i rn,j − ∥rn ∥2 δij ) shown together with the integrals over the bands (red curve). As the
n=1
integrals are all of similar magnitude, the system seems to be well-
− ∫ ρ(r)(3ri rj − ∥r∥2 δij ) d3 r, (24) equilibrated. Power spectra computed by TRAVIS already appeared
Mol many times in the literature.22,25,33,37,41,72,115–117
where qn and rn are the core charge and position of the nth atom
in the molecule, respectively, and δij is the Kronecker delta. We D. Bulk phase normal modes
have shown that the spectra obtained via Voronoi integration are
Based on an approach presented by Mathias et al.,118,119 we
very similar to those based on Wannier localization, but the prob-
implemented a bulk phase normal mode analysis in TRAVIS.115
lems with aromatic systems are no longer present.113 Furthermore,
By iteratively diagonalizing the matrix of all atomic velocity cross-
our integration algorithm is very fast (∼1 s per full frame on one
correlation functions, the power spectrum of the system is decom-
CPU core) so that a lot of computational time is saved when com-
posed into normal modes that are maximally localized in the fre-
pared to Wannier localization. Finally, also higher molecular mul-
quency space. To do so, a set of reference structures is required,
tipole moments such as the electric quadrupole moment (which
onto which the molecules in the trajectory are projected. Due to
is required to compute ROA spectra) can easily be obtained. The
this requirement, the approach works best for rigid molecules and
required volumetric electron densities along the trajectory can be
has severe limitations for flexible side chains longer than ethyl
stored in our highly efficient lossless compression format BQB7 (see
groups. In Fig. 13, an exemplary application of the normal mode
above). Therefore, we propose the Voronoi integration technique as
analysis is shown, where the individual C–H stretching motions of
a new standard method for obtaining molecular properties in bulk
[EMIm]+ cations are distinguished in bulk phase simulations of pure
phase systems.
[EMIm][OAc] and an [EMIm][OAc]/water mixture.115 There also
exist other articles in the literature where this feature of TRAVIS has
C. Power spectra been applied.37,120–122
The Fourier transform of the atomic velocity autocorrelation
function yields the power spectrum of a system, also known as
E. Infrared and Raman spectra
the vibrational density of states (VDOS). In contrast to the exper-
imentally accessible vibrational spectra, the power spectrum is not The infrared spectrum of a system can be computed as the

01 November 2023 13:58:33

based on selection rules but contains all motions of a system. If the Fourier transform of the molecular dipole autocorrelation function
units are chosen correctly, and if atom mass weighting is applied, along a simulation trajectory. Soon after the first ab initio molec-
the integral over a peak in the power spectrum directly yields the ular dynamics (AIMD) simulations were performed, the first pre-
temperature of the corresponding degree of freedom in Kelvins. dicted infrared spectrum based on AIMD was published in 1997,123
This is an interesting feature to determine if a simulation is well and several applications and further developments followed.124–126
equilibrated: According to the equipartition theorem, all degrees For Raman spectra, the molecular polarizability needs to be cor-
of freedom should possess the same temperature. If the integrals related instead. The first predicted Raman spectrum from AIMD
in the power spectrum indicate that this might not be the case, appeared in 2002.127 It has been demonstrated many times in the
a longer equilibration can be recommended. One example is pre- literature that these spectra from AIMD simulations agree very well
sented in Fig. 12, where a power spectrum of a system containing 20 with the experiment. A module for the computation of such spectra
was added to TRAVIS in 2013.128 While it initially made use of Wan-
nier centers, it is now based on our Voronoi integration method (see
above).113 Infrared and Raman spectra computed by TRAVIS have
been published many times in the literature.87,122,129–134

F. VCD and ROA spectra

The vibrational circular dichroism (VCD) spectrum of a molec-
ular system can be computed via cross-correlating the molecular
electric dipole moment pMol (t) and the molecular magnetic dipole
moment mMol (t). However, the magnetic dipole moment cannot
be easily obtained from a Born–Oppenheimer molecular dynamics
(BOMD) simulation, as magnetic moments arise from electric cur-
rents, and the Born–Oppenheimer approximation does not account
for electric currents. One solution is to apply some kind of pertur-
bation theory (e.g., NVPT135,136 ), which, however, requires much
computational resources. Recently, we developed a purely classical
approach that reconstructs the electric currents from the difference
FIG. 12. Power spectrum of a system containing 20 different harmonic diatomic in total electron density in subsequent steps along the trajectory.137
molecules with different atom masses and force constants (black curve) together Once the electric current density j(r) has been obtained on a grid,
with the integral depicting the temperature of each normal mode (red curve). integration over the molecular Voronoi cells yields the molecular

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-14

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

FIG. 13. With the normal mode decom-

positions of the power spectra in sim-
ulations of pure [EMIm][OAc] and an
[EMIm][OAc]/water mixture, the C–H
stretching bands of the [EMIm]+ cation
can be distinguished.115

magnetic dipole moments mMol via the classical expression where vn is the velocity of the nth atom in the molecule. Based
on this approach, which is implemented in TRAVIS, we published
1 NMol 1 the first ab initio prediction of a bulk phase VCD spectrum in
mMol = 3
∑ qn (rn × vn ) − ∫ r × j(r) d r, (25) 2016.137
2 n=1 2

01 November 2023 13:58:33

Mol

FIG. 14. Schematic approach to compute the electromagnetic moments for ROA spectra from standard BOMD simulations.138 Rows represent successive BOMD time steps.
Red font depicts quantities under an external electric field.

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-15

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

Computing Raman optical activity (ROA) spectra from molec-

ular dynamics simulations is slightly more involved. The required
correlation functions are based on the electric dipole–electric dipole
polarizability αMol , the electric quadrupole–electric dipole polariz-
ability AMol , and the magnetic dipole–electric dipole polarizability
G′Mol . It can be shown139 that
d Mol
αMol = p , (26)
dE

d Mol
AMol = AMol = Q , (27)
dE

d Mol T
G′Mol = −(G′
Mol T
) = −( m ) (28)
dE

with the electric dipole–electric quadrupole polarizability AMol and

the electric dipole–magnetic dipole polarizability G′
Mol
so that all
three quantities can be obtained via finite differences of an applied FIG. 16. Predicted resonance Raman spectra of uracil in water for all possible laser
external electric field E. The schematic approach to compute the wavelengths (vertical axis).144 Rows are normalized to show the relative intensity
required polarizabilities for ROA spectra is presented in Fig. 14. We ratio. These were the first predicted bulk phase resonance Raman spectra in the
implemented this approach in TRAVIS and were able to present literature.
the first ab initio prediction of a bulk phase ROA spectrum in the
literature138 in 2017, which shows a very good agreement with the
experiment140 (see Fig. 15). phase systems and therefore is able to take into account full sol-
vent influence as well as some anharmonic effects. Apart from that,

01 November 2023 13:58:33

G. Resonance Raman spectra our approach is not limited to generalized gradient approximation
(GGA) DFT but works with all real time electron structure meth-
Resonance Raman spectra can be computed via temporal cor- ods, including, e.g., DFT with hybrid functionals. In our article, we
relation of the dynamic polarizability tensor along a trajectory. presented the first ab initio prediction of a bulk phase resonance
One approach to obtain this tensor is based on real-time time- Raman spectrum in the literature on the example of uracil in water.
dependent density functional theory (RT-TDDFT).141 It has been We have also shown that the solvent effect is very important in this
applied several times in the literature,142,143 but only to static cal- case, as only the bulk phase spectrum with full solvent influence
culations of isolated molecules in vacuum. In 2019, we imple- agrees well with the experiment.
mented this method in TRAVIS and applied it to BOMD trajecto- Our newly developed method for predicting resonance Raman
ries in order to obtain the resonance Raman spectrum.144 In con- spectra does not require a certain laser wavelength as input but
trast to the existing methods, our implementation works for bulk yields the full set of resonance Raman spectra for all possible
laser wavelengths in one pass. The spectra in dependence on the
laser wavelength can be visualized as a contour plot, as shown
in Fig. 16 on the example of uracil in water. This feature is not
only helpful to understand the coupling between vibrational modes
and electronic excitations in the system investigated but it is also
beneficial in the design of new interesting experiments by pick-
ing laser wavelengths at which interesting resonance effects can be
expected.

VI. CONCLUSION
In this article, we have presented some methods and analyses
that have recently been introduced to the TRAVIS program pack-
age while also giving a glimpse of some of the underlying algorithms
that are essential for efficiency and accuracy of these analyses. While
some functions in TRAVIS are known since decades, others are
very recent and have been implemented directly after having been
developed. For example, TRAVIS has been used to compute the
FIG. 15. Predicted Raman optical activity (ROA) spectrum138 of liquid first ab initio predictions in the literature of bulk vibrational cir-
(R)-propylene oxide (black) together with the experimental spectrum (red).140 This cular dichroism (VCD) spectra,137 bulk phase Raman optical activ-
was the first predicted bulk phase ROA spectrum in the literature. ity (ROA) spectra,138 and bulk phase resonance Raman spectra144

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-16

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

20
within the last years. Since the release in 2011, the original TRAVIS M. Kohagen, M. Brehm, J. Thar, W. Zhao, F. Müller-Plathe, and B. Kirchner,
article1 has already been cited more than 400 times now, and more “Performance of quantum chemically derived charges and persistence of ion
than 800 different co-authors contributed to these studies. As the cages in ionic liquids. A molecular dynamics simulations study of 1-n-butyl-3-
methylimidazolium bromide,” J. Phys. Chem. B 115, 693–702 (2011).
authors of TRAVIS, we promise that we will continue to include new 21
M. Brüssel, M. Brehm, T. Voigt, and B. Kirchner, “Ab initio molecular dynamics
methods and features into the code while at the same time main- simulations of a binary system of ionic liquids,” Phys. Chem. Chem. Phys. 13,
taining old and existing functions and improving the quality of the 13617–13620 (2011).
documentation. 22
A. S. Pensado, M. Brehm, J. Thar, A. P. Seitsonen, and B. Kirchner, “Effect
of dispersion on the structure and dynamics of the ionic liquid 1-ethyl-3-
methylimidazolium thiocyanate,” ChemPhysChem 13, 1845–1853 (2012).
ACKNOWLEDGMENTS 23
M. Brehm, H. Weber, A. S. Pensado, A. Stark, and B. Kirchner, “Proton
We would like to thank all the TRAVIS users around the transfer and polarity changes in ionic liquid-water mixtures: A perspective on
hydrogen bonds from ab initio molecular dynamics at the example of 1-ethyl-
world for their trust (and sometimes, unfortunately, also patience). 3-methylimidazolium acetate-water mixtures—Part 1,” Phys. Chem. Chem. Phys.
M.B. acknowledges financial support from the DFG through Project 14, 5030–5044 (2012).
No. Br 5494/1-1. 24
M. Brehm, H. Weber, A. S. Pensado, A. Stark, and B. Kirchner, “Liquid
The data that support the findings of this study are available structure and cluster formation in ionic liquid/water mixtures—An extensive
from the corresponding author upon reasonable request. ab initio molecular dynamics study on 1-ethyl-3-methylimidazolium acetate/water
mixtures—Part 2,” Z. Phys. Chem. 227, 177–204 (2013).
25
REFERENCES F. Malberg, M. Brehm, O. Hollóczki, A. S. Pensado, and B. Kirchner, “Under-
standing the evaporation of ionic liquids using the example of 1-ethyl-3-
1
M. Brehm and B. Kirchner, “TRAVIS—A free analyzer and visualizer for Monte methylimidazolium ethylsulfate,” Phys. Chem. Chem. Phys. 15, 18424–18436
Carlo and molecular dynamics trajectories,” J. Chem. Inf. Model. 51(8), 2007– (2013).
26
2023 (2011). M. Thomas, M. Brehm, O. Hollóczki, and B. Kirchner, “How can a carbene be
2
W. Humphrey, A. Dalke, K. Schulten et al., “VMD: Visual molecular dynamics,” active in an ionic liquid?,” Chem. - Eur. J. 20, 1622–1629 (2014).
27
J. Mol. Graphics 14, 33–38 (1996). O. Hollóczki, M. Macchiagodena, H. Weber, M. Thomas, M. Brehm, A. Stark,
3
N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, “MDAnal- O. Russina, A. Triolo, and B. Kirchner, “Triphilic ionic-liquid mixtures: Fluori-
ysis: A toolkit for the analysis of molecular dynamics simulations,” J. Comput. nated and non-fluorinated aprotic ionic-liquid mixtures,” ChemPhysChem 16,
Chem. 32, 2319–2327 (2011). 3325–3333 (2015).
4 28
I. T. Todorov, W. Smith, K. Trachenko, and M. T. Dove, “DL_POLY 3: New M. Brehm, H. Weber, M. Thomas, O. Hollóczki, and B. Kirchner, “Domain

01 November 2023 13:58:33

dimensions in molecular dynamics simulations via massive parallelism,” J. Mater. analysis in nanostructured liquids: A post-molecular dynamics study at the
Chem. 16, 1911–1918 (2006). example of ionic liquids,” ChemPhysChem 16, 3271–3277 (2015).
5 29
S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” M. Brehm, G. Saddiq, T. Watermann, and D. Sebastiani, “Influence of small
J. Comput. Phys. 117, 1–19 (1995). fluorophilic and lipophilic organic molecules on dipalmitoylphosphatidylcholine
6 bilayers,” J. Phys. Chem. B 121, 8311–8321 (2017).
D. A. Case, T. E. Cheatham III, T. Darden, H. Gohlke, R. Luo, K. M. Merz, Jr.,
30
A. Onufriev, C. Simmerling, B. Wang, and R. J. Woods, “The Amber biomolecular S. Gehrke, M. von Domaros, R. Clark, O. Hollóczki, M. Brehm, T. Welton,
simulation programs,” J. Comput. Chem. 26, 1668–1688 (2005). A. Luzar, and B. Kirchner, “Structure and lifetimes in ionic liquids and their
7 mixtures,” Faraday Discuss. 206, 219–245 (2018).
M. Brehm and M. Thomas, “An efficient lossless compression algorithm for
31
trajectories of atom positions and volumetric data,” J. Chem. Inf. Model. 58, S. Pylaeva, M. Brehm, and D. Sebastiani, “Salt bridge in aqueous solution: Strong
2092–2107 (2018). structural motifs but weak enthalpic effect,” Sci. Rep. 8, 13626 (2018).
8 32
Grace Development Team, Grace, (c), see http://plasma-gate.weizmann.ac.il/ M. Brehm, M. Pulst, J. Kressler, and D. Sebastiani, “Triazolium-based ionic
Grace, 1996–2008. liquids: A novel class of cellulose solvents,” J. Phys. Chem. B 123, 3994–4003
9 (2019).
T. Williams, C. Kelley et al., Gnuplot 4.6: An Interactive Plotting Program,
33
http://gnuplot.sourceforge.net/ (2013). M. Brüssel, M. Brehm, A. S. Pensado, F. Malberg, M. Ramzan, A. Stark, and
10
Persistence of Vision Pty. Ltd., Persistence of Vision Raytracer (Version 3.6), B. Kirchner, “On the ideality of binary mixtures of ionic liquids,” Phys. Chem.
retrieved from http://www.povray.org/download/ (2004). Chem. Phys. 14, 13204–13215 (2012).
11 34
W. R. Inc., Mathematica, Version 12.0, Champaign, IL, 2019. A. Korotkevich, D. S. Firaha, A. A. H. Padua, and B. Kirchner, “Ab initio molec-
12 ular dynamics simulations of SO2 solvation in choline chloride/glycerol deep
B. Cordero, V. Gómez, A. E. Platero-Prats, M. Revés, J. Echeverría, E. Cremades,
F. Barragán, and S. Alvarez, “Covalent radii revisited,” Dalton Trans. 2008, 2832– eutectic solvent,” Fluid Phase Equilib. 448, 59–68 (2017).
35
2838 (2008). S. Zahn, “Deep eutectic solvents: Similia similibus solvuntur?,” Phys. Chem.
13
R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (CRC Chem. Phys. 19 (5), 4041–4047 (2017).
Press, 1981). 36
L. Gontrani, F. Trequattrini, O. Palumbo, L. Bencivenni, and A. Paolone, “New
14
H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, “GROMACS: A experimental evidences regarding conformational equilibrium in ammonium−
message-passing parallel molecular dynamics implementation,” Comput. Phys. bis(trifluoromethanesulfonyl)imide ionic liquids,” ChemPhysChem 19(20), 2776–
Commun. 91, 43–56 (1995). 2781 (2018).
15 37
C. A. Shelley and M. E. Munk, “Computer perception of topological symmetry,” M. Brehm and D. Sebastiani, “Simulating structure and dynamics in small
J. Chem. Inf. Comput. Sci. 17, 110–113 (1977). droplets of 1-ethyl-3-methylimidazolium acetate,” J. Chem. Phys. 148, 193802
16
H. L. Morgan, “The generation of a unique machine description for chemical (2018).
structures: A technique developed at chemical abstracts service,” J. Chem. Doc. 5, 38
L. M. Lawson Daku, “Spin-state dependence of the structural and vibra-
107–113 (1965). tional properties of solvated iron(II) polypyridyl complexes from AIMD simu-
17
C. A. Shelley, “Heuristic approach for displaying chemical structures,” J. Chem. lations: II. Aqueous [Fe(TPy)2 ]Cl2 ,” Phys. Chem. Chem. Phys. 21(2), 650–661
Inf. Comput. Sci. 23, 61–65 (1983). (2019).
18 39
M. Brehm and M. Thomas, The BQB format, see https://brehm-research.de/bqb. J. Thar, M. Brehm, A. P. Seitsonen, and B. Kirchner, “Unexpected hydrogen
19
K. A. F. Röhrig and T. D. Kühne, “Optimal calculation of the pair correlation bond dynamics in imidazolium-based ionic liquids,” J. Phys. Chem. B 113, 15129–
function for an orthorhombic system,” Phys. Rev. E 87, 045301 (2013). 15132 (2009).

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-17

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

40 60
M. Kohagen, M. Brehm, Y. Lingscheid, R. Giernoth, J. Sangoro, F. Kremer, U. Kapoor and J. K. Shah, “Globular, sponge-like to layer-like morphological
S. Naumov, C. Iacob, J. Kärger, R. Valiullin, and B. Kirchner, “How hydrogen transition in 1-n-alkyl-3-methylimidazolium octylsulfate ionic liquid homologous
bonds influence the mobility of imidazolium-based ionic liquids. A combined series,” J. Phys. Chem. B 122, 213–228 (2018).
theoretical and experimental study of 1-N-butyl-3-methylimidazolium bromide,” 61
M. Campetella, A. Mariani, C. Sadun, B. Wu, E. W. Castner, and L. Gontrani,
J. Phys. Chem. B 115, 15280–15288 (2011). “Structure and dynamics of propylammonium nitrate-acetonitrile mixtures: An
41
K. Wendler, M. Brehm, F. Malberg, B. Kirchner, and L. Delle Site, “Short intricate multi-scale system probed with experimental and theoretical techniques,”
time dynamics of ionic liquids in AIMD-based power spectra,” J. Chem. Theory J. Chem. Phys. 148(13), 134507 (2018).
Comput. 8, 1570–1579 (2012). 62
F. L. Celso, Y. Yoshida, R. Lombardo, C. J. Jafta, L. Gontrani, A. Triolo, and
42
O. Hollóczki, D. S. Firaha, J. Friedrich, M. Brehm, R. Cybik, M. Wild, A. Stark, O. Russina, “Mesoscopic structural organization in fluorinated room temperature
and B. Kirchner, “Carbene formation in ionic liquids: Spontaneous, induced, or ionic liquids,” C. R. Chim. 21(8), 757–770 (2018).
prohibited?,” J. Phys. Chem. B 117, 5898–5907 (2013). 63
L. S. Vermeer, B. L. de Groot, V. Réat, A. Milon, and J. Czaplicki, “Acyl chain
43
G. Bekçioǧlu, C. Allolio, and D. Sebastiani, “Water wires in aqueous solutions order parameter profiles in phospholipid bilayers: Computation from molec-
from first-principles calculations,” J. Phys. Chem. B 119, 4053–4060 (2015). ular dynamics simulations and comparison with 2H NMR experiments,” Eur.
44
M. L. S. Batista, G. Pérez-Sánchez, J. R. B. Gomes, J. A. P. Coutinho, and E. Biophys. J. 36, 919–931 (2007).
64
J. Maginn, “Evaluation of the GROMOS 56ACARBO force field for the calculation R. Giernoth, A. Bröhl, M. Brehm, and Y. Lingscheid, “Interactions in
of structural, volumetric, and dynamic properties of aqueous glucose systems,” ionic liquids probed by in situ NMR spectroscopy,” J. Mol. Liq. 192, 55–58
J. Phys. Chem. B 119, 15310–15319 (2015). (2014).
45 65
Q. R. Sheridan, S. Oh, O. Morales-Collazo, E. W. Castner, J. F. Brennecke, and A. B. W. Kennedy and H. R. Sankey, “The thermal efficiency of steam engines.
E. J. Maginn, “Liquid structure of CO2 -reactive aprotic heterocyclic anion ionic Report of the committee appointed to the council upon the subject of the defini-
liquids from X-ray scattering and molecular dynamics,” J. Phys. Chem. B 120, tion of a standard or standards of thermal efficiency for steam engines: With an
11951–11960 (2016). introductory note. (Including appendixes and plate at back of volume),” Minutes
46
F. Lo Celso, G. B. Appetecchi, C. J. Jafta, L. Gontrani, J. N. Canongia Proc. Inst. Civil Eng. 134, 278–312 (1898).
66
Lopes, A. Triolo, and O. Russina, “Nanoscale organization in the flu- J. D. Hunter, “Matplotlib: A 2D graphics environment,” Comput. Sci. Eng. 9,
orinated room temperature ionic liquid: Tetraethyl ammonium (trifluo- 90–95 (2007).
romethanesulfonyl)(nonafluorobutylsulfonyl)imide,” J. Chem. Phys. 148, 193816 67
G. Voronoi, “Nouvelles applications des paramètres continus à la théorie
(2018). des formes quadratiques. Deuxième mémoire. Recherches sur les parallélloèdres
47
K. Bernardino, T. A. Lima, and M. C. C. Ribeiro, “Low-temperature phase tran- primitifs,” J. Reine Angew. Math. 1908(134), 198 (1908).
sitions of the ionic liquid 1-ethyl-3-methylimidazolium dicyanamide,” J. Phys. 68
C. H. Rycroft, “Voro++: A three-dimensional Voronoi cell library in C++,”
Chem. B 123(44), 9418–9427 (2019). Chaos 19, 041111 (2009).
48 69
M. Zhao, B. Wu, S. I. Lall-Ramnarine, J. D. Ramdihal, K. A. Papacostas, E. C. H. Rycroft, G. S. Grest, J. W. Landry, and M. Z. Bazant, “Analysis of granular

01 November 2023 13:58:33

D. Fernandez, R. A. Sumner, C. J. Margulis, J. F. Wishart, and E. W. Castner, flow in a pebble-bed nuclear reactor,” Phys. Rev. E 74, 021306 (2006).
“Structural analysis of ionic liquids with symmetric and asymmetric fluorinated 70
B. J. Gellatly and J. L. Finney, “Calculation of protein volumes: An alternative to
anions,” J. Chem. Phys. 151, 074504 (2019). the Voronoi procedure,” J. Mol. Biol. 161, 305–322 (1982).
49
X.-Y. Guo, C. Peschel, T. Watermann, G. F. v. Rudorff, and D. Sebastiani, “Clus- 71
F. M. Richards, “The interpretation of protein structures: Total volume, group
ter formation of polyphilic molecules solvated in a DPPC bilayer,” Polymers 9, 488 volume distributions and packing density,” J. Mol. Biol. 82, 1–14 (1974).
(2017). 72
D. S. Firaha, M. Kavalchuk, and B. Kirchner, “SO2 solvation in the 1-ethyl-3-
50
M. H. Dokoohaki, A. R. Zolghadr, and A. Klein, “Impact of the chemi- methylimidazolium thiocyanate ionic liquid by incorporation into the extended
cal structure on the distribution of neuroprotective N-alkyl-9H-carbazoles at cation–anion network,” J. Solution Chem. 44, 838–849 (2015).
octanol/water interfaces,” New J. Chem. 44(4), 1211–1220 (2020). 73
T. D. N. Reddy and B. S. Mallik, “Heterogeneity in the microstructure and
51
L. K. Scarbath-Evers, R. Hammer, D. Golze, M. Brehm, D. Sebastiani, and dynamics of tetraalkylammonium hydroxide ionic liquids: Insight from classical
W. Widdra, “From flat to tilted: Gradual interfaces in organic thin film growth,” molecular dynamics simulations and Voronoi tessellation analysis,” Phys. Chem.
Nanoscale 12, 3834–3845 (2020). Chem. Phys. 22(6), 3466–3480 (2020).
52
O. Russina, F. L. Celso, and A. Triolo, “Pressure-responsive mesoscopic struc- 74
A. Mariani, R. Caminiti, M. Campetella, and L. Gontrani, “Pressure-induced
tures in room temperature ionic liquids,” Phys. Chem. Chem. Phys. 17, 29496– mesoscopic disorder in protic ionic liquids: First computational study,” Phys.
29500 (2015). Chem. Chem. Phys. 18 (4), 2297–2302 (2016).
53
T. C. Lourenço, Y. Zhang, L. T. Costa, and E. J. Maginn, “A molecular dynam- 75
D. O. Abranches, N. Schaeffer, L. P. Silva, M. A. R. Martins, S. P. Pinho, and J. A.
ics study of lithium-containing aprotic heterocyclic ionic liquid electrolytes,” P. Coutinho, “The role of charge transfer in the formation of type I deep eutectic
J. Chem. Phys. 148, 193834 (2018). solvent-analogous ionic liquid mixtures,” Molecules 24(20), 3687 (2019).
54 76
F. Lo Celso, G. B. Appetecchi, E. Simonetti, M. Zhao, E. W. Castner, U. T. Cosby, U. Kapoor, J. K. Shah, and J. Sangoro, “Mesoscale organization and
Keiderling, L. Gontrani, A. Triolo, and O. Russina, “Microscopic structural and dynamics in binary ionic liquid mixtures,” J. Phys. Chem. Lett. 10(20), 6274–6280
dynamic features in triphilic room temperature ionic liquids,” Front. Chem. 7, 285 (2019).
(2019). 77
R. Macchieraldo, S. Gehrke, N. K. Batchu, B. Kirchner, and K. Binnemans,
55
M. Macchiagodena, G. Mancini, M. Pagliai, G. Cardini, and V. Barone, “New “Tuning solvent miscibility: A fundamental assessment at the example of induced
atomistic model of pyrrole with improved liquid state properties and structure,” methanol/n-dodecane phase separation,” J. Phys. Chem. B 123, 4400–4407
Int. J. Quantum Chem. 118(9), e25554 (2018). (2019).
56 78
V. F. Sears, “Neutron scattering lengths and cross sections,” Neutron News 3, T. D. N. Reddy and B. S. Mallik, “Nanostructure domains, voids, and low-
26–37 (1992). frequency spectra in binary mixtures of N,N-dimethylacetamide and ionic liquids
57
E. Lorch, “Neutron diffraction by germania, silica and radiation-damaged silica with varying cationic size,” RSC Adv. 10(3), 1811–1827 (2020).
79
glasses,” J. Phys. C: Solid State Phys. 2, 229–237 (1969). S. Gehrke, R. Macchieraldo, and B. Kirchner, “Understanding the fluidity of
58
F. L. Celso, B. Aoun, A. Triolo, and O. Russina, “Liquid structure of dibutyl condensed phase systems in terms of voids—Novel algorithm, implementation
sulfoxide,” Phys. Chem. Chem. Phys. 18(23), 15980–15987 (2016). and application,” Phys. Chem. Chem. Phys. 21, 4988–4997 (2019).
59 80
M. Campetella, M. Montagna, L. Gontrani, E. Scarpellini, and E. Bodo, “Unex- X. Huang, C. J. Margulis, Y. Li, and B. J. Berne, “Why is the partial molar volume
pected proton mobility in the bulk phase of cholinium-based ionic liquids: New of CO2 so small when dissolved in a room temperature ionic liquid? Structure and
insights from theoretical calculations,” Phys. Chem. Chem. Phys. 19, 11869–11880 dynamics of CO2 dissolved in [BMIm+ ][PF6 − ],” J. Am. Chem. Soc. 127, 17842–
(2017). 17851 (2005).

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-18

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

81 103
N. N. Medvedev, V. P. Voloshin, V. A. Luchnikov, and M. L. Gavrilova, “An L. Van Hove, “Correlations in space and time and Born approxima-
algorithm for three-dimensional Voronoi S-network,” J. Comput. Chem. 27, tion scattering in systems of interacting particles,” Phys. Rev. 95, 249–262
1676–1692 (2006). (1954).
82 104
J. Ingenmey, S. Gehrke, and B. Kirchner, “How to harvest Grotthuss diffusion R. G. Gordon, “Molecular motion in infrared and Raman spectra,” J. Chem.
in protic ionic liquid electrolyte systems,” ChemSusChem 11, 1900–1910 (2018). Phys. 43, 1307–1312 (1965).
83 105
A. Einstein, “Über die von der molekularkinetischen theorie der wärme M. Thomas, “Theoretical modeling of vibrational spectra in the liquid phase,”
geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen,” Ann. Springer theses, Springer International Publishing, 2016.
Phys. 322, 549–560 (1905). 106
J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele, “CP2K: Atomistic
84
V. Calandrini, E. Pellegrini, P. Calligari, K. Hinsen, and G. R. Kneller, simulations of condensed matter systems,” Wiley Interdiscip. Rev.: Comput. Mol.
“nMoldyn—Interfacing spectroscopic experiments, molecular dynamics simula- Sci. 4, 15–25 (2014).
tions and models for time correlation functions,” École Thématique Soc. Française 107
M. Brehm and M. Thomas, Tutorial on Computing Vibrational Spectra, see
Neutronique 12, 201–232 (2011). https://brehm-research.de/spectroscopy.
85
R. Stefanovic, G. B. Webber, and A. J. Page, “Nanostructure of propylammo- 108
M. Borgerding, KISS FFT Library, see http://sourceforge.net/projects/kissfft/.
nium nitrate in the presence of poly(ethylene oxide) and halide salts,” J. Chem. 109
P. Kahlig, “Some aspects of Julius Von Hann’s contribution to modern clima-
Phys. 148(19), 193826 (2018). tology,” in Interactions Between Global Climate Subsystems (American Geophysi-
86
U. Cerajewski, J. Träger, S. Henkel, A. H. Roos, M. Brehm, and D. cal Union, 2013), pp. 1–7.
Hinderberger, “Nanoscopic structures and molecular interactions leading to a 110
L. Verlet, “Computer “experiments” on classical fluids. I. Thermodynamical
dystectic and two eutectic points in [EMIm][Cl]/urea mixtures,” Phys. Chem. properties of Lennard–Jones molecules,” Phys. Rev. 159, 98–103 (1967).
Chem. Phys. 20, 29591–29600 (2018). 111
87
G. H. Wannier, “The structure of electronic excitation levels in insulating
H. Scheiber, Y. Shi, and R. Z. Khaliullin, “Communication: Compact orbitals crystals,” Phys. Rev. 52, 191–197 (1937).
enable low-cost linear-scaling ab initio molecular dynamics for weakly-interacting 112
N. Marzari and D. Vanderbilt, “Maximally localized generalized Wannier
systems,” J. Chem. Phys. 148(23), 231103 (2018).
88
functions for composite energy bands,” Phys. Rev. B 56, 12847–12865 (1997).
M. S. Green, “Markoff random processes and the statistical mechanics of time- 113
M. Thomas, M. Brehm, and B. Kirchner, “Voronoi dipole moments for the
dependent phenomena. II. Irreversible processes in fluids,” J. Chem. Phys. 22,
simulation of bulk phase vibrational spectra,” Phys. Chem. Chem. Phys. 17, 3207–
398–413 (1954).
89
3213 (2015).
J. Wuttke, LMFIT—A C Library for Levenberg-Marquardt Least-Squares Min- 114
M. Mantina, A. C. Chamberlin, R. Valero, C. J. Cramer, and D. G. Truhlar,
imization and Curve Fitting, based on work by B. S. Garbow, K. E. Hillstrom,
“Consistent van der Waals radii for the whole main group,” J. Phys. Chem. A 113,
J. J. Moré, and S. Moshier, retrieved in 2010 from https://jugit.fz-juelich.de/mlz/
5806–5812 (2009).
lmfit. 115
90 M. Thomas, M. Brehm, O. Hollóczki, Z. Kelemen, L. Nyulászi, T. Pasinszki,
M. Macchiagodena, G. D. Frate, G. Brancato, B. Chandramouli, G. Mancini,

01 November 2023 13:58:33

and B. Kirchner, “Simulating the vibrational spectra of ionic liquid systems:
and V. Barone, “Computational study of the DPAP molecular rotor in various
1-ethyl-3-methylimidazolium acetate and its mixtures,” J. Chem. Phys. 141,
environments: From force field development to molecular dynamics simulations
024510 (2014).
and spectroscopic calculations,” Phys. Chem. Chem. Phys. 19(45), 30590–30602 116
M. C. Kirkegaard, J. L. Niedziela, A. Miskowiec, A. E. Shields, and B. B.
(2017).
91 Anderson, “Elucidation of the structure and vibrational spectroscopy of synthetic
M. H. Kowsari and S. Ebrahimi, “Capturing the effect of [PF3 (C2 F5 )3 ]− vs.
metaschoepite and its dehydration product,” Inorg. Chem. 58(11), 7310–7323
[PF6 ]− , flexible anion vs. rigid, and scaled charge vs. unit on the transport prop-
(2019).
erties of [BMIm]+ -based ionic liquids: A comparative MD study,” Phys. Chem. 117
T. Bauer, S. Maisel, D. Blaumeiser, J. Vecchietti, N. Taccardi, P. Wasser-
Chem. Phys. 20 (19), 13379–13393 (2018).
92 scheid, A. Bonivardi, A. Görling, and J. Libuda, “Operando drifts and DFT study
C. Dreßler and D. Sebastiani, “Effect of anion reorientation on proton mobility
of propane dehydrogenation over solid- and liquid-supported GaxPty catalysts,”
in the solid acids family CsHy XO4 (X = S, P, Se, y = 1, 2) from ab initio molecular
ACS Catal. 9(4), 2842–2853 (2019).
dynamics simulations,” Phys. Chem. Chem. Phys. (published online 2019). 118
93 G. Mathias and M. D. Baer, “Generalized normal coordinates for the vibra-
F. Stillinger, “Theoretical approaches to the intermolecular nature of water,”
tional analysis of molecular dynamics simulations,” J. Chem. Theory Comput. 7,
Philos. Trans. R. Soc., B 278, 97–112 (1977).
94
2028–2039 (2011).
D. C. Rapaport, “Hydrogen bonds in water: Network organization and life- 119
G. Mathias, S. D. Ivanov, A. Witt, M. D. Baer, and D. Marx, “Infrared spec-
times,” Mol. Phys. 50, 1151–1162 (1983). troscopy of fluxional molecules from (ab initio) molecular dynamics: Resolv-
95
A. Khintchine, “Korrelationstheorie der stationären stochastischen prozesse,” ing large-amplitude motion, multiple conformations, and permutational symme-
Math. Ann. 109(1), 604–615 (1934). tries,” J. Chem. Theory Comput. 8, 224–234 (2012).
96
A. Luzar and D. Chandler, “Effect of environment on hydrogen bond dynamics 120
B. Koeppe, S. A. Pylaeva, C. Allolio, D. Sebastiani, E. T. J. Nibbering, G.
in liquid water,” Phys. Rev. Lett. 76, 928 (1996). S. Denisov, H.-H. Limbach, and P. M. Tolstoy, “Polar solvent fluctuations drive
97
A. Luzar and D. Chandler, “Hydrogen-bond kinetics in liquid water,” Nature proton transfer in hydrogen bonded complexes of carboxylic acid with pyridines:
379, 55 (1996). NMR, IR and ab initio MD study,” Phys. Chem. Chem. Phys. 19(2), 1010–1028
98
A. Luzar, “Resolving the hydrogen bond dynamics conundrum,” J. Chem. Phys. (2017).
113, 10663–10675 (2000). 121
M. Gług, M. Z. Brela, M. Boczar, A. M. Turek, Ł. Boda, M. J. Wójcik, T.
99
S. Gehrke and B. Kirchner, “Robustness of the hydrogen bond and ion pair Nakajima, and Y. Ozaki, “Infrared spectroscopy and Born–Oppenheimer molecu-
dynamics in ionic liquids to different parameters from the reactive flux method,” lar dynamics simulation study on deuterium substitution in the crystalline benzoic
J. Chem. Eng. Data 65, 1146–1158 (2019). acid,” J. Phys. Chem. B 121(3), 479–489 (2017).
100 122
S. Ghahramani, F. Yousefi, S. M. Hosseini, and S. Aparicio, “High-pressure L. M. Lawson Daku, “Spin-state dependence of the structural and vibrational
behavior of 2-hydroxyethylammonium acetate ionic liquid: Experiment and properties of solvated iron(II) polypyridyl complexes from AIMD simulations:
molecular dynamics,” J. Supercrit. Fluids 155, 104664 (2020). Aqueous [Fe(bpy)3 ]Cl2 , a case study,” Phys. Chem. Chem. Phys. 20, 6236–6253
101
S. Gehrke and O. Hollóczki, “Hydrogen bonding of N-heterocyclic carbenes in (2018).
123
solution: Mechanisms of solvent reorganization,” Chem. - Eur. J. 24, 11594–11604 P. L. Silvestrelli, M. Bernasconi, and M. Parrinello, “Ab initio infrared spectrum
(2018). of liquid water,” Chem. Phys. Lett. 277, 478–482 (1997).
102 124
T. Stettner, S. Gehrke, P. Ray, B. Kirchner, and A. Balducci, “Water in protic M. P. Gaigeot, R. Vuilleumier, M. Sprik, and D. Borgis, “Infrared spectroscopy
ionic liquids: Properties and use of a novel class of electrolytes for energy storage of N-methylacetamide revisited by ab initio molecular dynamics simulations,”
devices,” ChemSusChem 12, 3827–3836 (2019). J. Chem. Theory Comput. 1, 772–789 (2005).

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-19

© Author(s) 2020
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

125 135
M.-P. Gaigeot, M. Martinez, and R. Vuilleumier, “Infrared spectroscopy in A. Scherrer, R. Vuilleumier, and D. Sebastiani, “Nuclear velocity perturbation
the gas and liquid phase from first principle molecular dynamics simulations: theory of vibrational circular dichroism,” J. Chem. Theory Comput. 9, 5305–5312
Application to small peptides,” Mol. Phys. 105, 2857–2878 (2007). (2013).
126 136
M.-P. Gaigeot, N. A. Besley, and J. D. Hirst, “Modeling the infrared and circular A. Scherrer, R. Vuilleumier, and D. Sebastiani, “Vibrational circular dichroism
dichroism spectroscopy of a bridged cyclic diamide,” J. Phys. Chem. B 115, 5526– from ab initio molecular dynamics and nuclear velocity perturbation theory in the
5535 (2011). liquid phase,” J. Chem. Phys. 145, 084101 (2016).
127 137
A. Putrino and M. Parrinello, “Anharmonic Raman spectra in high-pressure M. Thomas and B. Kirchner, “Classical magnetic dipole moments for the sim-
ice from ab initio simulations,” Phys. Rev. Lett. 88, 176401 (2002). ulation of vibrational circular dichroism by ab initio molecular dynamics,” J. Phys.
128
M. Thomas, M. Brehm, R. Fligg, P. Vöhringer, and B. Kirchner, “Computing Chem. Lett. 7, 509–513 (2016).
138
vibrational spectra from ab initio molecular dynamics,” Phys. Chem. Chem. Phys. M. Brehm and M. Thomas, “Computing bulk phase Raman optical activity
15, 6608–6622 (2013). spectra from ab initio molecular dynamics simulations,” J. Phys. Chem. Lett. 8,
129
M. Etinski and B. Ensing, “Puzzle of the intramolecular hydrogen bond of 3409–3414 (2017).
139
dibenzoylmethane resolved by molecular dynamics simulations,” J. Phys. Chem. A D. A. Long, The Raman Effect: A Unified Treatment of the Theory of Raman
122(28), 5945–5954 (2018). Scattering by Molecules (John Wiley & Sons Ltd., 2002).
130 140
E. O. Fetisov, D. B. Harwood, I.-F. W. Kuo, S. E. E. Warrag, M. C. Kroon, J. Šebestík and P. Bouř, “Raman optical activity of methyloxirane gas and
C. J. Peters, and J. I. Siepmann, “First-principles molecular dynamics study of a liquid,” J. Phys. Chem. Lett. 2, 498–502 (2011).
141
deep eutectic solvent: Choline chloride/urea and its mixture with water,” J. Phys. H. Chen, J. M. Mcmahon, M. A. Ratner, and G. C. Schatz, “Classical elec-
Chem. B 122(3), 1245–1254 (2018). trodynamics coupled to quantum mechanics for calculation of molecular optical
131 properties: A RT-TDDFT/FDTD approach,” J. Phys. Chem. C 114, 14384–14392
M. T. Ruggiero, J. Kölbel, Q. Li, and J. A. Zeitler, “Predicting the structures and
associated phase transition mechanisms in disordered crystals via a combination (2010).
142
of experimental and theoretical methods,” Faraday Discuss. 211, 425–439 (2018). M. Thomas, F. Latorre, and P. Marquetand, “Resonance Raman spectra of
132
G. Cassone, J. Sponer, S. Trusso, and F. Saija, “Ab initio spectroscopy of water ortho-nitrophenol calculated by real-time time-dependent density functional the-
under electric fields,” Phys. Chem. Chem. Phys. 21(38), 21205–21212 (2019). ory,” J. Chem. Phys. 138, 044101 (2013).
133 143
B. Dutta and J. Chowdhury, “Existence of dimeric hydroxylamine-O-sulfonic J. Mattiat and S. Luber, “Efficient calculation of (resonance) Raman spectra
acid: Experimental observations aided by ab initio, DFT, Car–Parrinello and and excitation profiles with real-time propagation,” J. Chem. Phys. 149, 174108
Born–Oppenheimer on the fly dynamics,” Chem. Phys. Lett. 732, 136645 (2019). (2018).
134 144
H. Hoshina, T. Kanemura, and M. T. Ruggiero, “Exploring the dynamics of M. Brehm and M. Thomas, “Computing bulk phase resonance Raman spec-
bound water in nylon polymers with terahertz spectroscopy,” J. Phys. Chem. B tra from ab initio molecular dynamics and real-time TDDFT,” J. Chem. Theory
124(2), 422–429 (2020). Comput. 15, 3901–3905 (2019).

01 November 2023 13:58:33

J. Chem. Phys. 152, 164105 (2020); doi: 10.1063/5.0005078 152, 164105-20

© Author(s) 2020