I am trying to simulate a large (I hope its large given how large systems these days can be : - ) ) system (~ 1 million particles).
The input script I am using is:
# -- initialise --
units real
atom_style full
boundary p p p
timestep 10
region simbox block 0 325 0 325 0 325
create_box 2 simbox bond/types 1 angle/types 1 extra/bond/per/atom 1 extra/angle/per/atom 2 extra/special/per/atom 5
info system
# -- force field --
include "ff.settings"
mass 1 14.0270
mass 2 18.015324
# -- generate box --
molecule pol polymer.mol
molecule mw mw.mol
create_atoms 0 random 1 123512 simbox mol pol 123135
create_atoms 1 random 1133740 123512 simbox mol mw 425641 overlap 2
mass 1 14.0270
mass 2 18.015324
set type 1 charge 0
set type 2 charge 0
write_data input_structure.data
quit
LAMMPS gets stuck when I run this script. It was not able to do anything for days on 96 cores.
Can someone tell me if I am doing something wrong or missing anything here?
For such large number of atoms to insert, it could be best not to use the random style, but instead place the atoms (or molecules?) on a regular grid.
Alternatively, with the random style, you can remove the “overlap 2” to speed up the insertion, and then perform an energy minimization to deal with the atoms that are too close to each other. This option is a bit dirty and would only work well for simple enough species, like simple particles, but will likely fail for more complex molecules.
You could also create a smaller isotropic sample and replicate it on a supercell. A bit of thermal annealing at fixed volume and (if it makes sense) with a simpler/faster pair potential can then be used to homogenize the supercell.
That is not very big by today’s standards. People have run LAMMPS with tens of billions of atoms.
Filling a volume with atoms in random places is intrinsically non-parallel, and in particular with overlap detection which requires additional communication, as @simongravelle noted.
More importantly, it is a hugely inefficient process to set up a system this way, as you need to run equilibration for a significant time. As @hothello noted preparing and simulating a smaller system until equilibration and then replicating it and equilibrating a little more with a dissipative thermostat to “heal” the intrinsic periodicity from the replication is the winning strategy. While starting from random geometries makes for a more unordered system, starting from grid positions is more problematic as it can take a very long time to get rid of the “memory” of the ordering. You have to ask yourself why are you running such a large system? and how much would each of those approaches mentioned above interfere with getting accurate results you desired.
All the advice so far in the thread has been really helpful and you should pay attention to it. I’m just giving followup advice based on having to do something similar recently – a lot of this is repeating what’s been contributed earlier, and remember that ultimately you are responsible for your own research.
For a million molecule situation, I would start with placing a thousand molecules using create_atoms in a much, much larger box than necessary, followed by replicate 10 10 10 to create that million molecule box, and then minimising and annealing.
Making an oversized box and replicating it before minimising helps the subsequent annealing somewhat because a configuration with very low mass density intrinsically has a bit more entropy, and we need every little bit of entropy we can get to make sure we’re not simulating an artificially over-ordered set of configurations.
However, bringing that configuration into a condensed-phase equilibrium gives off a lot of heat (the enthalpy of condensation!) which can “blow out” the subsequent simulation, even if there are no initial steric clashes. Energy minimisation is mandatory. In my particular case, I then “froze” my sample (npt integration at 30 K and 1 atm; running such low temperatures gives the thermostat and barostat more leeway to smoothly remove the enormous heat of condensation) for a few nanoseconds before bringing it gradually back to room temperature and pressure.
The simulation still needs much, much more annealing after that; I used the protocol from this paper: https://pubs.acs.org/doi/10.1021/acsapm.1c01798 and it seems to have given good results so far.