forces between two atoms using reax

hello!

I am trying to get the pairwise energy between two atoms, 1 carbon and another nitrogen kept 12 Å from each other initially. I displace the atom but I keep getting energy as 0. I want only pairwise energy so I haven’t provided any velocity/temperature. Can someone please advice me? Thank you. Here is my input script :

#Deposit file

clear
echo both

#-------------------------------INITIALISE SYSTEM-----------------------------
units real
atom_style charge
boundary p p p
timestep 0.25
#---------------------------DEFINE GROUPS AND AND VARIABLES-------------------
read_data input_equi.dat

group n2 type 2
group carbon type 1

#-------------------------------SET POTENTIAL-------------------------------
pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.reax C N
neigh_modify delay 1 check yes

velocity carbon set NULL NULL NULL units box
velocity n2 set NULL NULL NULL units box
compute reax all pair reax/c
variable ew equal c_reax[11]

thermo_style custom epair evdwl elong v_ew

#-------------------------------LOOP DEFINITIONS---------------------------
variable imax equal 1000
variable i loop ${imax}
label loopstart
displace_atoms n2 move 0 0 -0.05 units box

run 1
next i
jump SELF loopstart

Here is my data file:

#datafile

2 atoms
2 atom types

0 50 xlo xhi
0 50 ylo yhi
0 50 zlo zhi

Masses

1 12.0107
2 14.0067

Atoms

1 1 0 25 25 25
2 2 0 25 25 37

hello!

I am trying to get the pairwise energy between two atoms, 1 carbon and
another nitrogen kept 12 Å from each other initially. I displace the atom
but I keep getting energy as 0. I want only pairwise energy so I haven't
provided any velocity/temperature. Can someone please advice me? Thank you.
Here is my input script :

i ​can't reproduce it with the current (14 Mar 2016) version of LAMMPS.​

when extracting the potential energy from the log, i get the following
graph.
did you look at the _entire_ log file?

axel.​

scan.png

thank you for your reply. I am using the 16 feb version. I tried extracting ‘pe’ as well but I still get all zeros. In your graph what is the y axis? is it the distance? Also, is there something I am missing in my input script. I want only the energy between two atoms, van der Walls energy, which is why i didnt provide any temperature or velocity. I also don’t perform any time integration.

Here is my datafile:

#Deposit file

clear
echo both

#-------------------------------INITIALISE SYSTEM-----------------------------
units real
atom_style charge
boundary p p p
timestep 0.25
#---------------------------DEFINE GROUPS AND AND VARIABLES-------------------
read_data input_equi.dat

group n2 type 2
group carbon type 1

#-------------------------------SET POTENTIAL-------------------------------
pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.reax C N
neigh_modify delay 1 check yes

velocity carbon set NULL NULL NULL units box
velocity n2 set NULL NULL NULL units box
compute reax all pair reax/c
variable ew equal c_reax[11]

thermo_style custom epair evdwl elong v_ew pe

#-------------------------------LOOP DEFINITIONS---------------------------
variable imax equal 1000
variable i loop ${imax}
label loopstart
displace_atoms n2 move 0 0 -0.05 units box
run 0
next i
jump SELF loopstart

thank you for your reply. I am using the 16 feb version. I tried
extracting 'pe' as well but I still get all zeros. In your graph what is
the y axis? is it the distance?

​then something is wrong with your input or your reaxff forcefield file (i
used the one for RDX as that contains C and N)

​y-axis is pairwise potential energy (epair).

Also, is there something I am missing in my input script. I want only the
energy between two atoms, van der Walls energy, which is why i didnt
provide any temperature or velocity. I also don't perform any time
integration.

i didn't change anything in your script.​​ i just used grep to extract the
thermo output (the original plot was made from data without the last grep,
the one attached from only the non-zero data listed below):

lmp_serial -in in.scan | grep -B1 ^Loop | grep '^ \+[-0-9]' | grep -v '^
\+0 '

-0.01780097 -0.01780097 0 -0.01780097
-0.02330534 -0.02330534 0 -0.02330534
-0.02777709 -0.02777709 0 -0.02777709
-0.10290092 -0.10290092 0 -0.10290092
-0.10944439 -0.10944439 0 -0.10944439
-0.11608291 -0.11608291 0 -0.11608291
-0.12274591 -0.12274591 0 -0.12274591
-0.12934439 -0.12934439 0 -0.12934439
-0.13576734 -0.13576734 0 -0.13576734
-0.14187774 -0.14187774 0 -0.14187774
-0.14750772 -0.14750772 0 -0.14750772
-0.15245305 -0.15245305 0 -0.15245305
  -0.1564666 -0.1564666 0 -0.1564666
-0.15925087 -0.15925087 0 -0.15925087
  -0.1604491 -0.1604491 0 -0.1604491
-0.15963518 -0.15963518 0 -0.15963518
-0.15630167 -0.15630167 0 -0.15630167
-0.14984612 -0.14984612 0 -0.14984612
-0.13955503 -0.13955503 0 -0.13955503
-0.12458536 -0.12458536 0 -0.12458536
-0.10394301 -0.10394301 0 -0.10394301
  -4.8250689 -4.8250689 0 -0.040755598
  -4.7938862 -4.7938862 0 0.004776899
  -4.7643883 -4.7643883 0 0.062028816
   -4.744242 -4.744242 0 0.13321061
  -4.7460138 -4.7460138 0 0.22090537
  -4.7886743 -4.7886743 0 0.32812832
  -4.8989344 -4.8989344 0 0.45839554
  -5.1120894 -5.1120894 0 0.61580318
  -5.4720299 -5.4720299 0 0.80511904
  -6.0301284 -6.0301284 0 1.031888
  -6.8428367 -6.8428367 0 1.3025538
    -7.96802 -7.96802 0 1.6245991
  -9.4602586 -9.4602586 0 2.006707
  -11.365548 -11.365548 0 2.4589475
  -13.715962 -13.715962 0 2.9929908
  -16.524874 -16.524874 0 3.6223545
  -19.783293 -19.783293 0 4.3626868
  -23.457698 -23.457698 0 5.232092
  -27.489571 -27.489571 0 6.2515057
   -31.79668 -31.79668 0 7.4451241
  -36.276481 -36.276481 0 8.8408983
  -40.814216 -40.814216 0 10.471101
  -45.304465 -45.304465 0 12.372975
    -49.7031 -49.7031 0 14.58948
   -54.12485 -54.12485 0 17.170141
  -58.985774 -58.985774 0 20.172026
  -65.206116 -65.206116 0 23.660861
  -74.605522 -74.605522 0 27.712297
  -90.256063 -90.256063 0 32.413367
  -112.44587 -112.44587 0 37.86413
   -133.3331 -133.3331 0 44.179554
  -150.24749 -150.24749 0 51.491638
  -165.19515 -165.19515 0 59.951826
  -176.60607 -176.60607 0 69.733724
  -184.04063 -184.04063 0 81.036155
  -187.82522 -187.82522 0 94.086584
    -190.705 -190.705 0 109.14494
  -216.17583 -216.17583 0 126.50782
    -254.679 -254.679 0 146.51317
  -244.94509 -244.94509 0 169.54527
  -218.61507 -218.61507 0 196.04014
  -185.40254 -185.40254 0 226.49123
   -148.9953 -148.9953 0 261.4551
  -108.58195 -108.58195 0 301.55702
  -62.793852 -62.793852 0 347.4959
  -10.515872 -10.515872 0 400.04792
   49.242657 49.242657 0 460.06784
   117.45504 117.45504 0 528.48625
    195.1225 195.1225 0 606.30048
   283.27971 283.27971 0 694.55508
   382.96994 382.96994 0 794.30617
   495.18806 495.18806 0 906.5602
   620.78034 620.78034 0 1032.1724
   760.27903 760.27903 0 1171.6815
   913.63412 913.63412 0 1325.0415
   1079.7798 1079.7798 0 1491.1893
   1255.9307 1255.9307 0 1667.3411
   1436.4323 1436.4323 0 1847.8429
   1610.8376 1610.8376 0 2022.2482
   1760.3479 1760.3479 0 2171.7585
   1844.7527 1844.7527 0 2256.1634
   1760.3479 1760.3479 0 2171.7585
   1610.8376 1610.8376 0 2022.2482
   1436.4323 1436.4323 0 1847.8429
   1255.9307 1255.9307 0 1667.3411
   1079.7798 1079.7798 0 1491.1893
   913.63412 913.63412 0 1325.0415
   760.27903 760.27903 0 1171.6815
   620.78034 620.78034 0 1032.1724
   495.18806 495.18806 0 906.5602
   382.96994 382.96994 0 794.30617
   283.27971 283.27971 0 694.55508
    195.1225 195.1225 0 606.30048
   117.45504 117.45504 0 528.48625
   49.242657 49.242657 0 460.06784
  -10.515872 -10.515872 0 400.04792
  -62.793852 -62.793852 0 347.4959
  -108.58195 -108.58195 0 301.55702
   -148.9953 -148.9953 0 261.4551
  -185.40254 -185.40254 0 226.49123
  -218.61507 -218.61507 0 196.04014
  -244.94509 -244.94509 0 169.54527
    -254.679 -254.679 0 146.51317
  -216.17583 -216.17583 0 126.50782
    -190.705 -190.705 0 109.14494
  -187.82522 -187.82522 0 94.086584
  -184.04063 -184.04063 0 81.036155
  -176.60607 -176.60607 0 69.733724
  -165.19515 -165.19515 0 59.951826
  -150.24749 -150.24749 0 51.491638
   -133.3331 -133.3331 0 44.179554
  -112.44587 -112.44587 0 37.86413
  -90.256063 -90.256063 0 32.413367
  -74.605522 -74.605522 0 27.712297
  -65.206116 -65.206116 0 23.660861
  -58.985774 -58.985774 0 20.172026
   -54.12485 -54.12485 0 17.170141
    -49.7031 -49.7031 0 14.58948
  -45.304465 -45.304465 0 12.372975
  -40.814216 -40.814216 0 10.471101
  -36.276481 -36.276481 0 8.8408983
   -31.79668 -31.79668 0 7.4451241
  -27.489571 -27.489571 0 6.2515057
  -23.457698 -23.457698 0 5.232092
  -19.783293 -19.783293 0 4.3626868
  -16.524874 -16.524874 0 3.6223545
  -13.715962 -13.715962 0 2.9929908
  -11.365548 -11.365548 0 2.4589475
  -9.4602586 -9.4602586 0 2.006707
    -7.96802 -7.96802 0 1.6245991
  -6.8428367 -6.8428367 0 1.3025538
  -6.0301284 -6.0301284 0 1.031888
  -5.4720299 -5.4720299 0 0.80511904
  -5.1120894 -5.1120894 0 0.61580318
  -4.8989344 -4.8989344 0 0.45839554
  -4.7886743 -4.7886743 0 0.32812832
  -4.7460138 -4.7460138 0 0.22090537
   -4.744242 -4.744242 0 0.13321061
  -4.7643883 -4.7643883 0 0.062028816
  -4.7938862 -4.7938862 0 0.004776899
  -4.8250689 -4.8250689 0 -0.040755598
-0.10394301 -0.10394301 0 -0.10394301
-0.12458536 -0.12458536 0 -0.12458536
-0.13955503 -0.13955503 0 -0.13955503
-0.14984612 -0.14984612 0 -0.14984612
-0.15630167 -0.15630167 0 -0.15630167
-0.15963518 -0.15963518 0 -0.15963518
  -0.1604491 -0.1604491 0 -0.1604491
-0.15925087 -0.15925087 0 -0.15925087
  -0.1564666 -0.1564666 0 -0.1564666
-0.15245305 -0.15245305 0 -0.15245305
-0.14750772 -0.14750772 0 -0.14750772
-0.14187774 -0.14187774 0 -0.14187774
-0.13576734 -0.13576734 0 -0.13576734
-0.12934439 -0.12934439 0 -0.12934439
-0.12274591 -0.12274591 0 -0.12274591
-0.11608291 -0.11608291 0 -0.11608291
-0.10944439 -0.10944439 0 -0.10944439
-0.10290092 -0.10290092 0 -0.10290092
-0.02777709 -0.02777709 0 -0.02777709
-0.02330534 -0.02330534 0 -0.02330534
-0.01780097 -0.01780097 0 -0.01780097

scan.png

I used the same potential that you mentioned, reax.rdx and got the same curve. I guess there was a problem with my potential. Thank you for your help.

Good day!

As I mentioned in my previous emails, I am trying to output the van der Walls forces between one carbon and nitrogen atom with respect to distance. I noticed that using pairstyles reax and reax/c are giving me different values for lowest energies. I am attaching the graphs for the two as well as the graph I plot using the analytical formula that is stated in this paper “ReaxFF: A Reactive Force Field for Hydrocarbons” http://pubs.acs.org/doi/abs/10.1021/jp004368u. I used the parameters provided in my potential file.

I tried looking into the source files for USER-reax/c and was trying to understand how van der Walls energy is calculated and why I am getting such a high (negative) values when there is absolutely no velocity/temperature provided. The analytical result shows that a minimum of about -0.15 kcal/mol occurs at ~3.8 A. Why is there such a large difference between the analytical and lammps-obtained values? And also why is there a difference between the results using of reax and reax/c pair styles?

Attached below is my input and data files. Thank you.

#Input file

clear
echo both

#-------------------------------INITIALISE SYSTEM-----------------------------
units real
atom_style charge
boundary p p p
timestep 0.25
atom_modify map array
read_data input.dat

#---------------------------DEFINE GROUPS AND AND VARIABLES-------------------
group n type 2
group carbon type 1

#-------------------------------SET POTENTIAL-------------------------------
#pair_style reax/c NULL checkqeq yes
#pair_coeff * * ffield.reax C N

pair_style reax
pair_coeff * * ffield.reax 1 4

compute 1 all property/atom z
compute 2 all pe/atom pair

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

thermo_style custom epair evdwl elong ecoul
fix print all print 1 “((c_1[1]) - (c_1[2])) (evdwl)” file energy_r.dat screen no
fix print2 all print 1 “(evdwl) (ecoul) $(epair)” file vdw_coul_pair.dat screen no

#-------------------------------LOOP DEFINITIONS---------------------------

variable imax equal 5000
variable i loop ${imax}

label loopstart
displace_atoms carbon move 0 0 0.005 units box

run 1
next i
jump SELF loopstart

Date file:

#input.dat

2 atoms
2 atom types

0 50 xlo xhi
0 50 ylo yhi
0 50 zlo zhi

Masses

1 12.0107
2 14.0067

Atoms

1 1 0 25 25 25
2 2 0 25 25 25.1

analytical.png

pairstyle_comparision.png

Ray can probably answer the Q as to why pair_style reax vs reax/c
would be giving different answers. That should almost
never be the case.

Steve