Is there a clean way to combine two simulations' atoms in a third simulation?

LAMMPS LAMMPS Mailing List Mirror
1

Hi! I’m trying to simulate two initial states, one per simulation, and then combine them both in a third simulation where one simulation would go from (0,0,0) to (x,y,z) for example, and the other one would go from (0,0,z) to (x,y,z2), if that makes sense, as if placing one simulation box on top of the other one.

I’ve tried different combinations of the write/read data and write/read restart functions, but as the dump files have clearly shown, between the moment the atom mesh is created and the moment those files are written, the box varies in volume because of its thermostating, making what I want to do quite difficult.

The best solution I’ve found is to generate one box entirely, evaluate the [min,max] of its size and then using an appropriate fix deform on the second box to make it fit in the x and y boundaries (for the given example above), but even then, to make sure the vertical component of the box coincides with the first one without creating a big gap between the two boxes or straining it too much, keeping it from reaching equilibrium with its own thermostating, is quite difficult.

I’m guessing the cleanest way to do so would be to use a third party program or Python/C++ script that would ensure both boxes are in accordance with one another, generates a data file and then read by LAMMPS in the final simulation, but since that solution’s quite advanced, I wanted to ask if there was a way in LAMMPS I just hadn’t thought of in the first place!

I’ve also tried simply creating both boxes in the same simulation, segregating them in different regions and groups and thermostating both with different fixes (one has to have initially twice the temperature of the other one), but this resulted in both boxes reaching a somewhat intermediate temperature, which I figure is caused simply by the fact it would otherwise create an unphysical system because of the step in temperature through the mesh.

Thanks a lot for your input!
Have a good day,
Antoine.

2

Hi Antoine -

If I’m understanding you correctly, I think that the fixedpoint option in fix npt would solve your problem (assuming you meant barostating instead of thermostating in your reply. Fixedpoint allows you to specify each box where to do the dilation for the volume changes - you could set it for the upper edge of your cell (where your boxes connect) for the first box, and the lower edge of the cell for the second box so that you didn’t change the z location of where the boxes connect.

That being said, another issue that you are probably going to encounter as well is the fact that your x,y directions will likely reach different sizes if you equilibrate them differently. You might need to look at doing something like NPzAT (where you only adjust the box volume in the z direction during NPT moves to do this.

Without more examples/scripts it is hard for me to be more detailed, but hopefully something like this might work!

Zeke

3

Hi! I’m trying to simulate two initial states, one per simulation, and then combine them both in a third simulation where one simulation would go from (0,0,0) to (x,y,z) for example, and the other one would go from (0,0,z) to (x,y,z2), if that makes sense, as if placing one simulation box on top of the other one.

there are different ways. it all depends on the details.
there are effectively two main issues that need to be dealt with: force field and geometry.

some comments on the force field issue:

  • it is best to handle this separately and thus write out data file without force field parameters included and instead write out the force field parameters as a separate file with write_coeffs.
  • from there on it depends on the kind of force field and identity of atom types. if they are identical, you can just import one of the write_coeff files and be done. if they are completely different, you need to figure out a way to get the “mixed” parameters. if you use different pair styles, you need to use pair style hybrid and a third pair style to define the interactions between the two substyles. you would also have to renumber the atom types (and bond types etc.) of the second system (can be one on import of the data file).
  • if there is an overlap, things are the most complex. best to treat it as if those are two entirely different systems and just use the same parameters multiple times.

as for the system merging. a lot depends on whether you have periodic boundaries or not. if you have a fully periodic system, you will cut interactions apart and where you connect the two systems, you will have incompatible positions and a system that is equilibrated to a different state. best to use the change_box command to enlarge the box but do it with an empty group of atoms, so the atom positions are not scaled and thus you add a small space.

if both boxes have the same dimensions in x and y, you can easily combine the two systems. if they don’t, you should use the change_box command (but this time with the group all to scale atom positions with the box), so that the box lengths agree. if you have non-periodic boundaries, just enlarge the box with an empty group.

you can combine the two boxes using read_restart add, but for that to work well, you need to do some extra steps when loading the first box.
if you have different atom/bond/etc. types, you need to use the extra/*** keywords to add a sufficient number of atom/bond/etc types and special atoms, since those numbers are locked in when the simulation box is created. to avoid complications from image flags, you should use the change_box command after loading the first system, so that there is sufficient space in z when adding the second system. otherwise LAMMPS will enlarge the box as needed, but that part of the code does not know about image flags and thus if those are non-zero, you get a mess. then you can use read_restart add with the needed offset parameters and shift for displacing the coordinates. i would then write out the combined data file and check it carefully and visualize it.

please note that this can be done in steps and that you need to use the clear command multiple times, if you do this from within a single script.
but best is to put each step into a separate file (unless you need to process a very large number of files following a pattern).

HTH,
Axel.

P.S: for some cases, the route via external tools may be simpler. I have written a tool in VMD’s topotools plugin that merges systems. in VMD changing the box does not affect coordinates (since they are unwrapped in the first place). but then you need to do displacements manually and also change the box accordingly and there are no checks at all. so it requires more care, but may be simpler for specific systems.

4

Thanks a lot for both of your inputs and the detailed answers, this gives me quite a lot to work with that I hadn’t thought about at first! Both of your answers complete each other pretty well too, which is greatly appreciated, thanks again for your input!

Luckily for me, both boxes are of the same atom type and will use the same potential and pair_coeffs, but its periodic boundary conditions will definitely cause issues I hadn’t foreseen, I’ll try all of this out!

Have a great day!