lj/cut force default correction?

Hi all,

I’m using LAMMPS to compute the LJ (periodic) force between 2 atoms to compare with a code I’m writing. I’ve thoroughly checked all inputs to be identical. Two Ar atoms 3A apart.

In the PairLJCut::compute function, I printf’d the force and recompiled/re-ran: -13.337237. My code also gets this answer.

But when LAMMPS dumps a file with atom forces via
variable fx equal fx[1]
fix print2 all print 1 “${fx}” file f.txt screen no title "fx"

I get -13.3280278708918 (in that f.txt file)

Is there some default correction to the force for lj/cut ? It’s driving me crazy

Please assist. I’m copying below my input and data files.


LAMMPS input generated by MCMD on 10-10-2018 at 15:40:54

some variables for easy scripting later.

variable temperature equal 0
variable freq equal 1
variable nstep equal 5

some global params

units real
boundary p p p
atom_style full

the atoms file

read_data lammps.data

site (pseudo-atom) masses

mass 1 39.94000 #Ar

timestep 1.000000 # femptoseconds (1e-15 s)

methods for potential/force calculations

pair_style lj/cut 500 # this gets overwritten by pair param’s below

Atomic parameters (used for mixing)

#pair_coeff 1 1 0.238067 3.405000 8.512500
pair_coeff 1 1 0.23806705118 3.405000 498

Apply LB mixing rules.

pair_modify mix arithmetic


group moving molecule > 0 # for a MOF simulation, the MOF is molecule 1


#compute rdf all rdf 50 [pairs… e.g. 8 1 8 8…]

more variables etc. useful for MD

variable step equal step
variable time equal step2
variable timeNS equal time/1000000
compute themsd all msd
variable diffusion_coeff equal c_themsd[4]/(6
time)*0.1 # cm^2/s
variable te equal c_pe+c_ke
compute pe all pe
compute ke all ke
compute movingtemp moving temp
thermo_style custom step etotal ke pe evdwl ecoul
thermo ${freq}

fix printdiff all print 1 “{diffusion_coeff}" file D.txt screen no title "D_Coeff_cm^2/s" #fix printmsd all print 1 "{themsd[4]}” file msd.txt screen no title “MSD_A^2”

variable vx equal vx[1]
variable vy equal vy[1]
variable vz equal vz[1]
fix printmyatom all print 1 “{vx} {vy} ${vz}” file v.txt screen no title “vx vy vz”

variable fx equal fx[1]
variable fy equal fy[1]
variable fz equal fz[1]
fix print2 all print 1 “{fx} {fy} ${fz}” file f.txt screen no title “fx fy fz”

variable x equal x[1]
variable y equal y[1]
variable z equal z[1]
fix print3 all print 1 “{x} {y} ${z}” file x.txt screen no title “x y z”

variable rsq equal (x[1]-x[2])(x[1]-x[2])+(y[1]-y[2])(y[1]-y[2])+(z[1]-z[2])*(z[1]-z[2])
fix print4 all print 1 “${rsq}” file rsq.txt screen no title “rsq”

write files (dumps)

dump Dump all xyz ${freq} lammps_traj.xyz
dump_modify Dump sort id
dump_modify Dump element Ar

set NVT

velocity all create ${temperature} 12345 rot yes mom yes dist gaussian
fix nve moving nve


run ${nstep}

Data file:


2 atoms
1 atom types

-500 500 xlo xhi
-500 500 ylo yhi
-500 500 zlo zhi


1 1 1 0.00000 0.000 0.00000 0.00000 #Ar
2 2 1 0.00000 3.000 0.00000 0.00000 #Ar

I just realized,

The file output does not include the initially computed force before the 1st step happens. So the discrepancy is just that F(t=0) is not shown in the output file f.txt

You can also dump the force on each atom and it will include
the step 0 forces. For a 2-atom system that should
be what you want.