计算科学 - 在 LAMMPS 中实现基本示例 - 吾爱随笔录

我想在 LAMMPS 中实现“分子动力学模拟的艺术”一书中的一个基本示例。它描述了 Lennard-Jones 势中分子的二维运动，定义为

u (r_{i j}) = {\begin{matrix} 4 ε [{(\frac{σ}{r_{i j}})}^{12} - {(\frac{σ}{r_{i j}})}^{6}] + ε & r_{i j} < r_{c} = 2^{1 / 6} σ \\ 0 & r_{i j} \geq r_{c} \end{matrix}

$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.$ 使用温度和密度作为参数。二维区域的尺寸是

(20, 20)

$\left(20, 20\right)$ ，每个节点都有一个分子，产生 400 个分子。测量的值是总动能、总势能和总压力。

因此，我在 lamps 中的实现是

#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!"

该实现大致正确，还是完全错误？
执行文件时，执行停止并出现错误

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

我该如何调试呢？