Mechanism for saving/returning to previous state without writing to disk

LAMMPS LAMMPS Development
1

TLDR:

Is there a way to create a copy of the lammps C++ or Python Object? Can we deepcopy a Lammpy Object in Python (or is there a copy constructor in C++ bindings)
If not, can we reset the system back to a previous state without writing to file, e.g. by setting the positions and velocities through the Python interface?

More details
We are trying to implement a method using the Lammps Python interface which takes multil child trajectory off of a parent trajectory at intervals. The parent and child systems are identical in terms of positions, velocities, bonds, etc at the start (a mapping and force is then applied to the child). Here is an image that show the idea:

children
children1772×1078 41.5 KB

As this should run in parallel on a supercomputer (mpi4py) with lots of short runs, ideally we want to avoid writing to disk if possible. The most elegant way to do this would be to create a single parent object, evolve it in time and then copy/clone to make children runs at given intervals, run these child systems to get output and deallocate them, before running the parent forward again. Is this possible? I get an error with copy.deepcopy (which I wouldn’t expect to work for a wrapped C code anyway) and looking into the Lammps library.cpp C++ source code it seems there is no copy constructor (I assume this is not trivial to create for a general Lammps run)?

If not, is there a way to save the state of the system, e.g. the following psudo-code is a minimal example of what I want to do,

#Get the array of initial positions and velocity
nlmp = lmp.numpy
ids0 = np.copy(lmp.extract_atom("id"))
r0 = np.copy(lmp.extract_atom("x"))
v0 = np.copy(lmp.extract_atom("v"))
#Run a simulation forward
lmp.command("run 1000")
#Get pointers to variables
idsp = lmp.extract_atom("id")
rp = lmp.extract_atom("x")
vp = lmp.extract_atom("v")
#Set variables back to original values
for n in range(ids0.shape[0]):
    for ixyz in range(3):
        rp[idsp[n]][ixyz] = r0[ids0[n],ixyz]
        vp[idsp[n]][ixyz] = v0[ids0[n],ixyz]

#This run should then be identical to the one above
lmp.command("run 1000")

When I try this, I get either errors or different output for the two runs (depending on indexing choices). Do additional variables need to be reset here or is there are better way of doing this?

I’m aware that another option to acheive the above would be,

#Get the array of initial positions and velocity
#Run a simulation forward
lmp.command("write_data file.dat ")
lmp.command("run 1000")

#This run should then be identical to the one above
lmp.command("read_data file.out")
lmp.command("run 1000")

but ideally we want to avoid the overhead of writing to disk, especially for big systems. Another option here is if write_data could be piped to a string variable or stdout.

2

There is no copy constructor for the LAMMPS class, that would be too much of a problem with some of the data stored within (MPI communicator, file descriptors, threads), but you can use write_restart instead of write_data. This is much faster, much more compact, and stores far more data. With read_restart you can recover most of your system settings.

3

Hi akohlmey, thanks for the quick reply and for confirming no copy is possible. Looks like write/read_restart may be the way to go.
Would you expect setting positions and velocities back to previous values to work or too fragile an approach for general systems?

4

This is going to be tricky business the way you are doing it, because the “extract” functions work on atoms that are present as local and ghost on a given MPI rank only. So you have to make certain that you have the exact same domain decomposition which is essentially impossible. You would have to use “scatter/gather” functions to collect to and broadcast from the MPI root process. But at that point, using read/write restart is effectively the same. On most modern Linux system, you can write to a ramdisk at /tmp/ (make sure to delete after reading, so you reclaim the RAM) and then there is no more overhead than the MPI communication.

5

What about doing this from within LAMMPS? Store initial coordinates and velocities with fix store/state and after the end of a simulation restore it with set command.

6

That’s interesting, I didn’t know that was possible to use a ramdisk (I was looking at overloading the file object passed to lammps to achieve something similar).

For the proof of concept we are instantiating each lammps instance using MPI_COMM_SELF so all natom and its ghost copies would be on that processor. Appreciate this will not work in parallel without scatter/gather operations (which the write/read from disk essentially take care of).

I think a write to disk solution above is a more portable stratergy and if we suggest using the ramdisk if possible, the performance penalty should be minor.

Thanks again for your help!

7

Even in the serial case, you may have a different number of ghost atoms and local atoms may be at different positions and thus get memory corruption or a segfault. You would have to enforce a neighborlist rebuilt to have consistent data. This is - of course - enforced on the restart, but not automatic if you store (local) data at some point in the time integration. Between reneighboring, atoms may move outside of the box.

Also, re-neighboring may change the order of atoms (and atom sorting), so you really need the gather/scatter mechanism and use atom IDs (tags in the source code) to make sure you are restoring data to the correct atoms. I recently wrote a fix where two partitions have to be kept in perfect sync, so I can use mixed forces and that was really tricky to get right.

8

I didn’t know about fix store/state, interesting, I’ll look into it further. I suspect it will be susceptible to the same limitations as extracting the positions at a given time…

9

I thought by setting the atoms using the ids directly (tags) i.e.
rp[idsp[n]][0] = r0[ids0[n],0]
I would avoid these problems with ordering/sorting on a single process. This doesn’t solve the ghost copies which I assume change each timestep. If I limited the copies to real atoms only, then set only these real atoms on restart and forced a regeneration of the ghost/halos + rebuilt the lists then maybe it would work. I assume by the time I’ve done all this, maybe best to just read a binary input file.

10

The local atom indexing is not the same as atom IDs/tags.
From within the C++ code you would need to use the atom->map(ID) to get the corresponding index in the local atom->x, atom->vx arrays. When running in parallel you would also get indices to ghost atoms (index >= atom->nlocal) and -1 for atoms that are not present as local or ghost. For small systems, the same atom ID may be present multiple times (as ghost).

11

You would have to use either “extract” or “gather” in the python interface to get the data, or use atom style variables within LAMMPS input script and then the set command. The necessary mapping and communication will be done by LAMMPS internally in this case.

My suggestion is to use “clear” and then “read_restart” as initial implementation (possibly writing/reading to/from a ramdisk) and then observe if this has a significant performance impact.

This could just turn out to be a problem of “premature optimization™”.

12

Ah, sorry. I misunderstood your problem.

13

‘This could just turn out to be a problem of “premature optimization"’

Very possibly. I originally thought it might be simpler to just set these values directly rather than writing to disk.

It also potentially allowed some nice division of labour in parallel, with a parent on its own process sending updated ‘r/v’ values to each child. Much more scalable if no file system is involved.

Thank you again for your help.

14

BTW, the fastest way to broadcast a binary file over MPI is to open it read-only on the root rank, then mmap() it and just broadcast the resulting buffer as a char buffer to any other rank that needs it and and there just use a low-level write or a copy-on-write mmap(). That will minimize unnecessary buffering.