Different temperature/pressure behavior in small vs. large F-diamane systems using ReaxFF

aibrahim · September 12, 2025, 3:55am

Hello everyone,

I’m working on annealing simulations of fluorinated diamane (F-diamane) using LAMMPS and ReaxFF. I’ve run into some puzzling differences between a smaller and a larger replicated system, and I thought I’d seek feedback on the potential reasons behind this behavior, as I’ve viewed posts here that discuss how stable/“constant” temperature and pressure are attributed to big systems.

My work is as follows: an initial structure of an F-diamane slab with fluorine fully terminating both sides. A small system I created consisted of 384 atoms (256 C + 128 F), and I created a larger system repeating the smaller unit cell in the x- and y-directions, consisting of 3456 atoms (2304 C + 1152 F). Below is the script I use for annealing:

# 1) INITIALIZATION
# ------------------------------
variable T_start equal 10.0       # Starting temperature
variable T_end   equal 773.15     # Annealing target temperature

units real
atom_style charge
boundary p p p

read_data f_diamane_big.data

# ------------------------------
# 2) INTERACTION MODEL
# ------------------------------
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield.reax.FC.txt C F

# ------------------------------
# 3) GROUPS & CHARGES
# ------------------------------
group carbon type 1
group fluorine type 2
group all_atoms union carbon fluorine

fix qeq all_atoms qeq/reaxff 100 0.0 10.0 1.0e-6 reaxff

# ------------------------------
# 4) PRE-ANNEAL RELAXATION
# ------------------------------
fix boxrelax all_atoms box/relax x 0.0 y 0.0 vmax 0.001
minimize 1.0e-6 1.0e-9 1000 10000
unfix boxrelax

velocity all_atoms create ${T_start} 4928459 mom yes rot yes dist gaussian

# ------------------------------
# 5) PHASE 1: Temperature Ramp
# ------------------------------
fix temp_ramp all_atoms nvt temp ${T_start} ${T_end} 50
timestep 0.25
thermo 100
run 10000   # 2.5 ps temperature ramp

unfix temp_ramp

# ------------------------------
# 6) PHASE 2: Hold at T_end
# ------------------------------
fix anneal all_atoms nvt temp ${T_end} ${T_end} 50
fix reax_bonds all_atoms reaxff/bonds 100 bonds.reax

dump traj all_atoms custom 500 dump.lammpstrj id type x y z q
run 90000   # 22.5 ps hold at 773.15 K

# Total simulation time = 25 ps

Below are the respective temperature and pressure plots for the simulation for each system:

Small system clearly showing non-convergence, and so I thought of running the same tests on a bigger system:

Note that the pressure plot was truncated to better display the behavior at the beginning of the run. The rest of the run exhibited smooth trends like the ones shown after the initial fluctuations.

I’m happy for the most part with the temperature trends of the big system, as temperatures are more tightly packed or don’t fluctuate as much. The one concern I have is the initial temperature spike, which I believe is attributed to really high pressure:

   Step          Temp          E_pair         E_mol          TotEng         Press     
        22   10            -502971.03      0             -502868.04     -12120.994    
  -->  100   365.77929     -511393.37      0             -507626.32      48195.76     
       200   413.21375     -514338.02      0             -510082.46     -6067.4625    
       300   160.32654     -513229.38      0             -511578.22     -14845.023    
       400   204.45019     -514380.16      0             -512274.59      2087.0516    
       500   115.7196      -514030.76      0             -512839         6686.7493    
       600   106.61311     -514253.49      0             -513155.51      12875.709    
       700   60.446976     -513999.85      0             -513377.32     -34039.995    
       800   45.34627      -513956.07      0             -513489.07      19513.915    
       900   19.933844     -513774.81      0             -513569.52     -11997.832    
      1000   35.513314     -513879.51      0             -513513.77      20665.207    
      1100   43.236585     -513853.29      0             -513408.01     -41481.65     
      1200   66.191752     -513891.15      0             -513209.46      18514.136    
      1300   83.267452     -513866.3       0             -513008.76     -7292.2423    
      1400   105.54025     -513858.32      0             -512771.39      7358.7768

despite attempting to avoid this issue which I encountered with different simulations in the past by using fix box/relax to rescale the simulation box and resolve any potential in-plane strain.

Onto pressure trends, they seemed to inevitably settle down/converge but I’m curious about the cause of the initial fluctuations and if they’re attributed to my “big” system not being that big, bad practice/scripting on my part, or an alternate reason.

— LAMMPS version: 29 Aug 2024 - Update 1, on Windows

Thank you in advance and apologies for any oversights or misconceptions

aibrahim · September 12, 2025, 4:20am

Forgot to add that every run prints:

WARNING: Triclinic box skew is large. LAMMPS will run inefficiently.

My tilt factor for the respective systems (xy = 10.089 for small system, 30.266 for large system) is ~50% of the x-box length.
Not sure if this large skew is just inefficient, or if it’s affecting the stability of my runs and introducing artificial strain/instability?

Thank you

initialize · September 12, 2025, 11:19am

If you can share your full set of input files (so that others can reproduce it, ideally on a personal computer) then maybe you can get more people to help you. At the current moment I only have some general suggestions (base on the fact that E_pair drops significantly in initial steps):

make sure your minimization is actually converged (check lammps output);
check the trajectory in initial steps to see if there are large displacements;
increase the frequency of charge equilibration. Set Nevery to 1 in fix qeq/reaxff and see if it still happens.

Germain · September 12, 2025, 12:34pm

Hi @aibrahim,

In addition to @initialize comments, note that your potential energy goes down as the simulation runs. Your minimization only lasts for 1000 steps. This is very short. Are you sure your configuration reaches a local minimum before the minimization ends?

Your timestep is also in the high part of reaxff models. You might consider reducing it a bit unless you have some reference stating that this is a correct value for your model.

aibrahim · September 12, 2025, 7:36pm

Thanks for the insights @initialize, attached below are the files I’m using too:

f_diamane_big.data (300.7 KB)
ffield.reax.FC.txt (25.3 KB)

aibrahim · September 12, 2025, 8:06pm

Hi @Germain,

I only now just realized that I’ve been making the mistake of going over/taking for granted the LAMMPS output to check for the status of the minimization, but to answer your question here is the output:

OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
  using 1 OpenMP thread(s) per MPI task
Loaded 1 plugins from C:\Users\ahmed\AppData\Local\LAMMPS 64-bit 22Jul2025 with Python\plugins
Reading data file ...
  triclinic box = (0 0 -10) to (60.533997 52.423979 44.198049) with tilt (30.266999 0 0)
WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:221)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  3456 atoms
  read_data CPU = 0.032 seconds
Reading potential file ffield.reax.FC.txt with DATE: 2013-06-28
WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:300)
2304 atoms in group carbon
1152 atoms in group fluorine
3456 atoms in group all_atoms

Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 16 9 10
  2 neighbor lists, perpetual/occasional/extra = 2 0 0
  (1) pair reaxff, perpetual
      attributes: half, newton off, ghost
      pair build: half/bin/ghost/newtoff
      stencil: full/ghost/bin/3d
      bin: standard
  (2) fix qeq/reaxff, perpetual, copy from (1)
      attributes: half, newton off
      pair build: copy
      stencil: none
      bin: none
Setting up cg style minimization ...
  Unit style    : real
  Current step  : 0
WARNING: Energy due to 3 extra global DOFs will be included in minimizer energies
 (src/min.cpp:222)
Per MPI rank memory allocation (min/avg/max) = 361.2 | 361.2 | 361.2 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press          Volume
         0   0             -497604.85      0             -497604.85      62430.708      171993.88
        22   0             -502971.03      0             -502971.03     -12147.301      178956.67
Loop time of 45.5032 on 1 procs for 22 steps with 3456 atoms

99.1% CPU use with 1 MPI tasks x 1 OpenMP threads

Minimization stats:
  Stopping criterion = linesearch alpha is zero
  Energy initial, next-to-last, final =
     -497604.845836383  -502971.026322075  -502971.026322075
  Force two-norm initial, final = 342542.92 7492.2337
  Force max component initial, final = 242162.95 2837.2084
  Final line search alpha, max atom move = 1.6412657e-16 4.6566129e-13
  Iterations, force evaluations = 22 91

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 45.456     | 45.456     | 45.456     |   0.0 | 99.90
Neigh   | 0.03053    | 0.03053    | 0.03053    |   0.0 |  0.07
Comm    | 0.001656   | 0.001656   | 0.001656   |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0.000108   | 0.000108   | 0.000108   |   0.0 |  0.00
Other   |            | 0.01525    |            |       |  0.03

Nlocal:           3456 ave        3456 max        3456 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:           3894 ave        3894 max        3894 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:         935172 ave      935172 max      935172 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 935172
Ave neighs/atom = 270.59375
Neighbor list builds = 2
Dangerous builds = 0
Setting up Verlet run ...
  Unit style    : real
  Current step  : 22
  Time step     : 0.25
Per MPI rank memory allocation (min/avg/max) = 359 | 359 | 359 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press
        22   10            -502971.03      0             -502868.04     -12120.994
       100   365.77929     -511393.37      0             -507626.32      48195.76
       200   413.21375     -514338.02      0             -510082.46     -6067.4625

So it seems like the minimization didn’t in fact converge to a minimum. I’ll definitely give the minimization more time: minimize 1.0e-6 1.0e-9 10000 100000 and perhaps also tighten the tolerances to 1.0e-8 1.0e-10 for example?

Regarding the timestep, I’ve found timesteps under 0.5 to not be problematic with my ReaxFF work. For this specific simulation, I’ve tried 0.5, 0.25, and 0.1 but the behavior described above persisted for all of them. Not sure if that means the issue isn’t the timestep or if I should go lower.

initialize · September 12, 2025, 10:05pm

So that is the problem. Particularly note the force:
Force two-norm initial, final = 342542.92 7492.2337
The final force of 7492.2337 is pretty large and indicates that your system is far away from a “stable” configuration. And indeed you can notice large shifts of atom positions in the first steps of MD.

The real way of solving the problem should be modify your initial structure so that it is actually close to the real structure. If you just want to reach a minimum then you can try switching to a different minimization algorithm, e.g.

fix boxrelax all_atoms box/relax x 0.0 y 0.0 vmax 0.001
minimize 0 1.0e-9 20 10000
unfix boxrelax
min_style fire
minimize 1.0e-6 1.0e-9 1000 10000
min_style cg
fix boxrelax all_atoms box/relax x 0.0 y 0.0 vmax 0.001
minimize 0 1.0e-9 1000 10000

Now after minimization the force two-norm is 0.99270295 which is reasonably small, and the temperature does not skyrocket in the MD steps.

However, I have no idea if minimization in this way actually leads to a physically meaningful structure (which could only be determined by yourself).

aibrahim · November 4, 2025, 8:21pm

Hi again, thank you for your help! and apologies for the delay. I proceeded by adding the minimization steps recommended by you to my script:

# ------------------------------
# 1) INITIALIZATION
# ------------------------------
variable T_start equal 10.0       # Starting temperature
variable T_end   equal 773.15     # Annealing target temperature

units real
atom_style charge
boundary p p p

read_data f_diamane_big.data

# ------------------------------
# 2) INTERACTION MODEL
# ------------------------------
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield.reax.FC.txt C F

# ------------------------------
# 3) GROUPS & CHARGES
# ------------------------------
group carbon type 1
group fluorine type 2
group all_atoms union carbon fluorine

fix qeq all_atoms qeq/reaxff 100 0.0 10.0 1.0e-6 reaxff

# ------------------------------
# 4) PRE-ANNEAL RELAXATION
# ------------------------------
fix boxrelax all_atoms box/relax x 0.0 y 0.0 vmax 0.001

# 1) Quick relaxed-energy step to let box change slightly (no force tol)
#    (minimize parameters: etol ftol maxiter maxeval)
minimize 0 1.0e-9 20 10000

unfix boxrelax


min_style fire
minimize 1.0e-6 1.0e-9 1000 10000


min_style cg
fix boxrelax all_atoms box/relax x 0.0 y 0.0 vmax 0.001
minimize 0 1.0e-9 1000 10000
unfix boxrelax

# Printing final minimization diagnostics in the log (they will appear automatically).
# more informative thermo style so you can see box dims and energies during the run:
thermo_style custom step temp pe etotal press vol lx ly lz

velocity all_atoms create ${T_start} 4928459 mom yes rot yes dist gaussian

# Start recording trajectories and bonds before temperature ramp
fix reax_bonds all_atoms reaxff/bonds 100 bonds.reax
dump traj all_atoms custom 500 dump.lammpstrj id type x y z q

#short nve energy conservation test: Watch etotal/pe drift. Energy should be reasonably conserved (no explosive drift). If large drift, structure still has problems
fix test_nve all nve
thermo 100
run 8000
unfix test_nve


# ------------------------------
# 5) PHASE 1: Temperature Ramp
# ------------------------------
fix temp_ramp all_atoms nvt temp ${T_start} ${T_end} 200
timestep 0.25
thermo 100
run 40000   # 2.5 ps temperature ramp

unfix temp_ramp

# ------------------------------
# 6) PHASE 2: Hold at T_end
# ------------------------------
fix anneal all_atoms nvt temp ${T_end} ${T_end} 250

run 160000   # 22.5 ps hold at 773.15 K

# Total simulation time = 25 ps

I monitored Etotal for both the NVE run and the NVT run for phase 2 and it stayed consistent during each run, so the structure didn’t have problems as the energy didn’t drift

I then plotted the results for temperature and pressure:

The temperature is very good now, but pressure seems to still face issues with convergence, perhaps even worse than before. I’m unsure if this is due to the small size of my system or another potential cause? It looks like it started off relatively stable and then started to fluctuate much more during heating and afterwards too. I had a look at the log file and the violent fluctuation indeed starts with the temperature ramp and persists after that. The structure was also sound and remained intact at the end of the simulation and throughout it.

Are my pressure fluctuations reasonable for NVT on a slab? or should I run (fix npt) for example or switch to p p f and only barostat in-plane maybe?

Thank you again

initialize · November 5, 2025, 5:21am

There are several things you need to know:

the pressure value given by thermo output is not very meaningful if your system is not physically 3D (e.g. slabs);
the equilibrium box size at finite temperature (i.e. in MD) is generally different from the one at “zero temperature” (i.e. minimization);
the pressure value always has large fluctuations due to the limited system size in typical MD simulations (there are already lots of discussions on this topic in this forum), and this is usually not a serious problem by itself. Finite size effects does exist in MD but you should not determine its severity by the pressure fluctuations.

So due to the point 1 and 3, you should not need to worry about the pressure fluctuation itself.

However, to get a physically good simulation model, you probably need to deal with the point 2 above. This is honestly out of my knowledge (I don’t think npt is good here since the pressure output itself is not meaningful in this context). Maybe someone else with experience in similar systems can help you.

srtee · November 9, 2025, 1:44am

In your very first post you had pressure fluctuations ranging about 30,000 atm.

In your latest post the fluctuations range about 8,000 atm and include instantaneously negative values, which is typical for a solid at equilibrium at room temperature and pressure. As @initialize said, there are extensive discussions on pressure fluctuations in small systems on this forum. Not only are they normal and expected, they are related to the compressibility of your material, so you can get a rough order of magnitude of the “computational compressibility” and compare it to known values.

akohlmey · November 9, 2025, 2:23am

You should not look at the total pressure but at the components in the periodic dimensions.
For non-periodic dimensions what LAMMPS uses to compute pressure is the length of the simulation box, but this is for convenience reasons. In practice the pressure component for a non-periodic boundary is by definition 0.0 since what a LAMMPS non-periodic dimension really represents is an infinitely large vacuum around your slab.

To make it easier to see trends despite the large fluctuations (e.g. to determine if you get closer to equilibrium) people usually use fix ave/time to compute a running average over a sufficiently large window.

P.S.: BTW the y axis of your pressure vs. timestep plot is mislabeled. Also, shouldn’t the x axis be labeled “Time” and not “Timestep” or if it is the number of timesteps, then there should be no unit?