Implementation of basic example in LAMMPS [closed]

Question

I wanted to implement a basic example from the book "The Art of molecular dynamics simulation" in LAMMPS. It describes the 2d-movement of molecules in a Lennard-Jones-potential, defined as $$u(r_{ij})=\left\{\begin{matrix}4\varepsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]+\varepsilon&r_{ij}<r_c=2^{1/6}\sigma\\0&r_{ij}\geq r_c\end{matrix}\right. $$ using temperature and density as parameters. The dimensions of the 2d-area are $\left(20, 20\right)$, having a molecule sitting at each of the nodes, resulting in 400 molecules. The values measured are the total kinetic energy, the total potential energy and the total pressure.

My implementation in lammps thus is

#First test, if I am able to write a lammps script on my own
#------------------Init-----------------
clear
units metal
dimension 2
boundary p p p
atom_style body nparticle 2 6
atom_modify map array#Array or Hash

#------------Create Atoms----------------
lattice sq 1
region box block 0 20 0 20 0 1 units lattice
create_box 1 box

lattice sq 1 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
replicate 1 1 1

#-----------Interatomic potential---------
pair_style lj/cut 1.122466
pair_coeff * * 1 1
neighbor 0.4 bin
velocity all create 1.44 87287 loop geom

#------------Define Settings---------------
compute eng all pe/atom
compute keng all ke/atom
compute eatoms all reduce sum c_eng
compute keatoms all reduce sum c_keng
#compute peng all pressure

#-------------Run minimization--------------
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal "count(all)"
variable teng equals "c_eatoms"
variable length equals "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms) = ${length};"
print "Cohesive energy (eV) = ${ecoh};"
print "All done!"

Is that implementation approximately correct, or is there something completely wrong?
When executing the file, the execution stops with the error

ERROR: Illegal variable command (../variable.cpp:512)

How can I debug that?

Answer 1

It seems your questions is more on "How to debug that error" than on the error itself. If you wish me to look into the error specifically, I will do this, I promise :).

Otherwise, it is great to learn from your own mistakes and to learn how to debug yourself. The best way to debug a lammps script is to use the command:

echo both

at the beggining of the script. This command will output each line processed by lammps one by one to the terminal. You must understand that the script is read and interpreted line by line. So if you put echo both, each line that is interpreted will be displayed as it is interpreted. This way, when the code crash, the line above it will be the last one that has been processed and the error will most likely reside within that line. In your case, I believe your error is most likely related to one of your compute command.