Calculation of Diamagnetic Susceptibility Tensors of Organic Crystals: From Coronene to Pharmaceutical Polymorphs

Understanding why crystallization in strong magnetic fields can lead to new polymorphs requires methods to calculate the diamagnetic response of organic molecular crystals. We develop the calculation of the macroscopic diamagnetic susceptibility tensor, χcryst, for organic molecular crystals using periodic density functional methods. The crystal magnetic susceptibility tensor, χcryst, for all experimentally known polymorphs, and its molecular counterpart, χmol, are calculated for flexible pharmaceuticals such as carbamazepine, flufenamic acid, and chalcones, and rigid molecules, such as benzene, pyridine, acridine, anthracene, and coronene, whose molecular magnetic properties have been traditionally studied. A tensor addition method is developed to approximate the crystal diamagnetic susceptibility tensor, χcryst, from the molecular one, χmol, giving good agreement with those calculated directly using the more costly periodic density functional method for χcryst. The response of pharmaceutical molecules and crystals to magnetic fields, as embodied by χcryst, is largely determined by the packing in the crystal, as well as the molecular conformation. The anisotropy of χcryst can vary considerably between polymorphs though the isotropic terms are fairly constant. The implications for developing a computational method for predicting whether crystallization in a magnetic field could produce a novel or different polymorph are discussed.


INTRODUCTION
Though diamagnetism is a much weaker effect than paramagnetism or ferromagnetism, the advent of stronger continuous or pulsed magnets (currently up to 45.5 T continuous, 1 100 T pulsed) makes it much easier to study the response of diamagnetic materials, ranging from polymers to water, to an external magnetic field. Possibly, the most dramatic example is the levitation of a living frog in an inhomogeneous magnetic field, 2 as the (diamagnetic) water in the frog moves to the lower magnetic field region. The socalled Moses effect, the deformation of the surface of water or other diamagnetic liquids in a strong and steady magnetic field, has been explained as the balance between gravity and magnetic forces originating from inhomogeneity of the field. 3,4 This effect can also be exploited to separate diamagnetic particles by their densities. 5 Other applications in microfluidics use magnetic fields for the manipulation and direction of the movements of small diamagnetic particles. 6 Strong magnetic fields have been used in protein crystal growth in a microgravity environment. 7 More applications of applying a magnetic field to diamagnetic materials can be found in the reviews by Yamaguchi and Tanimoto, 8 and Maret and Dransfeld. 9 Diamagnetic molecular crystals can become oriented in magnetic fields, 10 and this is being exploited to aid the solution of crystal structures by aligning microcrystals in powder samples. 11 Crystals grown in magnetic fields have also been observed to have a different morphology and nucleation rate. 12 The much stronger magnetic fields around white dwarf stars, up to 10 5 T, have been shown to promote a new bonding mechanism, perpendicular paramagnetic bonding, in hydrogen. 13 When a diamagnetic molecule is put in a laboratory external magnetic field B, the motion of its electrons and nuclei is perturbed. The first-order response of the diamagnetic molecule to B is the creation of induced electron current density j (as illustrated for benzene in Supporting Information Figure S1), which subsequently induces a magnetic field, according to the Biot−Savart law, with a component opposite to B. The response of the molecule can be described as an induced magnetic dipole m, with its direction and magnitude determined by its molecular magnetic susceptibility (magnetizability), χ mol , a second-rank tensor. 14 When organic molecules crystallize, the magnetic response depends not only on the number of molecules and their magnetic response but also on the way the molecules are orientated within the crystal. Hence, different polymorphs 15 may have significantly different induced magnetic moments m when present in a magnetic field, which will be reflected in their different crystal magnetic susceptibilities (χ cryst ).
A key issue is whether a magnetic field can change which polymorph crystallizes, as reported for terpyridine. 16 Indeed, a previously unknown polymorph (β) of coronene was grown in a 1 T magnetic field 17 and is more thermodynamically stable than the readily crystallized γ form at low temperatures. Because finding all the stable polymorphs of a pharmaceutical is a key step in pharmaceutical development, 18,19 and the existence of the new coronene polymorph could be predicted by a crystal structure prediction (CSP) study, this raises the question of whether crystallization in a strong magnetic field could produce the first sample of a novel polymorph of a pharmaceutical. 20 CSP studies on pharmaceuticals usually generate more thermodynamically plausible structures than the polymorphs found, despite the extensive industrial screening that is usually carried out in pharmaceutical development, 21,22 and therefore unusual experiments can be undertaken to target a particular structure predicted by the CSP study to be more stable than the more readily crystallized forms. 23−26 In order to understand whether the crystallization of different polymorphs of a pharmaceutical could be affected by a strong magnetic field, we need to investigate the diamagnetic response of pharmaceutical polymorphs.
For diamagnetic systems, the induced magnetic moment m interacts with the external magnetic field, B, to change the total energy of the system by E B = −∫ m(B)·dB > 0, where m(B) denotes the dependence of m on B. This destabilizing interaction leads a diamagnetic molecule or crystallite to either move toward areas of lower B in a nonuniform magnetic field or orientate itself to align the smallest possible component of χ (and so the smallest m) along the direction of B in a uniform magnetic field. 9 Because the anisotropy of magnetic susceptibility determines the tendency of a molecule to reorientate in a magnetic field, 9 a full knowledge of the χ tensor is needed, not just the isotropic part (χ iso ). In order to understand how a magnetic field could affect the competitive nucleation and growth of polymorphs of organic crystals, we need to quantify the crystalline magnetic susceptibility tensor, χ cryst , and understand how this depends on the molecule, its conformation, and the molecular packing in the crystal. Reliable experimental measurements of magnetic susceptibility are difficult and most available measurements 27,28 are for the isotropic term only, measured in gases, liquids, or for powders, with the notable exception of the measurements of the anisotropy of χ of an unspecified coronene polymorph. 29 It is even more challenging to compare the crystal magnetic susceptibilities of different polymorphs because of the difficulty in growing and quantifying the size of suitable single crystals. Indeed, for molecules, extensive research in first principles calculations of molecular magnetizability over the past two decades 30 has reached the stage that the best calculated values are recommended for calibration of experimental measurements. 31 There are extensive calculations of molecular magnetic properties because of the link to aromaticity and correlations with other properties. 32−36 Although quantum chemical packages, for example, DALTON 37 and Gaussian, 38 are readily available for the calculation of the magnetic susceptibility of molecules, 39 most applications have focused on the magnetizability of small rigid molecules or fullerenes such as C 60 . 31,40,41 Virtually nothing is known about the magnetic susceptibility of pharmaceuticals, which generally comprise various aromatic rings and other πdelocalized groups with flexible linking groups between them. Such molecules can be observed in multiple conformations in their crystalline forms, leading to conformational polymorphs. 42 The molecular conformation in the most stable crystal form can have its aromatic groups in a substantially different relative orientation from the most stable conformation of the isolated molecule. 43 Thus, the diamagnetic response of any phase from which a pharmaceutical crystallizes, for example, a supersaturated solution, may be affected by changes in the conformation as the crystallization occurs.
This paper pioneers the calculation of crystal magnetic susceptibility tensors for organic pharmaceuticals from periodic density functional calculations and how they can be estimated from molecular magnetizabilities, provided that the conformational dependence is taken into account. The crystal magnetic susceptibility can be obtained from CASTEP 44 using the routines developed for calculating solid-state NMR chemical shifts. We then introduce the use of tensor addition (TAmol) as an approximation to calculate χ cryst from χ mol . Using periodic density functional theory (DFT) calculated χ cryst as a reference, the tensor addition method is assessed for 21 organic crystal structures covering the experimentally known polymorphs of nine diverse organic molecules. This set of calculations leads to the conclusion that different polymorphs of the same molecule can have significantly different magnetic anisotropies, Δχ an cryst , in the crystal. This anisotropy can be linked to the electron delocalization in the molecule, its conformation, and the crystal packing of the specific polymorph. In Section 4, we discuss how the tensor approximation allows CSP methods 45 to be extended to produce a magnetic property−crystal structure− energy landscape to establish whether polymorphs have sufficiently distinct magnetic properties that their crystallization could be affected by a magnetic field.

THEORY AND CALCULATION OF DIAMAGNETISM
OF ORGANIC PHARMACEUTICAL MOLECULES AND CRYSTALS As most organic pharmaceuticals are closed-shell molecules and contain only lighter atoms, only the orbital motions of the electrons couple with the magnetic field, which can be treated as a small perturbation. We only consider time-invariant and spatially homogeneous magnetic fields and use the Gaussiancgs unit system 45 (Supporting Information Section 9) for ease of comparison with the experimental data. The first-order response of a closed-shell, diamagnetic molecule to an external field, B, is the creation of the induced electron current density, j, which leads to an induced magnetic dipole m = χ·B that always has an antiparallel component to B. The components of magnetic susceptibility can be determined from the induced current densities j by 46 r r j r cB where c is the speed of light, B = |B| is the magnitude of the external magnetic field, and r is the real-space position vector relative to an arbitrary origin. j β (1) (r) is the first-order induced current density component in the Cartesian β direction, which can be calculated with coupled perturbed Hartree−Fock. 47 Calculations of the magnetic-field-induced current densities in organic molecules have been widely used to characterize ring currents 48 and aromaticity in large organic systems 32,49 including toroidal carbon nanotubes. 50 The Journal of Physical Chemistry A pubs.acs.org/JPCA Article The magnetizability χ can also be defined as the secondorder derivative of the total energy E of the system with regard to B The calculated diamagnetic susceptibility can be diagonalized to obtain three eigenvalues χ i (i = 1...3), which are labeled by magnitude with χ max ≤ χ mid ≤ χ min < 0. The tensor can be characterized by the isotropic average, χ iso = 1/3(χ max + χ mid + χ min ) < 0, and the anisotropy Δχ an = χ max − χ min < 0. The χ tensor can be visualized as an ellipsoid, 14 with its three eigenvectors parallel to and proportional to the three axes of the ellipsoid. Only when the external magnetic field coincides with the directions of one of the three eigenvectors of the χ tensor is the induced magnetic moment m antiparallel to B. We use χ for molar magnetic susceptibility, in units of 10 −6 cm 3 / mol or 1 cgs-ppm, defined per mole of molecules, to facilitate comparison between polymorphs with different numbers of molecules in the unit cell and isolated molecules. The dimensionless volume magnetic susceptibility, κ, is usually used for the magnetic response of a unit volume of material (gas, solution, or solid), where V m is the molar volume of the specific sample or the polymorph, ρ is the density, and M is the molecular weight. The calculation of magnetic properties requires consideration of gauge invariance. 30,51 Many methods have been proposed to address the need for magnetic susceptibility and nuclear shielding tensors to be independent of gauge origin, with earlier methods largely superseded by GIAO (gauge including atomic orbitals), 30,39,52−54 which explicitly uses gauge-including atomic orbitals (London atomic orbitals) as basis functions ξ p (B) = ξ p (0)exp[(−i/2c)(B × R p )·r], where ξ p (0) is the field-free basis function located on nucleus p and R p is the vector pointing from an arbitrary gauge origin to the nucleus. London atomic orbitals effectively shift any choice of global gauge origin back to the atomic centers, which, for a one-electron one-center problem, gives the correct solution to the first order. 55 This constrained choice of local origins removes all reference to the global gauge origin, thus ensuring gauge-origin independent results and fast basis set convergence. 51,56 The calculation of molecular magnetic susceptibilities of small molecules can be carried out with considerably more accurate charge densities, such as Møller−Plesset, 57 MCSCF, 58 coupled-cluster, 59 and full CI, 60 than could be afforded for larger pharmaceutical molecules. These accurate charge densities cannot be used for crystals, as they are not implemented in periodic electronic structure codes. The magnetizabilities χ mol calculated from many types of charge densities have been benchmarked against coupled-cluster and experimental values for cyclopropene and smaller molecules, showing the importance of including diffuse functions in the basis set. 57 Following the development of DFT methods, Lutnaes et al. established a benchmark database for a set of 28 small molecules at the CCSD(T) level, extrapolated to the complete basis set limit. 31 This showed that nearly all GGA and hybrid functionals perform similarly to the HF method with the range-separated hybrid CAM-B3LYP functional slightly better. 31 However, their work only considered the isotropic term of magnetizability, and it has been suggested that χ iso is not as sensitive to electron correlation as the anisotropy. 30 Our choice of method had to be one that would be affordable for pharmaceutical molecules and crystals, limiting the choice to HF and DFT methods. A simple analysis (Supporting Information, Section 4) examining the effects of choice of density functionals, size of basis sets, and the CGST or GIAO methods on the χ mol of pyridine using Gaussian 09 38 largely confirmed the findings in the literature. 31,40 Hence, χ mol was calculated for a range of molecules, ranging from rigid aromatic to flexible pharmaceutical molecules, using the GIAO method and 6-31G(d,p) basis set, with both the PBE and PBE0 functionals: the former for direct comparison of χ cryst calculated with the periodic PBE method and the tensor addition TAmol method and the latter to estimate the effect of the improved hybrid functional on χ.
2.1. Periodic DFT Calculations of Crystal Magnetic Susceptibility. For crystals, the use of plane wave basis sets changes the gauge invariance problem. Although formal expressions for χ cryst have been derived in the 1960s, 61,62 practical calculations were only first reported just over a decade ago for solid neon and diamond using the LDA functional 63,64 in developing the calculation of solid-state NMR nuclearshielding constants. This approach uses a sum-over-states perturbation expansion for the susceptibility for a magnetic field with a finite wavevector, with the macroscopic χ cryst being the limit for a field of infinite wavelength (i.e., uniform B). 63 Strictly speaking, this approach is only valid for all-electron calculations. In almost all practical periodic DFT codes, pseudopotentials are used to remove both the inner nodal structures of valence electrons and the electron contributions in the core area to greatly reduce the number of plane waves needed for calculations. However, for calculations of NMR nuclear shielding including any second-row elements, the inner nodal structure of the full electron density cannot be ignored. The gauge-including projector-augmented wave approach, 65−67 a variation of the projector-augmented wave approach, 68 adds the para-and diamagnetic contributions from core electrons to those from the valence electrons. The induced electronic currents j as calculated by CASTEP for a molecule of benzene in a large periodic box are shown in Supporting Information, Figure S1 which compares well with the molecular calculation in ref 32. The equivalence of the values of χ mol for aromatic molecules calculated by the molecular method (Gaussian 09) and by the periodic plane wave method for each molecule in a sufficiently large periodic box by CASTEP is demonstrated in Supporting Information, Figure S3.
All crystal structures considered in this paper were first optimized with CASTEP 44 using the PBE functional 69 and Tkatchenko and Scheffler's (TS) dispersion correction scheme, 70 a methodology that is widely used for modeling crystal structures of larger organic molecules, 71,72 particularly in CSP studies. 20 Plane wave cut-off energies and k-point grids were carefully selected for the polymorphs of each molecule to ensure tight convergence of the total energy (Supporting Information, Section 5) and converged to a maximum force of less than 0.001 eV/Å. The crystal diamagnetic susceptibilities χ cryst of the polymorphs were calculated using this PBE charge distribution. The effects of improved intramolecular electron correlation on χ cryst were roughly estimated using the PBE0 functional and the tensor addition method.  73 Hence, there is a need to investigate approximate methods that could be used to calculate the magnetizability of different polymorphs more efficiently. A first approximation to the crystal diamagnetic susceptibility χ cryst is the tensor addition of molecular diamagnetic susceptibilities χ mol , which would be exact if the intermolecular interactions within the crystal did not affect the magnetic response of the individual molecules. With this assumption, the magnetic response of a crystal to a field B is just the tensor sum of the individual response of every molecule in it and is determined by the relative orientations of all the Z molecules in the unit cell to the field, assuming the crystal to be perfect and infinite. The mathematical derivation in terms of the three-dimensional rotation matrices R i between the crystal-fixed and the ithmolecule-fixed coordinate systems is given in the Supporting Information, Section 4. The result for the crystal magnetic susceptibility in the crystal-fixed coordinate system (χ cryst | cryst-fixed ) is, in terms of the tensors or equivalently in terms of the tensor components as The calculation of χ mol-i of the ith molecule in the crystal in its molecule-fixed coordinate system should use the same conformation as in the crystal. We have chosen to use the same molecule-fixed coordinate system as the organic crystal modeling program DMACRYS, 74 as this provides the coordinates in a user-defined molecule-fixed axis system for each independent molecule in the crystallographic asymmetric unit cell along with the rotational matrices in eq 3a. The coordinates can be used directly to calculate χ mol of each molecule using Gaussian 09. However, it is important to note that coordinates of the molecules generated by mirror planes, glide planes, and inversion centers in the crystal differ from those of the molecule defining the asymmetric unit cell by a change of sign of the z-axis in order to maintain a right-hand axis system. Thus, χ mol for molecules generated by these symmetry operations also need transforming by changing the sign of χ xz (=χ zx ) and χ yz (=χ zy ). A schematic of the process of calculating the diamagnetic susceptibility of a crystal using the tensor addition method in contrast with direct calculation using periodic electronic structure methods is shown in Figure  1.
For each crystal structure, either PBE-TS optimized or experimental, Z′ calculations of molecular magnetizabilities were performed on each molecule in the asymmetric unit cell, extracted from the crystal structure with NEIGHCRYS, 74 which also generates rotational matrices for all molecules in the unit cell. These matrices and the χ mol tensors are combined using the specifically written code TAmol. For a direct comparison to χ cryst calculated with periodic DFT, χ mol was calculated with the PBE functional and a 6-31G(d,p) basis set using Gaussian 09 38 and also contrasted with χ mol calculated with PBE0/6-31G(d,p). The differences between the experimental and PBE-TS optimized structures were usually small (Supporting Information, Table S4), and the effects of these differences on calculated χ cryst are shown in Supporting Information, Tables S5 and S6. CCDC software Mercury (ver. 4.1.3) was used for visualizing and comparing crystal structures and generating crude morphology models for coronene polymorphs using the Bravais−Friedel−Donnay− Harker (BFDH) method which is based on the interplanar spacing. 75 rmsd n is the minimum root-mean-square (rms) difference between two crystal structures calculated by overlaying all nonhydrogen atoms in a crystalline cluster of n molecules. Figure 1. Diagram of the two approaches used for the calculation of χ cryst from a crystal structure: along the red arrow, χ cryst is calculated using the periodic DFT method; alternatively, it can also be calculated with much lower cost with the tensor addition approach that only involves Z′ calculations of χ mol . The most CPU-intensive steps are the electronic structure calculations in yellow boxes.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article 3. RESULTS 3.1. Diamagnetic Susceptibilities of Pharmaceutical Molecules. The physical origin of the diamagnetic susceptibility is the induced currents within the molecule, and this leads to additivity relationships for different atoms and functional groups with particularly large contributions for conjugated, delocalized groups, such as aromatic rings. 33,34,76 Empirical methods 28 have been used for estimating the isotropic magnetic susceptibility of molecules (Supporting Information, Section 3) from the contribution for each atom and additional terms for specific types of bonding. This principle extends to the pharmaceutical molecules in that the isotropic magnetizabilities for acridine, p-methyl chalcone, and carbamazepine ( Figure 2, Supporting Information, Table S5) are relatively similar as each molecule has two benzene rings linked by another conjugated system. However, the anisotropy of their magnetizabilities varies considerably more. For less symmetrical organic molecules with aromatic rings oriented in different directions, such as carbamazepine and the slightly nonplanar p-methyl chalcone, the anisotropy of χ mol is reduced. In a planar molecule such as acridine, χ max lies normal to the aromatic ring plane, whereas in carbamazepine, it lies in the plane where the multiple rings have the largest projection leading to a smaller χ max .
The results in Figure 2 suggest that the anisotropy of χ mol will be sensitive to the molecular conformation. Many organic molecules, in particular pharmaceuticals, comprise multiple aromatic rings which are flexibly linked and usually have conformations far from planarity. The fenamates exhibit a wide range of conformations in many of their polymorphs, 24 and the small energy differences with the change in molecular conformation are difficult to calculate accurately 77 because of the subtle balance between electron delocalization and steric clashes. The eigenvalues of χ mol for flufenamic acid (Figure 3a), calculated for the conformations generated by scanning the most flexible torsion angle ϕ and letting the rest of the molecule adapt, show that χ max and χ mid vary considerably with ϕ. The largest magnitude of χ max (∼−220 cgs-ppm) corresponds to the conformation with both aromatic rings nearly coplanar (ϕ ≈ 0 or 180°). When the two aromatic rings are perpendicular to each other (ϕ ≈ 90°), χ max is close in value to χ mid . In contrast, χ min varies little (∼−115 cgs-ppm) as it lies in the direction that remains roughly edge-on to the two aromatic rings in the molecule and so is barely affected by the change in ϕ between the two aromatic rings. A similar behavior is observed in the eigenvalues of χ mol for p-methyl chalcone, where there is a complex variation of the three components of χ mol (Figure 3b) due to the OC−CC group, which changes from a syn (ϕ ≈ 0°) to an anti (ϕ ≈ 180°) conformation. In both cases, χ iso is approximately constant with the conformational change.
3.2. Calculations of the Diamagnetic Susceptibility of Organic Crystals by Periodic DFT and Tensor Addition Methods. The diamagnetic susceptibility tensors of the polymorphs of coronene (Table 1) are clearly determined by the relative orientation of the aromatic rings, when superimposed on the crystal structure ( Figure 4). The calculated χ max cryst (per molecule) of each polymorphs is smaller than χ max mol for the isolated molecule, as in neither polymorph are the molecules stacked with their molecular symmetry axes perfectly aligned. The β form, which was first discovered by crystallization in a magnetic field, 17 has a far smaller angle between the ring planes than the γ form which usually crystallizes, and it therefore has a far larger χ max cryst and larger anisotropy. This is what would be expected from the tensorial addition of χ mol (Figure 4) if one considers the difference in angles between the normals to the molecular planes. Indeed, compared to the periodic PBE calculated χ cryst , the error in using tensor addition is small (Table 1 and Supporting  Information, Table S5). Using PBE0 χ mol in the tensor addition method consistently leads to insignificantly larger magnitudes of χ cryst (Table 1 and Supporting Information, Table S4).
The intermolecular interactions between coronene molecules are relatively weak compared with those involved in pharmaceuticals; therefore, the coronene polymorphs might be considered a favorable case for the neglect of the redistribution of the electron density among the molecules, which is inherent in the tensor addition method. Table 1 shows that the tensor addition method is also a good approximation for carbamazepine form III, which, like most carbamazepine polymorphs, contains the relatively strong double amide hydrogen bonding 78 that is considered a reliable synthon. 79 The tensor addition method also works well for a system with strong π···π stacking interactions, as shown in Table 1 for form IV of the benzene/hexafluorobenzene (C 6 H 6 /C 6 F 6 ) cocrystal in which benzene and hexafluorobenzene molecules alternate to form stacks. This cocrystal was originally described as a "chargetransfer solid," indicating a π−π* charge transfer between benzene and hexafluorobenzene, although spectroscopy only shows small frequency shifts from the two pure solids. 80 The tensor addition method also works well for many other organic crystal structures (Figures 5, 6, Supporting Information, Tables S5 and S6).
3.2.1. Variation of Anisotropy of χ cryst in Polymorphs. The calculated χ cryst for a range of organic crystals (Supporting Information, Tables S5 and S6), using both the periodic DFT and the tensor addition methods, shows that for polymorphs of the same molecule, χ iso cryst remains fairly constant and similar to χ iso mol , but there can be marked differences in Δχ an cryst between the polymorphs (Figures 5 and 6). In many cases, for example, acridine, the variation is completely the result of the differences in the crystal packing of the molecules, whereas in other cases, for example, flufenamic acid, different molecular conformations in each polymorph also plays an important role. Form II of acridine, which has the long axes of the molecules parallel, but the nitrogen positions alternating, has a Δχ an cryst of −150 cgsppm, close to that of χ mol , whereas the other polymorphs 81 The errors in using tensor addition or between using PBE0 or PBE for χ mol are far smaller than the differences in Δχ an cryst between different polymorphs of the same molecule as shown for the highly polymorphic flufenamic acid and acridine in Figure 6. Figure 6 also shows that using the experimental crystal structures and molecular conformations gives very similar χ cryst to those calculated using the lattice energy optimized structure (Supporting Information, Table S4). An exception is form IV (P2 1 2 1 2 1 , Z′ = 3) of acridine, where the PBE-TS optimized structure differs

DISCUSSION
We have shown that a reasonable estimate of the diamagnetic susceptibility tensor of an organic crystal can be calculated using periodic DFT with plane waves and pseudopotentials, using the methodology developed and successfully used for calculating solid-state NMR spectra. 65,66 This has been done using the PBE GGA functional, which is routinely used for modeling organic polymorphs, 82 estimating the differences in the magnetizability tensor χ cryst between polymorphs of a range of organic molecules, including flexible pharmaceuticals. We show that the anisotropy of the tensor, Δχ an cryst , can vary considerably between polymorphs, depending on the difference in the relative orientations of the conjugated molecular fragments, such as aromatic rings.
These calculated crystal diamagnetic susceptibilities are limited in applicability to magnetic fields of a strength where the response, characterized by the χ tensor, comes entirely from the response of the electronic orbital motion. Our calculated χ cryst tensors are also limited by the quality of the charge density: a GGA functional (PBE in this paper) is the highest one on Jacob's ladder of DFT functionals, 83 that could be afforded for covering a wide range of organic crystal structures. Just as an improved treatment of the exchange− correlation with hybrid functionals, such as PBE0, can affect the relative lattice energies of polymorphs, 73 our calculations using PBE0 for the same molecules suggest that this introduces a small change (<5%) in χ cryst . Work on calculating χ mol for small molecules with better descriptions of electron correlation and delocalization, 30,84 and including the effects of molecular vibrations, 85 suggests that the calculated values could change slightly when better theoretical treatments are available in periodic codes. However, these first estimates of the diamagnetic susceptibility tensors of pharmaceutical polymorphs are sufficiently realistic for considering the effect of a magnetic field on crystallization.
The calculation of magnetic susceptibility tensors for organic molecules and crystals is necessary because of the experimental difficulty in measuring this property accurately. Using the Faraday scale, Gouy's scale, or the improved Evans scale 86 and a fluid or powder sample, only the isotropic term of the magnetic susceptibility χ iso is readily accessible and still hard to measure accurately. 31 The assumption of random relative orientations in a gas or liquid is only valid for molecules, that are nearly isotropic in their interactions, and this is not the case for the majority of molecules with a large χ mol and will be further complicated if there is an ensemble of conformations present. It would take even more effort to measure the anisotropy of χ for a single crystal. 29 It is difficult to grow goodquality single crystals of sufficient size for experimental measurement of the χ tensor for many pharmaceutical polymorphs, let alone measuring the size of the crystal. The volume magnetic susceptibility, κ, is usually reported, which becomes challenging for nanoclusters or nanocrystals because of the difficulty in measuring the mass or size accurately. The practice of compacting powder samples to tablets limits the accuracy and reproducibility of all magnetic property measurements. 87 Figure 4. Molecular magnetic susceptibility χ mol of molecular coronene (left) and crystal magnetic susceptibility χ cryst of the β (middle) and γ (right) polymorphs plotted to the same scale and superimposed on the crystal structures and BFDH models of the morphology such that the tensor eigenvalues are correctly aligned. The angle between the normals to the molecular planes is indicated for each polymorph.  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article This paper calculates χ cryst , per mole of molecules, for an infinite crystal. The response of a real crystal to a magnetic field will also depend on its size and shape. Most organic molecules lack the symmetry to define the orientation of the χ tensor although χ max mol will usually be fairly perpendicular to the plane of the aromatic rings if the conformation of the molecule has such a plane (Figure 2). The relationship between the χ cryst tensor and the crystal structure depends on the space group symmetry but is rarely completely defined. 11 The most common space group for organic molecules (P2 1 /c or equivalently P2 1 /n) has one χ cryst tensor eigenvector along the unique crystal axis, but the other eigenvectors can lie anywhere in the perpendicular plane.
The relationship between χ cryst and the crystal morphology is even more uncertain, as the morphology also reflects the intermolecular interactions. In crystals that comprise layers of molecules containing the aromatic rings, the interlayer direction often corresponds to the slowest growth direction, as the interactions between the layers are only van der Waals dispersion, and so, the dominant face comprises aromatic rings. In this case, χ max cryst would be large and tend to be perpendicular to the dominant face in the plate-like crystals. Coronene is unusual in having π···π stacking as the strongest intermolecular interaction, and so, the stacking direction is the fastest growing direction, coinciding with χ max cryst ( Figure 4). In many pharmaceuticals, the crystal structure may involve stronger intermolecular interactions in all three directions and many different angles between the aromatic rings within a triclinic crystal, and so, there may be no (even approximate) symmetry relationship between the orientation of χ cryst and the morphology. Hence, calculations will be essential for interpreting any experimental measurements of χ cryst on real crystals which are finite in size. As most organic crystals are insulators, they are unlikely to have induced surface eddy currents in contrast to metals. However, with diamagnetism being such a weak phenomenon, even miniscule defects of paramagnetic or ferromagnetic nature might significantly alter the entire magnetic response of the whole crystal.
Although a major finding of this paper is that Δχ an cryst can differ significantly between polymorphs, it is not obvious by what mechanism this difference in Δχ an cryst would lead to the selective crystallization of a particular polymorph. For a diamagnetic molecule or crystallite, the anisotropy of χ leads to an energy difference ΔE an between the alignment of χ max cryst or χ min cryst parallel to a uniform magnetic field B However, for our sample of organic molecules with at least one and up to three aromatic rings, Δχ an mol is usually between −50 and −150 cgs-ppm, which even in a 20 T magnetic field would lead to a ΔE an about 10 −4 times smaller than kT at room temperature. Hence, the anisotropy of χ is not sufficient to orient individual pharmaceutical molecules crystallizing in a magnetic field. However, the reorientation energy for a nanocrystallite with a diameter of a few hundred nanometers could be large enough to overcome thermal motion at room temperature in a 1T magnetic field. The degree of alignment of a crystallite would be very dependent on the crystallite size and shape as well as the magnetic field and temperature. 11 A strong magnetic field could affect partially ordered molecular aggregates or prenucleation clusters 88 once they had grown beyond a similar size. Current experimental research seeks to determine at what stage the magnetic field appears to be affecting the crystallization. However, if polymorphs did not differ in χ cryst , no laboratory magnetic field would be able to affect their relative crystallization rates.
Tensor addition of χ mol provides a very good estimate of χ cryst for a wide range of crystals of neutral organic molecules ( Figures 5, 6, Supporting Information, Tables S5 and S6). Compared to χ cryst calculated directly with the more costly periodic DFT method, the rms difference for both the eigenvalues and the isotropic term is only about 5 cgs-ppm, and the maximum difference is around 10 cgs-ppm, when using the same PBE functional for χ mol . This shows that the magnetic response of each individual molecule in the crystal, to a good approximation, remains the same as if it was isolated but in its crystalline conformation. The tensor addition method assumes that there is no change in the charge distribution on crystallization, that is, the delocalization of the electrons, particularly those which play a significant role in Δχ an cryst , is not significant between the molecules in a crystal. This is a significant difference between organic crystals and metals or covalent crystals, where the electrons are delocalized over the crystal. The assumption that the charge distribution of a molecule is the same in isolation as in the crystal is often used for estimating electrostatic contributions to the lattice energy in CSP studies (the Ψ mol approach). 89 The tensor addition approximation may break down when the π electron density is more delocalized over adjacent molecules, for example, when a significant pressure is applied to a crystal. Another approximation inherent in the tensor addition method is the neglect of any local field effect, but as any induced magnetic fields generated by adjacent molecules are extremely small for diamagnetic crystals, this approximation should be valid in all cases. However, with these caveats in mind, the tensor addition method could be used with higher quality molecular susceptibilities than the PBE or PBE0 densities that we have used.
A significant implication of the success of the tensor addition method for organic crystals is that it should be equally applicable for any other form of aggregation of organic molecules with a specified structure, such as a crystallite or nucleating cluster of a finite size. Hence, we have established a method of estimating the difference in the diamagnetic response to a magnetic field between any defined size and shape of crystallite or cluster of molecules. Large, conformationally flexible molecules rarely form isotropic clusters.
The isotropic component χ iso of the diamagnetic susceptibility tensor appears to be very similar between phases of the same molecule as this differs little between polymorphs and the isolated molecule (Supporting Information, Tables S5 and S6). This is consistent with the observation by Pascal and others 28 that the isotropic molecular susceptibilities χ iso can be estimated by additivity methods. It is therefore reasonable that tensor addition can be applied for the full anisotropic tensor for crystals, clusters, or even larger molecules (Supporting Information, Section 3). The conformational flexibility of pharmaceuticals means that the anisotropy of the tensor can vary significantly within any phase. Although there is a general tendency for molecules to adopt extended conformations on crystallization, 90 many pharmaceuticals are not planar within their crystal structures even if their molecular diagrams may seem formally fully conjugated, as seen for flufenamic acid and p-methyl chalcone (Figure 3). The differences in the crystal packing often produce significant The Journal of Physical Chemistry A pubs.acs.org/JPCA Article differences in the tensor components between polymorphs, but this can be estimated quite readily by using tensor addition. One application of the calculation of crystal magnetic susceptibility is to complement CSP studies to determine whether there are unobserved putative structures that have contrasting magnetic properties to, and are also thermodynamically competitive with, the known polymorphs. The calculated magnetic properties could then be a part of the assessment of what experiments might crystallize a targeted putative polymorph for the first time. 20,24 The success of the tensor addition method means that a good estimate of the diamagnetic susceptibility can be readily calculated for many crystal structures which may be thermodynamically competitive in a CSP study. 20,91 An "energy−structure−magnetic property" map, an adaptation of the energy−structure− function maps used to design new functional molecular materials, 92 is shown in Figure 7. This shows the low-energy structures generated in a CrystalPredictor 93 /CrystalOptimizer 94,95 CSP study of the nonsteroidal anti-inflammatory drug diclofenac (details in Supporting Information, Section 8), classified by Δχ an cryst (Figure 7). It is clear that the unobserved structure with the second lowest lattice energy has a greater Δχ an cryst than the known forms. The χ iso values are not significantly different between any of the known and hypothetical structures, so the variation in κ is mainly from the density differences of the CSP structures (Supporting Information, Figure S14).
The variation in Δχ an cryst between the lowest energy unobserved and all the observed polymorphs merely shows that the crystallization of diclofenac might be affected by a strong magnetic field to favor a new form. A more accurate determination of the free energy differences 20 and studies of the crystallization behavior of the known forms would be required to help decide whether it was worth attempting crystallization of diclofenac in a very strong field. However, the development of efficient and reliable computational methods to determine the differences in the χ cryst tensors of pharmaceutical polymorphs is a first necessary step toward the goal of realizing the potential for the polymorphic control of a crystallization process by a magnetic field.

CONCLUSIONS
In this paper, we have developed a methodology for calculating the diamagnetic response of organic crystals with applications to pharmaceutical molecules. The magnitude of the diamagnetic susceptibility tensor, χ iso , in all the phases is dominated by the size of the molecule with significant contributions from the delocalized π orbitals. However, the components of the χ tensor and its anisotropy depend on the relative orientation of the aromatic rings, and therefore, they are affected by conformational changes. The diamagnetic susceptibility tensors for organic crystals χ cryst can be calculated using periodic density functional methods; though currently only with PBE or similar functionals that are widely used in modeling pharmaceutical polymorphs. We have shown that a good estimate of χ cryst can be obtained by the tensor addition of χ mol using an experimental or computer-generated crystal structure for the packing and molecular conformation.
The anisotropy of the diamagnetic susceptibility tensor can differ substantially between polymorphs when there is a significant difference in the relative orientations of the aromatic and other functional groups that can maintain an induced ring current in a magnetic field. Thus, it is now possible to estimate whether the magnetic properties of pharmaceutical polymorphs differ so much that their crystallization could be differentially affected if crystallized in a strong magnetic field. When we understand the mechanism by which a magnetic field can affect crystallization, the calculation of energy−structure− magnetic property maps may enable the targeted discovery of new polymorphs.
Illustration of induced ring currents in benzene; derivation of tensor addition for estimating the diamagnetic susceptibility tensor of an organic crystal; group additivity; experimental and optimized crystal structures considered in the study; tables of calculated molecular and crystal magnetic susceptibilities; crystal structure prediction study of diclofenac; and notes on different systems of units for magnetic properties (PDF)