A non-empirical intermolecular force-field for trinitrobenzene and its application in crystal structure prediction

An accurate potential energy surface (PES) for the intermolecular interaction between organic molecules is a key requirement for simulating the properties of the gas phase and the liquid, the solid, and other condensed phases including crystalline polymorphs and the transitions between them.¹ The methods of realistically evaluating many condensed phase properties require an analytical model so that huge assemblies of molecules can be simulated over long time periods. Such analytical models are generally pairwise intermolecular potentials; however, non-additive many-body contributions also need to be considered when studying condensed phases. Furthermore, potentials for flexible molecules may depend on molecular conformation. Many molecules change conformation when packed with other molecules, which can result in changes in the molecular charge distribution that affect the intermolecular interactions. Hence, while model intermolecular potentials that are able to describe the properties in all phases of the inert gases, such as argon, have been available for decades² and considerable progress has been made on a highly accurate pair potential for water (see Refs. 3 and 4 and references therein), most atomistic modeling of condensed phases of organic molecules has used empirical force-fields. Force-fields have been heavily developed for biomolecular simulations, but the fitting to experimental data inevitably absorbs system-specific errors and so limits accuracy if they are used for molecules and properties outside the scope of fitting and validation.⁵ Hence, advancing the computer simulation work requires a transferable approach to deriving non-empirical potentials, based on the theory of intermolecular forces, and building on the electronic level modeling of intermolecular interactions.

Since electron correlation is the cause of the dispersion, which is important and sometimes dominant to the long-range attractive force between organic molecules, even calculating dimer interaction energies requires high-level electronic structure treatments.^6,7 As the charge densities of closed shell organic molecules are only slightly perturbed by intermolecular interactions, their intermolecular forces are well described by perturbation methods. Indeed, perturbational methods such as symmetry-adapted perturbation theory (SAPT),^4,8,9 and its variant based on density-functional theory [SAPT(DFT)],^6,7 are ideally suited for this problem as they are not only accurate, and also, SAPT(DFT) is computationally efficient. Additionally, perturbation theory provides a rigorous match between the short-range energies and the long-range analytic expansions using multipoles, polarizabilities, and dispersion coefficients. The size of organic molecules means that a density-functional charge distribution of each molecule is the most practical basis for developing a model intermolecular potential using SAPT(DFT)^6,7 to provide the essential link between electronic and atomistic level modeling of the intermolecular forces.

In SAPT(DFT),^10,11 the intermolecular interaction energy between a pair of rigid molecules is defined through the sum of energy components,

Here, $Eelst(1)$ and $Eexch(1)$ are the first order electrostatic and exchange–repulsion energies, respectively. $EIND(2)$ is the second-order induction energy, which includes the exchange-induction energy, and $EDISP(2)$ is the total second-order dispersion energy, which includes the exchange-dispersion energy,

The term $δintHF$ estimates the main error in truncating the expansion to second-order¹² and is dominated by the higher order induction energy contributions. For condensed phases, the usual summation over the intermolecular pair potential between all pairs of molecules in the system is an approximation (the pairwise additive approximation). The induction energy includes contributions from quantum charge-delocalization and classical polarization effects.¹³ While both terms are likely to be non-additive, it is the classical polarization contribution that is dominant in the condensed phase as it is long-ranged. This term is difficult to evaluate and include in simulations as it requires numerical iteration to self-consistently account for the mutual polarization of all components in the lattice.¹⁴ 

SAPT(DFT) and related supermolecular partitioning methods have been the basis of various methods of deriving non-empirical intermolecular potentials (see Refs. 15 and 16). Non-empirical potentials are particularly valuable for modeling energetic materials, where the key functional groups, such as NO₂, have rather different charge distributions from the N and O atoms in well-parameterized biomolecular force-fields. The difficulties in handling explosives naturally leads to a lack of experimental data for empirical parameterization, and this has led to a desire for reliable computational prediction of condensed phase properties. Pioneering studies of the organic solid-state using analytical fits to SAPT(DFT) interaction energy calculations have concentrated on small energetic molecules,¹⁷ such as RDX^18–20 and FOX-7²¹ where predictive modeling requires reliable extrapolation into the repulsive region that is not sampled in fitting empirical potentials. These force-fields are in an isotropic atom–atom form for ease of use in molecular dynamics simulation programs.

This paper seeks a non-empirical model potential for use in Crystal Structure Prediction (CSP) studies of energetic molecules. The vision behind one of the earliest CSP methods was to predict whether an energetic molecule would crystallize in a sufficiently dense crystal to be a promising material prior to synthesis,²² although CSP studies can also be used to help solve the crystal structures.^23,24 The anisotropy of the atomic charge density arising from features such as lone pairs and π electron densities has proved to be important in realistic CSP studies,^25,26 with one of the leading family of methodologies being based on using a distributed multipole representation of the molecular charge density.²⁷ While most CSP studies have used empirical isotropic exp-6 repulsion-dispersion potentials in conjunction with an anisotropic atom–atom electrostatic model, an early non-empirical potential that additionally included anisotropy in the repulsion²⁸ for C₆Br₂FClH₂ was successful in an international blind test.²⁹ This potential completely neglected the induction contribution to the lattice energy.²⁸ A model for the use of atomic dipolar polarizabilities has since been implemented in the organic solid state modeling code DMACRYS,¹⁴ although the numerical iteration required means that the polarization energy cannot be included in structure optimization. A fully anisotropic non-empirical atom–atom model intermolecular potential for pyridine has been derived using the program CamCASP³⁰ and used within DMACRYS in CSP studies to study the polymorphs of pyridine.¹⁶ The non-empirical potential was successful in identifying the structure of a high-pressure polymorph³¹ unlike the empirical repulsion–dispersion potential FIT¹⁴ that is often used in CSP studies.¹⁶ Thus, this paper seeks to extend this approach to predicting the crystal structures of energetic molecules. The use of an atomistic intermolecular potential is far more computationally efficient than using periodic dispersion-corrected density functional (DFT-D) calculations on hundreds of low energy crystal structures, as it is being increasingly used in CSP studies.²⁷ Thus, further examples of evaluating the use of non-empirical intermolecular pair potentials in CSP are urgently needed.

The energetic molecule chosen for this study was the high-explosive trinitrobenzene (TNB), which is extremely volatile in its dehydrated powder form.³² It has a melting point of ∼400 K,^32,33 and although it was first discovered in 1883,³⁴ a crystal structure was not determined until 1972.³⁵ The polymorphism of TNB (Fig. 1) was first reported in 2004 after the addition of trisindane in an attempted cocrystallisation³⁶ that resulted in the serendipitous crystallization of forms II (Z′ = 2) and III (Z′ = 1). The experimental data on forms II and III are limited to the melting points and crystal structures, with form III having the highest melting point. This, and a force-field lattice energy estimate, suggests³⁶ that form III is more thermodynamically stable than kinetically favored form I. Form III is estimated to be more stable than form I by 0.6 kJ mol⁻¹ by periodic PBE-D3(BJ) lattice energy calculations.³⁷ 

The crystal structures of TNB (Fig. 1) show that the crystalline conformations do not differ much in the NO₂ torsion angles, and consequently, a rigid-molecule intermolecular potential is worth developing as a first approximation. TNB is highly symmetrical in its gas-phase optimized geometry so that there are only five distinct atomic types. This type of rigid-molecule approximation was often reasonably successful for substituted aromatic molecules in early CSP studies.^38,39 Incorporating flexibility into the force-field so that it accurately balances the inter and intramolecular forces is highly challenging, as shown by a contemporary study of α-RDX²⁰ as well as CSP studies on flexible pharmaceuticals.^40,41 Indeed, larger changes in the NO₂ torsion angles would affect the charge distribution of the nitro group to an extent that would require reparameterization of the model potential for accurate lattice energies.⁴² 

We can develop non-empirical derived intermolecular force-fields (DIFFs) for TNB with an atomic multipolar expansion of long-range components V_LR and an exponential model for the short-range terms V_SR. The functional form¹ for the interaction between molecules M and N containing atoms i and k of types ι and κ (types are C_H, C_N, N, O, or H) at a separation R_ik and orientation Ω_ik can be written as

The electrostatic and induction terms are implicitly summed over the permanent atomic multipoles $Qti$ (where t = lk; l is the total angular momentum, k is the z-component of angular momentum, and l ≤ 4), induced dipoles $ΔQti$ and the interaction tensors $Ttuik$ are defined using the axes (Fig. 1) on the two molecules.¹ The induced moments and the isotropic dispersion coefficients, $C6ικ$⁠, $C8ικ$⁠, and $C10ικ$⁠, are calculated from atomic polarizabilities that are derived from the molecular charge distribution using linear-response Kohn–Sham theory together with suitable distribution techniques. The dispersion and polarization expansions are damped using the Tang–Toennies⁴³ damping function f_n. All the short-range terms are modeled by an exponential term where the atomic repulsion anisotropy ρ(Ω_ik) is derived specifically for the molecule by fitting to SAPT(DFT) calculations, with guidance from the distributed density overlap model.¹⁶ The partitioning of the interaction energy and molecular properties into atomic contributions is based on the basis-space version⁴⁴ of the iterated stockholder atoms approach of Lillestolen and Wheatley.^45,46

The DIFF force-field functional form is complex. However, the complication of the anisotropy is likely to be important for many molecules as there is a significant difference in the reported quality of fits to SAPT(DFT) dimer intermolecular potential energy surfaces from using isotropic or anisotropic models. One automated generation of isotropic model potentials reported a typical error of about 0.8 kJ mol⁻¹ in the negative energy region,⁴⁷ whereas a different approach that included anisotropy exhibits average errors of 0.2 kJ mol⁻¹ in the attractive region for a different wide range of interacting dimers.⁴⁸ Given that lattice energy differences between polymorphs are often less than 2 kJ mol⁻¹, with about 80% calculated to be less than 4 kJ mol⁻¹,⁴⁹ the modeling of the organic solid state requires the best possible accuracy. The solid state also samples the potential energy surface at a much wider set of relative separations and orientations in the known and CSP generated crystal structures than those that are tested in dimer calculations. Thus, the ability to reproduce crystal structures and the relative energies of known and hypothetical structures generated in a CSP study is a severe test of the accuracy of a force-field.

In this paper, a set of DIFF potentials was developed using a combination of CamCASP,³⁰ Psi4,⁵⁰ NWChem,⁵¹ and ORIENT.⁵² These potentials have the same long-range terms but increase the cost of development, first by the inclusion of anisotropy in the short range term and then allowing the isotropic short range coefficients to be relaxed by fitting to the total interaction energy, calculated up to second-order SAPT(DFT) interaction energies.^30,52 The size of TNB meant that we needed to make some changes in the DIFF potential development method compared to that used for pyridine¹⁶ and water;³¹ however, the philosophy of the approach remains the same. The various DIFF model potentials were tested in a CSP (Z′ = 1) study to assess whether form III was found and also to study the relative stability of all observed polymorphs. The aim is to advance our ability to model the condensed phases of energetic molecules such as TNB by developing a workflow for deriving a sufficiently realistic non-empirical potential and implementing it in a CSP study.

The non-empirical potential derivation begins from the isolated molecular structure of TNB obtained by optimizing the isolated molecular conformation using the Gaussian09 program,⁵³ with the PBE0 functional^54–56 and the d-aug-cc-pVTZ Dunning basis.⁵⁷ This molecular conformation was held rigid for the development of the force-fields and in all simulations, except those used to assess this approximation. The three non-empirical DIFF potentials were constructed using the workflow in Fig. 2. These had the same long-range terms but differed in whether the short-range terms were isotropic, DIFF-srL0, or anisotropic with terms up to L = 2, DIFF-srL2. Additionally, the model parameters determined using first-order SAPT(DFT) energies were either allowed to relax against the total SAPT(DFT) intermolecular energies, indicated by “(rel),” or kept at their un-relaxed values, indicated by “(norel).” The comparisons were made with the non-empirical isotropic exp-6 model potential (FIT)¹⁴ whose parameters were fitted to organic crystal structures and have widely been used in CSP studies. The FIT model was either combined with multipoles computed using the Gaussian-Distributed Multipole Analysis (GDMA)⁵⁸ method with the PBE0/6-31G(d,p) charge density, denoted as FIT(GDMA,L4), or the same ISA distributed multipoles used in the DIFF potentials [FIT(ISA,L4)]. Here, “L4” denotes that these multipoles contained terms of maximum rank of 4 on all sites.

The molecular properties of TNB were computed using the modified basis-space iterated stockholder atoms (BS-ISA)^44,59 and the ISA-Pol⁵⁹ algorithms implemented in CamCASP.³⁰ The Kohn−Sham orbitals and orbital energies needed for these calculations were calculated using the NWChem. 6.6⁶⁰ code. We used the CS00 (Casida-Salahub) asymptotically corrected⁶¹ PBE0 functional and the augA-Sadlej basis set,⁶² a Sadlej-pVTZ basis that has been augmented with additional diffuse functions and functions of f-symmetry. This basis set was chosen as a cost-efficient yet accurate alternative to the larger d-aug-cc-VTZ basis set (see Fig. 1 of the supplementary material). The asymptotic shift (Δ) for TNB was calculated using the experimental ionization potential of 0.4028 a.u.⁶³ and the energy of the highest occupied molecular orbital, also computed with PBE0/augA-Sadlej, to give Δ = 0.0539 a.u. The linear-response density functional calculations were performed using the hybrid ALDA+CHF kernel in CamCASP 6.0.³⁰ 

The SAPT(DFT) dimer interaction energy calculations used to determine the induction energies and regularized induction energies for the polarization models were computed using CamCASP 6.0 with the above numerical details except that we used the aug-cc-pVTZ basis in the MC+BS format using a 3s2p1d1f midbond basis from the CamCASP basis library. Note that in CamCASP 6.0, the second-order exchange-induction energy is computed without the single-exchange (S²) approximation using a closed-shell implementation of the general energy expression derived by Schäeffer and Jansen.⁶⁴ The first-order SAPT(DFT) interaction energies were computed using the above basis but with CamCASP 5.9. However, for reasons of computational efficiency, the interaction energies computed up to second-order were performed using the CamCASP with DFT orbitals and orbital energies calculated with the Psi4 1.1 code. In this case, the asymptotic correction was changed to be the GRAC correction⁶⁵ as Psi4 does not support the CS00 scheme. The two asymptotic correction schemes use the same long-range form but differ in the manner in which the long-range functional is spliced with the PBE0 functional. Consequently, some differences are to be expected, but these are likely smaller than the overall errors of the approach

The molecular charge density was partitioned into atomic charge densities to obtain the permanent distributed multipole moments $Qti$ of rank l ≤ 4 using the most recent iterated stockholder atoms algorithm⁴⁴ implemented in using CamCASP 6.0, the ISA-A algorithm.⁵⁹ The ISA-A method used the augA-Sadlej main basis with PBE0/CS00. The additional basis sets needed for the ISA calculation included the aug-cc-pVTZ⁶⁶ basis for the density-fitting and the aug-cc-pVQZ-RI basis modified with the ISA Set2 s-functions from the CamCASP basis library. The permanent atomic multipoles $Qti$ of rank up to hexadecapole were obtained for each atom.⁵⁹ The ISA-Pol⁵⁹ algorithm as implemented in CamCASP was used to derive the anisotropic dipolar (L1) polarizabilities $αtt′ik$ for use in modeling the polarization energy and the isotropic L3 atom–atom dispersion model,⁶⁷ i.e., the $C6ικ,C8ικ, and C10ικ$⁠.^14,67 The long-range multipolar properties used in the DIFF potentials are given in the supplementary material.

The parameters in the short-range potential include the terms in the exponential part of the Born–Mayer potential and also the damping parameters for the dispersion and polarization models. The latter were determined by direct fitting to selected SAPT(DFT) second-order energies. However, for the Born–Mayer repulsive terms, we adopted the multi-stage fitting technique, which has been extensively described by Misquitta and Stone.¹⁶ In this approach, the initial parameter estimates are obtained via the distributed density-overlap parameterized against a dense set of first-order SAPT(DFT) energies. This provides the potential DIFF-srL2(norel) and a model with only isotropic repulsion DIFF-srL0(norel).

The first-order exchange, $Eexch(1)$⁠, and electrostatic, $Eelst(1)$⁠, interaction energies were computed using a combination of NWChem. and CamCASP¹⁶ for 2000 pseudo-random dimer configurations generated using the quasi-random Sobol sequence and Shoemake’s uniform random rotations algorithm,¹⁶ implemented in CamCASP and used previously for molecules as non-spherical as pyridine and C₆Br₂ClFH_2.²⁸ The distributed electrostatic energy, $Velst(1)$⁠, was then calculated for the 2000 dimer configurations using the multipoles $Qti$ and used to calculate the penetration energy (⁠$Epen(1)=Eelst(1)−Velst(1)$⁠) for each dimer configuration.

The fitting of the first-order short range energies to the sum of atom–atom exponential terms in Eq. (3), $Eexch(1)+Epen(1)∼∑i∈M,k∈NGexp−αικRik−ρικΩik$⁠, allowing for all the symmetry allowed shape-anisotropy terms of rank 1 or 2 (see Table VI of the supplementary material), is an ill-defined problem. Hence, the parameters are initially estimated by using the ISA atomic charge densities and assuming that $ESR(1)$ is proportional to the overlap of the charge distributions for each of the 15 different atom–atom interactions in TNB,

The isotropic model was derived by using the fits to the distributed density overlap model to provide the initial parameters. These parameters were relaxed in stages to fit the 2000 $Eexch(1)$ + $Epen(1)$ energies using constraints to keep the relaxed parameter values close to the initial overlap model values. The anisotropic model was derived in the same way, with the fitting to the distributed density overlap model being done systematically to determine which anisotropic terms made a significant improvement on the isotropic fits. Once the choice of anisotropic terms had been made, these parameter values were tightly constrained in the relaxation stage. Further details are in 2.4 of the supplementary material.

The final potential parameters required for Eq. (3) are those required for damping the multipolar induction and dispersion expansion at a short range to allow for molecular overlap using the Tang–Toennies damping function,⁴³ 

Note that, in principle, the damping parameters are dependent on the sites involved. Indeed, this has been demonstrated for the water dimer,⁴ but due to limitations of the DMACRYS¹⁴ code, here, we are restricted to single damping parameters for the polarization and dispersion models. This results in a degradation of the quality of the resulting models and the need for compromises. A preliminary estimate of $βdisp$ = 1.8 a.u. can be derived⁵⁹ from the ionization potential of TNB. This was investigated by comparing the second-order dispersion energies $EDISP(2)=Edisp(2)+Eexch−disp(2)$ and the multipolar dispersion energies $VDISP(2)$ of DIFF-srL2(norel) of both the configuration scan and crystal structure dimer datasets (defined below) for a range of values of β_disp. The estimate β_disp = 1.8 a.u. gave a good approximation of the dispersion energy at short-range, but a value of β_disp = 1.65 a.u. gave the best fit to the data (see 2.3.2 of the supplementary material) and so was used in all the DIFFs.

Defining the polarization damping parameter is more involved as this parameter is more site-dependent than the dispersion damping and also because the polarization damping needs to be determined by fitting not to the second-order induction energy, but to the second-order polarization energy defined by regularized SAPT(DFT). We compared $EPOL(2)=EIND2Reg$ with the multipolar polarization energies $VPOL(2)$ of DIFF-srL2(norel) for various values of β_pol. It has been suggested that $βpol=βdisp$ = 1.35 a.u.,⁶⁸ but this damping is insufficient. No particular value of β_pol performed perfectly across all separations, particularly for repulsive configurations (see 2.2.2 of the supplementary material); however, β_pol = 1.0 a.u. was selected as it gave the smallest errors around the equilibrium distances of the low energy dimers and for the various two-molecule contacts in the crystal (see Figs. 3 and 4 of the supplementary material). These values of β_disp = 1.65 a.u. and β_pol = 1.0 a.u. were used in Eq. (3) to complete the parameterization of the isotropic DIFF-srL0(norel) and anisotropic DIFF-srL2(norel) model intermolecular potentials.

The SAPT(DFT) energies to second-order include a charge-delocalization energy,

Both the $ECD(2)$ and $EPOL(2)$ include second-order exchange-induction effects.¹³ $ECD(2)$ is very small for TNB except for highly repulsive configurations (⁠$Eint(2)>150$ kJ mol⁻¹) and, in most dimers, is negligible in comparison with the polarization energy (see 2.6 of the supplementary material). Contributions to the interaction energy from third and higher-orders in the interaction operator are usually included via the $δintHF$ term; however, this term is known to overestimate these effects in complexes with little or no hydrogen bonding, as is the case for TNB. Furthermore, we do include higher-order polarization effects through the iterated polarization model. After investigation (see 2.6 of the supplementary material), we have not included the $δintHF$ term in this study.

Finally, three datasets of TNB dimer intermolecular energies, $EintSAPTDFT=Eint=Eelst(1)+Eexch(1)+EPOL(2)+EDISP(2)$⁠, were calculated using the dimer configurations described below. The configuration scans and the pseudo-random dataset were used for refining the isotropic short-range parameters $ρ00ικ$ and $α00ικ$ to absorb errors from the partitioning of the energies to give DIFF-srL2(rel) (see 2.5 of the supplementary material for details). A set of dimer orientations derived from the CSP-generated crystal structures was used as an independent test of the quality of all the model intermolecular potentials for TNB.

The gas phase dimers corresponding to minima on the DIFF-srL0(norel) were found using the basin-hopping algorithm, a two-part iterative global optimization treatment that employs a global stepping algorithm with local minimizations at each step⁶⁹ using ORIENT,⁵² as a good way of locating minima on a bumpy potential energy surface. For TNB, an energy corresponding to kT at a temperature of 500 K determined the acceptance criterion for up to 100 000 steps with a maximum displacement Δx = 0.10 and rotation Δθ = 30° for each step. The CLUSTER procedure in CamCASP was used to find the unique dimer geometries orientated by their axes of inertia.

The configuration scan dataset, used for fitting the damping parameters and illustrating the potential fits, comprised up to 33 points for each dimer orientation, with the center of mass separation scaled from the minimum energy structure by a factor between 0.76 and 1.4. The set of dimers in this dataset was considered too correlated to use in the relaxation process, so a truncated dataset of only 9 points per dimer orientation (45 points in all) was used in the relaxation of the short-range parameters.

A further 1000 random configurations were added to the 2000 generated for the first order energies by restarting the quasi-random Sobol sequence. The energies of these 3000 pseudo-random dimer configurations were evaluated using DIFF-srL2(norel). $Eint(2)[SAPTDFT]$ calculations using Psi4/CamCASP were initiated for the 200 most stable and 130 least stable (most repulsive) configurations. The final pseudo-random $Eint(2)[SAPTDFT]$ dataset comprised 190 random dimers between −10 kJ mol⁻¹ and −22 kJ mol⁻¹ (the energy of the most stable pseudo-random dimer) and 92 repulsive structures with energies ranging from 3 kJ mol⁻¹ to 65 kJ mol⁻¹. The results of an initial fit to the configuration scan dataset influenced this choice of points for the pseudo-random dataset that was added to it for the fitting.

The 100 most stable CSP structures obtained after lattice energy optimization with FIT(ISA,L4) were analyzed with the Crystal Similarity tool in Mercury⁷⁰ to give the geometries of all pairs of molecules in van der Waals contact. E_int[SAPT(DFT)] energies were calculated on the 66 most common distinct dimer configurations. This dataset was used to test the quality of the intermolecular potential models.

The $ρ00ικ$ and $α00ικ$ short-range parameters of DIFF-srL2(norel) were relaxed to the E_int[SAPT(DFT)] energies of the 45 of the 160 structures found in the configuration scans of the 5 gas-phase minima and the 282 pseudo-random structures, as described in 2.5 of the supplementary material.

The TNB CSP study employed CrystalPredictor 2.0.1^71,72 to use the optimized isolated molecular structure of TNB to generate Z′ = 1 crystal structures within the 59 most probable space groups⁷³ using the FIT potential and the point charges from the ISA multipole model [FIT(ISA,L0)]. This generated 13 956 unique structures that were lattice energy minimized using DMACRYS¹⁴ with the three non-empirical potentials [DIFF-srL2(rel), DIFF-srL2(norel), and DIFF-srL0(norel)] and the empirical isotropic exp-6 potential FIT¹⁴ with the two distributed multipole electrostatic models [FIT(ISA,L4) and FIT(GDMA,L4)]. Starting the lattice energy minimizations from the Crystal Predictor/FIT(ISA,L0), output structures ensured that more of the differences in the model potential energy surfaces were sampled than the usual hierarchical approach in CSP studies of minimizing from the lattice energy minima obtained with the previous model in the hierarchy of expected accuracy.

The structure optimizations were carried out without the polarization term. A Hessian check confirmed that each structure was a minimum in the lattice energy, and the structures clustered to remove duplicates. The polarization term was subsequently included as a single point energy evaluation, iterated so that the induced moments were maintained the crystal symmetry. Corresponding calculations were performed on the experimental crystal structures after the molecular conformation had been replaced by the optimized isolated molecular structure.

The crystal structure similarity tool in Mercury 3.6⁷⁴ was used to determine the number of molecules, n, in a maximum coordination cluster (30 for crystals, 2 for dimer comparisons, and 1 for molecular conformations), which matched within a 30% distance in intermolecular atom–atom distances and 30° in interatomic intermolecular angles, and the corresponding root mean square difference, RMSD_n, excluding hydrogen atoms.

The planar optimized structure has 5 unique atoms, all of which have anisotropic charge distributions within the ISA-A partitioning, as shown in Fig. 3. The corresponding long-range potential parameters are given in Tables I–III of the supplementary material. However, the derivation of the short-range anisotropy showed that only anisotropy on O and N made a significant difference, reducing the weighted root mean squared (rms) residuals in reproducing the $Eexch(1)$ + $Epen(1)$ energies from 1.16 kJ mol⁻¹ to 0.77 kJ mol⁻¹ (see Table V of supplementary material). The required repulsion and its anisotropy can be explained by considering that these are the only atoms where the repulsive wall is sampled over a range of orientations.

The dimer structures that correspond to minima in the intermolecular pair potential are fairly insensitive to the model, with there being four structures that have the planes of the molecules roughly parallel but significantly displaced and one where there is a large angle between the molecular planes (Table I). However, the relative energies and order of the minima differ significantly with the model potential.

The relaxation of the isotropic repulsion parameters to the total second-order energies can be very dependent on the choice of the dimer configurations included if these dimers are not chosen in a balanced way. If the fitting was just done to the configuration scan dataset (with 33 dimers for each of the 5 dimer orientations), there was an improvement representing those energies, but the energies of the crystal dimer dataset became significantly worse (see 2.5 of the supplementary material). This is a reflection of the heavily correlated dimers in the fitting set and that these orientations are not important in the crystalline phase. However, relaxation to the balanced dataset of pseudo-random dimers and some configuration scan points does improve the fit to the total energies, as is just detectable in Fig. 4 for the configuration scans. Each component is well represented by the model potential, apart from the polarization energy at a short range. The real test of the model potentials is in Fig. 5 for the independent test of the crystal dimer dataset. The relaxation has somewhat improved the model potential for the crystal dimers, although it was expensive to generate sufficient total intermolecular energies {E_int[SAPT(DFT)]} so that the pseudo-random dimers covered sufficient orientation space for the relaxation step to be at all worthwhile.

The overall relative success of the different model potentials in representing the second-order energies is shown in Table II. This confirms that the model potentials give a good representation of the SAPT(DFT) energies, including the relative orientations sampled in crystal structures, as well as to the points used in the weighted relaxation process (see 2.5 of the supplementary material), many of which are significantly up the repulsive wall.

The CSP lattice energy landscape for the final potential DIFF-srL2(rel) is shown in Fig. 6. The most stable and densest structures are all very similar to form III with the same 30-molecule coordination sphere. The global minimum structure is considerably more stable than the lattice energy minima reached by optimizing from the experimental form III structure, by 12 kJ mol⁻¹, but overlays 30 molecules with an RMSD₃₀ of 0.6 Å. It seems likely that many, if not all, structures with a 30 molecule overlay with form III are in shallow potentials energy wells and would become the same dynamic structure in a MD simulation under ambient conditions, i.e., correspond to form III if this lattice energy landscape (CSP_0) were reduced by allowing the molecular motions.⁷⁵ The lowest energy structure that does not approximate to form III (the Pbca structure, TNB No. 55 in Fig. 6) has two layers of molecules in common with form III but differs by the third layer (see Fig. 12 of the supplementary material).

The majority of computer-generated structures are less dense than form III, although denser than form II, and more stable than form I. A selection of the most diverse structures is shown in Table III, which emphasizes that only some of the crystal structures contain a dimer related to the gas phase minima (see Fig. 11 of the supplementary material), and none contain any of the three most stable dimer configurations. Hence, the CSP with DIFF-srL2(rel) is successful, within the limitations that only Z′ = 1 structures were generated, and so, forms I and II could not be found in the search.

The balance of contributions to the lattice energy varies considerably over the crystal structures (Table III), particularly the dominant dispersion energy, although this is partially offset by the variation in the short-range energies, as the denser structures are most stabilized by the dispersion. The polarization energy is significant, approaching half the electrostatic energy in some cases. Differences in the polarization energy account for a large proportion of the variation in lattice energies between the structures that approximate form III.

The relative stability of the CSP generated structures is very dependent on the potential energy model, as shown by the different crystal energy landscapes in Figs. 8–10 of the supplementary material and Fig. 7. The CSP study with the model potential commonly used for pharmaceuticals [FIT(GDMA,L4)] is almost useless, with all the known forms being unrealistically less stable than the hypothetical structures. Changing the electrostatic model to use the ISA partitioning together with a better quality charge density brings the structures within a more plausible energy range.

All the non-empirical DIFF potentials give the structures in a plausible energy range and have form I as metastable to forms II and III but differ in the relative rankings of forms II and III and the selected unobserved structures. The introduction of repulsion anisotropy introduces some reranking. The relaxation of the parameters by fitting to the second-order energies makes some difference, changing the relative stabilities of forms II and III and the structures that approximate form III from other hypothetical structures.

All the non-empirical models were lattice energy minimized without the polarization energy, which was then added as a single-point correction. Hence, the changes in relative energies include a contribution from small changes in the structure. The exception is the comparison between the final potential DIFF-srL2(rel) before (see Fig. 8 of the supplementary material) and after the polarization term (Fig. 6). This shows that the polarization energy does make a considerable difference in spreading out the energies, including the energy range of the structures that approximate form III. However, the improvement of the electrostatic model within the model potential has the greatest effect on the CSP study.

The final derived DIFF-srL2(rel) potential is highly successful in that it represents the intermolecular energies as calculated by SAPT(DFT) to a good accuracy over a wide energy range (Table II). The use of distributed multipoles, polarizabilities, and dispersion coefficients guarantees the correct long-range behavior, but the balance of the long and short-range terms is well represented in the potential energy wells for the low energy dimers (Fig. 4). This model potential has successfully translated to the solid state, providing a realistic CSP study for TNB. It is a significant improvement on a CSP study using the semi-empirical model potential [FIT(GDMA,L4)], which has been mainly used for pharmaceuticals. Indeed, early CSP studies of some chloronitrobenzenes³⁹ struggled to find an acceptable model potential. The improvement in the model potential was most dramatic on improving the distributed electrostatic model as the change in both the CSP (Fig. 7) and gas-phase dimer (Table I) results was greatest in going between FIT(GDMA,L4) and FIT(ISA,L4). However, the relative ranking of the minima in both phases was sensitive to the differences between the DIFF potentials.

This study has assumed that TNB is rigid and remains in the highly symmetrical planar conformation of the ab initio optimized isolated molecule upon crystallization. Analysis of the known polymorphs in Fig. 1 shows that this is an approximation, and the nitro groups twist somewhat to give non-planar structures. The most readily crystallized structure, form I, has a nitro group torsion that is more than 30° from planar in one of the two independent molecules. We estimated the effect of this conformational adjustment on crystallization by extracting the experimental molecular conformations from the crystal structures, recalculating the ISA multipoles and polarizabilities, and reoptimizing the crystal structures keeping the observed conformations rigid. This significantly improved the reproduction of form I and stabilized the lattice energies for forms II and III (see 3.1.1 of the supplementary material). However, even without including any energy penalty for the conformational change, the recently discovered forms II and III were more stable than form I. This sensitivity to the conformational change limits the conclusions that can be drawn from the CSP study of TNB. However, it seems very likely that form I is kinetically favored under most crystallization conditions, but the recently discovered forms II and III are more stable and that further polymorph screening and detailed crystallography may reveal more polymorphs or disorder in TNB. Ideally, a further CSP would be done in which the NO₂ groups of TNB could change conformation in response to the packing forces. This could be done by analytically rotating the multipoles during the structure relaxation,⁴² but a final revaluation of at least the multipoles for each crystalline conformation would be needed for the final lattice energies.⁴² 

A second approximation in the use of DIFF pair potentials in condensed phases is the modeling of the non-pairwise additive terms, i.e., the many-body terms that involve three or more molecules. The only non-additive term included in the DIFF potentials is the polarization energy and its exchange counterpart. The polarization contribution to the lattice energy is relatively large (Table III), several times the polarization energy per molecule for the dimers extracted from crystal structures $EPOL(2)[SAPTDFT]$ $<3.5 kJ mol−1$⁠; see Fig. 2 of the supplementary material}. This is in contrast to simple ionic solids, where the net electrostatic field at each ion in the crystal can be zero by symmetry so that the dipolar polarization energy is zero, although there would be a huge polarization energy for a pair of ions. The variation of the polarization energy between different crystals is significant (Table III). However, the fitting of a molecular damping parameter β_pol was not very satisfactory for TNB (see Fig. 2 of the supplementary material), and this could have made the polarization energy sensitive to the distance in the nearest neighbor contacts. This coupled with the uncertainty in the applicability of the $δintHF$ correction and not including the polarization energy in the lattice minimization could have led to an overestimation of the relative differences in the polarization contribution to the lattice energies of the form III type structures (Fig. 6).

The modeling of the polarization within the specific crystal structure via atomic polarizability tensors⁷⁶ is an advance on the isotropic polarizabilities within the AMEOBA force-field or the calculation of the distributed multipoles within a polarizable continuum model with a dielectric constant typical of those of organic crystals. The latter approach is often used in CSP studies,⁷⁷ as it often improves the relative energies between structures with different types of hydrogen bonding. We note that the molecular three-body Axilrod–Teller–Muto dispersion contribution was not included in modeling the crystals, and this would have provided a repulsive term that would share the sensitivity to the nearest neighbor atom–atom distances and at least partially cancel the attractive polarization term. Hence, it is important to extend our understanding of the many-molecule terms in condensed phase modeling, noting that the physics of many-body (non-additive) contributions limits the accuracy of the lattice energy evaluated from a very accurate two-body model such as the DIFF.

This CamCASP-derived TNB DIFF model is the first time the ISA partitioning of the molecular charge density into atoms has been used for all the components of the potential for a molecule larger than water.⁴ These atomic charge densities were used to derive the atomic multipoles,⁴⁴ polarizabilities,⁵⁹ and dispersion⁵⁹ parameters and, via the distributed density overlap model, to provide a first estimate of the repulsion potentials and form of anisotropy.¹⁶ The more rigorous ISA partitioning method provides maximally spherical atoms, thus minimizing the number of important anisotropy terms. From this approach, we were able to identify the significant anisotropy on the O and N atoms in TNB (Fig. 3) and thereby construct a model with short-range anisotropy on these two atom types to significantly improve the model potential (see Table V of the supplementary material). The ISA partitioning also allows the use of better wavefunctions and hence improved polarizabilities and $C6ικ$⁠, $C8ικ$⁠, and $C10ικ$ dispersion coefficients, which is particularly important given how the dispersion dominates the lattice energies (Table III). The use of the ISA partitioning to represent each individual component of the interaction energy is expected to reduce the errors in each contribution and hence reduce the need to relax the parameters to fit the total energies. Certainly, the unrelaxed model potential was already a good approximation for TNB and resulted in nearly the same relative energies of the crystal structures in the CSP study (Fig. 7). This is a promising finding for the development of potentials for large molecules, as the cost of deriving the DIFF-srL2(rel) potential was significantly more than the DIFF-srL2(norel). This is due to the need to calculate several hundred second-order dimer interaction energies, sampling a wide range of relative orientations, in order for the relaxation to be unbiased and hence an improvement (see 2.5 of the supplementary material). It is noteworthy that there are relatively few TNB gas-phase dimer structures (minima in the dimer potential energy surface, Table I), and so, the properties of TNB in the gas-phase would sample a very limited portion of the dimer intermolecular potential surface. The condensed phases sample a wide range of relative orientations (Table III) although the relative orientations in van der Waals contact are rarely those of the gas-phase dimers. Indeed, the non-empirical potentials have the further advantage of accurately representing the repulsive wall, which would be sampled when the crystals are compressed.⁵ The use of the theory of intermolecular forces in the DIFF approach produces model pair potentials that should describe the intermolecular pair potential surface over the configurations sampled in all phases. This is a significant advantage over methods that depend critically on the sampling of configurational space, such as machine learning or fitting to total energies.^78–80 

An anticipated advantage of the ISA partitioning of the charge density is that the changes in the atomic potential parameters with conformation may be limited to the physical redistribution of charge with conformation⁴² (i.e., the changes that go beyond the geometric changes that can be done using analytical tensor rotation). This could greatly improve the efficiency of modeling the changes in conformation that occur on crystallization. The CrystalPredictor^81,82 and CrystalOptimizer^82,83 approach to CSP models the effects of conformational changes by building up databases of conformational energies, associated forces, and distributed multipoles as a function of the torsion and bond angles that vary between polymorphs and using local approximate model interpolations. The costs of using the more accurate molecular charge densities found necessary in this study may be mitigated by the better transferability of the potential parameters with conformation. It is also possible that a DIFF-style intermolecular force-field and a separately parameterized atom–atom intramolecular force-field for the conformational energy penalty⁸⁴ could be used to efficiently generate and rank crystal structures in a CSP study. (The approximation that the same isotropic parameters can be used for intermolecular and intramolecular forces is a major limitation of traditional transferable force-fields, which has limited their usefulness in CSP studies.^84–86) There is a great need to have force-fields that are accurate enough to give realistic relative stabilities of observable and hypothetical crystal structures generated by CSP and can be used in molecular dynamics calculations to simulate the crystals under ambient temperature and determine the effect of varying the temperature and pressure.⁷⁵ The DIFF approach could provide such force-fields: DIFF-srL0(norel) potential could, in principle, be used in TINKER⁸⁷ after adjusting the atomic multipoles to have a maximum rank of 2.

We have derived a set of progressively more accurate and detailed non-empirical model intermolecular pair potentials for trinitrobenzene using the CamCASP and ORIENT programs. Central to these potentials is the Iterated Stockholder Atoms (ISA) partitioning of the molecular charge density into atomic charge densities using the basis-space ISA algorithm (BS-ISA) implemented in the CamCASP suite of codes. The potentials include analytic long-range parameters from the distributed atomic multipole tensors, distributed dipolar polarizability tensors, and distributed isotropic $C6ικ$⁠, $C8ικ$⁠, and $C10ικ$ dispersion coefficients, corresponding to a high-quality molecular charge density, and frequency-dependent molecular response functions. The long-range parameters were derived using the BS-ISA and ISA-Pol algorithms using the CamCASP code. The short-range potential parameters were determined through a hierarchical algorithm based on the distributed density-overlap model using the BS-ISA expansions for the atoms-in-a-molecule densities and SAPT(DFT) interaction energies. This algorithm involves multiple computationally well-defined steps that established a minimal set of molecule-specific atomic anisotropy terms. The initial short-range parameters in the analytical models were determined by fitting to an extensive set of computationally inexpensive SAPT(DFT) first-order exchange and penetration energies, $Eexch(1)$ + $Epen(1)$⁠. The Tang–Toennies short-range damping of the dispersion and polarization used a molecular parameter that was fitted to second-order dispersion and polarization energies, respectively, although these parameters could have been estimated. The resulting model potential, DIFF-srL2(norel), gave a useful crystal structure prediction (CSP) study for TNB and reproduced the dimer interactions energies very well.

The SAPT(DFT) interaction energies up to second-order were then considered, with the charge-delocalization term found to be negligible and the $δintHF$ term found to be insignificant and therefore neglected. A dataset of SAPT(DFT) dimer energies, $Eint=Eelst(1)+Eexch(1)+EIND(2)+EDISP(2)$⁠, was used to relax the isotropic short range parameters in order to absorb the errors in modeling the separate contributions and to also absorb the charge-delocalization energy. This final model, DIFF-srL2(rel), further improved the ability to model an independent dataset of total interaction energies of crystalline dimer configurations (not used in the potential derivation) with an RMSD of 0.5 kJ mol⁻¹. The improvement was sufficiently small to suggest that using ISA partitioning and SAPT(DFT) first-order calculations to produce a model for each significant contribution to the intermolecular potential could be extended to producing an accurate anisotropic atom–atom model intermolecular potential for much larger molecules than TNB.

The non-empirical DIFF intermolecular potentials were used in successful crystal structure prediction (CSP) studies of TNB, despite sampling a very different set of relative orientations of the molecules than the minima in the dimer potential energy surface. The relative energies of the observed and computer-generated crystal structures of TNB were reasonable, unlike those generated using an empirical model potential frequently used in CSP for pharmaceuticals. The study supports other evidence that the recently discovered form III of TNB is the most stable polymorph.

The chief weaknesses of the intermolecular potential model, from both theory and the results, were in the treatment of the non-pairwise additive (many-molecule) terms and in the assumption that the molecule was rigid. We were unable to include non-additivity in the dispersion model, and the use of a single damping parameter in the polarization model may have resulted in the many-body polarization energy being too sensitive to the crystal structure. However, the rigid molecule approximation was probably the more important limitation for modeling the crystal structures of the known polymorphs and the CSP study, as small changes in the NO₂ torsion angles can significantly improve the packing of the molecules in some crystal structures.

This study makes significant progress toward a methodology for deriving non-empirical anisotropic atom–atom force-fields for organic molecules that would be capable of realistic modeling of all phases in molecular dynamics simulations.

See the supplementary material for the justification for choice of basis set for SAPT(DFT), the TNB atomic properties used in the DIFF potentials, the definition of polarization energy and damping the multipolar polarization energy, the distributed dispersion model, the derivation of the repulsion anisotropy, the refinement of the potential parameters, and discussion of the omitted terms and for further details of the crystal structure modeling of TNB, reproduction of experimental polymorphs, effect of conformational flexibility, sensitivity to the force-field, CSP results with different force-fields, and the CSP generated structures for DIFF-srL2(rel).

This paper honors Ruth Lynden-Bell FRS, an outstanding theoretical chemist, who has mentored and inspired many students, including two of the authors, to pursue careers in academia during her career at the University of Cambridge and Queen’s University Belfast.

We thank Professors Claire Adjiman and Costas Pantelides for allowing us to use the Crystal Predictor code, and Professor Anthony J. Stone for his long-standing collaborations with both A.J.M. and S.L.P., which have resulted in many of the developments given in this paper, and also for the adaption of ORIENT. We also thank Mark Storr and other scientists at AWE and the M3S Centre for Doctoral Training (EPSRC Grant No. EP/G036675/1) for financial and general support to A.A.A. The CSP computational software is being developed under Grant No. EP/K039229/1.

The model potentials are available as input files for ORIENT and DMACRYs from https://app.ph.qmul.ac.uk/wiki/ajm:potentials:start. The CSP generated crystal structures are available from the corresponding author on reasonable request.