Hi all,
Context:
I’m conducting MD simulations of the fluorination of bilayer graphene using the ReaxFF force field. My goal is to study the system at various temperatures (700K, 500K, 300K…) to observe C-F bonding. The graphene bilayer was built in VMD, and periodicity was manually verified and visualized in OVITO. It appeared structurally sound when I replicated the unit cell in both x and y directions.
The Problem
When I begin the simulation, the system is first minimized to 0 K. I then apply a thermostat to bring it to 700 K. While the temperature ramps correctly initially, I consistently observe a sharp temperature spike to ~830 K shortly after velocity initialization — even when fluorine atoms are not added, and even at very small timesteps (e.g., 0.1 fs). The spike later relaxes and the system stabilizes, but I’m trying to resolve this transient behavior. Below are some output lines showing this behavior (sorry for the crowded view):
Step Temp E_pair E_mol TotEng Press
58 700 -1190540.1 0 -1175752.6 92958.971
100 522.85601 -1186389.9 0 -1175344.6 67898.317
150 732.60283 -1190359.5 0 -1174883.3 29355.424
200 803.10733 -1191565.6 0 -1174599.9 19041.264
250 777.73721 -1191451 0 -1175021.3 16333.987
300 816.93246 -1193152.1 0 -1175894.4 15345.358
350 829.22071 -1194330.6 0 -1176813.3 13201.88
400 796.83794 -1194624 0 -1177790.8 10883.634
450 770.03978 -1194958.6 0 -1178691.5 8697.1269
500 768.01306 -1195783.1 0 -1179558.8 8110.3273
550 766.29399 -1196605.2 0 -1180417.3 7558.6894
600 753.48731 -1197179.8 0 -1181262.4 6456.0001
650 743.74551 -1197760.6 0 -1182048.9 5839.4056
700 740.67396 -1198429 0 -1182782.2 5036.8335
750 728.1488 -1198858.7 0 -1183476.6 5000.4342
800 720.16754 -1199329.6 0 -1184116.1 5154.1862
850 715.07721 -1199800 0 -1184694 5559.0157
900 701.42814 -1200011.8 0 -1185194.2 5787.4037
950 695.25765 -1200282.6 0 -1185595.3 5893.9998
What I’ve Tried So Far
- Tested a wide range of timesteps: 0.1 fs, 0.25 fs, 0.5 fs, and 1 fs. The spike still appears, though slightly smaller at lower timesteps.
- Ran the simulation without adding fluorine atoms to try and isolate the source — the spike persists, so the F atoms aren’t the root cause.
- Verified the box bounds and structure replication manually and visually (VMD and OVITO). The geometry looks clean and periodic.
- Checked that my ReaxFF file is suited for my system (link), which resulted in a way softer spike instead of ~1280 K previously.
Suspected Causes / Next Steps I’m Considering
Based on a 2015 paper by Jensen et al. (J. Comput. Chem., 36, 1587–1596), which explores the effect of timestep and thermostats on ReaxFF behavior in carbon systems, I’m thinking this temperature spike may be related to:
- Lack of equilibration via Langevin damping or similar methods to smooth the transition from minimization to NVT
- Localized stress or energy pockets forming in pristine carbon structures under thermal activation (they report similar ReaxFF effects in graphene and diamond)
- Underlying graphene lattice issues that I should perhaps look into, maybe run tests using a simpler force field (e.g. AIREBO) to try and diagnose the problem?
So now I’m trying to decide my next step, if I should introduce a short equilibration phase before the NVT run? or are the other points listed worth looking into again?
Here is my script for reference:
# 1) Initialization
variable T equal 700 # Starting temperature
units real
atom_style charge
boundary p p f # periodic in x and y, fixed in z
# Read graphene structure
read_data graphene.data
# Use ReaxFF
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield.reax.FC.txt C F
# 2) Addition of Fluorine atoms to system
# Define groups for Gr and F
group graphene type 1
group gas type 2
# Define regions for F atoms, above and below the graphene sheets
region above_graphene block 3 93 3 93 15 25 units box
region below_graphene block 3 93 3 93 -25 -15 units box
# Create fluorine atoms in the above and below regions
create_atoms 2 random 200 12345 above_graphene overlap 1.5
create_atoms 2 random 200 34134 below_graphene overlap 1.5
# Charge equilibration
fix qeq all qeq/reaxff 10 0.0 10.0 1.0e-6 reaxff maxiter 400
# Prevent atoms from leaving box along the z axis
fix wall all wall/reflect zlo EDGE zhi EDGE
compute tgraphene graphene temp
compute tgas gas temp
# Minimization and initialization of temperature
minimize 1.0e-6 1.0e-9 1000 10000
velocity all create ${T} 4928459 mom yes rot yes dist gaussian
# Output simulation data
dump mydmp all custom 500 dump.lammpstrj id type x y z q
thermo 50
# Full-system NVT simulation: apply one thermostat to the entire system
fix nvt_sim all nvt temp ${T} ${T} 50
# Monitor bond formation (all bonds)
fix reax_bonds all reaxff/bonds 50 bonds.reax
timestep 0.5
run 40000
— LAMMPS version: 29 Aug 2024 - Update 1, on Windows
Thank you in advance! and apologies for any oversights
Some things with eight legs are not spiders, and some simulations with smoothly changing temperatures may not be valid. It is up to you to know what you expect to see in your simulation and understand the results.