Hi everyone,
I am trying to compare some custom Fortran code for a certain system where there are WCA interactions besides others with LAMMPS. I am looking at energy values to see the correctness of my code, but simplifying the system to see where the mismatch could be, I found that in LAMMPS for a simple case where I have two single atoms at a distance r=\\sigma and they only interact with the pair style lj/cut, the potential energy is not \\epsilon but \\epsilon/2. I don’t understand why this happens, because I thought at that distance the PE should be \\epsilon by definition. This is the script I am employing (minimal test).
units lj
dimension 2
boundary p p p
atom_style atomic
newton on
#Box and two atoms
region box block -5 5 -5 5 -0.5 0.5
create_box 1 box
create_atoms 1 single 0.0 0.0 0.0
create_atoms 1 single 1.0 0.0 0.0 # distance = 1 = sigma
mass 1 1.0
#Pair style: WCA (LJ truncated at 2^(1/6))
variable rcut equal 1.122462
pair_style lj/cut ${rcut}
pair_coeff * * 1.0 1.0 ${rcut}
pair_modify shift yes
#Thermo output
thermo_style custom step temp pe
thermo 1
#Minimal integration to trigger energy computation
velocity all set 0.0 0.0 0.0
fix 1 all nve
#Run a single step
run 1
And this is the output:
LAMMPS (29 Aug 2024)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
Created orthogonal box = (-5 -5 -0.5) to (5 5 0.5)
1 by 1 by 1 MPI processor grid
Created 1 atoms
using lattice units in orthogonal box = (-5 -5 -0.5) to (5 5 0.5)
create_atoms CPU = 0.000 seconds
Created 1 atoms
using lattice units in orthogonal box = (-5 -5 -0.5) to (5 5 0.5)
create_atoms CPU = 0.000 seconds
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info …
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.422462
ghost atom cutoff = 1.422462
binsize = 0.711231, bins = 15 15 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/2d
bin: standard
Setting up Verlet run …
Unit style : lj
Current step : 0
Time step : 0.005
Per MPI rank memory allocation (min/avg/max) = 3.048 | 3.048 | 3.048 Mbytes
Step Temp PotEng
0 0 0.5
1 0.014237099 0.49284091
Loop time of 2.9475e-05 on 1 procs for 1 steps with 2 atoms
Performance: 14656488.550 tau/day, 33927.057 timesteps/s, 67.854 katom-step/s
44.1% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
Pair | 1.257e-06 | 1.257e-06 | 1.257e-06 | 0.0 | 4.26
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 3.353e-06 | 3.353e-06 | 3.353e-06 | 0.0 | 11.38
Output | 1.1315e-05 | 1.1315e-05 | 1.1315e-05 | 0.0 | 38.39
Modify | 4.191e-06 | 4.191e-06 | 4.191e-06 | 0.0 | 14.22
Other | | 9.359e-06 | | | 31.75
Nlocal: 2 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1 ave 1 max 1 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1
Ave neighs/atom = 0.5
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00
I would appreciate some help with this because it is very strange to me.
Thanks,
Sohan.