Computing inelastic neutron scattering spectra from molecular dynamics trajectories

Inelastic neutron scattering (INS) provides a weighted density of phonon modes. Currently, INS spectra can only be interpreted for perfectly crystalline materials because of high computational cost for electronic simulations. INS has the potential to provide detailed morphological information if sufficiently large volumes and appropriate structural variety are simulated. Here, we propose a method that allows direct comparison between INS data with molecular dynamics simulations, a simulation method that is frequently used to simulate semicrystalline/amorphous materials. We illustrate the technique by analyzing spectra of a well-studied conjugated polymer, poly(3-hexylthiophene-2,5-diyl) (P3HT) and conclude that our technique provides improved volume and structural variety, but that the classical force field requires improvement before the morphology can be accurately interpreted.


Normal mode decomposition of Molecular Dynamics Trajectory
The aim of this section is to explore the physical meaning behind eqs. 3-6 in the main text. We start by postulating the existence of a set of j normal modes in the system, which generates an exact dynamical description of the atomic displacements according to, where u i (t) is the displacement of atom i at time t, u j (t) is the displacement along the normal mode coordinate j, and c i j is the unitary transformation matrix that connects the atomic and normal mode descriptions. In the harmonic approximation, the transformation matrix is time-independent, so taking the time derivative of both sides reproduces eq. 3 in the main text. With this in mind, we take the discrete Fourier transform of v i (t), We identify ω j as 2π j/t total and t n as n∆t, where t total is the length of the MD simulation, ∆t is the timestep of the simulation, j is an integer ranging from 0 to N − 1, and n is a similar integer. This identification puts the discrete Fourier transform in a more usual form: where N is the number of time-steps in the simulation. Since v i (t) must be a real number, the Fourier transform must be Hermitian, which means V i j = V * i,N− j . Thus, we can write the transform as, where the right-most exponential is always 1. If we write V i j as |V i j |e iδ j , then we can simplify the transform to, where we identify V i0 as the drift in the simulation box, which we set to zero. If we insert definitions of ω j and t n into the cosine argument, we get, Comparing to eq. 3 and eq. 4 in the main text, we see that, which means that our mode expansion in eq. 3 of the main text is analogous to the Fourier transform of the time series of the velocities of atom i. To be clear, this is an approximation in the general case, which becomes exactly true when all vibrational eigen-modes have non-degenerate frequencies. In this case, the Fourier transform will be able to pick up the dynamics of each mode independently. If there are modes that have degenerate eigen-frequencies, then constructive and destructive interference is possible between the degenerate Fourier components of the dynamics. As mentioned in the main text, we assume that each of the degenerate modes are decorrelated, which means the expectation value of the phase ( δ j ) is zero. We can see the impact of non-zero phases by writing out a set of degenerate modes at frequency ω j , where we observe that the degenerate modes can interfere with each other, and the interference is determined by the relative values of δ d . In a perfectly harmonic system, this is a problem, because the Fourier coefficients lose the information of the relative phases of the degenerate modes. However, we avoid this problem by noting that the system is weakly anharmonic, which decorrelates the modes, and we also use a weak thermostat to decorrelate the modes manually. In this case, the phase of each degenerate mode has a finite correlation time (e.g. δ d (0)δ d (t) ∼ e −t/τ ). If the total MD simulation time is much longer than the largest correlation time, then the Fourier transform becomes similar to an average over several Fourier transforms of randomized initial conditions (within the correct thermodynamic ensemble). The impact of this average is most easily observed when considering the Fourier transform of the velocity autocorrelation function. If we label the left-hand side of eq. S8 as v i (t), we find that its power spectrum is, where f j = ω j /(2π). Expanding the square leads to, We see that the d = d terms in the sum are always non-zero, and the d = d terms carry a phase e i(δ d −δ d ) that averages to zero if we average over a series of random phases as discussed above. Thus, we obtain an equation for degenerate modes in terms of the Fourier coefficients of v i (t) as, We remind the reader that V i j is the j−th Fourier component of the velocity function of atom i. As a result, we see that the relation that the square Fourier component |V i j | 2 is the proportional to the sum of the contributions of all modes that are degenerate at frequency ω j . To see how eq. 5 in the main text relates to the numerical objects used in the calculation, we insert the relationships between the mode expansion parameters and the Fourier transform into eq. 5, However, in a discrete basis, the delta function changes by δ ( f − f j ) → N∆tδ j j where δ j j is the Kronecker delta, and j becomes the free index of f in a discrete representation. Therefore, we can write the above equation as, where we have changed the summed index on the left-hand side to j to explicitly write the free index as j. The above equation demonstrates how eq. 5 in the main text is related to Fourier transforms of the atomic velocities, which shows how the approach in the main text is implemented in the code. Note that the left-hand side of above equation differs from eq. 5 in the main text because of the relation: 2πδ We derive eq. 8 in the main text from eqs. 5-7 by first rearranging eq. 7 into, which we insert into eq. 6 to get, Multiplying both sides by δ ( f − f j ), and summing over the index j, we get, where we identify the components of the velocity power spectrum from eq. 5 in the main text on the right-hand side of the above equation. Substituting eq. 5 in the above equation produces eq. 8 from the main text. We also see that the above equation can be written in terms of Fourier coefficents through eq. S13, which is how these equations are implemented in the code.

Deriving Intermediate Scattering Function from Velocity Correlations
We start by defining the convolution operation between B and any arbitrary function as an operator, which allows the general description of any overtone: The convolution operator is written into the scattering law as a power series of convolution operations. That power series simplifies to the exponential of the convolution operator: where v i (ω) is the Fourier transform of the classical velocity trajectory of atom i. By taking the Fourier transform of eq S19 we obtain: Using the identity of thermal averages, we show that eq S21 is equivalent to,

Correcting Anharmonic Effects in INS Spectra Using Classical Dynamics
One characteristic feature of anharmonic quantum oscillators is that they have decreasing energy-level spacings with increasing quantum number as shown in Figure S1b. In classical trajectories of isolated anharmonic oscillators, the power spectrum contains overtones in addition to the peak at the fundamental frequency. This differs from harmonic quantum oscillators, which always have the same energy difference between adjacent energy levels ( Figure S1a), and the classical power spectrum only has one peak at the fundamental frequency. We took advantage of this property in our formulation of neutron scattering spectra with the usage of Fourier convolutions to create overtone spectra. These convolutions implicitly assume a harmonic formulation because they create a set of peaks that are equally spaced in energy. While isolated oscillators can have anharmonic potential energy surfaces, it is also possible that anharmonicities couple two different harmonic modes together as in Figure S1c. The coupling between these states cause the transfer of energy between them, which means, when using phonon terminology, that the phonon populations diffuse between the different states. Once a phonon is created in an initial state, then it has a specific lifetime before it diffuses to a different state. This behavior manifests as a broadening of peak widths in the classical power spectrum. The coupling between different modes also decreases the energy level spacings. The most general description of anharmonicities in the solid state include coupled anharmonic potential energy surfaces as seen in Figure S1d. Since we seek to account for anharmonic effects in the INS spectra, we seek to augment the convolutions to account for the changes in energy level spacing. We find a starting description of anharmonic oscillators in Atkins 1 , in which the energy levels are given by a power series of the quantum number of the harmonic oscillator levels. If we truncate it to the lowest anharmonic order, E n =hω n + 1 2 + x e n + 1 2 2 (S25) then we calculate the difference between energy levels to be, We can see how the energy level spacing decreases with increasing n. We ignore the condition that the energy level difference becomes negative at large values of n because the important spectral features occur at low n, and the anharmonicity constant is expected to be small, particularly at low temperatures. We can correct our spectra by performing additional Fourier convolutions that shift the energy level differences. We define a function k i (ω) that represents the deviation from the harmonic energy level differences of the ith level. Using this set of functions, we can compute the nth overtone spectrum as, which is defined to be consistent with eq S26 above. This is not formally exact, and the approximation breaks down when there are multiple modes at the same frequency ω that have different energy level shifts. Thus, k i (ω) represents the average energy shift for modes at frequency ω. We need to define k i (ω) in terms of quantities calculated from classical trajectories. We start by assuming that the quantum mechanical energy differences can be described by second order perturbation theory, is equal to zero for the lowest order anharmonic operator (Γ i, j,k (â i +â † i )(â j +â † j )(â k +â † k ), whereâ i andâ † i are the annihilation and creation operators for mode i). Computing k i (ω) is equivalent to finding the set of E (2) i for all modes in the system.
The anharmonicities manifest as decay rates of periodic motions, which broadens peaks in the power spectrum (typically with a Lorentzian lineshape). If we model the classical dynamics as a set of decaying harmonic oscillators (as in the main text), then we can compute a set of decay constants (k j ) that correspond to frequencies (ω j ).
We aim to relate the relaxation rates (k j ) of the classical simulation to the second order energy correction (E (2) j ) found by perturbation theory to correct the overtones. We start with a first order correction to the time dependent wavefunction projected onto an arbitrary final state (assuming the initial state is the vibrational ground state for all oscillators).
where c f is the wavefunction coefficient for state f . In this equation, V is the anharmonic operator truncated to lowest order. From this, we can find the time-dependent population density of state f .
where V f i is the off-diagonal matrix element connecting states i and f . From this equation, we find that the time scale of the conversion from state i to f is π/ω. Now, we calculate the population transfer rate over this time scale. Defining the rate as, We can compute the time average from 0 to π/ω for each rate into all the possible final states. Since the rates can be combined additively to find the total rate leaving state i, we can define the total rate of disappearance from state i as, Carrying out the algebra we find a natural relationship between the average rate out of state i and E (2) i : This equation provides the quantum-classical correspondence that we need to correct for anharmonicities in our simulated INS spectra. However, since the dynamics were computed classically, the mode energies may be below the zero-point energy, which leads to inaccuracies in the observed anharmonic transfer rates. In addition, the rates computed from classical dynamics must be consistent with the thermodynamic ensemble of the simulation (we assume a NVT ensemble), so the RHS of eq S35 is the thermal average over the set of E (2) i . It may be sufficient to simply ignore the thermodynamic average, and assume that the anharmonic energy perturbations are equal the rates computed from low-temperature classical MD simulations. If one desires to correct for the thermodynamic average, then we suggest to perform multiple MD simulations at different temperatures, which will approximate the functional dependence of the set of rates on temperature, {R i (T )}. Finally, interpolate each of these functions to find the rate when the temperature is equal to¯h ω k B , where k B is the Boltzmann constant.

Impact of Thermostats and Barostats
For simulations of inelastic neutron scattering experiments, the choice of whether to use thermostats and barostats is difficult. Thermostats unphysically alter velocities on characteristic timescales, which can impact the observed velocity correlation functions use in our calculation of INS spectra. Barostats can also impact the observed velocity correlation functions, because the volume fluctuations also cause changes in interatomic distances, which impact potential energies and consequently the atomic dynamics. For both barostats and thermostats, we chose large time constants (τ = 10 ps), to minimize the impact on computed spectra; a 10 ps timescale corresponds to a frequency of ∼ 3 cm −1 , which is small enough such that it minimally impacts our computed spectra. We start by testing the impact of different thermostats (each with τ = 10 ps) in Figure S2a. In the main text, we chose to incorporate an Andersen-massive thermostat because our formalism required that all degenerate oscillators (degenerate meaning oscillating at the same frequency) be completely decorrelated. We chose to test a velocity-rescale, Andersen-massive, and no thermostat (NVE ensemble) in Figure S2a, and we see very little impact from both the presence of the thermostat, and type of thermostat used. The reason for choosing the two different thermostats is to show that complete decorrelation of all oscillators is not necessary at this scale; an Andersen-massive thermostat ensures complete oscillator decorrelation on the time scale τ, whereas the velocity-rescale thermostat will cause a decay of oscillator correlations on the time scale τ. If the phase of the independent oscillators in the simulation were heavily correlated, then there would be differences between peak heights computed from different thermostats (or no thermostat). The simulation results suggest that the oscillator phases are randomized by the large number of degenerate oscillators, or by anharmonicities present in the form of the classical forcefield. Each oscillator randomly starts at a phase between −π to π, which averages to 0 if there are many oscillators at the same frequency (law of large numbers). Since the simulation contains sixteen 20mer polymer chains, there are 320 identical chemical 6/8 groups, which is likely large enough to randomize the oscillator phases.
The impact of barostats may be less controllable than thermostats, but they may play a subtly important role in dynamical systems. In Figure S2b, we show the impact of using a Parrinello-Rahman barostat on the computed INS spectra. We see that the barostat dampens the low frequency modes in the system. To test if this effect depends on the simulation box size, we computed the INS spectrum for an NPT simulation of a 2 × 2 × 2 box of polymer chains (128 chains), and we see very little change between the computed spectra. This means that the lattice coordinates are coupled to the low energy modes of the molecules in the system, which is somewhat expected. However, the question of whether to use a barostat or not remains unanswered. It is important to note that a real crystalline material has phonon modes with wavelengths larger than the lattice dimensions used here, which are completely omitted from our simulations. This means that the barostat may mimic the action of these long range phonons, and thus may be physically appropriate to use for computing INS spectra. However, we did not choose to use a barostat in the main text because the mode dampening effect is strongly modulated by the compressibility value, which is unknown for P3HT. For our NpT simulations, we assumed that the compressibility of P3HT is 4.5 · 10 −5 bar −1 , which is the same as water at 300 K and 1 atm. Figure S3. Comparison between two amorphous simulated INS spectra; one with a velocity-rescale thermostat (1 picosecond time constant), and one without a thermostat. The inset expands the spectra at extremely low energies, revealing the position of a tall peak at ∼ 3 cm −1 .
We also compared the impact of the time constant of the thermostat on observed spectra. In Figure S2, the time constant was set to 10 picoseconds. In Figure S3, we demonstrate that a time constant of 1 picosecond using a velocity-rescale thermostat contributes a large, unphysical spike at ∼ 3 cm−1, which is eliminated when the thermostat is turned off. Interestingly, the frequency of an oscillator with a period of 1 picosecond is approximately 3 cm −1 . This is evidence that the thermostat actually has a measurable impact on the simulated INS spectrum. Going further, this implies that the thermostat contributes to the correlated dynamics in the system, which is somewhat unexpected because the thermostat acts stochastically on the system.

Experimental Inelastic Neutron Scattering Comparisons
In Figure S4, we compare the INS spectra for pristine P3HT, and P3HT doped with F4TCNQ at a 4% mole ratio. We see that the energies of the peaks do not change, only the heights of the peaks. This is important to show because the pristine measurement was performed during a different measurement cycle than the doped and regiorandom (RRa) measurements. Between the two measurement cycles, the VISION spectrometer was slightly modified, such that the two measurements are not directly comparable. The changes in peak heights in Figure S4a is the manifestation of the instrument changes. Ideally, we would like to directly compare RRa P3HT with pristine, undoped P3HT, but this is not a fair comparison. Since the peaks for 4% doped P3HT are at the same energies as pristine P3HT, we use it as a surrogate, while implicitly assuming that the changes in peak heights are completely due to the changes in instrumentation. As a result, we compare RRa P3HT with 4% doped P3HT to demonstrate the differences between the amorphous and semicrystalline phases of P3HT. In Figure S5 we show the comparison between the calculated INS spectrum derived from MD and DFT simulations over the full energy range. In the region between 1600 cm −1 −3000 cm −1 , there are no observed peaks, but the background signal comes from the overtones of the spectra. We show that the INS spectrum computed from MD displays the proper scatter signal strength in this region, which demonstrates that our methodology is accurately treating the overtone contributions to INS spectra.