eqeq in reax/c

Dear all,

I performed some checks on eqeq (charge equilibration) energy of reax/c. I have plotted eqeq vs system size at different temperatures. Plot attached. The result I get seems a bit suspicious to me. eqeq converges for about 30000 atoms which seems too big ‘requirement’. Manual says qeq does not work for the simulation box size < 10Å in any dimension, but the cell corresponding to 30000 atoms is much bigger than 10Å. Fluctuation is not ‘huge’ but still important and increases with temperature.

The check was done for Al2O3 with default control file. I’m using triclinic box. Could there be any problem related to it?

Script below:

log temp.log

variable t index 1 2 3 4 5 6 7 8 9 10 11

units real
atom_style charge
boundary p p p
dimension 3
atom_modify map array sort 0 2.0

read_data Al2O3.data #30 atom unit cell, triclinic

replicate $t $t $t

pair_style reax/c NULL
pair_coeff * * ffield.reax.Al_Al0_AlN 7 3

fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

neighbor 2.5 bin
neigh_modify every 1 delay 0 check yes

velocity all create 300 24577 dist gaussian
fix 2 all npt temp 300 300 100.0 iso 0.0 0.0 1000.0
#fix 3 all nve

print '-- t----' compute reax all pair reax/c variable ew equal c_reax[11] variable ep equal c_reax[12]/({t}^3)
variable eqeq equal c_reax[14]/({t}^3) variable coulsum equal (v_ep+v_eqeq)/({t}^3)
variable evdwl equal evdwl/(${t}^3)

thermo 1

thermo_style custom step etotal vol v_evdwl ecoul v_ew v_ep v_coulsum v_eqeq

fix out all print 1 ‘{t} {eqeq} {evdwl} {coulsum}’ append temp300.out title “”

run 1

next t
jump in.temp


Dear Manana,

I think you are seeing this behavior because of fix npt. Due to
differences in system sizes, it is normal for fix npt to give slightly
different equilibrium sizes. Size of the cell affects equilibrium
bond lengths then affects equilibrium charges then the eqeq energies.
This may be why there are more dramatic changes for higher temperature

To make sure fix qeq and pair_reax/c are working correctly for your
test, it would be better to start with 0K and nve. Verify first eqeq
does not change with system size under this condition. Then introduce
temperature, and then fix npt. If they do change, then it is just
temperature or system size effects.