Accurate force fields and methods for modelling organic molecular crystals at finite temperatures †

We present an assessment of the performance of several force fields for modelling intermolecular interactions in organic molecular crystals using the X23 benchmark set. The performance of the force fields is compared to several popular dispersion corrected density functional methods. In addition, we present our implementation of lattice vibrational free energy calculations in the quasi-harmonic approximation, using several methods to account for phonon dispersion. This allows us to also benchmark the force fields’ reproduction of finite temperature crystal structures. The results demonstrate that anisotropic atom–atom multipole-based force fields can be as accurate as several popular DFT-D methods, but have errors 2–3 times larger than the current best DFT-D methods. The largest error in the examined force fields is a systematic underestimation of the (absolute) lattice energy.


Introduction
Intermolecular interactions are central to many fields of science, 1 such as supramolecular chemistry, drug formulation and crystal engineering, 2 for which computational methods have been developed with the aim of ab initio crystal structure prediction (CSP). 3mputational modelling of molecular crystals with force fields has a long and rich history, [4][5][6][7][8] but recently we have seen a trend favouring density functional theory (DFT) calculations, perhaps most notably in the latest blind tests of crystal structure prediction, [9][10][11] where accurate lattice-and free energy calculations are a cornerstone of current prediction methods. 9,12,13ystal structure prediction is, in part, motivated by the need to anticipate polymorphism, where a molecule adopts different crystal structures under different crystallisation conditions.Differences in crystal packing between polymorphs can have important impacts on physical properties.The lattice energy difference between most polymorph pairs is smaller than 2 kJ mol −1 , 14,15 so that many-body dispersion, 16,17 induction, 18,19 lattice vibrations 14,20 and pressure 21,22 can affect the relative stability of polymorphs, in principle necessitating a calculation of the Gibbs free energy and very intricate computational modelling.Also, since the relevant region of predicted crystal energy landscapes typically contain tens, or even hundreds of plausible polymorphs, 23,24 the computational cost associated with evaluating a School of Chemistry, University of Southampton, Southampton, UK; Tel: +44 2380599174; E-mail: g.m.day@soton.ac.uk † Electronic Supplementary Information (ESI) available: Additional results, a description of co-prime splitting and Debye calculations, expressions for thermodynamic properties and phonon densities of states.See DOI: 10.1039/b000000x/ each crystal structure must be kept reasonable.
6][27][28][29][30][31] However, it is difficult to develop and rigorously validate computational models for molecular crystals, because of the limited amount of available high-accuracy experimental reference data on the properties of molecular crystals.Experimental sublimation enthalpies are typically limited in accuracy to ±4.9 kJ mol −1 . 32Consequently, the accuracy of computational methods cannot usually be assessed to better than this margin of error.
The benchmark study by Otero-de-la-Roza and Johnson 33 compared several common periodic DFT-D methods and a recently developed exchange-hole dipole moment dispersion correction method (XDM) using a set of 21 crystal structures (C21) of small organic molecules for which accurate crystal structures and sublimation enthalpies are known.Temperature-free benchmark lattice energies were obtained by calculating and subtracting thermal contributions.
The C21 benchmark set was corrected and extended by Reilly and Tkatchenko to form the X23 benchmark, 34 which has rapidly become quite popular.The benchmark set is currently restricted to single component crystal structures of small organic molecules; a future extension to multi-component crystals, such as salts and co-crystals would be useful to test methods on these systems, which have industrial relevance and where a different balance of interactions might be present.
Several DFT methods and dispersion correction schemes have been assessed using the X23 set.These include the GGA func- 5][36][37] In addition, the PBEh-3c functional and density functional tight binding (DFTB) have been benchmarked against X23. 37,38ince the interactions between organic molecules are relatively weak, temperature changes can have a significant influence on crystal structures and their properties.Therefore, it is desirable to be able to model molecular crystals at specific temperatures.The quasi-harmonic approximation 39 is commonly used to this end.However, when using force fields fitted to reproduce experimental structures, the intermolecular interaction parameters absorb some of the effects of lattice dynamics.Thus, including temperature explicitly can result in a double-counting of the thermal expansion and zero-point energy.In an effort to determine how large this error is, we also compare computationally thermally expanded crystal structures to finite-temperature experimental data.
In this study we use the X23 benchmark to assess the performance of several closely related empirically parameterised intermolecular atom-atom force fields, all derived from the work of Donald Williams, [40][41][42][43][44] when coupled with a DFT description of intramolecular energies.In essence, we present computationally very efficient methods for modelling crystal structures at finite temperatures with highly accurate force fields.

Force fields and crystal structures
In this benchmark study we have assessed the performance of the following force fields in reproducing the X23 benchmark lattice energies: • FIT: The FIT force field, 41 which is used with electrostatics described by distributed, atomic multipoles.
• W99rev6311: A revised version of the W99 force field, 45 which was parameterised to improve hydrogen bonding when using atomic multipolar electrostatics.
• W99rev6311P3: A revised W99 force field, 45 parameterised to be used with multipoles from a charge density calculated within a polarisable continuum model with dielectric constant ε r = 3.0.
• W99rev6311P5: The same as W99rev6311P3, but using ε r = 5.0 during the PCM calculation for the electrostatic model.
The force fields we benchmark here are all hybrid models that utilise a DFT model for the molecular geometry and energy, along with an intermolecular atom-atom model potential.The different force fields are strongly related and have the same functional form; an exp-6 potential for non-electrostatic intermolecular interactions and a separate treatment of electrostatics; either with atomic point charges or atomic multipoles.
All force fields included here use transferable, empirically parameterised repulsion-dispersion models.
Thus, no reparameterisation of the non-electrostatic interaction model is required when applied to a new molecule.The electrostatic model is generated from a quantum chemical calculation through either fitting to the molecular electrostatic potential or partitioning of distributed multipoles.In their original forms, FIT, W99 and W99rev provide repulsion-dispersion parameters for carbon, hydrogen, oxygen and nitrogen, though they have been used successfully with sulfur and halogen parameters taken from other sources. 468][49] Therefore, we have tested W99 with both multipoles (DMA) and molecular electrostatic potential derived atomic charges (ESP).
The newer revised versions of the W99 force field were partially re-parameterised to perform optimally with distributed multipoles. 45Parameterisation was performed against a large set of low temperature crystal structures to produce a force field whose parameters have absorbed as little thermal expansion as possible.Different versions of the force field were parameterised for use with multipoles derived either from a ground state charge density of the isolated molecule (W99rev6311), or from a polarised charge density obtained by embedding the molecular calculation in a polarising environment (W99rev6311P).It has been suggested 19 that the molecular DFT calculations should be performed within a polarisable continuum model 50,51 (PCM) with a dielectric constant of about 3, which is a typical permittivity of organic crystals.This should to some extent model an average polarisation of the molecular electron density in the crystal structure, and the accompanying induction energy.Polarised multipoles can then be derived from the polarised charge density.We can think of this as a mean field approximation of molecular polarisation in crystals, ignoring the local polarisation of functional groups arising from specific intermolecular interactions, particularly around hydrogen bonds.If a force field is parameterised specifically to be used with PCM-derived multipoles, the localized polarisation will be partly absorbed in the repulsion-dispersion parameters.This approach was attempted in parameterising the W99rev6311P force field.A relative permittivity of 3.0 was used in the parameterisation, but in our experience higher values can lead to improvements.Here we present results for ε r = 3.0 and 5.0.
We also include the older FIT potential, 41 which is a revision of the W84 force field, 40 again re-parameterised to perform optimally with atomic multipole electrostatics.
The force fields investigated here were parameterised with foreshortened hydrogen positions.Foreshortening moves the interaction centre of hydrogen atoms off the nucleus, along the bond towards the heavy atom to which the hydrogen is bonded.The FIT force field did not foreshorten hydrogen atom positions during parameterisation, but we find that foreshortening nevertheless improves its performance.We have therefore foreshortened the hydrogen positions of the DFT geometries by 0.1 Å during lattice energy calculations.The distributed multipoles and point charges were determined at the foreshortened H positions.The X23 benchmark consists of 23 crystals of small organic molecules.4][55] The X23 hexamine structure corresponds to the experimental structure HXM-TAM09 56 and the succinic acid structure is the monoclinic β polymorph, best represented by SUCACB02. 57For compatibility with the lattice energy minimiser DMACRYS, the symmetries of crystal structures with Z < 1 were modified to obtain whole molecules in the asymmetric unit.

Lattice energy calculations
The lattice energy consists of the intermolecular cohesive energy and an intramolecular relaxation energy, i.e. the difference in energy of the in-crystal molecular geometry and electronic structure relative to the gas phase ground state geometry.
Care was taken to use the correct lowest energy gas phase molecular conformers of oxalic acid and succinic acid.Cytosine has several tautomers and we have used the hydroxy(-amino) tautomer, in its 2b rotamer, in gas phase calculations, see Fig 1, which has been shown to be the lowest energy gas phase form. 58Molecular geometries and energies were calculated at the B3LYP/6-311G(d,p) level of theory using GAUSSIAN 09 59 revision D.01.For PCM calculations (W99rev6311P3 and W99rev6311P5), the relaxation energy contains the energy due to the geometric distortion of the molecule in the crystal and the electronic relaxation from the polarised (PCM) to unpolarised (in vacuo) electron density.
The intermolecular cohesive energy between any two molecules N and M is calculated as a sum over atom-atom pair interactions: where i and k are atoms of type ι and κ belonging to molecules M and N, respectively, separated by the distance R ik .The first two terms model the short range repulsive and attractive (dispersion) non-electrostatic intermolecular interactions, with force field parameters A ικ , B ικ and C ικ .Atomic partial charges were fitted to the molecular electrostatic potential using the program MULFIT 2.1. 60Distributed multipoles 61 were derived from the molecular charge density using GDMA 2.2.06. 62he crystal structures (except carbon dioxide) were first lattice energy-minimised using the program CRYSTALOPTIMIZER 2.4.2, 63,64 using the W99rev6311 force field.CRYSTALOPTIMIZER minimises the sum of intra-and intermolecular energies with respect to crystal packing degrees of freedom and a set of selected intramolecular degrees of freedom (dihedrals and bond angles), chosen automatically as described previously. 14The optimised structures were then re-optimised with the different intermolecular force fields, using DMACRYS 2.0.4, 46keeping the molecular geometries rigid.
For the W99rev6311P3 and W99rev6311P5 force fields, we recalculated the multipoles from a charge density obtained in an external-iteration polarisable continuum model with isotropic relative permittivity of 3.0 and 5.0, respectively.
Carbon dioxide was treated specially.A molecular geometry optimisation was performed in GAUSSIAN 09 and partial charges were calculated with the CHelpG routine. 65Multipoles and a molecular local axis system were constructed by manually adding non-interacting dummy atoms.
Charge-charge, charge-dipole and dipole-dipole interactions were calculated using Ewald summation, while repulsiondispersion interactions and all higher multipole-multipole interactions were calculated between whole molecules to a center of mass cutoff distance of 20 Å.

Lattice dynamics.
At high temperatures, thermal expansion significantly affects the crystal structure and for high-accuracy calculations of crystal properties, this cannot be neglected.We have implemented lattice vibrational free energy calculations in the quasi-harmonic approximation, so that the zero-point energy, vibrational frequencies and thermal pressures can be calculated using these force fields.This allows an explicit modelling of thermal effects.
All lattice dynamics calculations were performed in the rigidmolecule approximation, as implemented in DMACRYS, 46,66,67 with the FIT and W99rev6311P5 force fields.][70] Our method for sampling phonon frequencies and calculating the resulting thermodynamic properties have also been described in our previous work on polymorph energy differences. 14However, recent improvements in the treatment of phonon dispersion will be described here in some detail.
A great challenge in lattice dynamics calculations lies in the convergence of the phonon density of states g(ω).The anisotropic phonon dispersion in molecular crystals necessitate that several k-points are included in the calculations, and preferably in some configuration that efficiently samples the first Brillouin zone. 71,72e sample k-points by forming several linear supercells, each elongated along one lattice vector.The supercell expansion in each direction is chosen such that the distance between the kpoints is less than some target distance in reciprocal space, which we choose as 0.12 Å −1 in this work.This results in between 17 and 34 unique k-points being sampled for the crystals in the benchmark.However, sampling k-points close to the Γ-point is necessary to converge the lattice dynamical entropy and requires large supercells to yield long wavelength phonons at small k.Long supercells can be split into several shorter supercells with mutually co-prime expansion coefficients, which improves computational cost scaling, but does not improve the sampling near the Γ-point.To this end, we have implemented approximate methods to further improve the convergence, which we describe below.Results of convergence tests are included in the supplementary data †.

Debye approximation.
The Debye approximation interpolates the phonon dispersion between the Brillouin zone centre and the nearest explicitly sampled k-point. 73Debye frequency ω D is calculated from the elastic stiffness tensor C ik jl obtained from DMACRYS.66 From the elastic tensor we then calculate 3 × 3 Christoffel matrices Γ Γ Γ for 13 direction vectors km in the Brillouin zone (see ESI †).
From the eigenvalues of the Christoffel matrices, we obtain the phase velocity of sound ν in direction km .
The crystal's density is denoted ρ and the identity matrix I.We average over the three eigenvalues and obtain a mean velocity of sound ν( km ).Assuming a sinusoidal phonon dispersion around Γ, we calculate the Debye frequency ω D ( km ) as Averaging over the 13 km we obtain the Debye frequency ω D , which is an approximation of the average phonon frequency on an ellipsoid around the Γ-point.
Adding the Debye approximation to the long wavelength acoustic phonon contributions (see Eq. 8) improves the convergence of calculated thermodynamic properties with respect to the number of sampled k-points, see Fig S1 -S4.

Kernel density approximation.
Lattice dynamics calculations can only ever produce a list of discrete phonon frequencies, not the quasi-continuous density of states g(ω).However, a kernel density estimate 74,75 of the phonon spectrum rapidly converges to a distribution similar to the true density of states, 76 see Fig S7 -S10.Hence, we replace each discrete phonon frequency with a narrow Gaussian distribution and the phonon density of states is approximated with the kernel density KDE(ω).
The sum is over the six rigid-molecule vibrational modes of each molecule at each k-point, for a total of n = 6ZN k − 3 non-zero phonons for a crystal with Z molecules per unit cell, each with 6 degrees of freedom.N k is the number of sampled k-points.The kernel density is by definition normalized to unity.The choice of the kernel bandwidth h is somewhat arbitrary, but by testing we have found that a suitable bandwidth can be chosen as a fraction of the standard deviation of the phonon frequencies, h = std(ω i )/20, for the k-point sampling used here.The calculated vibrational energy varies negligibly with respect to the denominator in the range 10-100, see Fig S5 and S6.The kernel density estimate is not intended to smooth the density of states; we use it to model a small phonon dispersion around each explicitly sampled k-point.
The Helmholtz free energy is The vibrational contribution to the free energy F vib (T ) for one unit cell with rigid molecules, each with 6 degrees of freedom can then be calculated from the Debye frequency and the density of states g(ω) as: 73 where D(x) is the Debye function The kernel density and the integrals are conveniently evaluated numerically using SciPy.Naturally, carbon dioxide has only 5 intermolecular degrees of freedom, and was treated separately.

Thermal pressures.
At elevated temperatures, anharmonic vibrations lead to a thermal expansion of the crystal unit cell.To model this, we employ a calculated thermal pressure 33,77 for the temperature of interest as a finite difference in F vib (T ) between unit cells with slightly different volumes.
We believe that an isotropic scaling of the unit cell for the thermal pressure calculation could lead to errors, since molecular crystals can have strongly anisotropic thermal expansion.One option is to change the unit cell dimensions by directly scaling the lattice vectors in proportion to the elastic compliance constant in each direction.However, we found that the most consistent approach is to instead expand the unit cells by geometry optimising the crystal structure at a pressure different from the ambient pressure P 0 .We have chosen −300 MPa, taking a negative pressure as this should result in a crystal structure that is as similar to a thermally expanded structure as possible.The negative pressure causes a volume expansion between 2.4 and 7.7% in the benchmark crystal structures.
Provided that the same k-points are sampled in both the unexpanded and expanded structures, the vibrational energy varies linearly for volume changes up to about 15%, so the thermal pressure can be calculated as a finite difference between only two unit cells.The thermal pressure is then 77 Again geometry optimising to minimise E latt +PV at P = P 0 − P th results in a crystal structure close to the Gibbs free energy minimum at ambient pressure.The Gibbs free energy can then be calculated as 3 Results and Discussion

Lattice energies
In Figure 2 we show a comparison of the performance of the different force fields in reproducing the X23 benchmark lattice energies.Table 1 shows the mean absolute deviation (MAD), mean absolute relative deviation (MA%D), the systematic mean error (ME) and the random error as one standard deviation (SD).
Energy The force fields all systematically underbind the crystals.This is expected since the force fields were parameterised to thermally expanded experimental structures and sublimation enthalpies.Expressed as a percentage mean signed error, we have −7.9% for the FIT and −23.3% for the W99_ESP force field and the others in between.
The multipole electrostatics of W99_DMA improves on the point charge model in W99_ESP significantly, reducing errors by about a third.As was observed when re-parameterising W99, 45 W99rev6311 does not significantly improve the accuracy of lattice energies compared to W99_DMA; the improvement was largely seen in how well the geometries of the structures were reproduced.
The revised force fields utilizing polarised multipoles perform better than the non-polarised version since the polarisation generally increases the strength of directional electrostatic interactions, principally hydrogen bonds, and this reduces the systematic underbinding.It is clear that the PCM method improves both the systematic and random errors.The use of a dielectric constant in the PCM calculation of 5 rather than the value 3 applied during parameterisation improves the energies further and is still in the range of reasonable dielectric constants for organic materials, so can be physically justified.
The best models in terms of lattice energy are the revised force fields with multipoles derived from a PCM model with a dielectric constant of 5 (W99rev6311P5) and the FIT force field.The latter is the most accurate of the tested force fields, and it performs remarkably well both with respect to systematic and random errors.
The accuracy of the force field results depends partly on the level of theory used for the intramolecular energies.However, apart from oxalic acid, where the crystalline geometry is calculated to be between 14.7 and 16.4 kJ mol −1 above the gas phase minimum, the molecular relaxation energies are small (less than 3.2 kJ mol −1 in the unpolarised models and up to 6.8 kJ mol −1 with PCM, ε r = 5.0).Given these magnitudes, we expect that errors in the intramolecular energies are a small contribution to overall errors.
In Table 2 and Figure 3 we compare the two best force fields with several DFT-D methods previously benchmarked against the X23 or C21 sets.For published C21 results, we have recalculated the errors relative to the X23 benchmark lattice energies and we have performed PBE-XDM and PBE-D2 calculations on hexamine and succinic acid using QUANTUMESPRESSO. 78ergy model MA%D MAD ME SD TPSS-D3 35 5.21 The FIT and W99rev6311P5 force fields reproduce lattice energies better than PBE-TS and B3LYP (without dispersion correction).The random error is of the same magnitude for most of the benchmarked methods.
Most solid state calculations on molecular crystals are performed using GGA functionals, such as PBE.This is due to the cost of evaluating exact exchange, which is necessary for hybrid functionals, using plane-wave basis sets.Only with the most accurate dispersion correction methods: exchange-hole dipole (XDM), many-body dispersion (MBD) and Grimme's D3 method, is the PBE functional significantly better than the force fields we have tested here.We note that the force field-like dispersion correction used in these DFT methods has a more complex form than what is used in the current force fields: XDM uses a sum of C 6 , C 8 and C 10 dispersion terms, while MBD includes non-pairwise additive dispersion contributions.It seems reasonable to expect that the future development of force fields including such terms could reduce the gap between the current force fields and the best DFT-D methods.
The hybrid functionals HSE06 and PBE0 and the meta-GGA functional TPSS together with D3 or many-body dispersion methods are also significantly better than the force fields. 35However, even using the best DFT-D methods, the errors are only about 2 or 3 times smaller than the force field methods.

Lattice parameters and geometries
In Table 3 we compare the performance of the force fields in reproducing experimental lattice parameters, without any thermal adjustments.Symmetry-independent lattice vector lengths for 22 of the crystal structures are used, excluding urea.With the force field models, we find that the experimental tetragonal P42 1 m urea structure (UREAXX02) is a vibrationally averaged structure corresponding to a saddle point on the potential energy surface. 66Breaking the symmetry results in a stable orthogonal P2 1 2 1 2 structure.This prevents a direct comparison of lattice parameters for this structure.The FIT force field performs very well in reproducing the experimental lattice parameters.It is comparable in accuracy to PBE-XDM.The B86b-XDM functional, PBE-TS and PBE-D2 are better than the force field methods in reproducing lattice parameters.The W99rev6311P5 force field also reproduces the crystal structures well, with the exception of benzene and naphthalene.The W99-based force fields do not perform well for crystals with aromatic ring systems in T-stacked configuration, 31 see Table S3.However, systems in π-stacked configurations are well reproduced. 28The point charge force field W99_ESP performs significantly worse than the other methods and the errors in lattice parameters are about twice as large as the multipole-based force fields.

Energy model MAD [Å] MA%D
The errors are not systematic; the mean errors in lattice vector lengths are merely −0.03 Å and +0.04 Å for the FIT and W99rev6311P5 force fields respectively.PBE-TS reproduces crystal geometries very well, 27,33 but not lattice energies. 34Of course, an energy model must perform well on both energies and geometries to be useful and a benchmark study really should consider both.Unfortunately, this is not always done.
A common way to specify crystal structure similarity is the root mean square deviation in atomic coordinates in a cluster of molecules taken from a crystal structure. 9,79Taking 20-molecule cluster comparisons (RMSD 20 ), the FIT and W99rev6311P5 force fields reproduce the experimental crystal structures with an average RMSD 20 of 0.234 Å and 0.258 Å, respectively, see Table S3.

Vibrational free energies
An accurate method should reproduce the shape of the lattice energy surface well enough to give accurate vibrational energy contributions and thermal expansion.Therefore, for benchmarking purposes, we argue that these contributions should be calculated with the method being benchmarked.
Lattice dynamics calculations involve a much higher computational cost than lattice energy minimisation.We and others have previously shown that it is crucially important to adequately sample phonons across the Brillouin zone to achieve reliable free energy contributions. 14,20Given the high computational expense of some of the methods currently being developed for lattice energy calculations, methods for accelerating the convergence of vibrational energies with respect to reciprocal space sampling could enable otherwise impractical free energy calculations.Furthermore, even using computationally efficient force field methods, free energy calculations can become the computational bottleneck in large studies of polymorphism or crystal structure prediction, where many hundreds of structures sometimes need to be evaluated.
Figure 4 shows the convergence of the calculated vibrational energy, F vib , with respect to the target distance between sampled k-points for two of the benchmark crystal structures, benzene and acetic acid (results for urea and cyanamide are shown in the supplementary information †).For the purposes of testing, we treat calculations with a maximum k-point spacing of 0.1 Å −1 as converged.This corresponds to between 25 and 37 unique k-points for these four crystal structures.
The magnitude of F vib increases (F vib becomes more negative) as the k-point sampling density is increased, which is mainly due to including the contributions of the low frequency acoustic vibrations.For this reason, adding a Debye approximation of the acoustic modes between the Brillouin zone centre and the nearest sampled k-point improves the convergence significantly.
The KDE approach to approximating phonon dispersion around each sampled k-point also reduces errors in F vib (Fig. 4), but by less than the Debye correction.While the energetic correction when using the KDE is sometimes small, its use has the effect of smoothing the convergence of F vib by reducing the quantization effects resulting from discrete sampling of k-points (e.g.near 0.8 Å −1 in the benzene results and 1.0 Å −1 for acetic acid).
Used together, the Debye and KDE approximations improve the F vib calculations dramatically, allowing accurate calculations of the vibrational energy contributions at much lower computational cost; in the cases studied here, sampling only a few k-points explicitly is sufficient.

Thermal pressures
Calculating the thermal expansion when using empirically fitted force fields could double-count the expansion, as some thermal effects have already been absorbed by the parameters of the force field.However, since modelling of thermal expansion can be important, we compare computationally thermally expanded crystal structures to the benchmark data, to see how large the errors are from this double-counting.
In this work, we have used the computational efficient approach of modelling thermal expansion by including a thermal pressure during lattice energy minimization.Thermal pressures using FIT and W99rev6311P5 are reported in Table 4.Our results do not agree with the values reported by Otero-de-la-Roza and Johnson 33 .The discrepancy is most likely due to the different k-point sampling, the unit cell expansion method, and our methods for improving convergence of the vibrational contributions to the free energy.We find that including phonon frequencies from a single k-point leads to inaccurate thermodynamic properties. 14

Structure
Temp  For the FIT force field, the deviations in lattice vector lengths for the thermally expanded structures are MAD = 0.188 Å and MA%D = 2.82%.The W99rev6311P5 force field yields MAD = 0.264 Å and MA%D = 3.89%.After thermal expansion, experimental crystal structure geometries are reproduced with an average RMSD 20 of 0.234 Å and 0.314 Å for the two force fields respectively.
Since the force fields were parameterised to reproduce experimental crystal structures, the explicit modelling of thermal expansion here constitutes a double counting of some of the thermal expansion, although low temperature (< 100 K) crystal structures were deliberately chosen in re-parameterising W99 to minimise this effect.This should cause a systematic positive error in the lattice vector lengths.However, the mean errors are only +0.14 Å and +0.05 Å for the W99rev6311P5 and the FIT force fields respectively.Again the FIT force field performs very well and does not have a large systematic overestimation of unit cell These results confirm the soundness of the computationally very efficient thermal pressure method and that the error caused by the double-counting of implicit and explicit thermal expansion is negligible.

Conclusions
The FIT and W99rev6311P5 force fields perform very well for molecular crystals.Both have a relatively large systematic error in that they underestimate absolute intermolecular energies by 8 and 15% respectively, while the mean absolute deviations and random errors are comparable to several popular GGA DFT-D methods.The performance of force fields with atomic multipole electrostatics greatly reduces errors compared to an atomic partial charge model.
For the application of crystal structure prediction, or other computational studies of polymorphs, the systematic error is less of a concern, as it is the stability difference between polymorphs that must be reproduced accurately.The force fields described here have random errors comparable to PBE-XDM, PBE-D2 and PBE-TS.In fact, only hybrid or meta-GGA functionals, or PBE with sophisticated (many-body or XDM) dispersion corrections outperform the force fields in reproducing experimental lattice energies and geometries.Even then, the errors are only 2 or 3 times smaller.These findings could help guide the development of better force fields, where developments in the functional form have thus-far focused largely on accurate electrostatics.More elaborate forms of dispersion in intermolecular force fields should be developed along the of the successful DFT dispersion correction schemes.
We have implemented methods by which it is possible to very rapidly model the thermal expansion of molecular crystals.This allows us to estimate Gibbs lattice vibrational free energies of crystals of small rigid molecules in a matter of a few CPU-hours.The accuracy of vibrational energy contributions is improved by including a Debye correction to the phonon density of states, along with a kernel density estimate of phonon dispersion around each sampled k-point.The error caused by the double-counting of some thermal contributions in empirical force fields is very small, and accurate modelling of crystals at finite temperatures can be performed with these force fields.
For future benchmarks, we encourage the consideration of both lattice energies and crystal geometries, since generally useful methods must give good accuracy in both.This must be done with care: thermal effects are large enough to significantly affect benchmarking results.In particular, 0 K calculated geometries cannot be directly compared to experimental crystal structures.Thermal pressures can be used to correct for this, but the thermal pressure (or any other method chosen to model thermal expansion) should be calculated with the method being benchmarked and phonon dispersion must be accounted for sufficiently well for accurate vibrational energies.We have described efficient methods for Brillouin zone integration that makes the calculated thermodynamic properties virtually independent of the exact choice of sampled k-points.
Since the impressive crystal structure prediction results by Neumann et al. 11 in the fourth and fifth crystal structure prediction blind tests 9,12 there has been a major focus on DFT-D methods in the computational crystallography community.However, other DFT-D methods have performed less well and there is little correlation between computational expense and predictive ability in currently used methods. 13One should not assume that density functional methods are always more accurate than intermolecular force fields, so long as the force field is constructed with a respect for the necessary functional form to describe the required physical interactions.When it comes to large-scale studies of organic molecular crystals, particularly crystal structure prediction, anisotropic atomic multipole based force fields are competitive with the best alternatives and should continue to be developed.

Fig.
Fig. Relative errors in lattice energies for the X23 set of crystal structures calculated with six different intermolecular force fields.Negative values represent underbinding of the crystal.

Fig. 3
Fig.3Relative error in lattice energies calculated with FIT, W99rev6311P5 and several PBE-D methods.PBE-D results reproduced from Reilly and Tkatchenko34 and Otero-de-la-Roza & Johnson.33 Negative values represent underbinding.

Fig. 4
Fig.4 Convergence of F vib with respect to the maximum distance between sampled k-points, for (top) benzene and (bottom) acetic acid.Neat refers to the vibrational free energy calculated with neither the Debye nor KDE corrections for phonon dispersion.The leftmost data points correspond to k = 0 phonons only.
* Anthracene and naphthalene not included.

Table 3
Mean absolute deviation (MAD) and mean absolute relative deviation (MA%D) in lattice vector length for the different force fields and a selection of DFT-D methods.The urea structure was not included in the force field calculations.
* Hexamine and succinic acid not included.

Table 4
Thermal pressures in MPa calculated with the FIT and W99rev6311P5 force fields at the experimental temperatures.
* Supercells of carbon dioxide were unstable during thermal expansion with FIT.