Geometric Dumping

Generating geometrically spaced trajectory data with geo-rev11.py and geometricdump.com, and analyzing it in the cyclic modes of 2_time.py and s4_tb_fft_ww.py.

← / → or space · d toggles dark mode · click the slide counter to jump

The problem

The linear modes of 2_time.py and s4_tb_fft_ww.py work on linearly spaced dumps. That is fine until the dynamics spread over many decades: a correlation function can only be evaluated at time separations that exist in the saved trajectory, and a single linear spacing cannot cover all the decades at once.

Numbers for the production region used throughout these slides (106 steps ≈ 100 τα, 100k particles, ~1.2 MB per frame):

Dump scheduleFrames / runSize / runLag coverage
linear, every step106~1.2 TBeverything
linear, every 1000 steps1,0001.2 GBnothing below 1 LJ; the first 3 decades are missing
geometric (25 dumps/cycle)2,7753.3 GB1 step to full span, uniform in log t

Why repeated short sequences

A single geometric series from step 1 to the end of the run would solve storage. It fails for a different reason. Correlation functions are averaged over initial times t0: the same separation t is measured from many starting frames. The translational overlap, for example:

wj(t0+t, t0; a) ≡ Θ(a − |rj(t0+t) − rj(t0)|)
Fo(t) ≡ (1/Nt₀) Nt₀Σt₀=1 (1/N) NΣj=1 ⟨wj(t0+t, t0; a)⟩

The schedule is therefore built in two levels:

Cost: separations between the cycle length Λ and its multiples are coarser. In practice the set of available separations stays dense enough (slide 13 shows it exactly).

The schedule

Generated live with the geo-rev11.py --small algorithm (geometric runtimes). Each tick is one frame. The time axis is fixed at 160 LJ, so changing -t or the cycle count changes the coverage you see.

Top: the dumps in real (linear) time. Bottom: one row per cycle, each showing the time of its dumps measured from that cycle's own start, on a log axis. The rows are identical: every cycle repeats the same pattern.

geo-rev11.py

Prints three text files that a LAMMPS loop reads line by line. One frame per DCD file.

FileContents
run_times.txtsteps to run before each dump (argument to run)
dump_times.txtcumulative timestep of each dump (dump trigger / bookkeeping)
file_names.txtoutput filename for each frame
OptionMeaning
-a [float](LJ time units) first geometric term. Without --old, as used here, this is the first runtime; with --old it would be the first dump offset. Values below one timestep floor to one timestep.
-t [float](LJ time units) endpoint used with -a and -n to set the geometric ratio. Without --old, as used here, it is the largest single runtime and the span is whatever the runtimes sum to; the site solves this backwards from the span you ask for.
-n [int](count) number of geometric terms. Without --old every runtime is at least one timestep, so no term is dropped and this is exactly the dump count.
-d [float](LJ time units) simulation timestep, used to convert LJ times to integer timesteps.
-r [int](count) number of times to repeat the generated sequence back to back.
-c [int](file index) initial filename counter before each generated file increments it. Set this to the last already-used index; the first generated file is -c + 1.
-m [int](timesteps) offset added to every cumulative dump time, used to continue the timestep counter after a previous sequence.
--small(mode flag) generate the simple geometric-series schedule used by this workflow. Without it, geo-rev11.py uses the older long-series mode.
--old(mode flag) treat the geometric terms as cumulative dump offsets and compute runtimes as their differences. This workflow does not use it: the runtimes are the geometric terms instead.
--test(diagnostic flag) print diagnostics instead of loop-file output. Do not use this output as LAMMPS schedule files.
Output files are opened in append mode. Delete them before re-running, or two schedules get concatenated.

The LAMMPS loop

Choose the version that matches how the job will run. The ZIP includes both files for every long+short schedule.

Legacy: use this compact loop for a fresh run that is expected to finish without resuming from a checkpoint. It resets the timestep and does not check whether trajectory files already exist.
# --- Looping Parameters Setup --- # variable loopvar loop 3366 variable eruntime file run_times.txt variable filename file file_names.txt variable dumpint file dump_times.txt reset_timestep 0 dump 1 all dcd 2 traj1.dcd dump_modify 1 unwrap yes run 1 undump 1
  • The initial one-step run writes the timestep-zero frame to traj1.dcd.
  • The three file-style variables then consume one line per trajectory file.
  • This version is intentionally unchanged from the original site workflow.
label loopstart dump 1 all dcd ${dumpint} ${filename} dump_modify 1 unwrap yes run ${eruntime} undump 1 # next loopvar next eruntime next filename next dumpint jump SELF loopstart
  • Each iteration advances by the runtime in run_times.txt and writes one DCD.
  • There is no restart catch-up logic. Starting this loop again at timestep zero can overwrite earlier data.

default vs --old: what is geometric?

default: geometric runtimes

runtimek = ⌊t0 rk⌋; dump time = running sum. -t sets the largest single gap. The sequence ends wherever the sum lands, so the site solves -t backwards from the span you ask for. No runtime can floor to zero, so -n is exactly the dump count.

Used for both sequences. Repeating the same runtime list still reproduces identical offsets in every cycle, which is what the analysis scripts verify.

--old: geometric dump times

dump offsetk = ⌊t0 rk⌋ exactly; runtime = difference from the previous dump. The sequence ends exactly at -t, but early terms can floor onto each other and disappear.

Not used here. Shown for contrast: it is the alternative reading of the same geometric series.

Top: 27 dumps spanning 9 LJ with geometric runtimes, one sequence. Bottom: -a 0.001 -t 9 -n 27 --old -r 3, on the same linear axis. Both production commands use --small; the non---small path is the legacy long-series mode and predates the current workflow.

Two ways to use the schedule

A: short sequences only

System is equilibrated before recording. Every cycle start is an equivalent initial time. Analyze with 2_time.py -G … and s4_tb_fft_ww.py -G ….

B: long sequence, then short sequences

During aging, correlations depend on the waiting time tw, and the evolution is slow in log tw, so one geometric series samples it with a few hundred frames. Analyze the long region with 2_time_aging.py (explicit initial frames); analyze the short region in the unified scripts’ cyclic modes, starting at the right file index (-z).

Don't point either script’s cyclic mode at the stitched long+short set: the long region doesn't repeat, and the cycle check fails (correctly).

geometricdump.com

You can hand-pick geo-rev11.py flags, but they interact: the ratio depends on -a, -t, -n, and the stitch values -c and -m depend on the long sequence's output. The site takes system-level inputs and derives the flags.

You provide
  • dt1 (long) and dt2 (short)
  • τα
  • τequil
  • expressions like 100*tau_alpha work
It derives
  • long: runtimes geometric, ratio fixed by matching the largest gap and summing to the seam, -n from the first run time, companions added after
  • short: length 0.9 τα, -n 27, -r = total/length
  • stitch: -c, -m from the long output (+ the step-0 file)
You get
  • both commands with per-flag explanations
  • consistency checks
  • DCD size estimate
  • plots
  • downloadable schedule files
  • LAMMPS loop snippet

The reason for the specific defaults is dump density at τα: next slide.

Why cycle length ≈ 0.9 τα

A schedule can look reasonable and still measure τα badly. The available separations get sparse toward the end of a cycle and dense again right after the cycle boundary. The cycle length decides which region τα falls in.

τα = 10 LJ, dt = 0.001, -n 27

At the default Λ = 0.9 τα, τα sits just past the cycle boundary, in dense spacing. Push Λ above τα and τα lands in the sparse tail of the cycle, where the nearest available separations are far apart.

Worked example

dt1 = dt2 = 0.001, τα = 10 LJ, τequil = 200 LJ, production 100 τα, companions on.



python geo-rev11.py -a 0.001 -t 2.32607519009 -n 665 -d 0.001 -c 1 --small -r 1 -m 0

Long sequence: one sequence whose runtimes are geometric, ending at the seam so the largest gap matches the short sequence's final 2.66 LJ gap, with 665 terms (no runtime can floor to zero, so all 665 become schedule dumps), first dump one timestep in, filename counter starting at 1 (the step-0 dump took traj1.dcd), timestep counter starting at 0. The site then adds companions at +3 and +5 timesteps, bringing 665 dumps to 1,801:



python geo-rev11.py -a 0.001 -t 2.32413132057 -n 27 -d 0.001 -c 1802 --small -r 111 -m 200000

Short sequences: 111 repeats of a sequence whose runtimes are geometric, each spanning 9 LJ = 0.9 τα with a goal of 27 dumps and the first dump 1 step in, continuing the filename counter from the long sequence and the timestep counter from its last dump:




Long (with companions)Short (production)
dumps665 schedule + 1,136 companions = 1,801111 cycles × 27 = 2,997
filestraj2.dcdtraj1802.dcdshort1803.dcdshort4799.dcd
timesteps1 – 200,000 (= 200 LJ)200,001 – 1,199,000 (= 1,000 LJ past the seam)

One short sequence (production cycle), in steps. Dump times, measured from the start of the sequence (the last is 9 LJ / 0.001 = 9,000). Offsets 4 and 5 are the default companion separations:

1 2 4 5 8 11 16 23 33 47 66 94 134 191 271 385 546 775 1100 1562 2217 3147 4467 6340 9000

and the corresponding runtimes (the gaps between dumps):

1 1 2 1 3 3 5 7 10 14 19 28 40 57 80 114 161 229 325 462 655 930 1320 1873 2660

Sizing the long sequence

The long sequence has one job: sample waiting time after the quench. Two things about it are already decided before you choose anything.

The seam is fixed

The last long dump lands on τequil, where production starts. Call it S in timesteps.

The coarsest spacing is fixed

The widest gap in the long sequence is matched to the widest gap of a production cycle, G, so the sampling density does not jump at the seam.

Those two facts pin the ratio completely. The last gap of a geometric sequence ending at S is S(1 − 1/r), so setting it equal to G gives

r = S / (S − G)

Nothing is left to tune there. The only free choice is a1, the first dump time, and once you pick it the number of terms follows:

terms = ln(S / a1) / ln(r)

So a1 is the whole knob. It is the earliest waiting time you will ever be able to report, and it is also what sets the file count. The default is one timestep, the earliest possible.

Starting later saves less than you would expect. Once adjacent terms are closer together than one timestep, which happens below a1 = 1/(r − 1), they floor onto the same integer and collapse into a run of consecutive timesteps. Those dumps were going to be taken anyway, so pushing a1 right mostly discards early waiting times without buying much storage back.

The requested term count and the number of files are not the same number when flooring is in play. The site reports both, and the second one is what the analysis scripts need.

A value, or a rate?

One frame at waiting time T, paired with a later frame, gives a correlation value. But some quantities are correlations of the relaxation rate: how fast the correlation is falling at T. A rate is a difference, so it needs two frames close together at T, not one.

The production region gets this for free: its cycles repeat, so short separations exist everywhere. The aging region does not repeat, so if the schedule only records single frames, no rate can ever be recovered from it. The fix is to add a companion dump a few timesteps after each scheduled dump.

worked example: τα = 10 LJ, τequil = 200 LJ, dt = 0.001

A rate at waiting time T needs a main dump δ0 plus nearby companions δ1, δ2; the next main dump is δ0′. Top: t < tx (companions land on scheduled dumps). Middle: t > tx (companions are extras in the gap). Bottom: gap versus waiting time with tx.

The separations are taken from the production cycle's own dump offsets, indexed by i. That is not decoration: it guarantees the same physical separation exists on both sides of the seam, so a rate means the same thing in the aging region and in production. It is also why the companions that would fall past the seam can simply be dropped, since the first production dumps already supply them.

Watch the plot as you move i. The red line is the crossover tx: left of it the schedule is already denser than the companion separation, so δ0′ arrives at or before the widest companion and those frames land on dumps that were being taken anyway; right of it every dump owns its own companions in the gap. Raising i pushes tx later; raising a1 throws away early waiting times.

When dt1 and dt2 differ, a separation must be a whole number of timesteps on both grids. If dt1/dt2 is p/q in lowest terms, only offsets divisible by p survive, and the site advances i to the first pair that does.

What cyclic mode does before computing

The analysis scripts don't read dump_times.txt directly. The schedule is recovered from the DCD metadata. The cyclic modes of 2_time.py and s4_tb_fft_ww.py both start with these steps (L = frames per cycle, the -G argument):

  1. Measure the cycle. The times of the first L+1 frames give the L time gaps that make up one cycle.
  2. Verify it. Every frame time in the whole trajectory has to match that repeating pattern. If one frame doesn't, the script stops and reports which frame in which file.
  3. Rotate. The starting point within the cycle is shifted so analysis begins at the smallest gap.
  4. List the available separations. Adding the repeating gaps end to end gives every time separation that exists in the data. Requested lags are then moved to the closest available separation (compared on a log scale).

Step 2 is why --small matters at generation time, and why a wrong -G or -z fails immediately instead of producing wrong numbers.

Lag selection: -g and -U

-g N asks for N lag values, geometrically spaced between a first lag and the longest separation in the data. The data only contains specific separations, so each requested value is moved to the closest one available. When two requests land on the same separation it is kept once, which is why fewer values come out than went in.

-U sets the first requested lag. After the cycle is rotated to begin with its shortest gap, -U n (positive) starts at n times that gap. -U 0, the old default, starts at 1.0 LJ, so no separations shorter than 1 LJ get requested. Try it below.



worked-example data: 2,775 frames, 25-frame cycle (actual progression.py algorithm)

2_time.py: linear vs cyclic mode

One script now handles both frame layouts. Supplying -G selects cyclic mode; omitting it selects linear mode.

linear: omit -G

  • defaults to traj<n>.dcd; count -n, prefix -N
  • evenly spaced dumps
  • initial times every -d frames
  • lags used directly as frame counts

shared

  • MSD, overlap (-a), Fs (-q), orientational correlations (--polyatomic*)
  • run averages and across-run standard deviations (-r runs)
  • frame range -s/-m, last initial time -k
  • lag request flags -f / -g / -w, --low-interval / --high-interval
  • -P run-level workers, normally equal to -r
  • output column 1 is the real lag time, so results plot together

cyclic: give -G

  • defaults to short<n>.dcd; count -n, prefix -N, first file -z
  • requires -G, the true frames per cycle
  • verifies the cycle from DCD metadata; file tests first*
  • initial times at cycle starts
  • requested lags moved to available separations; -U


python 2_time.py -r 5 -P 5 -n 2997 -N short -z 1803 -G 27 -g 50 -U 1 > out.txt

For the worked example: average over 5 runs with 5 run-level workers, reading 2,997 short files per run starting at short1803.dcd, with 25 frames per cycle, requesting 50 geometric lags with the first lag one frame in:

* file tests: before computing, the script checks that every expected DCD file exists and opens readably (--skip-file-tests skips this part; the cycle verification always runs). Cyclic mode needs at least 2 full cycles of frames. --negvals is forced off in cyclic mode.

S4 analysis

S4(q, tb): spatial correlations of mobility. For each interval length tb, particles get a mobility weight w(t0, t0+tb), the weights are binned on an FFT grid, and the structure factor of that field is computed. Output per tb: total, self, and distinct parts, with errors from across-run scatter. At least 2 runs required.

s4_tb_fft_ww.py: linear vs cyclic mode

One script now handles both layouts. Supplying -G selects cyclic mode and adds the cycle steps from slide 12; omitting it selects linear mode.

linear: omit -G

  • count -n; use -N traj for standard traj<n>.dcd files
  • tb values in frame units
  • -d initial-time spacing in frames (default 10)
  • -M margin in frames

shared

  • S4 total / self / distinct with across-run errors
  • -x FFT grid, -y box size (both required)
  • ta (--ta), tc (--tc), usually 0; W1/W2 flags
  • qshell flags -q / -v / -l
  • -k last initial time; -i one output file per tb; -r ≥ 2
  • -P run-level workers, normally equal to -r

cyclic: give -G

  • prefix -N (default short), count -n, first file -z
  • requires -G, the true frames per cycle
  • tb values moved to available separations; -U, --geom-terms-factor
  • -d in cycles (default 1); -M in cycles
  • realized (ta, tb, tc) at each initial time must match the first occurrence within -S; mismatches are skipped


python s4_tb_fft_ww.py -r 5 -P 5 -n 2997 -N short -z 1803 -G 27 -g 30 -x 24 -y 47.0 --W1-theta 0.25

For the worked example: average over 5 runs with 5 run-level workers, reading 2,775 files named short*.dcd starting at index 593, with 25 frames per cycle, 30 requested tb values, a 24³ FFT grid for a box of side 47.0, and a Heaviside mobility function (translational overlap) with threshold radius 0.25:

Things to be careful about

IssueBehavior / fix
re-running geo-rev11 without deleting outputsschedules concatenate silently (append mode)
set_len is the dumps per cyclewith geometric runtimes no term is dropped, so -n 27 gives 27 dumps and -c 27 / -G 27. The wrong value fails the cycle check
forgetting the step-0 dumpit consumes the first filename; first generated file is index 2, and the first production file in the example is short1803.dcd (1 + 1,801 + 1)
wrong -z in long+short datalong-sequence frames enter the cycle check and it fails
--skip-file-testsskips only the existence/readability tests; the cycle check always runs
-U < 0unsupported; the script warns that this path is incorrect, so leave -U nonnegative
S4 with 1 runhard error; errors come from across-run scatter
ta or tc ≠ 0 without --W2-*hard error
too few framescyclic 2_time.py needs at least 2 cycles; cyclic S4 warns if -M leaves too little room and may fail at the last tb

Command quick reference

The worked example end to end: dt1 = dt2 = 0.001, τα = 10 LJ, τequil = 200 LJ, production 100 τα, 5 runs.

# 1. schedule files (or use geometricdump.com) python geo-rev11.py -a 0.001 -t 2.32607519009 -n 665 -d 0.001 -c 1 --small -r 1 -m 0 # companions at +3 and +5 are added by the site, not by geo-rev11 python geo-rev11.py -a 0.001 -t 2.32413132057 -n 27 -d 0.001 -c 1802 --small -r 111 -m 200000
# 2. LAMMPS: loop over run_times.txt / dump_times.txt / file_names.txt # use the loop snippet from slide 6
# 3. two-time correlations on the production region python 2_time.py -r 5 -P 5 -n 2997 -N short -z 1803 -G 27 -g 50 -U 1 > out.txt
# 4. S4 on the production region python s4_tb_fft_ww.py -r 5 -P 5 -n 2997 -z 1803 -N short -G 27 -g 30 -x 24 -y 47.0 --W1-theta 0.25

-G takes the true dumps per sequence (25), not the generation -n goal (27). -z 1803 = 1 step-0 file + 1,801 long files (with companions) + 1.

Summary

geo-rev11.pygenerates the schedule files (run_times, dump_times, file_names) for the LAMMPS loop
geometricdump.comderives the geo-rev11 flags from dt1/dt2, τα, τequil; emits both commands and the files
2_time.pyMSD, overlap, Fs, orientational correlations; -G selects cyclic dumps
s4_tb_fft_ww.pyS4(q, tb); -G selects cyclic dumps
2_time_aging.pytwo-time correlations vs (tw, tw+t) for the long-sequence region

Repo docs: files_readme/, geometric-dumping/README.md.

slide 1