I have made a triangular lattice. And all the points are 1 unit of distance apart and I have set the bond is harmonic with rest length of 1 so without adding any other potential or force the force on each particle should be 0 and I have calculated it using a md code and it comes out to be zero. But then I run lammps for the same configuration the forces in dump file is not zero. Here is the basic input script I am using.
dimension 2
units lj
atom_style bond
# Read initial configuration
read_data init_config.txt
# Define bonds (springs)
bond_style harmonic
bond_coeff 1 200.0 1
# Define the neighbor list
neighbor 0.3 bin
neigh_modify delay 5
# Output atom properties
dump 1 all custom 100 trajectory/atoms.lammpstrj id type x y z
dump 2 all custom 100 forces/forces id x y z fx fy fz
# Specify simulation settings
# thermo 100
# Run the simulation for the required duration
run 1000
# Write a restart file
write_data restart.txt
And here is the output of dump file.
ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
100
ITEM: BOX BOUNDS pp pp pp
-5.0000000000000000e+00 1.4500000000000000e+01
-5.0000000000000000e+00 1.2794200000000000e+01
-5.0000000000000000e-01 5.0000000000000000e-01
ITEM: ATOMS id x y z fx fy fz
1 0 0 0 0 0 0
2 1 0 0 0 0 0
3 2 0 0 0 0 0
4 3 0 0 0 0 0
5 4 0 0 0 0 0
6 5 0 0 0 0 0
7 6 0 0 0 0 0
8 7 0 0 0 0 0
9 8 0 0 0 0 0
10 9 0 0 0 0 0
11 0.5 0.866025 0 0 0 0
12 1.5 0.866025 0 0 0 0
13 2.5 0.866025 0 0 0 0
14 3.5 0.866025 0 0 0 0
15 4.5 0.866025 0 0 0 0
16 5.5 0.866025 0 0 0 0
17 6.5 0.866025 0 0 0 0
18 7.5 0.866025 0 0 0 0
19 8.5 0.866025 0 0 0 0
20 9.5 0.866025 0 0 0 0
21 0 1.73205 0 0 0 0
22 1 1.73205 0 0 0.00300001 0
23 2 1.73205 0 0 0.00300001 0
24 3 1.73205 0 0 0.00300001 0
25 4 1.73205 0 0 0.00300001 0
26 5 1.73205 0 0 0.00300001 0
27 6 1.73205 0 0 0.00300001 0
28 7 1.73205 0 0 0.00300001 0
29 8 1.73205 0 0 0.00300001 0
30 9 1.73205 0 0 0 0
I don’t know what I am doing wrong here.