Analysis of the conformational profiles of fenamates shows route towards novel, higher accuracy, force-fields for pharmaceuticals †

In traditional molecular mechanics force fields, intramolecular non-bonded interactions are modelled as intermolecular interactions, and the form of the torsion potential is based on the conformational profiles of small organic molecules. We investigate how a separate model for the intramolecular forces in pharmaceuticals could be more realistic by analysing the low barrier to rotation of the phenyl ring in the fenamates (substituted N -phenyl-aminobenzoic acids), that results in a wide range of observed angles in the numerous fenamate crystal structures. Although the conformational energy changes by significantly less than 10 kJ mol (cid:2) 1 for a complete rotation of the phenyl ring for fenamic acid, the barrier is only small because of small correlated changes in the other bond and torsion angles. The maxima for conformations where the two aromatic rings approach coplanarity arise from steric repulsion, but the maxima when the two rings are approximately perpendicular arise from a combination of an electronic eﬀect and intramolecular dispersion. Representing the ab initio conformational energy profiles as a cosine series alone is ineﬀective; however, combining a cos2 x term to represent the electronic barrier with an intramolecular atom–atom exp-6 term for all atom pairs separated by three or more bonds (1–4 interactions) provides a very effective representation. Thus we propose a new, physically motivated, generic analytical model of conformational energy, which could be combined with an intermolecular model to form more accurate force-fields for modelling the condensed phases of pharmaceutical-like organic molecules.


Introduction
Huge advances in our understanding of biomolecular behavior have been made using molecular mechanics force-fields such as AMBER, 1 GROMOS, 2 CHARMM 3 and OPLS-AA, 4 repaying the immense effort that has gone into parameterizing these force-fields.However, even for protein and nucleic acid polymers, which are well suited to assumptions of transferability of parameters for specific residues or bases, the search for increasing accuracy for more demanding energetic predictions continues, 5 with more complex forms such as AMOEBA 6 or the addition of numerical grids to model dihedral cross-terms more accurately. 7The traditional force-field includes explicit periodic torsional potentials, and applies the non-bonded terms to all intramolecular atom pairs separated by three or more covalent bonds (1-4 interactions and above).Thus, traditional force-fields are built on the reasonable assumption that, in biomolecules, the intramolecular non-bonded interactions are the same as intermolecular non-bonded interactions.The explicit torsional potential helps the force-field to model the low energy conformations of the peptide interacting with itself, ligands and solvents.The importance of an accurate balance of inter and intramolecular forces for challenging applications involving molecular recognition of pharmaceuticals, such as computer-aided drug design, cannot be overemphasized. 8owever, it is clear that the traditional force-field model is limited in accuracy, 9,10 for example, the electrostatic models that are successfully used in modelling intermolecular forces are not valid at some 1-4 distances because of charge cloud overlap 11 leading to significant penetration effects. 12,13The customary approximation in many force-fields of just halving the atom-atom non-bonded intermolecular functions for the 1-4 intramolecular energy terms cannot be very accurate.
Applying traditional force-fields to pharmaceutical-like molecules, with flexible bonds linking multiple functional groups, often relies on poorly justified transferability assumptions for lack of data for empirical fitting 14,15 and hence are often unsuccessful for simulating the properties of pharmaceuticals materials. 16or example, there have been no successes based on the use of force-fields for final lattice energy evaluations in the blind tests of organic crystal structure prediction (CSP). 17This can be due to the force-field giving a qualitatively wrong conformation even for the isolated molecule, as in the case of aspirin. 18Lattice energy minimizations with force-fields can change the molecular conformation within the crystal so much that the relative positions of the functional groups leads to a qualitatively different crystal structure. 19In other cases, the structures may be reproduced adequately but the failure to rank the energies properly has been traced to the use of the same charges and van der Waals interactions for the intermolecular and intramolecular forces. 20onsequently, the most successful approaches to CSP rely on expensive electronic structure calculations of the molecular conformational energy.In one CSP approach, a specific tailormade force-field is parameterized for the molecule from dispersion corrected density functional (DFT-D) calculations and used to generate crystal structures, but then the most promising crystal structures require refinement by periodic DFT-D lattice energy minimizations. 21In another CSP approach, ab initio calculations on a single molecule either in isolation 22,23 or in a polarizable continuum, 24 are used to evaluate the energy penalty for changes in conformation (DE intra ), and provide the distributed multipole representation of the molecular charge density 11,25,26 used for the electrostatic contribution to the intermolecular lattice energy.Thus, the crystal structure prediction methodologies that have advanced to aid solid form screening in pharmaceutical development [27][28][29][30][31] require a very large number of electronic structure calculations to define the conformational potential energy surface of the molecule.Hence, the approach needed to give the relative energies of different possible crystal structures is far too computationally demanding to be used in Molecular Dynamics simulations.Such simulations are highly desirable for calculating the relative free energies of organic polymorphs [32][33][34] as ambient conditions are rather too close to the melting temperatures of organic solids to rely on the harmonic approximation.Therefore, we need to model the molecular flexibility of typical pharmaceutical molecules by a force-field that accurately reproduces the relative energies of the known and thermodynamically competitive crystal structures and yet can be evaluated sufficiently quickly for realistic Molecular Dynamics simulations.Such a force-field would be used for assessing crystal stability at ambient temperatures, calculating the relative free energies of known and potential polymorphs, and for simulating nucleation and other molecular recognition processes of the molecules. 35he success of CSP studies suggest that it is worth investigating decoupling the models for the intermolecular forces from those for the intramolecular forces (i.e.conformational profiles) as a route to more accurate force-fields for pharmaceutical molecules.The intermolecular forces could be modelled by anisotropic atom-atom potentials, as currently used in CSP studies, so the requirement is for an analytical form for the intramolecular energy changes (DE intra ).Hence, for a preliminary investigation into how we could model conformational energy differences of pharmaceuticals, we investigate a single torsion angle that both exemplifies the challenges of conformational flexibility in crystal structure prediction, and has long been seen as a key determinant in the pharmacological activity of a family of analgesics. 36The fenamates (Fig. 1) are so prone to conformational polymorphism, 37 that the fenamate unit has been termed a polymorphophore. 38The non-steroidal antiinflammatory drug flufenamic acid (FFA) holds the current record for having crystal structures determined for nine polymorphs, 39 and another, tolfenamic acid (TA), has at least five polymorphs. 38In the majority of fenamate crystal structures, the intramolecular and intermolecular (carboxylic acid dimer) hydrogen bonds are preserved, and it is the torsion angle (x) defining the orientation of the (substituted) phenyl ''paddle wheels'' that varies, leading to the marked differences in the crystal packing.Crystal structure prediction studies on fenamic (FA) and tolfenamic acid (TA) 40 show that it is the subtle compromise between the packing of the substituted phenyl rings and the small conformational energy penalty that leads to the polymorphism of TA and monomorphism of FA.Mefenamic acid (MA) differs from TA only by a chloro/methyl exchange at a position (R 2 Fig. 1) that would not be expected to affect the torsional profile, and the chloro/methyl substitution can lead to isostructural crystal structures, [41][42][43][44][45] yet the crystal energy landscape of MA is distinct from that of TA. 46 Many studies have emphasized the difficulty in evaluating the relative energies of fenamate polymorphs, [47][48][49] or controlling the polymorphic outcome by varying the crystallization conditions. 50A study of the distribution of the fenamate-like torsion angle x in the organic crystal structures within the Cambridge Structural Database 51 shows that a wide range of angles can be adopted, but these correlate with very low conformational energies. 40his study analyzes the one dimensional conformational profile for the fenamate torsion, x, contrasting FA, TA, ClFA, and MA (Fig. 1), which differ in substituents that would be expected to change the conformational profile (R 1 ) and those sufficiently distant (R 2 ) to be expected to have little effect.To find an analytical model that can reproduce these lowenergy conformational profiles accurately requires a physically justified functional form for effective parameterization.The cause of rotational barriers has long been controversial, 52 with the ethane rotation barrier still generating discussion as to whether the origin is steric, hyperconjugation [53][54][55] or an electrostatic effect. 56The quantitative distinction between electronic effects from changes in the molecular orbitals, as opposed to steric ''non-bonded'' effects, is dependent on the precise definition and type of charge density calculation.We investigate the qualitative issue for the fenamate molecules by constructing model molecules with minimal steric effects so that conformational profile is dominated by the change in the electronic effects, loosely termed conjugation.This type of effect would be expected to be represented by a few terms of the traditional explicit torsional term of the general form: where g defines the phase shift, and n the periodicity.For FA, x defined in Fig. 1 requires that g = 0 and symmetry dictates that n is an even integer.For substituted fenamates, odd values of n could contribute, though the conformational profile should be symmetric about x = 0.If the origin of the conformational profile is predominantly steric, caused by the varying repulsion between the overlapping charge distributions of 1-4 atoms with x, then it should be well represented by non-bonded atom-atom interactions.An exp-6 atom-atom model would be expected to give a better representation than a Lennard-Jones (R À12 ) model of the variation of the repulsion with distance, given the success of the overlap model in parameterizing intermolecular repulsion potentials. 57,58hus a crude starting point for our investigation of the ''nonbonded'' contribution to the torsion potential is the exp-6 atom-atom model potential with a parameterization that has been developed for modelling the intermolecular forces between organic molecules in crystals: [59][60][61][62] EðxÞ ¼ X where atoms i and k of atomic types i and k are separated by intramolecular distances R ik , calculated from the molecular conformation with torsion angle x.Not including an explicit electrostatic term in the intramolecular potential considerably simplifies the implementation and extension to larger molecules.][61][62] The aim of this study is to establish the physical basis for the variation in conformational energy of the fenamates to determine what might be a reasonable analytical model.This tests whether the understanding of torsional potentials that has been developed for small molecules, such as ethane, need amending for larger organic systems.The work concentrates on fenamic acid (FA) and tolfenamic acid (TA) but the analysis is extended to related molecules to assess generality.If we can find an appropriate analytical functional form for the conformational energy of a family of molecules, then the combination of separate intermolecular and intramolecular potentials would provide more accurate analytical force-fields for pharmaceutical molecules.

Ab initio conformational energy profiles
The conformational profiles used throughout this study, unless otherwise specified, were relaxed conformational energy scans at PBE0/6-31+G(d) level of theory, carried out using GAUSSIAN03. 63e found that the results of TA were sensitive to the atoms used to define the torsion angle, (i.e.Fig. S1 of the ESI † shows that defining the torsion by H 6 -N 1 -C 8 -C 9 or H 6 -N 1 -C 8 -C 13 could double the height of the maximum at x = 01, and even using C 7 -N 1 -C 8 -C 13 could show differences at high energies as x approached 1801).This common observation that torsional scans depend on which 4-atom set is used to describe a torsion about a rotatable bond complicates the analysis of rotamer distributions. 64The profiles were also dependent on the starting points.Hence, to determine the starting geometry, we performed a full optimization near each symmetry independent potential minimum.The grid consisted of the optimized structures and points x = 5n1, with the scans going from the potential minima.Highly repulsive points as x approached 1801 for the substituted fenamates were omitted.

Investigation of electronic versus steric effects
To attempt to separate out the steric effects from electronic effects, a series of model molecules where the steric effects had been minimized were studied.When the benzoic acid group was replaced by a series of smaller molecular fragments, such as hydrogen atom, the three bonds to the nitrogen atoms were constrained to be coplanar by fixing an improper torsion angle relating the nitrogen to the three bonded atoms.This prevents the pyramidalization at the nitrogen in the torsional potential of phenylamine 56 and the major rearrangement of the second hydrogen that occurs in a relaxed scan of H 7 -N 1 -C 8 -C 9 .

Atom-atom modelling of torsion potentials
A starting point for considering intramolecular steric interactions is the exp-6 atom-atom intermolecular parameters derived by Gavezzotti by fitting to crystal structures and heats of sublimation of hydrocarbons, oxahydrocarbons, azahydrocarbons, chlorohydrocarbons and nitro compounds. 59This provides the parameters for all intramolecular interactions involving C, N, O, Cl and H.The same parameters are used for all C and H atoms, whether aromatic or in the methyl or carboxylic acid groups, hence the atom typing is crude compared with current force-fields e.g.Sybyl typing. 65There is a polar hydrogen type, HB, which we use for both polar hydrogens (H 1 and H 6 in Fig. 1) available from the extension of the exp-6 parameterization to hydrogen-bonded crystals. 60The HBÁ Á ÁO/N exp-6 potentials have particularly deep wells as they have absorbed the electrostatic effects in intermolecular hydrogen bonding. 60Since the intermolecular parameter set does not have parameters for the HÁ Á ÁHB and CÁ Á ÁHB interactions, these were fitted in this study, using the parameters for HÁ Á ÁH and CÁ Á ÁH as a crude starting point.The intermolecular exp-6 parameters are given in Table S1 of the ESI, † with other types available from the scheme in Table S2 (ESI †).
The atom-atom interactions are summed over all 1-4 and higher bond-paths in the entire molecule (i.e. this explicit intramolecular force-field does not distinguish 1-4 from the other intramolecular interactions, as the 1-4 interactions are not always the shortest intramolecular atom-atom distances Tables S3 and S4, ESI †).The maximum bond-path is 1-11 for all three fenamates, and most pharmaceutical molecules are sufficiently small that there is no need to define a summation limit in terms of intramolecular distance or bond-path length.The use of 1-4 distances as the shortest intramolecular interactions included in the atomatom summation is traditional, although we note that 1-3 interactions are used elsewhere, e.g. in the CSP code in GRACE where they have sometimes been found to be problematic, including in the case of a bulky side group attached to an aromatic ring which required specific scaling down. 66he atom-atom formulation results in many virtually constant terms, such as the HÁ Á ÁH, CÁ Á ÁC and HÁ Á ÁC contributions from within the aromatic rings.These terms contribute to the baseline energy, E base , defined as the minimum energy found in the scan with a specific parameterized model.

Fitting analytical models to the conformational profiles
For this preliminary investigation of suitable functional forms, we have not applied any weighting to the conformational profile beyond restricting the points to conformational energies below 10 kJ mol À1 , i.e. not seeking to accurately represent the steric barrier above 1451 for TA and MA.This gives N p = 37 energy data points for the FA fit and N p = 30 for TA and MA.
Since there is considerable correlation between the atomatom coefficients, particularly the two repulsion parameters (A ik and B ik ), 67,68 we seek to rescale selected repulsion and dispersion coefficients, giving a linear model that can be combined with an appropriate cosine term and an approximation to the baseline constant, c B E base , to give A generalizable method of deriving analytical models for conformational profiles was developed during this work (Fig. S4, ESI †), which includes the FORTRAN code and NAG 69 library routines for systematically comparing the ability of various selections of the linear parameters (a, b ik , g ik , c) to represent the ab initio data by least squares using a general linear regression model. 70 Results

Conformational energy profile of the fenamates
The conformational profiles for four of the fenamates (Fig. 2) show that there are two distinct minima, which are only close in energy for the symmetric FA and R 2 -substituted ClFA.There is one potential maximum that varies a little between the molecules for the planar conformation (as drawn in Fig. 1), another maximum when the aromatic rings are approximately perpendicular, and a third where there is a significant steric clash for the R 1 -substituted fenamates TA and MA as the other planar conformation is approached.An analysis of the observed values of this torsion angle in crystal structures containing the fenamate fragment (Fig. 2 of ref. 40) shows that the observed angles are clustered around the two minima, consistent with the expectation that most molecules adopt low energy conformations in crystal structures. 37Hence, the two low energy barriers (around 5-9 kJ mol À1 ) clearly have a major effect on the crystal packing and are large compared with most measured polymorphic energy differences, including those of TA which cover less than 2 kJ mol À1 . 40esting the sensitivity of this conformational profile to the choice of ab initio method (Fig. 3) shows that even obtaining a conformational profile in qualitative agreement with that derived from experimental crystal structures is sensitive to method.The HF scan has only one minimum, at a conformation that is not observed in the crystal structures of fenamates; a CSP study based on this conformational profile would generate qualitatively incorrect crystal structures.Only the ab initio methods that include some description of electron correlation produce a maximum at around 901.Nevertheless, there is fair agreement in the conformational barriers for these methods.However, it is notable that for TA and MA evaluating the energy using PBE0 geometries at the MP2 level swaps the relative energy of the two minima, with the PBE0 calculations being in better agreement with the analysis of the crystalline conformations of fenamate-like fragments with a substituent at C 13 .Evaluating the conformational profile within a polarizable continuum model 71 with e = 3, a typical dielectric constant of organic crystals, 72 showed a reduction in energy penalty around This journal is © the Owner Societies 2015 the energy maxima (x B 901) for all the fenamates (Fig. 3).The application of this polarizable continuum model (PCM) has been shown to improve the relative energy ranking of some conformational polymorphs. 72The results in Fig. 3 confirm that we cannot obtain definitive ab initio conformational energy scans, but the scans in Fig. 2 are adequate for the purposes of this study.
To establish the importance of changes in the other torsion angles, bond angles and bond lengths during the relaxed scans in Fig. 2, the conformation scans were repeated, starting from the fully optimized structure at the PBE0/6-31+G(d) level of theory and only allowing the angle x to change.Thus each conformation is identical when a single point energy is evaluated for each ab initio method.The energy differences (Fig. 4) are very marked for all the fenamates, even in the lowest energy regions.The extent to which the positions of the other atoms relax to lower the conformational energy barrier is very marked around the energy maxima (Fig. 5), even though these maxima correspond to an energy of less than 7 kJ mol À1 in the relaxed scans (Fig. 4).Even for FA, there is a significant change in the internal hydrogen bond with the reorientation of the benzoic acid around all maxima.Both methyl torsions play a role in reducing the energy around the 901 maximum for TA (Fig. 5).Indeed, even for FA, the changes in the other conformational variables produce the slight asymmetry between the two minima (Fig. S2 of the ESI †).Hence, the contrast between the relaxed and rigid (Fig. 4) scans and the corresponding conformations (Fig. 5) confirm that the changes in the other conformational variables play a major role in lowering the conformational barrier over a wide range of x angles, including those sampled within the crystal structures. 40

Splitting electronic from steric contributions to the barrier to rotation
To establish the importance of the intramolecular steric clashes, the rigid and relaxed scans were repeated with model molecules in which the benzoic acid group was replaced with a smaller fragment and the bonds around the nitrogen constrained to be planar so as to avoid the pyramidalization of the amine (Fig. 6).The potential energy scan has a very large maximum at 901 for planar-N-constrained phenylamine (PA), approximately 6 times higher than the barrier in the fenamates.(Note that Fig. 6 has the same scale as the rigid scans in Fig. 4, covering a larger energy range than the relaxed scans in Fig. 2 and 3).Replacing one constrained hydrogen with a methyl (Fig. 6b) produces a slightly larger barrier but very little asymmetry, strongly suggesting that this is an electronic effect of conjugation between the lone pair on the nitrogen and the aromatic ring.Adding a double bond to PA reduces the barrier to rotation by 11.10 kJ mol À1 (Fig. 6c), implying that the conjugation with the benzoic acid ring of the fenamates will have contributed significantly to reducing the electronic barrier.Adding a carboxylic acid that forms an intramolecular hydrogen bond to the N-H group further reduces this barrier by almost 7 kJ mol À1 (Fig. 6d).The intramolecular hydrogen bond in this model molecule (Fig. 6d) varies in length from 1.94 to 1.92 Å as x changes from 0 to 1801, in comparison with the FA hydrogen bond varying from 1.87 to 1.84 Å, and so we may infer that the intramolecular hydrogen bonding in the fenamates will similarly reduce the electronic barrier.In contrast, substituting Cl and CH 3 at the meta position of PA (Fig. 6e and f) shows only a small change to the barrier height, h, with 3-chloroaniline (h = 37.0 kJ mol À1 ) and 3-methylaniline (h = 35.6 kJ mol À1 ) having only slightly larger barriers than PA (h = 35.2kJ mol À1 ).
The scans in Fig. 6a, b, e and f clearly have no maxima at 0 or 1801 confirming that the steric clash between the aromatic C 6 -H and the C 9 -H or C 13 -R 1 groups of the phenyl ring are responsible for these maxima.These curves are very well reproduced by (h/2)(1 À cos(2x)) where h is the potential maximum.As the nitrogen substituents get larger (Fig. 6c  and d), there are signs of additional steric effects at 0 and 1801 and a larger difference in the barrier height at 901 between MP2 and HF calculations, suggesting there is more change in intramolecular dispersion.The difference between a rigid and relaxed scan is small (Fig. 6a and b), and the difference in the curves with type of calculation are relatively minor compared with the qualitative difference between the HF and correlated methods for the fenamates (Fig. 3).The overriding conclusion from contrasting the conformational scans of model molecules with minimal steric effects (Fig. 6) with those of the fenamates (Fig. 2) is that there is an electronic contribution to the torsional barrier at x = 901, which can be represented by a (h/2)(1 À cos(2x)) term.This can be rationalised as resulting from the changing conjugation between the nitrogen lone pair and the phenyl ring.(The simple idea of conjugation, that a fenamate with x = 0 is stabilized by having a p orbital delocalized over both rings is inappropriate as this conformation is not planar (Fig. 5).)The changes in the phenyl molecular orbitals with conformation are similar for PA and the fenamates. 73owever, the electronic effects are not solely responsible for the maxima at x B 901.The intramolecular atom-atom dispersion contribution (ÀC ik /R ik 6 ) from eqn (2) using the intermolecular parameters produces a significant maximum in this region (ESI, † Section 1.3).The repulsion component (ÀB ik R ik ) gives maxima at 0 and 1801 and minima at 901, consistent with the expectation that these two maxima occur because of steric clashes.Thus, the analysis of the torsional potentials of the fenamates reveals that there is an electronic effect from the change in the orbital interactions such as ''conjugation'' between the aromatic rings and the nitrogen lone pair  This journal is © the Owner Societies 2015 destabilizing the non-planar conformations and steric effect from the variation in overlap of non-bonded atoms.For these larger molecules, in contrast to the well-studied small molecule torsional potentials (e.g. of ethane), intramolecular dispersion and small changes in the other conformational variables also make a very significant contribution to the torsional profile.
where N k and N p are the number of fitted coefficients and data points respectively, must include cos 4x to have the correct number of minima for the fenamates.This term alone gives a poor position of the minima for FA and is qualitatively wrong for TA (Fig. 7).A least squares fit including the lower cosines (N k = 5, eqn ( 4)) gives a qualitatively reasonable representation (Fig. 7), but further improvement is slowly converging (Fig. S5 in the ESI †).This demonstrates that the cosine series is a fitting exercise, not reflecting the physics.It is effectively modelling the relaxed scans by a functional form that assumes the scan is rigid (i.e.only the torsion angle changes), whereas there are (Fig. 4) significant relaxation effects that reduce the conformational barrier.Poor convergence of the cosine model expansion has also been reported for a biphenyl torsion within a dye which required 7 terms, 74 while polynorbornene 75 required 6 and 15 terms for the meso and racemic dimer respectively.3.3.2.Rescaling the repulsion model.It is possible to get an excellent fit to the torsion potentials by summing the exp-6 potential (eqn ( 2)) over all 1-4 and higher intramolecular atomatom distances (ESI, † Section 3.2), provided that a few of the intermolecular repulsion coefficients parameters A ik are rescaled by a factor b ik .Rescaling just two repulsion contributions for FA (CÁ Á ÁH and HÁ Á ÁH) and five for TA (CÁ Á ÁH, CÁ Á ÁHB, CÁ Á ÁN, HÁ Á ÁH and HÁ Á ÁHB) produces a model that reproduces the torsional profile well (Fig. S6 in the ESI †).However, some of the fitted rescaling parameters were negative (b ik o 0), which implies that the exponential steric repulsion had become attractive.Thus, again, this appears to be an unphysical fitting exercise.
3.3.3.Combined physical model.If we assume a model that describes both the electronic effects and allows rescaling of the atom-atom interactions: then there are a huge number of ways of finding a satisfactory fitting of the data (ESI, † Tables S8 and S9).However, only when b and g are positive do the contributions retain the repulsive and attractive (dispersion) nature respectively, and a negative value of a is required to give a maximum around x = 901 corresponding to ''conjugation'' (cf.Fig. 6).Fitting a and rescaling only a few atom-atom interactions gives a qualitatively accurate fit (Fig. 8).It is not surprising that the parameters involving H and HB require significant rescaling as the intermolecular values were not well defined. 59Although virtually perfect fits can be obtained (Table S8 in ESI †), the variation in the fitted parameters is significant, which is not surprising given the exponential sensitivity of the repulsion to changes in atom-atom distances.These changes can be substantial, for example the two minima in the FA scan at x = 38.941and 144.711 correspond to conformations that differ by 0.61 Å in the 1-4 distance between the amide proton and C 13 (or C 9 ) although the minima only differ in energy by 1.88 Â 10 À3 kJ mol À1 .Even the C 8 -C 9 and C 8 -C 13 aromatic bond lengths differ by AE0.0037 Å for x = 0 or 1801, but only by AE0.0004 Å for x = 80 or 1001, with larger changes in the bonded hydrogen positions.In contrast to FA, TA gives a qualitatively acceptable fit (Table S9 in ESI †) only when at least three types of atom-atom parameters are rescaled (N k = 8), including CÁ Á ÁN, as shown in Fig. 8. TA differs from FA in having many more intramolecular distances that change significantly with x including some within the same aromatic ring, such as methyl-chloro interactions (Table S5 in ESI †).Thus, we are able to obtain a variety of analytical models of the form of eqn ( 5) that can reproduce the torsional profiles of FA and TA with a high degree of accuracy, despite the significant variation in many atom-atom distances during the relaxed torsional scans.The ease with which this physically justified model could reproduce the conformational dependence of FA and TA shows that the approach is promising.

Transferability
We can further investigate the physical applicability of the atomatom plus electronic functional form (eqn (5)) by testing whether the models can describe the torsional potentials of related molecules.The difference between the energy scans when the substituents are far from the varying torsion angle, for example the change of H to Cl atom i.e. from FA to ClFA or between TA and MA (Fig. 2), should not change the main steric interactions, and will only have a small effect on the electronic term (Fig. 6e and f).
Transferring a set of parameters fitted to FA and adding Cl parameters, does indeed (Fig. 9) give some of the asymmetry in the well depths seen for ClFA, and with fitting the electronic term the higher central barrier (Fig. 2) is also reproduced.Similarly, a set of parameters fitted to TA can reproduce the lower barrier at x = 901 and higher barrier at x = 01 in MA (Fig. 9), despite the conformational relaxation of the methylmethyl interaction being somewhat different from that of the methyl-chloro geometry (Fig. 5).Further examples of the transferability of the parameters are given in the ESI † Section 4.

Physical origins of torsional potentials
By analyzing the low energy torsional barrier in the fenamates (Fig. 2), it is clear that larger organic molecules retain the contributions identified for small model molecules, such as ethane, in that there is both an electronic and a steric component.However, as the molecules become larger, the effect of small correlated changes in the other bond angles and the dispersion contribution become very significant.The difference between the low energy torsional barrier and the one calculated holding other conformational variables constant is surprisingly large (Fig. 4).This means that attempts to represent the conformational barrier by a cosine series, ignoring the position of the other atoms beyond 1-4, degenerates into an ineffective fitting exercise (Fig. 7 and as shown in ref. 74 and 75).An atom-atom formulation can directly reflect the geometric changes in the relaxation.
We could not use a definitive ab initio torsional potential for each fenamate, so used a set of consistent, qualitatively realistic, potentials, because of the variation in the relative energies within the affordable methods (Fig. 3).As molecules increase in size, there is an increasing contribution to the conformational profile from the intramolecular equivalent of the intermolecular dispersion.Since dispersion is an electron correlation effect, this makes converging to an accurate ab initio torsional profile very demanding of the type and quality of electronic structure calculation 76,77 because of the importance of electron correlation and intramolecular basis set superposition error. 78The Tyr-Gly peptide conformational minima, 79 alanine dipeptide fc energy maps 7 and the barriers to torsional rotation in p-conjugated polymers 80 have also been shown to vary significantly with choice of post-Hartree-Fock theoretical approach.Electron correlation plays a critical role in what we can qualitatively recognize as through space intramolecular dispersion effects and changing conjugation of the molecular orbitals.This is in addition to the variation in the repulsion and electrostatic interactions within the molecule that would be expected from the sensitivity of intermolecular interactions to the ab initio method. 81Separating the ''through space'' intramolecular dispersion from the other electron correlation effects that contribute to the (h/2)(1 À cos(2x)) electronic barrier from ''conjugation'' or delocalization between the two aromatic rings is probably not quantitatively meaningful when using a quantum mechanical method that approaches the quantitative accuracy needed.The challenge of extending the reliability and accuracy of electronic structure methods to larger molecules, which are more typical of pharmaceuticals and realistic biological molecules in isolation or condensed phases, is the subject of much active research. 76,82This study emphasizes the risk in using affordable but approximate electronic structure methods to provide a large data set of conformational energies for fitting, as HF methods would provide qualitatively misleading results for the fenamates.
The most generalizable analytical models for torsion potentials will represent the physical origins of the contributions.For the fenamates, the electronic term is well represented by a cos(2x) contribution, and an atom-atom model is appropriate for representing the steric and dispersion contributions and automatically includes the effect of relaxation of the rest of the molecule.However, the simple exp-6 model used here is only a first approximation for the intermolecular forces 83,84 and could not be expected to translate accurately to the shorter intramolecular distances that vary with the torsion angle x.Pairs of atoms of atomic types that would rarely, if ever, be found in van der Waals intermolecular contact can be at very short and varying 1-4 distances within a molecule.Modelling conformational energies using a simple atom-atom exp-6 form is effective, but the interactions involving hydrogen atoms, and the methyl carbon nitrogen interaction for TA and MA, were described by significantly different parameters from those empirically fitted for modelling intermolecular forces.These atomic types are involved in some of the atom-atom distances that change most with x.The original intermolecular parameters appear to be able to capture the smaller changes from molecular relaxation adequately.The net result is that it is possible to obtain an analytical expression in the form of eqn ( 5) that can model the conformational curves of the individual fenamates extremely well (Fig. 8) and could be transferable (Fig. 9).
The range of the sets of atom-atom parameters that can reproduce the limited ab initio data on the torsional profiles of the fenamates shows that much more extensive sets of ab initio calculations with greater variations in the other degrees of freedom would be required for fitting eqn (5) to provide a robust analytical model.It would be helpful to have more stringent constraints on what would constitute a physically reasonable range of parameter values based on more careful characterization of intramolecular ''steric'' interactions.Nonetheless, the functional form appears promising for the ability to represent the complex interactions that lead to the low energy torsional potentials in fenamates.

Towards more accurate force-fields for pharmaceutical molecules
This approach to modelling conformational energies of the fenamates could be extended to many pharmaceuticals that comprise approximately rigid molecular fragments joined by flexible linkages that allow the molecule to adopt a wide range of conformations.We can envisage a general scheme for determining such potentials for a given molecule following a crystal structure prediction (CSP) study 22,85 which involves the calculation and storage of a large database of ab initio conformational energies and forces for the pharmaceutic al molecule. 27,28,30This database will cover most of the range of conformations that are likely to be sampled in a Molecular Dynamics study of the molecules in condensed phases, with a strong bias towards the conformations that occur in low energy crystal structures, including known and possible polymorphs. 86This database could be used to parameterize the analytical conformational energy model, adapting the fitting routines written for this study.
Using a physically motivated analytical functional form ensures that the extrapolation to other high energy conformations will be realistic.Building the analytical force-field in conjunction with a CSP study for a specific molecule would have the advantage that the analytical intramolecular forcefield could be validated by ensuring that it reproduced the crystal energy landscape, i.e. that the energies of different packing, hydrogen bonding and stacking modes were correctly balanced with the accompanying conformational changes.The use of a physically based functional form is more conceptually pleasing in its generality than fitting the ab initio data by a This journal is © the Owner Societies 2015 molecule-specific force-field defined relative to the lowest energy conformation, 87 or constructing a neural network potential. 88he application of separate analytical potentials for both intra-and intermolecular terms will require adaptation of molecular modelling codes; however, the coding for the energies, forces, and second derivatives of the proposed intramolecular force-field (eqn ( 5)) is already in most codes.The calculation of inter-and intramolecular terms would need to be separated in programs that use traditional force-fields such as DL_POLY 32,89 for Molecular Dynamics simulations.However, the greater accuracy of the intramolecular forces is most needed in combination with the more accurate anisotropic atom-atom intermolecular potentials for organic molecules.The analytical intramolecular potential models could be incorporated in the rigid-molecule codes DMACRYS 25 and DL_MULTI 90 which use distributedmultipole electrostatic models for static lattice and Molecular Dynamics modelling of organic crystals respectively.The conformation dependence of the distributed multipoles would need to be considered, but new methods of partitioning the charge density 91 may reduce the conformation dependence, or it could be represented by an analytical model 92 or interpolation scheme. 93lthough this change in approach to pharmaceutical forcefields is envisaging a specific model fitted for each molecule, the physical basis of the current model (eqn ( 5)) and the results in Fig. 9 and Fig. S7 (ESI †) suggest that a reasonably transferable set of atom-atom intramolecular exp-6 potentials could be fitted for families of molecules.Deriving a transferable model would require a very large dataset of ab initio conformational profiles of many molecules calculated at an appropriate accuracy.The transferability of the electronic term (a coefficients) would also need investigating.However, using separate atomatom models for the forces within and between molecules could provide a significant improvement in accuracy on current force-fields, whilst maintaining the advantages of transferability for families of flexible pharmaceuticals.

Conclusions
The torsional potentials of organic molecules not only include short-range electronic ''conjugation'' effects and steric interactions, but also have a significant contribution from the intramolecular dispersion and small concerted changes in other conformational variables.This has two important consequences.Firstly the ab initio determination of organic molecule conformations is very sensitive to the treatment of electron correlation.Secondly, it is not possible to view a torsion as being simply transferable (i.e. the potential is defined by just the atomic types involved in 1-4 interactions) or expect it to be effectively modelled as a cosine series However, we have shown that an appropriate cosine term for the short-range electronic effects plus an isotropic atom-tom exp-6 intramolecular potential can model the conformational profiles of the fenamates well, provided that some of the coefficients are fitted to ab initio torsion potentials.It is clear that discarding the assumption that the same atom-atom models can be used for inter-and intramolecular forces is a route forward to more accurate force-fields.

Fig. 2
Fig. 2 Relaxed conformational scans at PBE0/6-31+G(d) level of theory for the fenamates.The minima were at x = 38.941and 144.711 for FA, 40.631 and 111.861 for TA, 44.081 and 110.481 for MA, and 35.861 and 148.381 for ClFA.

Fig. 4
Fig. 4 Comparison of rigid (dotted line) scans of FA (left) and TA (right) at HF, MP2 and PBE0 methods with 6-31+G(d) basis set as a function of torsional angle x.The relaxed PBE0 scans from Fig. 3 are shown for comparison as a solid line.

Fig. 6 3 . 3 .
Fig. 6 The relaxed (solid lines) torsional scan of planar-N-constrained models for the phenyl rotation, where the benzoic acid group of fenamic acid has been replaced by (a) hydrogen atom (PA), (b) methyl, (c) vinyl and (d) prop-2-enoic acid, and the hydrogen in the meta position of PA has been replaced by (e) chlorine, and (f) methyl using HF, PBE0 and MP2 methods with the 6-31+G(d) basis set.Plots of (h/2)(1 À cos(2x)) where h is the height of the barrier of relaxed PBE0/6-31+G(d) scans are shown in green.For (c)-(f) the PBE0/6-31+G(d) relaxed scans of PA from (a) are shown in grey for comparison.In (a) and (b) a rigid scan at the PBE0/6-31+G(d) level of theory is shown by a dotted line.

Fig. 7
Fig.7Comparison of the ab initio relaxed scan at PBE0/6-31+G(d) level of theory with linearly fitted least square cosine series model for (a) FA and (b) TA, with N k = 5 (red lines), the optimal cos 4x terms (blue lines), and N k = 21 (green lines, for which the quality of fit s intra = 0.17 kJ mol À1 for FA and coincidentally for TA).

Fig. 8 Fig. 9
Fig. 8 Comparison of the relative energies from ab initio calculations (solid black lines) of (a) FA and (b) TA with selected physical models.For FA, the selected N k = 6, s intra = 0.26 kJ mol À1 model has b C-HB = 8.01, g C-HB = 16.54,b H-H = 10.74,g H-H = 58, a = À3.39, and c = 123.69kJ mol À1 , whilst the N k = 6, s intra = 0.28 kJ mol À1 model has b C-H = 0.28, g C-H = 7.83, b H-HB = 11.09,g H-HB = 34.79, a = À1.17, and c = 75.59kJ mol À1 .For TA, the parameters of the selected fits N k = 8, s intra = 0.15 kJ mol À1 and N k = 10, s intra = 0.08 kJ mol À1 are highlighted in Table S8 of the ESI.† Fig. 9 Comparison of the relative energies from ab initio calculations (solid black lines) with models using transferred b ik and g ik parameters from FA and TA to (a) ClFA and (b) MA respectively.The red curve has only had the baseline adjusted, whereas the blue curve has the a parameter refitted to ClFA or MA respectively.The transferred FA parameters are those N k = 6, s intra = 0.28 kJ mol À1 in Fig. 8, while those of TA N k = 10 s intra = 0.26 kJ mol À1 are highlighted in Table S8 of the ESI.† The grey dotted lines give the conformational profiles for (a) FA and (b) TA from these parameters.