Adding thermalized supercell to running simulation

Hi All,

I have read in the literature about a ‘moving window’ molecular dynamics technique which involves observing some dynamic phenomenon over long time scales by feeding in fresh material upstream and removing it downstream (and observing said phenomena in some steady-state sample region in between). This has been done for metals previously, and I am interested in repeating this for more complex materials that have a longer equilibration time. As a first set of ‘proof of concept’ experiments, I am attempting to join two thermalized supercells of a crystalline material during NPT relaxation using the following scheme.

  1. Add 4x4x20 supercell, assign random velocities, NPT equilibrate for 2 ps
  2. Use read_data add append shift to add another 4x4x20 supercell to the first along the z-axis, assign random velocities to these atoms corresponding to EQ temp
  3. Continue NPT equilibration for 8 ps

Code is attached at the bottom. From what I have read in the documentation, it seems like this should work (given the description of the read_data append add commands), but I am running into a segmentation fault (11) on my local cluster shortly after LAMMPS attempts to insert the new atoms. I have varied the shift distance to ensure that this problem isn’t due to atoms overlapping and continue to run into the same problem.

I am now testing a version where the NPT is unfixed before adding the atoms (since N is no longer constant), and then refixed, but I am not optimistic that this will work. The sense I get from the documentation is that the primary use case for adding data files to a system is to prepare a complex initial state, not modify one that has already been running.

Has anybody done something similar, and if so, how did you work around the SEGFAULT?


log test_relax.log

units           real
atom_style      full
boundary        p p p
timestep 1

bond_style      harmonic
angle_style     harmonic
dihedral_style  hybrid harmonic multi/harmonic
improper_style  harmonic
#pair_style      buck/coul/long 11.0 11.0
pair_style              hybrid/overlay table linear 2550 coul/long 15.0
kspace_style    pppm 1.0e-6
#kspace_modify   slab 3.0
special_bonds   lj/coul 0.0 0.0 1.0


pair_coeff 1 1 table HMX_X6table.txt HH
pair_coeff 1 2 table HMX_X6table.txt HN1
pair_coeff 1 3 table HMX_X6table.txt HN2
pair_coeff 1 4 table HMX_X6table.txt HO
pair_coeff 1 5 table HMX_X6table.txt HC
pair_coeff 2 2 table HMX_X6table.txt N1N1
pair_coeff 2 3 table HMX_X6table.txt N1N2
pair_coeff 2 4 table HMX_X6table.txt N1O
pair_coeff 2 5 table HMX_X6table.txt N1C
pair_coeff 3 3 table HMX_X6table.txt N2N2
pair_coeff 3 4 table HMX_X6table.txt N2O
pair_coeff 3 5 table HMX_X6table.txt N2C
pair_coeff 4 4 table HMX_X6table.txt OO
pair_coeff 4 5 table HMX_X6table.txt OC
pair_coeff 5 5 table HMX_X6table.txt CC

pair_coeff * * coul/long

dihedral_coeff          1 multi/harmonic 8.452 0.000 -5.28998 0.00 -3.16003
dihedral_coeff          2 harmonic -0.08 -1 3
dihedral_coeff          3 multi/harmonic 0.095 -1.485 1.61000 -0.22 0.0

neighbor          2.0 bin
#neigh_modify  every 1 delay 5
neigh_modify    delay 0 every 1 check yes one 10000 page 100000

thermo_modify  line one flush yes

velocity all create 300.0 2151658 dist gaussian mom yes rot yes

compute		PPA all stress/atom NULL
compute 	vorvol all voronoi/atom peratom yes
compute		vz all property/atom vz

variable 	meanP atom -((c_PPA[1]+c_PPA[2]+c_PPA[3])/(3*c_vorvol[1]))*0.000101325 #const added for GPa conversion

thermo 10
thermo_style custom step vol temp press pe ke etotal epair evdwl ecoul elong emol ebond eangle edihed eimp enthalpy density lx ly lz pxx pyy pzz pxy pxz pyz
thermo_modify line one

shell mkdir addDumps

dump                    d1 all custom 500 ./addDumps/Dump.*.dump id mol type xu yu zu vx vy vz fx fy fz
dump_modify             d1 sort id

fix 2 all npt temp 300.0 300.0 100.0 aniso 1 1 1000

run 2000

variable zh equal zhi-20+zlo+3

unfix 2

read_data add append shift 0 0 ${zh}

velocity all create 300.0 2151654 dist gaussian mom yes rot yes

fix 3 all npt temp 300.0 300.0 100.0 aniso 1 1 1000

compute dz all chunk/atom bin/1d z 10 6 units box

fix pressure_profile all ave/chunk 10 5 100 dz v_meanP

run 8000

There are two issues here:

  • the fact that you get a segmentation fault
  • the model that you want to implement

To track down the segmentation fault, we need to have more information and - ideally - enough information (and the inputs) to reproduce it. That is assuming that the segfault is a “real” error and not just a consequence of LAMMPS aborting somewhere and the MPI library not stopping properly. For that we need to know what version of LAMMPS you are running with and on what platform with what settings (plus the inputs) and a representative log file or screen output.

Some comments on the second item. What you are after is a complex thing and I don’t see your input coming even close to achieving what you describe. For example, it cannot work with fully periodic boundary conditions. You need to have non-periodic boundaries in the direction where you “feed” the simulation new material and remove “processed” material.
You would have something that induces a flow of material (e.g. by using fix move on the inserted atoms, then switching to thermostatting them while having enough space between that zone and the zone where the point of interest is and then a sufficient buffer of “processed” material before it will be discarded). You also need to make certain that that addition of atoms does not cause a shock (atoms that previously had no neighbors suddenly have some and that is the same as if they were hit very hard with a blunt object).
Thus, as an alternative, I would urge you to consider a periodic flow simulation setup where you instead have two regions, one where the point of interest happens and a second where you “repair/restore” the system to its original state. That is much easier to achieve in a simulation because it does not require adding/removing material and thus avoids a lot of complications.

There is one case where LAMMPS explicitly supports the situation you describe and that would be a traveling shockwave. For that the SHOCK package has the fix append/atoms command — LAMMPS documentation. If you need to set up something similar, I would suggest you try to make that kind of process work first, as it should be much easier to append new atoms on fixed lattice positions than adding atoms from a data file.

1 Like

Thank you for the detailed response and helpful feedback! What I am going for is complex indeed, and what I have presented shouldn’t even be characterized as a facsimile of the intended end product, rather it is a fraction of a baby step in the right direction to characterize the process of removing the free surface in the upstream direction and the pressure wave that results (i.e. to see if thicker additions or thinner additions have a weaker influence on the sample region of interest). I have stuck with PBC’s here to continue my use of the PPPM k space style, though I suppose that I could switch to a fixed BC and use the slab option in k_space_modify.

That said, your pointer to the append atoms command is very helpful and think is very suitable for my use case. My only concern at this point is the fact that I am working with a molecular solid where a molecule is many atomic layers thick. This is why I was trying the ‘append a slab of unit cells in the direction of shock propagation’ approach. If I start adding on an atom-by-atom layer basis I think this would cause problems with huge energies near the surface when only parts of a molecule have been added, exacerbating the issues you mentioned with changes in the thermodynamic properties in the region of interest. I may be misinterpreting the documentation, but it appears that additional material can only be added atom by atom. Or, if I specify the correct z parameter in append/atom with the right basis, can I get a layer of unit cells appended?

Your suggestion of a periodic flow simulation is very intriguing, and I think is probably the right way to go about my set of experiments. However, I don’t think I will be able to avoid the issues related to adding/removing material, as there needs to be some initial impulse from a piston to get the shock going initially and keeping the piston material would result in a nightmare of reflected shocks in the sample region (if the piston is not fixed, which is the only way I could imaging going about the flow simulation you describe). I don’t have must experience with MSST, but it seems like this could work for my use case, but I would like to keep my simulations entirely atomistic without bringing the continuum approximations into it.

You may find some success locating prior literature if you search for “open boundary conditions” – as in, the boundaries of your simulation are open to a reservoir, instead of being either closed or periodic. For example, this recent paper looks very interesting and may have solved or addressed questions you’re asking: Cookie Absent

1 Like