Designing a new algorithm for crystal structure prediction in magnetic materials

(1)

KTH - Royal Institute of Technology

Master of Science Thesis

Designing a new algorithm for crystal structure prediction in

magnetic materials

Supervisors:

prof. Artem Oganov, Skoltech prof. Anatoly Belonoshko, KTH

Candidate:

Michele Galasso

June 2017

(2)

ii

(3)

Abstract

The USPEX code is a powerful and reliable methodology for crystal structure prediction in solid materials. It consists of an evolutionary algorithm which starts from a population of symmetrically generated structures and improves it at each step by making use of the best individuals. Two different methods are used in order to generate the first population: Topological Random and Symmetric Random.

This thesis extends both of the above methods by implementing generation of magnetic structures in USPEX. For Topological Random, magnetic structures are generated by placing different spin states in different Wyckoff positions of the unit cell, chosen with the right multiplicity. For Symmetric Random, instead, symmetry group decomposition is used in order to split magnetic atoms into two groups to which different spin states are assigned.

For a better fulfillment of this task, also a new algorithm for magnetic structure recognition has been introduced in USPEX. Finally, a computational test has been performed in order to test the new developed method. By confronting statistics on 1000 newly generated structures, the test shows that the implemented initialization improves the quality and the variety of structures if compared to random initialization.

iii

(4)

iv

(5)

Acknowledgements

First and foremost, I would like to thank my supervisor Prof. Artem Oganov for bringing me into his group with such open arms and for giving me constant support throughout the whole project. During these few months have been working in a wonderful environment, I met extraordinary people and I understood what it does really mean to work with science.

I am deeply thankful to the whole USPEX team and in particular to Pavel Bushlanov, for giving me immediate and precious help all the countless times when I was stuck with some problem and for helping me a lot to improve my MATLAB and Python skills.

I also would like to thank Prof. Anatoly Belonoshko for opening to me this thesis opportunity abroad and for giving me important help and feedback during the project. Thanks to him I have been not only working for five months with top level scientists, but I have also experienced the contact with a new and extremely interesting culture and I have become fluent in my fourth language.

Last but not least, I would like to thank my family for making me who I am and for giving me constant support throughout all of the important choices I had to make in my life.

v

(6)

vi

(7)

Introduction

This project consists in designing and implementing a series of new functions for the USPEX code, a powerful evolutionary algorithm for crystal structure prediction in solid materials. In particular, we aim to improve initialization of magnetic structures in order to make searches for new magnetic materials more efficient and robust. Chapter 2 gives a general presentation of the USPEX code and of its working principles. It also contains a number of references to scientific publications in which the various algorithms implemented in USPEX are explained in more detail. Chapter 3, instead, focuses on the available tools for structure generation, Topological Random and Symmetric Random, and explains their working principles. After this preliminary treatment, chapter 4 tackles the problem of magnetic structure generation. It starts with an introduction to the different possible spin states of magnetic atoms and then it explains how we implemented initialization of magnetic structures, first in Topological Random and then in Symmetric Random. In order to check if the symmetrically generated structures are of the desired magnetic type, we have also implemented in USPEX a new algorithm capable of determining the magnetic type of a general structure. This algorithm is introduced and described in chapter 5. Finally, we measured the performance of the newly implemented generation of magnetic structures by means of a computational test. This test, together with its detailed results, is described in chapter 6. It shows that the newly implemented algorithm improves the quality and the variety of magnetic structures if compared to random spin initialization.

In appendix A we report the source code of all the new functions that have been introduced in USPEX with this project. The explanation of the code and references to the different subsections of appendix A are provided

1

(11)

2 Chapter 1. Introduction

throughout the whole text of the thesis.

(12)

Chapter 2

The USPEX code

The USPEX code [1] is an efficient and reliable methodology for crystal structure prediction originally conceived by Artem R. Oganov and Colin W.

Glass at the ETH Zurich in 2004. Its main purpose is the prediction of the stable crystal structure for a given chemical compound at any pressure and temperature condition, without requiring any experimental input. In addition, USPEX predicts for the given compound also a number of low- energy metastable structures and it calculates several other material properties. USPEX is mostly written as a set of MATLAB functions, to which some Python scripts are mixed. This combination of programming languages exists because, at the time of this project, the USPEX developing team was working on the new version 10.1 which will be released completely rewritten in Python code.

The knowledge of the crystal structure is the most important piece of information about crystalline solids, since many physical properties can be calculated from it. However, the problem of crystal structure prediction on fully theoretical basis is an old riddle of science and despite many different approaches have been proposed [2–4], their success rates and their limits of applicability have always revealed very harsh if compared to USPEX [5].

The purpose of any algorithm that attempts to predict the crystal structure of some material is to find the arrangement of atoms in the unit cell which corresponds to a minimum of the Gibbs free energy. The mathematical function that links all possible structures of a particular compound with their Gibbs free energy is called free energy landscape of the material. It is immediately evident that this energy surface is extremely complex and multidimensional and a point by point exploration of the whole landscape is

3

(13)

4 Chapter 2. The USPEX code

not viable with the computing power available today¹.

Unlike its predecessors, which focused mainly on exploring one or a few good regions of the free energy landscape, USPEX tries to scan it as whole looking for the most promising regions and it does it by merging ab initio total energy calculations with a specifically devised evolutionary algorithm.

The term ab initio in physics is used to denote a type of calculation made from well established physical laws, which does not use any empirical as- sumption. An evolutionary algorithm, instead, is an algorithm which selects the best solutions for a given problem by starting from an initial population of possible solutions and by improving this population at each step, simu- lating the natural process of Darwinian evolution. Evolutionary algorithms are very general and they can be applied to a big class of scientific problems.

In particular, through the USPEX code, they have revealed to be successful also in the field of crystal structure prediction.

The success of the evolutionary approach for chemical systems can be explained in terms of the properties of their free energy landscape. As already mentioned the whole exploration of this surface is not viable due to its enormous size, but luckily only a tiny fraction of it corresponds to reasonable physical structures. In addition, due to crystal symmetry, to one good structure in general correspond many equivalent minima of the free energy surface and due to the similarity of good structures most or all of the low energy minima are located in the same region of the landscape. Therefore, the aim of USPEX is to locate these small and promising regions and to con- centrate there the search for the stable crystal structure. Its evolutionary algorithm consists of the following steps:

1. Initialization: the first generation is produced randomly.

2. The structures are relaxed to the nearest local minimum by means of ab initio energy calculations. For this step USPEX can use several first-principles simulation packages, and in this project we have used VASP [6].

3. Among the locally optimized structures, a certain number of the worst ones are rejected and the remaining structures participate in creat- ing the next generation through the action of variation operators. A specified number of the best structures (usually, one) in the current generation survive and compete in the next generation.

1The dimensionality of the free energy surface is 3N + 3, where N is the number of atoms in the unit cell. In detail, this is due to the 3N − 3 degrees of freedom in the atomic positions plus the 6 degrees of freedom which describe the whole unit cell.

(14)

5

4. The new generation is relaxed (by going back to step 2) or the simulation is stopped if some halting criterion is met.

Four variation operators are used in step 3: heredity, atomic permutation, mutation and transmutation. Heredity produces new structures by matching slices of two or more parent structures; these slices are chosen in random directions and with random relative size. Atomic permutation consists in switching the identities of two or more atoms in the unit cell. Mutation consists in the random change of the cell vector of a structure and/or of the atomic positions, while transmutation consists in transmuting a randomly selected atom into another chemical species present in the system.

(15)

6 Chapter 2. The USPEX code

(16)

Chapter 3

Structure generating tools

At the starting point of this project, among the limitations of the USPEX code there was the impossibility to fully take into account magnetic moments when finding the ground state of magnetic materials. A first set of MAT- LAB functions for handling magnetic calculations had already been written, but some important issues still needed to be addressed. First of all, the production of initial structures for the first generation of the evolutionary algorithm did not take magnetic moments into account. Moreover, some criteria needed to be implemented for handling spins during structure relaxation and for determining the magnetic state of a given structure. Our goal is to implement magnetic calculations in USPEX by adapting as much as possible the existing tools and by writing a few new MATLAB functions.

The first issue that we address is the generation of good initial structures for the first generation. For non-magnetic materials, two different methods were already present in USPEX for this purpose, namely Topological Random [7]

and Symmetric Random [8].

3.1 Topological Random

Topological Random is a newly conceived algorithm for crystal structure generation that has been recently implemented in USPEX. It combines data min- ing techniques with the purely mathematical tool of symmetry groups and it does it by using databases of known crystal structures. These databases are then coupled with an algorithm capable of generating new structures not present in any database and corresponding to arbitrary stoichiometries.

More in detail, the Random Topological approach consists of two steps: first the known crystal structures are abstracted leading to a small set of ideal

7

(17)

8 Chapter 3. Structure generating tools

nets, and second the new crystal structures are generated from the ideal nets.

After the abstraction each structure is represented as a graph, made of nodes and edges. Atomic types and actual atomic positions are ignored, therefore many different structures happen to have the same topology (but they may very well correspond to different space groups). Even more general than the concept of topology is the concept of underlying topology of a structure. This last is defined as the usual topological graph where all the nodes with coordination numbers 1 or 2 have been removed. In this way it is possible to obtain a list of underlying topologies for all the crystal structures observed in nature. This is a very important piece of information about crystalline solids since it isolates the underlying topologies chosen by nature and it also includes a huge (and in principle infinite) number of structures with good underlying topologies that, indeed, have never been observed. At this point, only one last piece of information about geometry is needed before proceeding with structure generation. This geometrical information consists in choosing, for each of the underlying topologies, only its most symmetrical spatial representation. In other words, from all the possible mathematical space groups associated with a given topology, only the most symmetric one is considered (namely the space group with the highest number of symmetry elements). Luckily this combination of topological and geometrical information has already been collected in a database used by the TOPOS package [9].

From this database the Topological Random algorithm makes its structure generation, by using symmetry group decomposition in order to associate the nodes of the topological graph with atoms of different type. More in detail, the algorithm is the following:

1. One topology is randomly chosen from the database.

2. The most symmetric space group related to the chosen topology is decomposed in several subgroups. This means that a few less symmetrical space groups are obtained from the initial one by discarding one or more symmetry elements.

3. Wyckoff positions¹for atoms of different type are chosen from different subgroups and their multiplicity is selected according to the desired stoichiometry (for example for a stoichiometry of the type A2B₃, two different Wyckoff positions of multiplicity 2 and 3 are chosen from different subgroups).

1For a detailed explanation of the concept of Wyckoff position, see the next section.

(18)

3.2. Symmetric Random 9

In this way each topology in the database generates an infinite family of structures, each having the same topology but with different group symme- tries. The symmetry of a given structure in each family is a subgroup of the space group which describes the base topology.

3.2 Symmetric Random

Symmetric Random, on the other hand, is a fully theoretical algorithm which does not use any experimental database. It generates each crystal structure according to a mathematical space group, randomly selected among the 230 possible symmetry groups for 3D crystals. Each group is constituted by a set of symmetry operations which link the spatial points of the unit cell with each other, describing in this way the symmetry of the resulting crystal. In the International Tables for Crystallography this symmetry is described in an equivalent way in terms of Wyckoff positions: a set of positions which implicitly contain the symmetry elements of the chosen space group and in this way they specify where the atoms can be found in a crystal with that symmetry. Consider, for example, a space group whose only symmetry element is a symmetry plane at y = 0. The first Wyckoff position of this group will be the arbitrary point (x, y, z) with y 6= 0. This position has multiplicity 2, since a second equivalent point is generated at (x, −y, z) by the action of the symmetry plane. The second Wyckoff position, instead, will be (x, 0, z)and its multiplicity in this case is only 1 since it lies on the symmetry plane and no additional points are generated by its action. The first Wyckoff position that we have seen in this simple example is called general since it does not reside upon any symmetry element of the space group, while the second is called special. Special Wyckoff positions correspond to atoms which lie upon one or more symmetry elements, and for this reason they always have a smaller multiplicity than the general positions. Furthermore, one or more of their coordinates must be fixed in order to force their lying upon symmetry elements.

In the Symmetric Random algorithm, after the creation of a cell shape consistent with the chosen space group, an atom is randomly placed on a general Wyckoff position and it is duplicated in other points of the unit cell through the action of space group symmetry operators. If after this step two symmetrically related atoms happen to be too close to each other, the algorithm merges them in a single atom placed on a special Wyckoff position which is obtained by averaging the coordinates of the two initial atoms. In order to get the desired stoichiometry, some atoms may have to

(19)

10 Chapter 3. Structure generating tools

be placed on special Wyckoff positions with the right multiplicity. Once the right stoichiometry has been reached for the first atomic type, the algorithm repeats the same process for the subsequent atomic types. It can also happen that at some point two atoms of different type are placed in positions that are too close to each other, and in this case the structure in process of generation is immediately discarded without relaxation. Finally, to allow the possibility of subsequent symmetry break by means of the variation operators, a small random displacement is applied to all the atomic positions during structure relaxation. This procedure for the construction of the first generation has the advantage of providing very diverse initial structures and to better cover the search space.

In USPEX it is possible to generate structures with Topological Random and Symmetric Random also in subsequent generations, in order to provide a few new symmetric structures at each step of the evolutionary algorithm.

(20)

Chapter 4

Generation of magnetic structures

The main idea behind this chapter is to adapt the already existing tools Topological Random and Symmetric Random to the generation of magnetic structures. For this purpose, we used two different methods in order to es- tablish a geometrical arrangement among spin states in the newly generated structures. In Topological Random, different spin states of the same atomic species are placed in different Wyckoff positions as they were atoms of different type, while in Symmetric Random a symmetry element of the space group associated to the structure is broken and spin states are placed according to the remaining symmetry elements. Before talking explicitly about magnetic initialization, we will introduce the different spin states that have been considered in this work.

4.1 Spin states

Since magnetism is generated by electron spin, the number of unpaired electrons in a specific atom or compound indicates how magnetic that system is. Among the different atomic species present in nature, transition metals are those which tend to have a large number of unpaired electrons. In this work we have considered only the first two rows of transition metals (see fig.

4.1 [10]) as atoms of magnetic type by default. The reason for this choice is that the transition metals with greater atomic number, due to the chemical nature of their orbitals, almost never give rise to materials with magnetic properties.

In addition to the ordinary spin up and spin down states we have also 11

(21)

12 Chapter 4. Generation of magnetic structures

Figure 4.1: Periodic Table of the Elements

(22)

4.1. Spin states 13

Figure 4.2: Low-spin crystal field diagram

Figure 4.3: High-spin crystal field diagram

considered high-spin and low-spin configurations of transition metals, which arise from the different ∆0 splitting of the d orbitals. If the energy cost of placing an electron into an already singly occupied orbital is less than ∆0 a low-spin state will occur, otherwise we will have a high-spin state (see figs.

4.2 and 4.3 [11]). Moreover, if these two energy costs are comparable, some atoms will prefer a high-spin state and some other will prefer a low-spin state resulting in a low/high-spin mixed configuration.

(23)

4.2 Initialization

The input parameters of USPEX are specified in the INPUT.txt file. In particular the optional parameter magRatio allows the user to specify a vector with seven entries, which contains the relative abundances of different spin polarized structures in the initial population. The meaning of each entry in the magRatio vector are the following:

1. NM – non-spin-polarized

2. FM-LS – low-spin ferromagnetic 3. FM-HS – high-spin ferromagnetic 4. AFM-LS – low-spin anti-ferromagnetic 5. AFM-HS – high-spin anti-ferromagnetic

6. FM-HSLS – low/high-spin mixture ferromagnetic 7. AFM-HSLS – low/high-spin mixture anti-ferromagnetic

4.2.1 Topological Random

The Topological Random algorithm generates structures through the MAT- LAB function random_topology.m and a set of other functions connected to it. It takes as input parameter the vector numIons, which specifies the number of different atomic species in the unit cell. In order to generate a magnetic structure in which atoms with different spin states are treated as different atomic types, we substituted the input parameter numIons in ran- dom_topology.mwith the new vector numIonsMagnetic. This vector has the same structure of numIons, but it lists also the abundances of every particular spin state for each atomic type. For example, if our unit cell is Fe3C1

and in it 2 iron atoms occur with high-spin up and 1 with low-spin up, the numIonsvector will be 3 1 and the numIonsMagnetic vector will be 2 1 1.

The vector numIonsMagnetic is created by the MATLAB function mag- netize_numIons.m, together with the vectors magneticFlag and the variable magType. The complete code of this function is reported in appendix A.1.

The vector magneticFlag has the same size of numIonsMagnetic and it indicates the spin magnitude and direction for each slot in numIonsMagnetic, while magType is a simple number which indicates the magnetic type of each initialized structure according to table 4.1.

(24)

4.2. Initialization 15

magType meaning

1 NM

-2 FM-LS

2 FM-HS

-3 AFM-LS

3 AFM-HS

4 FM-HSLS

5 AFM-HSLS

Table 4.1: Possible values of the variable magType

In magnetize_numIons.m a random seed is used in order to select which type of magnetic structure has to be generated, according to the user-defined abundances memorized in the vector magRatio. For the initialization of high- spin (HS) and low-spin (LS) states we used, respectively, a magnetic moment of 4µB and 1µB, where µB is the Bohr magneton¹. Unless otherwise stated, all the magnetic moments in this work are expressed in units of the Bohr magneton. More in detail, the spin initialization for each magnetic type is done in the following way:

NM All the atoms are assigned spin 0.

FM-LS All the magnetic atoms are assigned spin +1, while the non-magnetic atoms are assigned spin 0.

FM-HS All the magnetic atoms are assigned spin +4, while the non-magnetic atoms are assigned spin 0.

AFM-LS The magnetic atoms of each type are split in two halves (in pres- ence of an odd stoichiometry, one of the two halves will have an element more than the other). The first half is assigned spin +1 and the second half spin −1. The non-magnetic atoms are assigned spin 0.

AFM-HS Same as in AFM-LS, but here the two halves are assigned, respectively, spins +4 and −4.

FM-HSLS The magnetic atoms of each type are split in two halves. The first half is assigned spin +1 and the second half spin +4. The non- magnetic atoms are assigned spin 0.

1In SI units, 1 µB≈ 9.274 × 10⁻²⁴J T⁻¹

(25)

AFM-HSLS The magnetic atoms of each type are split (if possible) in four groups, that are assigned with the four different spin states. If there are not enough atoms of the same type to form four groups but more magnetic atoms of another species are present in the unit cell, then these last atoms are used to fill in the missing polarizations. As always, the non-magnetic atoms are assigned spin 0.

If the magnetic atoms are not enough to set up an AFM configuration, a warning message is displayed and a FM configuration is chosen. The impor- tance of the warning message lies in the fact that this change of magnetic configuration can be a signal that the original proportions in the vector magRatiowill not be fully respected by the generated structures. If no magnetic atoms are present at all, a NM configuration is chosen and a different warning message is displayed. The pictures in figure 4.4 are examples of different magnetic structures created by Random Topology.

The so obtained variables numIonsMagnetic, magneticFlag and magType are then used by the function initialize_magMom.m (listed in appendix A.2). This simple function puts them together in a new vector called mag- mom_ions, which has length equal to the total number of atoms in the unit cell (call it N) plus one. The first slot in magmom_ions is nothing but the old variable magType, while the other N slots contain the value of the magnetic moments for each atom in the unit cell. This vector will then be used by subsequent MATLAB functions in order to pass the magnetic moments to the VASP package. Lastly, the VASP code will take care of the ab initio magnetic calculations and it will determine the magnetic moment of each atom after structure relaxation.

4.2.2 Symmetric Random

For Symmetric Random the trick of considering different spin states as different atomic types could not be used, since it would have considerably increased noise in the generating algorithm and it would not have produced structures of good quality. Therefore we decided to leave the structure generating algorithm as it is also for magnetic materials and then to place spins on magnetic atoms by breaking one of the symmetry elements of the structure’s space group.

The MATLAB function getSpinGuide.m, which is reported in appendix A.3, does the task of dividing magnetic atoms in the unit cell into two groups, that are then initialized with different spin states by initialize_magMom.m.

This function makes use of the two variables numSites and Operation, given

(26)

(a) NM

(b) FM-LS (c) FM-HS

(27)

(d) AFM-LS (e) AFM-HS

(f) FM-HSLS (g) AFM-HSLS

Figure 4.4: Structures of different magnetic types generated by Random Topology for the supercell Fe12C4. The small brown spheres denote the carbon atoms, while for the iron atoms different colours are used for different spin states.

(28)

in output by the Symmetric Random algorithm. The vector numSites has length equal to the number of Wyckoff positions chosen for the atoms in the unit cell, and each of its slots specifies the multiplicity of the Wyckoff position corresponding to the slot’s index. For example, if for the system Mn8Bi8 we have a numSites vector equal to 4 4 8, this means that the 8 manganese atoms are disposed in two different Wyckoff positions each of multiplicity 4, while all of the bismuth atoms are disposed in a unique Wyckoff position of multiplicity 8. The variable Operation, instead, is a collection of 3x4 matrices which describe all the symmetry operators of the structure’s space group. These 3x4 matrices are composed of a 3x3 matrix of rotation and a 3x1 vector of translations along the Cartesian axes.

In general, if G is a space group, it can be written in terms of its symmetry elements as

G = {I, g₂, g₃, . . . , g_M} (4.1) where I is the identity and the number M of operators is called rank of the space group. Every space group is closed, since if we take any couple of operators gi, gj ∈ G, the operator given by their product still lies in G. If we remove one of the gi from the space group in order to obtain a less symmetrical subgroup, we have to remove also a number of other elements in order to make sure that the new group is also closed and that it is not possible to obtain the removed element by any combination of the remaining operators. In this way we obtain a subgroup of rank M/n, where nis the rank of the symmetry element that has been broken². However, the choice of the particular subgroup is not trivial since usually there are many alternatives and it is hard to write an algorithm which works for the general case, independently of the particular space group chosen. For this reason we implemented group decomposition in a slightly different way, which is still able to produce a random subgroup from a general space group and it ensures that the broken element is of rank 2. Later this will allow us to split each Wyckoff position into two groups of coordinates of equal size, one for each spin state. The implemented algorithm is the following:

1. Two random operators are chosen from the space group G.

2. The operators are multiplied by themselves and with each other in each possible combination, in order to verify if the chosen subgroup is closed. If new elements arise, they are included in the group and the operation is repeated until a closed group is found.

2The rank of a symmetry operator is defined as the number of times that operator has to be applied to a position in the unit cell in order to obtain back the starting point.

(29)

3. If the number of elements in the subgroup is M/2 the symmetry elements are saved, otherwise they are discarded and the algorithm restarts from step 1.

In order to improve performance and to avoid the infinite cycle, we excluded from this search the space groups with rank M ≤ 3, since for all of them the immediately less symmetric subgroup is P1 ³ which contains only the identity operator. In these cases atomic positions for different spin states are just divided at random. In addition, we noticed that the space groups P23, F23 and I23, which have 6 symmetry elements, cannot be reduced in any subgroup with 3 elements. Therefore also in this case Wyckoff positions are randomly split.

For all the other groups where an operator of rank 2 can successfully be broken, the operators of the obtained subgroup are used in order to split Wyckoff positions. For each slot in numSites, we take the coordinates of the first atom corresponding to the selected Wyckoff position and we apply to it all the operators of the subgroup. Usually after this operation we obtain half of the atomic coordinates relative to that Wyckoff position. For example if a mirror plane has been broken, we obtain only the atomic positions on one side of the plane and we miss the operators which give the positions on the other side. The obtained coordinates are then marked to be assigned with a spin value, while the remaining ones will be assigned with a different spin value. However, it can be the case that the operators of the subgroup are enough to obtain all the atomic coordinates corresponding to the chosen Wyckoff position. This happens when we have a special Wyckoff position and the atomic coordinates all lie on the symmetry element that has been removed. In this case there is no way to symmetrically split positions in two groups, therefore the spin states are assigned randomly.

Got the vector spinGuide, the function initialize_magMom.m does the task of assigning spin values to the two groups of atomic positions depending on the structure’s magnetic type. For the case AFM-LS, the two groups will be given spins +1 and −1, for the case AFM-HS they will be given spins +4 and −4, and for the case FM-HSLS, they will be given spins +1 and +4. More complex is the case AFM-HSLS, where we would have needed four groups of symmetrically related coordinates for each Wyckoff position in order to assign the four possible spin values. However, splitting each Wyckoff position into four groups would have been problematic, given that in general Wyckoff positions have relatively low multiplicity and many space groups would not

3Here and in the following, we refer to space groups according to their international short symbols.

(30)

have allowed this splitting, increasing the number of exceptions. So what we did in this case is to use the two symmetrically obtained groups for spin-up and spin-down states and to assign HS and LS at random, keeping as always an even balance between them. In general, we believe that the initialization of HS and LS states is not as important as the initialization of spin-up and spin-down states, since HS and LS can more easily transform into each other during structure relaxation. Different types of magnetic structures generated with Symmetric Random are reported in figure 4.5.

The function getSpinGuide.m makes use of the simple function ap- plySymmetry.m, used to let symmetry operators act on an atomic positions.

The whole code of this function is reported at the end of appendix A.3.

(31)

(a) NM

(b) FM-LS (c) FM-HS

(32)

(d) AFM-LS (e) AFM-HS

(f) FM-HSLS (g) AFM-HSLS

Figure 4.5: Structures of different magnetic types generated by Symmetric Random for the unit cell Mn8Bi8. The big purple spheres denote the bismuth atoms, while for the manganese atoms different colours are used for different spin states.

(33)

(34)

Chapter 5

Identification of magnetic types

Two main reasons make it important to implement in USPEX identification of magnetic type for a general structure. First, among the several output files that USPEX creates there is a file named Individuals_magmoment, which collects the spin values of all the atoms in each structure after relaxation and the new magnetic types of the structures. This last information is not trivial to obtain, given the possibility that the magnetic type changes during relaxation. The second reason, instead, is about the use of magnetic type identification for improving the algorithm which generates the initial magnetic structures. This is very important in order to give consistency to the whole algorithm which handles magnetic materials in USPEX and it will be recalled and better explained at a later stage.

The MATLAB function magTypeIdentify.m, reported in appendix A.4, takes care of investigating the magnetic type of a general structure. First of all, two thresholds are defined in this function in order to distinguish between the states spin zero, low-spin and high-spin. In order to get some reasonable values for these parameters, we run a few tests in which a population of 50 magnetic structures for the unit cell Fe12C4was generated and then relaxed.

Once obtained the spin values after relaxation, we did a statistics of the final magnetic moments of all the magnetic atoms (fig. 5.1). From this statistics we chose to consider as spin zero all the spins with absolute value below 0.2, as low-spin states all the spins with absolute value between 0.2 and 1.2 and as high-spin states all the spins with absolute value greater than 1.2.

More in detail, our algorithm for magnetic type identification consists of the following procedure:

• The spins in the unit cell are screened. If all of them have values under the defined threshold 0.2 the structure is classified as non-magnetic,

25

(35)

26 Chapter 5. Identification of magnetic types

Figure 5.1: Statistics of the final magnetic moments for 50 structures of Fe12C4 after a 5 step VASP relaxation.

otherwise more steps are needed in order to determine its magnetic type.

• For each magnetic atom in the unit cell, we examine its spin value in relation to the spins of all its nearest neighbours of magnetic type (see fig. 5.2 [12]). A variable criterion is created and initialized to zero.

Then its value is incremented by one every time that two magnetic atoms with concordant spins happen to be nearest neighbours and it is decremented by one every time that two magnetic atoms with discordant spin happen to be nearest neighbours. If at the end of this calculation we have criterion ≥ 0 the structure is considered FM, otherwise it is considered AFM.

• The number of high-spin states in the unit cell are counted. If they are less than 1/3 of all the magnetized atoms the structure is considered to be in LS state; if they are more than 2/3 of all the magnetized atoms the structure is considered to be in HS state; otherwise the structure is classified as a mixed HSLS state.

As discussed in chapter 4, the initialization of magnetic structures follows

(36)

5.1. Neighbour determination 27

Figure 5.2: Nearest and next-nearest neighbours are displayed for the central atom of a body centred cubic lattice (bcc).

a slightly different procedure in which the neighbours of each atom are not taken into account. This procedure has no important consequences for FM structures, where the spins are all concordant and the variable criterion will have the maximum possible value, but some inconsistencies can arise in the case of AFM structures. When an AFM structure is initialized, the magnetic atoms in the unit cell are just divided into two groups (for spin up and spin down states) according to some symmetry criterion. In some cases, a structure with two separate spin-up and spin-down domains can be generated and the variable criterion for it can be greater than zero even if the user wanted that structure to be AFM. For this reason, all the newly generated structures are checked with magTypeIdentify.m just after their initialization and if the automatically identified magnetic type does not correspond to the magnetic type that wanted to be generated, the structure is discarded and the system tries to generate a new one of the same magnetic type.

5.1 Neighbour determination

A non trivial step for implementing magnetic type identification has been to identify the nearest neighbours for each magnetic atom in the unit cell, in a general structure. In order to solve this problem, we wrote the MAT- LAB function determineNeighbors.m, reported in appendix A.5. This function, in turn, uses the two Python scripts coordination_numbers.py and

(37)

28 Chapter 5. Identification of magnetic types

dists.pythat have been connected to the MATLAB part of USPEX by two simple adapters. The first script already existed in the previous versions of USPEX: it uses a powerful empirical model in order to calculate the effective coordination number CN of each atom in the unit cell [13]. The second script, instead, has been written on the basis of the first one and it returns a matrix which contains all the distances between pairs of atoms in the unit cell, including the atoms in the 26 copies of the cell around itself in three di- mensions. The code of these two Python scripts is also reported in appendix A.5.

In order to restrict our attention only to atoms of magnetic type, the function determineNeighbors.m starts with a screening of the unit cell which removes all the non-magnetic atoms before passing the necessary information to the Python scripts. In this way we get in output from the scripts the magnetic coordination number CN M and the inter-atomic distances only between atoms of magnetic type. Obtained this output, the main function determineNeighbors.m considers as neighbours of each magnetic atom all the atoms within a distance

R =r 3V³ _cell

4πN_M (C_{N M}+ 1) (5.1)

where Vcell is the volume of the unit cell and NM is the total number of magnetic atoms. In other words, R is derived from the following relation

4

3πR³ = V_cell

N_M (CN M+ 1) (5.2)

as the radius of a sphere whose volume contains, on average, CN M+ 1atoms given the density of magnetic atoms in that particular cell. These atoms are nothing but the central atom plus its CN M nearest neighbours of magnetic type. Given that the density of atoms in the unit cell is in general non uniform, this is not an exact solution and it can happen that more or less atoms than CN M are considered to be nearest neighbours. However, this solution has proven to be useful for our purposes and not too computationally expensive.

(38)

Chapter 6

Test of the new algorithm

In order to investigate the performance of our algorithm, we performed a computational test by generating 1000 crystal structures for the unit cell Mn8Bi8, half with Topological Random and half with Symmetric Random.

We relaxed each structure to the nearest local minimum with a 5 step VASP relaxation and then we performed energy and spin statistics.

As a control group for this test, we repeated the same simulation with random spin initialization and with no identification of magnetic type. In this case, each structure has been generated with Topological Random or Symmetric Random as it was non-magnetic and the spins on the Mn atoms have been assigned by giving them a random value in the set {−4, −1, 0, 1, 4}.

6.1 Energy statistics

In the experimental group, we found that 723 structures were correctly relaxed while the remaining were discarded during relaxation, due to their diverging energy or to the non-fulfilment of some distance constraint. The average energy of the relaxed structures is

E_E = −96.9708 eV (6.1)

In the control group, instead, we found that only 548 structures correctly completed relaxation while the remaining were discarded. This high per- centage of structures which failed relaxation is an index of a lower quality of the initial structures in this case, given that VASP parameters have been left equal to those of the experimental group. The average energy of the relaxed structures in the control group is

EC = −97.1621 eV (6.2)

29

(39)

30 Chapter 6. Test of the new algorithm

The two average energies have almost the same value, however from figure 6.1 we see that the structures in the experimental group exhibit a variety of energies that is bigger than the variety of energies in the control group. This information is even more evident if we compare the two energy variances of the experimental group and of the control group:

V (E_E) = 2.6680 eV² V (E_C) = 1.2858 eV² (6.3) Therefore, in the experimental case, very good or very bad structures have a greater probability of being generated and this enhances the possibility of obtaining a structure closer to the ground state of the magnetic material.

6.2 Spin statistics

The statistics of spin values in the relaxed structures for the experimental group and for the control group are displayed in figure 6.2. From these statistics we can conclude that the Mn8Bi8 cell is a system in which the magnetic atoms prefer the high-spin configuration, since low-spin states are almost absent. We also notice that, despite the fact that more atoms have been initialized to spin zero in the control group (1/5 versus 1/7 of the total magnetic atoms), the abundance of spin zero states is much higher in the experimental group. This is due to the absence of completely non-magnetic structures in the control group, where the more abundant but sparse spin zero atoms are influenced by their spin-polarized neighbours which in turn make them assume a non-zero spin during relaxation.

The higher amount of spin zero states in the relaxed structures is a good feature of the experimental treatment, which has the ability to preserve some ordered non-magnetic configurations that would have been lost in the case of random spin initialization.

(40)

6.2. Spin statistics 31

(a) Experimental group

(b) Control group

Figure 6.1: Energy statistics for relaxed structures in the experimental group and in the control group.

(41)

32 Chapter 6. Test of the new algorithm

(a) Experimental group

(b) Control group

Figure 6.2: Spin statistics for relaxed structures in the experimental group and in the control group.

(42)

Chapter 7

Conclusions

The work presented in this thesis, despite the impression of completeness that the reader could have, is only a small step in the full implementation of magnetism in USPEX. During the project we focused mostly on the designing and coding part and, as it can be seen from the appendices, the quality of the new code has become better and better over time. However, the computational test described in chapter 6 is the first test of the new algorithm with a consistent amount of structures and more extensive tests are needed in order to refine all the details and to tune the newly introduced parameters.

In addition, magnetism needs to be implemented in the subsequent generations of the evolutionary algorithm by extending the existent variation operators. New rules are needed in order to decide how to assign spin values to child structures from their parents and how to implement purely magnetic variation operators such as spin mutation (the creation of a new structure by randomly changing the magnetic type of an existent structure). Not less important is the translation of the new MATLAB functions into Python for the upcoming USPEX 10.1, which will lead to a consistent improvement in the efficiency of the code.

We believe that, when all the conditions to take full advantage of our work will be reached, the algorithm described in this thesis will hopefully lead to the discovery of new magnetic properties of materials and it will contribute to the already big list of successes of the ambitious USPEX project.

33

(43)

34 Chapter 7. Conclusions

(44)

Bibliography

[1] Artem R Oganov and Colin W Glass. Crystal structure prediction using ab initio evolutionary techniques: Principles and applications. The Journal of chemical physics, 124(24):244704, 2006.

[2] J Pannetier, J Bassas-Alsina, J Rodriguez-Carvajal, and V Caignaert.

Prediction of crystal structures from crystal chemistry rules by simu- lated annealing. Nature, 346(6282):343–345, 1990.

[3] Roman Martoňák, Alessandro Laio, Marco Bernasconi, Chiara Ceri- ani, Paolo Raiteri, Federico Zipoli, and Michele Parrinello. Simula- tion of structural phase transitions by metadynamics. Zeitschrift für Kristallographie-Crystalline Materials, 220(5/6):489–498, 2005.

[4] Stefan Goedecker. Minima hopping: An efficient search method for the global minimum of the potential energy surface of complex molecular systems. The Journal of chemical physics, 120(21):9911–9917, 2004.

[5] Colin W Glass, Artem R Oganov, and Nikolaus Hansen. Uspexevo- lutionary crystal structure prediction. Computer Physics Communica- tions, 175(11):713–720, 2006.

[6] G Kresse. Comput. matter sci. 6, 15 (1996);(d) kresse, g., and furth- muller. Phys. Rev. B, 54:11–169, 1996.

[7] Pavel V Bushlanov, Vladislav A Blatov, and Artem R Oganov.

Topology-based crystal structure generator, (unpublished).

[8] Andriy O Lyakhov, Artem R Oganov, Harold T Stokes, and Qiang Zhu.

New developments in evolutionary structure prediction algorithm uspex.

Computer Physics Communications, 184(4):1172–1182, 2013.

[9] Vladislav A Blatov, Alexander P Shevchenko, and Davide M Proser- pio. Applied topological analysis of crystal structures with the program package topospro. Crystal Growth & Design, 14(7):3576–3586, 2014.

35

(45)

36 Bibliography

[10] http://sciencenotes.org/printable-periodic-table.

[11] G.L. Miessler and D.A. Tarr. Inorganic Chemistry. Pearson Prentice Hall, 2004.

[12] http://www.physics-in-a-nutshell.com/article/12.

[13] Rudolf Hoppe. Effective coordination numbers (econ) and mean fictive ionic radii (mefir). Zeitschrift für Kristallographie-Crystalline Materials, 150(1-4):23–52, 1979.

(46)

Appendix A

MATLAB and Python functions

A.1 magnetize_numIons.m

1 f u n c t i o n [ n u m I o n s M a g n e t i c , m a g n e t i c F l a g , m a g T y p e ] = m a g n e t i z e _ n u m I o n s ( n u m I o n s )

2

3 g l o b a l O R G _ S T R U C

4

5 % n u m I o n s M a g n e t i c m e m o r i z e s LS - up , LS - down , HS - up and HS - do w n as d i f f e r e n t a t om t y p e s

6 % m a g n e t i c F l a g c o l l e c t s i n f o r m a t i o n a b o u t the s pi n d i r e c t i o n and m a g n i t u d e of e v e r y s lo t in n u m I o n s M a g n e t i c

7 % d e v e l o p e r - G a l a s s o M i c h e l e

8

9 a t o m T y p e = O R G _ S T R U C . a t o m T y p e ;

10 m a g R a t i o = O R G _ S T R U C . m a g R a t i o ;

11 r a n d S e e d = r an d() ;

12 m a g T y p e = 1; % d e f a u l t = NM

13

14 for i = 2:l e n g t h( m a g R a t i o )

15 m a g R a t i o ( i ) = sum( m a g R a t i o ( i -1: i ) ) ;

16 end

17 18 % NM

19 if r a n d S e e d <= m a g R a t i o (1)

20 n u m I o n s M a g n e t i c = n u m I o n s ;

21 m a g n e t i c F l a g = z e r o s (1 ,l e n g t h( n u m I o n s M a g n e t i c ) ) ;

22 23 % FM

24 e l s e i f r a n d S e e d > m a g R a t i o (1) && r a n d S e e d <= m a g R a t i o (3)

37

(47)

38 Appendix A. MATLAB and Python functions

25 for i = 1:l e n g t h( n u m I o n s )

26 if ( a t o m T y p e ( i ) >= 21 && a t o m T y p e ( i ) <= 30) || ( a t o m T y p e ( i )

>= 39 && a t o m T y p e ( i ) <= 48) % if the a to m is m a g n e t i c

27 n u m I o n s M a g n e t i c ( i ) = n u m I o n s ( i ) ;

28 m a g n e t i c F l a g ( i ) = 1;

29 e ls e % if the a t om is non m a g n e t i c

30 n u m I o n s M a g n e t i c ( i ) = n u m I o n s ( i ) ;

31 m a g n e t i c F l a g ( i ) = 0;

32 end

33 end

34

35 % FM - LS

36 if r a n d S e e d > m a g R a t i o (1) && r a n d S e e d <= m a g R a t i o (2)

37 m a g T y p e = -2;

38

39 % FM - HS

40 e ls e

41 m a g T y p e = 2;

42 for i = 1:l e n g t h( m a g n e t i c F l a g )

43 m a g n e t i c F l a g ( i ) = 4* m a g n e t i c F l a g ( i ) ;

44 end

45 end

46

47 if sum(abs( m a g n e t i c F l a g .* n u m I o n s M a g n e t i c ) ) == 0

48 d is p(’ No m a g n e t i c a t o m s in the u n it ce l l . I m p o s s i b l e to set up m a g n e t i c c a l c u l a t i o n s . ’) ;

49 m a g T y p e = 1;

50 end

51

52 % AFM

54

55 j = 1; % i n i t i a l i z e s the c o u n t e r

56

57 % t hi s s w i t c h e r a l l o w s to ev e n as mu c h as p o s s i b l e the s p i n s up and do w n in the AFM ca s e .

58 % If for e x a m p l e we h av e Fe3Ni3 , in the f i r s t sl o t it wi l l set 2 s p i n s up and 1 do w n and in the s e c o n d 2 d ow n and 1 up .

59 s w i t c h e r = tr u e ;

60

63

64 if n u m I o n s ( i ) == 1 && s w i t c h e r == tr u e % ca s e w he n t h e r e is on l y 1 m a g n e t i c at o m in the sl o t

65 n u m I o n s M a g n e t i c ( j ) = 1;

66 m a g n e t i c F l a g ( j ) = 1;

(48)

A.1. magnetize_numIons.m 39

67 j = j +1;

68 end

69

70 if n u m I o n s ( i ) == 1 && s w i t c h e r == f a l s e

72 m a g n e t i c F l a g ( j ) = -1;

73 j = j +1;

74 end

75

76 if n u m I o n s ( i ) ~= 1 && s w i t c h e r == t r ue

77 n u m I o n s M a g n e t i c ( j ) = c ei l( n u m I o n s ( i ) /2) ; % h al f ( or h al f +1) sp i n up

78 m a g n e t i c F l a g ( j ) = 1;

79 j = j +1;

80 n u m I o n s M a g n e t i c ( j ) = f l o o r( n u m I o n s ( i ) /2) ; % the r e m a i n i n g sp i n d o wn

81 m a g n e t i c F l a g ( j ) = -1;

82 j = j +1;

83 end

84

85 if n u m I o n s ( i ) ~= 1 && s w i t c h e r == f a l s e

86 n u m I o n s M a g n e t i c ( j ) = f l o o r( n u m I o n s ( i ) /2) ; % ha l f ( or half -1) sp i n up

87 m a g n e t i c F l a g ( j ) = 1;

88 j = j +1;

89 n u m I o n s M a g n e t i c ( j ) = c ei l( n u m I o n s ( i ) /2) ; % the r e m a i n i n g sp i n d o wn

90 m a g n e t i c F l a g ( j ) = -1;

91 j = j +1;

92 end

93 s w i t c h e r = ~ s w i t c h e r ;

94 el s e % if the at o m is non m a g n e t i c

95 n u m I o n s M a g n e t i c ( j ) = n u m I o n s ( i ) ;

96 m a g n e t i c F l a g ( j ) = 0;

97 j = j +1;

98 end

99 end

100

101 % AFM - LS

102 if r a n d S e e d > m a g R a t i o (3) && r a n d S e e d <= m a g R a t i o (4)

103 m a g T y p e = -3;

104

105 % AFM - HS

106 el s e

107 m a g T y p e = 3;

108 for i = 1:l e n g t h( m a g n e t i c F l a g )

109 m a g n e t i c F l a g ( i ) = 4* m a g n e t i c F l a g ( i ) ;

110 end

111 end

(49)

112

113 if l e n g t h(f in d( m a g n e t i c F l a g .* n u m I o n s M a g n e t i c ) ) == 1

114 d is p(’ O nl y one m a g n e t i c at o m in the un i t cell , i m p o s s i b l e to set up an AFM s t a t e . Set as FM i n s t e a d . ’) ;

115

116 if m a g T y p e == -3

117 m a g T y p e = -2;

118 e ls e

119 m a g T y p e = 2;

120 end

121

122 e l s e i f sum(abs( m a g n e t i c F l a g .* n u m I o n s M a g n e t i c ) ) == 0

124 m a g T y p e = 1;

125 end

126

127 % FM - HS L S

129

131 s w i t c h e r = tr u e ; % i n i t i a l i z e s the sw i t c h e r , to e ve n the n u m b e r of HS v e r s u s LS

132 m a g T y p e = 4;

133

136

137 if n u m I o n s ( i ) == 1 && s w i t c h e r == tr u e % t r e a t s the c a se wh e n t h e r e is on l y 1 m a g n e t i c at o m in the sl o t

139 m a g n e t i c F l a g ( j ) = 1;

140 j = j +1;

141 end

142

143 if n u m I o n s ( i ) == 1 && s w i t c h e r == f a l s e

145 m a g n e t i c F l a g ( j ) = 4;

146 j = j +1;

147 end

148

149 if n u m I o n s ( i ) ~= 1 && s w i t c h e r == tr u e

150 n u m I o n s M a g n e t i c ( j ) = c e il( n u m I o n s ( i ) /2) ; % ha l f ( or h a lf +1) LS

151 m a g n e t i c F l a g ( j ) = 1;

152 j = j +1;

153 n u m I o n s M a g n e t i c ( j ) = f l o o r( n u m I o n s ( i ) /2) ; % the r e m a i n i n g HS

(50)

A.1. magnetize_numIons.m 41

154 m a g n e t i c F l a g ( j ) = 4;

155 j = j +1;

156 end

157

158 if n u m I o n s ( i ) ~= 1 && s w i t c h e r == f a l s e

159 n u m I o n s M a g n e t i c ( j ) = f l o o r( n u m I o n s ( i ) /2) ; % ha l f ( or half -1) LS

160 m a g n e t i c F l a g ( j ) = 1;

161 j = j +1;

162 n u m I o n s M a g n e t i c ( j ) = c ei l( n u m I o n s ( i ) /2) ; % the r e m a i n i n g HS

163 m a g n e t i c F l a g ( j ) = 4;

164 j = j +1;

165 end

166 s w i t c h e r = ~ s w i t c h e r ;

167 el s e % if the at o m is non m a g n e t i c

168 n u m I o n s M a g n e t i c ( j ) = n u m I o n s ( i ) ;

169 m a g n e t i c F l a g ( j ) = 0;

170 j = j +1;

171 end

172 end

173

174 if l e n g t h(fi n d( m a g n e t i c F l a g .* n u m I o n s M a g n e t i c ) ) == 1

175 di s p(’ On l y one m a g n e t i c a t om in the u n it cell , i m p o s s i b l e to set up a m i x e d FM - H SL S s t a t e . Set as FM - LS i n s t e a d . ’) ;

176 m a g T y p e = -2;

177 end

178

180 di s p(’ No m a g n e t i c a t o m s in the un i t c e ll . I m p o s s i b l e to set up m a g n e t i c c a l c u l a t i o n s . ’) ;

181 m a g T y p e = 1;

182 end

183

184 % AFM - HS L S

187 m a g T y p e = 5;

188

>= 39 && a t o m T y p e ( i ) <= 48) % if the at o m is m a g n e t i c

191

192 w h i l e j <= sum( n u m I o n s (1: i ) )

193 if e x i s t(’ n u m I o n s M a g n e t i c ’) == 0

194 n u m I o n s M a g n e t i c = [];

195 end

196

197 k = 1 + mod ( j -1 ,4) + 4*( i -1) ;

(51)

198

199 if l e n g t h( n u m I o n s M a g n e t i c ) < k % c r e a t e s the sl o t in n u m I o n s M a g n e t i c

200 n u m I o n s M a g n e t i c ( k ) = 0;

201 end

202

203 if mod ( k ,4) == 1

204 m a g n e t i c F l a g ( k ) = 1;

205

206 e l s e i f mod ( k ,4) == 2

207 m a g n e t i c F l a g ( k ) = -4;

208

209 e l s e i f mod ( k ,4) == 3

210 m a g n e t i c F l a g ( k ) = -1;

211

212 e l s e i f mod ( k ,4) == 0

213 m a g n e t i c F l a g ( k ) = 4;

214 end

215

216 n u m I o n s M a g n e t i c ( k ) = 1 + n u m I o n s M a g n e t i c ( k ) ;

217

218 j = j +1;

219 end

220 e ls e % if the a t om is not m a g n e t i c

221 k = 1 + mod ( j -1 ,4) + 4*( i -1) ;

222 n u m I o n s M a g n e t i c ( k ) = n u m I o n s ( i ) ;

223 m a g n e t i c F l a g ( k ) = 0;

224 j = j + n u m I o n s ( i ) ;

225 end

226 end

227

228 if l e n g t h(f in d( m a g n e t i c F l a g .* n u m I o n s M a g n e t i c ) ) == 1

229 d is p(’ O nl y one m a g n e t i c at o m in the un i t cell , i m p o s s i b l e to set up an AFM s t a t e . Set as FM i n s t e a d . ’) ;

230 if sum(abs( m a g n e t i c F l a g ) ) == 1

231 m a g T y p e = -2;

232 e ls e

233 m a g T y p e = 2;

234 end

235 end

236

239 m a g T y p e = 1;

240 end

241 end

242

243 i = 1;

Designing a new algorithm for crystal structure prediction in magnetic materials

KTH - Royal Institute of Technology

Master of Science Thesis