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.
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 why the structure broke down before the structure relaxation in the NPT ensemble. Is it just because the Al2O3 is too long or the system is not orthogonal or something else?
Would anyone please give me some suggestions on this issue?
Thank you in advance.
