NEMD temperature spiking to 1000K at grain boundary

Ankit213910 · February 23, 2026, 3:29am

Dear LAMMPS experts,
I am using heat add/subtract method to get a temperature profile across LiAlO2 using ReaxFF potentials for single crystal and bicrystal with grain boundary. Here, 0.0005 eV heat is added at the two ends/step and -0.001 ev/step heat is subtracted from the center. The shape is long in x direction as shown below:
LAMMPS data file written by OVITO Basic 3.5.4

   38400  atoms
       3  atom types

  0.000000000000     418.410640000000  xlo xhi
  0.000000000000      31.380798000000  ylo yhi
  0.000000000000      31.633710000000  zlo zhi

Masses

        1   6.94000000              # Li
        2   26.98153860             # Al
        3   15.99900000             # O

The code is below

#---------initialization---------
echo screen
units metal
dimension 3
boundary p p p
atom_style charge

read_data single_p1_x

variable t equal 300
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s

##############################

Force Field Parameters

##############################

pair_style reaxff lmp_control safezone 2.0 mincap 100
pair_coeff * * reaxff Li Al O

timestep 0.001

1st equilibration run

velocity all create 300 87287

region hotl block 0 10.41 INF INF INF INF
region hotr block 408 418.41 INF INF INF INF
region cold block 200 221 INF INF INF INF
compute Thot all temp/region hotl
compute Tcold all temp/region cold

fix 1 all nvt temp $t $t 0.001
fix 3 all qeq/reax 1 0.0 10.0 1e-6 param.qeq

thermo 100
run 10000

unfix 1
unfix 3

velocity all scale 300

2nd equilibration

reset_timestep 0

fix 1 all nve
fix 3 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
fix hotl all heat 1 0.0005 region hotl
fix hotr all heat 1 0.0005 region hotr
fix cold all heat 1 -0.001 region cold

thermo_style custom step temp c_Thot c_Tcold
thermo 1000
run 10000

temperature profile

compute layers all chunk/atom bin/1d x lower 25 units box

fix 2 all ave/chunk 10 100 1000 layers temp file profile.heat
thermo 100

run 100000

For the single crystal in x direction, I get a pretty linear profile shown below.

However, as soon as I put 2 grain boundaries at the 1/4 length and 3/4 length (to maintain symmetry), the same exact input file gives a very unexpected temperature profile with temperature spikes at the location of grain boundaries as shown below

May I ask if I am doing something wrong here? I tried different grain boundary types like sigma 5 and sigma 41 and both are showing similar behavior. Is there any suggestions/corrections I can try please?

I highly appreciate any inputs from you.

Sincerely,

Ankit

Germain · February 23, 2026, 3:03pm

Hi @Ankit213910,

First it appears that you know how to triple quote your text to past it as plain text in your messages. Please think about doing so with input files since this makes them easier to read. Symbols can be interpreted in the Markdown rendering and pollute your message.

As for your issue there is very little that people could tell you without more information on your grain boundaries systems. Have you only looked if the grain boundaries where stable or the final state of your simulations? The experimental melting temperature of LiAlO2 is close 1900K. Is the model stable between 800K and 1400K? It is possible you have melted parts of your system with void in between. Note that if you have diffusion in your system due to heat flux of stress, you should use dynamic groups if you want to keep your heat sink and sources at the same locations.

Ankit213910 · February 24, 2026, 5:04am

in.stablegb (2.1 KB)

last (2.0 MB)

final (2.6 MB)

log_output.out (31.0 KB)

Dear Germain,

Thank you for taking the time to look through the problem.

I looked back at the simulation and realized that the system does have relaxation problem.

The temperature is constantly tending to shoot up to high values of 2000K.

For that I tried to relax in NVT ensemble using very small damping of 0.2. But as soon as I used a NPT even with 0.5 ps damping the temperature shot up again.

All files are attached

input file- in.stablegb

starting structure- final

last output structure- last

log output- log_output.out

Also, if tried to skip the NPT ensemble relaxation and moved straight to the heat add/subtraction in NVE ensemble but the system temperature shoots again to high values.

Do you have any suggestions to treat the GB to stabilize it?

Sincerely,

Ankit

Germain · February 24, 2026, 6:58pm

You’re using isotropic NPT on a solid system that is highly anisotropic which a plane dislocation in its center.

This is doomed to fail.

Try using the aniso keyword of fix NPT. Please read the documentation of the NPT command to understand the issue.

Ankit213910 · March 20, 2026, 1:13am

Dear Germain,

Thank you for your thoughtful suggestions.

after your inputs, I tired the following methods.

Construction of structure file- I made 3 crystals (length along x axis) and attached them along x axis. Crystal 1 is 104 Angstrom X 31 A X 31 A, crystal 2 is 208 A X 31 A X 31 A but is rotated by just 1 degree around the x axis using atomsk, crystal 3 is same size as crystal 1 with dimensions 104 Angstrom X 31 A X 31 A. All 3 crystals are individually charge neutral with no missing atoms or bonds.

All 3 crystals are then attached together (attached file is called combine1.6) and I intentionally put a 1.6 Angstrom gap in between crystals as shown below in the image. I put the gap so that atoms dont overlap with each other. Most bond lengths in LiAlO2 are 1.6 A or larger so I put a gap of 1.6 A.

Then I minimize and relax the structure in NVT ensemble. After that I do the NEMD by supplying heat at left and right ends and extracting heat from the center. The location of the grain boundary still shows temperature spike as shown in the right plot below. Previously without the 1.6 angstrom gap, I was observing temperature spikes of up to 1000 degrees which has currently come down to just 50 degree but is still hotter than the heat source bin.

The single crystal gives a nice profile matching with the experimental thermal conductivity.

To avoid the temperature spike, I tried to relax the system again in NPT ensemble as shown in the input file and using this snippet

fix npt_check all npt temp 300 300 0.5 aniso 0 0 1.0

fix 3 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
thermo 200
thermo_style custom step temp press pxx pyy pzz pe etotal vol density
run 10000
unfix npt_check
unfix 3

but the system with grain boundary crashes immediately saying “bond check failed” and pressure extremely high positive in all 3 directions.

Do you recommend any other way to construct the grain boundary to avoid this problem? Any other suggestion to obtain the kaptiza resistance at the grain boundary will be helpful.

The structure file and input files are attached.

looking forward to your response.

combine1.6 (1.8 MB)

input.txt (2.8 KB)

Germain · March 20, 2026, 1:54pm

If this is the “bondchk” error this is a reaxff issue.

Please have a look at the memory management issues detailed in the manual. The original version if known to be very sensitive to those issues. You can either try to increase memory management parameters of try using the Kokkos version which is more robust on this regard.

akohlmey · March 21, 2026, 3:30am

That usually means that your geometry is not correct.

Ankit213910 · March 22, 2026, 2:17am

Thank you Germain and Axel for your thoughts.
May I ask how I should correctly make the grain boundary?

Currently I am doing this:

I made 3 crystals (length along x axis) and attached them along x axis. Crystal 1 is 104 Angstrom X 31 A X 31 A, crystal 2 is 208 A X 31 A X 31 A but is rotated by just 1 degree around the x axis using atomsk, crystal 3 is same size as crystal 1 with dimensions 104 Angstrom X 31 A X 31 A. All 3 crystals are individually charge neutral with no missing atoms or bonds.

All 3 crystals are then attached together (attached file is called combine1.6) and I intentionally put a 1.6 Angstrom gap in between crystals as shown below in the image. I put the gap so that atoms dont overlap with each other. Most bond lengths in LiAlO2 are 1.6 A or larger so I put a gap of 1.6 A.

Would you recommend any other method please?

I will start looking into the memory issues of reaxff, once I am able to correct the grain boundary construction.

I highly appreciate you giving this a thought.

Sincerely,

Ankit

initialize · March 22, 2026, 2:38am

Apparently your input structure has some atoms abnormally close to each other: