Green Kubo results processing

[ashwani12-cmd]

Hi everyone,

I’m currently computing lattice thermal conductivity using Equilibrium Molecular Dynamics (EMD) with the Green–Kubo (GK) formalism in GPUMD, employing a NEP (Neural Equivariant Potential) model. I would really appreciate feedback from the community on whether my simulation protocol and post-processing workflow look physically and numerically sound, or if there are aspects that should be improved.

I want to calculate the x , y and z thermal conductivity in separately.

my run.in file:

potential ./nep_y2025_m12_d20_h16_m08_s16_generation100000.txt

minimize fire 1.0e-5 100000 1

ensemble nve
time_step 0
dump_xyz -1 0 1 relaxed.xyz
run 1

Equilibration

velocity 300
time_step 1
ensemble npt_ber 300 300 100 0.0 0.0 0.0 0.0 0.0 0.0 60.84 60.84 73.12 17.82 17.82 10.89 1000
dump_thermo 1000
run 1000000

Production

ensemble nve
dump_thermo 10000
compute_hac 10 100000 1
run 10000000

My post processing python script :

import sys
sys.path.append('/home/IITB/multiscale-mechanics/ashwani12/software/tools/gpyumd')

import numpy as np
import matplotlib.pyplot as plt
from gpyumd.load import load_hac
import os

num_runs = 1

============================================================

Plot style (journal quality)

============================================================

plt.rcParams.update({
"font.size": 16,
"axes.linewidth": 2.0,
})

def style(ax):
ax.tick_params(which='major', length=8, width=2,
direction='in', top=True, right=True)
ax.tick_params(which='minor', length=4, width=2,
direction='in', top=True, right=True)
ax.minorticks_on()

============================================================

GK time window (plateau region)

============================================================

tau_range = range(30000, 40000) # adjust carefully!

============================================================

Load HAC data

============================================================

hac_dict = {}
for run_num in range(num_runs):
hac_dict[f'run{run_num}'] = load_hac(
num_corr_points=100000,
output_interval=1,
directory="./"
)['run0']

Time array

gk_tau = hac_dict['run0']['t']

============================================================

Initialize accumulators

============================================================

kxx_ave = np.zeros_like(hac_dict['run0']['kxi'])
kyy_ave = np.zeros_like(hac_dict['run0']['kyi'])
kzz_ave = np.zeros_like(hac_dict['run0']['kz'])

============================================================

Loop over runs

============================================================

for runkey in hac_dict.keys():

# X direction (in + out)
kxx_ave += 1 * (hac_dict[runkey]['kxi'] +
                  hac_dict[runkey]['kxo'])

# Y direction (in + out)
kyy_ave += 1 * (hac_dict[runkey]['kyi'] +
                  hac_dict[runkey]['kyo'])

# Z direction (already full)
kzz_ave += hac_dict[runkey]['kz']

============================================================

Average over runs

============================================================

nruns = len(hac_dict)
kxx_ave /= nruns
kyy_ave /= nruns
kzz_ave /= nruns

============================================================

Plateau averages

============================================================

kxx = np.mean(kxx_ave[tau_range])
kyy = np.mean(kyy_ave[tau_range])
kzz = np.mean(kzz_ave[tau_range])

print(f"kxx (GK) = {kxx:.6f}")
print(f"kyy (GK) = {kyy:.6f}")
print(f"kzz (GK) = {kzz:.6f}")

============================================================

Plot GK convergence

============================================================

fig, ax = plt.subplots(figsize=(5, 4))

ax.plot(gk_tau, kxx_ave, label=r"$\kappa_{xx}$", lw=2)
ax.plot(gk_tau, kyy_ave, label=r"$\kappa_{yy}$", lw=2)
ax.plot(gk_tau, kzz_ave, label=r"$\kappa_{zz}$", lw=2)

Plateau window

ax.axvspan(gk_tau[tau_range.start],
gk_tau[tau_range.stop-1],
color='gray', alpha=0.2, label="Plateau")

ax.set_xlabel("Correlation time (ps)")
ax.set_ylabel("Thermal conductivity (W/mK)")
ax.legend()
style(ax)

============================================================

Save & show

============================================================

plt.savefig("gk_kxx_kyy_kzz_convergence.pdf",
bbox_inches="tight", pad_inches=0)

plt.show()
plt.close()

print("✅ Plot saved: gk_kxx_kyy_kzz_convergence.pdf")

gk_kxx_kyy_kzz_convergence.pdf

[brucefan1983]

It's not easy to check if all the steps are perfect.

I try to point out one thing:

You are studying a material with low thermal conductivity, so it is better to use a small sampling interval, to make the GK integration more accurate. Also you don't need so large correlation time. So I would change to use:

compute_hac 1 100000 1

[brucefan1983]

And if you have H atom in your system, perhaps you need to reduce the time step to 0.5 fs.

[ashwani12-cmd]

My system has only Mg and Bi , that is Mg3Bi2 and can you comment on this part of my post processing script. This only for single run

============================================================

Loop over runs

============================================================

for runkey in hac_dict.keys():

# X direction (in + out)
kxx_ave += 1 * (hac_dict[runkey]['kxi'] +
                  hac_dict[runkey]['kxo'])

# Y direction (in + out)
kyy_ave += 1 * (hac_dict[runkey]['kyi'] +
                  hac_dict[runkey]['kyo'])

# Z direction (already full)
kzz_ave += hac_dict[runkey]['kz']

============================================================

Average over runs

============================================================

nruns = len(hac_dict)
kxx_ave /= nruns
kyy_ave /= nruns
kzz_ave /= nruns

So, here # X direction (in + out) kxx_ave += 1 * (hac_dict[runkey]['kxi'] +hac_dict[runkey]['kxo']) , need 1 or divide by 2. which one will be correct.

[brucefan1983]

There is no need to divide by 2.

You are studying bulk materials, so the in- and out-components for the x direction (similarly for the y direction) should be summed up.

See documentation:

https://gpumd.org/dev/gpumd/output_files/hac_out.html#hac-out