Potential energy of granular; modification advice

I know that none of the granular system() functions bother to return a potential energy since the potentials are dissipative.

For my situation of a 1D chain of non-rotating spheres using gran/hertz/history, where k_t, gamma_t and gamma_n are zero, I want the potential energy to make sure that total energy is conserved.

That was easy enough to fix in the code. For pairs, this:

pe = ccel*(radsum-r)*2.0/5.0;

was added to file pair_gran_hertz_history.cpp, and returned from the system() function and added to the pair energy given to ev_tally_xyz() in the compute() function.

For the walls,

pe = ccel*(radius-r)*2.0/5.0;

was added in file fix_wall_gran.cpp, and the value added to the atom's energy in the post_force() function;

double *eatom = force->pair->eatom;
...
eatom[i] += pe/2.0;

The total energy is now conserved, however two issues remain;

1. The grain-grain pair potential energy shows up correctly in the thermo and dump (see input file below), but the grain-wall energy only shows up in the dump of compute c4.

How can I get the thermo output to recalculate the pairs energy to report the total potential energy correctly? Am I just out of sequence?

2. The potential is off by an exact multiplicative factor, and I can't figure out why. Any ideas?

Thanks!

Thad

I know that none of the granular system() functions bother to return a
potential energy since the potentials are dissipative.

For my situation of a 1D chain of non-rotating spheres using
gran/hertz/history, where k_t, gamma_t and gamma_n are zero, I want the
potential energy to make sure that total energy is conserved.

That was easy enough to fix in the code. For pairs, this:

pe = ccel*(radsum-r)*2.0/5.0;

was added to file pair_gran_hertz_history.cpp, and returned from the
system() function and added to the pair energy given to ev_tally_xyz()
in the compute() function.

For the walls,

pe = ccel*(radius-r)*2.0/5.0;

was added in file fix_wall_gran.cpp, and the value added to the atom's
energy in the post_force() function;

double *eatom = force->pair->eatom;
...
eatom[i] += pe/2.0;

The total energy is now conserved, however two issues remain;

1. The grain-grain pair potential energy shows up correctly in the
thermo and dump (see input file below), but the grain-wall energy only
shows up in the dump of compute c4.

How can I get the thermo output to recalculate the pairs energy to
report the total potential energy correctly? Am I just out of sequence?

energy contributions from fixes are usually not added to the total
energy. you need to use: fix_modify energy yes.
http://lammps.sandia.gov/doc/fix_modify.html

axel.

And there is a line like this you need to include
in the Fix::set mask() method:

mask |= THERMO_ENERGY;

Steve

Thanks Steve and Alex, but I just ran a test with this advice and no luck.

Below, c_c4 has the potential energy I added to the fix code, but c_c5 does not.
Not a total loss, I can dig the correct pe out of the dumped data.

Thad

Thanks Steve and Alex, but I just ran a test with this advice and no luck.

since you are using modified source code, it is difficult to debug
this from remote. if you are doing it correctly, it should work as
advertised. what you put in your input script is only a part of the
story.

axel.

In FixWallGran::post_force(int vflag) I add:

double pe;
double *eatom = force->pair->eatom;

and inside the loop of local atoms,

pe = hertz_history(...)

now returns the calculated pe. The next line is

eatom[i] += pe;

That's all. I may have added the pe to the wrong energy variable, I just don't know. It was the simplest way I could figure to add it to the atom's energy.

In FixWallGran::post_force(int vflag) I add:

double pe;

double *eatom = force->pair->eatom;

and inside the loop of local atoms,

pe = hertz_history(...)

now returns the calculated pe. The next line is

eatom[i] += pe;

That's all. I may have added the pe to the wrong energy variable, I just
don't know. It was the simplest way I could figure to add it to the atom's
energy.

yes. you tally the fix energy into the per atom *pair* energy. the
"clean" way to do this if you want the per atom energy, would be to
return this into a per_atom vector computation.

otherwise, if you simply want the total, you need to accumulate the
energy and return it as a global scalar computation. see for example
fix_spring.cpp. it is also recommended to have a careful look at fix.h
for the description of the available global flags.

axel.