Too many time steps in fix/ave for RDF?

LAMMPS LAMMPS Mailing List Mirror
1

Greetings,

I am performing simulations of point-charge mixtures with coul/long forces. For my application, it is important to have “nice” radial distribution functions, which I am trying to achieve by averaging over every single time step as follows

compute gij all rdf 100 1 1 1 2 2 1 2 2
fix gofr all ave/time 1 {nt} {nt} c_gij[*] file rdf.dat mode vector

where ${nt} is the number of time steps in the run. I find that when nt=10000, all goes as expected – the file ‘rdf.dat’ has all my g®’s and they look lovely. However, if I increase nt to 15000 or 20000, leaving all other simulation parameters the same, the resulting ‘rdf.dat’ contains only the header comments – no data. The log file shows no errors or warnings.

Any leads on why I lose RDF output when increasing the number of time steps? Are there internal memory caps in LAMMPS that I might be running up against?

Regards,

Nathaniel Shaffer
University of Iowa
Physics and Astronomy

2

Greetings,

I am performing simulations of point-charge mixtures with coul/long
forces. For my application, it is important to have "nice" radial
distribution functions, which I am trying to achieve by averaging over
every single time step as follows

​this reasoning doesn't make much ​sense. the configurations of subsequent
MD steps are highly correlated and thus have extremely little statistical
relevance and thus do not contribute much to the g(r) computation. for
liquids, there is little improvement doing g(r) computations more often
than every 1000 steps. increasing the total number of samples to reduce the
noise, through doing longer simulations, will indeed improve smoothness and
correctness. How large the statistical noise is, also depends on the width
of intervals over which the g(r) is accumulated. the wider the interval,
the smaller the noise, but also the less precise the result (e.g. in
absolute peak location and height).

compute gij all rdf 100 1 1 1 2 2 1 2 2
fix gofr all ave/time 1 \{nt\} {nt} c_gij[*] file rdf.dat mode vector

where ${nt} is the number of time steps in the run. I find that when
nt=10000, all goes as expected -- the file 'rdf.dat' has all my g(r)'s and
they look lovely. However, if I increase nt to 15000 or 20000, leaving all
other simulation parameters the same, the resulting 'rdf.dat' contains only
the header comments -- no data. The log file shows no errors or warnings.

Any leads on why I lose RDF output when increasing the number of time
steps? Are there internal memory caps in LAMMPS that I might be running up
against?

​this could happen when the value of ${nt} is larger than the number of
time steps in the "run" command.

​axel.​

3

hello, developer,

I want to know the meaning of mask of atom and the reason for using the keyword. What is the difference with type?

4

hello, developer,

I want to know the meaning of mask of atom and the reason for using the
keyword. What is the difference with type?

​the Atom::mask variable contains the information about which group(s) an
atom belongs to.​

​axel.​

5

Hi Axel,

Thanks for the reply.

        Greetings,
        
        I am performing simulations of point-charge mixtures with
        coul/long forces. For my application, it is important to have
        "nice" radial distribution functions, which I am trying to
        achieve by averaging over every single time step as follows
        
​this reasoning doesn't make much ​sense. the configurations of
subsequent MD steps are highly correlated and thus have extremely
little statistical relevance and thus do not contribute much to the
g(r) computation. for liquids, there is little improvement doing g(r)
computations more often than every 1000 steps. increasing the total
number of samples to reduce the noise, through doing longer
simulations, will indeed improve smoothness and correctness. How large
the statistical noise is, also depends on the width of intervals over
which the g(r) is accumulated. the wider the interval, the smaller the
noise, but also the less precise the result (e.g. in absolute peak
location and height).

Sure, point taken. For this discussion, let's say I'm more interested in
the technical question of how the averaging fix is being applied.

        Any leads on why I lose RDF output when increasing the number
        of time steps? Are there internal memory caps in LAMMPS that I
        might be running up against?
        
​this could happen when the value of ${nt} is larger than the number
of time steps in the "run" command.

I should have mentioned that I am running for ${nt} time steps, so that
the RDF-related input lines are, together:

variable nt equal 10000 # or 15000 or 20000
...
compute gij all rdf 100 1 1 1 2 2 1 2 2
fix gofr all ave/time 1 \{nt\} {nt} c_gij[*] file rdf.dat mode vector
...
run ${nt}

so I think I've eliminated the possibility of a mismatch between the
averaging time and the run time. To be safe, I also tried using 'run
$(v_nt+1)' to guard against possible off-by-one counting errors. The
fact that nt=10000 produces RDF output, while nt=15000 or nt=20000 do
not, is the basis for my hunch that it's a memory issue.

On the side, I also tried the following, taking to heart your comment on
not over-sampling:

variable nt equal 100000
...
compute gij all rdf 100 1 1 1 2 2 1 2 2
...
fix gofr all ave/time 1000 \(v\_nt/1000\) {nt} c_gij[*] file rdf.dat mode
vector
run $(v_nt+1)

which by my reading of the docs amounts to: "Run for 100001 time steps,
collecting an averaged g(r) on the 100000th time step based on 100
samples taken every 1000 steps." This script also yielded no RDF output.
This throws a wrench into my "memory limit" hypothesis, since it has the
fewest samples of any of my attempts. It makes me think I'm
misunderstanding the ave/time fix at some very basic level.

-- Nathaniel

6

i have no problems for any value of ${nt} when i use the following modified version of the melt input.

lattice fcc 0.8442

region box block 0 10 0 10 0 10
create_box 2 box
create_atoms 1 box
region half block 0 10 0 10 EDGE 4.9
mass * 1.0
set region half type 2

velocity all create 3.0 87287

pair_style lj/cut 2.5

pair_coeff * * 1.0 1.0
neighbor 0.3 bin

neigh_modify every 5 delay 5 check yes

compute gij all rdf 100 1 1 1 2 2 1 2 2
fix gofr all ave/time 1 {nt} {nt} c_gij[*] file rdf.dat mode vector

fix 1 all nve

thermo 500
run ${nt}

7

Hi Axel,

Thanks for sticking we me on this.

i have no problems for any value of ${nt} when i use the following
modified version of the melt input.

lattice fcc 0.8442

region box block 0 10 0 10 0 10
create_box 2 box
create_atoms 1 box
region half block 0 10 0 10 EDGE 4.9
mass * 1.0
set region half type 2

velocity all create 3.0 87287

pair_style lj/cut 2.5

pair_coeff * * 1.0 1.0
neighbor 0.3 bin

neigh_modify every 5 delay 5 check yes

compute gij all rdf 100 1 1 1 2 2 1 2 2
fix gofr all ave/time 1 \{nt\} {nt} c_gij[*] file rdf.dat mode vector

fix 1 all nve
thermo 500
run ${nt}

​however, there i will get the header-only rdf.dat file you are
getting, if i insert:

run 100

*before* defining the compute and the fix ave/time.

this is because with the additional 100 steps, the fix ​ave/time will
start collecting data only on step \{nt\}, that is the {nt}-100 th
step in the second run.

so i suspect, you left out the crucial information, that you are not
starting your production simulation at step 0, but some other step.
and then because of the mismatch, averaging will start much later, as
it is based on the *time step number* and not the time step count in
the current run.

You nailed it -- I'm running off a restart file read in at the top of
the script. Didn't read the ave/time docs closely enough to realize
that it operates on the absolute time step.

to prove this, all you need to do is to replace

run ${nt}

with

run \{nt\} run {nt}

​
and you should have some output in your rdf.dat

however, the clean solution would be to either put the command

reset_timestep 0

at some strategic location before the production run command.

Yep, both of these check out.

For any curious future readers, here's why I was getting the behavior I
was getting:

I was performing my runs by loading a restart from a 50000-time step
equilibration. When I run and average with 10000 time steps, I do get
output because of luck -- the starting time step (50000) is a multiple
of 10000, so I have a full averaging period (50000-60000). However,
with 15000 time steps of running/averaging, the ave/time fix doesn't
really start until step 60000 (the first multiple of 15000 after the
restart), and then when the run ends on step 65000, it hasn't completed
a full averaging period. Using 20000 time steps of running and
averaging fails for the same reason, but 25000 succeeds because 25000
is a multiple of 50000.

Cheers and thank you,
Nathaniel