comparison of lammps and reaxff stand alone energies during minimization

Hi,

I am minimizing a Ni - O system (4000 Ni atoms with O in a tetrahedral or octahedral interstitial) using a reaxff potential. The total potential energy and other partial contributions of step 0 are listed below, both from lammps and reaxff stand alone code (4000 Ni atoms with O in a octahedral site).


EPOT
EBOND
EATM
EVAL
EVDW
ECOUL
ECHARGE
ReaxFF -405982.0919
-1250487.34
684519.52
50.02
159938.38
-2.66
-0.02














Lammps -405982.3542
-1250487.336
684519.5235
49.7545
159938.38
-2.6559
-0.0203

The final energies upon minimization are as follows:


EPOT
EBOND
EATM
EVAL
EVDW
ECOUL
ECHARGE
ReaxFF -405990.9952
-1250484.56
684517.62
48.45
159930.13
-2.54
-0.1














Lammps -406003.5147
-1250485.42
684518.8921
44.9878
159920.619
-2.3819
-0.2122

In the final energies, the maximum difference is coming from eval and echarge. The energies for pure nickel system are equal from simulation run with lammps as well as reaxff stand alone code. The above mentioned energy difference does translate in to something significant as the values from lammps suggest that putting O in either octahedral or tetrahedral interstitial costs almost equal energy where as the energies from reaxff stand alone code results in octahedral interstitial O costing less energy. The result from reaxff stand alone code is what has been predicted earlier from DFT calculations. Can anyone give more insight on this error?

regards,
Karthik Guda, PhD
Penn State University

This is likely due to the version of LAMMPS, control parameters and minimization settings. Was control file used? If yes, what were the control parameters? How was minimization set up? Was it minimized with fix box/relax? How was the minimization stopped? Check these and maybe you can track the origin of inconsistency.

Ray

The version I am using is 28 Mar 2012. I did not use any control file, no box/relax has been used and the minimization stopped by saying “linesearch alpha is zero”. Below is my input file.

units real
atom_style charge # corase-grain liquids, solids, metal
boundary p p p # boundary conditions

read_data Ni_octa1.data

neighbor 2 bin
neigh_modify every 10 delay 0 check no
pair_style reax/c NULL checkqeq yes
pair_coeff * * ffield.reax 6 3
compute gv all pair reax/c
variable eb equal c_gv[1]
variable ea equal c_gv[2]
variable ev equal c_gv[5]
variable ep equal c_gv[12]
variable eqeq equal c_gv[14]
variable evdw equal c_gv[11]
thermo_style custom etotal vol press pe ke temp ecoul evdwl elong etail pxx pyy pzz pxy pxz pyz lx ly lz xy xz yz v_eb v_ea v_evdw v_eqeq v_ev v_ep fnorm fmax
thermo_modify line multi flush yes
thermo 1
dump 1 all custom 10 ini_test3.dump id type q xu yu zu
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
minimize 0.0 1.0e-8 1000 100000

The relaxed structure from the above simulation has a fnorm value of 0.017 where as the RMSG (fnorm) of forces of the relaxed structure from reaxff stand alone code is about 0.48. I reran lammps minimization run with ftol = 0.5, and the new relaxed structure has an fnorm value of 0.44. Below is the comparison between the energies of the two strutures (lammps fnorm 0.44 and reaxff 0.48):


EPOT
EBOND
EATM
EVAL
EVDW
ECOUL
ECHARGE
ReaxFF -405990.9952
-1250484.56
684517.62
48.45
159930.13
-2.54
-0.1














Lammps -406003.5099
-1250485.365
684518.8496
44.9965
159920.603
-2.3832
-0.2109

The discrepancy still exists and the difference in energies translate into disagreement in insertion energy of O in octahedral and tetrahedral positions

regards,
Karthik Guda

Since you are using a March version, one source of inconsistency arose from valence bond order cutoff.

Subroutine srtang, reac.f, stand alone code:
if (bo(i1)*bo(i2).lt.0.00001) goto 50.

Function Valence_Angles, reaxc_forces.cpp, LAMMPS user-reaxc (March version):
…bo_ij->BO * bo_jk->BO > 0.001) ) …

These “0.00001” and “0.001” are the cause of your different EVal. Please update to the current distribution, with which the 0.001 is changed to a default value of 0.00001.

Ray

Hi,

I am now using the current distribution of lammps (22 Aug 2012). The difference in energies still exists in the final structures and the difference mainly comes from Eval,Echarge and Evdw.

Below are the energies:


EPOT
EBOND
EATM
EVAL
EVDW
ECOUL
ECHARGE
ReaxFF -405990.9952
-1250484.56
684517.62
48.45
159930.13
-2.54
-0.1














Lammps -406003.5098
-1250485.364
684518.8494
44.9967
159920.6025
-2.3818
-0.2124

regards,
Karthik