Hey everyone,
I’m working on a LAMMPS simulation to calculate thermal conductivity using the Green-Kubo method with an MLIP (machine-learned interatomic potential) for a Mg-Bi system. I’m encountering an issue where my calculated thermal conductivity is about two orders of magnitude higher than expected (e.g., ~100-1000 W/mK instead of ~1-10 W/mK for similar materials). I’ve attached my LAMMPS input script and Python processing script below for reference.
Mg3Bi2 simulation
clear
units metal
dimension 3
boundary p p p
atom_style atomic
read_data lammps_pattern.data
timestep 0.001
neighbor 2.0 nsq
variable boxx equal xhi
variable boxy equal yhi
mass 1 24.305
mass 2 208.980
MTP potential
pair_style mlip mlip.ini
pair_coeff * *
thermo 100
thermo_style custom step temp pe ke etotal press pxx pyy pzz pxy
Compute and print desired thermodynamic information
dump 1 all custom 10000 dump_min.lammpstrj id type x y z
dump_modify 1 sort id
thermo_style custom step pe ke temp pxx pyy pxy pzz xlo xhi ylo yhi fmax fnorm
thermo 1000
min_style fire
minimize 0 1e-10 100000 100000
undump 1
reset_timestep 0
###################################################################################################
Equilibration Stage 1
Assign initial velocities
variable temp equal 300.0
variable temp2 equal 2*{temp}
variable nfreq equal 10000
variable runtime equal 100*{nfreq}
velocity all create ${temp} 41663
Do NVT
dump 1 all custom 10000 dump_NVT.lammpstrj id type x y z
dump_modify 1 sort id
fix 1 all nvt temp {temp} {temp} 1.0
run ${runtime}
undump 1
write_restart nvt.restart
reset_timestep 0
###################################################################################################
Store heat flux
variable nev equal 1
variable nrep equal 1
variable nfreq equal 1
variable nfreq0 equal 10000
variable runtime equal 600*${nfreq0}
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all centroid/stress/atom NULL virial
compute myStress all stress/atom NULL virial
compute 5 all heat/flux myKE myPE myStress
fix JJ all ave/time {nev} {nrep} ${nfreq} c_5[1] c_5[2] c_5[3] file J0Jt.dat ave one
run
thermo {nfreq0} run {runtime}
Python post processing script
gk_publishable_plots_parallel.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import correlate
from concurrent.futures import ProcessPoolExecutor
import inspect
import json
from scipy.io import savemat
import time
step_times =
start_time = time.time()
def log_step(msg):
current_time = time.time()
elapsed_time = current_time - start_time
step_times.append((msg, elapsed_time))
print(f"Line {inspect.currentframe().f_back.f_lineno}: {msg} (Elapsed: {elapsed_time:.2f} s)")
def compute_k_for_lag(lag, corr_xx, corr_yy, corr_zz, nS):
norm = np.arange(nS, nS - lag, -1)
mid = nS - 1
kx = np.sum(corr_xx[mid:mid + lag] / norm)
ky = np.sum(corr_yy[mid:mid + lag] / norm)
kz = np.sum(corr_zz[mid:mid + lag] / norm)
return kx, ky, kz
def gk_publishable_plots():
log_step(“Setting constants”)
deltaMD = 0.001
vol = 133.5058665981581515*10
kb = 8.617343e-5
theta = 300.0
log_step("Reading J0Jt.dat")
with open('J0Jt.dat', 'r') as f:
_ = f.readline(); _ = f.readline()
data = np.genfromtxt(f)
log_step("Transposing and shaping input data")
A1 = data.T
nS = A1.shape[1]
log_step(f"Using {nS} time steps from J0Jt.dat")
log_step("Creating lag vector")
I = np.concatenate((np.arange(2000, 5001, 500), np.arange(10000, 100001, 500)))
tau_ps = I * deltaMD
log_step(f"Correlation window = {tau_ps[-1]:.1f} ps, step = {deltaMD} ps")
#log_step("Creating lag vector")
I = np.concatenate((np.arange(100, 5001, 500), np.arange(6000, 100001, 500)))
tau_ps = I * deltaMD
log_step("Calculating autocorrelations")
corr_xx = correlate(A1[1], A1[1], mode='full')
corr_yy = correlate(A1[2], A1[2], mode='full')
corr_zz = correlate(A1[3], A1[3], mode='full')
log_step("Running parallel conductivity computation")
results = []
with ProcessPoolExecutor() as executor:
futures = [executor.submit(compute_k_for_lag, lag, corr_xx, corr_yy, corr_zz, nS) for lag in I]
for future in futures:
results.append(future.result())
log_step("Converting results and applying factor")
kcond = np.array(results)
factor = (deltaMD / (vol * kb * theta**2)) * 1602.176
kcond *= factor
log_step("Computing in-plane average")
kinplane = np.mean(kcond[:, 0:2], axis=1, keepdims=True)
kcond = np.hstack((kcond, kinplane))
log_step("Calculating standard deviation")
std_k = np.std(kcond[:, 3])
log_step("Extracting plateau region (15-25 ps)")
idx_low = np.where(tau_ps >= 15)[0][0]
idx_high = np.where(tau_ps <= 25)[0][-1]
k = {
'kxx': np.mean(kcond[idx_low:idx_high+1, 0]),
'kyy': np.mean(kcond[idx_low:idx_high+1, 1]),
'kzz': np.mean(kcond[idx_low:idx_high+1, 2]),
'k': np.mean(kcond[idx_low:idx_high+1, 3]),
'std': np.max(kcond[idx_low:idx_high+1, 3]) - np.min(kcond[idx_low:idx_high+1, 3])
}
print(f'kxx={k["kxx"]:.2f}, kyy={k["kyy"]:.2f}, kzz={k["kzz"]:.2f}, '
f'k={k["k"]:.2f} ± {k["std"]:.2f} W/mK')
log_step("Saving conductivity data to thermal_conductivity.txt")
out_k = np.column_stack((tau_ps, kcond))
header_k = 'tau_ps[ps] k_xx[W/mK] k_yy[W/mK] k_zz[W/mK] k_inplane[W/mK]'
np.savetxt('thermal_conductivity.txt', out_k, header=header_k, fmt='%.6e', comments='')
with open('thermal_conductivity_summary.json', 'w') as jf:
json.dump(k, jf, indent=4)
log_step("Saved summary dictionary to thermal_conductivity_summary.json")
log_step("Plotting and saving thermal conductivity figure to PDF")
plt.figure()
labels = ['k_xx', 'k_yy', 'k_zz', 'k_inplane']
for i in range(4):
plt.plot(tau_ps, kcond[:, i], marker='o', label=labels[i])
plt.xlabel('Correlation Time τ_I (ps)')
plt.ylabel('Thermal Conductivity (W/mK)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('thermal_conductivity.pdf')
plt.show()
log_step("Computing normalized HCACF")
Icorr = np.arange(1, 100001, 100)
mid = nS - 1
normHC = corr_xx[mid + Icorr] / corr_xx[mid]
log_step("Saving normalized HCACF to normalized_HCACF.txt")
tau_I = Icorr * deltaMD
out_hcacf = np.column_stack((tau_I, normHC))
header_hcacf = 'tau_ps[ps] normalized_HCACF'
np.savetxt('normalized_HCACF.txt', out_hcacf, header=header_hcacf, fmt='%.6e', comments='')
log_step("Saving normalized HCACF to normalized_HCACF.mat (backup)")
savemat('normalized_HCACF.mat', {'Icorr': Icorr, 'normHC': normHC})
log_step("Plotting and saving normalized HCACF figure to PDF")
plt.figure()
plt.plot(tau_I, normHC, 's-')
plt.xlabel('τ (ps)')
plt.ylabel('Normalized HCACF')
plt.grid(True)
plt.tight_layout()
plt.savefig('normalized_HCACF.pdf')
plt.show()
log_step("Writing timing log to timing_log.txt")
with open("timing_log.txt", "w") as f:
f.write("Step Description | Time (s)\n")
for msg, t in step_times:
f.write(f"{msg} | {t:.2f}\n")
return k, normHC
if name == “main”:
log_step(“Calling gk_publishable_plots()”)
gk_publishable_plots()