Classical and quantum simulations of dislocations in crystals

Classical and Quantum simulations of dislocations in crystals Author: Santiago Sempere? Supervised by: Dr. Claudio Cazorla† Co-supervised by: Prof. Jordi Boronat♣ University of New South Wales†and Universitat Politecnica de Catalunya?,♣ October 17, 2016 Acknowledgements Above all, I would like to express my most sincere gratitude to my supervisor of the thesis, Dr. Claudio Cazorla. From the very first moment, he tried his best to help me out with all the inconveniences. Aside from guiding and counseling me in all the technical aspects, he effectively transmitted me courage and optimism in those tough times when all you have are non-sense results. I also want to thank Prof. Jordi Boronat for giving me the opportunity of living this one-in-a-lifetime experience and giving me some advice on everything I needed. Concerning the technical aspects, I would like to thank the invaluable help offered by Guillem Ferr´e for providing and helping me with the PIMC code to be able to perform the quantum simulations. Also, thanks to prof. Anna Serra for giving us a beam of light on what was going on with our simulations at the time we were more lost than ever. Last but not least, I would like to thank my parents for making possible this experience of studying and living in Australia. Abstract Dislocations are known to play a key role in the plastic behavior of materials. At the quantum level, experimentalists working on the archetypical bosonic quantum solid 4He have observed unusual material properties such as giant plasticity and superfluid mass transport. Although the theoretical explanation for these observations remains elusive, their interpretation has involved the role of dislocations unquestionably. In this thesis, we aim to fulfill the lack of theoretical support for these experiments through atomistic simulations of dislocations. As a first approach, we have characterized the dislocations under classical conditions for a hcp system of Xe with an LJ interatomic potential. Our results reveal the dissociation of the dislocation into two Shockley partial dislocations bounding a broad region of stacking fault. Also, our findings show a very small Peierls Stress τpwhich results in the absence of lattice resistance to dislocation motion at finite temperature. This provides a key insight into a behavior thought to be exclusive to quantum systems. To assess the features of the dislocation at the quantum regime, we employ Path Integral Monte Carlo simulations. We have applied on-the-fly and a posteriori methods of analysis to compute the behavior of the dislocation, but none of them have provided conclusive results. Nevertheless, clear evidence for superfluid-like behavior in either the dislocation cores or stacking fault region is not observed in our simulations. Contents 1 Introduction 4 1.1 A bit of history and motivation . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Objectives and outline of the thesis . . . . . . . . . . . . . . . . . . . . . . 6 2 Hcp crystals and theory of dislocations 8 2.1 Hexagonal Close-Packed crystals . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 CrystallineDefects ............................... 11 2.2.1 Burgers circuit and Burgers vector . . . . . . . . . . . . . . . . . . 13 2.2.2 Types of dislocations . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.3 Partial dislocations. Shockley partial dislocations . . . . . . . . . . 15 2.2.4 Dislocationmotion........................... 16 2.3 Dislocations in hcp structures . . . . . . . . . . . . . . . . . . . . . . . . . 18 3 Simulation methods 22 3.1 Classicalmethods................................ 22 3.1.1 Staticsimulations............................ 23 3.1.2 Dynamic simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Quantummethods ............................... 26 3.2.1 Path Integral Monte Carlo (PIMC) . . . . . . . . . . . . . . . . . . 26 3.3 Analysismethods................................ 33 3.3.1 Burgersvector ............................. 33 3.3.2 Angular distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.3 Differential displacement analysis . . . . . . . . . . . . . . . . . . . 35 3.3.4 Nearest neighbors analysis . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.5 Common Neighbor Analysis (CNA) . . . . . . . . . . . . . . . . . . 37 3.3.6 Neighbor-Common Parameter (NCP or CNP) analysis . . . . . . . 38 4 Results 40 4.1 Classicalregime................................. 40 4.1.1 Zero temperature properties . . . . . . . . . . . . . . . . . . . . . . 41 4.1.2 Finite temperature properties . . . . . . . . . . . . . . . . . . . . . 52 4.1.3 About the boundary conditions . . . . . . . . . . . . . . . . . . . . 56 4.2 Quantumregime ................................ 57 5 Conclusions 62 5.1 Classicalregime................................. 62 5.1.1 Relaxation................................ 62 Author: Santiago Sempere Simulation of dislocations in crystals 5.1.2 Peierlsstress .............................. 63 5.1.3 Finite temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Quantumregime ................................ 64 6 Further work 65 3 Chapter 1 Introduction 1.1 A bit of history and motivation In the 1970’s Alexander Andreev and Ilya Lifshitz at Moscow [1] and Geoffrey Chester at Cornell University [2] suggested the existence of a state of matter in which crystalline order and Bose-Einstein condensation coexist, from their investigations in solid 4He. They called it supersolidity and it was based on the presence of a measurable number of vacancies which at very low temperatures permit a part of a solid to flow on top of the other. This hypothesis aroused big expectations within the experimentalists and some of them started trying to reach this state of matter employing several mass flow and torsional oscillator experiments [3] at temperatures of few dK. The whole idea of the experiment is that if some of the atoms could flow on top of the others, i.e. decouple, this would result in a shift of the oscillation period. Firsts experiments found no evidence of this ability to flow, revealing the non-existence of the so-called supersolidity and killing the enthusiasm of most experimentalists on the topic for over two decades. It was not until Moses Chan and Eusong Kim [4], re-reproducing the experiments previously proposed by Legget [5], the pioneer in the torsional oscillator, but introducing a higher number of vacancies, observed that there was a change in the period when cooling the solid below 0.15 K. These results were interpreted as a clear evidence of what was likely to be the discovery of a new state of matter, however, remained unclear. From the theory, we know vacancies have a too large formation energy (∼15K) so as to be present in a relevant number in the ground state of solid 4He. Meanwhile, many researchers followed the experiments conducted by Chan and Kim [4]. In 2007, James Day and John Beamish [6] studied the elastic and mechanical properties of solid 4He, concluding they were originated by the presence of dislocations and their interplay with impurities of 3He. They found an intriguing coincidence between the temperature at which the shift of period occurred, i.e. 0.15K, with the temperature at which the shear modulus increased sharply. Then, the following question arose: was the change in the period a matter of a change in the structural properties instead of the existence of a new state of matter? This idea gained a major weight afterward, when other experimentalists [7,8] showed a correlation between the moment of inertia and the structure of the solid. Finally, in 2012, a better-designed torsional oscillator was built by Chan and Kim [9] and no more evidence of supersolidity, within the experimental Author: Santiago Sempere Simulation of dislocations in crystals errors, was found. These new observations made clear that the behavior observed in the first experiments was due to the same causes underlying the mechanical and elastic properties of 4He observed by Day and Beamish (i.e. dislocations and 3He impurities). As a result of the frustrated searches of a new state of matter and following these works, Haziot et al. began to study how the stiffness of the quantum solid 4He was affected by the temperature [10], deriving the conclusions illustrated in Fig. 1.1. They were able to show how the stiffness, i.e. the shear modulus, drops when cooling down the temperature until it reaches a point, T≃0.15 K, where it drastically grows again. This behavior was named as ”Giant plasticity”. Experimentalists attributed this behavior to the free gliding of dislocations in a particular direction. The fact that below a critical temperature the stiffness grows, what could be counterintuitive, is thought to be due to the pinning of the dislocations by the 3He impurities, which at such low temperatures are more likely to condensate and mix with the pure solid 4He. However, this hypothesis lacks either theoretical or experimental evidence to support it. Figure 1.1: Shear modulus of a solid crystal of 4He. At a temperature of T≈0.2K the shear modulus is sharply reduced to 72 bar, much less than the normal value of 127 bar. Figure taken from [10] Recently, Boninsegni et al. [11], via atomistic simulations, predicted the existence of superfluidity in the core of a screw dislocation in 4He, in fact, a 1D Luttinger-liquid system. This prediction led to many following hypotheses involving superfluidity phenomena along dislocation cores, resulting in the prediction of other quantum mechanical properties such as dislocation superclimbing or the syringe effect [12,13]. This work preceded other experiments [14–18] carried out by Ray, Hallock et al. showing the existence of ”superfluid mass transport”. The intepretations of these experiments infer the existence of superfluid dislocation cores, as predicted by Boninsegni et al.. However, the very recent work published by Borda et al. [19] came to show that the cores were not superfluid, in either the edge or the screw dislocations, invalidating Boninsegni’s prediction. Far from a clear conclusion, further tests on this topic are required. These experiments and theoretical predictions show how very unexpected phenomena can be observed at the quantum crystals and explains why there exists a growing interest in the quantum field from materials’ scientists. From this perspective, we conclude that a theoretical explanation must be done, since many yet-to-be-understood results emerging from the experiments exist and a lack of knowledge in the quantum dislocation field has to be fulfilled. Thus, the aim of our work is to be able to give a proper theoretical explanation of this phenomena by means 5 Author: Santiago Sempere Simulation of dislocations in crystals of atomistic simulations at a quantum level. A first approach from the classical point of view is strictly necessary to set the foundations of the work and to have a reference on what to expect in the quantum regime, although, as we know, everything can happen within it. 1.2 Objectives and outline of the thesis The primary purpose of this thesis is to give a theoretical explanation of the phenomena observed in the experiments mentioned previously, through the simulation and analysis of atomistic models of 4He at the quantum regime. In these simulations, we would like to see the existence, or in its defect, the nonexistence, of superfluidity and the appearance of free gliding of dislocations under the application of a minuscule shear stress, what would lead to the so-called giant plasticity. Moreover, an ambitious study of the behavior of the dislocations in the presence of 3He impurities, to show the role these impurities play in the quantum regime would be liked to carry it out if possible. Several experimentalists claim they have a pinning effect on the dislocations, disabling the free gliding. In order to accomplish this specific purpose, we have to make a first approach to dislocations from the classical regime to make the path smoother into the quantum dislocations. Furthermore, these first classical simulations could provide us with valuable data to compare with the quantum one. In this sense, we study the properties, both dynamic and static, of an edge dislocation at the classical regime, i.e. using classical potential interactions and at a classical finite temperature. The outline of the thesis is the following: 1. In chapter 2, we describe from a theoretical point of view the basic players of our game, the hcp crystal and the dislocations. We deeply explain the properties of each of them, doing a brief review in crystallography and showing step by step the features that characterize the dislocations. 2. In chapter 3 we present the different techniques used for the study of the dislocations at the classical and quantum regimes. We mainly introduce the technique of Molecular Dynamics and give rather a simple explanation of how the Feynman’s formalism of path integral is implemented into a Monte Carlo method to perform a quantum simulation. Also, we expose the tools used -or triedto characterize and monitor the dislocation during or after the simulation. 3. Chapter 4 is devoted to showing the results obtained in each of the stages of the project. Consequently, we can distinguish two main parts: the classical regime and the quantum regime. Also, some interconnections and interrelations are done among both methods analyzing differences and similarities between them. By doing this, we can have an idea if whether we can rely on the results obtained at the quantum regime, since not much literature is available, or not. This is the most relevant part of the thesis since we reveal the results of the simulations performed under different conditions. Moreover, we study the relevance of size errors in our system, which can be found to be considerably high and incur wrong 6 Author: Santiago Sempere Simulation of dislocations in crystals interpretations. This fact enhances the importance of a good classical background, in which we can reach significant systems in our simulations, in order to estimate the effect of the size and bear it in mind for the quantum simulations. In the classical regime, we have used atoms of Xe and a LJ potential, since it is not possible to simulate solid 4He at such high temperatures without applying very high pressure to give cohesion to the system. 4. In chapter 5, we expose the achieved conclusions from our work and discuss the validity of our experiments due to countless simulation problems, such as system sizes or erroneous interpretations. Also, we will validate some of the hypothesis done by the experimentalists and mentioned in the first chapter in section 1.1. 5. Finally, in chapter 6, we propose a path to follow to improve our work and develop further investigations on this intriguing field, such as the simulation of a crystal of solid 4He containing dislocations and 3He impurities. 7 Author: Santiago Sempere Simulation of dislocations in crystals possible jumps v: up v1= [01], down v2= [01], right v3= [10] and left v4= [10]. Fig. 2.5 shows how to perform this analysis in a 2D system. The technique can be easily extended to a more complex three-dimensional system. Furthermore, accordingly to what has been said, the resulting Burgers vector bis: b=rN−r0= N X i=1 ∆ui(2.7) where, ∆ui:= ri−(ri−1+vi) (2.8) is the difference between where the atom is in the actual, distorted, system and where it should be with respect to the perfect lattice. This value is computed at every step and summed up at the end. This way of computing the Burgers vector will be especially useful when handling with dissociated partial dislocations. Next on, we will explain how we dealt with the issues related to the partial dislocations. Section 3.3.1 is devoted to the implementation of an algorithm to find the Burgers vector and its position. Figure 2.5: Three Burgers circuit depicted in an atomic plane perpendicular to the edge dislocation line in a simple cubic lattice. Notice how circuit 1 and circuit 2 are non-closed since they enclose the dislocation, whereas circuit 3 is a closed loop. Ei and Siare the corresponding starting and finishing atoms of the corresponding Burgers circuit. The sense of the vector ξis defined to point out from the paper so that all the circuits flow in the counterclockwise direction. Figure taken from [21] 14 Author: Santiago Sempere Simulation of dislocations in crystals 2.2.2 Types of dislocations There exist two type of dislocations: screw and edge dislocations. To better understand the geometry and creation of a dislocation let’s introduce them through a SC lattice. In Fig. 2.6 we can visualize the concept of both dislocations. Fig. 2.6a is an edge dislocation that has been created by the insertion of a half of a plane in the solid. On the other hand, Fig. 2.6b shows an screw dislocation that has been created by ’cutand-slip’ operation. The plane defined by the Burgers vector band the dislocation line ξ(section 2.2.1) is called the glide plane. For the edge dislocation, where the Burgers vector is perpendicular to the dislocation line, the glide plane is well-defined, whereas there is not a well-defined glide plane for the screw dislocation because its Burgers vector and dislocation line are parallel to each other. (a) Edge Dislocation (b) Screw Dislocation Figure 2.6: Types of dislocations. Figure taken from [21] 2.2.3 Partial dislocations. Shockley partial dislocations As introduced before, a partial dislocation is a dislocation whose Burgers vector bhas a modulus smaller than a lattice vector, for instance, bp= [1 20]. Then, to compute its value, a small correction is made to Eq. (2.6) as follows: b= N X i=1 ∆uiH(−|∆ui|) (2.9) with H(x) being the Heavyside function (Step function). is a parameter that chooses when the distortion is too big to be accepted in the sum or not, it is calibrated by trial-and-error. Computing the difference between the end and the beginning positions would lead us to an incorrect result, either a whole lattice Burgers vector or a null one. Partial dislocations are found in the presence of a stacking fault. In fact, each of the partial dislocations places at the end of the stacking fault, where the stacking fault merges with the perfect lattice structure. The formation of a low energetic stacking fault in an hcp with a dislocation is very common; therefore, the dislocation tends to split up into two Shockley partial dislocations. 15 Author: Santiago Sempere Simulation of dislocations in crystals AShockley partial dislocation is a kind of partial dislocation. These are the ones associated with a slip, and its formation can be compared to the one of an edge dislocation in an elastic model. The creation of these partial dislocations is due to a low energetic stacking fault that is more energetically favorable than a single edge dislocation. If we depict the stacking fault energy surface (γ), we would appreciate that the energy of the system after a certain displacement fshows a minimum, corresponding to a metastable state. In Fig. 4.2 in section 4.1.1 is represented the stacking fault energy surface (γ) for α-Zr and we can observe the minima that explains the dissociation into two Shockley partial dislocations. 2.2.4 Dislocation motion We can distinguish between two main types of dislocation movement. •Glide or conservative motion: the dislocation moves in the plane formed by the dislocation line and the Burgers vector. This kind of movement is the most typical one and it is even more predominant at low temperatures. This kind of mobility is very anisotropic with respect to the glide plane in the non-screw dislocations. •Climb or nonconservative motion: the dislocation moves out of the glide surface, and, therefore, normal to the Burgers vector. Climbing only occurs at high temperatures when the insertion or emission of atoms is possible. It is called non-conservative motion due to the change in the number of atoms contained in the crystal. In the screw dislocations the glide plane is not defined; therefore, the only possible movement is gliding. As mentioned before, the motion of a dislocation is highly related to the mechanical properties of the crystal, in other words, the motion of a dislocation is the mechanism for a crystal to deform plastically. In most cases, to do so we are required by an external force, the driving force, that pushes the dislocation until it moves. This force is applied as a shear stress τ=F/A. The dislocation mobility, which can be expressed as a function of the applied forces on the dislocation M(f), is mainly influenced by two types of forces: extrinsic and intrinsic forces. The extrinsic forces are those provoked by obstacles or impurities found in the crystal or externally. On the other hand, the intrinsic forces are those that arise from the interatomic interactions and due to the lattice resistance. Next section explains the latter type of forces in more detail. 2.2.4.1 Intrinsic resistance to dislocation motion Within a crystal there exists an energy barrier, when no external forces or energy are added, that oppose the dislocation movement, due to the periodicity of the lattice. In this case, it is called the Peierls energy barrier, named after the theoretical physicist Rudolph Peierls, and is the energy needed to break and create the bonds between the atoms in the dislocation core. Hence, this energy depends sensitively on the form of 16 Author: Santiago Sempere Simulation of dislocations in crystals the force-distance relation between individual atoms, i.e. on the interatomic potential. Fig. 2.7 illustrates the process of breaking and creation of bonds. Figure 2.7: Dislocation gliding. A shear stress has been applied and the dislocation moves from left to right, resulting on a plastic deformation in the end. Figure taken from [21]. Another vital concept in dislocation motion is the Peierls stress. This is the minimum shear stress required to move the dislocation one lattice vector at T= 0. The calculation of this value is essential in all the studies about dislocations since it gives a measure of the lattice resistance of the crystal, which is directly linked to its plasticity. The value of the Peierls stress is related to the disregistry of the atoms across the slip plane. In Fig. 2.8 we can appreciate how the atoms are displaced to accommodate the dislocation in the crystal. The disregistry is characterized by the displacement difference ∆u=u(B)−u(A) between two atoms on adjacent sites above(A) and below(B) the slip plane. The width wof the dislocation is defined as the region where the disregistry is greater than a half of its maximum. In the 1940’s Peierls and Nabarro calculated the dislocation energy per unit of length and found it to oscillate with period of b/2 and a maximum fluctuation, the Peierls energy, given by [22] Ep=Gb2 π(1 −ν)e−2πw b(2.10) where Gis the shear modulus and νis the Poisson’s ratio. The maximum slope of the energy function is the critical shear stress to move the dislocation through the crystal. Dividing by bwe get the Peierls stress τp=2π bEp=G (1 −ν)e−2πw b(2.11) This simple model agrees much better with the experimental results than Eq. 2.6 , that describes the theoretical shear strength in a crystal. From the expression, notice that the τpof a dislocation scales as the negative exponential of the width of the dislocation. Hence, a very wide dislocation core will be very mobile. In general, edge dislocations have smaller Peierls stresses than screw dislocations. 17 Author: Santiago Sempere Simulation of dislocations in crystals Figure 2.8: Scheme of the disregistry suffered by the atoms when we introduce the dislocation. The colored atoms indicate the positions after the insertion of the dislocation, whereas the empty circles indicate the positions in the perfect, undistorted configuration. Figure obtained from [22]. 2.3 Dislocations in hcp structures As shortly introduced previously, dislocations in crystals with an hcp lattice structure undergo specific processes that make the dislocation behavior more complex than in other lattices. Making use of the Thompson tetrahedron, Fig. 2.9 for hcp structures (it is normally used for the fcc lattice), we can describe the most important and common dislocations found in hcp structures. •Perfect dislocation with the Burgers vector being one of the vectors in the basal plane. Regarding the thetraedron this would be either AB, BC, CA, BA, CB or AC. •Perfect dislocation with the Burgers vector perpendicular to the basal plane, represented by ST and TS. In this case the glide plane is non-basal, is the prism plane (1010). •Perfect dislocation with one of twelve Burgers vector of the type SA/TB •Imperfect basal dislocation of the Shockley kind, regard 2.2.3, with Burgers vector Aσ,Bσ,Cσin either one or the other senses. •Imperfect dislocations with the Burgers vector perpendicular to the basal plane, represented by σS,σTand the counter-sense vectors. •Imperfect dislocations which are a combination of the two latter cases, namely, AS,BS, etc. Although these vectors go from one atomic site to another one, are still considered imperfect dislocations since their surroundings at each atomic site are not identical and they are not a lattice vector. The crystallographic representation of these dislocations via the Miller-Bravais indices is found in Table 2.1 Many of the imperfect dislocations are a result of the existence of a low energetic stacking fault. We can distinguish three primary basalplane stacking faults that do not affect the nearest neighbor arrangements of the perfect stacking ABABAB... , [23] 18 Author: Santiago Sempere Simulation of dislocations in crystals Two of them are intrinsic and called I1and I2, and the third one is extrinsic and called E. The change in the stacking is as follows: ABABABAB...→ABABBABA...→ABABCBCBCB...(I1) (2.12) where the middle stage is produced after the removal of a basal layer, what produces a high energy stacking fault. The system goes to the next stage by a slip of 1 3h1010i, arriving at a low energy stacking fault. ABABABA...→ABABCACA...(I2) (2.13) Fault I2is a result of a slip of 1 3h1010iin a perfect crystal ABABABA...→ABABCABAB...(E) (2.14) The extrinsic fault (E) is a consequence of the insertion of an extra plane. These faults introduce a region in the space where the stacking is different than in the rest of the crystal, thus, in this region a face-centered cubic stacking (ABC) can be found and so have a characteristic stacking-fault energy γ. The way to explore the stacking-fault energy surface γ(f) is to displace a half of a crystal a vector fwhile the other half remains in the same position and letting the system go to the equilibrium position. Repeating this procedure all over the space and mapping for x∈[0, a) and z∈[0,1.6a) we can draw the energy surface. From the data of the energy surface, we can predict and understand why the system goes to one state or the other since it always will try to go to the minima. In chapter 4, where the results are displayed, we will take a look to the stacking energy surface of α-Zr [24] and relate it to the stacking fault and the dissociation of the dislocation observed. Figure 2.9: Burgers vector in the HCP structure. Taken from [25] Dissociation of a perfect dislocation into two Shockley partial dislocations 19 Author: Santiago Sempere Simulation of dislocations in crystals Type AB TS SA/TB Aσ σS AS b1 3h1120i[0001] 1 3h1123i1 3h1100i1 2[0001] 1 6h2203i b a c (c2+a2)1/2a/√3c/2 (a2 3+c2 4)1/2 b2a2c2=8 3a211 3a21 3a22 3a2a2 Table 2.1: Dislocations in Hexagonal Close-Packed lattices The shortest lattice vector of the HCP is 1 3h1120iand the most common slip planes are (0001) and {1100}, which correspond to the basal and the first order prism planes, respectively. The preference of the glide plane is determined by the energy and stability of the stacking fault. If the stacking fault I2with vector 1 3h1120iexists, then the perfect dislocation AB dissociates into two Shockley partial dislocations bounding a ribbon of stacking fault, which has a fcc-like stacking. The reaction of this process is as follows: AB→Aσ+σB(2.15) which in crystallographic notation is: 1 3[1120] →1 3[1010] + 1 3[0110] (2.16) The geometry of the dissociation is two partial dislocations lying on the basal plane at ±30◦to the perfect vector (in some texts it refers to the dislocation line, then it would be ±60◦) and the reduction in energy given by b2is 1/3, as shown in table 2.1. Notice how after the dissociation the partial dislocations are not a pure edge dislocation anymore, since they have a screw component, although the average screw component for the complete dislocation remains to be null. The method of the Differential displacement introduced afterward in section 3.3 is especially useful to differentiate between the edge and the screw components of the Burgers vector in a dislocation. The schematic representation of this dissociation is illustrated in 2.10 Figure 2.10: Schematic representation of the I2stacking fault bounded by two Shockley partial dislocations (Aσand σB). The arrows indicate the two errors in the stacking. Notice how it varies to ABC, corresponding to a face-centered cubic structure. Figure taken from [22]. 20 Author: Santiago Sempere Simulation of dislocations in crystals Despite the existence of other low-energetic stacking faults that led to other dislocation dissociations, we are not going to go into the detail of these, since the principal and most commonly observed one, actually, the only observed in our work, is the previous one. 21 Chapter 3 Simulation methods This chapter is devoted to the description of the simulation methods used during the project to understand and simulate the behavior of the dislocations in a crystal under multiple conditions. Firstly, we will introduce the classical methods, in particular, the Molecular Dynamics. It will allow us to analyze the energy and the structural and dynamical properties of the dislocations in classical crystals. We will use this piece of information as a reference set to which compare the results obtained in quantum crystals. Secondly, we introduce the Quantum methods and, especially, the one employed in our simulations, the Path Integral Monte Carlo. Why do we use atomistic simulations? The interatomic interactions are the fundamental basis underlying the mechanical properties of materials. Therefore, to understand the behavior of dislocations, and, hence, the plastic deformation of a material, it is necessary and sufficient to study the collective behavior of the atoms in a crystal containing a dislocation. Indeed, other researchers have used continuum models, as Pessoa et al. [26], to perform a simulation of the mobility of a dislocation in a quantum crystal. In these cases, the assumption of many classical values for the parameters ruling the model makes the simulation less rigorous and reliable from the physical point of view. 3.1 Classical methods Quantum mechanical motion and the interaction between the electrons can be important and relevant in the interatomic forces. This fact makes, sometimes, very complicated to describe the interaction among atoms, since a solution of the Schroedinger’s equation to describe the electronic interaction is needed to be strictly rigorous. The methods that use this solution are known as first principles methods or ab initio methods. The intrinsic difficulty of the Schroedinger’s equation makes these methods very costly from the computational point of view, until the point it is unfeasible to simulate a system of more than a few thousands of atoms. Commonly, to describe well the behavior of a dislocation we need much bigger systems that can only be approached by a less sophisticated model. In the following chapters, we explain how this model is. Author: Santiago Sempere Simulation of dislocations in crystals 3.1.1 Static simulations 3.1.1.1 Relaxation The probability to find a system in an state µcharacterized by {ri,pi}in the phase space and in thermal equilibrium at a temperature T is defined by the Boltzmann’s law, as follows: p(µ) = 1 Ze−H({ri,pi}) kbT(3.1) where H({ri,pi}) = N X i=1 |pi|2 2m+V({ri}) (3.2) is the Hamiltonian of the system, Zis the partition function, which ensures the proper normalization of the probability density, defined as: Z=ZN Y i=1 dridpie−H({ri,pi}) kbT(3.3) and kbis the Boltzmann’s constant. According to this probability distribution, the odds to find the system in a determined state, for instance µ, decreases exponentially with increasing energy H({ri,pi}). Eventually, at the low-temperature limit, the most likely state of the system is at the global minimum of the energy surface H({ri,pi}). The global minimum of the energy, from Eq.3.2, is found when pi= 0,∀i. Furthermore, the minimum of the potential energy V({ri}) gives a good description of the system at low temperatures. Then, a relaxation of a system consists on searching for the minimum of the potential energy V({ri}) in order to find the most stable state of the system at T= 0K. Searching for minima is a widely studied field, and still an active area of research in computational sciences. There exist several algorithms to seek for the minima of a function, such as the steepest descent algorithm, the Hessian-free truncated Newton algorithm and the conjugate gradient algorithm. In the present work, we have used the Polak-Ribiere version of the more general conjugate gradient relaxation (CGR) algorithm [27]. The CGR algorithm relies on the computation of the atomic forces to displace the atoms in such directions. The algorithm works iteratively until the condition |F|<  is reached. It is worth mentioning the CGR algorithm does not guarantee to arrive at a global minimum of the energy, only a local one. Fig. 3.1 illustrates why this can happen. V(x) has three relative minima, two local ones, and the global one. If we use the CGR method for the relaxation, we may end up at a local minimum, since the forces (F=−∂V (x)/∂x) will be zero. To ensure we arrive at the global minimum, we can use brute-force and run the simulation starting from many randomly selected initial configurations. Clearly, though, it is a very inefficient way to proceed. 23 Author: Santiago Sempere Simulation of dislocations in crystals and, hence, the expectation value is: hˆ Oi=ZdRρ(R,R;β)O(R)≈ZM Y j=1 dRjρPA(Rj+1,Rj;)O(Rj) (3.26) applying the boundary condition RM+1 =RM. The main feature we must notice in this equation is the fact that the product p(R1, ..., RM+1) = QM j=1 ρPA(Rj+1,Rj;) is positive definite and its integral over the whole space is equal to 1. Therefore, can be interpreted as a probability function and is suitable to work as the probability distribution that rules the sampling of the degrees of freedom in, for instance, the Metropolis algorithm. Thus, many observables that give important physical properties can be computed at any temperature in a rather simple way [31]. Moreover, in virtue of the Trotter formula, as we get closer to the limit M→ ∞ the integral becomes exact and we obtain the exact value of the expectation value ˆ O. This is the reason why this method is often referred as an exact method. The Classical Isomorphism The Path Integral Monte Carlo (PIMC) method describes a system of N particles by means of M different configurations RMthat constitute what we have called the ”path” in the imaginary time in the space of configurations. However, this can be interpreted as a classical system made of N×Mparticles, each of these called bead. Somehow, each of the Ninitial particles becomes a polymer formed by Mbeads. From the identity M X m=1 (Rm−1−Rm)2= M X m=1 N X i=1 (ri,m−1−ri,m)2= N X i=1 M X m=1 (ri,m−1−ri,m)2(3.27) we can see the kinetic interaction (Eq. 3.19) of each of the beads and interpret it as if each of the beads is connected to the next one by a spring-like bond. The condition rM+1 =rMintroduced in Eq. 3.26 guarantees that the first bead is connected to the last bead, forming a closed polymer. Now, adding the potential interaction, i.e the interaction between beads of different polymers, does not change this picture. In the primitive approximation, the potential contribution is given by Eq. 3.20 and the interaction between classical beads is identical to the two-body potential of quantum atoms. Nevertheless, these classical polymers interact in a very particular way, since only the beads at the same ”time”, i.e. imaginary time, in the path interact with each other. This makes the computation so much simpler. Regard Fig.3.3 for a visual understanding of the presented case. In summary, the PIMC method consists of considering the quantum particles as classical ring polymers interacting between them in a particular way. Each polymer consists of a chain of beads connected by ideal springs. This representation enables us to introduce the delocalisation of the particles due to its zero-point motion. The action is defined as minus the logarithm of the density matrix. For a given j, this is: S(Rj+1,Rj;) = −ln[ρ(Rj+1,Rj;)] (3.28) 30 Author: Santiago Sempere Simulation of dislocations in crystals Figure 3.3: Scheme of two classical ring polymers representing the quantum atoms in the PIMC formalism. Each of the numbered circles is a bead and its index is an instant in the imaginary time. The springs connecting each of the beads in each chain shows the kinetic action and the red dashed line represents the potential action acting on beads in the same imaginary time, i.e. same index so specifies the interaction between the beads in the classical analogy of the quantum system. It strongly depends on the choice of the thermal density, see the review [35] for how to get a ”good” action. Particularly, the kinetic action derived from Eq. 3.19 is Sk=MkBT 4λ((Rj+1 −Rj)2(3.29) Such expression gives an intuitive idea of how the spring-like interaction behaves in the polymers. At high temperatures, the quadratic term is big and, hence, the harmonic potential between beads is strong. Hence, the beads are close to each other and does not allow the polymer to spread in the space, reducing the delocalization of the particle. As the temperature decreases the polymer becomes bigger and the particle is more delocalized in the space, increasing its zero-point motion. The number of chosen beads Mto represent each polymer is essential to get a reliable and well-behaved system and it strongly depends on the temperature. For high T, where no delocalization is found at all, we are able to represent the quantum system by a small number of beads per polymer, being M=1 for the classical limit where the polymer shrinks into a single particle. In the other hand, for low temperatures, we need to increase the number of beads M to make a proper representation of the system. Summarizing, the number of beads scales with the inverse of the temperature, i.e. M∝1/T. This fact represents one of the weaknesses of this method, since when we approach the most interesting part from the quantum level, this is, at very low temperatures, the computation becomes more and more costly due to the very low efficiency of the sampling of the long chains involved. The way to improve its efficiency is to make a better approximation for the thermal density matrix, in order to work with greater values of . 31 Author: Santiago Sempere Simulation of dislocations in crystals The Takahashi-Imada approximation [36,37] and the Chin approximation [38] are among the approximations that make the algorithm more efficient and capable of performing simulations of systems at lower temperatures. We will not get into the detail of these approximations in this work, see each of the mentioned references for a very detailed explanation or [31] for a general overview of these approximations. The permutation sampling The previous sections highlight the most important features of the PIMC algorithm, but it is not entirely complete. Indeed, the previous algorithm only holds for systems made up of distinguishable particles. When attempting to deal with quantum many-body systems, we have to take into account the quantum statistics of particles. To recover the right expression for the thermal density matrix ρ, when dealing with the Bose or Fermi statistics, we have to sum over all the possible permutations of the particle labels in one of the two arguments, deriving the following expression: ρB F=1 N!X P (±1)Pρ(R1,PR2;β) (3.30) where Pis one of the N! permutations of the particle labels, P is the number of transpositions of the permutation Pand the term (±1) is either + or −if we are dealing with bosons or fermions, respectively. Analogously to the distinguishable particles case, we can recognize a probability distribution and perform a Monte Carlo procedure to calculate the sum over the permutations. A remarkable detail that changes from the previous sections is the symmetrization of the thermal matrix density, where the boundary condition RM+1 =R1becomes RM+1 =PR1.PR={rp(1),rp(2), ..., rp(N)}, p(i) being the label of the particle permutated with the i-th atom. This new condition indicates that the last bead of the i-th polymer is no longer compulsory connected to the first bead of the same chain, but it is connected to the first bead of the p(i)-th polymer. This makes a substantial difference in the mapping onto the classical system, since now exists the possibility that the system is not composed of Nidentical closed polymers made up of Mbeads, but it is possible to find out polymers formed by L×M beads, which represents the permutation cycle between Lbosons. To sum up, the explained so far, taking into account the indistinguishability of the particles can lead some of the polymers to open their chains and join other opened polymers to form a longer chain. Concerning the physics of the problem, the formation of longer polymers increases the delocalization of the particles, decreasing the kinetic energy of the particle. If we decrease the temperature enough, the polymers are very spread out and the beads of different polymers are close together making it very plausible that they ”collapse” and form a longer chain. In this latter case, we can arrive at a limit for which a macroscopic size polymer is formed and, consequently, a phase transition occurs and we have a fraction of the atoms in a superfluid phase. In a periodic system, the order parameter is measured by the winding number: the number of times a polymer wraps around the periodic boundary conditions. From a practical point of view, the implementation of these movements are difficult and, thus, makes the computation of the superfluidity unfeasible for systems greater 32 Author: Santiago Sempere Simulation of dislocations in crystals than N∼100 atoms. However, the worm algorithm, first implemented by Prokof’ev et al. and extended to the continuous case by Boninsegni et al. (see Ref. [39–41] for further details), provides a tool that allows us to give an efficient description of the thermodynamic properties connected to the bosonic statistics of the quantum systems, such as superfluidity. 3.3 Analysis methods In this section, we describe the techniques used to estimate the structural and dynamical properties of dislocations in the simulations. Different estimators were tested to determine which one describes more accurately the location and features of the dislocation during the simulation. We performed extensive classical molecular dynamics tests since the computational expense associated with this type of simulations is low in comparison to that of quantum simulations. In view of our MD test outcomes, we decided which estimator was more suitable to be implemented in the PIMC code, both in terms of accuracy and computational affordability. 3.3.1 Burgers vector As explained in section 2.2.1, the Burgers vector is the property that fully characterizes the dislocation, therefore, it is logic to begin with it in the search of a suitable estimator that allows us to monitor the dislocation. The implementation of the algorithm to compute the Burgers vector and its position is the following: 1. The first step contains three simultaneous steps (a) Set the dimension of the Burgers circuit. It is important to bear in mind that there is a trade-off in the choice of the dimension of the Burgers circuit since the choice of a large circuit will correspond in a much less costly computation, but the precision in the position of the Burgers vector will be very low. Notice that this method only gives us an intuitive idea of where is the dislocation and it is strongly conditioned by the dimension of the Burger’s circuit. In our case, we decided to make a circuit of 2 ×2. (b) Choose a determined plane perpendicular to the dislocation line. The Burger’s circuit will run over this plane. (c) Choose the initial atom. It is important to recall that we have to explore the whole plane, these means we will have to sample the plane and do many Burgers circuit on it. Fig. 3.4 illustrates the exposed idea. 2. From the initial atom r0move to the nearest atom (ri) to the desired position, ri+vi. 33 Author: Santiago Sempere Simulation of dislocations in crystals Figure 3.4: Schematic representation of the implementation of the algorithm to compute the position and value of the Burgers vector. We have sampled the chosen plane in 6 different circuits each of them starting from the atom with the big dot in it. The circuit is 2 ×2 and we observe that the circuit number 5 is not closed, since it encloses the dislocation and we can find the burgers vector. 3. Compute the difference between the actual position of the atom and the place it should have been, as explained in section 2.2.1, ∆ui=ri−(ri+vi) 4. Repeat the same procedure for all the circuit, this is for i= 2,3,4. 5. Compute the Burgers vector from the relative displacements of each step ∆ui, as b= N X i=1 ∆uiH(−|∆ui|) 6. Repeat the same procedure spanning all over the plane to find the dislocation core This method of analysis relies on the fundamentals of the dislocation theory, which makes the results very descriptive and complete. The method lacks robustness since it is very sensitive to small fluctuations of the atoms that can lead the algorithm to construct wrong Burgers circuits. Also, the position of the dislocation is not very accurate, since we have a certain range in which the dislocation core must be, delimited by the size of the circuit. Concerning the computational cost, the algorithm computes the nearest neighbors of all the atoms, which is quite costly. 3.3.2 Angular distribution The introduction of a dislocation in a crystal changes its periodicity and, therefore, there is a structural variation from the perfect crystal. One of the properties that changes is the angular distribution of each atom in comparison with its neighbors. In this context, we have created a parameter per atom that measures the structural variation using the computation of the difference between the angle formed by the atom and its four closest neighbors, considering just one plane perpendicular to the dislocation line, and the angle they would form in the perfect crystal. The parameter is calculated as 34 Author: Santiago Sempere Simulation of dislocations in crystals follows: φi= 4 X j=1 (θj i−θ0,j i)2(3.31) where θj iis each of the four angles formed by the atom and its four neighbors and θ0,j i are the angles formed by the atoms in the perfect lattice, i.e. 0o, 90o, 180o, 270o. The strength of this method lies in the fact of its simple implementation and interpretation; however, it is very sensitive to the atoms fluctuating from their equilibrium positions due to thermal energy. Moreover, it does not provide much information about the dislocation or the stacking fault such as its width or the orientation of the Burgers vector. Once again, the most costly part of the analysis is the computation of the nearest neighbors. In Fig. 3.5 we can appreciate the result given from this analysis. As observed, the biggest dispersion in the angle is shown in the dislocation cores, allowing us to identify its position. It is a very visual way of localizing the dislocations, however, unfortunately, it does not give us a clear idea about the width of the stacking fault or the Burgers vector. Figure 3.5: Example of the calculation of the angular distribution in a arbitrary yplane of a system containing 2240 atoms (Lx= 12a). 3.3.3 Differential displacement analysis As introduced previously in section 2.2.4, in order to accommodate the dislocation, the atoms across the slip plane are displaced. We can compute the spread of the disregistry associated with a planar edge dislocation by means of the distribution of 35 Author: Santiago Sempere Simulation of dislocations in crystals the components [bx, by] of the Burgers vector in the glide plane. These components are given by the following expressions: ρbi=d(∆ui) dx (3.32) where i=x, y and ∆uiis described by: ∆ui=uabove,i −ubelow,i (3.33) and uabove,i and ubelow,i are defined: uabove,i =rdisloc above,i −rperf above,i (3.34) where rdisloc above,i are the positions of the atoms above the glide plane in the system with the dislocation and rperf above,i are the positions of the atoms in the non-distorted lattice. The definition of the displacement for the atoms below the glide plane is analogous. Therefore, ∆uiis the difference between the displacement of the atoms above the glide plane and the atoms below the glide plane with respect to the perfect lattice positions. The xcomponent of the Burgers vector corresponds to the edge component of the dislocation, whereas the ycomponent refers to the screw component. The total value of the Burger vector can be computed by integrating over all the values of xalong the glide plane. This is: bi=ZLi 0 ρbidx (3.35) The computation of the differential displacement gives us much information about the dislocation since from its analysis we can derive, for instance, the presence of a stacking fault ribbon bounded by two partial dislocations, what is the width of the dislocation core and what is the separation between the partial dislocations. In chapter 4, we show many examples of this method of analysis, which has shown up to be extremely useful. Concerning the computational cost of the method, it does not require the information about the atom-atom distances, which in most cases is the most costly part of the algorithms. For small systems and at finite temperature, it is very imprecise. 3.3.4 Nearest neighbors analysis An analysis of the nearest neighbors can be very useful to locate the core of the dislocation and, thus, monitor the position of the dislocation during the evolution of the simulation. Basically, we count how many atoms are within a distance rcutoff of each atom. The rcutoff is prescribed to be a value between the distance of the first and second neighbors. A possible cutoff for an hcp system is [42]: rhcp cutoff =1 2(1 + r4+2x2 3)a(3.36) where x= (c/a)/1.633 and ais the lattice parameter. 36 Author: Santiago Sempere Simulation of dislocations in crystals Therefore, if some of the atoms have a value of nearest neighbors different than 12, which is the coordination number for an hcp lattice, this would mean they are in a region of the crystal where the distortion is high and they can be identified as a part of the dislocation core. The algorithm calculates the distances between atoms, which is computationally costly. The major strength of this method lies on its simplicity, both of implementation and interpretation. In the other hand, it gives limited information about the features of the dislocation such as, dislocation core width, stacking fault width (although, it is given implicitly) and Burgers vector components. 3.3.5 Common Neighbor Analysis (CNA) The CNA method was firstly introduced by Honeycutt and Andersen [43] to study the local structure environment. The method consists of creating a diagram for a pair of atoms, αand β, formed by a set of four indexes (1,2,3 and 4). The criterion of the indexes is the following 1. indicates if αand βare nearest neighbors or not. It is 1 if they are and 2 if they are not. Two atoms are nearest neighbors if they are closer than a prescribed rcutoff , which is, in general, the first minimum of the pair distribution function g(r). 2. is the number of common nearest neighbors shared by αand β. In a perfect fcc or hcp lattice is 4. 3. indicates the number of bonds among the common neighbors 4. differentiates between diagrams with same indexes 1, 2 and 3 and different bonding among the neighbors. Table 3.1 shows the distribution of diagrams in each of the perfect lattices. From this analysis, we can determine the local structure in a crystal and distinguish and recognize the different existing defects. CNA diagram fcc bcc hcp 1421 1 0 0.5 1422 0 0 0.5 1441 0 3/7 0 1661 0 4/7 0 Table 3.1: Relative presence of each diagram in the fcc, bcc and hcp crystal structures. Table taken from [44] In our case, this method has shown up to be very useful, so it is a very effective way to localize the dislocations cores and the stacking fault. In section 4.1.1, we show an image where this method of analysis has been used to depict where the stacking fault is since it is a region of fcc structure, and which atoms form the dislocation core. Also, during the simulations at a certain finite temperature, we have used this parameter to monitor the motion of the dislocation. 37 Author: Santiago Sempere Simulation of dislocations in crystals The salient advantage of utilizing this method is that it is already implemented in the LAMMPS code we employ in the classical simulations. On the other hand, its implementation is not as easy as other estimators, what makes it less suitable to use it for the analysis of the quantum simulations. 3.3.6 Neighbor-Common Parameter (NCP or CNP) analysis This method of analysis was introduced by Tsuzuki et al. [44] and it combines the strengths of two well-known methods, the common neighbor analysis (CNA), presented previously, and the centrosymmetry parameter. The three of them are used to characterize the structure of a crystal and to differentiate structural defects such as dislocations, stacking faults, grain boundaries, cracks and surfaces. However, this new method of analysis was chosen as the most appropriate for our case since, on one hand, the CNA has a complex way of implementation and it may be intricate to interpret, and, on the other hand, the CSP is only fully defined for centrosymmetric crystals, i.e. fcc and bcc. Although this does not prevent us from using them, we would incur a loss of accuracy compared to the CNP. Also, another advantage of this method over the CNA is that it does not require an explicit knowledge about what are the opposite neighbors, what for the hcp structure can be a tedious and a not well-behaved numerical task. The CNP consists on the computation of a single parameter Qiper atom. Depending on the value given by this parameter, that atom will be arranged in a determined structur. The definition of Qifollows: Qi=1 ni ni X i=1 | nij X k=1 Rik +Rjk| 2 (3.37) where the index jsums over the ninearest neighbors of atom iand index kgoes over the nij common neighbors between atom iand it’s nearest neighbor j.Rik is the vector pointing from atom ito atom k. Respectively, Rjk is the vector pointing from atom jto atom k. For the bcc and fcc lattices this parameter is null and for the hcp lattice, its value depends on the lattice parameter. In the Path-Integral Monte Carlo algorithm, we must make use of the beads formalism to compute this parameter, thus the definition of the parameter Qiresults to be: Qi,j =1 ni ni X i=1 | nil X k=1 (rk,j −ri,j)+(rk,j −rl,j)| 2 (3.38) where Qi,j is the NCP parameter for the bead jof atom i. The vector (rk,j −ri,j) is analogous to the vector Rik of the classical expression and indicates the vector pointing from bead jof atom kto the bead jof atom i. Now, to obtain the NCP parameter per atom we have to sum over all the beads of each atom, Qi=1 nb nb X j=1 Qi,j (3.39) This would give us the value of the parameter Qithat allows us to identify the local lattice structure within a crystal. 38 Author: Santiago Sempere Simulation of dislocations in crystals The main strengths of this method are its simpleness of implementation and interpretation. Computationally, it is costly since it requires the computation of the atom-atom distances. 39 Author: Santiago Sempere Simulation of dislocations in crystals Figure 4.6: Schematic representation of the system employed to retrieve the Peierls stress in Method 1. We can differentiate three parts: U the upper part, L the lower part and A the region of mobile atoms. The letter P indicates the periodic images. 4.1.1.2 Peierls stress The Peierls stress is a key quantification of the lattice resistance to dislocation motion in a crystal. This value is commonly referred as the critical shear stress (CRSS) for glide at T= 0K. In the present work, we have used two methods to compute this value, which should retrieve similar results. 1st Method This procedure is the most reported one in the literature and it is the most used among researchers to compute the Peierls stress τp[24,47]. This method consists of dividing the solid into three different regions, U -the upper part-, A -the region of mobile atomsand L -the lower partalong the zaxis and displacing U a certain increment while relaxing the region of mobile atoms. Fig.4.6 shows an scheme of the division of the system. In our case, since the dislocation line is parallel to the yaxis, we displace the upper slab of the crystal an increment ∆uin the xaxis. Hence, we can depict the xz −σxz graph, where the strain xz is xz = ∆u/Lz and the shear stress is σxz =Fx/LxLz, where Fxis the sum of the forces exerted on the atoms of the upper part. After each displacement, we minimize the potential energy until the A region of the system reaches the equilibrium. The atoms in the L and the U regions are freeze and fixed to their positions after the displacement. The results obtained are found to be very sensitive to the size of the system and to the thickness of the upper region we pick. Also, it is important to choose a proper strain increment. In Fig. 4.7 we can appreciate the effects of the size of the upper region in the computation of the Peierls stress. For the narrower slab, we obtain much more inaccurate results and they are not consistent since for a slight increment of the strain we obtain substantial differences than for the others. The results correspond to a system of 18424 atoms, contained in a box of Lx = 26band Lz = 12c. The stress profiles obtained from this simulations are not very satisfying. We would have expected to get a first part of the plot showing a linear profile, exhibiting the 46 Author: Santiago Sempere Simulation of dislocations in crystals (a) Evolution of the stress as a function of the strain in a crystal of 18424 atoms and where the upper slab is formed by a region of vertical size of ∼1c (b) Evolution of the stress as a function of the strain in a crystal of 18424 atoms and where the upper slab is formed by a region of vertical size of ∼2.5c Figure 4.7: The lower plots of both figures are a zoom in of the upper plot, in order to be able to observe the elastic behavior of the crystal. elastic behavior of the crystal, until the dislocation is strained enough to move to the next valley of the energy surface, indicated by a sharp drop in the shear stress; and this repeated periodically. A hint of this behavior is shown in our simulations; however, it is not as regular as others have reported, for instance, see [47] or [24]. Although, we have made an estimation of the Peierls stress from Fig. 4.8. We have averaged the stress value since the point of the strain at which the dislocation starts to move, which consistently coincides with the maximum of the stress value obtained, until the end of the simulation. This estimation gives us a reference value of τp= 1.759 ±.17116MPa To prove the accuracy of the method, the Orowan formula gives a theoretical relationship between the plastic shear strain , the mean distance of the dislocation motion x 47 Author: Santiago Sempere Simulation of dislocations in crystals Figure 4.8: Evolution of the stress as a function of the strain in a crystal of 18424 atoms and the dislocation position as a function of the same strain. and the dislocation density ρD[22] =xbρb=xb LxLz (4.2) which substituting the values obtained in our simulations, x= 6b, we get = 0.01001 (4.3) In our simulation we have that at the point the dislocation starts to move ini = 0.0031 and the strain for which the dislocation has moved over 6 Burgers vector is 6b= 0.0143. Therefore, the difference in the strain is ∆=6b−ini = 0.0112, which is in good agreement with the obtained value. Osetsky and Bacon [47] report that accurate results for this analysis are only obtained when simulating a box with Lx≥130band Lz≥80c(These are the lengths of the sides of the plane perpendicular to the dislocation line). Therefore, we have constructed a system of 354455 atoms and the mentioned box sizes to perform the same simulation. The results obtained are shown in Fig. 4.9. Notice the regular pattern obtained for the stress in Fig. 4.9. The Peierls stress obtained in this case is σp= 0.6415 MPa. This value is far from the obtained with the system of 18424 at, meaning that the size effects are very significant, but it allows us to estimate the order of magnitude of the Peierls stress to be around ∼1MPa. In the latter case, the displacement every step is larger, ∆u=b/10. It is important to stand out the boundary conditions used in this attempt to obtain the Peierls stres. In this case, we have used PBC in the glide plane directions, but non PBC in the +z direction, since it will make another dislocation appear above and under the system. This is not a problem for the classical case, since the simulation of a system without PBC is within the range of our possibilities, however, for the quantum computations we strictly need the use of PBC to be able to get reliable results, since we are very limited on the size of the system. Moreover, in this method a part of the system, U and L, does not take part in the relaxation of the system since they are forced 48 Author: Santiago Sempere Simulation of dislocations in crystals Figure 4.9: Evolution of the stress as a function of the strain in a crystal of 354455 atoms and the dislocation position as a function of the same strain. to stand still. Once again, this procedure is unfeasible for the quantum simulations. This is the reason why we need to seek for an alternative way of computing the Peierls stress. A possible way to solve the boundary problem is to introduce a system with a dipole of dislocations on it, or even better, with a quadrupole of dislocations, which allow us to apply boundary conditions and it reduces the effect of the image dislocations. 2nd Method This second method attempts to reproduce the same results achieved by the other method by studying the energy of the system at every step. In other words, in this approach, we simply modify the box by adding a tilt (xz, xy, yz), remapping the atoms according to this tilt and relaxing the system to reach its minimum in the potential energy surface. After the tilt, the box is not anymore orthogonal but triclinic. The lattice vectors that describe it are: A= (a, 0,0), B = (xy, b, 0) and C= (xz, yz, c). This deformation will introduce a shear stress into the system that can be computed by means of the derivative of the energy with respect to the strain. In the present case, we will only introduce a xz tilt. The exact expression is as follows: σxz =1 V ∂E ∂nxz (4.4) where σxz is the shear stress (MPa) applied to the system by the deformation, Eis the potential energy of the system, i.e. the potential energy since we are at T= 0K, and nxz =xz/Lzis the displacement strain defined from the tilt factor xz. In this method we use PBC in all three directions, to make a more realistic approach, despite the creation of other dislocations at the boundaries. For the previous system of 18424 at. and 2240 at. we compute this analysis. In Fig. 4.10 we plot the energy of the system, the stress of the system and the position of the dislocation as a function of the strain for the system of 18424 atoms. For the smaller system, we obtain roughly a similar profile. The energy profile is consistent with the expectations since it shows a certain periodicity distorted by the presence of 49 Author: Santiago Sempere Simulation of dislocations in crystals a dislocation. From section 2.2.4 we know that the Peierls stress, τp, is defined as the maximum slope of the energy profile. Therefore, for these cases we obtain: τ18242at. p= 3.3756694 MPa;τ2240at. p= 5.2336159 MPa (4.5) These values are rather small and agree with the observed in method 1 and in the finite temperature simulations, where almost no lattice resistance opposes to the movement of the dislocations. However, an intriguing fact arising from the dislocation position plot makes us doubt on these results. Fundamentally, the Peierls stress is defined as the shear stress at which the dislocation starts to move at T= 0K. In the observed plots, the dislocation moves for values of the shear stress ∼1.1MPa. This value is much smaller than the global maximum, but it ”coincides” with the local maximum. Moreover, it is in better agreement with method 1, so we take this value as the Peierls stress obtained from method 2. In Wang’s et al. work [48], they fit the profile of the energy to a mathematical function related with the elasticity theory. In our case, although the profile of the energy exhibits a periodic curve, we have not been able to fit the data to a suitable function. Choicely, as explained, we have relied the observations on the fundamentals of the Peierls stress relating the shear stress values to the motion of the dislocation. Figure 4.10: Representation of the Energy (eV), the shear stress σxz (MPa) and the position of the dislocation (˚ A) as a function of the strain for a system of 18424 atoms. Analoguely to Method 1, we calculate, via the Orowan formula, if the theoretical expectations match the experimental results. In this case the mean dislocation motion x= and the box sizes are the same as before, Lx= 26band Lz= 14c= 14 ·1.633b. Therefore, =xbρb=xb LxLz (4.6) 50 Author: Santiago Sempere Simulation of dislocations in crystals where substituting the values obtained in our simulations we get = 0.02523 (4.7) In our plot, we obtain a strain of 0.0629 to move the mentioned distance, which is far from the theoretical value. The strain for which the dislocation starts to move does coincide with the strain observed in Method 1. A possible inaccuracy that led us to wrong results might be the forces of the periodic images and the fact of having a dislocation at the boundary. In this sense, to compute the size and the image force effects we have performed the same analysis for three different systems, the first one consisting of a system of 1368 atoms with a single dislocation, the second one is a quadrupole of dislocations, as shown in Fig. 4.11a, with 5472 atoms and a third much bigger one with 247680 atoms containing a quadrupole too. The introduction of a quadrupole is the unique way to vanish the periodic image forces in the dislocation cores. Nonetheless, after the first tilt (a) View of the initial system containing the quadropole from the (1010) plane. The introduction of the quadropole is due to the minimization of the image effects [49]. Notice the arrangement of the four dislocations, which can not be reproduced by means of a single dipole. (b) View of the relaxed system containing the quadrupole from the (1010) plane. The dislocations have recombined and formed a wide stacking fault where atoms are arranged in a FCC lattice. The initial separation between the dislocations is very small and when we apply a small stress the dislocations tend to recombine. Figure 4.11: Quadrupole of dislocations in a box containing 5472 atoms and with Lx= 20a. of the box, the dislocations seem to merge into one big stacking fault, recombining and disappearing, as shown in Fig. 4.11b. Obviously, this deceiving behavior disables us to compute the value of the Peierls stress. If we wanted to perform the computation this way, we would need a huge system, where the partial dislocations are too far from each other to interact. In Fig. 4.12, we plot the Differential displacement analysis for the quadrupole case. We would have expected four peaks in the differential displacement in the xdirection, 51 Author: Santiago Sempere Simulation of dislocations in crystals Figure 4.12: Plot of the relative displacement and differential displacement in the planes where the dislocations were introduced. The green curve depicts the DD and relative displacements of the planes of the dislocations of the upper part of the box, whereas the blue one describes the distortion of the atoms in the lower dislocation since there are four partial dislocations, but only two appear. We assume this is due to the interaction between the dislocations, that somehow balance the displacement of the atoms. In the ydirection we obtain the expected results, with four peaks and a value of zero for the screw component of the Burgers vectors (the integration is null). To get an estimation of the size error it would have been interesting to perform the same simulations but introducing a dipole rather than a quadrupole. We know this scheme will be more sensitive to the image forces than the quadrupole, but less than a single dislocation. Once again, this has been imposible due to time constraints. 4.1.2 Finite temperature properties In this section, we show the results obtained from the simulations of a system with a single dislocation at a finite temperature. The system is relaxed from the initial configuration, as introduced previously. The boundary conditions applied to the simulations are PBC in the glide plane axes, i.e xand y, and in the zdirection. However, the simulation box is enlarged in the +zand −zedges to create a vacuum gap between the vertical images. This procedure is done to avoid the creation of a second dislocation in the upper and lower edges that would interact with the principal dislocation. We have performed different simulations using different ensembles and different box sizes, Lx= 12a(2240 atoms) and Lx= 26a(18424 atoms). The temperatures at which we have carried out the simulations are rather low, 25Kand 50K. For the NVT and NPT ensembles, the temperature is fixed through the simulation and, in the case of the NPT, the pressure is set and fixed to zero as well. Also, the shear stresses are forced 52 Author: Santiago Sempere Simulation of dislocations in crystals to be zero in the latter case. In the NVE simulation, we set the initial velocities to be the corresponding ones to a kinetic energy of 25K, by the equipartition theorem. Then, at the beginning of the simulation, half of this kinetic energy transforms in potential energy, since the total energy (E) of the system must be conserved in this ensemble. Therefore, the decrease in the kinetic energy provokes the temperature to drop to a half of its initial value, 12.5Kand 25K, respectively. The simulations are performed with a timestep of ∆t= 1 fs and 800000 steps long, resulting in simulations of 800ps. The results obtained for each simulation are exhibited in Figs.4.13, 4.14, 4.15 and 4.16. Figure 4.13: Representation of the position of the dislocation versus the time expressed in simulation time steps. The system is formed by 2240 atoms and the temperature of the system is fixed to 25 K through all the simulation. The simulation has been performed using three different ensembles, NVT,NPT and NVE. In the case of the NVE ensemble, the system reaches the equilibrium at a temperature of 12.5 K Figure 4.14: Representation of the position of the dislocation versus the time expressed in simulation time steps. The system is formed by 2240 atoms and the temperature of the system is fixed to 50 K through all the simulation. The simulation has been performed using three different ensembles, NVT,NPT and NVE. In the case of the NVE ensemble, the system reaches the equilibrium at a temperature of 25 K 53 Author: Santiago Sempere Simulation of dislocations in crystals Figure 4.15: Representation of the position of the dislocation versus the time expressed in simulation time steps. The system is formed by 18424 atoms and the temperature of the system is fixed to 25 K through all the simulation. The simulation has been performed using three different ensembles, NVT,NPT and NVE. In the case of the NVE ensemble, the system reaches the equilibrium at a temperature of 12.5 K Figure 4.16: Representation of the position of the dislocation versus the time expressed in simulation time steps. The system is formed by 18424 atomsaverageand the temperature of the system is fixed to 50 K through all the simulation. The simulation has been performed using three different ensembles, NVT,NPT and NVE. In the case of the NVE ensemble, the system reaches the equilibrium at a temperature of 25 K The position of the dislocation has been monitored with three different estimators to give consistency and robustness to the analysis. For the case of the differential displacement (3.3.3), due to its sensitiveness to the thermal fluctuations, we average the positions of the atoms over five timesteps, achieving very decent results. The other two methods of analysis are the study of the nearest neighbors (3.3.4) and the computation of the common neighbor analysis (3.3.5). The analysis of the nearest neighbors gives the blue curve in the plots, which is in the middle of the position of the partial dislocation cores. The reason for this is that, for the sake of simplicity, we have averaged the position of all the atoms that exhibited a different number of nearest 54 Author: Santiago Sempere Simulation of dislocations in crystals neighbors than twelve, obtaining the center of the stacking fault as a result. Due to the variation of the volume in the NPT simulation, the computation of the nearest neighbors gives erronous results, so in these cases is better to look at the other two estimators. Notice the movement of the dislocation when applying no shear stress at all. Cautiously, we regard the residual stresses that might affect the dislocation. In Fig. 4.17 we represent the temporal evolution of the σxz which is the candidate to be responsible for the movement of the dislocation. As observed, for the NVT and NVE ensembles it reaches values of 1.5MPa in Fig.4.17a and of almost 4MPa in Fig. 4.17b. Comparing these maximum values to the approximations of the Peierls stress obtained in section 3.1.1 it is consistent that there is movement under these residual stresses. Moreover, we want to emphasize the case of the NPT ensemble, where the stresses are very close to being null but we still observe some motion. This manifests the very low resistance of the lattice to dislocation movement since only with the energy coming from the thermal influence is sufficient to make the dislocation move. (a) System of 18424 atoms (b) System of 2240 atoms Figure 4.17: Temporal evolution of the σxz component of the shear stress matrix at T= 25K In addition, the stress values for the small (Fig. 4.17b) system are greater than in the big system (Fig. 4.17a). This is in good agreement with the fact that the periodic image forces are greater in the small box than for the bigger system. The stacking fault bounded by the dislocations appear to attain an equilibrium width after a transient period. The mean values for each of the ensembles, box sizes and temperatures are shown in tables 4.1 and 4.2 55 Chapter 5 Conclusions In this thesis, we have performed atomistic simulations in the same system by using both classical and quantum approaches. We have dealt with an hcp lattice structure containing a dislocation -or in specific cases more than one-, which is a line defect that mediates the plastic behavior of materials. Firstly, we have studied the behavior of the dislocation at the classical regime, both at zero and at finite temperatures. To do so, we have made used of an already built code (LAMMPS) that has permitted us to focus more on the results rather than in the Molecular Dynamics code implementation. Applying several methods of relaxation, we have achieved the most stable, and thus the most likely, configuration of the system at zero temperature. From this state, we have applied and developed some techniques to study the plastic features of the solid containing the dislocation employing the study of the Peierls stress. Secondly, we have performed numerical simulations of a similar configuration in the quantum regime by means of a Path Integral Monte Carlo method (PIMC). We have been able to reproduce some of the results obtained by other researchers in very recent reports such as Borda et al. [19]. Due to time constraints, we have not been able to fully extend our simulation study to the regime of low temperatures, in which interesting phenomena involving dislocations are expected to occur. We leave this piece of research as future work (see Chapter 6). Next on, we report the main conclusions obtained at each of the parts of the thesis 5.1 Classical regime 5.1.1 Relaxation As a first contact with the dislocations, we have relaxed the system utilizing several techniques in order to obtain the most stable state, i.e. the state of minimum energy and forces, after the extraction of a semiplane to create the dislocation in the lattice. We have tested how the final system can be quite sensitive to some details in the relaxation procedure such as the shifting of the potential or the relaxation of the box to minimize the forces and stresses. At the relaxed state, we have recovered a well- Author: Santiago Sempere Simulation of dislocations in crystals known structure for hcp lattices with a dislocation, the dissociation of the dislocation into two Shockley partial dislocations bounding a stacking fault. Consequently, the atoms in the stacking fault region show up to be arranged in a fcc lattice. These results have been possible by the development of an algorithm to study the disregistry of the atoms, what it is often known as the relative displacement. As explained, the derivative of this function gives the corresponding edge and screw components of the Burgers vector, in the present case, (bx, by). Therefore, we observe how the initial dislocation splits up into two partial dislocations as follows: 1 3[1120] →1 3[1010] + 1 3[0110] (5.1) where 1 3[1120] is the initial Burgers vector of an edge dislocation with no screw component, and it splits up into two dislocations with individual non-zero screw component. The width of the stacking fault, which depends on the elastic properties of the crystal, has shown up to be very broad. Its theoretical width is out of our scope and, thus, we have observed that every time we enlarged the size of the box that contained the Burgers vector (in our case, Lx) we achieved a different, greater, width for the stacking fault region. 5.1.2 Peierls stress We have implemented two different techniques to retrieve the value of the Peierls stress for the classical system of Xe. We have observed that comparing both methods we have not been able to achieve a final estimation of the Peierls stress, but we are capable of giving a guess in the order of magnitude of it, which seems to be enough to explain some of the results obtained at finite temperature. A consistent order of magnitude with the simulations would be τp∼1MPa. This value is rather low and it is the cause of the motion of dislocations at finite temperature in the absence of external forces, only due to the thermal fluctuations. The computation of the shear stress by means of the first method developed is very sensitive to the box size and the upper slab size. In our simulations, it has shown up to give a very irregular pattern for the small boxes, when we expected a regular periodical profile as reported in many cases, for instance, see [24]. This is not a problem for classical simulations, since we have computational resources to overcome these difficulties, and, therefore, this method can provide a good estimation for the Peierls stress, as shown in Fig. 4.9 for the system of 300000 atoms. Nonetheless, for the quantum simulations, due to its intrinsic complexity, it seems unfeasible to perform simulations of more than a few thousand of atoms. On the other hand, we have developed a method to compute the Peierls stress that, from a qualitative point of view, offers better results than the former one. In this approach, we recover a very smooth profile for the energy when tilting the box more and more. We certainly believe this method, when studied the errors induced by the size effects, is the most suitable to study the Peierls stress at zero temperature in the quantum regime. 63 Author: Santiago Sempere Simulation of dislocations in crystals 5.1.3 Finite temperature In chapter 4.1.2, we have simulated the previously relaxed system of Xe at certain finite temperatures, T= 25Kand T= 50K, with no stresses or external forces applied. We had the purpose of simply confirming our expectations of observing no motion at all. However, we observe movement of the dislocation under no stresses, other than the residual ones, manifesting the absence of almost any intrinsic lattice resistance for the system. We have distinguished two sources for this motion, depending on the ensemble used for the simulation. In the NVT and NVE ensembles, we have observed that a large residual stress, of the same order of magnitude of the Peierls stress τp, fluctuates in time, causing the dislocation to fluctuate as well. In the other hand, in the NPT ensemble, where we force the stresses to be as close to zero as possible, we conclude the motion is due to thermal fluctuations, which led the dislocation to move to one side or the other alternatively. This is an important conclusion since gives an evidence that the absence of lattice resistance to dislocation motion is not an exclusive feature of quantum systems, but it has more to do with the kind of interaction, i.e. the strength of the interatomic potential, between particles. Moreover, we observe the width of the dislocation remains constant throughout the simulation and its value coincides with the initial width, obtained from the relaxation. 5.2 Quantum regime Motivated by recent experiments done in the field, we have computed the behavior of a single dislocation in a 4He solid at very low temperatures, i.e T= 1.5Kand T= 3K. Although not being able to study thoroughly the dislocation itself, we have been able to compute the properties of the system with a dislocation. We have observed no sign of superfluidity, either in the cores of the dislocation, as reported by Boninsegni et al. [11], or the stacking fault region. Although, a possible cause why we have not observed either of those might be the high temperature at which we performed the simulations or their small length. In Boninsegni’s et al. work, they span a temperature range from 0.2 to 1K, whereas, in Borda’s et al., they simulate the system at a temperature of T= 0.267K. These values are far from ours, 1.5Kand 3K. Due to many doubts regarding the structural analysis of the system at T= 3K, we questioned if the solid had melted into a crystal or not. From the analysis of the radial distribution functions and the static structure factor, we concluded that it was still a solid for that temperature. These issues prevented us from doing a much longer simulation that could give us reliable and meaningful results. The results provided by the on-the-fly estimators have shown up to be non-sense. Therefore, we have concluded that to develope a proper analysis of the dislocation motion and to be able to track the dislocation throughout the simulation, we have to average the positions of the beads over several time steps to reduce the zero-point motion effect. 64 Chapter 6 Further work Throughout the thesis we have been able to cover most of the objectives proposed related to the classical regime, however, due to the appearance of unexpected results in this first part, we have not been able to perform an exhaustive analysis in the quantum regime. Therefore, in this section, we present the aspects that require further work and a path to follow in order to achieve the objectives we proposed at the beginning. Concerning the quantum simulations, we have been able to make an approach to the simulation of a single dislocation in a solid of 4He at low temperatures. Our first objective was to make a study of the behavior of the dislocation under quantum effects and to clarify the controversy of the supersolidity; however, many improvements must be done to extract conclusions from our simulations. Firstly, we should be able to reimplement the estimators in the PIMC code in order to make them work every some iterations after we have averaged the beads positions. The zero-point motion of the particles distorts the structure of the crystal, disabling us to compute the properties of the dislocation. Secondly, we should review the simulation of a single dislocation with PBC. It would be logical, if possible, to better simulate a dipole of dislocations to avoid problems at the boundaries. With these improvements and performing a larger simulation, we would be able, in principle, to achieve the same results as Borda et al. concerning the dislocation behavior, but better results upon the superfluidity of the stacking fault and the dislocation cores. In our PIMC code, it has been implemented the Worm algorithm that allows us to properly calculate the superfluidity for systems greater than a few hundreds of atoms. In Borda et al., they do not make use of this algorithm (at least they do not mention it), so that is why we doubt in some of the technical aspects of their procedure and we think we could make a better study, at least concerning the quantum aspects of the simulation. Moreover, it would be interesting to span a greater range of temperatures, trying to go under a few dK. Bear in mind as a useful reference that the giant plasticity phenomenon has been reported to appear at T≃0.15K, so it would be interesting to be able to achieve this temperature to contrast the experiments with the simulations. Moreover, it would be important to perform a quantum static analysis of the system with a dislocation. This is, to calculate the properties at T= 0K. The Path Integral Ground State method allows us to perform such analysis, that would be very useful to assess the Peierls stress for the quantum system and to depict the stacking fault energy Author: Santiago Sempere Simulation of dislocations in crystals surface γas well. The study of both properties would give us an idea of the lattice resistance to motion of dislocation, and, therefore, the resistance to plastic deformation and the stability of the dissociation of the edge dislocation into two Shockley partial dislocations. In the end, we would be able to give a proper theoretical explanation to the facts observed in the experiments, such as Giant Plasticity, from a fundamental point of view. 66 Bibliography [1] A. F. Andreev and I. M. Lifshtiz. Sov. Phys. JETP , 29 , 1107, 1969. [2] G. V. Chester. Phys. Rev. A 2, 256–258, 1970. [3] E.L Andronikashvili. Zh. Eksp. Teor. Fiz. 16, 780, 1946. [4] E. Kim and M.H.W. Chan. Nature 427, 225, 2004. [5] A. Leggett. Physica Fennica, 8, 125, 1973. [6] J.Day and J. Beamish. Nature 450, 853, 2007. [7] J. Reppy. Condensed matter: The supersolid’s nemesis. E. S. Reich, Nature 468, 748 (2010). [8] J. Reppy. Phys. Rev. Lett. 104, 255301, 2010. [9] D.Y. Kim and M. H. W. Chan. Phys. Rev. Lett. 109, 155301, 2012. [10] A. Haziot, X. Rojas, A. D. Fefferman, and J.R. Beamish ans S.Balibar. Phys. Rev. Lett. 110, 035301, 2013. [11] M. Boninsegni, A. B. Kuklov, L. Pollet, N. V. Prokof’ev, B. V. Svistunov, and M. Troyer. Phys. Rev. Lett. 99, 035301, 2007. [12] A.B. Kuklov. Phys. Rev. B 92, 134504, 2015. [13] A.B. Kuklov, L. Polllet, N.V. Prokof’ev, and B. V. Svistunov. Phys. Rev. B 90, 184508, 2014. [14] M.W Ray and R.B. Hallock. Phys. Rev. Lett. 100, 235301, 2008. [15] M.W Ray and R.B. Hallock. Phys. Rev. B 79, 224302, 2009. [16] M.W Ray and R.B. Hallock. Phys. Rev. B 84, 144512, 2011. [17] M. W. Vekhov, W. Mullin, and R.B. Hallock. Phys. Rev. Lett. 113, 035302, 2014. [18] M. W. Vekhov and R.B. Hallock. Phys. Rev. B 90, 134511, 2014. [19] E. Borda, W. Cai, and M. de Koning. Phys. Rev. Lett. 117, 045301, 2016. [20] A Serra, D.J Bacon, and R.C Pond. Acta Metallurgica 36, Pag. 3183-3203, 1988. [21] V. V. Bulatov and W. Cai. Computer simulations of dislocations. Oxford University Press, Oxford, 2006), 2006. Author: Santiago Sempere Simulation of dislocations in crystals [22] D.Hull and D.J Bacon. Introduction to Dislocations. BH, 2011. [23] D.J. Bacon and M. H. Liang. Philosophical Magazine A, 53:2, 163-179, 1986. [24] H.A. Khater and D.J. Bacon. Acta Materialia, 2978-2987, 2010. [25] A Berghezan, A Fourdeux, and S Amelinckx. Acta Metallurgica, 9, 464, 1961. [26] R. Pessoa, S. A. Vitiello, and M. de Koning. Phys. Rev. Lett. 104, 085301, 2014. [27] E. Polak and G. Ribi`ere. Revue Francaise d’Informatique et de Recherche Op´erationnelle, 16, 35–43., 1969. [28] W. Cai, Li J., and Yip S. Comprehesive Nuclear Materials, pages 249–265, 2012. [29] D.J. Bacon and J.W. Martin. Philosophical Magazine A, 43:4, 883-900, 1981. [30] Murray S. Daw and M. I. Baskes. Phys. Rev. B 29, 6443, 1984. [31] Riccardo Rota. Path Integral Monte Carlo and Bose-Einstein condensation in quantum fluids and solids. PhD thesis, Universitat Polit`ecnica de Catalunya, 2011. [32] Bernard Bernu and David M. Ceperley. Path integral monte carlo. In J. Grotendorst, D. Marx, and A. Muramatsu, editors, Quantum Simulation of Complex Many-Body Systems: From Theory to Algorithms, pages 51–61. John von Neumann Institute for Computing, 2002. [33] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. J. Chem. Phys , 21 , 1087, 1953. [34] H. Trotter. Proc. Am. Math. , 10 , 545, 1959. [35] D. Ceperley. Rev. Mod. Phys. 67, 280, 1995. [36] M. Takahashi and M. Imada. J. Phys. Soc. Jpn. , 53 , 963, 1984. [37] M. Takahashi and M. Imada. J. Phys. Soc. Jpn. , 53 , 3765, 1984. [38] S. A. Chin and C. R. Chen. J. Chem. Phys , 117 , 1409, 2002. [39] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn. Physics Letters A, 238, 253, 1998. [40] M. Boninsegni, N. V. Prokof ’ev, and B. V. Svistunov. Phys. Rev. E , 74 , 036701, 2006. [41] M. Boninsegni, N. Prokof ’ev, and B. Svistunov. Phys. Rev. Lett. , 96 ,070601, 2006. [42] A. Stukowski. Modelling Simul. Mater. Sci. Eng. 20, 4, 2012. [43] Honeycutt and Andersen. J. Phys. 91 (19), pp 4950–4963, 1987. [44] H. Tsuzuki, P.S. Branicio, and J. P. Rino. Comp. phys. Comm. 177, 518-523, 2007. [45] S. Plimpton. Fast parallel algorithms for short-range molecular dynamics. J Comp Phys, 117, 1-19, 1995. http://lammps.sandia.gov. [46] K. Momma and F. Izumi. J. Appl. Crystallogr., 44, 1272-1276, 2011. 68 Author: Santiago Sempere Simulation of dislocations in crystals [47] Yu N. Osetsky and D. J. Bacon. Modelling Simul. Mater. Sci. Eng. 11, 427-446, 2003. [48] G. Wang, Alejandro Strachan, T. C¸agin, and W. A. Goddard. Modelling Simul. Mater. Sci. Eng. 12, S371-S389, 2004. [49] S. K. Yadav, R. Ramprasad, A. Misra, and X. Liu. Acta Materialia 74, pag. 268-277, 2014. [50] C. Cazorla and J. Boronat. Phys. Rev. B 88, 224501, 2013. [51] M. C. Gordillo and D. M. Ceperley. Phys. Rev. Lett. 79, 3010, 1997. [52] R. A. Aziz, F. R. McCourt, and C. C. Wong. Mol. Phys. , 61 , 1487, 1987. 69