Looking forward to your reply(some questions)!

[pppyk]

Checklist

I ran your code, but there are still some questions, I hope you can help answer them, thank you.

❓ Question

1. Is the data in Figure 2 a rapid heating process?

2. Can you share the data in Figure 2?

3. For the last lindemann value output in the program, why is it divided by nframes and not divided by n-1

📎 Additional context

[github-actions[bot]]

Hello @pppyk, thank you for your interest in our work!

If this is a bug report, please provide screenshots and minimum viable code to reproduce your issue, otherwise we can not help you.

[Cloudac7]

It could be shown that .

Well, from definition of Lindemann index, the upper term below the root symbol (i.e. np.devide(array_var, nframes)) is just a term yielding the variance of r_ij through the nframes according to Welford Algorithm, thus should be nframes instead of N-1 in the original equation.

However, there seems a bug in the original code yielding Lindemann index per frame, where nframes should be substituted into iframe - 1 (i.e. number of frames from the first frame to current frame), for it would give wrong variance for each frame except the first (of course, 0) and the last one for it ignored the change of timescale. @N720720

# lindemann.index.per_frame
def calculate(frames: npt.NDArray[np.float32]) -> npt.NDArray[np.float32]:
    for q, coords in enumerate(frames):
        ......
        if first:
            lindemann_indices = 0
            first = False
        else:
            np.fill_diagonal(array_mean, 1)
            lindemann_indices = np.zeros((natoms), dtype=dt)  # type: ignore[assignment]
            lindemann_indices = np.divide(np.sqrt(np.divide(array_var, nframes)), array_mean)  # type: ignore[assignment]
            lindemann_indices = np.mean(
                np.asarray([np.mean(lin[lin != 0]) for lin in lindemann_indices])  # type: ignore[attr-defined]
            )

        lindex_array[q] = lindemann_indices
    return lindex_array

[N720720]

Thanks for your contribution. I will look into it.

[N720720]

I added unit tests to reproduce this false behavior

[N720720]

Checklist

I ran your code, but there are still some questions, I hope you can help answer them, thank you.

Question

1. Is the data in Figure 2 a rapid heating process?
2. Can you share the data in Figure 2?
3. For the last lindemann value output in the program, why is it divided by nframes and not divided by n-1

Additional context

1. Yes, we are looking at a heating ramp with 50 thousand steps. However, due care must be exercised, when one interprets such simulated phase transitions. Since 50 thousand simulation steps can lead to the impression that in fact we are looking at a slow process. Quite the opposite is the case since one simulation step translates to 2 femtoseconds. Thus 50 thousand time steps portray only a “real” time of 0.1 nanoseconds (given our timestep length). For a heating ramp that goes from 400K to 1200K this translates to a heating rate of 8 · 10¹² K/s for 50 thousand steps. The heating rate is probably unrealistically fast, only comparable to a short intense laser impulse hitting a cluster, but the constraint lies here in the computational time and cost. To get a heating rate of 4 kelvin per second for a heating ramp difference of 800 kelvin 100 quadrillion (1 · 10¹⁷ ) timesteps would be needed. With the actual simulation setup used, around 8 simulation steps per second where calculated, this would take close to 4 · 10⁸ years to calculate.

2. There are two 50k heating ramps of 459 Atoms clusters in the folder for the unit tests.

3. I fixed this in version 0.6.0