6
$\begingroup$

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:

####################################################################################################
#                                      Lennard-Jones fluid                                         #
####################################################################################################


#--------------------------------------Dordoni Luca 0880068----------------------------------------#
#----------------------------------------------13.5a-----------------------------------------------#


#-------------------------------------------Parameters---------------------------------------------#

# Number of integration timesteps
variable steps equal 1e5

# 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 ${seed} 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-----------------------------------------------#


# minimize to improve the random distance at which atoms are placed and avoid errors.
minimize 1e-4 1e-7 10000 10000

reset_timestep 0


#---------------------------------------------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 10 ./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

# thermo_style custom time etotal
# thermo 200
 

#-----------------------------------------------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.

Edit. Here is the log file:

LAMMPS (29 Sep 2021 - Update 1)
####################################################################################################
#                                      Lennard-Jones fluid                                         #
####################################################################################################


#--------------------------------------Dordoni Luca 0880068----------------------------------------#
#----------------------------------------------13.5a-----------------------------------------------#


#-------------------------------------------Parameters---------------------------------------------#

# Number of integration timesteps
variable steps equal 1e5

# 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 V equal 10*$L*$L
variable V equal 10*10*$L
variable V equal 10*10*10
variable N equal ${rho}*$V
variable N equal 0.3*$V
variable N equal 0.3*1000

# 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)
region sim_box block -5 $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L)
region sim_box block -5 5 -$(0.5 * v_L) $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L)
region sim_box block -5 5 -5 $(0.5 * v_L) -$(0.5 * v_L) $(0.5 * v_L)
region sim_box block -5 5 -5 5 -$(0.5 * v_L) $(0.5 * v_L)
region sim_box block -5 5 -5 5 -5 $(0.5 * v_L)
region sim_box block -5 5 -5 5 -5 5

# Create simulation box based on region
create_box 1 sim_box
Created orthogonal box = (-5.0000000 -5.0000000 -5.0000000) to (5.0000000 5.0000000 5.0000000)
  1 by 1 by 1 MPI processor grid


#-----------------------------------------Create Atoms---------------------------------------------#


# Place particles in box randomly
create_atoms 1 random $N ${seed} NULL
create_atoms 1 random 300 ${seed} NULL
create_atoms 1 random 300 123 NULL
Created 300 atoms
  using lattice units in orthogonal box = (-5.0000000 -5.0000000 -5.0000000) to (5.0000000 5.0000000 5.0000000)
  create_atoms CPU = 0.000 seconds


#---------------------------------------------Masses-----------------------------------------------#


mass 1 $m
mass 1 1


#----------------------------------------Pair Coefficients-----------------------------------------#


# Truncated LJ interactions for repulsion among particles
pair_style lj/cut ${Rc}
pair_style lj/cut 2.5

# 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}
pair_coeff 1 1 1 ${sigma} ${Rc}
pair_coeff 1 1 1 1 ${Rc}
pair_coeff 1 1 1 1 2.5

# neighbor $(2.*v_Rc) bin        # not needed but may be useful


#-------------------------------------------Minimize-----------------------------------------------#


# minimize to improve the random distance at which atoms are placed and avoid errors.
minimize 1e-4 1e-7 10000 10000
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:188)
Neighbor list info ...
  update every 1 steps, delay 0 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 2.8
  ghost atom cutoff = 2.8
  binsize = 1.4, bins = 8 8 8
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.214 | 4.214 | 4.214 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0            0     31424599            0     31424599     37710155 
     503            0   -5.2799826            0   -5.2799826  -0.12954156 
Loop time of 0.106709 on 1 procs for 503 steps with 300 atoms

99.8% CPU use with 1 MPI tasks x no OpenMP threads

Minimization stats:
  Stopping criterion = energy tolerance
  Energy initial, next-to-last, final = 
      31424598.7313598  -5.27947676473381  -5.27998264715565
  Force two-norm initial, final = 7.6269644e+11 7.3211830
  Force max component initial, final = 4.5240368e+11 1.1862632
  Final line search alpha, max atom move = 0.017969632 0.021316714
  Iterations, force evaluations = 503 1005

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.065804   | 0.065804   | 0.065804   |   0.0 | 61.67
Neigh   | 0.030718   | 0.030718   | 0.030718   |   0.0 | 28.79
Comm    | 0.0043468  | 0.0043468  | 0.0043468  |   0.0 |  4.07
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 0.005841   |            |       |  5.47

Nlocal:        300.000 ave         300 max         300 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:        840.000 ave         840 max         840 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:        7939.00 ave        7939 max        7939 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 7939
Ave neighs/atom = 26.463333
Neighbor list builds = 146
Dangerous builds = 1

reset_timestep 0


#---------------------------------------------Fixes------------------------------------------------#


# Initialize velocities to MB distribution and zero total momentum
velocity all create $T ${seed} mom yes
velocity all create 2 ${seed} mom yes
velocity all create 2 123 mom yes

# Microcanonical ensemble (total energy constant)
fix 1 all nve


#------------------------------------------Measurements--------------------------------------------#


# Include this if for visualisation with vmd.
dump movie all atom 10 ./movie_13_5_a.lammpstrj

# Current timestep
variable t equal step*${dt}
variable t equal step*0.005


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

# thermo_style custom time etotal
# thermo 200


#-----------------------------------------------Run------------------------------------------------#


# set timestep for integration
timestep ${dt}
timestep 0.005

# 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}
run 100000
ERROR on proc 0: Substitution for illegal variable Etot (../input.cpp:617)
$\endgroup$
11
  • $\begingroup$ +1. Welcome to our new community and thank you very much for contributing your question here. We hope to see much more of you in the future!!! $\endgroup$ Jan 2, 2022 at 22:26
  • 2
    $\begingroup$ Running the same input script as provided by you, the simulation does not give any error with my LAMMPS installation. Can you upload the log file printed by LAMMPS as well? $\endgroup$
    – PBH
    Jan 3, 2022 at 1:34
  • 2
    $\begingroup$ It seems that this could be a problem with your build (29 Sep 2021). I suggest that you try with a different version of LAMMPS. The calculation ends without any errors in my case. $\endgroup$
    – PBH
    Jan 3, 2022 at 13:38
  • 2
    $\begingroup$ for debugging, consider this (which will print KE, PE, and Etotal to the log file every 100 steps). First command: thermo_style custom step ke pe etotal Second command: thermo 100 $\endgroup$ Jan 14, 2022 at 7:47
  • 1
    $\begingroup$ I'm not sure what the deal is with the print, possibly syntax. But now that I'm thinking about this again, you could also try (a) a small modification to include the variable flag v_ to "\$t \${v_Etot}" in the print command, or (b) calculating the total energy directly from the two computes inside the print statement, or (c) using the fix ave/time command with the KE and PE computes and/or the total energy variable and a file output instead of the print. $\endgroup$ Jan 20, 2022 at 23:13

0