Simulating the Damaged Tungsten for Thermal Conductivity

In my previous blog post, I mentioned using classical molecular dynamics implemented in LAMMPS. Let’s dig a little into the procedure and learn how we calculate the thermal conductivity of damaged Tungsten.
Although there are a massive collection of things you can do with LAMMPS, our procedure involve three main steps, and these are:
- Introducing damage to simulate irradiated Tungsten / Creation-Relaxation Algorithm(CRA)
How would you imagine a neutron bombardment to a material? Neutrons hit the surface, and the impact will spread through the material, displacing the atoms from their natural sites to a higher energy position. This is one of the ways to simulate the damage. However, the time required to observe the damage on a broader scale and for high doses is pretty lengthy. This also means that the use of computational resources is extremely high. Luckily, Creation-Relaxation Algorithm is a shortcut to simulate high-dose damage swiftly. In this iterative method, randomly selected atoms are displaced such that there is a vacancy in the initial position of the atom and an interstitial(Figure 1). These pairs are also called Frenkel pairs. Then, energy minimization is performed, and the atoms settle to their new positions. Finally, these steps are repeated until the desired damage level is reached. In our study, we applied damage to the middle 10% of the simulation box(Figure-2).

2. Thermally equilibrating the system:
Since we will investigate the thermal conductivity, it’s time to introduce the temperature and pressure factor to our damaged system.
- Isothermal-Isobaric Ensemble (NPT): In this thermodynamical ensemble keeping the number of particles(N), pressure(P), and temperature(T) constant, we increase the temperature of the system to the desired value. While this happens, the volume of the simulation box changes until it converges to a value for a given pressure and temperature.
- Canonical Ensemble (NVT): After obtaining the well converged system from the results of NPT, we put the system into canonical ensemble where the number of particles(N), volume(V), and temperature(T) are constant. The main reason we apply this step after thermally equilibrating the system is to obtain a statistically probable state of the system under given conditions.
- Introducing temperature gradient in microcanonical ensemble(NVE): In this step, we construct a heat source(hot region) and heat sink(cold region) within the simulation box. This is achieved by withdrawing kinetic energy from the heat sink and supplying the same amount of kinetic energy to the heat source(Figure 2). To obtain a steady state temperature gradient between “hot” and “cold” regions, we let the energy be distributed along the simulation box in microcanonical ensemble. As you can interpret the constants of the system from its abbreviation, the number of particles(N), volume(V), and total energy of the system(E) are held constant.

In step-2 and step-3, the interstitial atoms can fall into the vacant sites, or new atoms can leave their sites and create a Frenkel pair. You can see an animation in the featuring image taken from one of our calculations in step-3. If you want to learn more about how we calculate thermal conductivity, you can check out this blog post from Paolo, who was one of the last year’s participants in the project!
In my final blog post, I will share my results with you!
Leave a Reply