Calculation of thermal conductivity of a graphene layer

[Fishpack]

Hello, everyone. I have tried to figure out the thermal conductivity of a graphene layer with the unit of “W/mK” (not W/K) using NEMD methods.

The structure is a graphene nanosheet which has the size of 12 nm × 9 nm. The total number of the carbon atoms is approximately 2240.

For the equilibrium, the structure was optimized and simulated using Nose-Hoover thermostat at 300 K. In this simulation, the time step was 1 fs and the whole process time was 4 ns. The thermostat timescale was 100 fs. The reservoir temperature and final temperature were 300 K. The tersoff_C_2010 potential was applied.

The non-equilibrium simulation was then performed using Nose-Hoover thermostat at 300 K with the same condition of the equilibrium simulation (ex. time and temperature). Also, The heat source and heat sink were set at the left side and right side of the structure. Using the non-equilibrium momentum exchange in ATK, I want to make the temperature gradient in the area between the heat source and heat sink. Then, I will calculate the thermal conductivity of the graphene layer.

I have several questions about the simulation.
1)    I know that the nanosheet is 2D and periodic in the y and z directions. Is it required to set the periodic boundary condition for the top and bottom side additionally? If the boundary condition is periodic, there are no scattering effects related on the phonon or electrons because the top and bottom side are not terminated anymore, Right? I expect for the case; there are no dependencies on the width of the structure for the thermal conductivity. I want to know about this.

2)    For the non-equilibrium momentum exchange, should we set the fixed boundary at the left and right side of the structure besides the heat sink and source?

3)     For the Nose-Hoover thermostat, the temperature gradient can be generated using reservoir temperature. Is there are no problem to obtain the thermal conductivity of the graphene layer although some synchronization issues are between the heat source and sink.

4)    I read a several papers about NEMD for the calculation of the thermal conductivity by LAMMPS. In that case, the method is similar methods kinds of setting the different reservoir temperature on both sides for the Nose-Hoover thermostat.

5)    Please comment about the procedure of simulation, some notes and missing point and so on.

[Petr Khomyakov]

You may have a look at our tutorial on the thermal conductance calculations at http://docs.quantumwise.com/tutorials/interfacial_thermal_conductance/interfacial_thermal_conductance.html. There is also our recent webinar video on this topic, posted on youtube, see https://www.youtube.com/playlist?list=PLkOUx2yxHrqxH5szv2HUX0JwqD25Y2Ab8&utm_source=ZohoCampaigns&utm_campaign=Send+links+to+MD+webinars+participants_2017-03-06&utm_medium=email. 

You may also be interested in watching the webinar video on molecular dynamics simulations in the ATK, see https://www.youtube.com/playlist?list=PLkOUx2yxHrqweTQGgrgiH5wC_WlHKwq5i&utm_source=ZohoCampaigns&utm_campaign=Send+links+to+MD+webinars+participants_2017-03-06&utm_medium=email.

[Julian Schneider]

1)    I know that the nanosheet is 2D and periodic in the y and z directions. Is it required to set the periodic boundary condition for the top and bottom side additionally? If the boundary condition is periodic, there are no scattering effects related on the phonon or electrons because the top and bottom side are not terminated anymore, Right? I expect for the case; there are no dependencies on the width of the structure for the thermal conductivity. I want to know about this.

With ATK-Classical you cannot switch off the periodic boundary conditions, so you should make sure that you have a sufficiently large vacuum buffer on top of th graphene layer if you want to simulate an isolated layer.

2)    For the non-equilibrium momentum exchange, should we set the fixed boundary at the left and right side of the structure besides the heat sink and source?

You have two options how to set up the boundaries in transport direction:
a) Add some vacuum buffer in transport direction and fix the outermost terminating atoms. This way you'd only have periodicity in the in-plane direction perpendicular to the transport direction.
b) You treat the system periodic in the transport direction, as well. Then the thermal current runs in both directions. Here you in principle don't need to fix anything.

3)     For the Nose-Hoover thermostat, the temperature gradient can be generated using reservoir temperature. Is there are no problem to obtain the thermal conductivity of the graphene layer although some synchronization issues are between the heat source and sink.

In principle you could just set two thermostats on heat sink and source to obtain the gradient. However, in this case you'd have no straightforward way to measure the magnitude of the thermal flux. So, I'd recommend using the NonEquilibriumMomentumExchange method.

4)    I read a several papers about NEMD for the calculation of the thermal conductivity by LAMMPS. In that case, the method is similar methods kinds of setting the different reservoir temperature on both sides for the Nose-Hoover thermostat.

According to recent papers, the alternative way to using NonEquilibriumMomentumExchange  not really setting reservoir temperatures, but to set achieve a specified thermal by heating and cooling the heat and sink with a constant heating power, which at steady state has the same magnitude as the thermal curent in the system. This method will be available in the upcoming ATK-2017 release.

5)    Please comment about the procedure of simulation, some notes and missing point and so on.

To obtain reliable bulk conductivity values, you may have to run a series of simulations at increasing system sizes and then extrapolate to infinite lengths, as described e.g. in http://aip.scitation.org/doi/abs/10.1063/1.4767516.
As, Petr said, for more explanations and details, you may want to look at our tutorials and the webinar.

[Fishpack]

Now, I've simulated the NEMD simulation based on your posts and the tutorials. That was helpful for me. Thank you.
By the way,  Is it acceptable that the kinetic energy and the temperature increases as the time step goes up under the non-equilibrium momentum exchange?
We want to calculate and obtain the thermal conductivity of a graphene layer at 300 K. But, when we increase the time step to find the steady-state condition, the kinetic energy and the temperature also increases linearly.
I attached my python script for the simulation.
Thank you

[Julian Schneider]

The kinetic energy should not increase significantly  during the simulation. By increasing the time step, do you mean the total number of MD steps or the actual time step  (e.g. 1 fs)?
The latter should not be increased. Increasing it to e.g. 5 fs can in fact lead to energy conservation being violated.
If you just increased the number if MD steps and the energy increases, then the system may not be sufficiently equilibrated before the NEMD simulation. Or your temperature gradient has become so large that part of the system experiences such a high temperature that it runs out of equilibrium.

[Fishpack]

I meant that the number of MD steps not the actual time step. In my simulation, I used the 1 fs time step. And I think that the system is equilibrated enough because the temperature is almost constant at 300 K.

I want to ask something

1. What do you think the suitable momentum exchange interval? I used 150 intervals for the simulation. When we used 150 intervals, the temperature was increased.

2. How can we decide the convergence or steady state of the temperature profile? Are there any criteria for that? For my case, the temperature gradient is changed in a short time, and it is hard to figure out the convergence of the temperature profile.

3. In the tutorial of the calculating interfacial thermal conductance using molecular dynamics, you compared the two temperature profiles. Are there a proper time range for the integration of the temperature profile? The temperature gradient is changed when we modify the time range.

4. Can you tell me the reason why the increase of the system temperature during the simulation of the non-equilibrium momentum exchange? In the case, the system is isolated, and no energy can be provided.

5. What is the meaning of "truncated content"? I saw the warning in my log file of the simulation.

[Julian Schneider]

Sorry for the late answer, I took me some to look into your system.

1. What do you think the suitable momentum exchange interval? I used 150 intervals for the simulation. When we used 150 intervals, the temperature was increased.

150 should be fine in general, although it depends a bit on your system. Generally, the exchange interval allows you to modify the magnitude of the thermal current. For a system with a small cross section, the same thermal current can lead to a larger thermal flux and thus a larger temperature gradient. The temperature gradient should ideally not be too large, on the other hand large enough so that it can be measured with sufficient accuracy. 

2. How can we decide the convergence or steady state of the temperature profile? Are there any criteria for that? For my case, the temperature gradient is changed in a short time, and it is hard to figure out the convergence of the temperature profile.

You should plot your temperature profile at different time intervals during the simulation. If you don't see a difference between the temperature profiles of subsequent time intervals then you can assume that your temperature profile has reached a steady-state.

3. In the tutorial of the calculating interfacial thermal conductance using molecular dynamics, you compared the two temperature profiles. Are there a proper time range for the integration of the temperature profile? The temperature gradient is changed when we modify the time range.

At first, I would recommend you to use a smaller log interval than 10 000, e.g. something like 500 or 1000. This will make the temperature profile a lot smoother, as you will have more snapshots to average over. Remember, for your system you have to calculate the profile along the B-direction, not the default C-direction.
Yes, the temperature profile will depend on the chosen time interval. On the one hand, as for a small time interval your average might not be very good and the profile can become fuzzy. On the other hand, as long as you haven't reached a steady state, changing the length of the time interval will also include snapshots from different phases of the simulation, which will result in a different profile, that's just the nature of the non-equilibrium simulation. So, as a rule of thumb, you could use intervals of 50-100 ps to check for convergence, and once you know that you have reached a steady-sate, you should average the profile over several 100 ps, to get a smooth profile.
If your profile is too fuzzy you can also try and increase the bin width a bit.

4. Can you tell me the reason why the increase of the system temperature during the simulation of the non-equilibrium momentum exchange? In the case, the system is isolated, and no energy can be provided.

I can reproduce your problem, and it actually seems that the total energy increases during time. Although I'm not 100% sure, I would guess a major part of the problem is that the time step is too large and the conservation of total energy is not given any more. In this case a smaller time step, maybe 0.5 fs, might at least fix the problem to some extent. Furhtermore, the system is very soft and buckles quite a lot during the simulation which might make it additionally difficult.

5. What is the meaning of "truncated content"? I saw the warning in my log file of the simulation.

That is just a message that the output line in the log file is longer than the width of the terminal. You can ignore that.

[Fishpack]

Thank you for your reply.

I performed the non-equilibrium momentum exchange simulation for 1 ns.

The time step was modified with 0.5 fs. And the log interval was 2000.

Still, we could see the increase of the temperature and the total energy.
(I think that the slope of the increase of the temperature is reduced when we compared to the 1 fs time step)

Please, see my attached files that contain the graph of the temperature and the total energy.
One of them is the case that we applied the 200 momentum exchange interval.
and the other is the case of the 400 momentum exchange interval.

1) Is it acceptable increase of the temperature in the graph? If not, I want you to give me some advise to fix them.

2) I expected the unbalance of the energy exchange between the heat source and the heat sink might cause this problem. Is there are method to obtain the heat flux at the heat source and the heat sink? and are there options at the heat source and the heat sink to modify the heat flux of each part?