Structure breaks down in the simulation with LAMMPS

LAMMPS LAMMPS General Discussion
1

I used the following input file to simulate the thermal conductance of the Al2O3/Bilayer Graphene/Al2O3 structure along the z axis.

variable         kb equal 8.617333262e-5
variable         t equal 300.0
variable         th equal 400.0
variable         tc equal 200.0

units            metal
dimension        3
newton           on
boundary         p p p
atom_style       atomic
neighbor         0.3 bin
timestep         0.0005

read_data        atomic_structure_atomic

pair_style       hybrid vashishta lj/cut 12.0 airebo 3.0
pair_coeff       * * vashishta ./Al2O3.vashishta O Al NULL NULL
pair_coeff       * * airebo ./CH.airebo NULL NULL C C
pair_coeff       1 3 lj/cut 0.006199133907248657 3.692748434432002
pair_coeff       1 4 lj/cut 0.006199133907248657 3.692748434432002
pair_coeff       2 3 lj/cut 0.042256836241015484 3.1764162195782846
pair_coeff       2 4 lj/cut 0.042256836241015484 3.1764162195782846

compute          ke1 all ke/atom
variable         tem atom c_ke1/1.5/${kb}

thermo           1000
thermo_style     custom step temp pxx pyy pzz press vol
velocity         all create $t 898758 mom yes rot yes dist gaussian

min_style        cg
minimize         1.0e-8 1.0e-8 1000000 1000000

#1st euqilibrium in NPT ensamble
thermo           2000
thermo_style     custom step temp pe ke etotal vol pxx pyy pzz press
fix              1 all npt temp $t $t 0.05 aniso 0.0 0.0 0.5
dump             xyz_npt all xyz 1000 dump_npt.xyz
dump_modify      xyz_npt element O Al C1 C2
run              500000
unfix            1
undump           xyz_npt
reset_timestep   0
write_data       atomic_structure_atomic_npt.dat

#2nd equilibrium in NVT ensamble
thermo           2000
thermo_style     custom temp pe ke etotal vol pxx pyy pzz press
fix              2 all nvt temp $t $t 0.05
dump             xyz_nvt all xyz 1000 dump_nvt.xyz
dump_modify      xyz_nvt element O Al C1 C2
run              500000
unfix            2
undump           xyz_nvt
reset_timestep   0

#Dividing the zone
variable         zbl equal zlo
variable         zbh equal zhi
variable         nla equal 100
variable         sca equal 1/${nla}
variable         len equal ${zbh}-${zbl}
variable         dz equal ${len}/${nla}
variable         zll equal ${zbl}+${dz}*2
variable         zlh equal ${zbl}+${dz}*6
variable         zrl equal ${zbh}-${dz}*6
variable         zrh equal ${zbh}-${dz}*2

region           reg_whole block INF INF INF INF ${zbl} ${zbh} units box
region           lef_fixed block INF INF INF INF INF ${zll} units box
region           hot block INF INF INF INF ${zll} ${zlh} units box
region           cold block INF INF INF INF ${zrl} ${zrh} units box
region           rig_fixed block INF INF INF INF ${zrh} INF units box
region           fixed_both union 2 lef_fixed rig_fixed

group            fixed region fixed_both
group            hot region hot
group            cold region cold
group            move subtract all fixed

reset_timestep   0
change_box       all boundary p p f

velocity         fixed set 0 0 0
fix              fixed fixed setforce 0 0 0

fix              3 hot langevin ${th} ${th} 0.05 14565 tally yes
fix              4 cold langevin ${tc} ${tc} 0.05 16576 tally yes
fix              5 all nve

variable         time equal step
variable         el equal f_3
variable         er equal f_4

compute          6 all chunk/atom bin/1d z lower ${sca} units reduced
fix              7 all ave/chunk 1 10000 10000 6 v_tem file profile.heat
thermo           2000
thermo_style     custom step temp ke pe etotal press vol
dump             xyz_nve all xyz 1000 dump_nve.xyz
dump_modify      xyz_nve element O Al C1 C2
run              500000
unfix            7
undump           xyz_nve
reset_timestep   0

fix              8 all ave/chunk 1 10000 10000 6 v_tem file time_average_temperature.dat
fix              eou all print 2000 "${time} ${el} ${er}" file energy.dat title "time e1 e2" screen no

thermo           20000
thermo_style     custom step temp ke pe etotal press vol
dump             fin all xyz 20000 fin.xyz
dump_modify      fin element O Al C1 C2
run              500000

The length of the Al2O3 on each side of the graphene is 60 nanometers to include the phonon scattering effect. The z axis is perpendicular to both x and y axes and the angle between x and y axes is 120 degree.

After I launched the calculation, I found that the Al2O3 structure broke down into several pieces along the z axis before the structure relaxation in the NPT ensemble, shown in the figure below.

image
image1920×1200 45.1 KB

This is minimization part before the NPT ensemble in the log.lammps file and I see no problem with them.

......
......
......
  (6) neighbor class addition, perpetual, trim from (5)
      attributes: full, newton on, cut 6.3
      pair build: trim
      stencil: none
      bin: none
Per MPI rank memory allocation (min/avg/max) = 10.31 | 10.63 | 12.28 Mbytes
   Step          Temp           Pxx            Pyy            Pzz           Press          Volume    
         0   300           -1076017.5     -1076023.9     -730882.13     -960974.51      615437.33    
      1000   300           -688351.21     -196985.71     -289204.07     -391513.66      615437.33    
      2000   300           -579771.31     -172489.76     -190150.25     -314137.11      615437.33    
      3000   300           -537753.47     -161814.31     -113940.49     -271169.42      615437.33    
      4000   300           -501933.82     -150357.51     -69332.576     -240541.3       615437.33    
      5000   300           -482425.4      -141388.12     -45320.842     -223044.79      615437.33    
      6000   300           -474994.54     -139769.04     -31335.682     -215366.42      615437.33    
      7000   300           -466052.17     -136965.52     -15117.402     -206045.03      615437.33    
      8000   300           -462875.62     -137443.51     -8425.4562     -202914.86      615437.33    
      9000   300           -460286.1      -137468.76     -5082.0359     -200945.63      615437.33    
      9705   300           -458378.44     -138243.2      -3363.5665     -199995.07      615437.33    
Loop time of 1113.12 on 96 procs for 9705 steps with 75400 atoms

99.4% CPU use with 96 MPI tasks x no OpenMP threads

Minimization stats:
  Stopping criterion = energy tolerance
  Energy initial, next-to-last, final = 
      -757842.22766864  -826743.919712325  -826743.927871976
  Force two-norm initial, final = 7411.7974 1.2164204
  Force max component initial, final = 902.49618 0.21573228
  Final line search alpha, max atom move = 0.24797873 0.053497017
  Iterations, force evaluations = 9705 18380

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 4.9564     | 221.08     | 334.1      | 426.3 | 19.86
Neigh   | 75.142     | 499.64     | 703.98     | 587.1 | 44.89
Comm    | 8.8415     | 161.43     | 962.61     |1534.9 | 14.50
Output  | 0.0035638  | 0.0036056  | 0.0036534  |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 231        |            |       | 20.75

Nlocal:        785.417 ave        1050 max           0 min
Histogram: 4 3 3 2 3 1 4 12 45 19
Nghost:        10156.7 ave       12907 max        1716 min
Histogram: 3 0 2 4 5 7 5 7 48 15
Neighs:              0 ave           0 max           0 min
Histogram: 96 0 0 0 0 0 0 0 0 0

Total # of neighbors = 0
Ave neighs/atom = 0
Neighbor list builds = 3452
Dangerous builds = 9

#1st euqilibrium in NPT ensamble
......
......
......

I wonder the reason why the structure broke down before the structure relaxation in the NPT ensemble. Is it just because the Al2O3 is too long or the potential issue?

Would anyone please give me some suggestions on this issue?

If the error comes from the potential, would anyone please recommend some proper potentials to describe Al2O3 for the thermal conductance calculation?

Thank you in advance.

2

Maybe try a simpler bulk system first and try reproducing equilibrium properties of the materials in the literature to verify your potential.

3

Thank you for your suggestions.

I will do the test.

Can I ask another question regarding this calculation? I set the higher temperature (400 Kelvin) on one side of Al2O3 and lower temperature (200 Kelvin) on the other side of Al2O3. According to this setting, is it possible to obtain a plateau figure for the time averaged temperature v.s. the distance along the thermal transport direction as follows?

Untitled
Untitled792×612 4.98 KB

Do you think such time averaged temperature v.s. the distance along the thermal transport direction relation makes sense?

Thank you again.

4

I searched the literature (Janak Tiwari and Tianli Feng, Accurate prediction of thermal conductivity of Al2O3 at ultrahigh temperatures, Phys. Rev. B 109, 075201, 2024), and found that the phonon mean free path of Al2O3 is 50 nanometers.

In my simulation, I set the length of Al2O3 as 60 nano meters on each side of graphene layer. Do you think that the length of Al2O3 in my simulation is long enough to include the phonon scattering effect and avoid the plateau in the temperature v.s. distance figure?

I would really appreciate that you could provide any suggestions/comment.

5

These questions are mostly related to science and not specific to LAMMPS software, and therefore not likely to be answered on the LAMMPS forum because we aren’t experts on this particular science problem. I would suggest you talk to your local advisor/mentor on the science part.