[lammps-users] Energy Conservation Question

LAMMPS LAMMPS Mailing List Mirror
1

I am running a simple simulation with about 4000 atoms using LJ pair
potential. My command script is below. When I look at the total energy
(TotEng) of the system from the log file, there appears to be
fluctuations that are on the order of 0.4% of KE of the system if I have
the fix temp/rescale commented out (lines with #*** following them). Is
this sort of fluctuation in the energy typical? When I uncomment the
fix temp/rescale commands (lines with #*** following them), I get even
larger fluctuations that approach 2% of the KE of the system. My
understanding is that the "fix_modify fixID energy yes" would account
for all the energy added or subtracted from the group of atoms and add
is to the potential energy of the system so that energy conservation of
the entire system can be monitored. Is this correct or am I missing
something? If it is correct, why might I be getting such large
fluctuations in my total system energy (TotEng)? Could it be that the
energy added from the fix_modify command is only adding the energy from
the last timestep, opposed to all the energy added since the last log
printing?

Also if I use "fix_modify fixID themo yes", where the fixID is the ID
of the fix temp/rescale command, is the output to my log file the energy
the fix added to the system since the last log print or since the last
timestep?

A side question, is there a way to introduce a random seed into the
velocity command within in LAMMPS or is this something I must do
external to the program?

Thank you,

Rob

#3D LJ simulation for TBR

echo both

variable d equal 010307 #Change this to be date

variable s equal 2 #Change this to be scan number

log simulation$dscan$s.dat

#Initialization, units are Angstrom, ps, grams/mole, eV, K

units metal

boundary p p f

atom_style atomic

pair_style lj/cut 8.5125

#Atom definitions

lattice fcc 5.2865

region afixed block 0 5 0 5 -20 -18.01

region aheat block 0 5 0 5 -18.0 -16.01

region abox block 0 5 0 5 -16.0 -.01

region bbox block 0 5 0 5 0 15.99

region bcool block 0 5 0 5 16 17.99

region bfixed block 0 5 0 5 18.0 20

region aside block 0 5 0 5 -20 -0.0001

region bside block 0 5 0 5 0 20

region worldbox block 0 5 0 5 -20 20

region activebox block 0 5 0 5 -18 18

create_box 2 worldbox

create_atoms 1 region aside

create_atoms 2 region bside

mass 1 40

mass 2 80

group gaheat region aheat

group gaside region abox

group gafixed region afixed

group gbcool region bcool

group gbside region bbox

group gbfixed region bfixed

group gactivebox region activebox

velocity gaheat create 46 4932345 # mom yes rot yes dist
gaussian

velocity gaside create 46 8932345 # mom yes rot yes dist
gaussian

velocity gbside create 34 6932345 # mom yes rot yes dist
gaussian

velocity gbcool create 34 2932345 # mom yes rot yes dist
gaussian

#Settings

pair_coeff * * .0103 3.405 8.5125

neighbor 2.0 bin

neigh_modify every 20 delay 0 check yes

fix heat gaheat temp/rescale 1 23 23 0.01 1 #***

fix cool gbcool temp/rescale 1 17 17 0.01 1 #***

timestep 0.005

fix 1 gactivebox nve

fix_modify heat thermo yes energy yes #***

fix_modify cool thermo yes energy yes #***

temperature 1 gaheat full

temperature 2 gaside full

temperature 3 gbside full

temperature 4 gbcool full

thermo 50

run 100000

2

Hi,
I am trying to print out the potential (energy and forces) as a function
of interatomic distance using the pair_write command. I am using the EAM
potential (cuu3.eam). My question is whether LAMMPS prints out the total
energy (both embedded term and the interatomic term) or only the
interatomic term? I am confused because when I plot the printed
energies, there is no attractive region.

Thanks,
Arun

3

You may be getting fluctutations b/c you are using LJ with a cutoff.
If you use pair_modify shift yes, you can eliminate the energy
discontinuity at the cutoff.

I think the energies printed by fix temp/rescale are only instantaneous
(on the timestep printed), not cummulative. So they may be
less useful for monitoring eng conservation.

Generally, if your system is at equilibrium, then Eng fluctuations
are OK, but you shouldn't get drift in E or T, even if the temp/rescale
is turned off. If you do, try decreasing the timestep. Also the
fluctuations should decrease with increasing N. 4000 atoms is not
a lot.

Re: velocity random seeds - I'm not sure what you're asking. The
velocity create command requires a random number seed entered
by the user. This is to allow for reproducibility from run to
run and across machines. If you are asking if the seeds can be
generated internally so they are invisible to the user, then LAMMPS
doesn't currently do that.

Steve

4

For pair_write, only the pair term (not the embedding term) is
computed for the EAM potentials. This is b/c the embedding
term typically comes from the collection of surrounding atoms, not
from a single i-j interaction.

However, when you print per-atom energy (via the dump custom command),
you will get the total eng for each atom (embedding + pair).

Steve