Simultaneous Prediction of the Magnetic and Crystal Structure of Materials Using a Genetic Algorithm

We introduce a number of extensions and enhancements to a genetic algorithm for crystal structure prediction, to make it suitable to study magnetic systems. The coupling between magnetic properties and crystal structure means that it is essential to take a holistic approach, and we present for the first time, a genetic algorithm that performs a simultaneous global optimisation of both magnetic structure and crystal structure. We first illustrate the power of this approach on a novel test system—the magnetic Lennard–Jones potential—which we define. Then we study the complex interface structures found at the junction of a Heusler alloy and a semiconductor substrate as found in a proposed spintronic device and show the impact of the magnetic interface structure on the device performance.


Introduction
In order to meet the challenges posed by modern and emerging technologies, it is increasingly necessary to look beyond existing, known materials. Many fields, from solar cells to spintronic devices, call for materials with unprecedented performance characteristics, or even entirely new behavior. Searching for new materials experimentally is expensive and time-consuming, but the advent of efficient, accurate computational materials modeling offers a potential way forward. Magnetic materials are of particular interest, with applications from fast, high-density data storage such as magnetic RAM devices [1] and heat-assisted magnetic recording (HAMR) [2], to new spintronic and quantum devices, such as spin-valves [3,4]. Magnetic materials include conventional ferromagnets along with more exotic structures, such as antiferromagnets, ferrimagnets, and spin glasses. These materials are already at the heart of many important technologies, but play an increasingly important role in developing and future technologies.
Many of the strong permanent magnets in use today rely on rare-earth elements [5] and there is concern over the sustainability of these elements' availability [6]. Therefore, developing ferromagnets made from more readily available materials is desirable.
In order to determine a material's properties, it is not sufficient to know its chemical composition, it is also vital to understand the crystal structure of the material. In general, atoms will arrange themselves in a material so as to minimise the total (free) energy. One method to determine crystal structure is to estimate an initial atomic configuration, compute the atomic forces and lattice cell stress according to a suitable model, and then adjust the atomic positions and lattice parameters in order to minimise the energy (where the forces and stresses are zero). This procedure works well if the initial geometry is sufficiently close to the true ground state, and this workflow has been the backbone of computational materials studies for many decades. However, the resultant optimised structure only represents the local energy minimum; that is, this procedure finds the lowest energy structure which is closest to the initial geometry. It is possible that a lower-energy structure exists, but that it is separated from the initial geometry by an energy barrier which the local optimisation procedure will not overcome. This situation is depicted in Figure 1. In fact, many materials do possess a range of low-energy phases, and determining the stable phase(s) under given environmental conditions is extremely challenging, even for well-studied classes of materials. In general, finding the full range of stable phases requires a global optimisation method, usually involving a wide range of simulations, encompassing atomic configurations which span all the possible phases. Developing global optimisation methods is an active research field in materials modelling, and existing methods include basin hopping [7][8][9], random structure searching [10], meta-dynamics [11], swarm optimisation [12] and evolutionary algorithms. The latter class includes genetic algorithms [13][14][15][16], which have met with great success in materials modelling and are the focus of this work.

E(x)
Conventional global optimisation methods for crystal structure prediction often work in the basis of the atomic positions and lattice vectors. However when considering systems with nontrivial magnetic structures, this representation leads to a multi-valued energy landscape to search, even when the crystal structure itself is relatively well-established. As a result of this, conventional algorithms can struggle when it comes to searching the combined phase space of magnetic-and crystal-structures for magnetic materials.
Therefore, in this work, we present a novel genetic algorithm (GA) for the simultaneous global optimisation of both the magnetic structure and the crystal structure of a material. We shall explain the key ideas, and then illustrate with two example studies. The first is a new model system, suitable for the study of magnetic systems at the atomic scale-the magnetic Lennard-Jones potential. Whilst the conventional Lennard-Jones potential is a well-known test system for many structural and dynamical algorithms, it has not yet (to our knowledge) been extended to study magnetic systems. We shall introduce this system and explore some of its fundamental behaviour as a means to test our new magnetic GA. The second system we shall study is more complex, and is inspired by recent experiments on spintronic systems, and is an interface between a Heusler alloy and a Ge substrate. We shall use the magnetic GA to find the optimal interface structure and show how the novel predicted structure can explain the experimental data.

Materials Modelling
The first step in finding a stable crystal structure is to have a reasonable model for the energy of any particular atomic configuration. There are two main classes of models in this context: classical forcefields; and quantum mechanical methods.
Classical forcefields are multivariate functions which take the atomic species and positions as input, and return the internal energy of the system and the atomic forces (the first derivative of the energy with respect to the atomic positions). Since atomic interactions depend on the positions of the atoms relative to each other, rather than the absolute positions themselves, forcefields are usually expressed in terms of bond lengths and, depending on the particular forcefield, bond angles and bond torsions. One of the earliest and simplest interatomic forcefields was proposed by Lennard-Jones [17], and expresses the energy E LJ of an atomic configuration as where r ij is the distance between atoms i and j and ǫ ij and σ ij are parameters of the potential which control the equilibrium bond energy and bond-length, respectively. Forcefields are generally computationally cheap and simulations of millions of atoms may be carried out routinely. The main drawbacks of forcefield methods are the need to select and parameterise an appropriate forcefield for the material, and the inability of most forcefields to model dynamic changes in chemistry such as bond breaking or formation.
Quantum mechanical approaches centre on solving the many-body Schrödinger equation for the electrons and nuclei which comprise the material. These approaches have two principal advantages over forcefield methods: they are parameter-free since the form of the Schrödinger equation is known exactly; and by modeling the electrons explicitly, the approaches handle seamlessly chemical complexities such as bond formation and breaking, or environmentally-dependent oxidation states. There is one significant drawback, however, in that this is an extremely computationally intensive task and solving the many-body Schrödinger equation in its original form is unfeasible for materials. Most practical quantum-mechanical materials simulations are based on a reformulation of quantum mechanics known as density functional theory (DFT).
DFT is founded on the Hohenberg-Kohn theorem [18], which proves that the ground state energy of a many-body quantum system is a unique, universal functional of the particle density and has no explicit dependence on the many-body wavefunction. Kohn and Sham [19] used this to demonstrate a formal link between the ground state of the many-body Schrödinger equation, and the ground state of a related auxiliary set of coupled single-particle equations. Crucially, these auxiliary equations may be solved with far less computational resources than the original many-body expressions [20,21].
In principle the mapping between the many-body system and the auxiliary system is exact, but there is one term in the auxiliary equations which is unknown: the exchange-correlation potential. The exchange-correlation potential must be approximated in practical calculations, with the most common approximations being built on known limits and numerically-exact results for weakly interacting electron gases, e.g., the approximation of Perdew, Burke and Ernzerhof [22]. Despite these approximations, DFT has enjoyed popularity and great success in a wide range of materials simulations [23,24].

Genetic Algorithms
A common solution to the need for global energy optimisation is to use an ensemble, or 'population', of starting atomic geometries. Each of these candidate geometries can then be optimised independently with a local optimisation method, in an attempt to map out the possible local minima. Since the initial geometries are unlikely to span enough of the configuration space to find all the local minima, the geometries are typically updated in order to explore other regions of the configuration space. Due to the independent nature of these members, population-based methods often have more scope for parallelism than other methods.
One popular class of population-based methods for structure prediction is that of genetic algorithms (GAs) [13][14][15][16]. In these algorithms, the 'fitness' of members of a population is evaluated in some appropriate manner (e.g., from the binding energy per atom) in an attempt to decide which of them are most likely to find solutions of interest. The population is then updated by selecting favorable (e.g., low-energy) members of the population, and generating new population members by combining these members to produce 'child' members, a process known as 'crossover', followed by some random mutation(s). Depending on the problem at hand, a variety of methods exist for this fitness evaluation, member recombination, and selection stages.
When considering how to represent an optimisation problem for a genetic algorithm, a choice needs to be made about how the search vector is to be represented for these operations. A simple GA may treat the input vector simply as a 1D array of real numbers or as a bit string representing that vector. However, it is often the case that a more physically meaningful representation of the search vector can improve the convergence characteristics of the algorithm. In addition to pure GAs, a number of hybrid GAs exist, incorporating ideas from other algorithms. For example, memetic algorithms allow each member of the population to perform local searches.
In this work, we started from the GA developed by Abraham and Probert [14], which was created specifically for the prediction of crystal structures in periodic boundary conditions. In this GA, the crossover is performed in real-space by using a pair of periodic cuts to select material from each parent structure. This material is combined to form child structures, as illustrated in Figure 2. Real-space crossover of two unit cells, each containing six atoms across two different atomic elements (black and white circles) in two dimensions using a pair of periodic cuts (dashed lines). The cuts define two sets of atoms for each of the parents (left): the 'inner' set comprises atoms lying between the two cuts; and the 'outer' set comprises the remainder of the atoms. Child structures (right) are formed by combining atoms from the 'inner' set of one parent, with those from the 'outer' set of the other parent. Following crossover, the child structures are mutated before a local optimisation method is performed. Several different types of mutations are allowed: deformation of the cell vectors; perturbations of the atomic positions; and inter-atomic swaps (see Figure 3). Mutations are important to GAs, as they are the only way of introducing fundamentally new structural features which are not present in previous populations. Once the local optimisation method has optimised all of the population structures, the fitness of the structures is evaluated. Based on the fitness, some structures are removed from the population (preferentially those with low fitness) and some structures are chosen for crossover (preferentially those with high fitness).
One potential weakness of the method outlined thus far is the tendency of population members to become more and more similar, a process known as 'stagnation'. In extreme cases, every single population member may converge to the same structure; this structure usually represents a stable, low-energy phase but it is not necessarily the global minimum-energy structure, and once a population has stagnated it is almost impossible for it to diversify again significantly.
In order to prevent stagnation, the GA was extended in Ref. [25,26] to incorporate a structure-factor based fingerprinting technique, which enables the GA to differentiate structures in order to penalise structures that are too similar, lowering their fitness and encouraging diversity within the population. For a given structure x containing N atoms, the (non-magnetic) fingerprint is defined as: where r i is the position of atom i, ρ i is the charge density of the nucleus of atom i, defined to be Z i at r i and zero elsewhere, and j 0 is the spherical Bessel function of the first kind. For each structure, the value of Λ is calculated for a range of values of wavenumber k. Using this fingerprint, the difference between two structures can then be quantified, for example using an R-factor inspired by the Pendry R-factor. Using this, the difference between two structures x and x ′ is defined as:

Extending the GA for Magnetic Materials
A GA may be extended to optimise magnetic systems by the simple expedient of including magnetic effects in the calculation of energies and forces, and allowing the local optimisation method to minimise the energy with respect to the magnetic degrees of freedom as well as the crystal structure. However, this naïve approach suffers from two severe problems: firstly, the magnetic energy landscape itself has multiple minima, necessitating a global energy optimisation method; secondly, the magnetic structure and crystal structure of a material are often coupled, and may not be treated independently of each other. Elemental iron may be considered a prototypical example of both effects, possessing an antiferromagnetic face-centred cubic phase in addition to the familiar ferromagnetic body-centred cubic phase (see Figure 4). For these reasons it is important to consider the magnetic-and atomic-structure of a material on an equal footing, and to be able to predict both simultaneously.  In order to account for magnetic effects within the GA itself, the magnetism needs to be included within the representation of the structure and evolved using updated GA operations such as crossover and mutations. For the representation of magnetism in the system, this work assigns an atomic spin to each atom, which can be either an additional degree of freedom for collinear spin systems or an additional three degrees of freedom for non-collinear spin systems. When considering DFT simulations, spin information will typically be represented as the electronic spin density across all space. In order to project this onto the atomic representation, Hirshfeld analysis [28] was used to partition the electronic spin density into regions of space associated with each atom. This information can then be transmitted from parent structures to offspring during the crossover operation.

Perturbation/Permutation Operations
In order for a GA to optimise the spins efficiently, mutation operations need to be extended to also affect the spin degrees of freedom. This work defines two spin mutations: perturbation and permutation of the atomic spins.
In the case of perturbations, a distinction needs to be made between collinear and non-collinear spin systems. For collinear spin systems, there is only one scalar value per atom to optimise. The perturbation therefore takes the same role as the atomic position perturbations, where the role of the perturbation is to move the system from one basin of attraction to another. Since there is a maximum value these spins can take, i.e., the total number of electrons associated with the atom, issues may arise trying to add an additional perturbation to an existing spin. For example, an atom with saturated spin given a positive perturbation will become unphysical. As a result, the spin perturbation is evaluated as a uniform random number between -spin_max and spin_max, where spin_max is a user parameter, set to the maximum spin value expected on any atom. For example, for a d-block transition metal, spin_max should be set to 5h/2, since the d shell is capable of having a maximum spin of 5h/2. This perturbation allows atomic spins to spontaneously magnetise, demagnetise or flip, irrespective of the initial value.
In the case of non-collinear spins, a similar procedure is performed. In this case, the atomic spin is set to a random vector within the sphere of radius spin_max. Again, this allows any spin state to be found by the perturbation within the range specified, without the risk of being stuck in a state based on the structure's history.
For permutations, the situation is analogous to swapping the atomic positions of different species. However, since atomic spin is closely related to the crystal structure, it only makes sense to swap atomic spins of atoms of the same species. For example, swapping the spin of a nickel and oxygen atom in NiO makes no physical sense, since all the spin exists on the Ni atoms, and the O atoms generally have zero spin.  We shall perform all our calculations using the general purpose DFT code CASTEP [29], which relaxes the electronic charge density and spin (in either collinear or non-collinear form) to find the electronic ground state. There may be multiple local minima with different spin configurations.

Start
Hence the input atomic spin from the GA is used to initialise the electron density, and CASTEP then acts as a local optimiser for the atomic spin. This process may then be repeated as CASTEP performs a geometry optimisation to find the (local) minimum energy configuration of atomic positions and cell vectors ( Figure 5).

Magnetic Fingerprinting
To encourage diversity in the population, it is important to be able to quantify the differences between structures. In Abraham and Probert [25,26], a fingerprint was introduced that was translationally and rotationally invariant, and could successfully distinguish different structures. However, it will fail to distinguish two magnetic systems that have the same crystal structure but a different arrangement of spins. Hence we propose an augmented fingerprint Λ aug , inspired by magnetic scattering experiments, that can differentiate magnetic structures [30]: where Λ(k) is the original crystal structure fingerprint as defined in Equation (2) and is the magnetic structure fingerprint term. Here S i is the spin on atom i and q is some scaling parameter which allows the overall fingerprint to be more or less sensitive with respect to differences in magnetic structure and crystal structure. The R-factor in Equation (3) can be used with this augmented fingerprint to quantify the difference between two magnetic structures. As an example, the effect of perturbations on a 6-atom Fe unit cell is shown in Figure 6. The black crosses represent the R-factor difference of structure with perturbed atomic positions with respect to the original reference structure. The red squares and blue circles represent similar perturbations with an additional perturbation to the atomic spins. Here a value of q = 5 was used to make significantly different spin structures (∆S ≥h/2) appear as different as significantly different crystal structures (∆r > 0.2 Å), both of which can be represented by an R-factor difference of R xx ′ > 0.03.

Fictional Magnetic Potential: LJ + S
In order to test the GA on magnetic systems, a pair-potential model including magnetic effects was defined. Since a significant portion of these effects are intrinsically quantum mechanical or many-body in nature, it is difficult to capture all of these effects in empirical potentials without limiting them to some specific regime. Magnetic moments on isolated atoms tend to arise from partially filled electronic states, as dictated by Hund's rules. As a result, atoms with nonzero magnetic moments tend to be polyvalent, causing an additional challenge to empirical pair potentials.
Since magnetic effects are often too complex to include explicitly in pair-potentials, their effects are usually included implicitly through the mechanical and thermal properties used to parameterise them [31]. There are some empirical potentials however that try to explicitly include magnetic effects [31][32][33] and some machine learning models have been trained to deal with magnetic systems [34]. These models tend to require far more complex effects than simply pair-wise interactions.
Given all this, attempting to describe real magnetic materials using a pair-potential is beyond the scope of this work. Instead, the spatial dependence of the energy from the standard (non-magnetic) Lennard-Jones potential was combined with the magnetic dependence of the energy from the lattice-based Heisenberg model in order to create a potential which behaves somewhat like a magnetic material. This work defines a generalised magnetic Lennard-Jones potential as: Dipole term Anisotropy term (6) In general, a magnetic material may have contributions to the potential from exchange between sites i and j (with spatial dependence given by the f ex (r ij ) term; a dipole-dipole interaction of the given form; and a symmetry-dependent anisotropy term.  (3)) comparing the magnetic fingerprints (Equation (4)) of a 6-atom Fe unit cell to the same cell with a range of perturbations applied to the atomic positions. In addition to just performing perturbations to the positions (black crosses), the processes was repeated with an additional perturbation to the atomic spins (red squares and blue circles).
In this study, the magnetic GA was tested on a reduced form of this potential where B = C = 0 and f ex = exp(−αr ij ). However this form of f ex could be modified to include more complex forms of the exchange interaction such as the RKKY interaction. The form of the potential used is shown in Figure 7.

CFAS/n-Ge Interface
One material that is of interest for its magnetic properties is the full Heusler alloy Co 2 FeAl 0.5 Si 0.5 (CFAS). CFAS is a half-metal and has been proposed as a candidate material for novel spintronic applications such as spin valves and magnetic tunnel junctions. In order to use Heusler alloys for device applications, the material needs to be attached to electrical contacts. For a number of spintronic applications of Heusler alloys, both metallic and semiconducting contacts are required.
It is known that the half-metallic nature of Heusler alloys at the interface can be sensitive to the exact configuration of the interface atoms, and this can affect the performance of the device. For example, it is known that Co 2 MnSi (CMS) grown on a silver surface changes significantly depending on whether it terminates with a Co layer or a Mn/Si layer [35]. For the CFAS/semiconductor interface, germanium provides an extremely good lattice match with only a 0.2% mismatch, compared to around 4.5% for CFAS/Si. However, when CFAS is grown on n-type germanium, significant mixing at the interface can occur, as evidenced by energy dispersive X-ray spectroscopy (EDS) [36] (see Figure 8). In addition, annealing the sample at 575 K forms a plateau in the cobalt intensities and, to a lesser extent, the iron and germanium intensities, as seen in the EDS results in Figure 9. This suggests that a stable phase is forming at this point. Since this is in the centre of the interface region and all the EDS intensities are about half of their respective bulk values, it is proposed that this structure contains half the number of atoms in the bulk CFAS and half the number of atoms in a germanium unit cell. This results in a potential Ge 4 Co 4 Fe 2 AlSi structure, containing 12 atoms in the formula unit, although the crystal structure of this interface phase is unknown. Since the electronic properties of this phase will affect the half-metallicity of the CFAS/n-Ge interface, the structure of this phase was investigated using the new magnetic GA.

LJ+S
For real-world magnetic materials, it is often the exchange term which dominates the magnetic interaction. Because of this, the work presented in this section will use a simplified version of the LJ+S potential, ignoring the dipole and anisotropy terms (i.e., choosing B = C = 0) and using the simple exchange interaction form f ex = exp(−αr ij ).
In this case, only 2 parameters need to be chosen, A and α. As a result of this, the effect of the magnetic modifications to the Lennard-Jones potential are to raise and/or lower the energy of aligned and anti-aligned spins, depending on the sign of A. If A is positive, pairs of aligned spins will act to lower the energy and anti-aligned spins will raise it. In this case, ferromagnetic magnetic structures are expected to have the lowest energy, since aligned spins would act to decrease the total energy of the system. If A is negative, pairs of anti-aligned spins will lower the energy and pairs of aligned spins will raise it. In this case, antiferromagnetic structures are expected to have the lowest energy, since aligned spins would act to increase the total energy of the system.
In order to get physically sensible values for A and α, they can be parameterised in a number of ways. The way which has been chosen here is to parameterise the values to the equilibrium distance r min of an Fe dimer, along with the difference in energy at r min and 2r min between aligned and anti-aligned dimers, as demonstrated in Figure 7. These values have been found by performing a DFT calculation using the CASTEP code and fitting the results to the V LJ+S potential form (Equation (6)).

Algorithm Performance
The magnetic GA was run on two different six-atom systems of the LJ+S potential. The first used the parameters listed in Table 1 and the other had the same parameters but with the sign of A inverted. This provides both a ferromagnetic and antiferromagnetic system on which to test the algorithm. The GA was allowed to optimise the atomic positions, spins and lattice parameters. The GA was run with population size of 24 members per generation, and ran for a maximum of 50 generations for the FM case and 80 generations for the AFM case. Each member was locally optimised using the TPSD algorithm for 50 iterations of the local optimisation to get into the quadratic region, then further converged with the BFGS algorithm. The structures were converged to a tolerance of 1 meV. Since the LJ+S system spins could not be locally optimised by CASTEP, the GA spin mutation operation was modified to be a normally distributed rotation around the surface of the unit sphere, centring on the spin's previous position.
A mutation amplitude of 2 Å was used for the ions, at a rate of 0.03 mutations per atom. This was chosen to allow on average one atom every member to mutate, but significantly enough to hop lattice sites. The spin mutation rate was 0.06 per atom, since the mutations in the spin would be less dramatic. Figure 10 shows the enthalpies of each structure as a function of the sum of the z component of the spins for both the AFM and FM parameterisations. It can be seen that, as the calculation progresses, both the enthalpies and spins of each system converge to their stable states, i.e., ∑ S z = 0h/2 for the AFM case and ∑ S z = 6h/2 for the FM case. Since the magnitude of the spins on each atom is kept constant at | S i | = 1h/2, these states correspond to fully AFM and fully FM structures respectively.  Figure 10. Convergence of the enthalpy and spin (∑ S z ) for each structure found by the GA relative to the ground state for the two parameterisations of the LJ + S potential: AFM (left) and FM (right). The color represents the generation in which that structure was found.

Final Structures
The final lowest enthalpy structures are shown in Figure 11. Common neighbour analysis by the software visualisation package Ovito [37] reveals the ground state structure for the AFM parameterisation to be face-centred cubic (FCC) and the ground state structure of the FM system to be body-centred cubic (BCC). It can be seen that as well as finding the crystal structures, the GA was able to correctly align and anti-align the spins for the FM/AFM structures respectively. This illustrates the power of our new magnetic GA to perform global optimisation in both spin and crystal structure spaces simultaneously.

Heusler/Ge Interface
For the CFAS/Ge interface structure search, the calculation was performed on the 12-atom unit cell containing four germanium, four cobalt, two iron, one aluminium and one silicon atom. The GA was allowed to fully optimise the atomic positions. The lattice parameters were constrained to that of the full Heusler alloy, a = b = c = 5.676 Å, since there is very little lattice mismatch to the germanium lattice and it is not expected that this would change much over such a short interface region. The GA was run with a population size of 20 members per generation for a maximum of 100 generations.
An ionic position mutation amplitude of 3 Å was used with a rate of 0.03 mutations per atom. This was chosen to allow on average one atom every two members to mutate, but significantly enough to hop lattice sites. A slightly higher ionic permutation rate of 0.04 was chosen since, for the bulk Heusler, atomic swaps on the lattice sites can be comparable in energy to the ground state. If a similar crystal structure is found, it is likely that exploring this kind of disorder would be a good idea. The spin mutation rate was 0.06 per atom. This was chosen to be higher than the atomic mutation rate because not all atoms would be expected to magnetise, meaning that some mutations would get negated by the local optimisation.
Each member was optimised using the BFGS algorithm to a tolerance of 1 meV. The energies were computed using CASTEP, with a cutoff energy of 650 eV and k-point spacing of 0.04 Å −1 . The PBE functional [22] was used in conjunction with Hubbard U values of 2.1 eV on the d-orbitals of iron and cobalt, in keeping with previous CASTEP studies of the related Heusler alloys [35]. Figure 12 shows the range of enthalpies for CFAS/Ge structures found at each generation of the GA. It can be seen that, on average, the enthalpy of structures in decreases significantly over the first 10 generations, and remains below that of the random structures searched in generation 0 for the remainder of the calculation. This implies that the GA successfully explores more favourable regions of the configuration space. In addition, it can be seen that every generation retains a significant distribution of structures over a range of about 2.5-3 eV. This suggests that the population is remaining significantly diverse and not converging to a single low enthalpy solution. Finally, the inset to Figure 12 shows how the lowest enthalpy value found changes during the calculation. The overall lowest enthalpy structure of the run was found in generation 72. The calculation proceeded for a further 28 generations without finding a lower enthalpy structure and hence this can be taken as a reasonable candidate for the global lowest enthalpy structure. Figure 13 shows the distribution of magnetic structures found, plotted against their energy. Two clusters of structures are found in the energy range which would be thermally accessible during growth. The first of these are ferromagnetic, with a total spin of about 8h/2 and a modulus spin of about 9.9h/2. The second cluster is antiferromagnetic, with a total spin of about 0h/2 and a modulus spin of about 7.5h/2.  Figure 12. Total enthalpy per 12-atom cell of CFAS/Ge structure relative to the overall lowest enthalpy structure found. The lines represent the highest, lowest and mean enthalpies, along with the overall best structure found so far for each generation. The inset shows the best structures found in the 0-0.05 eV range.  Figure 13. Total spin of the Co 2 FeAl 0.5 Si 0.5 (CFAS)/Ge GA structures, plotted against enthalpy of the 12 atom cell, relative to the overall lowest enthalpy structure found. Thermally accessible energies are marked for 575 K (blue dots) and 300 K (red dashes). The color of the points represents the generation in which the structure was found. Figure 14 shows these regions in more detail. It can be seen that each cluster contains a number of structures. These are atomic swaps of the lowest enthalpy structure in each magnetic state. Candidate structures are labelled A1-A3 for the AFM structures and F1-F4 for the FM structures.

Resultant Structures
The lowest enthalpy structure F1, shown in Figure 15 turns out to be a half Heusler alloy, where the X 1 site contains cobalt, the X 2 site contains a vacancy, the Y site contains the iron, silicon and aluminium, and the Z site contains germanium. The F2 structure has the same crystal structure as F1 but with a slightly different spin structure. The F3 and F4 structures have the same crystal structure except one germanium atom is swapped with the silicon atom from the Y site. The results are summarised in Table 2. We also performed a Hirshfeld charge and spin analysis, which showed that in each of these structures the spin is localised on the iron and cobalt, with a spin of around 3h/2 on the iron atoms and 0.5h/2 on the cobalt atoms.  Unlike the full CFAS structure, none of these structures are half-metallic, with a significant number of states around the Fermi energy. It does not appear that the atomic disorder separating F1-4 has a significant effect on the electronic structure of the structures. Figure 16 shows the density of states for bulk CFAS, along with the lowest enthalpy FM and AFM structures found by the GA. It can be seen that, unlike the bulk CFAS structure, these structures are not half metallic, as there is no band gap in the minority spin channel. This provides fundamental insight into how this interface structure might degrade device performance. E-E f (eV) Figure 16. Density of states (DoS) for bulk CFAS (left) and two candidate half-Heuslers: F1 (middle) and A1 (right). The red and blue lines show the DoS for the majority and minority spins respectively relative to the Fermi energy. These are characteristic of the FM and AFM structures discussed in this paper. It can be seen that, unlike bulk CFAS, neither of the half-Heusler structures have a band gap in the minority spin channel.
This illustrates the power of our new magnetic GA in an ab initio context, to fully optimise spin and structure simultaneously, even in the absence of any experimental information on a candidate structure.

Conclusions
We have presented an enhanced GA for the structure prediction of magnetic materials, such that the magnetic and crystal structure may be predicted simultaneously. We have introduced new operations for atomic spin such as crossover, mutation and permutation operations. In addition, the idea of structural fingerprinting using the crystallographic structure factor was extended, based on ideas from magnetic neutron scattering theory, and it was demonstrated that the new fingerprint could identify both similar and distinct crystal and magnetic structures.
To test the magnetic GA, we introduced a novel pair potential V LJ+S which included magnetic-like effects. While this potential does not accurately model any particular real-life material, it provided a computationally efficient way of exploring magnetic materials without the cost of fully quantum calculations. Two parameterisations were given for this potential. The first was parameterised by DFT simulations of Fe dimers in both aligned and anti-aligned configurations. Since the Fe dimers preferred to be aligned, it was expected that this parameterisation would yield a FM structure. This was observed as the lowest energy structure found by the GA was FM FCC, and there was a clear trend in the energies of all the structures searched towards FM structures. The second parameterisation was the same as the first, except for a reversed sign on the exchange-like interaction. This was expected to result in AFM structures since pairs of aligned spins would raise the energy of the system. Indeed, the lowest energy structure discovered by the GA was an AFM BCC structure. In addition, there was a clear trend amongst the other structures towards AFM alignment.
Finally, as an example of studying novel magnetic structures of significant experimental interest, the interface between the Heusler alloy CFAS and n-doped germanium was investigated. Experimentally, the two materials showed significant mixing and, when annealed, a new phase formed at the interface as shown by EDS measurements. The structure of this phase was investigated with the new magnetic GA, and a half-Heusler structure was predicted.
To conclude, we have presented a new GA for structure prediction capable of optimising both magnetic and crystal structures simultaneously.