Notes on the computational method being developed for crystal structure and polymorph prediction>
CPOSS is developing a suite of computational methods to tackle the crystal structure prediction of an increasing range of molecules. The appropriate method to generate a worthwhile estimate of the crystal energy landscape depends primarily on the flexibility of molecule and range of functional groups present. Salts, hydrates and other systems with more than one molecule in the asymmetric unit cell require specific methodologies. This crystal energy landscape will show the range of different ways that the molecule can pack with itself, which often provides considerable insight into the experimental crystal structures.
For some purposes, the relative energies of the low energy crystal structures are needed more accurately. This can generally be done by reminimising the structures on the crystal energy landscape with a more accurate (and computationally expensive) method, which again depends on the molecular system in question.
The set of molecules in the database and the corresponding publications give a sense of the range of molecules that we are able to tackle, and what methods were found appropriate. We are very happy to discuss individual molecules and do a preliminary analysis of what combination of our current range of methods could give a worthwhile crystal energy landscape for the particular scientific problem. The following notes show the currently available methods and limitations.
Limitations on molecules
We use ab initio calculations to predict the structure of the molecule and its charge density. This is typically done at the MP2 6-31G** level, though SCF optimisations can be used for larger molecules, and better wavefunctions with larger basis sets and the GDMA2 method give more accurate electrostatic models. Hence, we are limited in the size of molecule to less than about 100 atoms depending on the number of many-electron atoms, such as Cl, S, Br.
Intermolecular forces
The ab initio wavefunction is analysed to give a distributed multipole model of the molecular charge density. This is used to calculate the electrostatic contribution to the lattice energy, essentially to the accuracy of the wavefunction, and far more accurately than can be achieved with atomic charges for most molecules. Thus, our methods are particularly advantageous for polar and hydrogen-bonding molecules, and those with p electron density, where the electrostatic directionality of features such as lone pair and p electrons, play a major role in the crystal packing.
All the other intermolecular forces are usually modelled as a exp-6 repulsion-dispersion model. We have plenty of experience in choosing suitable parameters for molecules involving C/H/N/O atoms from empirically fitted isotropic atom-atom potentials. We can also use anisotropic atom-atom repulsion models, and have developed specific potentials almost non-empirically for molecules containing Cl and Br, allowing high accuracy studies for chlorobenzenes, chlorothalonil etc. Current work exploiting recent developments in the theory of intermolecular forces, developed a non-empirical potential for C6Br2ClFH2, which correctly predicted its crystal structure.
Intramolecular Forces
Our routine search tool MOLPAK and lattice energy minimisation program DMAREL work with rigid molecules. Hence we can readily compute a crystal energy landscape for a molecule that is rigid. Molecules that have a few different low energy conformations, with high and steep energy barriers between them, are readily studied by doing a MOLPAK search for each conformer, and adding the ab initio energy penalty for adopting the alternative conformers to the lattice energy for evaluating the relative stability.
The relative energies of the crystal structures found in such rigid molecules searches can be refined by DMAflex to allow for small changes in conformation, such as small changes hydroxyl torsions or amine pyramidalisation, or other torsion angles, that may improve the hydrogen bonding geometry within the crystal structures.
For molecules where there are a wide range of conformations that are sufficiently low in energy that they might be observed as conformational polymorphs, then searching through this range becomes considerably more challenging. This can be done by doing many MOLPAK searches, each with a different conformation, to find which range of conformations have favourable crystal packings. The most promising structures can then have the low energy torsion angles refined by DMAflex. Hence the degree of flexibility is a major limitation on the molecules that can be studied.
Beyond single molecule predictions.
We have performed searches for monohydrates, diastereomeric salts, and co-crystals, which are all systems where there are two molecular fragments in the asymmetric unit. This can be done by doing a series of MOLPAK searches using likely relative positions of the two molecular fragments, particularly when these fragments will be hydrogen bonded. If there are a huge range of likely relative positions of the molecules within the asymmetric unit, then a more general search method, such as Crystal Predictor, is needed.
Limitations on crystal structures
MOLPAK can generate densely packed crystal structures with one molecule in the asymmetric unit cell (Z’=1), as staring structures for lattice energy minimisation. We routinely use 38 packing types covering 18 of the most common spacegroups (or far fewer for enantiomerically pure molecules), though other are possible. We use this method rather than more extensive search methods because of the relative expense of using accurate models for the electrostatic forces. It may miss some crystal structures within these space groups, if they form an unusual packing type. Since we check whether the minimisations correspond to a genuine minimum rather than a saddle point, and in the latter case, remove symmetry elements until a true minimum is found, we can generate some structures which correspond to Z’=2. Similarly, symmetric molecules can minimise to produce a crystal structure in a higher space-group which includes the molecular symmetry. However, in general the search only considers perfect crystal structures in common space groups with Z’=1.
We can also use Crystal Predictor for more extensive searches, covering a larger range of space groups and explicitly searching for Z’=2 structures, but this is far more demanding of computational resources.
Accuracy of relative energies
Once the set of low energy structures has been established, it is possible to re-evaluate their relative stability using more expensive and accurate methods. This ranges from using alternative model potentials, estimating the effect of the induction energy, using more accurate ab initio models for the intramolecular energy, and refining the molecular conformation by DMAflex. For certain molecules, an electronic structure based method of optimising the positions of all atoms and the unit cell may be worthwhile: we are evaluating the on-going development of such methods for organic crystal structures.
However, the above methods only provide increasingly accurate lattice energies, i.e. relative stabilities of the static crystal at 0 K. We really should be comparing room temperature free energies for most problems. We use the harmonic approximation to estimate the relative free energies at room temperature for molecules that are rigid. Such estimates suggest that the entropy differences generally only make a small difference, particularly if the difference in density between the structures is small. However, there are many cases where molecular flexibility or the onset of significant anharmonic motions means that there will be more significant entropy differences between molecules. Certainly the computational methods do not yet seem reliable for predicting transformation temperatures for polymorphs which change relative stability with temperature. We are investigating the use of meta-dynamics to give relative free energies.
Hence, CPOSS is still working hard to improve the accuracy of the relative energies of the different structures. However, the accuracy required is very dependent on the energy difference between the structures on the crystal energy landscape, the specific molecule and the aims of the study.
Energy Range of Possible Polymorphism
This has not been established experimentally, with a typical estimate being of order 10 kJ mol-1. However kinetics will play a major role, with for example, dehydrated hydrates being likely to be significantly more metastable than the less stable of two polymorphs which crystallise under the same conditions. Although the larger the energy difference, the greater the thermodynamic driving force towards the most stable structure, the kinetic barrier to transformation will be very dependent on the specific structures involved. Hence we tend to look at the energy distribution to decide which the low energy structures of interest are, looking for natural cutoffs and also considering the likely accuracy of the relative energies. We are storing unique structures within a generous energy range, for reconsideration as our understanding of polymorphism and ability to calculate the relative energies improves.
© UCL Chemistry Department 2022. This page was last updated on 17 August, 2022. If you have any problems with this page please email the WebMaster