First of all, I’m very new to LAMMPS so I’m sorry if my question is trivial.
I’m trying to model a gas of atoms that float in a square box, they obey an NVE integration and only feel a shifted Lennard-Jones potential, no other interaction is introduced. What I’m trying to do is to fill a file with energy values (alongside times) in order to plot them and see if energy is conserved. To do so, I use the fix print function, where previously I defined a variable Etot which should contain the energy I need, the full fix print line I’m using is:
fix print_energy all print 1000 "$t ${Etot}" file data_13_5_a.txt screen no
Which should create the file data_13_5_a.txt with times alongside energies. The problem is that I receive an error:
ERROR on proc 0: Substitution for illegal variable Etot (../input.cpp:617
I think the simulation itself is working because the dump file (see code below) .lammpstrj gives a result in Ovito.
In the end I also tried with another fix print (see below) but it gives me a similar error (Substitution for illegal variable K).
I don’t know what I’m doing wrong, for reference and errors, here is the full code I wrote:
#-------------------------------------------Parameters---------------------------------------------#
# Number of integration timesteps
variable steps equal 1e6
# Integration timestep
variable dt equal 5e-3
# Truncation radius and LJ parameters
variable epsilon equal 1
variable sigma equal 1
variable Rc equal 2.5
# Length of cubic simulation box
variable L equal 10.
# Number of particles
variable rho equal 0.3
variable V equal $L*$L*$L
variable N equal ${rho}*$V
# Temperature
variable T equal 2.
# Particles' mass
variable m equal 1.
# Seed for the RNG
variable seed equal 123
#-----------------------------------------Define Styles--------------------------------------------#
# Atom style used
atom_style atomic
#---------------------------------------Create cubic box-------------------------------------------#
# Define cubic region
region sim_box block -$(0.5 * v_L) $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L)
# Create simulation box based on region
create_box 1 sim_box
#-----------------------------------------Create Atoms---------------------------------------------#
# Place particles in box randomly
create_atoms 1 random $N 123 NULL
#---------------------------------------------Masses-----------------------------------------------#
mass 1 $m
#----------------------------------------Pair Coefficients-----------------------------------------#
# Truncated LJ interactions for repulsion among particles
pair_style lj/cut ${Rc}
# Shift the LJ potential to 0.0 at the cutoff (thus no discontinuities). The coefficients should
# be the ones of the LJ interaction formula: epsilon and sigma, also another for the shift.
pair_modify shift yes
# Coefficients of the pair interactions
pair_coeff 1 1 ${epsilon} ${sigma} ${Rc}
# neighbor $(2.*v_Rc) bin # not needed but may be useful
#-------------------------------------------Minimize-----------------------------------------------#
# minimise to improve the random distance at which atoms are placed and avoid errors.
minimize 1e-4 1e-7 1000 1000
#---------------------------------------------Fixes------------------------------------------------#
# Initialize velocities to MB distribution and zero total momentum
velocity all create $T ${seed} mom yes
# Microcanonical ensemble (total energy constant)
fix 1 all nve
#----------------------------------------Measurements----------------------------------------------#
# Include this if for visualisation with vmd.
dump movie all atom 100 ./movie_13_5_a.lammpstrj
# Current timestep
variable t equal step*${dt}
#--------------------------------------------Energies----------------------------------------------#
# group particles type 1 # test
# calculate total kinetic energy of the system, store it in the compute variable c_Kenergy
compute Kenergy all ke
# calculate total potential energy of the system, store it in the compute variable c_Venergy
compute Venergy all pe
# Store the energies (kinetic, potential and total) in three variables which can be accessed
# by the print command. Don't set Etot equal $K+$V, because LAMMPS will complain...
variable K equal c_Kenergy
variable V equal c_Venergy
variable Etot equal c_Kenergy+c_Venergy
#---------------------------------------------Run--------------------------------------------------#
# set timestep for integration
timestep ${dt}
# Write to a file
fix print_energy all print 1000 "$t ${Etot}" file "./data_13_5_a.txt" screen no
# fix energy_print all print 10 "$t $K $V $E" file energy.dat screen no title "#Timestep Kinetic Potential Total"
run ${steps}
Furthermore, I’m using Mac, and to run the simulation I wrote in the terminal lmp_serial -in input_13_5_a where input_13_5_a is the file containing my script.