Crystal effects in the vibrational spectra of one-dimensional molecular spin crossover crystals using molecular dynamics simulations

We use molecular dynamics simulations to model the iron phonon density of states (pDOS) of a spin crossover (SCO) crystal built from 39-nuclear [Fe(atrz)3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}Cl2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}] chains (atrz = 4-amino-1,2,4-triazole). This is possible by using an atomistic potential that depends on the spin state of the Fe atoms. While quantum-chemical methods are able to provide data on isolated SCO molecules, the computational costs necessitate the use of classical simulations to calculate crystal effects on the spectra. Data are provided for a high-spin (HS) crystal, in which all Fe atoms are in their HS state, and a low-spin (LS) crystal. HS and LS crystals show distinctively different spectra that are in agreement with experimental data. Compared to the single-molecule spectra, the crystal spectra show several distinctive features: (i) The spectra are smoother, and the main peaks are slightly reduced in amplitude. (ii) The effect of intermolecular interactions shows up in additional structure at long-wavelengths. Good agreement between experimental and simulated phonon spectra is observed. The lattice-dynamics analysis of the low-energy pDOS allows to extract the temperature dependence of the vibrational entropy of the LS and HS crystals. Our study thus provides novel information about the crystal effects in SCO crystals.


Introduction
The effects of lattice dynamics lie at the heart of the spin crossover (SCO) phenomenon [1][2][3][4][5]. Firstly, the metalligand vibrations (optical phonons) bring about the main vibrational entropy contribution of the spin transition. Secondly, the low-frequency vibrations (acoustic phonons) transmit the long-range elastic interactions determining the cooperativity of the spin transition in thermal equilibrium [6]. Thirdly, the lattice vibrations drive the fast structural trapping of the photoinduced high-spin state [7][8][9]. The effects of lattice-dynamical properties in the spin crossover material have been recently exhaustively reviewed by Molnár et al. [10].
In order to set up a computational model of these effects, a method is needed to provide insight into the vibrational properties of spin crossover solids. Yet, already the modeling of intermolecular effects-i.e., the calculation of electronic energies with or without geometry optimization-with quantum-chemical methods is a complex issue [11,12]. There are only few reports on the solid-state calculations of vibrations with density functional theory (DFT) methods [13][14][15][16][17][18]. This approach, surely being the future standard, has to accept a compromise between the computational demands and the choice of an appropriate functional and basis set. Another approach that yields a good reproduction of the vibrational patterns in the far-IR region, crucial for the SCO properties, is to model the 1D or 3D polynuclear complexes using their oligomeric models containing 7-15 Fe(II) centers [19][20][21]. In this case the use of B3LYP based functionals allows the modeling of SCO material in which the Fe(II) ions are bound via a system of covalent bonds rather than intermolecular interactions typical for molecular crystals. Interestingly, also the solid state DFT calculations seem to give better accuracy for the 3D polynuclear systems [17,18] than for the molecular crystals of the monomeric complexes [14][15][16].
In an attempt to model the larger systems, we have recently performed molecular dynamics (MD) simulations 1 3 345 Page 2 of 9 of the 1D polynuclear model of the spin crossover Fe(atrz) 3 Cl 2 complex (atrz = 4-amino-1,2,4-triazole) [22]. We determined force constants for bond stretching and bending of the high-spin (HS) and low-spin (LS) Fe atoms for use in atomistic molecular dynamics simulations by fitting to DFT data [20]. Our simulations allow to determine the phonon density of states (pDOS) of the iron atoms; this quantity describes the frequency distribution of the vibrational modes in the crystal in which Fe takes part and is experimentally accessible via the method of nuclear inelastic scattering (NIS) [23]. The Fe pDOS obtained within this approach for both high-spin and low-spin 25-nuclear 1D model of Fe(atrz) 3 Cl 2 fitted well to the experimental data obtained with NIS [20].
With this validated method of modeling of large chains we have decided to model crystals of such spin crossover 1D chains based on the [Fe(atrz) 3 Cl 2 ] complex [24]. MD offers the possibility of calculating the pDOS of models containing millions of atoms [25] and can handle crystals as easily as molecular structures [25,26]. For such models one may expect the ability to obtain the low-energy part of the pDOS corresponding to acoustic phonons from which the principal parameters of the lattice dynamicsentropy, specific heat and normalized force constants, among others-can be derived for both spin isomers [27]. In a very recent paper of the Toulouse group [28], a similar approach -however, with no DFT calibration of the force constantsfor a 3D SCO complex yielded a successful reproduction of the thermodynamic and mechanical parameters. Thus, the MD modeling would provide insight into lattice-dynamics effects. In addition, the MD modeling may be performed for crystals of different sizes which may shed light on the important problem of size effects observed for the 1D [29,30] and 3D [30] polymeric complexes.

Methods
The Fe(atrz) 3 Cl 2 system studied here was previously characterized by means of nuclear inelastic scattering and DFT calculations [20], yielding very good fits of the experimental and theoretical pDOS in the low-spin state. The force fields presented in the following are based on DFT modeling of the heptanuclear models of Fe(atrz) 3 Cl 2 [20]. Due to the large size of the systems studied here -containing 1275 atoms per molecule and up to 326 400 atoms for the largest crystal studied-a DFT calculation of the vibronic spectra is not possible and we resort to MD.
We built molecular crystals with a chain length of 39 iron cores, basing on the structure of [Fe(atrz) 3 (NO 3 ) 2 ⋅ 2 H 2 O] obtained in Ref. [31]. The molecule was constructed using the software Mercury 4.2.0 [32] from a Crystallographic Information File (CIF) obtained from the Cambridge Crystallographic Data Centre [33]. This data is based on XRD results by Grosjean et al. [31]. The chain length (14.25 nm) was chosen to be maximum for this software. Each molecule carries a charge of 78e, which is balanced out by 78 periodically arranged NO 3 ions, each with a charge of −1e . The packing mechanism of Mercury is able to reproduce this adjusted structure throughout the molecule and the subsequent periodic arrangement of the molecules in a molecular crystal. This structure is then saved in a Protein Data Bank (PDB) file. In the following process, we erased the water molecules from the PDB file and replaced NO 3 ions with Cl anions. The arrangement of the anions relative to the 1D chain was retained.
A representation of the structural model is shown in Fig. 1 with every third ligand and the chlorine ions omitted for clarity reasons. The Fe atoms terminating the chain on either end are only connected to three ligands each and are terminated with coordinated water. These atoms are aptly called terminal Fe atoms, Fe t , the inner atoms are designated Fe L or Fe H , depending on their spin state, i.e. low-spin or high-spin state.
In order to take the intramolecular interactions in the detailed geometry of the chain-like molecules into account, we employed the CHARMM potential [34]. The existing interaction parameters in the CHARMM database were extended by the parametrization by Meyer et al. [22] for bond-stretching and angle-bending interaction terms for Fe-ligand bonds in chain-like triazole SCO molecules. Here, special care was given to the fact that HS and LS Fe bonds differently to the N atoms. The CHARMM force field also includes non-bonded interaction terms in the form of the Coulomb and the van-der-Waals interaction, for which the CHARMM General Force Field (CGenFF) program [35] can supplement the relevant parameters. Partial charges were adopted from those calculated for nonanuclear molecules in Ref. [22]. SCO crystals can exist in a variety of spin states, since each internal Fe atom can be either in the HS or LS state; the terminal Fe atoms correspond to the HS state and were optimized in DFT as such. In the following, we focus on the two extreme spin configurations, in which all internal Fe atoms are high-spin (HS molecule) or all internal Fe atoms are low spin (LS molecule).
The molecules are arranged in a triclinic crystal structure with space group P 1 , determined by X-ray crystallography [31]. The molecules are aligned with their chains pointing in the same direction, denoted as [100], cf. The simulations were performed with the NAMD software [26,36] under periodic boundary conditions with a time step of 1 fs. After setting up the crystals according to the crystal structure, the system was relaxed for 50 ps at 100 K in an NVT environment. Temperature was controlled by a Langevin thermostat and fixed at T = 100 K.
The pDOS was determined analogous to the approaches used in [22,37]. We calculated the velocity autocorrelation function ( ) for all N Fe iron atoms in the molecular crystal where v n is the velocity vector of Fe atom n and ⟨… ⟩ t 0 denotes the average taken for different reference times t 0 . The autocorrelation function was calculated for times up to 50 ps. The Fourier transform of the velocity autocorrelation function then yields the pDOS g( ) , which is normalized to The calculation sampled all iron atoms present in the crystal and uses an average with equal weights. This method has proven useful for different simulation environments, such as single molecules or crystalline materials, as well as different temperatures. Figure 3 shows that the pDOS of the single 39-nuclear molecule has a rich band structure from 300 to 450 cm −1 in the LS state. This is due to the large iron movements within modes with strong Fe-N bending and stretching character [11,20]. In the HS state, the Fe-N distances are increased by about 10 % with respect to the LS state [1] which causes a downshift of the band structure to 200-280 cm −1 [38]. The shift in the peaks between the HS and LS spectrum is the main characteristics of these SCO molecules. Figure 4 shows how the spectrum changes if the molecules are arranged in a crystal structure of increasing size. The spectra become considerably smoothened. The main peaks (at 250 cm −1 for the HS, and at 370 cm −1 for the LS) are somewhat reduced in amplitude but not shifted in frequency. Figure 4 also demonstrates that our results are already well converged for the crystal sizes used here. In comparison to the 4 × 4 crystallite displayed in Fig. 2, which is elongated along the molecular axis, the 16 × 16 crystal has slightly larger extensions in the (100) plane perpendicular to the molecular axis, while the 12 × 12 crystal is nearly of a cubic cross section. Also this difference does not influence the

Results
spectra. In the following, we discuss the spectra obtained for the 16 × 16 crystal; a comparison of the lattice-dynamics parameters derived for the 16 × 16 and 12 × 12 models is given in the Supplementary Material. The most remarkable change between Figs. 3 and 4 appears at long wavelengths, resulting in an additional peak at a wavenumber of = 40 cm −1 . Fig. 5 focuses into the low-frequency part of the spectra. For crystals of a size bigger than 4 × 4, the pDOS can be approximated by at small wavenumbers, < 40 cm −1 , with a constant c. Figure 5 shows that the expected smooth Debye-like behavior of the pDOS for small wavenumbers is not sufficiently reproduced by the MD simulations of the 4 × 4 crystal. This finite-size effect disappears in the pDOS obtained for the 12 × 12 and 16 × 16 crystallites, which almost agree with each other. Figure 6 compares the converged crystal spectra of the HS and LS states. Several major differences with respect to the single-molecule spectra, Fig. 3 are noticeable. First of all, the frequency shift between LS and HS molecules is observable for the crystals as well: While the HS pDOS drops to zero at around 350 cm −1 , the LS crystal still produces signals at higher wavenumbers. Besides that, the major peaks in the pDOS did not change their position and shape: a flat single peak for the HS state and a narrow triple peak for the LS state. Both crystals feature the longwavelength phonons with an additional peak around 40 cm −1 .
The Fe partial density of states obtained with the MD approach presented in this study may be compared with those obtained experimentally by means of nuclear inelastic scattering (NIS) of synchrotron radiation [38]. The comparison of the MD data for the 16 × 16 crystal and experimentally obtained data for HS and LS is shown in Fig. 7.
For both the LS and HS phases, the MD crystal simulation correctly reproduces the broad maxima in the pDOS at 300-450 cm −1 for the LS crystal (corresponding to Fe stretching modes [38]) and the broad asymmetric peak at around 175-325 cm −1 for the HS phase, respectively. For the LS phase, additionally, the appearance of the two weak peaks at 150-250 cm −1 (Fe-N bending [38]) is reproduced for the crystal model. Interestingly, the Debye behavior of the phonons below around 50 cm −1 is also reproduced for the crystal simulation; this constitutes the main result of this work. The MD simulation predicts no less-intense peaks observed experimentally above 450 cm −1 . Instead, a weak cluster of peaks in the pDOS is predicted by MD simulations in the 550-600 cm −1 area. Also for the HS phase, the Debye region below 50 cm −1 is well reproduced, The remaining differences between experiment and MD may be attributed to two sources. (i) The force-field parameters used in the MD simulation were obtained from fitting to DFT calculations. However, the position of iron ligand modes can be simulated by DFT in our experience with an accuracy of ±30 cm −1 only. This accuracy is inevitably enhanced when using the fitted force fields in the MD simulations. (ii) In addition, the observed bands are a superposition of several molecular modes. Any inaccuracies in the mode frequencies and/or amplitudes will hence show up as further deviations between MD simulation and NIS experiment.
Further, the obtained pDOS of the LS and HS crystals were analyzed [39] to yield the lattice-dynamics parameters, including the Lamb-Moessbauer factor, f LM , total vibration amplitude x, mean internal energy U, specific heat C v , entropy S vib and normalized mean force constant D as a function of temperature. Here, f LM is a parameter of relevance to NIS experiments as it gives the fraction of -quanta emitted and absorbed without recoil during nuclear resonance experiments. A comparison of the relevant parameters obtained experimentally for the LS and HS   Table 1.
The results in Table 1 reveal that MD simulations yield a good agreement to the experimentally determined values of specific heat, vibrational entropy and mean internal energy. The estimation of the Lamb-Moessbauer parameter and total vibrational amplitude are more reasonable for the LS phase than for the HS phase. The values of the normalized force mean force constant follow only qualitatively the experimental results. The MD simulations lead to an overestimation of the mean force constant D both in the HS and the LS states. This might possibly be caused by an overestimation of packing interactions in the MD modeling.
Using the calculated values of the vibrational entropy for both spin phases, the vibrational spin transition entropy is displayed in Fig. 8. It has a value of around 10 J/(K mol) at the transition temperature (ca. 340 K) that is somewhat lower than the values 43-68 J/(K mol) determined calorimetrically for similar 1D systems [40].

Conclusions
We studied the vibrational spectra of one-dimensional molecular SCO crystals using molecular dynamics simulations. While quantum-chemical methods are able to provide data on isolated SCO molecules, the computational costs necessitate the use of classical simulations to calculate crystal effects on the spectra. We obtained the following findings.
1. We were able to calculate phonon densities of states of molecular SCO crystals in the low spin and the high spin state. The calculation of spectra of such large crystalline SCO materials was only possible using molecular dynamics in combination with a spin-state dependent atomistic potential. 2. The pDOS converges with increasing crystal size to a crystalline spectrum, both for the low-spin and the highspin state. As the crystalline structure is built up, the pDOS subsequently shows crystalline characteristics. 3. Increasing crystal size leads to a smoother spectrum as well as a broadening and subsequent lowering of the characteristic peaks in the spectrum. This can be attributed to bigger statistics with increasing system size as well as the crystalline structure. 4. The increase of the crystal size also leads to the appearance of a low-frequency peak in the spectrum with a quadratic increase at lower frequencies. This effect predominantly occurs at larger crystal sizes, when the spectrum is converging. 5. Good agreement between experimental and simulated phonon spectra is observed.
6. The lattice-dynamics analysis of the low-energy pDOS allows to extract the temperature dependence of the vibrational entropy of the LS and HS crystals.
In future work, it will be interesting to extend the studies to other SCO crystals, for which the experimental pDOS has been determined. These can be related one-dimensional SCO systems, such as [Fe(1,2,4-triazole) 2 (1,2,4-triazolato)] (BF 4 ) [21]; however, also more ambitious systems such as Acknowledgements We acknowledge support by the Deutsche Forschungsgemeinschaft -project number 268565370-SFB/TRR 173. Access to the computational resources provided by the compute cluster 'Elwetritsch' of the University of Kaiserslautern is appreciated. We also thank the Allianz für Hochleistungsrechnen Rheinland-Pfalz (AHRP) for providing CPU time within the projects TUKL-SPINPLUSVIB and TUK-SCO and Tim Hochdörffer for the calculation of the lattice-dynamics parameters.
Author Contributions All authors contributed to the study conception and design. Simulations and data analysis were performed by JH and RM. The first draft of the manuscript was written by HMU. All authors commented on the manuscript and read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interests
The authors have no competing financial or nonfinancial interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.