Electronic Excitations in Molecular Solids: Bridging Theory and Experiment

We evaluate the performance of continuum models in accounting for environmental influences on the component species in molecular crystals, using the [Ni(Et 4 dien)( η 2 -O,ON)( η 1 -NO 2 )] linkage-isomer system as a test case. As the spatial and temporal resolution accessible to experiment and theory converge, computational chemistry is an increasingly powerful tool for modelling and interpreting spectroscopic data. However, the study of molecular processes, in particular those related to electronic excitations 10 (e.g. photochemistry), frequently pushes quantum-chemical techniques to their limit. The disparity in the level of theory accessible to periodic and molecular calculations presents a significant challenge when modelling molecular crystals, since accurate calculations require a high level of theory to describe the molecular species, but must also take into account the 15 influence of the crystalline environment on their properties. In this article, we briefly review the different classes of quantum-chemical techniques, and present an overview of methods that account for environmental influences with varying levels of approximation. Using a combination of solid-state and molecular calculations, we quantitatively evaluate the 20 performance of implicit-solvent models for the [Ni(Et 4 dien)( η 2 -O,ON)( η 1 NO 2 )] linkage-isomer system as a test case. We focus particularly on the accurate reproduction of the energetics of the isomerisation, and on predicting spectroscopic properties to compare with experimental results. This work illustrates how the synergy between periodic and molecular 25 calculations can be exploited for the study of molecular crystals, and forms a basis for the investigation of more challenging phenomena, such as excited-state dynamics, and for further methodological developments. description of the valence wavefunctions. A kinetic-energy cutoff of 944.5 eV was used for the plane-wave 30 basis, and the electronic wavefunctions were evaluated at the Γ point; these

(e.g. photochemistry), frequently pushes quantum-chemical techniques to their limit. The disparity in the level of theory accessible to periodic and molecular calculations presents a significant challenge when modelling molecular crystals, since accurate calculations require a high level of theory to describe the molecular species, but must also take into account the 15 influence of the crystalline environment on their properties. In this article, we briefly review the different classes of quantum-chemical techniques, and present an overview of methods that account for environmental influences with varying levels of approximation. Using a combination of solid-state and molecular calculations, we quantitatively evaluate the 20 performance of implicit-solvent models for the [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] linkage-isomer system as a test case. We focus particularly on the accurate reproduction of the energetics of the isomerisation, and on predicting spectroscopic properties to compare with experimental results. This work illustrates how the synergy between periodic and molecular 25 calculations can be exploited for the study of molecular crystals, and forms a basis for the investigation of more challenging phenomena, such as excited-state dynamics, and for further methodological developments.

Introduction
Continuous improvements in spectroscopic techniques, such as the advent of next- 30 generation X-ray free-electron laser (XFEL) facilities, are allowing the structural dynamics of molecular crystals to be studied in unprecedented detail. These advances bring with them, however, a number of significant challenges in interpreting the data in terms of processes occurring at the molecular level. The length-and timescales accessible to these advanced experimental techniques are 35 rapidly converging on those which are accessible to theoretical study, and as such computational chemistry, e.g. within the ubiquitous density-functional theory (DFT) framework, can be a powerful complement to experiment.
A particular problem arises with the study of molecular solids, namely that the calculations need to capture both the chemistry of the molecular species, and the 40 influence of the crystalline environment. At present, there exists a large disparity between the level of theory which is accessible to molecular and solid-state (periodic) calculations, due to the higher complexity of the latter. However, accurate descriptions of molecular properties, particularly those related to electronic excitations, e.g. absorption spectra and excited-state dynamics, frequently require high-level methods. Through continual advances in computing power and software efficiency, the gap between molecular and solid-state calculations is closing, but at present molecular crystals still frequently stretch the limits of what is possible in 5 periodic calculations, mainly due to their large unit cells, and to the spectrum of weak non-bonding interactions which play an important role in defining their structure and properties.
To study electronic excitations in molecular solids, therefore, approximate methods of accounting for the effect of the crystalline environment in molecular 10 calculations are required. In this work, we review several different approaches, and evaluate quantitatively the performance of the simplest one for computing various properties of the well-studied [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] linkage-isomerisation system. The key focus of our work is to explore how best to exploit the synergy between solid-state and molecular calculations, and to provide new theoretical 15 insight to complement ongoing experimental work on this and related systems.

Towards Chemical Accuracy: The "Jacob's Ladder" of Approximations
Quantum chemistry aims to model the properties of quantum-mechanical systems from first principles by solving the Schrödinger equation, typically within the Born-Oppenheimer (clamped-nuclei) approximation. The ultimate goal is to be able to 20 predict properties with "chemical accuracy", that is, accuracy on the same scale as state-of-the-art experimental techniques. For example, historically, chemicallyaccurate energies are usually taken as implying an uncertainty of less than 1 kcal mol -1 (4.18 kJ mol -1 ), although tighter criterion are often required in practice. Due to Perdew and coworkers, 1 quantum chemistry techniques can be classified according 25 to a "Jacob's ladder" of approximations, extending to the "heaven" of chemical accuracy.
DFT, which recasts the many-body Hamiltonian in the electronic Schrödinger equation into a functional of the spatial electron density ݊(), 2, 3 is perhaps the most widely-used theoretical "workhorse" at present. Leaving aside implementation 30 details, such as the mathematical basis used to express the electron orbitals, the key parameter in DFT is the form of the exchange-correlation (XC) functional used to calculate the contributions to the total energy from quantum-mechanical electron exchange and correlation. The simplest functional form is the local-density approximation (LDA), in which the exchange and correlation energy for a given This journal is © The Royal Society of Chemistry [year] structure calculations on periodic systems.
Finally, there exist a growing number of "beyond DFT" and "post Hartree-Fock" methods, which are becoming increasingly popular and accessible with advances in computing power. Examples of this include GW theory, 4 a perturbative approach to treating many-body physics, and second-order Møller-Plesset perturbation theory 5 5 (MP2) and coupled-cluster theory, 6 which both aim at a more accurate description of electron correlation.
In practice, (meta-)GGA functionals typically offer a good balance between computational cost and accuracy for a number of properties, including energetics and forces, and hence equilibrium geometry and vibrational frequencies. 7 Hybrids 10 usually give more accurate electronic structures (i.e. orbital energies in molecular systems) than semi-local functionals, and so are often a prerequisite for the computation of optical properties and for studying electronic excitations in timedependent DFT (TD-DFT) calculations. Post-DFT methods are useful in cases where even higher accuracy is required, e.g. to construct benchmark datasets against which 15 functionals can be tested, and are sometimes necessary for systems where manybody effects are prominent, for which the more approximate approaches frequently fail.

Environmental Effects in Molecular Calculations: A Spectrum of Approaches
The most common methods for accounting for the solid-state environment in 20 calculations on the component species in molecular crystals may be divided into three classes, with differing computational cost and accuracy.
At one end of the scale, full periodic calculations can be performed. The molecular environment is treated explicitly, but the computational cost of using high-level theories is prohibitive for large systems, and lower-level theories may be 25 insufficient to describe certain properties with the required level of accuracy. At the other end, continuum models 8 (e.g. the polarisable-continuum model (PCM) 8 or COnductor-like Screening MOdel (COSMO) 9 schemes) assume that the most important environmental effects can be captured by a simple dielectric screening. This is likely to be a good approximation when the intermolecular interactions are 30 minimal, and is commonly used to implicitly model the effect of a solvent. Forming a "middle ground" between these techniques are embedding methods, in which a large system is separated spatially into different regions which are then treated at different levels of theory. 35

Explicit Periodic Calculations
With current hardware and software algorithms, it is feasible to perform periodic calculations on medium-sized molecular crystals (up to a few hundred atoms) using DFT with semi-local (meta-)GGA functionals, or, in some cases, with hybrids. [10][11][12][13] This is usually sufficient for studying equilibrium geometry and lattice dynamics 40 (e.g. vibrational frequencies), but semi-local functionals are often not able to describe electronic structure quantitatively, as is required, for example, for accurate prediction of optical properties.
A particular issue with molecular crystals, as opposed to many simple bulk materials, is that weak interactions (e.g. van der Waals' dispersion forces) often play 45 an important role in defining the structure. [14][15][16] Dispersion forces are, in principle, a non-local electron-correlation effect, and thus a first-principles quantum-chemical description requires a non-local functional. To circumvent this, several approximate This journal is © The Royal Society of Chemistry [year] methods have been developed to correct GGA energies and forces, e.g. the DFT-D2 17 /D3 18 and DFT-TS 19 /TS+SCS 20 methods, and non-local correlation functionals such as vdW-DF 21 /DF2. 22 However, these approximate forms may not account for more complex many-body dispersion interactions, which appear to be significant in some molecular systems. 16 5 A good compromise to obtain accurate electronic properties is to optimise structures at a moderate level of theory, and to then perform more accurate singleshot calculations with hybrid functionals. In many cases, this approach gives good results, but for systems where many-body effects are prominent, e.g. excitonic effects in optical spectra, 23 the required higher levels of theory sometimes cannot be 10 used due to computational limitations. Also, since geometry relaxation with nonlocal functionals is not feasible, it is not realistically possible to study excited-state dynamics with periodic calculations, which precludes the investigation of photochemical reactions. 15

Embedding
In typical embedding methods, a large bulk system is divided into regions, and the core (e.g. a single molecule, or a dimer, etc.) is treated with a quantum-mechanical (QM) method, while a surrounding "shell" is treated with an empirical molecularmechanics (MM) method such as a parameterised force field, reverting to a 20 continuum model at larger distances. The layers are then interfaced together to account for the interaction energy between them. The MM atoms are treated as point charges or multipoles, which can polarise the wavefunction in the QM region, while the QM atoms can likewise exert forces on the MM atoms. This allows for an explicit atomistic treatment of large systems, with high accuracy in a region of 25 interest, at a manageable computational cost. This QM/MM approach is commonly used to model biological systems, for example to look at redox processes at the active centres of enzymes.
This idea behind QM/MM embedding schemes can be extended to analogous QM/QM methods, where the core region is treated with a high-level quantum- 30 chemical method (e.g. MP2 or coupled-cluster) and the shell with a lower-level one such as DFT with a GGA functional. [24][25][26][27] This class of methods also encompasses techniques in which the periodic wavefunctions obtained from an extended-system calculation are converted to localised orbitals, and selected local states then treated with higher-level theories. 28, 29 35 Embedding schemes have successfully been applied to a variety of molecular crystals, [30][31][32] and the particularly ambitious study by Kochman et al. 32 utilised a scheme where a single molecule in a periodic DFT calculation was treated with a molecular code, allowing exploration of its excited-state potential energy surface using TD-DFT. 40

Continuum Models
At the lowest level of approximation, environmental effects can be included using continuum models, which assume that the main influence of surrounding molecules is a dielectric-screening effect. The key parameter in these models is the static 45 dielectric constant of the medium, ߝ ௦௧௧ , which captures its ability to screen an applied electric field. These methods are routinely used to model solvent effects in molecular calculations, but could also feasibly be used for molecular crystals, provided that the interaction between molecular units was negligible.

This journal is © The Royal Society of Chemistry [year]
The molecule is treated as being a solute within a cavity, surrounded by a dielectric continuum of the solvent. The charge-density cloud of the molecule polarises the medium, which responds by generating screening charges on the cavity surface. In the COSMO model, the screening charges are obtained by modelling the response of an ideal solvent with ߝ ௦௧௧ = ∞, which is then scaled to account for the 5 ߝ ௦௧௧ of a non-ideal medium. Some continuum models additionally contain the shape of the solvent molecule as a parameter, which is then taken into account when defining the cavity around the solute.
Continuum models are widely used in the study of solution chemistry (e.g. solution thermodynamics 33 and reactions 34,35 ), and to model biological systems. 36 10 They have also been used to model environmental effects on electronic excitations, 37, 38 e.g. within a PCM/TD-DFT formalism. 38 Continuum models are generally used in quantum chemistry to model solvents, although there are a few cases where they have been employed to model the dielectric environment of a crystal in a molecular calculation. 39, 40 15

Test System: Linkage Isomerisation in [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )]
Solid-state linkage isomerisation is a topical example of a reversible phase transition generating one or more metastable structural isomers in response to external stimuli, e.g. illumination or temperature changes. 41  complexes are perhaps the most well-known. Linkage-isomer systems are typically studied using photocrystallography, where a sample, usually a single crystal, is irradiated in situ on the diffractometer, allowing the steady-state populations of the different isomers to be quantified as a function of 25 temperature and of the wavelength of radiation used. 52 By performing pseudosteady-state experiments, where the crystal is continuously irradiated during the data-collection process, 41 additional short-lived metastable species can also be identified. 46 The [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] linkage-isomer system is an ideal test case 30 for continuum models. The isomerisation takes place within a large reaction cavity, which allows the ligand binding to be switched without inducing a large stress on the crystal. 39,40,43,51 Furthermore, previous computational modelling 40 and experimental measurements of the transition kinetics 46 both suggest that, to a good approximation, the molecules in the solid behave as isolated units, whose properties 35 and behaviour are influenced by the dielectric environment of the crystal. This system exists as three isomers, in which the η 1 -bound NO 2 ligand coordinates to the Ni centre via either N or O (Fig. 1). The N-bound ground-state (GS) isomer is obtained on cooling the crystal in the absence of illumination, and is the thermodynamically-stable form. The GS can be converted to the O-bound MS1 40 isomer with excellent yield on photoactivation, and is also formed thermally in significant population at moderate temperatures. A second O-bound isomer, MS2, was recently observed as a transient species close to the so-called metastable limit, the temperature at which MS1 begins to decay thermally, and can interconvert with MS1 with a relatively small energy barrier.  A number of open questions remain on linkage isomerisation. Firstly, although 5 the isomers themselves have been characterised, understanding the reaction paths which connect them is a work in progress. [53][54][55] Numerical modelling has allowed mechanisms to be proposed for some systems, 39,40 although explicit study of the excited-state potential energy surfaces, and establishing the photochemicalisomerisation pathways, remains an important undertaking. A second key challenge 10 is to understand the role that the crystalline environment plays in controlling the kinetics and energetics of the process, which would, in principle, provide valuable input to crystal-engineering approaches for tuning the properties of different systems for specific applications.
In the present work, we focus on the GS and MS1 isomers, since these have been 15 well characterised experimentally, 44,46 providing a good set of reference data against which to optimise computational parameters and to evaluate the performance of continuum models.

Methods
All computational work was performed within the Kohn-Sham DFT formalism. 3 20 Different codes were used for the periodic and molecular calculations, as outlined below.

Periodic Calculations
Periodic plane-wave pseudopotential calculations were performed using the Vienna ab initio simulation package (VASP) code. 56 Projector augmented-wave (PAW) 25 pseudopotentials 57, 58 were used to treat core electrons, with the outer s and p electrons of H, N and O, and the 4s, 3d and semicore 3p states of Ni, being treated as valence states. Where available, we used "hard" pesudopotentials, with a smaller core region to allow for more flexibility in the description of the valence wavefunctions. A kinetic-energy cutoff of 944.5 eV was used for the plane-wave 30 basis, and the electronic wavefunctions were evaluated at the Γ point; these parameters were found to be sufficient to converge the total energies to within less than 1 meV per atom. Spin-polarisation was used in all calculations, with each Ni complex constrained to have a triplet magnetic configuration, as was established to be the lowest-energy state in previous work. 40 A selection of DFT exchange-correlation (XC) functionals were used in these 5 calculations, as described in the text: we tested the PBE 59 and PBEsol 60 GGA functionals, the TPSS 61 meta-GGA, the dispersion-corrected PBE-D2, 17 PBE-D3 18 and TPSS-D2 functionals, and the PBE0 hybrid. 62 Convergence criteria of 10 -6 and 10 -5 eV were employed during electronic minimisation and atom position/unit-cell parameter optimisations, respectively. 10 During the calculation of the dielectric functions, as described in the text, the number of electronic bands was increased to 1,136, triple the default, to ensure convergence of the sum over empty electronic states. To obtain a high resolution for the calculated function, the number of grid points used to evaluate the electronic density of states (DOS) was increased to 5,000 or 10,000, depending on whether or 15 not the resulting function was to be used to compute optical properties.

Molecular Calculations
Molecular calculations were carried out using the NWChem 63 and Gaussian 09 64 codes. As described in the text, several Pople split-valence 65 and Dunning 66 basis sets were tested for the light atoms, together with the Los Alamos effective-core 20 pseudopotential and corresponding double-zeta basis set (LANL2DZ) for the Ni atom. 67 As in the VASP calculations, the molecular complex was constrained to have a triplet magnetic configuration. For comparison with the periodic calculations, we tested PBE, PBE-D2 and TPSS, and we also used the M06 68 and B3LYP 69 functionals, which are both popular choices for molecular calculations. Convergence 25 criteria of 10 -6 and 10 -5 were used for the total energy and density, respectively, while geometry optimisations were performed to a force tolerance of 5 x 10 -4 Hartree Bohr -1 . For the continuum calculations, we tested the COnductor-like Screening MOdel (COSMO) method, 9 as implemented in NWChem, and the polarisable-continuum model (PCM) in Gaussian. 30

Vibrational Spectroscopy
Room-temperature infrared (IR) and Raman spectra of crystalline [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] were collected for comparison with calculated vibrational spectra. IR spectra were obtained from a single crystal using a Perkin Elmer Spectrum 100 FT-IR spectrometer with the ATR accessory. A spectral resolution of 1 cm -1 between 35 4000 and 600 cm -1 was available with this instrument, and acquisition and processing was performed using the Perkin Elmer Spectrum software.
Raman spectra were collected from a single crystal using a Renishaw inVia Raman spectrometer with the Renishaw WiRE 4.1 software. A 532 nm excitation laser (Renishaw diode laser, 380 mW at source/1.9 mW after attenuation) was used 40 with a 20x objective lens (LEICA), and the instrument was calibrated to the 520 cm -1 line of Si. Spectra were obtained between ~100 and 4000 cm -1 at a resolution of 1.2 cm -1 .
It is worth noting that, since the GS-to-MS1 transition is thermally activated, at ambient temperature, and in particular under laser irradiation at 532 nm during the 45 Raman experiments, the spectra are expected to contain bands due to both the GS and MS1 forms.

This journal is © The Royal Society of Chemistry [year]
3 Results

Energetics
Among the most straightforward property to compute with DFT is (relative) total energy. The per-molecule enthalpy difference between the GS and MS1 crystals has been determined from photocrystallographic measurements of the temperature 5 dependence of the isomer populations, 44 providing a benchmark against which to quantitatively compare calculated energy differences.
In our recent work, 40 we found that the energetics were sensitive to the computational parameters employed, in particular the DFT functional used and, for molecular calculations, the choice of basis set. As a foundation for testing 10 continuum models, in this section we systematically investigate these dependencies, and quantify the effect of the crystalline environment on the energy differences between the isomers.
To test the effect of the choice of functional on the predicted energetics, we carried out plane-wave pseudopotential calculations on the reported GS and MS1 15 crystal structures 44 with six commonly-used functionals, viz. PBE, PBEsol, PBE+D2, PBE+D3, TPSS and TPSS+D2. For each, the atomic positions and unitcell parameters were fully relaxed, and the per-molecule GS-MS1 energy differences were then calculated from these optimised models; some structural parameters from the models (e.g. lattice parameters, cell volume and bond lengths) are given in the 20 supporting information. To investigate the effect of the crystalline environment on the energetics, we also carried out an equivalent set of calculations on the individual GS/MS1 complexes in the gas phase. These were performed by placing single molecules from the crystal structures in periodic cells, with an (initial) vacuum gap of 10 Å between images, and optimising the geometries. 25 Of the selection of functionals we employed, PBE represents a typical choice for many solid-state problems, although PBEsol, a variant of PBE revised to better reproduce the properties of solids, is also used routinely. PBE+D2 and PBE+D3 both add a semi-empirical correction to the PBE energies to account for dispersion forces, and are frequently used for systems where weak interactions are expected to be 30 significant. TPSS is a typical meta-GGA functional, and perhaps the most routinely used functional of this type, and we also tested it in combination with the D2 dispersion correction (TPSS+D2). Fig. 2 compares the per-molecule energy differences in the crystal and in the gas phase. Compared to the experimental enthalpy difference of 9.69 kJ mol -1 , all of the 35 functionals apart from TPSS significantly overestimate the energy difference between the GS and MS1 crystals. The PBE energy difference of 14.57 kJ mol -1 is around 1.5× larger than the experimental one, while PBEsol and the three dispersion-corrected functionals overestimate the difference by a factor of two. On the other hand, the TPSS value of 7.75 kJ mol -1 is in remarkably good agreement, 40 suggesting that the improved accuracy of meta-GGA functionals provides a good description of this system.  Comparing the gas-phase energy differences to those in the solid, all six 5 functionals predict that the difference between the isomers is substantially reduced. However, in all cases it remains positive, indicating that the GS isomer is the ground state in the gas phase as well as in the molecular crystal.
Having found that TPSS gives improved energetics compared to the GGA functionals, we performed single-point energy calculations on our six optimised GS 10 and MS1 structures using the PBE0 hybrid functional, and recalculated the energy differences (see supporting information). We observed some small variation between the different starting geometries, although the calculated values appear to be relatively insensitive to this. The best PBE0 value (7.05 kJ mol -1 ) is obtained when the atomic positions, but not the cell parameters, are optimised with TPSS as a 15 starting point; it is possible that the agreement might improve further if a geometry optimisation could be performed with the hybrid.
The tight convergence criteria employed in the plane-wave calculations on the molecular species means that the calculated energy differences should be close to the basis-set convergence limit for the molecular calculations. We therefore used 20 these values as a reference to optimise the basis set for these simulations. We carried out single-point energy calculations on the optimised gas-phase GS and MS1 structures with a number of different basis sets, and compared the calculated energy differences to the plane-wave values. Table 1 shows the results for a subset of the basis sets tested and the PBE functional; additional data for PBE, and a 25 corresponding set of data for PBE+D2, are given in the supporting information.
In general, the double-zeta 6-31G family of basis sets predict a qualitatively incorrect energy ordering. Among the triple-zeta basis sets, the accuracy of the predicted energy differences improves when diffuse and polarisation functions are added (6-311G+/++ and 6-311G(d)/(d,p), respectively). This appears to be more 30 important for the heavier atoms, with the 6-311+G(d) basis, which includes diffuse functions and d-type polarisation functions on the heavy atoms, giving very similar results to 6-311++G(2d,2p), which includes diffuse functions and two d/p This journal is © The Royal Society of Chemistry [year] polarisation functions on both light and heavy atoms. Table 1 Calculated single-molecule GS-MS1 energy differences with the PBE functional and various basis sets. The structures of the complexes are those obtained from the gas-phase planewave calculations, and are not re-optimised in these tests. 5 Basis Set ∆EMS1-GS / kJ mol -1 molec -1 We note that the results in Table 1 were all obtained using the LANL2DZ pseudopotential and corresponding double-zeta basis set to treat the Ni atom; during testing, we found that varying the Ni basis made relatively little difference to the energetics when the better-converged basis sets were used for the other atoms (see 10 supporting information), but that, without the pseudopotential, the computational cost was significantly increased.  Table 1, the structures of the complexes 15 are those obtained from the corresponding gas-phase plane-wave calculations, and are not reoptimised. For comparison, the reference plane-wave values are overlaid as dashed lines.
In contrast to the Pople basis sets, the Dunning bases are designed to give more systematic convergence with respect to the number of basis functions, although they are typically more computationally expensive than the Pople ones. Like the 20 equivalent Pople double-zeta basis sets, the double-zeta cc-pVDZ basis predicts an incorrect energy ordering, while the augmented cc-pVDZ (aug-cc-pVDZ) basis set and the bare and augmented triple-zeta (aug-)cc-pVTZ sets give the best results overall. For the augmented cc-pVDZ and the two cc-pVTZ basis sets, we observed near-quantitative agreement with the plane-wave values for energetics calculations with PBE, PBE+D2 and TPSS (Fig. 3), suggesting that the energy difference with these basis sets is close to the convergence limit. 5 We note that these basis sets are computationally expensive to work with, and we found them unwieldy for geometry optimisation; however, the results in Table 1 suggest that a good compromise would be to optimise with the more efficient 6-311+G(d) basis, and to then perform further calculations with one of the triple-zeta Dunning bases where possible, which is the approach taken in the following 10 sections.

Dielectric Properties
In this section, we calculate the macroscopic dielectric constants of the molecular crystals from solid-state calculations, and evaluate the performance of the implicitsolvent COSMO model in including the effect of the crystalline environment in the 15 molecular calculations.
As discussed in Section 1, the static dielectric constant (permittivity), ߝ ௦௧௧ , of a material measures its ability to screen an applied electric field, and as such can be used to model approximately the effect of a dielectric continuum (e.g. a solvent) on 20 various properties of a molecular system. ߝ ௦௧௧ can be decomposed into the sum of two contributions, viz. electronic polarisation, and ionic relaxation: 70 can be calculated either from the response of the system to a finite electric field (e.g. using density-functional perturbation theory (DFPT)), or obtained as the zero-frequency component of the real part of the frequency-dependent dielectric function ‫)ܧ(ߝ(‬ = ߝ ‫)ܧ(‬ + ݅ߝ ‫.))ܧ(‬ Calculating ߝ requires additionally the vibrational modes of the system to be evaluated, which for large systems can be 30 very time consuming.
In the current version of the VASP code, ߝ ௦௧ can be computed using DFPT, or by calculating the dielectric function via the linear optical response, while both components of ߝ ௦௧௧ can be calculated by analysing the vibrational modes using DFPT. However, VASP currently does not support DFPT calculations with 35 semi-empirical and meta-GGA functionals, and thus we were only able to calculate both components of ߝ ௦௧௧ with PBE and PBEsol. For the other functionals, we obtained the value of ߝ ௦௧ by calculating the dielectric function. Table 2 lists the values of ߝ ௦௧ , ߝ and ߝ ௦௧௧ for the GS and MS1 crystals obtained with PBE and PBEsol, and the corresponding values of ߝ ௦௧ computed with a 40 selection of other functionals are listed in Table 3.
There is a clear overall trend visible in this data. For both crystals, ߝ ௦௧௧ is around 3-4, with the major contribution being from electronic polarisation rather than ionic relaxation. This value is comparable to a solvent such as benzene, toluene or diethylamine (2.27, 2.38 and 3.58, respectively), rather than a more polar solvent 45 such as ethanol or water (24.5/80.1). On this scale, the variation between the This journal is © The Royal Society of Chemistry [year] functionals considered here is relatively small, which, in the absence of experimental measurements, lends a reasonable degree of confidence to the computed range. Having obtained a reliable estimate of the dielectric constants of the GS and MS1 crystals, we then attempted to reproduce the solid-state energetics in molecular calculations using COSMO. We calculated the GS-MS1 energy differences for molecular complexes optimised with a selection of functionals, viz. PBE, PBE-D2, 15 TPSS, B3LYP and M06, and with dielectric constants of ߝ ௦௧௧ = 3, 3.5 and 4; as discussed at the end of the previous section, the geometry optimisations were performed with the 6-311+G(d) basis set, and the final energy differences were obtained from single-point energy calculations with the cc-pVTZ basis.
The calculated energy differences are compared in Fig. 4. PBE, PBE-D2 and 20 TPSS were also used in the periodic calculations, allowing for a quantitative comparison. For all three functionals, the presence of a continuum increases the energy difference relative to the gas phase, mirroring the plane-wave calculations, and in general the correspondence between the energy differences obtained from the full periodic and continuum calculations is very good. The only exceptions to this 25 are the PBE-D2 results, the reasons for which are not clear. B3LYP and M06 are hybrid and meta-hybrid functionals, respectively, and thus should in principle yield more accurate energetics than the others. M06 performs similarly well to TPSS in these calculations, yielding energy differences fairly close to the experimental values. On the other hand, B3LYP significantly underestimates 30 the energy differences, yielding values lower than the TPSS gas-phase results, which suggests that this functional does not provide a good description of this system. We note in passing that we encountered problems with spin contamination in the This journal is © The Royal Society of Chemistry [year] B3LYP/ߝ ௦௧௧ = 3.5 single-point calculation with the cc-pVTZ basis set, and for the same reason we were not able to obtain energy differences with PBE0. Since this problem did not appear to affect the gas-phase calculations with the same functionals, and the small dielectric constant of the continuum is not expected to lead to large perturbations to the orbital structure, this is most likely an issue with 5 the current implementation of COSMO in NWChem.  10 included for comparison, and the experimentally-determined energy difference is overlaid as a dashed black line. As discussed at the end of the section on energetics, for each calculation, the geometries of the GS and MS1 molecules were optimised with the 6-311+G(d) basis set, and the energies then recalculated with the cc-pVTZ basis set. There is no B3LYP value for ߝ ௦௧௧ = 3.5, due to issues with this particular calculation (see text). 15 Overall, the results suggest that a polarisable-continuum model, using the dielectric constant obtained from periodic calculations, is an effective way to combine the strengths of solid-state and molecular calculations for this system. In the following sections, we use our optimised parameters to model various spectroscopic properties, in particular IR/Raman frequencies and UV-visible 20 absorption profiles.

Vibrational Spectra
Since (meta-)GGA functionals typically yield good forces, it is possible to perform periodic vibrational-frequency calculations, although this can be expensive for large systems. However, while calculating the IR intensity for a vibrational mode is fairly 25 straightforward, Raman intensities are more computationally demanding to model, and so for molecular crystals this property is ideally best obtained from molecular, rather than periodic, calculations.
The IR spectra of the molecular crystals were obtained with the PBE and PBEsol functionals as a by-product of the calculations of ߝ in the previous section. To 30 further quantify the performance of the polarisable-continuum model, we computed the spectra with PBE and a 6-311+G(d) basis set using the Gaussian 09 code, using the polarisable-continuum model 8 with diethylamine (DEA) as a solvent (ߝ ௦௧௧ = This journal is © The Royal Society of Chemistry [year] 3.6). The spectra are compared in Fig. 5. The overall correspondence between the two sets of spectra is very good. Although there are minor variations in the calculated intensities, the spectra have similar overall forms, and most of the calculated peak positions are likewise very 10 similar. The only notable exception is in the position of a pair of bands around ~1400 cm -1 , which correspond to a mixture of vibrations of the isomerisable ligands and parts of the Et 4 dien backbone, and which lie on top of each other in the periodic calculations, but are visibly separate in the continuum ones. This suggests a slight difference in the bonding of the ligand in the two sets of calculations, which may 15 well be related to the small energy differences between the periodic and continuum calculations evident in Fig. 4.
It is also worth noting that, despite the use of the continuum model, as only a single molecule is present, the molecular calculations would not be able to reproduce phonon modes corresponding to collective motions of two or more of the four 20 molecules in the crystallographic unit cell, which could potentially lead to discrepancies, in particular in the low-wavenumber ("fingerprint") region of the spectra. The good correspondence between the periodic and molecular spectra in Fig. 5 therefore indicates that, to a very good approximation, the vibrational modes of the molecular crystal correlate with those of the molecular species, lending 25 further support to the treatment of the molecules in the solid as independent entities.
Given the good performance of the continuum model with the calculated dielectric constant, we opted to compute IR and Raman spectra for both isomers with the more accurate M06 functional, which, as a hybrid functional, is in particular more likely to predict Raman polarisabilities more accurately than (meta-)GGAs. The calculated spectra are compared against room-temperature IR and Raman data collected from single crystals in Fig. 6.   Fig. 6 Simulated IR and Raman spectra of the GS (blue) and MS1 (red) molecules, compared to 5 room-temperature (~300 K) spectra taken from single crystals (black). The simulated spectra were obtained using the M06 functional and a polarisable continuum of diethylamine (DEA; ߝ ௦௧௧ = 3.6). As in Fig. 5, these spectra were broadened using a Lorentzian function with a full-width at half-maximum (FWHM) of 7.5 cm -1 .
There is a generally good correspondence between the simulated and experimental 10 spectra. There are some notable mismatches in the predicted intensities, particularly in the case of the Raman spectra, although this could be due in part to the (arbitrary) choice of the broadening width for the calculated spectra. However, the positions of the peaks appear to be fairly well modelled, enabling assignment of some of the IR bands between ~1000 and 1500 cm -1 to the GS and MS1 isomers. This agreement 15 nicely illustrates a potential utility of computation for supporting characterisation.

Optical Properties
Optical properties arise from electronic excitations, and are a response property of the system. Modelling them formally requires solution of the time-dependent Schrödinger equation, or, more practically, modelling a perturbation to the ground- 20 state electronic wavefunctions in response to an external field (e.g. as in the TD-DFT formalism).
The method implemented in VASP uses a perturbative linear-response method to obtain the frequency-dependent dielectric function, ‫,)ܧ(ߝ‬ which, in essence, entails an enumeration over transitions between filled and empty electronic states. 71  extinction coefficients (ߢ), and absorption coefficients (ߙ), the latter of which can be compared to experimentally-recorded UV-visible spectra. Fig. 7 compares these quantities for the GS and MS1 structures, computed from the TPSS-optimised crystal structures with a PBE0 ground-state reference. 5 Fig. 7 Optical properties of the GS (blue) and MS1 (red) isomers, computed from the TPSSoptimised crystal structures using a PBE0 ground-state reference. The plots compare the real (solid lines) and imaginary (dashed lines) parts of the frequency-dependent dielectric function ‫,)ܧ(ߝ‬ ߝ /ߝ , the refractive index, ݊, the extinction coefficient, ߢ, and the absorption coefficient, ߙ. Each of the latter three plots are annotated with the equations used to calculate the quantities from ‫.)ܧ(ߝ‬ 10 The UV-visible spectrum of the GS molecular crystal published previously 46 does not extend into the deep UV, but the key features in the visible region are a long absorption tail around 400 nm, plus a second, broad isolated peak at ~650 nm. The most likely comparison between this and the spectrum in Fig. 7 is that the absorption tail is blue-shifted by ~200 nm, while the peak at 650 nm is not reproduced. 15 For comparison, we calculated spectra for the GS and MS1 molecules using the linear-response TD-DFT implemented in Gaussian. As for the calculations in the previous section, we used the M06 functional, together with a 6-311+G(d) basis set and a polarisable continuum of diethylamine. To compute the spectra, we analysed the 30 lowest-energy triplet transitions. The calculated spectra are shown in Fig. 8. 20 The overall structure of both spectra is similar to the 150-250 nm region of the solid-state one, consisting of primary high-and low-intensity peaks with some fine structure. The red shift of ~200 nm compared to the periodic spectra suggests that the more sophisticated TD-DFT method implemented in the molecular code gives a better description of the optical-absorption profile, and the position of the absorption 25 tail agrees much better with the experimental data.
Based on our previous calculations, 40  delocalised metal-to-ligand transitions, with large components on the Ni d orbitals, while the majority of the states from ~450-350 nm can be assigned as metal-toligand charge-transfer bands. We note that, while the computational parameters used in these and in the present calculations are very similar, the continuum solvent used in the previous study (H 2 O) has a considerably larger dielectric constant (ߝ ௦௧௧ = 5 80.1). However, we found that, while the positions of the absorption bands were sensitive to the dielectric constant, their nature and relative ordering were not, with the assignments being the same both in the gas phase and in a continuum of H 2 O. The broad peak at 650 nm is again not reproduced in these spectra. However, considering the individual calculated transitions (see supporting information), there are several states with lower excitation energies, but with negligible oscillator 15 strengths. The orbitals involved have a significant Ni d component, suggesting them to be (symmetry forbidden) d-d transitions, which perhaps become weakly allowed when coupled to vibrations. This could, in principle, be modelled by performing TD-DFT calculations on structures displaced along the vibrational modes, although this is much more computationally demanding. 20 These results, and also those in our previous work, 40 suggest that further optimisation of the computational parameters may be required to reproduce accurately the experimental UV-visible spectra for this system. However, the superior performance of the continuum TD-DFT calculation over the periodic linearresponse method highlights the potential benefits of being able to use molecular 25 calculations, and their excited-state functionality, while accounting approximately for the influence of the crystalline environment.

Discussion and Conclusions
The results in the previous section demonstrate quantitatively the performance of polarisable-continuum models in approximately accounting for the influence of the dielectric constants obtained from periodic calculations, were able to reproduce a broad range of properties of the [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] system, including energetics, vibrational frequencies, and, to a reasonable extent, optical properties. The correspondence between the continuum and periodic calculations was very good, and the comparison with available experimental data was also generally 5 favourable. Moreover, and particularly with regard to modelling electronic excitations, the benefits of the higher levels of theory accessible in molecular calculations were clearly demonstrated.
The method explored in this work is fairly straightforward, and effectively combines some of the strengths of periodic and molecular calculations to allow for a 10 more complete theoretical characterisation of molecular solids. For its broader applicability, the most important proviso is the assumption inherent in the polarsisable-continuum approach, i.e. that the interactions between the molecular units in the solid are well represented by a dielectric-screening effect. It remains to be discussed, however, for which sorts of system this is a good approximation, and 15 for which systems it should be expected to fail, e.g. crystals with longer-range intermolecular interactions such as π-stacking. For systems where this method works well, however, we anticipate that theoretical calculations should be a valuable tool for modelling and understanding a broad spectrum of state-of-the-art experimental work. 20 While the main focus of this study has been on the methodology, the calculations also provide some interesting insight into the chemistry of [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )], reinforcing the findings of our previous study. 40 It is quite clear that the crystalline environment significantly influences the relative stabilities of the different linkage isomers, in this case increasing the energy difference between the 25 GS and MS1 forms relative to the gas phase. To explore this further, as a theoretical exercise we calculated the GS-MS1 energy difference using COSMO for a range of dielectric constants, spanning the two-orders-of-magnitude range covered by common solvents (Fig. 9). This plot clearly illustrates that, by hypothetically tuning the dielectric constant of the crystal, the energy difference between the two isomers can be adjusted over quite a large range. In general, a stronger dielectric screening appears to increase the energy difference, with the effect saturating beyond values of ߝ ௦௧௧ around 30. Up to this limit, however, there is potentially scope for controlling the kinetics of the 5 isomerisation process through the chemistry of the molecular species. This adds a new perspective to the established practice of using crystal engineering to create a "reaction cavity" for the isomerisation to occur within. 39,40,43,51 To investigate this finding in more detail, it would be interesting to calculate the dielectric constants of other known Ni-NO 2 linkage-isomer systems, and to relate 10 these to the photoconversion yield and/or the presence of absence of thermalisomerisation pathways. It may perhaps also be of interest to relate the dielectric constant to other properties, such as the optical-absorption profile, and the energy difference between magnetic configurations.
To conclude, polarisable-continuum models represent an effective means of 15 including the influence of the crystalline environment in molecular calculations on the [Ni(Et 4 dien)(η 2 -O,ON)(η 1 -NO 2 )] complexes. We have compared quantitatively a number of properties, and observed good overall agreement both with periodic calculations and experimental data, as well as obtaining new insight into the properties of this system, in particular the effect of dielectric environment on the 20 linkage isomerism. This work represents a step towards developing more general approaches to the theoretical characterisation of molecular solids, and, ultimately, to working on more complex topics, such as performing high-level electronic-structure calculations and exploring excited-state potential-energy surfaces.
This journal is © The Royal Society of Chemistry [year]