Problems reducing the size of the simulation box

I’m trying to compress the size of a simulation box in lammps, I’m using Reaxff as the force field because I have a large system of clusters of aromatic layers (C,H,O,N)

To compress the box, I’m using this command after an NPT run at 298K, 1 atm


variable a loop 10

label loop

change_box all x scale 0.9 y scale 0.9 z scale 0.9

run 1000

next a

jump in.data loop

The simulation box initially has the following dimensions

-38.9259 561.074008 xlo xhi

-374.371 225.628502 ylo yhi

8.173006 858.173006 zlo zhi

I need to compress that to fit the density of a solid, so the box would need to be at least (60x60x60)

The output file after the NPT remained in this step for a long time, and nothing else is happening


variable a loop 10

label loop

change_box all x scale 0.9 y scale 0.9 z scale 0.9

Changing box …

orthogonal box = (-2.41712 -367.139 18.4184) to (524.565 218.396 847.928)

orthogonal box = (-2.41712 -337.863 18.4184) to (524.565 189.12 847.928)

orthogonal box = (-2.41712 -337.863 59.8939) to (524.565 189.12 806.452)

run 1000

Setting up Verlet run …

Unit style : real

Current step : 500100

Time step : 1

“biochar.out” 287L, 140

Any thoughts about what the problem would be??

There is not much useful information here to base any suggestions on.

Some (random?) thoughts:

  • The change of the box volume by a factor of nearly 1500.0 seems pretty drastic. What phase is the system in originally?
  • The segment of code you present only induces a change of the volume by a factor of less than 2.5
  • You seem to be running with fix npt between each volume rescaling, which makes little sense. What should keep the system from expanding back again to the equilibrated volume?
  • The quoted output suggests you are using a timestep of 1fs. That seems extremely large for a system handled by ReaxFF. I would expect something in the neighborhood of 0.25fs or even less
  • ReaxFF is very sensitive toward large changes in the environment of atoms. It would be much easier to do the compression with a classical force field and then switch to ReaxFF after that system has equilibrated (and the explicit bonded interactions removed).
  • Have you tried/tested this procedure with a simpler input? and/or simpler force field?