2d energy comparison vs HOOMD


I am using the December 15, 2015, intel complied version of lammps and have come into a pickle with comparing potential energies with HOOMD.

My inputs and data file are shown below after this section.

I am switching from HOOMD to LAMMPS because of some core pinning issues and want to get the correct energies for a 2D simulation of a simple model. The potential energies from the HOOMD run and lammps run are not consistent(-4.7 reduced energy units vs -0.5). I calculated the potential energy by hand to see which energy it was nearest too and it was hoomd. I feel confident in my calculation of the potential energy( ~ -4.6) so I am wondering what I am doing wrong in my process in LAMMPS. I checked both end results which give me the roughly the same final configuration and distance between atoms (important in verifying my potential energy calculation). I have interconverted between lammps and hoomd before in 3D and had no problems but that was with units real in LAMMPS in which I constructed similar units in HOOMD. Wondering if there might be a bug or if I am not accounting for something right.


Michael Goytia

units lj
atom_style full
boundary p p p
dimension 2
variable tempi equal 0.75
variable tempf equal 0.0123

read_data random.data

pair_style lj/cut 2.5
pair_coeff 1 1 0.182097851198 1.04308220858 2.60770552144
pair_coeff 1 2 0.702854451025 0.866639911779 2.16659977945
pair_coeff 1 3 0.973685740843 0.777251743894 1.94312935973
pair_coeff 2 2 0.169736804547 0.690197614981 1.72549403745
pair_coeff 2 3 0.570151175955 0.600809447096 1.50202361774
pair_coeff 3 3 0.24318312159 0.600809447096 1.27855319803

#pair_modify shift yes
neigh_modify exclude molecule all
timestep 0.005

fix 2d all enforce2d

thermo 1000
thermo_style custom step temp pe ke etotal press vol

dump myDump all dcd 1000 random.dcd
write_dump all custom random.lammpstrj id mol type x y z
restart 100000 2d.restart
fix 1 all rigid molecule langevin \{tempi\} {tempf} 500 10
run 4000000

.data file
#Goytia's random molecule generator

6 atoms
3 atom types

-6.0 6.0 xlo xhi
-6.0 6.0 ylo yhi
-5.0 5.0 zlo zhi


1 5.00065548514
2 1.44874411347
3 0.589396705833


1 1 1 0 -5.01247611056 0 0
2 1 2 0 -3.5 0 0
3 1 3 0 -3.83528420926 0.66167837418 0
4 2 1 0 1.98752388944 0 0
5 2 2 0 3.5 0 0
6 2 3 0 3.16471579074 0.66167837418 0

Is the energy the same between the 2 codes on timestep 0?

You are running a simulation for a 6-atom system? How

many bodies?

Do you get the same answer if you don’t define rigid bodies?


what happens, when comment out the following line?

neigh_modify exclude molecule all


what happens, when comment out the following line?
neigh_modify exclude molecule all

Axel and Steve,

1)The energy goes down but still not near to 4.7(~0.98)

2)Both trajectories have the same energy at time step 0,

3)The system contains 2 rigid bodies.

4)In order to not define rigid bodies and integrate are you suggesting adding bonds, angles, and dihedrals and then run with another integrator?


do you see the same issue with the 17 Nov 2016 version?


We currently don’t have installed at the moment but I will try that and get back to you.

Michael Goytia

Hi all,

I think I found a bug in dump_modify. When the option pbc is set to yes and the box is changed to triclinic the dump file does not contain the right coordinates. I am using version 9 Nov 2016. I am attaching an input that reproduces the error.


in.small_bcc (875 Bytes)

ATT00001.txt (4.33 KB)