Available methods and codes for materials discovery / crystal structure prediction?

Question

There are many open databases and projects that allow you to access computed crystal structures of experimentally known compounds, registered in the Inorganic Crystal Structure Database (ICSD). These structures are usually then studied using a given materials modeling method, e.g. DFT.

Are there any methods and codes available to predict completely new materials, their crystal structures and properties with no experimental references?

Answer 1

Yes, there are! One early example that's still in use today is the Universal Structure Predictor: Evolutionary Xtallography (USPEX) method. You can find many "success cases" on their website if you're curious. The First-Principles-Assisted Structure Solution (FPASS) and Prototype Electrostatic Ground States (PEGS) methods are a couple of other codes used for structure prediction. Then there are many that are specific to a given application area. For instance, to predict the structures of metal-organic frameworks from specified molecular building blocks, there is the Topologically Based Crystal Constructor (ToBaCCo).

Each code out there works differently. Some try to use brute-force approaches while accounting for material-specific design principles. The lowest energy structure is then often take. You could imagine taking such structures from the Open Quantum Materials Database, for instance. Other methods use evolutionary algorithms, as summarized in this review. More recently, these crystal structure prediction algorithms are taking advantage of machine learning and artificial intelligence for the design and prediction of crystals with tailored properties for a given application of interest, such as a material that can separate $\mathrm{CO}_{2}$ from natural gas.

There are several reviews on the topic you may be interested in: 1, 2, 3.

Answer 2

stk (Supramolecular Toolkit)

Disclaimer: I am the author of stk.

stk, https://github.com/lukasturcani/stk, is an open-source Python library, which, among other things, is able to carry out automated material design through evolutionary algorithms.

stk allows you to customize the evolutionary algorithms, by specifying your own mutation and crossover operations you want to use, which are appropriate for your material, but it also provides some default and built-in ones.

The framework by which stk constructs materials is flexible and extendable, so you can define your own materials, if the built-in ones are not the ones you are interested in. At time of writing organic cages, metal-organic cages, metal complexes, rotaxanes, macrocycles, covalent organic frameworks, linear polymers are all built-in. If you ask him, the author is also quite likely to add any materials you might be interested in, or if you file an issue on GitHub.

stk also let's you deposit any calculated properties and constructed molecules into a MongoDB database for easy future retrieval.

stk documentation, https://stk.readthedocs.io/, contains tutorials on how to run and customize the evolutionary algorithm.

Note that stk does not include computational chemistry packages. Normally you will define a fitness function of the form

def get_fitness_value(molecule):

    # get fitness value somehow
    # for example, by using rdkit or ASE to get physical properties of the molecule
    # or writing the molecule to a file and then executing a computational chemistry package
    # and extracting the result

    return the_calculated_fitness_value

You would then execute the evolutionary algorithm in the following way

import stk
import pymongo

ea = stk.EvolutionaryAlgorithm(
    fitness_calculator=stk.FitnessFunction(get_fitness_value),

    # Other options for the evolutionary algorithm go here, such
    # the mutation, crossover, selection operations you want to use
    # and the initial population
)

# Using a MongoDB is optional, but very useful.
db = stk.ConstructedMoleculeMongoDb(pymongo.MongoClient())

# Run the EA for 50 generations.
for gen_index, generation in enumerate(ea.get_generations(50)):
    # Do stuff with the molecules in the generation.
    for mol_index, molecule_record in enumerate(generation.get_molecule_records()):
        # Like putting them into a MongoDB database.
        db.put(molecule_record.get_molecule())
        # Or writing them to a file
        molecule_record.get_molecule().write(f'{gen_index}_{mol_index}.mol')

Answer 3

Predicting a crystal structure without experimental input is often called "ex nihilo" (out of nothing) structure prediction. Most materials modelling programs have optimisation methods to take a user-supplied initial structure and optimise it (i.e. find a low energy, zero force/stress configuration). These methods are known as "local optimisation methods" because they find the low-energy, zero force/stress structure near the starting configuration, not necessarily the lowest energy structure that is possible. As an example, if your initial structure is a carbon diamond crystal then the software will find the optimal lattice constants of diamond and move the atoms to the high symmetry sites (if they weren't already there).

In order to predict crystal structures, you need to be able to find the structures with the lowest possible energies. Methods to do this are known as "global optimisation methods" and there are many different flavours. A select few are:

Simulated annealing

This method is inspired by the physical processes of crystallisation as matter cools. If your modelling software has molecular dynamics capabilities, then you can do this method yourself. Essentially you heat your initial system up, let the atoms wander around in configurational space, then cool it down and see where they go. Repeat... repeat... repeat... It is not very efficient, so I don't recommend it for "real research" in general, but it's very easy to do.

Basin Hopping

A related, much more sophisticated and reliable method is "basin hopping" which can actually be quite efficient for some classes of problem. In essence, this method tries to identify the "basins of attraction" as it proceeds (a "basin" is the set of configurations which all optimise to the same structure) and to force the system out of the known "basins" and, therefore, into a new one.

See, for example,

"Global structure search for molecules on surfaces: Efficient sampling with curvilinear coordinates", K. Krautgasser et al, J. Chem. Phys. 145(8), 084117 (2016); https://doi.org/10.1063/1.4961259

"Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms", D.J. Wales and J.P.K. Doye, J. Phys. Chem. A 101(28) 5111 (1997); https://doi.org/10.1021/jp970984n

There's even a basin hopping algorithm in scipy, though I've never used it.

Genetic algorithms (GA)

These methods are inspired by the development of populations of animals. They work by generating a large group of candidate structures (the "population"), evaluating how good they are (their "fitness") and, based on that, choosing some of the structures to mix together to create new candidate structures (called "breeding").

The "breeding" usually involves picking a pair of "parent" structures from the population, and creating one or more new "child" structures by (a) taking parts of each of the parent structures (called "crossover"); and (b) adding some random changes ("mutations"), such as perturbing the atomic positions. The new child structures are often filtered to remove some unlikely ones, and the remainder are evaluated for their fitness and added to the population. This forms the next "generation", who can then have children of their own, which can have children of their own etc. You also need to remove poor structures from the population ("death").

In the materials modelling context, it is common practice to optimise each child structure with the usual local optimisation methods in order to make them as good as they can be.

There are a huge number of different genetic algorithms because there are lots of choices, for example:

• How do you evaluate "fitness"?
• How do you use the "fitness" to choose "parent" structures?
• How do you generate the "children" from the "parents"?
• Which mutation changes do you allow?
• Do you let all the children enter the population?
• When do you remove structures from the population?

Genetic algorithms can be very efficient for materials prediction, however you should beware the huge number of choices! Even once you've answered all the questions above, you still have to decide how many structures should be in each generation, how many generations to run, what the actual mutation probabilities should be, which properties you wish to include in your fitness measure...

"Improved real-space genetic algorithm for crystal structure and polymorph prediction", N. L. Abraham and M. I. J. Probert, Phys. Rev. B 77, 134117 (2008); https://doi.org/10.1103/PhysRevB.77.134117

"Crystal structure prediction using ab initio evolutionary techniques: Principles and applications", A. R. Oganov and C. W. Glass, J. Chem. Phys. 124 244704 (2006); https://doi.org/10.1063/1.2210932

"XtalOpt: An open-source evolutionary algorithm for crystal structure prediction", D. C. Lonie and E. Zurek, Comp. Phys. Comm. 182 (2) 372-387 (2011); https://doi.org/10.1016/j.cpc.2010.07.048

And, if you'll forgive me advertising my own student's work:

"Simultaneous Prediction of the Magnetic and Crystal Structure of Materials Using a Genetic Algorithm", E. J. Higgins, P. J. Hasnip and M.I.J. Probert, Crystals 439 (2019); https://doi.org/10.3390/cryst9090439

There's an implementation of a genetic algorithm in the popular Atomistic Simulation Environment (ASE; https://wiki.fysik.dtu.dk/ase/), as well as XTalOpt and USPEX.

Bio-inspired algorithms

There are a huge number of methods drawn from the behaviour of living creatures, many based on how creatures find resources. Particle swarm optimisation, ant colony optimisation, bird-flocking... For example:

"Crystal Structure Prediction via Particle Swarm Optimization", Y. Wang, J. Lv, L. Zhu, and Y. Ma, Phys. Rev. B 82 094116 (2010); https://doi.org/10.1103/PhysRevB.82.094116

(The above work is implemented in the CALYPSO package.)

Ab Initio Random Structure Searching (AIRSS)

So far, all of the methods described have been general global optimisation methods which can be applied to pretty much any optimisation problem. In the context of materials modelling, this means they can be used with continuum models, coarse-grained molecular models, atomistic forcefields and quantum mechanical methods.

In recent years a much simpler scheme has demonstrated great success across a range of materials problems. The Ab Initio Random Structure Searching (AIRSS) method is, at its heart, a Monte Carlo-style random sampling of configurational space, coupled with the usual local optimisation of each new structure. It exploits the observation that, when the electrons are treated quantum mechanically, the low energy material structures lie within broad basins, which are relatively easy to find. This is not generally true of optimisation problems, and in fact even within materials modelling it is not true of many forcefields, which often have jagged energy landscapes; it does appear to be true of most "real" materials though, and this is reflected in the results of quantum mechanical models.

"Ab initio random structure searching", C.J. Pickard and R.J. Needs, J. Phys.: Condens. Matter 23 053201 (2011) ; https://doi.org/10.1088/0953-8984/23/5/053201

The AIRSS method is implemented in the AIRSS software: https://www.mtg.msm.cam.ac.uk/Codes/AIRSS

Answer 4

Structure prediction involves two steps:

1. Propose candidate structures. This step is necessary to explore configuration space, and different codes use different strategies. Examples of methods that have been extensively used in the solid state include stochastic, evolutionary, or particle swarm algorithms to generate candidate structures. Features that can help when generating candidate structures include the use of symmetry, knowledge about the chemistry of the material (e.g. stoichiometry), or partial knowledge about the structure (e.g. lattice constants but no atomic positions). These features effectively reduce the high dimensionality of the configuration space that is being explored, thus greatly reducing the computational cost. However, even when absolutely nothing is known about a material, these methods can provide good candidate structures.
2. Rank the candidate structures according to their energy. The thermodynamically stable structure of the system is the one with the lowest energy (or enthalpy at finite pressure, or Helmholtz free energy at finite temperature, or Gibbs free energy at finite pressure and temperature). In structure prediction, this problem is reduced to finding the lowest energy structure from the subset of candidate structures identified in the first step. Calculations are typically performed at the DFT level for solid state structure predictions, but of course can in principle be performed at any level of theory, with the only limitation being computational resources avaible (there may be many candidate structures).

How reliable are these methods? It is impossible to prove purely theoretically that a given structure is the ground state of a material. You could always have missed the real structure in your pool of candidate structures (e.g. due to computational limitations you may have searched for structures up to 100 atoms in the primitive cell, but the correct structure may have 200 atoms). This is why this is such a challenging problem. However, in practice structure prediction methods have been quite successful at predicting correct structures.

For actual codes, the Wikipedia page has a reasonable list, and I would encourage everyone to add to the list.