How to compute the rate of change of a global scalar in LAMMPS?

Hi LAMMPS developers and users,

I want to compute the rate of change (derivative) of a global scalar, i.e., the number of atoms in a group. LAMMPS has fix store/state to store the per-atom attributes, so that one can easily compute the derivatives of those per-atom quantities. However, I have not found such a command for global scalars. So I wonder if there is a convenient way to obtain this derivative on the fly without modifying the source code. Thanks in advance!

You can use fix vector and then the slope() variable function on that vector. Fix vector has an “nmax” keyword that can limit over how many values you will compute the slope. To reduce noise, it is often better to not use only the direct delta, but to compute the slope over a larger number of values.

Thanks for your advice. But what I want to obtain is the instantaneous rates. The slope function assumes that the value changes linearly, which is not the case in my simulation. Now I output it to a yaml file via fix print, and compute the derivatives during post-processing.

If you limit the vector to two elements, the result is the same.

1 Like

Wow, that’s true. I will have a try!

But as I already mentioned, such data is generally very noisy and computing the derivative is enhancing this, thus some averaging is usually beneficial and that you get from increasing the vector length.

1 Like

Yeah, now I understand why you suggested using the Nmax keyword. Perhaps I should specify Nmax as three or four for smoothing.

I use fix evaporate to implement the outflow boundary condition, and want to compute the outflow rate. Becuase fix evaporate computes a global scalar, I try to add it to fix vector:

fix         outflow fluid evaporate 1 10000 outlet 2332
fix         ns_in_out all vector 1 f_outflow nmax 3

But was given the error:

ERROR: Must set ‘extvector’ or ‘extlist’ when setting ‘vector_flag’ for fix vector. Contact the developer. (src/fix.cpp:138)

I don’t understand what I should do according to the error message. But when I store f_outflow in a variable and add the latter to fix vector, the error is fixed:

fix         outflow fluid evaporate 1 10000 outlet 2332
variable    n_outflow equal f_outflow
fix         ns_out all vector 1 v_n_outflow nmax 3

When reporting unexpected errors, please always report which LAMMPS version you are using and what platform you are on. Also, please see the forum guidelines post

LAMMPS version: 29 Aug 2024
Platform: both CentOS 7 and window 10 report this error.
Here’s a minimal working example:

variable npart  equal 300
units    lj
dimension 3
atom_style  atomic
boundary        p p p
neighbor        0.5   bin
region box block -10 10 -10 10 -10 10
create_box 1 box
region top block INF INF INF INF 3 INF
region bot block INF INF INF INF INF 3
create_atoms 1 random ${npart} 324523 bot
mass   1  1
pair_style soft  1.0
pair_coeff * * 10.0

velocity all create 2.0  34234123 dist gaussian
minimize 1e-4 1e-4 1000 1000
reset_timestep 0
fix 1 all nve
fix 2 all evaporate 100 10 top 49892
fix ns_out all vector 1 f_2 nmax 3
run 1

Please try making the following change, recompile, and try again:

  diff --git a/src/fix_vector.cpp b/src/fix_vector.cpp
  index 7c75f93a3a..de8fa83dbe 100644
  --- a/src/fix_vector.cpp
  +++ b/src/fix_vector.cpp
  @@ -114,9 +114,9 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
           error->all(FLERR, "Fix for fix {} vector not computed at compatible time", val.id);
   
         if (val.argindex == 0)
  -        value = ifix->extvector;
  +        value = ifix->extscalar;
         else
  -        value = ifix->extarray;
  +        value = ifix->extvector;
         val.val.f = ifix;
   
       } else if (val.which == ArgInfo::VARIABLE) {

P.S.: If you need an updated Windows installer package for Windows, I can probably produce one as well, but we are planning to release update2 for the stable version next week, and that will include this bugfix and new packages will be provided for that anyway.

1 Like

Thanks for your effort. The error is fixed now after the source code is modified.

You are so kind, but I will try building it on Windows by myself first.

It’s also worth noting that fix controller (documentation) tracks the time derivative of its process variable as an output scalar. Thus, for example, the following combination:

variable dummy internal 0.0
fix      outflow fluid evaporate 1 10000 outlet 2332
fix      ns_out  all   controller 5 1 0 0 1 f_outflow 0.0 dummy
fix      ns_ave  all   ave/time   5 20 100  f_ns_out[3]

should track the difference between f_outflow’s current value and its 5-steps-ago value, every five steps, in f_ns_out[3], and the fix ns_ave will further collect 20 values of ns_out and average over them (albeit somewhat superfluous in this example – sampling a difference every 5 steps, and averaging over 20 samples, is exactly the same as sampling a difference every 100 steps).

You should be able to use this without patching the source code and recompiling.

1 Like

Thank you for providing a new perspective on derivative computation. It is indeed very convenient. I will try it.