TRAVIS (“Trajectory Analyzer and Visualizer”) is a program package for post-processing and analyzing trajectories from molecular dynamics 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.

Within the last few decades, the methods of molecular dynamics (MD) and Monte Carlo simulations have become a central part of computational chemistry and physics. In contrast to static calculations, they implicitly include the effect of entropy at finite temperatures. Most experimentally accessible observables are only defined as ensemble averages, rendering it indispensable to rely on ensemble averaging also in computations if experimental results shall be reproduced or predicted. The direct result of a simulation is a trajectory that contains the positions and velocities of the particles in equidistant discrete time steps. Mathematically speaking, such a trajectory is a path through a 6n-dimensional space if n particles have been simulated. Such a high-dimensional object is not suitable for direct evaluation—methods for reducing the dimensionality are required. These methods will be termed trajectory analyses here.

It is a commonly known phenomenon that many working groups from the field of computational science use their own in-house scripts and codes to evaluate and analyze simulation trajectories. These codes all have their specific scopes, strengths, and limitations, and despite being developed in independent groups, many methods and algorithms are for sure almost equivalent in all these tools. This comes with several disadvantages. On the one hand, it wastes a lot of valuable scientific manpower by “reinventing the wheel” hundreds of times. On the other hand, it increases the risk of subtle scientific or coding problems not being detected because the user base of the code is often only one working group and therefore rather small. Hence, it is highly desirable to have common tools for analyzing trajectories, which are utilized by a large group of scientists.

In 2011, we published the TRAVIS (“Trajectory Analyzer and Visualizer”) program package1 as an open source free software, which specifically aims at providing as many trajectory analyses as possible within a unified program so that many scientists are hopefully no longer forced to develop their own in-house solutions to analyze trajectories. There already exist several other software packages for analyzing trajectories, e.g., the Visual Molecular Dynamics (VMD) program package2 or Python frameworks such as MDAnalysis.3 However, by its huge number of available analyses and by its interactive command-line user interface, TRAVIS is in a certain sense unique. Since the release in 2011, the original TRAVIS article1 has already been cited more than 400 times by more than 800 different co-authors (according to the journal homepage), and around 95% of these publications actually used the software to produce results in the manuscript. This gives us the impression that our code is somehow useful for the community.

In this article, we want to highlight some of the recently added methods and features in TRAVIS. At the same time, we are going to shortly introduce some of the underlying algorithms that are typically invisible for users but are, of course, very important to achieve the desired efficiency and accuracy of trajectory analyses. This manuscript is structured as follows: After a discussion of general features of TRAVIS in Sec. II, we will discuss specific features from the fields of structural, dynamical, and spectroscopic analyses in Secs. III–V. In this context, structural analyses are defined as analyses that are invariant under any re-shuffling of the trajectory frames, while dynamical analyses somehow depend on the ordering of the frames. Of course, spectroscopic analyses are also dynamical analyses, but they are discussed separately. We end our article with some conclusions.

Being licensed under the GNU GPL v3 license, TRAVIS is an open source free software and can be downloaded from http://www.travis-analyzer.de/.

Apart from the source code, Windows executables can also be obtained from the homepage. New versions, including bug fixes and additional features, are regularly supplied. TRAVIS is written in C++, and the developers’ version currently contains around 270 000 lines of source code. It does not require any external libraries or dependencies (except from a C++ compiler) and is platform independent. TRAVIS is currently not parallelized, i.e., it runs on one single central processing unit (CPU) core. However, we are working on introducing optional OpenMP-based parallelization for some suitable modules.

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-orthorhombic cell geometries and variable cell vectors (NpT ensemble) are supported. The cell geometry can directly be read from the trajectory if present. Very large systems with several hundred thousand of atoms can be easily handled.

TRAVIS can be used as an interactive program, asking questions to the user in a text-mode console interface. Most of the questions are more or less self-explanatory, and reasonable default values are provided in most cases. On the other hand, TRAVIS can also be used for batch processing with the help of input files. When creating plots and visualizations, TRAVIS does not directly output image files, but rather input files for visualization programs such as VMD,2 xmgrace,8 Gnuplot,9 Povray,10 and Mathematica.11 By supplying these files to the corresponding programs, high-quality images can be easily obtained without requiring much knowledge of the visualization program.

TRAVIS only reads atom labels, positions, and optionally velocities from the input trajectory; any topological information (molecules, bonds, etc.) is ignored. Instead, a distance-based bond detection is executed in one trajectory frame, and molecules are defined on the basis of the detected covalent bonds. A pair of atoms is defined to be covalently bonded if its distance is smaller than the sum of the covalent radii of the two atom types (taken from the literature)12 multiplied with a constant to account for thermal vibrations (by default 1.15, but can be chosen differently). To detect bonds, a straight-forward (“canonical”) approach of checking the distance of all pairs of atoms in the system is very inefficient, as its computational time scales with On2. TRAVIS uses a more sophisticated approach, which is based on a spatial domain decomposition. Domain decompositions and similar approaches such as the linked list method13 are commonly used in computational chemistry (in particular, in parallelized force field molecular dynamics codes such as LAMMPS5and Gromacs14) but have rarely been applied to trajectory post-processing. First, the longest possible bond distance is determined by finding the element in the system with the largest covalent radius, and taking this radius twice, multiplied by the constant mentioned above. Then, the simulation cell is tessellated into cubes of equal size and shape. The number of cubes along each cell vector is chosen as large as possible while still satisfying the constraint that the edge length of the cube is larger than the longest possible bond distance. All atoms in the system are then binned into these cubes. To find all covalent bond partners of a certain atom inside of such a cube, only atoms within the same cube and within the 26 neighboring cubes need to be considered. This yields a very fast bond detection algorithm that essentially scales with On.

The efficiency of this approach is easily visible in Table I, where the molecule recognition times for similar systems (containing ion pairs of the ionic liquid [EMIm][OAc], i.e., 1-ethyl-3-methylimidazolium acetate) of increasing size are shown. While the fourth column gives the molecule recognition time with the canonical algorithm, the fifth column depicts the timing of the domain decomposition approach. For relatively small systems, both methods are acceptable. However, for a system with around 500 000 atoms in a 17 nm cubic cell, the canonical approach takes almost 1 h for the molecule recognition, while the implementation based on domain decomposition only runs for 1 s.

TABLE I.

Execution times of the molecule recognition for different system sizes with canonical implementation (tCan) and domain decomposition implementation (tDom). The snapshots contain ion pairs of the ionic liquid [EMIm][OAc].

Ion pairsAtomsCell (nm)tCan(s)tDom (s)
36 936 2.121 0.003 0.001 
288 7 488 4.242 1.036 0.001 
2 304 59 904 8.485 83.009 0.011 
18 432 479 232 16.969 3240.000 1.002 
Ion pairsAtomsCell (nm)tCan(s)tDom (s)
36 936 2.121 0.003 0.001 
288 7 488 4.242 1.036 0.001 
2 304 59 904 8.485 83.009 0.011 
18 432 479 232 16.969 3240.000 1.002 

After the molecule recognition, topological atom codes are assigned to the atoms by using Shelley’s extension15 of the Morgan algorithm.16 The algorithm guarantees to give identical codes to topologically equivalent atoms and different codes otherwise. Atom numbers are assigned by descending atom codes. As a by-product, TRAVIS maintains a list of topologically identical atoms, which is helpful for many analyses. Furthermore, a full ring system detection and refinement17 is carried out, as some of the analyses are well-suited for cyclic structures. Finally, all individual molecules with identical sets of atom codes are considered to be topologically equivalent and grouped together as molecule types. We would like to point out again that this full protocol only takes 1 s on one CPU core for a system with around 500 000 atoms and 40 000 molecules, as shown in Table I.

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 (“losslessindeed means that the input Gaussian Cube file is bit-wise reproduced). The methods in TRAVIS for computing vibrational spectra (see below) require the total electron density on a grid along the trajectory, easily amounting to several terabytes of raw data. With our novel compression scheme, these data can be compressed significantly without losing any information. Figure 1 gives some details on the performance of the method. While multi-purpose compression 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 compression ratio of around 40:1 while still being faster than the former two methods in compression and decompression time. Furthermore, the BQB format allows for random access to the frames, which is, of course, not the case with multi-purpose compression formats. The high compression efficiency is achieved by applying a sequence of methods, including multi-dimensional polynomial extrapolation, Hilbert curve re-ordering, and canonical multi-table Huffman encoding.7 

FIG. 1.

Compression times (orange), decompression times (green), and compressed sizes (blue) of a volumetric electron density frame with different lossless compression algorithms. The BQB format7 discussed here is shown on the right.

FIG. 1.

Compression times (orange), decompression times (green), and compressed sizes (blue) of a volumetric electron density frame with different lossless compression algorithms. The BQB format7 discussed here is shown on the right.

Close modal

While our original aim was to compress volumetric data trajectories, we noted that the same methodology is also able to losslessly compress standard position trajectories. When starting from a XYZ trajectory file with a frame interval of 0.5 fs, the BQB format reaches a compression ratio of around 20:1 here. This is much more efficient than other existing trajectory formats such as DCD or Gromacs XTC14 while not sacrificing accuracy. More details can be found in Ref. 7.

The BQB format and the methods described above are implemented in TRAVIS, which can be used to compress or decompress volumetric or position trajectories, or directly read and analyze such trajectories. Apart from that, there is also a stand-alone tool for compressing and decompressing trajectories, called the “bqbtool.” A library (“libbqb”) to add support for the BQB format to other program packages will be released soon. For more information, refer to the website of the BQB format.18 

In this section, some structural analyses implemented in TRAVIS will be presented. As mentioned above, structural analyses are defined as analyses that are invariant under any re-shuffling of the trajectory frames in this context. Due to this property, structural analyses can also be applied to Monte Carlo simulations because they do not require that a physical time elapses along the trajectory.

Probably, the most common kinds of trajectory analyses are particle density histograms. In particular, the radial distribution function (RDF), also known as gr, is found in countless publications 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 function that automatically recognizes the first maximum and the first minimum of the RDF and prints out the distance as well as the intensity of these two points in the log file together with the particle number 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 computed 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 angular distribution functions (ADFs) and dihedral distribution functions (DDFs), respectively. Furthermore, the distance between a point and a line can be monitored, where the point can be any atom in some molecule and the line can be defined by any other two atoms in the same or a different molecule. The corresponding analysis is called point–line distribution function (LiDF). Similarly, the distance between a point and a plane can be observed, where the plane can be defined either by three points in some molecules or by specifying 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. As the resulting particle density is a volumetric dataset, only iso-surfaces of these data can be visualized. TRAVIS offers to output the resulting volumetric data as a particle density (nm−3) or relative to the uniform density, similar to RDFs. Furthermore, there are some advanced functions implemented such as smoothing SDFs, inverting SDFs, creating cut planes, and forcing SDFs to be symmetric to certain planes to improve sampling. There are many examples where SDFs computed by TRAVIS have been published in the literature.20,21,24,25,32,39–48

If the simulation cell is non-isotropic (e.g., systems with a surface or lipid bi-layer), it can be interesting to compute density profiles (DProf) along a fixed vector in the simulation cell. The result is a histogram that gives the particle density of a selected particle type (either in nm−3or relative to uniform density) in thin slices of the system perpendicular to the chosen vector. Density profiles computed by TRAVIS have already been published.29,49,50

Recently, we implemented a function called geometric density analysis (GeoDens) to observe the temporal development of the average particle density within a geometric region along the trajectory. In the specific application, we defined a cylinder between the two amino acid molecules involved in a salt bridge.31 The cylinder has a fixed radius, but the starting and ending points are tied to the molecules so that the cylinder moves together with the molecules. By monitoring the particle density of anions within this cylinder relative to the uniform density, we were able to conclude that the anions are repelled from the salt bridge (see Fig. S-18 of Ref. 31).

If two observables are correlated, the individual histograms of the observables do not give a complete picture of the situation, and a multi-dimensional histogram is required to understand the correlation effects (for example, see Fig. 2). Within a liquid phase trajectory,23,24 the length of a certain hydrogen bond is depicted on the horizontal axis, while the hydrogen bond angle is shown on the vertical axis. The 2D contour plot possesses a very strong signal, indicating a hydrogen bond, at low distances and angles of around 180°. The plot gives the information that at small distances, the angle is preferentially found to be close to 180°. This is a correlation between both quantities, and the standard 1D histograms would not be able to give this information.

FIG. 2.

Combined distribution function (CDF) depicting the O⋯H distance and the O⋯H–C angle in the ionic liquid [EMIm][OAc].23,24 A possible geometric hydrogen bond criterion is shown by the yellow rectangle.

FIG. 2.

Combined distribution function (CDF) depicting the O⋯H distance and the O⋯H–C angle in the ionic liquid [EMIm][OAc].23,24 A possible geometric hydrogen bond criterion is shown by the yellow rectangle.

Close modal

TRAVIS combines a multitude of scalar quantities into multi-dimensional histograms, which we call combined distribution functions (CDFs) here. By using combinations of distances, angles, dihedral angles, point–plane distances, and point–line distances, thousands of interesting CDFs can be computed from a single trajectory. There are many examples in the literature of CDFs computed with TRAVIS.21–26,30,33,37,40–42,46,51–54 It is even supported to combine three scalar quantities into a three-dimensional histogram in TRAVIS. The resulting histogram can be visualized by means of iso-surfaces, as it has been shown in the literature a few times.24,33,55

An even more direct way for understanding correlation between quantities is the so-called correlation plot in TRAVIS. Consider Fig. 3 where a CDF with two distances is shown in the left panel. While the horizontal axis depicts the C–O bond length in an acetate anion within a bulk phase simulation,23,24 the vertical axis shows the length of the other C–O bond in the same acetate anion. If both quantities would be uncorrelated, the 2D histogram would simply be the Cartesian product of the two 1D bond length histograms. To reveal the correlation, TRAVIS internally creates this Cartesian product of the two histograms and subtracts it from the actual CDF. The resulting correlation plot is shown in the right-hand panel of Fig. 3. Positive values indicate that the probability of finding this configuration is larger (i.e., positive correlation) than if the two quantities would be uncorrelated, while negative values depict situations less probable than in the uncorrelated case (negative correlation). We find that configurations with one elongated and one shortened bond are more likely than expected from the individual bond length statistics, while cases where both bonds are either long or short are less likely—a clear correlation effect.

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).

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).

Close modal

A very important experimental technique for condensed phase systems is measuring structure factors, also called Sq. The structure factor is related to gr of the system by a Fourier transform,

Sq=1+ρeiqrgr1dr.
(1)

Despite this easy relation, it is slightly more involved to compute Sq 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 the corresponding element types, and finally summed up. TRAVIS does this automatically and directly outputs the total Sq, together with all partial structure factors (i.e., contribution to the total structure factor from a certain pair of elements).27 Several commonly used normalization factors are available. Furthermore, it is possible to compute individual intramolecular and intermolecular contributions to the structure factor. Neutron scattering cross sections for the most important isotopes56 as well as x-ray scattering cross sections (given by Cromer–Mann coefficients) are included in TRAVIS, so both neutron and x-ray structure factors can be computed.27 Recently, the option to apply a Lorch-type window function57 to the structure factor calculation was implemented. Structure factors computed by TRAVIS have been published in the literature several times.27,52,55,58–62 TRAVIS can also compute dynamic structure factors Sq,ω, which is discussed in Sec. IV E.

There exists a large variety of so-called order parameters to determine perturbed or ordered structures. TRAVIS only contains a few of them, which are suitable to characterize lipid bi-layers that we have been simulating.29 The first one is the tail vector tilting analysis. If a lipid bi-layer is in the sol phase, the fatty acid tail vectors are oriented in random directions and constantly fluctuating. However, when the bi-layer undergoes a transition into the gel phase, the tail vectors are all tilted with the same angle into the same direction. By projecting all tail vectors into the X–Y plane (where the Z direction is the normal vector of the bi-layer), it can be easily seen if the bi-layer is in the sol or gel state (see Fig. 3 of Ref. 29).

The second order parameter that is implemented in TRAVIS is SCD, also called the deuterium order parameter, as it can be directly measured by nuclear magnetic resonance (NMR) experiments of deuterated samples.63,64 It is computed by the following equation:

formula
(2)

where k depicts the position in the alkyl chain, N is the number of lipid molecules, T is the total number of time steps, rikCHt is one of the C–H (or C–D) bond vectors of the kth carbon atom of the chain in time step t, and n is the membrane normal vector. SCD can take values between −0.5 and 1, where 0 corresponds to the statistical distribution (i.e., no ordering) and −0.5 indicates that the bond vectors are strictly within the membrane plane. An application of this analysis can be found in Fig. 4 of Ref. 29.

Finally, we implemented a function to compute the gauche ratio along alkyl chains, i.e., the average fraction of C–C–C–C dihedral angles at a given chain position being close to 60°, which is called the gauche configuration. This configuration is energetically less favorable than the anti configuration with dihedral angles of around 180° and therefore also statistically less populated. For an application of this analysis and interpretation of the results, see Fig. 5 and Table 2 of Ref. 29.

In condensed phase simulations of solutions, it is often very interesting to understand how the solvent molecules are arranged and oriented around the solute. Information on the solvent shells can be gained from RDFs and SDFs, but no details on the orientation of the individual solvent molecules are gained from these analyses. CDFs involving distance and angle can give information on the solvent orientation, but it is still not apparent what the preferred solvent orientation at a specific position around the solute molecule is. To overcome this, we implemented the plane projection analysis (PlProj). Consider Fig. 4 where an [EMIm]+ cation from a simulation of one ion pair [EMIm][OAc] in water is fixed. The color scale of the plot depicts the average particle density of water at the given position relative to uniform density (considering only a thin slice of the simulation cell defined by the [EMIm]+ring plane). A value larger than 1 means that water is more frequently found here than at some arbitrary position, while a value below 1 indicates a depletion of water molecules at that point. It can be seen that there exists a well-defined first solvation layer with particle densities above 1, which is separated from the water bulk phase by a clear region of depletion with particle densities below 1. The vectors in the plot depict the average orientation of water dipole vectors at each point. As expected, it can be seen that water forms hydrogen bonds to the three protic ring hydrogen atoms, and therefore, the water dipole vector points away from these atoms on average. Next to the non-polar regions of the molecule, there is no clear preferential orientation of the water molecules, resulting in very short average vectors.

FIG. 4.

Average particle density of water molecules around an [EMIm]+ cation (color scale) together with the average orientation of the water dipole vectors (arrows) in a simulation of one ion pair of [EMIm][OAc] in water.

FIG. 4.

Average particle density of water molecules around an [EMIm]+ cation (color scale) together with the average orientation of the water dipole vectors (arrows) in a simulation of one ion pair of [EMIm][OAc] in water.

Close modal

Many molecular liquids form hydrogen bonds, which often have a significant influence on the properties (just consider water). If there are several different hydrogen bond donors and acceptors in a complex system, a subtle balance of interactions is found, and even small changes in hydrogen bond donor or acceptor strengths can lead to large effects. Therefore, it is very important to understand the hydrogen bond topology and the competition for donor/acceptor positions in such systems. A simple estimate for the strength of a hydrogen bond is the distance and height of the first maximum in the X⋯H RDF, together with the number integral up to the first minimum, which gives the average number of hydrogen bonds of atom X. Let us consider a system where a cellulose oligomer is dissolved in a mixture of the ionic liquid [EMIm][OAc] and water. There are three hydrogen bond donors (cellulose alcohol H, water H, and [EMIm]+ ring protons) and three acceptors (cellulose O, water O, and [OAc] O). The hydrogen bond estimates from the pairwise RDFs could be written as a matrix of numbers, but this is not very illustrative. Recently, we found a much more appealing way to represent such data, which is based on Sankey diagrams.65 The concept is more than 120 years old now, but it is rarely found in computational chemistry and never has been used to visualize hydrogen bond topology to the best of our knowledge. Therefore, despite its early origin, we consider this a “modern” visualization technique. Some common plotting libraries such as MatPlotLib66 are able to produce Sankey diagrams, but not the circular variant described below. For a visualization of the hydrogen bond topology of the mixture described 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 correspond 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 represent such topologies, and we hope to see them in many future computational chemistry papers from now on. As no common plotting program is able to produce both linear and circular Sankey diagrams as shown, TRAVIS directly outputs the images as SVG vector graphics.

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.

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.

Close modal

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 acceptors. For example, if it would be of interest which of the 16 different 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 in the system, while the columns represent the different hydrogen bond donors. For each of the matrix elements, TRAVIS internally computes an RDF and extracts the distance and height of the first maximum. If the RDF does not indicate the presence of a hydrogen bond, the matrix element is simply filled with a black cross. If a hydrogen bond exists, the square is filled with a color according to a two-dimensional color scale that is shown on the right-hand side of Fig. 6. Both the distance and the height of the first RDF maximum are encoded in color. For example, red indicates a hydrogen bond with a very small distance and a very large maximum height (i.e., a very strong hydrogen bond). Yellow, on the other hand, corresponds to a hydrogen bond with a large maximum height but large distance, which is therefore not considered very strong. From looking at the matrix, it is immediately visible that the acetate anion forms strong hydrogen bonds to many donors, while many cellulose oxygen atoms are not involved in strong hydrogen bonds. There are two strong intramolecular hydrogen bonds clearly visible in cellulose, namely, O3⋯H25 and O8⋯H31. Obtaining this information just with standard RDF analyses would have required to compute and manually look at many hundred RDFs. We are convinced that 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).

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).

Close modal

Despite being more than 110 years old now, the Voronoi tessellation67 is still a fascinating mathematical formalism. When applied to a simulation cell with the atoms as Voronoi sites, in simple words, it partitions the box into Voronoi cells belonging to each atom in 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 cell. TRAVIS uses the free Voro++ library68,69 (integrated in the TRAVIS source code so that it does not have to be provided externally) to perform a periodic radical Voronoi tessellation.70,71 From such a tessellation, many interesting properties can be derived. On the one hand, Voronoi cells can be directly visualized by a 3D rendering via Povray.10 On the other hand, one can create statistics on the Voronoi cells’ volumes, surface areas, isoperimetric quotients (measures for the “sphericity” of a cell), and so on. It is also possible to monitor Voronoi surface coverage, i.e., which percentage of the faces of a molecule’s Voronoi cell is shared with other molecules of a certain type. Both the visualization and the statistics can be performed with TRAVIS, as already demonstrated in the literature.37,42,72,73

Another interesting result from a Voronoi tessellation is a simple and unbiased neighborhood criterion: If the Voronoi cells of two atoms share a common face, these two atoms are said to be neighbors. In contrast to other neighborhood criteria, e.g., based on the first minimum of the RDF, this approach convinces through its great simplicity, as it does not possess any empirical parameters. It is therefore an unbiased criterion and transferable to all kinds of systems. In Fig. 7, the neighborhood probabilities for all pairs of atoms in a bulk phase simulation of [EMIm][OAc]23,24 are visualized as a matrix. White corresponds to a neighborhood probability of zero, i.e., these two atoms never share a common Voronoi face, while purple represents a neighborhood probability of one, i.e., the one atom type always has at least one shared Voronoi face with an atom of the other type. The strong hydrogen bonds between acetate oxygen atoms and the protic ring hydrogen atoms can be directly seen, while some interaction between the hydrogen atoms of the aliphatic side chains is also visible. Such a neighborhood matrix captures the full neighborhood relationship of the system in one image and can therefore be considered a fingerprint of the liquid structure.

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 share a common Voronoi face at a time in a simulation of [EMIm][OAc].

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 share a common Voronoi face at a time in a simulation of [EMIm][OAc].

Close modal

To study microheterogeneity and microphase separation in complex liquid systems, we developed the Voronoi-based domain analysis (DomA) in 2015.28 First, the atoms of the system are categorized into different classes. This can be, e.g., polar, non-polar, side chain, aromatic, fluorinated, etc. Then, a radical Voronoi tessellation of the simulation cell is performed using all atoms of the system as Voronoi sites. The Voronoi cells of atoms from the same class, which are neighbors, i.e., share a common Voronoi face, are merged. This leads to domains of different sizes and shapes. A two-dimensional illustration of this concept is presented in Fig. 8, where the different fill colors represent the different atom classes and the bold black lines correspond to the domain boundaries. While in a perfectly homogeneous system, there will be many small domains, a certain degree of microheterogeneity will lead to a smaller number of larger domains because the atoms from the same class tend to form clusters. From performing this analysis in all trajectory frames and looking at the domain statistics (average domain count, volume, surface area, isoperimetric quotient, etc.), one can gain much insight into the microphase behavior of the system. This has already been successfully applied many times in the literature.27,60,73–78

FIG. 8.

Schematic representation of the Voronoi-based domain analysis.28 Colors represent different classes of atoms. Domains are highlighted by bold black lines.

FIG. 8.

Schematic representation of the Voronoi-based domain analysis.28 Colors represent different classes of atoms. Domains are highlighted by bold black lines.

Close modal

In a solution, the exclusion volume of a solute perturbs the structure of the pure liquid. Obviously, the required amount of energy for this perturbation is decreased if there already exist void cavities in the pure liquid, which are suited to accommodate the solute. Therefore, it can be expected that the size, shape, and mobility of void cavities play a crucial role in many fundamental processes of solvation and mass transport in fluid systems.

With the aim to quantify this effect, we have implemented a Voronoi-based void analysis in TRAVIS, where a set of spheres is grown into the voids between the molecules in a simulation cell (seeFig. 9).79 In contrast to other algorithms,80,81 we enhance this by the definition of a void domain: The spheres are handled as pseudo-atoms and the corresponding Voronoi cells are combined (cf. the domain analysis). With the aid of the isoperimetric quotient and autocorrelation functions, information on the shape and mobility of the void space can be extracted.

FIG. 9.

Example of the detected void spheres in a simulation box containing N-methylimidazole and acetic acid (see Ref.82for details). Several single spherical voids are combined into a much larger void space domain.

FIG. 9.

Example of the detected void spheres in a simulation box containing N-methylimidazole and acetic acid (see Ref.82for details). Several single spherical voids are combined into a much larger void space domain.

Close modal

Most of the particle density histograms discussed above only involve two molecules at a time. Any RDF is based on the distance between two atoms. In CDFs, there can be more atoms involved (e.g., when observing a distance and an angle, as shown inFig. 2), but there are still only two molecules involved at a time. One can easily think of analyses that involve more than two molecules.

For example, let us consider a simulation of a concentrated 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 atom is not further than 250 pm away from a Na+ ion. This analysis involves three particles at a time (one water molecule, one Na+, and one Cl). A straightforward approach would loop over all possible 3-tuples of these particles, skipping all tuples that do not fulfill the conditions and using the remaining tuples to populate the multi-dimensional histogram. However, this is an inefficient algorithm, as most of these tuples will not fulfill the conditions, and it is slow to compute so many unnecessary distances. A more optimized approach would first define a maximum observation range. All distances beyond this range are neglected. Then, it would utilize a spatial domain decomposition, as explained for the molecule recognition in Sec. II A, using the maximum observation range as minimal cube diameter. Now, for a particular water molecule, only the ions within the central cube and the 26 surrounding cubes need to be considered for forming three-tuples with that water molecule, excluding all the ions that would not fulfill the conditions anyway. This leads to an implementation that essentially scales linearly with system size.

We are currently working on an unified formulation for all such analyses, which we call the generalized tuple analysis (GenTup). Based on a central particle, as many particles as required can be 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 maximum observation range around the central particle (the maximum observation range can be chosen by the user). It will be possible to use many different condition relations between the particles (distances, angles, dihedral angles, triangle surface area, tetrahedron volume, etc.). To the best of our knowledge, such a generalized analysis has not yet been published in the literature. Our implementation in TRAVIS will be available in some months.

In this section, some dynamical analyses contained in TRAVIS will be shortly discussed. As already described above, dynamical analyses are defined as analyses that are not invariant under re-shuffling of the trajectory frames, i.e., they depend on the ordering of the frames in some way. Most dynamical analyses should not be applied to Monte Carlo simulations, as there exists no physical time in standard Monte Carlo approaches.

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 valuable property for validating simulations. The diffusion coefficient D of some molecule can be computed from the mean square displacement (MSD) via the Einstein relation83 

D=12dlimτptpt+τ2tτ,
(3)

where pt is the position vector of a particle at time t and d represents 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 coefficient. Some years ago, a fast method for computing MSDs via Fourier transform has been proposed in the literature;84 we are currently 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

Another approach to diffusion coefficients of simulations is the Green–Kubo diffusion relation,88 which relates the integral of the velocity autocorrelation function to the diffusion coefficient. It reads

D=1d0vtvt+τtdτ,
(4)

where vt is the velocity vector of the considered particle at time t and d again stands for the dimensionality of the system. In addition, the computation of velocity autocorrelation functions is implemented in TRAVIS so that results from both approaches can be compared (they are typically not identical for finite-length trajectories, and it is unclear which one is “better”).

Another very interesting dynamical property is the reorientation dynamics of certain vectors in a simulation, also known as rotational relaxation time. It can directly be determined by certain NMR and electron paramagnetic resonance (EPR) techniques and therefore offers a good opportunity for comparison between experiment and simulation. Based on a trajectory, this quantity can easily be expressed as an integral over the corresponding vector autocorrelation function Cτ. For comparison with some experiments, it is desirable to apply a Legendre polynomial Pn of order n to the vector dot product in the correlation, then termed nth order correlation function Cnτ and nth order reorientation time Tn. This leads to the equations

Cnτ=0Pnutut+τutut+τdt,
(5)
Tn=0Cnτdτ,
(6)

where ut is the chosen vector in the system at time t. As the first-order Legendre polynomial P1 is the identity, the reorientation dynamics for n = 1 equals the case without the Legendre polynomial. For typical simulation lengths, the vector autocorrelation function Cτ does not decay to zero so that the integration to obtain T cannot be directly performed. To overcome this, TRAVIS automatically fits a poly-exponential function to Cτ via a least-squares procedure (using a Levenberg–Marquardt minimizer),89 which is then analytically integrated to determine T as a result. Vector reorientation times computed with TRAVIS have already been published in the literature.20,31,58,90–92

In computer simulation, it is often desirable to determine the average lifetimes of certain aggregates (e.g., ion pairs or hydrogen bond complexes). It is for sure not a valid approach to simply observe the lifetime of such a complex once in a trajectory—a molecular dynamics simulation is a chaotic system, and by chance, any arbitrary lifetime can be obtained from such an “experiment” with some non-zero probability. A valid approach is to run a simulation until the investigated aggregate has formed and broken sufficiently many times and then use the autocorrelation formalism to determine its lifetime.93,94 To do so, first a geometric criterion for aggregation is chosen. This can be a simple distance criterion or a more complex variant. For example, when defining hydrogen bonds via a geometric criterion, it has shown to be effective to define a criterion for both the distance and angle (see the yellow rectangle inFig. 2for a reasonable hydrogen bond criterion). Then, an auxiliary function

βijt1if criteria fulfilled0otherwise
(7)

is defined, where βijt=1 exactly when an aggregate between molecules i and j at time t is intact according to the chosen geometric criterion. From this, the autocorrelation

cτ=1N2i=1Nj=1N0βijtβijt+τdt
(8)

is computed. The lifetime T of the aggregate can then be obtained as twice the total integral of this autocorrelation function

T=20cτdτ.
(9)

The factor of 2 can be understood from the following simple consideration: Imagine an aggregate that existed only once for a time interval a. The corresponding correlation function cτ is then a piecewise-linear function that runs through the two points c0=1 and ca=0. The integral of this function will be a2. A factor of 2 is required to retain the original lifetime of a.

In practice, the correlation function cτ will often not decay to zero within the simulation time. As explained in Sec. IV B, TRAVIS automatically fits a poly-exponential function to the correlation function via a least-squares procedure (using a Levenberg–Marquardt minimizer),89 which is integrated analytically to obtain the lifetime of the aggregate. Furthermore, in finite-size simulation cells, the autocorrelation function never decays to zero because any two molecules that formed an aggregate will eventually meet again. To obtain correct lifetimes, it is important to subtract the ensemble average value from the autocorrelation function before performing the exponential fit and the integration. TRAVIS automatically performs this correction.

The above-description corresponds to so-called intermittent autocorrelation functions. TRAVIS can also compute continuous autocorrelation functions, which are not discussed here. Aggregate lifetimes and autocorrelation functions computed by TRAVIS have been published in the literature several times.20,22,25,37,39,40,43,47,77,78

1. Efficient implementation

Computing aggregate lifetimes as described above would be highly inefficient. Imagine a liquid water force field trajectory with 1000 molecules and 106 trajectory steps. Storing the value of the βijt function for all pairs and all time steps would require 1012 memory cells, i.e., 1 TB of RAM if every function value is stored in 1 byte. A much better approach is to keep a list of aggregates. For any aggregate between two molecules, it is only stored which two molecules are involved in that aggregate and from which to which time step the aggregate was intact. Such a list requires only a tiny fraction of the memory described above, and for every pair i and j, the function βijt can be reconstructed from the list. By this approach, only the data for one pair of i and j have to be stored in 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 βijt can only take the values 0 and 1 and can therefore be written as the sum of rectangular functions. Let

Rabt1at<b0otherwise
(10)

be a rectangular function. The cross-correlation of two such functions Rabt and Rcdt has a simple piecewise-linear form (a trapezoid)

RabtRcdt+τt=τEFE,Eτ<F1,Fτ<GτHGH,Gτ<H0,otherwise
(11)

with the substitutions

Ecb,Fminca,db,Gmaxca,db,Hda.
(12)

As βijt can be written as a sum of rectangular functions (by knowing the aggregate lifetime intervals from the list), the resulting product in the autocorrelation can be expanded (“multiplied out”) to a sum of trapezoid functions as shown above. Then, the autocorrelation 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.

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 definition, 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

n(τ)=0τH(t)[1h(t)]ḣ(0)hdt,
(13)

which gives the probability of finding a pair that was initially interacting 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

ċ(t)=kdc(t)kfn(t).
(14)

The inverse of the rate constant of decay kd is defined as the average lifetime τlt=kd1 of the interaction.

The current implementation in TRAVIS30,99 enables the analysis for hydrogen bonding—based on a geometrical criterion as visualized in the left part of Fig. 10—as well as for ion pair dynamics based on a next neighbor criterion and a distance, as illustrated in the right part of Fig. 10. The analysis of the interaction dynamics allows a deeper insight into the molecular cosmos of chemical systems,100 especially in organocatalysis,101 electrolytes,102 and metal extractants.77 

FIG. 10.

Criteria for the reactive flux analysis. (Left) Hydrogen bond analysis: aggregation if distance rAH and angle α are below a cutoff, vicinity defined via the distance rAD. (Right) Ion pair analysis: aggregation defined via the next neighbor condition, vicinity defined via the distance rIC of the first solvation shell.

FIG. 10.

Criteria for the reactive flux analysis. (Left) Hydrogen bond analysis: aggregation if distance rAH and angle α are below a cutoff, vicinity defined via the distance rAD. (Right) Ion pair analysis: aggregation defined via the next neighbor condition, vicinity defined via the distance rIC of the first solvation shell.

Close modal

As described above, an RDF is computed by taking the positions of two particles in the same trajectory frame, measuring the distance between the two positions, and adding the result to a histogram. It is interesting to take the positions of the two particles from different trajectory frames with a “lag time” of τ between them. If this is performed for different lag times and represented as a 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 correlation function for the O⋯H distance in a simulation of liquid water is presented. It can be seen that the standard RDF is reproduced at τ = 0, but the features in the RDF quickly are blurred out for an increase in lag time. Van Hove correlation functions are particularly interesting because their two-dimensional Fourier transform is the dynamic structure factor Sq,ω, which can be measured by neutron scattering experiments. This offers another opportunity to directly compare experimental and computational results.

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 function corresponds to the dynamic structure factor Sq,ω.

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 function corresponds to the dynamic structure factor Sq,ω.

Close modal

It is known since a long time that vibrational spectra can be computed from molecular dynamics simulations by time correlation functions.104 This comes with several advantages over the widely used static–harmonic approach to spectra. Automatic conformer sampling takes place during the MD simulation, the full solvent effect can be included, realistic band shapes are directly obtained, and even some anharmonic effects (such as temperature-dependent frequency shifts as well as approximate overtones and combination bands) are covered.105 During the last few years, many functions for predicting vibrational spectra have been implemented in TRAVIS. This section will give an overview on these methods and techniques. A detailed step-by-step tutorial on computing vibrational spectra with CP2k106 and TRAVIS can be found on the Internet.107 

This section describes some general techniques used in TRAVIS, which are relevant for the calculation of all types of vibrational spectra. While most of these techniques are basic knowledge of signal processing, they are rarely discussed in articles on computer simulations.

1. Fast correlation

All prediction methods for vibrational spectra from molecular dynamics simulations are based on either autocorrelation or cross-correlation functions of some time series along the trajectory. Typically, these correlations are computed for molecular properties, not for the properties of the whole simulation box. The reason is that correlating properties of the whole system is equivalent to taking into account all cross-correlations between different molecules. However, the motion of very distant molecules is not correlated at all so that these cross-correlations mostly introduce noise into the spectrum. Computing all the molecular correlation functions requires much computational power because each such calculation scales with Omn, where n is the length of the time series and m is the chosen correlation depth. Fortunately, there exists the Wiener–Khintchine theorem,95 which states that computing the autocorrelation Cτ of a time series ft is equivalent to Fourier transforming the time series, taking the absolute square, and transforming back,

Cτ=ft¯ft+τdt=F1Ffτ2.
(15)

As the Fourier transform of a discrete dataset of length n can be efficiently computed by the fast Fourier transform (FFT) method that only scales with Onlogn, this saves a lot of computational time. Equation (15) is only a special case of the cross-correlation theorem so that a similar approach can also be used to speed up the calculation of cross-correlation functions via FFT. TRAVIS computes all autocorrelation and cross-correlation functions based on these approaches using the KISS FFT library108 (integrated in the TRAVIS source code).

2. Window function and zero padding

When considering Fourier transforms of data, it is very important to choose a suitable window function. There exists no Fourier transform without a window function—using no window function means using a rectangular window, which is a very poor choice. Applying a window function simply means element-wise multiplication of the discrete input data with the window function before Fourier transforming. By default, TRAVIS uses the Hann window function109 for computing spectra, which is given by

Wn=cos2πn2N1,
(16)

but other window functions (exponential, Gaussian) are also implemented. The total integral of the Fourier transform of a discrete dataset ai depends only on the first value a0. Therefore, any window function that fulfills the condition W0 = 1 does not modify the total integral at all. In other words, window functions alter the peak shapes of the spectrum but maintain the integral intensities of the bands.

By adding zeroes to the tail of the correlation function before Fourier transforming, the resolution of the resulting spectrum can be increased, which is called zero padding. No new spectral information is gained by this technique, it is just an interpolation between the data points of the original spectrum, i.e., a cosmetic measure. By default, TRAVIS uses a zero padding factor of 4.

3. Finite difference correction

For most spectra, the time derivative instead of the quantity itself enters the correlation functions. For a cosine function cosωt+φ, it can be shown105 that the second-order central finite difference derivative is given by

Dcosωt+φ=cosωt+Δt+φcosωtΔt+φ2Δt=ωsinωt+φsinωΔtωΔt,
(17)

i.e., the exact value of the derivative is multiplied by a sinus cardinalis (sinc) function sinωΔtωΔt. Since the product of two such time derivatives is formed in the correlation functions, the final spectra need to be divided by a factor of sinωΔtωΔt2 to obtain the correct intensities. This correction is implemented in TRAVIS and active by default. Due to the correction, high-quality spectra with accurate intensities can be obtained even if the dipole moment only was computed every 4.0 fs.

4. Improved sampling via time reversibility

Similar to most other properties, vibrational spectra should be invariant under reversal of the time direction: The trajectory with reversed time samples the same thermodynamic ensemble as the forward trajectory. As vibrational spectra are ensemble averages, they cannot differ for these two cases. Consider the cross-correlations of two time series at and bt along the trajectory. Let ãtaTt and b̃tbTt be the two time series of the trajectory with reversed time, where T is the total trajectory length. We find that

ãt+τb̃tt=ãt+τb̃tdt=aTtτbTtdt=at̃bt̃+τdt̃=bt̃+τat̃t̃
(18)

with the substitution t̃Ttτ. Thus, the reversal of the time direction is equivalent to swapping the first and second argument of the cross-correlations. If the spectra are invariant under time reversal, the underlying cross-correlation functions should also be invariant. Therefore, swapping the two arguments will yield identical correlation functions for an infinite trajectory length. However, in practice, sampling is not perfect, and the two correlation functions differ due to statistical noise. By simply taking the arithmetic average of both correlation functions (with forward time and reversed time) [see Eq. (19)], the sampling is effectively increased by a factor of 2. We call this approach the “commutator trick,”

CCFa,b12at+τbtt+bt+τatt.
(19)

It only improves sampling if at and bt are different time series, i.e., in the case of true cross-correlations. If a and b are identical (i.e., autocorrelation, as with correlating the dipole moment for infrared spectra), nothing is gained.

Another subtlety arises if one of both properties a and b is the magnetic dipole moment [which is required for vibrational circular dichroism (VCD) and Raman optical activity (ROA) spectra; see below]. This quantity depends linearly on the atomic velocities. A time reversal inverts all atomic velocities and therefore swaps the sign of the magnetic dipole moment. In this case, taking the average of both correlation functions would only yield noise, as both contributions would cancel out. The sign of the correlation function with reversed time has to be changed to obtain the desired effect, which is given in the following equation:

CCF*a,b12at+τbttbt+τatt.
(20)

5. Correcting the Verlet frequency shift

Describing the dynamics of a system with finite time steps does unfortunately not leave the vibrational frequencies of harmonic oscillators invariant. If a harmonic oscillator with “true” frequency ω0 is solved by a Verlet integration scheme110 with integrator time step Δt, it can be derived105 that the observed frequency ω̃ is given by

ω̃=1Δtarccos112ω02Δt2.
(21)

Fortunately, the effect is rather small. When simulating an oscillator with ω0 = 3000.0 cm−1 and Δt = 0.5 fs, the observed frequency is ω̃=2990.0 cm−1, which is a deviation of around 0.3%. When computing spectra from a trajectory, the observed frequencies ω̃ are given and the true frequencies ω0 are sought after. The inverse of Eq. (21) can be easily written as

ω0=22cosω̃ΔtΔt.
(22)

By using Eq. (22), the frequency axis of the spectra can be corrected for the frequency shift of the Verlet integrator. The correction is implemented in TRAVIS.

When vibrational spectra shall be computed from molecular dynamics simulations, molecular electric properties (e.g., dipole moments) are required. The standard approach to obtain those is the so-called Wannier localization.111,112 A unitary transformation that minimizes the spatial spread is applied to the set of molecular orbitals, yielding localized Wannier orbitals. For each such orbital, the center of mass (called “Wannier center”) can be considered an integer point charge. By assigning the Wannier centers to individual molecules, the molecular dipole moment can be computed. In the case of one isolated molecule, this is not even an approximation. However, this approach comes with several limitations. First, it 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 that the iterative approach to determine the Wannier centers sometimes does not converge at all. Apart from that, there arise problems with aromatic systems. The Wannier localization enforces an alternating single bond–double bond pattern in aromatic systems, which can introduce artificial spectral bands, as we have recently reported.113 When computing molecular polarizabilities for Raman spectra of aromatic systems, the numerical noise introduced by the Wannier localization is several orders of magnitude larger than the spectrum itself. Finally, to the best of our knowledge, it is not possible to derive molecular quadrupole moments from Wannier centers; however, these quadrupole moments are required for the computation of ROA spectra.

To overcome these issues when computing molecular properties, we developed the Voronoi integration in 2015.113 After a periodic radical Voronoi tessellation67,70,71 of the system (with the atoms as Voronoi sites and atomic van der Waals radii114as Voronoi radii) has been computed by using the Voro++ library,68,69 the total electron density ρr, which is supplied on a grid, is directly integrated within the Voronoi cell of a molecule to yield the desired molecular moments (e.g., electric dipole momentpMoland electric quadrupole momentQijMol),

pMol=n=1NMolqnrnMolρrrd3r,
(23)
QijMol=n=1NMolqn3rn,irn,jrn2δijMolρr3rirjr2δijd3r,
(24)

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 have shown that the spectra obtained via Voronoi integration are very similar to those based on Wannier localization, but the problems with aromatic systems are no longer present.113 Furthermore, our integration algorithm is very fast (1 s per full frame on one CPU core) so that a lot of computational time is saved when compared to Wannier localization. Finally, also higher molecular multipole moments such as the electric quadrupole moment (which is required to compute ROA spectra) can easily be obtained. The required volumetric electron densities along the trajectory can be stored in our highly efficient lossless compression format BQB7 (see above). Therefore, we propose the Voronoi integration technique as a new standard method for obtaining molecular properties in bulk phase systems.

The Fourier transform of the atomic velocity autocorrelation function yields the power spectrum of a system, also known as the vibrational density of states (VDOS). In contrast to the experimentally accessible vibrational spectra, the power spectrum is not based on selection rules but contains all motions of a system. If the units are chosen correctly, and if atom mass weighting is applied, the integral over a peak in the power spectrum directly yields the temperature of the corresponding degree of freedom in Kelvins. This is an interesting feature to determine if a simulation is well equilibrated: According to the equipartition theorem, all degrees of freedom should possess the same temperature. If the integrals in the power spectrum indicate that this might not be the case, a longer equilibration can be recommended. One example is presented in Fig. 12, where a power spectrum of a system containing 20 different diatomic molecules with harmonic bonds (black curve) is shown together with the integrals over the bands (red curve). As the integrals are all of similar magnitude, the system seems to be well-equilibrated. Power spectra computed by TRAVIS already appeared many times in the literature.22,25,33,37,41,72,115–117

FIG. 12.

Power spectrum of a system containing 20 different harmonic diatomic molecules with different atom masses and force constants (black curve) together with the integral depicting the temperature of each normal mode (red curve).

FIG. 12.

Power spectrum of a system containing 20 different harmonic diatomic molecules with different atom masses and force constants (black curve) together with the integral depicting the temperature of each normal mode (red curve).

Close modal

Based on an approach presented by Mathias et al.,118,119 we implemented a bulk phase normal mode analysis in TRAVIS.115 By iteratively diagonalizing the matrix of all atomic velocity cross-correlation functions, the power spectrum of the system is decomposed into normal modes that are maximally localized in the frequency space. To do so, a set of reference structures is required, onto which the molecules in the trajectory are projected. Due to this requirement, the approach works best for rigid molecules and has severe limitations for flexible side chains longer than ethyl groups. In Fig. 13, an exemplary application of the normal mode analysis is shown, where the individual C–H stretching motions of [EMIm]+ cations are distinguished in bulk phase simulations of pure [EMIm][OAc] and an [EMIm][OAc]/water mixture.115 There also exist other articles in the literature where this feature of TRAVIS has been applied.37,120–122

FIG. 13.

With the normal mode decompositions of the power spectra in simulations of pure [EMIm][OAc] and an [EMIm][OAc]/water mixture, the C–H stretching bands of the [EMIm]+ cation can be distinguished.115 

FIG. 13.

With the normal mode decompositions of the power spectra in simulations of pure [EMIm][OAc] and an [EMIm][OAc]/water mixture, the C–H stretching bands of the [EMIm]+ cation can be distinguished.115 

Close modal

The infrared spectrum of a system can be computed as the Fourier transform of the molecular dipole autocorrelation function along a simulation trajectory. Soon after the first ab initio molecular dynamics (AIMD) simulations were performed, the first predicted infrared spectrum based on AIMD was published in 1997,123 and several applications and further developments followed.124–126 For Raman spectra, the molecular polarizability needs to be correlated instead. The first predicted Raman spectrum from AIMD appeared in 2002.127 It has been demonstrated many times in the literature that these spectra from AIMD simulations agree very well with the experiment. A module for the computation of such spectra was added to TRAVIS in 2013.128 While it initially made use of Wannier 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

The vibrational circular dichroism (VCD) spectrum of a molecular system can be computed via cross-correlating the molecular electric dipole moment pMolt and the molecular magnetic dipole moment mMolt. However, the magnetic dipole moment cannot be easily obtained from a Born–Oppenheimer molecular dynamics (BOMD) simulation, as magnetic moments arise from electric currents, and the Born–Oppenheimer approximation does not account for electric currents. One solution is to apply some kind of perturbation 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 in total electron density in subsequent steps along the trajectory.137 Once the electric current density jr has been obtained on a grid, integration over the molecular Voronoi cells yields the molecular magnetic dipole moments mMol via the classical expression

mMol=12n=1NMolqnrn×vn12Molr×jrd3r,
(25)

where vn is the velocity of the nth atom in the molecule. Based on this approach, which is implemented in TRAVIS, we published the first ab initio prediction of a bulk phase VCD spectrum in 2016.137 

Computing Raman optical activity (ROA) spectra from molecular 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 polarizability AMol, and the magnetic dipole–electric dipole polarizability GMol. It can be shown139 that

αMol=ddEpMol,
(26)
AMol=AMol=ddEQMol,
(27)
GMol=GMolT=ddEmMolT
(28)

with the electric dipole–electric quadrupole polarizability AMol and the electric dipole–magnetic dipole polarizability GMol so that all three quantities can be obtained via finite differences of an applied external electric field E. The schematic approach to compute the required polarizabilities for ROA spectra is presented in Fig. 14. We implemented this approach in TRAVIS and were able to present 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).

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.

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.

Close modal
FIG. 15.

Predicted Raman optical activity (ROA) spectrum138 of liquid (R)-propylene oxide (black) together with the experimental spectrum (red).140 This was the first predicted bulk phase ROA spectrum in the literature.

FIG. 15.

Predicted Raman optical activity (ROA) spectrum138 of liquid (R)-propylene oxide (black) together with the experimental spectrum (red).140 This was the first predicted bulk phase ROA spectrum in the literature.

Close modal

Resonance Raman spectra can be computed via temporal correlation of the dynamic polarizability tensor along a trajectory. One approach to obtain this tensor is based on real-time time-dependent density functional theory (RT-TDDFT).141 It has been applied several times in the literature,142,143 but only to static calculations of isolated molecules in vacuum. In 2019, we implemented this method in TRAVIS and applied it to BOMD trajectories in order to obtain the resonance Raman spectrum.144 In contrast to the existing methods, our implementation works for bulk phase systems and therefore is able to take into account full solvent influence as well as some anharmonic effects. Apart from that, our approach is not limited to generalized gradient approximation (GGA) DFT but works with all real time electron structure methods, including, e.g., DFT with hybrid functionals. In our article, we presented the first ab initio prediction of a bulk phase resonance Raman spectrum in the literature on the example of uracil in water. We have also shown that the solvent effect is very important in this case, as only the bulk phase spectrum with full solvent influence agrees well with the experiment.

Our newly developed method for predicting resonance Raman spectra does not require a certain laser wavelength as input but 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 picking laser wavelengths at which interesting resonance effects can be expected.

FIG. 16.

Predicted resonance Raman spectra of uracil in water for all possible laser wavelengths (vertical axis).144 Rows are normalized to show the relative intensity ratio. These were the first predicted bulk phase resonance Raman spectra in the literature.

FIG. 16.

Predicted resonance Raman spectra of uracil in water for all possible laser wavelengths (vertical axis).144 Rows are normalized to show the relative intensity ratio. These were the first predicted bulk phase resonance Raman spectra in the literature.

Close modal

In this article, we have presented some methods and analyses that have recently been introduced to the TRAVIS program package 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 first ab initio predictions in the literature of bulk vibrational circular dichroism (VCD) spectra,137 bulk phase Raman optical activity (ROA) spectra,138 and bulk phase resonance Raman spectra144 within the last years. Since the release in 2011, the original TRAVIS article1 has already been cited more than 400 times now, and more than 800 different co-authors contributed to these studies. As the authors of TRAVIS, we promise that we will continue to include new methods and features into the code while at the same time maintaining old and existing functions and improving the quality of the documentation.

We would like to thank all the TRAVIS users around the world for their trust (and sometimes, unfortunately, also patience). M.B. acknowledges financial support from the DFG through Project No. Br 5494/1-1.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
M.
Brehm
and
B.
Kirchner
, “
TRAVIS—A free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories
,”
J. Chem. Inf. Model.
51
(
8
),
2007
2023
(
2011
).
2.
W.
Humphrey
,
A.
Dalke
,
K.
Schulten
 et al, “
VMD: Visual molecular dynamics
,”
J. Mol. Graphics
14
,
33
38
(
1996
).
3.
N.
Michaud-Agrawal
,
E. J.
Denning
,
T. B.
Woolf
, and
O.
Beckstein
, “
MDAnalysis: A toolkit for the analysis of molecular dynamics simulations
,”
J. Comput. Chem.
32
,
2319
2327
(
2011
).
4.
I. T.
Todorov
,
W.
Smith
,
K.
Trachenko
, and
M. T.
Dove
, “
DL_POLY 3: New dimensions in molecular dynamics simulations via massive parallelism
,”
J. Mater. Chem.
16
,
1911
1918
(
2006
).
5.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
6.
D. A.
Case
,
T. E.
Cheatham
 III
,
T.
Darden
,
H.
Gohlke
,
R.
Luo
,
K. M.
Merz
, Jr.
,
A.
Onufriev
,
C.
Simmerling
,
B.
Wang
, and
R. J.
Woods
, “
The Amber biomolecular simulation programs
,”
J. Comput. Chem.
26
,
1668
1688
(
2005
).
7.
M.
Brehm
and
M.
Thomas
, “
An efficient lossless compression algorithm for trajectories of atom positions and volumetric data
,”
J. Chem. Inf. Model.
58
,
2092
2107
(
2018
).
8.
Grace Development Team
, Grace, (c), see http://plasma-gate.weizmann.ac.il/Grace, 1996–2008.
9.
T.
Williams
,
C.
Kelley
 et al, Gnuplot 4.6: An Interactive Plotting Program, http://gnuplot.sourceforge.net/ (
2013
).
10.
Persistence of Vision Pty. Ltd.
, Persistence of Vision Raytracer (Version 3.6), retrieved from http://www.povray.org/download/ (
2004
).
11.
W. R. Inc.
, Mathematica, Version 12.0, Champaign, IL, 2019.
12.
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
2838
(
2008
).
13.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
CRC Press
,
1981
).
14.
H. J. C.
Berendsen
,
D.
van der Spoel
, and
R.
van Drunen
, “
GROMACS: A message-passing parallel molecular dynamics implementation
,”
Comput. Phys. Commun.
91
,
43
56
(
1995
).
15.
C. A.
Shelley
and
M. E.
Munk
, “
Computer perception of topological symmetry
,”
J. Chem. Inf. Comput. Sci.
17
,
110
113
(
1977
).
16.
H. L.
Morgan
, “
The generation of a unique machine description for chemical structures: A technique developed at chemical abstracts service
,”
J. Chem. Doc.
5
,
107
113
(
1965
).
17.
C. A.
Shelley
, “
Heuristic approach for displaying chemical structures
,”
J. Chem. Inf. Comput. Sci.
23
,
61
65
(
1983
).
18.
M.
Brehm
and
M.
Thomas
, The BQB format, see https://brehm-research.de/bqb.
19.
K. A. F.
Röhrig
and
T. D.
Kühne
, “
Optimal calculation of the pair correlation function for an orthorhombic system
,”
Phys. Rev. E
87
,
045301
(
2013
).
20.
M.
Kohagen
,
M.
Brehm
,
J.
Thar
,
W.
Zhao
,
F.
Müller-Plathe
, and
B.
Kirchner
, “
Performance of quantum chemically derived charges and persistence of ion cages in ionic liquids. A molecular dynamics simulations study of 1-n-butyl-3-methylimidazolium bromide
,”
J. Phys. Chem. B
115
,
693
702
(
2011
).
21.
M.
Brüssel
,
M.
Brehm
,
T.
Voigt
, and
B.
Kirchner
, “
Ab initio molecular dynamics simulations of a binary system of ionic liquids
,”
Phys. Chem. Chem. Phys.
13
,
13617
13620
(
2011
).
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
).
23.
M.
Brehm
,
H.
Weber
,
A. S.
Pensado
,
A.
Stark
, and
B.
Kirchner
, “
Proton 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-3-methylimidazolium acetate-water mixtures—Part 1
,”
Phys. Chem. Chem. Phys.
14
,
5030
5044
(
2012
).
24.
M.
Brehm
,
H.
Weber
,
A. S.
Pensado
,
A.
Stark
, and
B.
Kirchner
, “
Liquid structure and cluster formation in ionic liquid/water mixtures—An extensive ab initio molecular dynamics study on 1-ethyl-3-methylimidazolium acetate/water mixtures—Part 2
,”
Z. Phys. Chem.
227
,
177
204
(
2013
).
25.
F.
Malberg
,
M.
Brehm
,
O.
Hollóczki
,
A. S.
Pensado
, and
B.
Kirchner
, “
Understanding the evaporation of ionic liquids using the example of 1-ethyl-3-methylimidazolium ethylsulfate
,”
Phys. Chem. Chem. Phys.
15
,
18424
18436
(
2013
).
26.
M.
Thomas
,
M.
Brehm
,
O.
Hollóczki
, and
B.
Kirchner
, “
How can a carbene be active in an ionic liquid?
,”
Chem. - Eur. J.
20
,
1622
1629
(
2014
).
27.
O.
Hollóczki
,
M.
Macchiagodena
,
H.
Weber
,
M.
Thomas
,
M.
Brehm
,
A.
Stark
,
O.
Russina
,
A.
Triolo
, and
B.
Kirchner
, “
Triphilic ionic-liquid mixtures: Fluorinated and non-fluorinated aprotic ionic-liquid mixtures
,”
ChemPhysChem
16
,
3325
3333
(
2015
).
28.
M.
Brehm
,
H.
Weber
,
M.
Thomas
,
O.
Hollóczki
, and
B.
Kirchner
, “
Domain analysis in nanostructured liquids: A post-molecular dynamics study at the example of ionic liquids
,”
ChemPhysChem
16
,
3271
3277
(
2015
).
29.
M.
Brehm
,
G.
Saddiq
,
T.
Watermann
, and
D.
Sebastiani
, “
Influence of small fluorophilic and lipophilic organic molecules on dipalmitoylphosphatidylcholine bilayers
,”
J. Phys. Chem. B
121
,
8311
8321
(
2017
).
30.
S.
Gehrke
,
M.
von Domaros
,
R.
Clark
,
O.
Hollóczki
,
M.
Brehm
,
T.
Welton
,
A.
Luzar
, and
B.
Kirchner
, “
Structure and lifetimes in ionic liquids and their mixtures
,”
Faraday Discuss.
206
,
219
245
(
2018
).
31.
S.
Pylaeva
,
M.
Brehm
, and
D.
Sebastiani
, “
Salt bridge in aqueous solution: Strong structural motifs but weak enthalpic effect
,”
Sci. Rep.
8
,
13626
(
2018
).
32.
M.
Brehm
,
M.
Pulst
,
J.
Kressler
, and
D.
Sebastiani
, “
Triazolium-based ionic liquids: A novel class of cellulose solvents
,”
J. Phys. Chem. B
123
,
3994
4003
(
2019
).
33.
M.
Brüssel
,
M.
Brehm
,
A. S.
Pensado
,
F.
Malberg
,
M.
Ramzan
,
A.
Stark
, and
B.
Kirchner
, “
On the ideality of binary mixtures of ionic liquids
,”
Phys. Chem. Chem. Phys.
14
,
13204
13215
(
2012
).
34.
A.
Korotkevich
,
D. S.
Firaha
,
A. A. H.
Padua
, and
B.
Kirchner
, “
Ab initio molecular dynamics simulations of SO2 solvation in choline chloride/glycerol deep eutectic solvent
,”
Fluid Phase Equilib.
448
,
59
68
(
2017
).
35.
S.
Zahn
, “
Deep eutectic solvents: Similia similibus solvuntur?
,”
Phys. Chem. Chem. Phys.
19
(
5
),
4041
4047
(
2017
).
36.
L.
Gontrani
,
F.
Trequattrini
,
O.
Palumbo
,
L.
Bencivenni
, and
A.
Paolone
, “
New experimental evidences regarding conformational equilibrium in ammonium−bis(trifluoromethanesulfonyl)imide ionic liquids
,”
ChemPhysChem
19
(
20
),
2776
2781
(
2018
).
37.
M.
Brehm
and
D.
Sebastiani
, “
Simulating structure and dynamics in small droplets of 1-ethyl-3-methylimidazolium acetate
,”
J. Chem. Phys.
148
,
193802
(
2018
).
38.
L. M.
Lawson Daku
, “
Spin-state dependence of the structural and vibrational properties of solvated iron(II) polypyridyl complexes from AIMD simulations: II. Aqueous [Fe(TPy)2]Cl2
,”
Phys. Chem. Chem. Phys.
21
(
2
),
650
661
(
2019
).
39.
J.
Thar
,
M.
Brehm
,
A. P.
Seitsonen
, and
B.
Kirchner
, “
Unexpected hydrogen bond dynamics in imidazolium-based ionic liquids
,”
J. Phys. Chem. B
113
,
15129
15132
(
2009
).
40.
M.
Kohagen
,
M.
Brehm
,
Y.
Lingscheid
,
R.
Giernoth
,
J.
Sangoro
,
F.
Kremer
,
S.
Naumov
,
C.
Iacob
,
J.
Kärger
,
R.
Valiullin
, and
B.
Kirchner
, “
How hydrogen bonds influence the mobility of imidazolium-based ionic liquids. A combined theoretical and experimental study of 1-N-butyl-3-methylimidazolium bromide
,”
J. Phys. Chem. B
115
,
15280
15288
(
2011
).
41.
K.
Wendler
,
M.
Brehm
,
F.
Malberg
,
B.
Kirchner
, and
L.
Delle Site
, “
Short time dynamics of ionic liquids in AIMD-based power spectra
,”
J. Chem. Theory Comput.
8
,
1570
1579
(
2012
).
42.
O.
Hollóczki
,
D. S.
Firaha
,
J.
Friedrich
,
M.
Brehm
,
R.
Cybik
,
M.
Wild
,
A.
Stark
, and
B.
Kirchner
, “
Carbene formation in ionic liquids: Spontaneous, induced, or prohibited?
,”
J. Phys. Chem. B
117
,
5898
5907
(
2013
).
43.
G.
Bekçioǧlu
,
C.
Allolio
, and
D.
Sebastiani
, “
Water wires in aqueous solutions from first-principles calculations
,”
J. Phys. Chem. B
119
,
4053
4060
(
2015
).
44.
M. L. S.
Batista
,
G.
Pérez-Sánchez
,
J. R. B.
Gomes
,
J. A. P.
Coutinho
, and
E. J.
Maginn
, “
Evaluation of the GROMOS 56ACARBO force field for the calculation of structural, volumetric, and dynamic properties of aqueous glucose systems
,”
J. Phys. Chem. B
119
,
15310
15319
(
2015
).
45.
Q. R.
Sheridan
,
S.
Oh
,
O.
Morales-Collazo
,
E. W.
Castner
,
J. F.
Brennecke
, and
E. J.
Maginn
, “
Liquid structure of CO2-reactive aprotic heterocyclic anion ionic liquids from X-ray scattering and molecular dynamics
,”
J. Phys. Chem. B
120
,
11951
11960
(
2016
).
46.
F.
Lo Celso
,
G. B.
Appetecchi
,
C. J.
Jafta
,
L.
Gontrani
,
J. N.
Canongia Lopes
,
A.
Triolo
, and
O.
Russina
, “
Nanoscale organization in the fluorinated room temperature ionic liquid: Tetraethyl ammonium (trifluoromethanesulfonyl)(nonafluorobutylsulfonyl)imide
,”
J. Chem. Phys.
148
,
193816
(
2018
).
47.
K.
Bernardino
,
T. A.
Lima
, and
M. C. C.
Ribeiro
, “
Low-temperature phase transitions of the ionic liquid 1-ethyl-3-methylimidazolium dicyanamide
,”
J. Phys. Chem. B
123
(
44
),
9418
9427
(
2019
).
48.
M.
Zhao
,
B.
Wu
,
S. I.
Lall-Ramnarine
,
J. D.
Ramdihal
,
K. A.
Papacostas
,
E. D.
Fernandez
,
R. A.
Sumner
,
C. J.
Margulis
,
J. F.
Wishart
, and
E. W.
Castner
, “
Structural analysis of ionic liquids with symmetric and asymmetric fluorinated anions
,”
J. Chem. Phys.
151
,
074504
(
2019
).
49.
X.-Y.
Guo
,
C.
Peschel
,
T.
Watermann
,
G. F. v.
Rudorff
, and
D.
Sebastiani
, “
Cluster formation of polyphilic molecules solvated in a DPPC bilayer
,”
Polymers
9
,
488
(
2017
).
50.
M. H.
Dokoohaki
,
A. R.
Zolghadr
, and
A.
Klein
, “
Impact of the chemical structure on the distribution of neuroprotective N-alkyl-9H-carbazoles at octanol/water interfaces
,”
New J. Chem.
44
(
4
),
1211
1220
(
2020
).
51.
L. K.
Scarbath-Evers
,
R.
Hammer
,
D.
Golze
,
M.
Brehm
,
D.
Sebastiani
, and
W.
Widdra
, “
From flat to tilted: Gradual interfaces in organic thin film growth
,”
Nanoscale
12
,
3834
3845
(
2020
).
52.
O.
Russina
,
F. L.
Celso
, and
A.
Triolo
, “
Pressure-responsive mesoscopic structures in room temperature ionic liquids
,”
Phys. Chem. Chem. Phys.
17
,
29496
29500
(
2015
).
53.
T. C.
Lourenço
,
Y.
Zhang
,
L. T.
Costa
, and
E. J.
Maginn
, “
A molecular dynamics study of lithium-containing aprotic heterocyclic ionic liquid electrolytes
,”
J. Chem. Phys.
148
,
193834
(
2018
).
54.
F.
Lo Celso
,
G. B.
Appetecchi
,
E.
Simonetti
,
M.
Zhao
,
E. W.
Castner
,
U.
Keiderling
,
L.
Gontrani
,
A.
Triolo
, and
O.
Russina
, “
Microscopic structural and dynamic features in triphilic room temperature ionic liquids
,”
Front. Chem.
7
,
285
(
2019
).
55.
M.
Macchiagodena
,
G.
Mancini
,
M.
Pagliai
,
G.
Cardini
, and
V.
Barone
, “
New atomistic model of pyrrole with improved liquid state properties and structure
,”
Int. J. Quantum Chem.
118
(
9
),
e25554
(
2018
).
56.
V. F.
Sears
, “
Neutron scattering lengths and cross sections
,”
Neutron News
3
,
26
37
(
1992
).
57.
E.
Lorch
, “
Neutron diffraction by germania, silica and radiation-damaged silica glasses
,”
J. Phys. C: Solid State Phys.
2
,
229
237
(
1969
).
58.
F. L.
Celso
,
B.
Aoun
,
A.
Triolo
, and
O.
Russina
, “
Liquid structure of dibutyl sulfoxide
,”
Phys. Chem. Chem. Phys.
18
(
23
),
15980
15987
(
2016
).
59.
M.
Campetella
,
M.
Montagna
,
L.
Gontrani
,
E.
Scarpellini
, and
E.
Bodo
, “
Unexpected proton mobility in the bulk phase of cholinium-based ionic liquids: New insights from theoretical calculations
,”
Phys. Chem. Chem. Phys.
19
,
11869
11880
(
2017
).
60.
U.
Kapoor
and
J. K.
Shah
, “
Globular, sponge-like to layer-like morphological transition in 1-n-alkyl-3-methylimidazolium octylsulfate ionic liquid homologous series
,”
J. Phys. Chem. B
122
,
213
228
(
2018
).
61.
M.
Campetella
,
A.
Mariani
,
C.
Sadun
,
B.
Wu
,
E. W.
Castner
, and
L.
Gontrani
, “
Structure and dynamics of propylammonium nitrate-acetonitrile mixtures: An intricate multi-scale system probed with experimental and theoretical techniques
,”
J. Chem. Phys.
148
(
13
),
134507
(
2018
).
62.
F. L.
Celso
,
Y.
Yoshida
,
R.
Lombardo
,
C. J.
Jafta
,
L.
Gontrani
,
A.
Triolo
, and
O.
Russina
, “
Mesoscopic structural organization in fluorinated room temperature ionic liquids
,”
C. R. Chim.
21
(
8
),
757
770
(
2018
).
63.
L. S.
Vermeer
,
B. L.
de Groot
,
V.
Réat
,
A.
Milon
, and
J.
Czaplicki
, “
Acyl chain order parameter profiles in phospholipid bilayers: Computation from molecular dynamics simulations and comparison with 2H NMR experiments
,”
Eur. Biophys. J.
36
,
919
931
(
2007
).
64.
R.
Giernoth
,
A.
Bröhl
,
M.
Brehm
, and
Y.
Lingscheid
, “
Interactions in ionic liquids probed by in situ NMR spectroscopy
,”
J. Mol. Liq.
192
,
55
58
(
2014
).
65.
A. B. W.
Kennedy
and
H. R.
Sankey
, “
The thermal efficiency of steam engines. Report of the committee appointed to the council upon the subject of the definition of a standard or standards of thermal efficiency for steam engines: With an introductory note. (Including appendixes and plate at back of volume)
,”
Minutes Proc. Inst. Civil Eng.
134
,
278
312
(
1898
).
66.
J. D.
Hunter
, “
Matplotlib: A 2D graphics environment
,”
Comput. Sci. Eng.
9
,
90
95
(
2007
).
67.
G.
Voronoi
, “
Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire. Recherches sur les parallélloèdres primitifs
,”
J. Reine Angew. Math.
1908
(
134
),
198
(
1908
).
68.
C. H.
Rycroft
, “
Voro++: A three-dimensional Voronoi cell library in C++
,”
Chaos
19
,
041111
(
2009
).
69.
C. H.
Rycroft
,
G. S.
Grest
,
J. W.
Landry
, and
M. Z.
Bazant
, “
Analysis of granular flow in a pebble-bed nuclear reactor
,”
Phys. Rev. E
74
,
021306
(
2006
).
70.
B. J.
Gellatly
and
J. L.
Finney
, “
Calculation of protein volumes: An alternative to the Voronoi procedure
,”
J. Mol. Biol.
161
,
305
322
(
1982
).
71.
F. M.
Richards
, “
The interpretation of protein structures: Total volume, group volume distributions and packing density
,”
J. Mol. Biol.
82
,
1
14
(
1974
).
72.
D. S.
Firaha
,
M.
Kavalchuk
, and
B.
Kirchner
, “
SO2 solvation in the 1-ethyl-3-methylimidazolium thiocyanate ionic liquid by incorporation into the extended cation–anion network
,”
J. Solution Chem.
44
,
838
849
(
2015
).
73.
T. D. N.
Reddy
and
B. S.
Mallik
, “
Heterogeneity in the microstructure and dynamics of tetraalkylammonium hydroxide ionic liquids: Insight from classical molecular dynamics simulations and Voronoi tessellation analysis
,”
Phys. Chem. Chem. Phys.
22
(
6
),
3466
3480
(
2020
).
74.
A.
Mariani
,
R.
Caminiti
,
M.
Campetella
, and
L.
Gontrani
, “
Pressure-induced mesoscopic disorder in protic ionic liquids: First computational study
,”
Phys. Chem. Chem. Phys.
18
(
4
),
2297
2302
(
2016
).
75.
D. O.
Abranches
,
N.
Schaeffer
,
L. P.
Silva
,
M. A. R.
Martins
,
S. P.
Pinho
, and
J. A. P.
Coutinho
, “
The role of charge transfer in the formation of type I deep eutectic solvent-analogous ionic liquid mixtures
,”
Molecules
24
(
20
),
3687
(
2019
).
76.
T.
Cosby
,
U.
Kapoor
,
J. K.
Shah
, and
J.
Sangoro
, “
Mesoscale organization and dynamics in binary ionic liquid mixtures
,”
J. Phys. Chem. Lett.
10
(
20
),
6274
6280
(
2019
).
77.
R.
Macchieraldo
,
S.
Gehrke
,
N. K.
Batchu
,
B.
Kirchner
, and
K.
Binnemans
, “
Tuning solvent miscibility: A fundamental assessment at the example of induced methanol/n-dodecane phase separation
,”
J. Phys. Chem. B
123
,
4400
4407
(
2019
).
78.
T. D. N.
Reddy
and
B. S.
Mallik
, “
Nanostructure domains, voids, and low-frequency spectra in binary mixtures of N,N-dimethylacetamide and ionic liquids with varying cationic size
,”
RSC Adv.
10
(
3
),
1811
1827
(
2020
).
79.
S.
Gehrke
,
R.
Macchieraldo
, and
B.
Kirchner
, “
Understanding the fluidity of condensed phase systems in terms of voids—Novel algorithm, implementation and application
,”
Phys. Chem. Chem. Phys.
21
,
4988
4997
(
2019
).
80.
X.
Huang
,
C. J.
Margulis
,
Y.
Li
, and
B. J.
Berne
, “
Why is the partial molar volume of CO2 so small when dissolved in a room temperature ionic liquid? Structure and dynamics of CO2 dissolved in [BMIm+][PF6]
,”
J. Am. Chem. Soc.
127
,
17842
17851
(
2005
).
81.
N. N.
Medvedev
,
V. P.
Voloshin
,
V. A.
Luchnikov
, and
M. L.
Gavrilova
, “
An algorithm for three-dimensional Voronoi S-network
,”
J. Comput. Chem.
27
,
1676
1692
(
2006
).
82.
J.
Ingenmey
,
S.
Gehrke
, and
B.
Kirchner
, “
How to harvest Grotthuss diffusion in protic ionic liquid electrolyte systems
,”
ChemSusChem
11
,
1900
1910
(
2018
).
83.
A.
Einstein
, “
Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen
,”
Ann. Phys.
322
,
549
560
(
1905
).
84.
V.
Calandrini
,
E.
Pellegrini
,
P.
Calligari
,
K.
Hinsen
, and
G. R.
Kneller
, “
nMoldyn—Interfacing spectroscopic experiments, molecular dynamics simulations and models for time correlation functions
,”
École Thématique Soc. Française Neutronique
12
,
201
232
(
2011
).
85.
R.
Stefanovic
,
G. B.
Webber
, and
A. J.
Page
, “
Nanostructure of propylammonium nitrate in the presence of poly(ethylene oxide) and halide salts
,”
J. Chem. Phys.
148
(
19
),
193826
(
2018
).
86.
U.
Cerajewski
,
J.
Träger
,
S.
Henkel
,
A. H.
Roos
,
M.
Brehm
, and
D.
Hinderberger
, “
Nanoscopic structures and molecular interactions leading to a dystectic and two eutectic points in [EMIm][Cl]/urea mixtures
,”
Phys. Chem. Chem. Phys.
20
,
29591
29600
(
2018
).
87.
H.
Scheiber
,
Y.
Shi
, and
R. Z.
Khaliullin
, “
Communication: Compact orbitals enable low-cost linear-scaling ab initio molecular dynamics for weakly-interacting systems
,”
J. Chem. Phys.
148
(
23
),
231103
(
2018
).
88.
M. S.
Green
, “
Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids
,”
J. Chem. Phys.
22
,
398
413
(
1954
).
89.
J.
Wuttke
,
LMFIT—A C Library for Levenberg-Marquardt Least-Squares Minimization and Curve Fitting
, based on work by
B. S.
Garbow
,
K. E.
Hillstrom
,
J. J.
Moré
, and
S.
Moshier
, retrieved in 2010 from https://jugit.fz-juelich.de/mlz/lmfit.
90.
M.
Macchiagodena
,
G. D.
Frate
,
G.
Brancato
,
B.
Chandramouli
,
G.
Mancini
, and
V.
Barone
, “
Computational study of the DPAP molecular rotor in various environments: From force field development to molecular dynamics simulations and spectroscopic calculations
,”
Phys. Chem. Chem. Phys.
19
(
45
),
30590
30602
(
2017
).
91.
M. H.
Kowsari
and
S.
Ebrahimi
, “
Capturing the effect of [PF3(C2F5)3] vs. [PF6], flexible anion vs. rigid, and scaled charge vs. unit on the transport properties of [BMIm]+-based ionic liquids: A comparative MD study
,”
Phys. Chem. Chem. Phys.
20
(
19
),
13379
13393
(
2018
).
92.
C.
Dreßler
and
D.
Sebastiani
, “
Effect of anion reorientation on proton mobility in the solid acids family CsHyXO4 (X = S, P, Se, y = 1, 2) from ab initio molecular dynamics simulations
,”
Phys. Chem. Chem. Phys.
(published online 2019).
93.
F.
Stillinger
, “
Theoretical approaches to the intermolecular nature of water
,”
Philos. Trans. R. Soc., B
278
,
97
112
(
1977
).
94.
D. C.
Rapaport
, “
Hydrogen bonds in water: Network organization and lifetimes
,”
Mol. Phys.
50
,
1151
1162
(
1983
).
95.
A.
Khintchine
, “
Korrelationstheorie der stationären stochastischen prozesse
,”
Math. Ann.
109
(
1
),
604
615
(
1934
).
96.
A.
Luzar
and
D.
Chandler
, “
Effect of environment on hydrogen bond dynamics in liquid water
,”
Phys. Rev. Lett.
76
,
928
(
1996
).
97.
A.
Luzar
and
D.
Chandler
, “
Hydrogen-bond kinetics in liquid water
,”
Nature
379
,
55
(
1996
).
98.
A.
Luzar
, “
Resolving the hydrogen bond dynamics conundrum
,”
J. Chem. Phys.
113
,
10663
10675
(
2000
).
99.
S.
Gehrke
and
B.
Kirchner
, “
Robustness of the hydrogen bond and ion pair dynamics in ionic liquids to different parameters from the reactive flux method
,”
J. Chem. Eng. Data
65
,
1146
1158
(
2019
).
100.
S.
Ghahramani
,
F.
Yousefi
,
S. M.
Hosseini
, and
S.
Aparicio
, “
High-pressure behavior of 2-hydroxyethylammonium acetate ionic liquid: Experiment and molecular dynamics
,”
J. Supercrit. Fluids
155
,
104664
(
2020
).
101.
S.
Gehrke
and
O.
Hollóczki
, “
Hydrogen bonding of N-heterocyclic carbenes in solution: Mechanisms of solvent reorganization
,”
Chem. - Eur. J.
24
,
11594
11604
(
2018
).
102.
T.
Stettner
,
S.
Gehrke
,
P.
Ray
,
B.
Kirchner
, and
A.
Balducci
, “
Water in protic ionic liquids: Properties and use of a novel class of electrolytes for energy storage devices
,”
ChemSusChem
12
,
3827
3836
(
2019
).
103.
L.
Van Hove
, “
Correlations in space and time and Born approximation scattering in systems of interacting particles
,”
Phys. Rev.
95
,
249
262
(
1954
).
104.
R. G.
Gordon
, “
Molecular motion in infrared and Raman spectra
,”
J. Chem. Phys.
43
,
1307
1312
(
1965
).
105.
M.
Thomas
, “
Theoretical modeling of vibrational spectra in the liquid phase
,” Springer theses,
Springer International Publishing
,
2016
.
106.
J.
Hutter
,
M.
Iannuzzi
,
F.
Schiffmann
, and
J.
VandeVondele
, “
CP2K: Atomistic simulations of condensed matter systems
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
15
25
(
2014
).
107.
M.
Brehm
and
M.
Thomas
, Tutorial on Computing Vibrational Spectra, see https://brehm-research.de/spectroscopy.
108.
M.
Borgerding
, KISS FFT Library, see http://sourceforge.net/projects/kissfft/.
109.
P.
Kahlig
, “
Some aspects of Julius Von Hann’s contribution to modern climatology
,” in
Interactions Between Global Climate Subsystems
(
American Geophysical Union
,
2013
), pp.
1
7
.
110.
L.
Verlet
, “
Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard–Jones molecules
,”
Phys. Rev.
159
,
98
103
(
1967
).
111.
G. H.
Wannier
, “
The structure of electronic excitation levels in insulating crystals
,”
Phys. Rev.
52
,
191
197
(
1937
).
112.
N.
Marzari
and
D.
Vanderbilt
, “
Maximally localized generalized Wannier functions for composite energy bands
,”
Phys. Rev. B
56
,
12847
12865
(
1997
).
113.
M.
Thomas
,
M.
Brehm
, and
B.
Kirchner
, “
Voronoi dipole moments for the simulation of bulk phase vibrational spectra
,”
Phys. Chem. Chem. Phys.
17
,
3207
3213
(
2015
).
114.
M.
Mantina
,
A. C.
Chamberlin
,
R.
Valero
,
C. J.
Cramer
, and
D. G.
Truhlar
, “
Consistent van der Waals radii for the whole main group
,”
J. Phys. Chem. A
113
,
5806
5812
(
2009
).
115.
M.
Thomas
,
M.
Brehm
,
O.
Hollóczki
,
Z.
Kelemen
,
L.
Nyulászi
,
T.
Pasinszki
, and
B.
Kirchner
, “
Simulating the vibrational spectra of ionic liquid systems: 1-ethyl-3-methylimidazolium acetate and its mixtures
,”
J. Chem. Phys.
141
,
024510
(
2014
).
116.
M. C.
Kirkegaard
,
J. L.
Niedziela
,
A.
Miskowiec
,
A. E.
Shields
, and
B. B.
Anderson
, “
Elucidation of the structure and vibrational spectroscopy of synthetic metaschoepite and its dehydration product
,”
Inorg. Chem.
58
(
11
),
7310
7323
(
2019
).
117.
T.
Bauer
,
S.
Maisel
,
D.
Blaumeiser
,
J.
Vecchietti
,
N.
Taccardi
,
P.
Wasserscheid
,
A.
Bonivardi
,
A.
Görling
, and
J.
Libuda
, “
Operando drifts and DFT study of propane dehydrogenation over solid- and liquid-supported GaxPty catalysts
,”
ACS Catal.
9
(
4
),
2842
2853
(
2019
).
118.
G.
Mathias
and
M. D.
Baer
, “
Generalized normal coordinates for the vibrational analysis of molecular dynamics simulations
,”
J. Chem. Theory Comput.
7
,
2028
2039
(
2011
).
119.
G.
Mathias
,
S. D.
Ivanov
,
A.
Witt
,
M. D.
Baer
, and
D.
Marx
, “
Infrared spectroscopy of fluxional molecules from (ab initio) molecular dynamics: Resolving large-amplitude motion, multiple conformations, and permutational symmetries
,”
J. Chem. Theory Comput.
8
,
224
234
(
2012
).
120.
B.
Koeppe
,
S. A.
Pylaeva
,
C.
Allolio
,
D.
Sebastiani
,
E. T. J.
Nibbering
,
G. S.
Denisov
,
H.-H.
Limbach
, and
P. M.
Tolstoy
, “
Polar solvent fluctuations drive proton transfer in hydrogen bonded complexes of carboxylic acid with pyridines: NMR, IR and ab initio MD study
,”
Phys. Chem. Chem. Phys.
19
(
2
),
1010
1028
(
2017
).
121.
M.
Gług
,
M. Z.
Brela
,
M.
Boczar
,
A. M.
Turek
,
Ł.
Boda
,
M. J.
Wójcik
,
T.
Nakajima
, and
Y.
Ozaki
, “
Infrared spectroscopy and Born–Oppenheimer molecular dynamics simulation study on deuterium substitution in the crystalline benzoic acid
,”
J. Phys. Chem. B
121
(
3
),
479
489
(
2017
).
122.
L. M.
Lawson Daku
, “
Spin-state dependence of the structural and vibrational properties of solvated iron(II) polypyridyl complexes from AIMD simulations: Aqueous [Fe(bpy)3]Cl2, a case study
,”
Phys. Chem. Chem. Phys.
20
,
6236
6253
(
2018
).
123.
P. L.
Silvestrelli
,
M.
Bernasconi
, and
M.
Parrinello
, “
Ab initio infrared spectrum of liquid water
,”
Chem. Phys. Lett.
277
,
478
482
(
1997
).
124.
M. P.
Gaigeot
,
R.
Vuilleumier
,
M.
Sprik
, and
D.
Borgis
, “
Infrared spectroscopy of N-methylacetamide revisited by ab initio molecular dynamics simulations
,”
J. Chem. Theory Comput.
1
,
772
789
(
2005
).
125.
M.-P.
Gaigeot
,
M.
Martinez
, and
R.
Vuilleumier
, “
Infrared spectroscopy in the gas and liquid phase from first principle molecular dynamics simulations: Application to small peptides
,”
Mol. Phys.
105
,
2857
2878
(
2007
).
126.
M.-P.
Gaigeot
,
N. A.
Besley
, and
J. D.
Hirst
, “
Modeling the infrared and circular dichroism spectroscopy of a bridged cyclic diamide
,”
J. Phys. Chem. B
115
,
5526
5535
(
2011
).
127.
A.
Putrino
and
M.
Parrinello
, “
Anharmonic Raman spectra in high-pressure ice from ab initio simulations
,”
Phys. Rev. Lett.
88
,
176401
(
2002
).
128.
M.
Thomas
,
M.
Brehm
,
R.
Fligg
,
P.
Vöhringer
, and
B.
Kirchner
, “
Computing vibrational spectra from ab initio molecular dynamics
,”
Phys. Chem. Chem. Phys.
15
,
6608
6622
(
2013
).
129.
M.
Etinski
and
B.
Ensing
, “
Puzzle of the intramolecular hydrogen bond of dibenzoylmethane resolved by molecular dynamics simulations
,”
J. Phys. Chem. A
122
(
28
),
5945
5954
(
2018
).
130.
E. O.
Fetisov
,
D. B.
Harwood
,
I.-F. W.
Kuo
,
S. E. E.
Warrag
,
M. C.
Kroon
,
C. J.
Peters
, and
J. I.
Siepmann
, “
First-principles molecular dynamics study of a deep eutectic solvent: Choline chloride/urea and its mixture with water
,”
J. Phys. Chem. B
122
(
3
),
1245
1254
(
2018
).
131.
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 of experimental and theoretical methods
,”
Faraday Discuss.
211
,
425
439
(
2018
).
132.
G.
Cassone
,
J.
Sponer
,
S.
Trusso
, and
F.
Saija
, “
Ab initio spectroscopy of water under electric fields
,”
Phys. Chem. Chem. Phys.
21
(
38
),
21205
21212
(
2019
).
133.
B.
Dutta
and
J.
Chowdhury
, “
Existence of dimeric hydroxylamine-O-sulfonic acid: Experimental observations aided by ab initio, DFT, Car–Parrinello and Born–Oppenheimer on the fly dynamics
,”
Chem. Phys. Lett.
732
,
136645
(
2019
).
134.
H.
Hoshina
,
T.
Kanemura
, and
M. T.
Ruggiero
, “
Exploring the dynamics of bound water in nylon polymers with terahertz spectroscopy
,”
J. Phys. Chem. B
124
(
2
),
422
429
(
2020
).
135.
A.
Scherrer
,
R.
Vuilleumier
, and
D.
Sebastiani
, “
Nuclear velocity perturbation theory of vibrational circular dichroism
,”
J. Chem. Theory Comput.
9
,
5305
5312
(
2013
).
136.
A.
Scherrer
,
R.
Vuilleumier
, and
D.
Sebastiani
, “
Vibrational circular dichroism from ab initio molecular dynamics and nuclear velocity perturbation theory in the liquid phase
,”
J. Chem. Phys.
145
,
084101
(
2016
).
137.
M.
Thomas
and
B.
Kirchner
, “
Classical magnetic dipole moments for the simulation of vibrational circular dichroism by ab initio molecular dynamics
,”
J. Phys. Chem. Lett.
7
,
509
513
(
2016
).
138.
M.
Brehm
and
M.
Thomas
, “
Computing bulk phase Raman optical activity spectra from ab initio molecular dynamics simulations
,”
J. Phys. Chem. Lett.
8
,
3409
3414
(
2017
).
139.
D. A.
Long
,
The Raman Effect: A Unified Treatment of the Theory of Raman Scattering by Molecules
(
John Wiley & Sons Ltd.
,
2002
).
140.
J.
Šebestík
and
P.
Bouř
, “
Raman optical activity of methyloxirane gas and liquid
,”
J. Phys. Chem. Lett.
2
,
498
502
(
2011
).
141.
H.
Chen
,
J. M.
Mcmahon
,
M. A.
Ratner
, and
G. C.
Schatz
, “
Classical electrodynamics coupled to quantum mechanics for calculation of molecular optical properties: A RT-TDDFT/FDTD approach
,”
J. Phys. Chem. C
114
,
14384
14392
(
2010
).
142.
M.
Thomas
,
F.
Latorre
, and
P.
Marquetand
, “
Resonance Raman spectra of ortho-nitrophenol calculated by real-time time-dependent density functional theory
,”
J. Chem. Phys.
138
,
044101
(
2013
).
143.
J.
Mattiat
and
S.
Luber
, “
Efficient calculation of (resonance) Raman spectra and excitation profiles with real-time propagation
,”
J. Chem. Phys.
149
,
174108
(
2018
).
144.
M.
Brehm
and
M.
Thomas
, “
Computing bulk phase resonance Raman spectra from ab initio molecular dynamics and real-time TDDFT
,”
J. Chem. Theory Comput.
15
,
3901
3905
(
2019
).