Accurate and efficient representation of intramolecular energy in ab initio generation of crystal structures. I. Adaptive local approximate models

This article describes an important improvement in the CrystalPredictor II code: adaptive Local Approximate Models (LAMs). This improvement allows the most efficient use of computational effort to cover a flexible molecule’s conformational space, and is illustrated with a crystal structure prediction (CSP) investigation into the sixth blind test molecule 26.


Introduction
The primary aim of crystal structure prediction (CSP) techniques is to produce a ranked list of all the potential crystal structures for a molecule or set of molecules. Because of the significant effect that crystal structure has on solid-state properties, such as colour, solubility and hygroscopicity, such a ranked list offers a wealth of information and many opportunities to improve the development of new crystalline materials (Price et al., 2016;Neumann et al., 2015). In the case of the pharmaceutical industry, the appearance of a new or unexpected form or polymorph can have major legal and economic ramifications, particularly if solubility/bioavailability are affected, as illustrated by the cases of the appearance of Ritonavir form II (Chemburkar et al., 2000) and the Zantac litigation (Seddon, 1999). Furthermore, the ability to tune a molecule's solid-state properties through predictive approaches would be very useful to industries that rely on crystalline materials. Therefore, significant benefits are offered by the possibility of predicting a molecule's crystal structure(s), especially when this is possible via ab initio techniques that rely only on molecular structure information.
Whilst a relatively new field, CSP methods for organic molecules have undergone considerable improvements over the past few years, as seen in the increasing size and complexity of molecular targets in the blind tests organized by the Cambridge Crystallographic Data Centre, as well as the increasing level of success achieved in the tests (Day et al., 2005(Day et al., , 2009; Bardwell et al., 2011;Motherwell et al., 2002;Lommerse et al., 2000). The targets that CSP groups are being asked to investigate as a matter of routine are becoming more industrially relevant, with larger, more flexible molecules that could be seen as drug analogues now being considered. Indeed, in the case of molecule (XXIII) in the sixth blind test (Reilly et al., 2016), the target represents a former drug candidate for the treatment of Alzheimer's disease. All the targets were chosen so as to present challenges that test the theories and computational capabilities currently available.

Global search in CSP
The central tenet of CSP is that the crystal structures that are most likely to form will be low-energy minima on the freeenergy surface, with respect to structural variables, namely: cell lengths and angles, the molecular position and orientation, and the molecule's internal degrees of freedom (Pantelides et al., 2014;Brandenburg & Grimme, 2014;Woodley & Catlow, 2008;Price, 2008;Cruz-Cabeza et al., 2015). Thermodynamically, the most stable crystal structure (at given temperature and pressure) is the global minimum on the Gibbs free-energy surface; however, given the cost inherent in free-energy calculations and the comparatively small energetic contributions arising from entropic effects, most CSP methods use lattice energy/enthalpy rather than free energy in order to rank the predicted crystal structures.
A major factor in successfully identifying all likely polymorphs is the trade-off made between the accuracy of the model used to describe the differences in energy between a molecule's possible structures (often less than 5 kJ mol À1 ), and the extent of the search for low-energy minima across the entire free-energy surface. In view of this, most CSP techniques use a broadly two-stage methodology: a first-stage global search that is used to search for low-energy structures on the lattice energy surface using a relatively low-cost, less accurate lattice energy model; and a second-stage refinement that takes the most promising structures from the first stage and re-ranks them via local energy minimization, using a more accurate and computationally demanding lattice energy model. All the successful predictions in the sixth blind test (Reilly et al., 2016) used some variant of this multi-stage methodology.
In order to identify all potential low-energy polymorphs, the first stage must perform an extensive search (typically involving hundred of thousands of points) of the lattice energy surface over sufficiently wide ranges of the lattice energy model variables (cell lengths and angles, conformational degrees of freedom etc.); therefore, the efficiency of the lattice energy model is very important. Moreover, since only a relatively small proportion (typically a few hundreds) of the lowest-energy structures identified will be passed for refinement to the second stage, the lattice energy model employed by the first stage also needs to be sufficiently accurate not to exclude any potential polymorphs from further consideration.
Overall, achieving the right trade-off between the efficiency and accuracy of the first-stage lattice energy model is a key challenge for CSP. If the accuracy of the lattice energy model can be increased at moderate computational cost, the risk of missing low-energy structures can be decreased. Furthermore, the cost of the second stage can be reduced significantly as increased confidence in the ranked list of structures generated in the first stage typically allows a decrease in the number of structures that must be taken through to the computationally intensive refinement stage, opening the possibility for the latter to employ even higher-accuracy lattice energy models.
Before describing specific advances, we give a brief overview of the algorithm to the extent necessary for the purposes of this paper.

The CrystalPredictor algorithm
The CrystalPredictor algorithm (Karamertzanis & Pantelides, 2007, 2005Habgood et al., 2015) is a global search algorithm based on a large number of gradient-based local minimizations starting from crystal structures generated by a Sobol sequence (Sobol, 1967), a low-discrepancy technique that ensures the best coverage of the space of the variables that uniquely define a crystal structure.
The original version of the algorithm CrystalPredictor I (Karamertzanis & Pantelides, 2007, 2005 was used successfully in several CSP studies to produce initial ranked lists of crystal structures. However, in order to ensure that all experimentally known structures are identified by the CSP, it was often found to be necessary to refine the 1000 to 1500 lowest-energy structures in these initial lists, which resulted in very significant computational costs (see, for example, Vasileiadis et al., 2012Vasileiadis et al., , 2015.
The above issue with CrystalPredictor I was partly caused by the insufficiently accurate description of the effects of molecular conformation on both the intramolecular and intermolecular contributions to the lattice energy. This realisation led to an improved version of the algorithm, CrystalPredictor II (Habgood et al., 2015), using a more accurate energy function that utilizes local approximate models (LAMs; Kazantsev et al., 2010). LAMs allow the efficient and accurate calculation of intramolecular energy as a function of flexible torsion angles ('independent' conformational degrees of freedom, ). Moreover, LAMs also allow the values of those degrees of freedom that are not explicitly treated as flexible in the minimization (the 'dependent' degrees of freedom, " , including bond lengths, bond angles and any torsion angles that are not included in ) to be determined as functions of the independent conformational degrees of freedom, .
CrystalPredictor assumes that the lattice energy of a crystal is given as a function of the cell lengths and angles, collectively denoted as X, as well as the positions and orientations of the molecules in the asymmetric unit, collectively denoted as , and the molecules' independent conformational degrees of freedom, . The optimization then seeks to minimize a lattice energy function, U latt of the form where the intermolecular energy is separated into (a) an electrostatic term, U e ðX; ; Þ, evaluated by the Coulombic attraction between atom centres, based on point charges obtained using isolated molecule ab initio calculations, and (b) a repulsion/dispersion term, U rd ðX; ; Þ, described by Buckingham potentials whose parameters have been fitted to experimental data, typically the FIT potential (Cox et al., 1981;Williams, 1984;Coombes et al., 1996). We note that both U e and U rd are functions of molecular conformation as it affects intermolecular distances. In general, the electronic charge distribution within the molecule is also a function of molecular conformation, and therefore the atomic charges used in U e may also depend on .
The intramolecular energy contribution, U intra is given by where U intra " ; À Á is the intramolecular energy of an isolated molecule at conformation " ; À Á and U gas is the minimum energy of the unconstrained isolated molecule (i.e. in vacuo, with all internal degrees of freedom allowed to vary). To avoid expensive repeated ab initio calculations for the evaluation of the terms U intra and U e during the global search, a set of reference calculations at values of the on a regular grid are performed before the start of the global search, and are subsequently used in CrystalPredictor to obtain a low-cost approximation of these energies at any point. The two versions of CrystalPredictor differ in how the approximation is constructed. In CrystalPredictor II, the intramolecular energy at some value of the independent conformational degrees of freedom is calculated from the LAM with the closest matching conformation, ref , on the grid using an approximation of the form whilst the set of dependent degrees of freedom " is obtained by a linear approximation of the form where the matrices A and C and the vector b are given by (Kazantsev et al., 2011) where n is the total number of independent conformational degrees of freedom and min i and max i , i ¼ 1; . . . ; n, are selected to include all areas of practical interest, typically where ÁU intra is below 20 to 30 kJ mol À1 . LAMs are calculated at grid points whose location depends on the size of the search domain and a userspecified grid spacing . The conformational space is therefore partitioned into hyper-rectangles of the form The LAM validity range, , needs to be small enough to ensure that expressions (3) and (4) provide sufficiently good approximations for the intramolecular energy and dependent degrees of freedom within a certain conformational distance from ref .
The adoption of a regular LAM grid has been found to be effective in CSP for several molecules, such as -d-glucose, ROY and a pharmaceutical compound, BMS-488043 (Habgood et al., 2015). However, the number of LAM points needed to achieve a desired coverage grows exponentially with the number of degrees of freedom. For highly flexible molecules, where the number of independent conformational degrees of freedom is large, and the range of flexibility, [ min i¼1;...n , max i¼1;...n ], to be searched is wide, the number of LAMs to be calculated incurs a high computational cost. In such cases, the choice of an appropriate has a significant impact on both the accuracy and computational cost, and its determination requires substantial analysis of the molecule of interest prior to computing the grid points.

Aims
In this paper we propose improvements to the algorithm that address the issues identified above, leading to reduced computational cost and improved accuracy. In particular, we seek to achieve this by introducing an adaptive LAM implementation so that the LAM points no longer need to be placed on a regular grid.

research papers
A motivating example for the development of an improved algorithm is introduced in x2, based on molecule (XXVI) from the sixth blind test (Reilly et al., 2016). The adaptive LAM placement algorithm is described in x3, and the reduction in computational cost that it offers is analysed. Finally, in x4, we revisit the motivating example, applying to it the improved CrystalPredictor II algorithm in the context of a complete CSP study of molecule (XXVI) from the sixth blind test.

Motivating example: molecule (XXVI) of the sixth blind test
The recent blind test on crystal structure prediction methods, organized by the Cambridge Crystallographic Data Centre, sought to evaluate the capabilities of current computational methods in predicting the crystal structures of organic molecules. Five targets were chosen, representing challenges to the crystal structure prediction community. The two versions of CrystalPredictor were deployed by two of the participating groups, in combination with CrystalOptimizer, to identify Z 0 = 1 structures. This approach resulted in the identification of the known experimental structures within the predicted energy landscapes in most cases. However, in the case of molecule (XXVI), shown in Fig. 1, the multiple flexible torsion angles present particular difficulties, which are discussed here and motivate the development of an improved version of Crystal-Predictor II.
Molecule (XXVI) contains the common 1,1 0 -binaphthalene fragment, which can feature axial chirality, although no chiral precursors were present in the synthesis. As reported by Reilly et al. (2016), there are currently two known pure experimental forms: form (1) is a Z 0 = 1 structure crystallized in the P " 1 1 space group, while form (11) is a structure discovered through polymorph screening by Johnson Matthey (Pharmorphix), after the conclusion of the blind test, for which no structural data are currently available. In addition, there are nine reported solvates. Unusually for 1,1 0 -binaphthalene-based molecules, one O atom is unsatisfied in terms of hydrogen bonds, and the angles and torsions in the amide group and phenyl rings are somewhat outside expected ranges. This is a result of the bulkiness of the 1,1 0 -binaphthalene and phenyl groups, as well as the internal hydrogen bond occurring between the chlorine in one half of the molecule, and the amide group in the opposite half. The number and unusual values of the independent degrees of freedom, as well as its sheer size, contribute to the difficulties posed by this molecule.
In the sixth blind test (Reilly et al., 2016), the use of Crys-talPredictor I and CrystalOptimizer by the Price et al. group successfully led to the identification of form (1) as the lowest energy structure in the final landscape. The use of Crystal-Predictor I, however, required making severe assumptions on flexibility to limit the computational cost; as is usually done with CrystalPredictor I when there are many flexible degrees of freedom, the flexible torsion angles were divided into three groups (group 1 containing T1, group 2 containing T3-5, and group 3 containing T7). This approach has been successful in other investigations (Vasileiadis et al., 2012) and relies on the assumption that flexible torsions in distinct parts of the molecule can rotate independently, with their effect on U intra being largely unaffected by the values of the flexible torsions in the other torsion groups. However, in the case of molecule (XXVI), since the benzene groups that rotate in the different halves of the molecule are in close proximity to each other, such an assumption may not be fully justified in this case. The loss of accuracy arising from this treatment was acknowledged by the Price et al. team and was countered by applying the second-stage refinement to a wider than usual range of the low-energy crystal structures predicted in the CrystalPredictor I landscape. More specifically, a single iteration of Crystal-Optimizer was applied to each of the 9400 structures identified by CrystalPredictor I within 40 kJ mol À1 of the global minimum, thereby resulting in re-ranking of the structures. The full CrystalOptimizer calculation was then performed for the 1322 lowest-energy structures. Although in this case the experimental form was successfully identified, this decom-   position approach is not generally applicable to all molecules, and a more accurate method of covering the conformational space is needed (see Habgood et al., 2015) for a more complete discussion). Our research group's submission for molecule (XXVI) made use of CrystalPredictor II, with all seven torsional angles shown in Fig. 1 being treated as independent degrees of freedom, . The domain of each angle that was deemed to be relevant for CSP purposes was initially decided by analysis of crystal structures in the Cambridge Structure Database (CSD) and the results of one-dimensional scans through each independent degree of freedom. A LAM validity range Á of AE15 was used for all torsional angles, resulting in a grid comprising 2592 LAM points (see Table 1). The computation of the latter at the HF/6-31G(d,p) level of theory required approximately 200 000 CPU h [typically on Intel(R) Xeon(R) CPU E5-2660 v2 running at 2.20 GHz].
Our normal practice in the applications of CrystalPredictor II is to also evaluate the intramolecular energies at the edges of the search space; if these are found to be lower than a userspecified threshold (typically 20-30 kJ mol À1 ), the search space is expanded. In the case of molecule (XXVI), this investigation identified energies lower than 10 kJ mol À1 on the boundaries of the domains for torsions T3 and T5, and therefore these domains would normally have to be expanded quite significantly. However, a larger regular grid with the domain of the two key torsions extended by the necessary 120 would involve 11 858 LAMs, and their construction would require approximately 910 000 CPU h. As this was impracticable within the time constraints of the blind test, it was decided not to extend the search beyond the domains indicated in Table 1.
As indicated in Table 1, the experimental value for the torsion angle T5 is 222.46 , which unfortunately lies outside the search domain [95 , 185 ]. As a result, our search failed to identify form (1) of molecule (XXVI). This illustrates the importance of developing new techniques that would allow large conformational spaces to be covered efficiently by LAMs, which in turn provides the motivation for the present work.

An algorithm for adaptive LAM placement
This section presents an adaptive algorithm that automatically positions LAMs at points in the search domain of the independent degrees of freedom, where necessary to ensure the required degree of accuracy. Firstly, the revised algorithm for generating new LAMs is summarized in x3.1, with examples of its implementation given in xx3.2 and 3.3.

Adaptive generation of LAMs
The basic idea of the adaptive LAM placement algorithm proposed in this paper is to take an existing set of LAMs placed over the search domain of the independent conformational degrees of freedom , and to try to identify a point at which these LAMs may not attain the required accuracy. If such a point is found, then a new LAM is then generated at that point. The procedure is repeated until no new point is found to be necessary.
Establishing the exact error of the approximation provided by a LAM at a particular point would require performing the corresponding quantum mechanical calculation and comparing its results with the LAM predictions. As this would defeat the purpose of an efficient LAM placement algorithm, we choose to use an approximate criterion based on the difference in the predictions at a point between two neighbouring LAMs.
In particular, we assume that the maximum discrepancy between the predictions of two LAMs generated at points A . and B . respectively, is likely to occur around the mid-point M ¼ ð A þ B Þ=2. Using equation (3), we can then easily compute the quantities ÁU intra ð M Þ using the two LAMs. If we denote these by ÁU intra A ð M Þ and ÁU intra B ð M Þ, then a new LAM is generated at point M only if these quantities differ by more than a certain specified threshold, Á Ã , in absolute value, i.e.

ÁU intra
However, before deciding whether to generate a new LAM at M, there are two additional conditions we need to consider. First, it is unnecessary to generate a LAM at point M if the latter is unlikely to be inside the region which would be relevant for the purposes of CSP, i.e. if ÁU intra ð M Þ exceeds a given threshold, Á ÃÃ . Of course, the exact value of ÁU intra ð M Þ is not known, but it can be approximated by the values obtained by the two LAMs. Conservatively, we choose to consider the lower of these two values; therefore, another necessary criterion for a LAM to be generated at point M is A second consideration that needs to be taken into account is that the above reasoning is valid only if the LAMs A and B are indeed those nearest to point M. If there exists a third LAM C which is nearer to M than either A or B, then of course the accuracy of the approximations provided by the LAMs at A and B at point M is irrelevant: neither of those would be used during the search to determine the quantity ÁU intra ð M Þ. Therefore, a third necessary criterion for a LAM to be generated at point M is For each and every existing LAM k other than A and B, where the norm jj:jj is the Euclidean norm in conformational space. The above ideas provide the basis of the new adaptive algorithm for LAM generation. Given any set of LAMs, we consider each and every pair (A, B), determine its midpoint M, and test criteria (9)-(11). If all of those are found to be true, then a new LAM is generated at point M, and the procedure is repeated until no more new LAMs are found to be necessary.
In our current implementation, the algorithmic parameters Á Ã and Á ÃÃ are set by default at 1 and 20 kJ mol À1 , respectively. Using a smaller value of Á Ã leads to increased consistency between LAMs, but also results in the addition of a greater number of LAM points and hence higher computa-research papers tional cost. We have found the value of 1 kJ mol À1 to give an appropriate balance between cost and consistency. The default value of Á ÃÃ is chosen based on the assessment of Thompson & Day (2014) of the maximum energetic cost of molecular distortion away from gas phase conformation in naturally occurring polymorphs. Here, using a larger value of Á ÃÃ increases the reliability of the LAMs for higher-energy conformations, but this again comes at the cost of adding more LAMs. The norm in criterion (11) is based on the Euclidean distance The initial set of LAMs is constructed over a relatively coarse regular grid which is then subsequently refined according to the algorithm presented here, resulting in a complete set of LAMs prior to the start of the global search. A flowchart describing this process is provided in the supporting information. As in the previous implementation of CrystalPredictor II (Habgood et al., 2015), during the search, equations (3) and (4) are applied using the LAM that is nearest, in the Euclidean distance sense, to the current point .

Illustrative example 1: benzoic acid
In order to better understand the concept of adaptive LAM placement, we first consider a molecule with a single independent degree of freedom, namely benzoic acid (see Fig. 2). The chemically relevant domain for torsion angle T1 is initially covered by four LAMs based at the points T1 = À90, À30, +30 and +90 , at the M06/6-31G(d,p) level of theory.
As can be seen in Fig. 3(a), there is clearly a significant mismatch (9.2 kJ mol À1 ) in the intramolecular energy contribution predicted by adjacent LAMs at T1 = AE 60 . This can be corrected by inserting two LAMS at these positions, as illustrated in Fig. 3(b). On the other hand, there is no such mismatch at the boundary between the original second and third LAMs at T1 = 0 , and therefore no new LAM needs to be inserted there. This consistency check, in which different LAM predictions are compared to each other, ensures that the intramolecular energy is described consistently by the LAMs at the given boundary. It does not, however, guarantee that that ab initio accuracy is achieved, although we note that LAMs have been shown to represent ab initio results very well in their locality (Kazantsev et al., 2011). In the case of a symmetric molecule such as benzoic acid, the consistency of the LAMs at T1 = 0 could be attributed to the symmetric placement of LAM points and does not imply agreement with the ab initio energy value. This could be addressed through the manual addition of LAMs by the user to break symmetry where appropriate. Overall, achieving the same level of accuracy with a regular grid would require a grid spacing of = AE 15 , i.e. 7 LAMs overall (starting with one based at T1 = À90 ), as opposed to the 6 LAMs shown in Fig. 3(b). Whilst only a small saving is achievable in this simple case, much more marked efficiencies can be achieved for molecules involving multiple independent degrees of freedom, as illustrated by the next example.

Illustrative example 2: the ROY molecule
The adaptive algorithm is further illustrated for the ROY molecule (5-methyl-2-[(2-nitrophenyl)amino]-3-thiophenecarbonitrile) (Yu, 2010), which is considered here to involve two independent conformational degrees of freedom, T1 and T2, as shown in Fig. 4. These two degrees of freedom have broad ranges of flexibility, with T1 2 ½À20 ; 180 and T2 2 ½100 ; 260 , but within this overall conformational space there are large regions that are characterized by high intramolecular energy which are unlikely to be of relevance to CSP.
Starting with an initial uniform grid generated with = AE20 and comprising 20 LAMs, at the B3LYP/6-31G(d,p) level of theory, the application of the LAM generation algorithm results in the final set of 41 LAMs shown in Fig.  5(b). The minimum spacing between these LAMs is 14 ; a regular grid constructed over the original domain would require about 163 LAMs to achieve the same minimum spacing ( ' AE7 ). However, many of these LAMs would be unnecessary: for example, we note that the adaptive algorithm does not introduce any new LAM points in the region ½À20 ; 40 Â 100 ; 160 ½ . Fig. 5(b) also shows the positions of the six known experimental forms of ROY (Yu, 2010). This demonstrates that the algorithm does indeed focus computational effort on relevant areas of conformational space.
The intramolecular energy predictions by the original and final sets of LAMs are shown in Figs. 6(a) and (b), respectively.
It is clear that the low conformational energy regions are not rectangular, i.e. there is significant interaction between the two torsional angles. It can also be seen that the adaptive LAM placement leads to a smoother intramolecular energy surface in these key regions. The intramolecular energy contribution is also computed ab initio over the same range of degrees of freedom at 5 increments and shown in Fig. 6(c). Visual comparison of the three energy landscapes show that key qualitative features are captured by both LAM-based approximations. A more quantitative comparison is presented in Figs. 7(a) and (b), where the differences between the LAM approximation and the ab initio energies are computed at 5 intervals. The average absolute deviation for the regular coarse grid scheme is 0.75 kJ mol À1 , while for the adaptive scheme it is 0.56 kJ mol À1 . More importantly, it is evident that with the regular grid, there are many areas in which the error is more than 5 kJ mol À1 , particularly at the edges of LAM    Intramolecular energy (in kJ mol À1 ) as predicted by LAMs in 0.5 scan across conformational space; (a) under a regular coarse grid (Á = AE20 ), (b) using the adaptive LAM placement of Fig. 5(b). Crosses represent regular LAMs, circles non-uniform/adaptive LAMs and (c) Ab initio intramolecular energy based on a 5 scan.
validity. This can lead to the generation of a low-accuracy energy landscape during the global search, in which some structures are found to have unrealistically low or high lattice energy. Finally, it can be seen that in the areas surrounding the experimental structures (black triangles), improved accuracy is achieved.

CSP investigation of molecule (XXVI)
The proposed algorithm is now applied to molecule (XXVI) from the sixth blind test (cf. x2). As shown in Fig. 1, the molecule has 7 flexible degrees of freedom, several of which have broad ranges of flexibility. In this new CSP study, the search domains for all 7 torsion angles are extended until the intramolecular energies ÁU intra at the edges of the search space exceed 15 kJ mol À1 . While a larger cutoff value has been used in previous work (Habgood et al., 2015), 15 kJ mol À1 is a practical value given the computational cost, and the low likelihood of a molecular distortion with an energetic cost greater than 15 kJ mol À1 (Thompson & Day, 2014). As can be seen by a comparison of the search domains listed in Tables 1  and 2 this now results in much wider domains for T3 and T5 than those used in our original CSP study (Reilly et al., 2016).

Generation of an appropriate LAM set
In applying the algorithm proposed in this paper, we start with a relatively coarse regular grid with increments of 60 for T1, T3, T5 and T7, and 30 for T2, T4 and T6 (see Table 2), again, at the HF/6-31G(d,p) level of theory. Overall, this initial grid requires the generation of 1152 LAMs. We then apply the adaptive LAM placement algorithm of x3 in a single-pass mode, i.e. simply considering all pairs of points in the initial grid and deciding whether to place a LAM at their mid-point. Overall, this results in the generation of an additional 2491 LAMs.
The accuracy gain achieved by the judicious placement of new LAM points is illustrated in Fig. 8

Figure 7
Absolute difference between ab initio and LAM predicted intramolecular energies (kJ mol À1 ) based on a and T7 (indicated by a filled triangle). The figures show the differences between the value of ÁU intra predicted using the nearest LAM and the corresponding ab initio value. The underlying data are generated by varying T1 and T7 in 2 increments, while keeping the other 5 torsional angles constant at the values T2 = 180.0 , T3 = 170.0 , T4 = 70.0 , T5 = 230.0 and T6 = 180.0 . Fig. 8(a) shows results obtained using the initial LAM set on a regular grid. The four nearest LAMs used for this purpose are outside the domain shown. As can be seen, the values of ÁU intra involve non-negligible errors, with a maximum of 5.15 kJ mol À1 across the sub-region and a value of 1.01 kJ mol À1 at the experimental values of T1 and T7. On the other hand, Fig. 8(b) shows results obtained with the final LAM set which now includes a new LAM placed at the position indicated by the open circle. It can clearly be seen that the addition of this single new point in this sub-region results in very substantial reduction in the error in ÁU intra . The maximum error across this sub-region is now 0.27 kJ mol À1 , with the error at the experimental values of T1 and T7 being just 0.09 kJ mol À1 .
As has already been noted in x2, a regular grid that would cover the entire domain of interest at the required accuracy would have to incorporate 11 858 LAMs, whose construction would require approximately 910 000 CPU h on the computing hardware used for this study. Instead, the LAM set determined by the new adaptive LAM placement algorithm requires only 3643 LAMs, an overall reduction of just under 70%.

Global search using CrystalPredictor II
A global search over 1 000 000 candidate structures is performed using CrystalPredictor II, making use of the LAM set determined above. As shown in Fig. 9, this results in 81 unique structures being identified within 10 kJ mol À1 of the global minimum, with 465 and 1413 unique structures being identified within, respectively, 20 and 30 kJ mol À1 . The experimental form is identified as the 130th lowest energy structure, with a lattice energy 12.27 kJ mol À1 greater than the global minimum, and a good reproduction of the experimental geometry (RMSD 20 = 0.595 Å ).

Refinement of low-energy crystal structures using CrystalOptimizer
CrystalOptimizer minimizations are performed on the 1413 unique structures that were identified within 30 kJ mol À1 from the global minimum (cf. Fig. 9). The approach followed was identical to that in our original investigation carried out in the context of the sixth blind test (see supporting information in the blind test paper by Reilly et al. (2016). In particular, intramolecular energy and conformational multipoles were determined using quantum mechanical calculations at the PBE1PBE 6-31G(d,p) level of theory, and an extended set of independent conformational degrees of freedom was considered as seen in Fig. 10. The use of a different level of theory from CrystalPredictor implies that it is not possible to re-use the LAMs generated at the global search stage. If the same level of theory were used, this would result in a reduction of the number of quantum mechanical calculations at the refinement stage.
The resulting energy landscape is presented in Fig. 11. The experimental form is found at the global minimum, with another 17 structures having lattice energy within 10 kJ mol À1 from the global minimum, and 92 within 20 kJ mol À1 . We note that these numbers are significantly lower than the corresponding numbers of structures determined at the end of the global search (81 and 465, respectively); thus, refinement using a more accurate lattice energy model and taking account of a higher degree of conformational flexibility has resulted in substantial clarification of the polymorphic landscape. We also note that the geometry of the experimental structure is reproduced with good accuracy (RMSD 20 = 0.330 Å ), as illustrated in Fig. 12 and Table 3.  Table 3 Structural information for the predicted crystal structure for molecule (XXVI).

Figure 9
CrystalPredictor II energy landscape for molecule (XXVI) based on 1000 000 minimizations and adaptive LAM placement. The square denotes the experimental form, the solid line is the 10 kJ mol À1 cut-off from the global minimum, and the heavy and light dashed lines are the 20 and 30 kJ mol À1 cut-offs from the global minimum, respectively.
The computational cost of the CSP study is summarized in Table 4. The generation of LAM points remains the most significant cost but is now tractable given the high-dimensionality of this molecule.

Concluding remarks
The 2016 blind test (Reilly et al., 2016) revealed that achieving an appropriate balance between computational cost and accuracy in the global search for crystal structures remains a challenge for large molecules. The algorithm presented in this paper addresses this issue by introducing the adaptive placement of LAMs within the CrystalPredictor II algorithm, an improvement on the uniform grid scheme which had proved too computationally demanding to apply to molecule (XXVI). A higher density of LAM points is automatically achieved in chemically interesting areas of conformational space, thereby resulting in a more efficient use of expensive ab initio calculations. This, in turn, allows the CrystalPredictor II algorithm to handle larger molecules and to explore larger areas of conformational space, through an effective global search methodology. The successful application of this new approach to molecule (XXVI) realises one of the aims of the blind tests, namely to drive innovation in CSP by providing unique and challenging molecular systems.
Throughout the CrystalPredictor II calculations, the lattice energy for any given molecular conformation is computed by making use of the LAM that is nearest to . One undesirable effect of this approach is that discontinuities in both the lattice energy and its partial derivatives may occur at points that lie on the boundaries between adjacent LAMs. Such discontinuities may cause numerical difficulties for CrystalPredictor's gradient-based optimization algorithm in cases in which the path of the optimization iterations crosses one or more LAM boundaries. In practical terms, this is usually exhibited by the algorithm reaching a point from which it cannot achieve any Lattice energy landscape following CrystalOptimizer results for molecule (XXVI). The structures generated following the refinement of the 1413 structures generated by the global search stage within 30 kJ mol À1 of the lowest-energy structure are shown. The lowest 100 unique structures span a lattice energy range of 20.8 kJ mol À1 from the global minimum, whilst only 17 unique structures are identified with lattice energies within 10 kJ mol À1 of the global minimum. The square denotes the experimental form, the solid line is the 10 kJ mol À1 cut-off from the global minimum, and the heavy and light dashed lines are the 20 and 30 kJ mol À1 cut-offs from the global minimum, respectively.

Figure 12
Overlay of the global minimum predicted structure (green tubes) generated in the CrystalOptimizer energy landscape, and the experimental structure (grey tubes = C atoms, red = O, blue = N, white = H).

Figure 10
Independent conformational degrees of freedom used in the Crystal-Optimizer investigation of molecule (XXVI). Curly arrows represent torsions, block arrows represent angles. Table 4 Computational cost of CSP for molecule (XXVI).
Step No. of calculations CPU h (approximate) Step 0: construction of LAM regular grid 3643 280 000 Step 1: CrystalPredictor II minimizations 1 000 000 20 000 Step 2: CrystalOptimizer refinements 1413 80 000 Total -380 000 further reduction in lattice energy, despite the fact that the mathematical optimality conditions are not yet strictly satisfied. In the calculations reported in this paper and in our previous work, we have chosen to adopt a conservative approach whereby such points are still considered as candidates for further refinement. However, this may result in much additional computation: for example, in the case of molecule (XXVI), 1283 of the 1413 structures that underwent final refinement (cf. x4.3) actually belonged to this category. Part II of this paper (Sugden & Adjiman, 2016) is concerned with addressing this problem in a more fundamental manner by removing the discontinuities at the LAM boundaries. Data statement: Data underlying this article can be accessed on Zenodo at https://doi.org/10.5281/zenodo.56731, and used under the Creative Commons Attribution licence.