Energy difference at hot/cold region in NEMD simulation

Dear all,

I’m trying to estimate thermal conductivity of 2D layered MOF by setting hot/cold region in the simulation cell. When temperature gradient is applied in the cross-plane direction, there was a large difference between the absolute values of in/out potential energy at the two regions. The increasing rates were also different and the extracted energy at the cold region is always larger.
This problem is not observed in the in-plane direction, and I guess the limited energy transfer in the cross-plane direction could be the reason as only LJ potential is set between layers. However, I could not solve the problem while I changed several parameters like cut off or pppm/disp accuracy.
I’d like to ask if anyone encountered same situation and solve the problem.

Below is my simulation setting.

units real
atom_style full
boundary p p p
timestep 1

pair_style lj/long/coul/long long off 15
pair_modify tail yes mix arithmetic
dielectric 1.0
special_bonds lj/coul 0.0 0.0 1.0
neighbor 3.0 bin
neigh_modify every 1 delay 0 check yes
kspace_style pppm/disp 1.0e-7
kspace_modify disp/auto yes

fix fi_NVE all nve
fix fi_heat MOF_heat langevin 320 320 100 \{random\_seed\} tally yes fix fi\_cool MOF\_cool langevin 280 280 100 {random_seed} tally yes
fix_modify fi_heat energy yes
fix_modify fi_cool energy yes

Thank you,
Shingi Yamaguchi