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
clear
next t
jump in.temp
