Possible bug in compute heat/flux for inhomogeneous geometry?

LAMMPS LAMMPS Mailing List Mirror
1

Dear developers and users,

First of all, I followed the example on compute heat/flux page to study a single sheet of graphene without any problem. Then I'm trying to compute the heat flux in a system composed of a graphene sheet on a sio2 substrate in equilibrated constant energy ensemble. However, I found that the heat flux in this case is about several times larger than that case when there is only a single graphene sheet. This is unexpected because if Green-Kubo method is used based on such heat flux value, due to a large autocorrelation function, an unphysically large thermal conductivity value is resulted. I closely followed the example provided on the compute heat/flux page and I did define compute ke/atom etc. using the same group as used by compute heat/flux. I checked the formulation used to compute heat current, which is

J= sum e_i v_i + 0.5 sum (f_ij * ( v_i + v_j)) r_ij

I can correctly plot the phonon DOS, so I think the velocities are fine, and I cannot see why e_i, f_ij and r_ij could be wrong. So I'm thinking if there could be a bug when computing the heat flux in a inhomogeneous system like the one described here. I can post my input and other results if needed. Thanks!

The heat current I obtained for combined system (not divided by volume). For instance J_graphene_y is about 1000 rather than around 100~200
      Step T_total T_graphene T_sio2 J_total_x J_total_y J_graphene_x J_graphene_y J_sio2_x J_sio2_y
      980 306.39639 307.12495 306.26429 617.35114 -1044.5349 74.679789 -929.30782 542.67135 -115.22708
      990 303.30044 293.0246 307.97336 465.86349 -1173.2185 61.494599 -951.29552 404.36889 -221.92303
     1000 301.35906 284.70745 308.81584 243.9928 -1492.084 57.352149 -979.79229 186.64065 -512.29172
     1010 301.05119 288.73513 306.61394 274.21896 -1619.9197 99.706495 -1006.6044 174.51247 -613.31525
     1020 301.68731 297.34578 303.76705 334.78362 -1487.9385 179.08974 -1015.0053 155.69388 -472.93321
     1030 302.48324 300.79699 303.40359 534.94672 -1218.0079 247.80235 -1011.7078 287.14437 -206.3001
     1040 301.49924 298.45163 303.01366 625.72718 -972.88291 264.97314 -1028.5007 360.75403 55.61781
     1050 302.3801 298.5795 304.22397 701.26153 -1075.2168 231.08803 -1080.6749 470.17351 5.4580833
     1060 306.34253 304.64485 307.27023 860.98445 -1058.2901 184.32805 -1149.2712 676.6564 90.981065
     1070 309.61732 306.94997 310.97056 910.03291 -918.57712 161.21415 -1197.1053 748.81877 278.52813
     1080 307.01076 298.23791 311.0294 889.24083 -836.32682 162.17057 -1193.2289 727.07026 356.90206
     1090 303.89042 288.36854 310.85526 613.77726 -1077.6778 159.67712 -1141.8684 454.10014 64.190661
     1100 302.79404 289.6409 308.72351 227.62664 -1229.6045 134.10506 -1086.97 93.521584 -142.63447

The heat current I obtained for single graphene sheet for comparison
      Step T_graphene J_graphene_x J_graphene_y
        0 277.2783 -99.55638 -81.354693
       10 292.16742 24.749244 -185.91829
       20 294.96592 -10.440488 -163.69911
       30 288.76916 -218.61152 -151.89496
       40 288.98167 -105.2722 -155.36314
       50 293.3518 -116.39535 -3.4446361
       60 283.10224 -93.979776 -147.81794
       70 302.96888 -58.062824 -233.84752
       80 284.44289 -180.45572 -124.62181
       90 297.89118 -137.1921 -99.936134
      100 282.53282 -113.83959 -107.37939
      110 294.81692 -148.54612 -82.443137
      120 278.73482 -78.765497 -112.74231

Thanks,
Bo

2

I noticed that in last year's posts German commented on the difference between energy flux and heat flux, I'm wondering if what I'm having here is related to that difference. Thanks!

Bo

3

I imagine what are calling a "bug" has nothing
to do with the compute heat/flux code. It does
what it says it does, which is pretty simple. Whether it
makes sense to use this methodology in your
non-homogenous system or whether you
are interpreting the output correctly, is another question.

Steve

4

I agree with Steve. It's a problem to apply GK method to calculation of
thermal conductivity of interfaces. You probably may use non-equilibrium
method. like Muller-Plathe or direct non-equilibrium.

Best,

German

5

Dear German,

May I ask one question not quite related to Lammps? Regardless of Green-Kubo method and thermal conductivity, I'm just wondering why the fluctuation of energy flux in the graphene + low thermal conductivity sio2 combined system is larger than that of a system with only graphene? Shouldn't the energy flux in low thermal conductivity sio2 be much smaller? At least in NEMD direct method it can be found the energy flux in sio2 is much smaller. Thanks a lot!

Bo