How to equilibrate amorphous silica model


I’ve been trying to equilibrate a simulation box containing amorphous silica (12000 atoms) using reax/c force field. This system was simulated using a protocol reported by Fogarty, et al. (JCP, 132, 174704 (2010)), in which a box of SiO2 molecules was heated to 4000 K and then quenched to 300 K at a rate of 100 K/4 ps. An initial heat-quench simulation was performed in the NVT ensemble (‘fix nvt’; TDAMP = 25.0 fs) and then a second heat-quench simulation the the NPT ensemble (‘fix npt’ TDAMP = 25.0 fs; PDAMP 250.0 fs) was also performed.

It’s been a challenge to equilibrate this system. The amorphous silica structure resulting from the two heat-quench simulations was equilibrated for ~2ns in the NPT ensemble (target pressure of 0.0 atm; TDAMP = 25.0 fs; PDAMP 250.0 fs). Over the last 20 ps of this equilibration stage, the average pressure (taken every 10th timestep) was 0.2030 atm (std. dev 796.6 atm). As a quick check of whether or not this system was nearing equilibration, I performed a 20 ps-long NVE simulation on the structure produced by the NPT equilibration. The average temperature of the system during the NVE simulation was 298.5 K (std. dev 1.5 K) and the average pressure was 229.4 atm (std. dev 747.4 atm)- which indicates to me that this system is still far from equilibrated.

I understand that for small amorphous solids, equilibration can be a tricky and time-consuming and 2 ns isn’t necessarily all that long of an equilibration. However, I am concerned that this approach will effectively never yield a equilibrated system due to a poorly chosen barostat and/or PDAMP parameter. FWIW, the density of this system after 2ns of equilibration is comparable to reported values for amorphous silica simulated using ReaxFF- as are the coordination number, bond lengths and bond angle distributions. So, in these respects, this simulated system appears to be reasonable. Does anyone have any advice on what I should do differently to equilibrate this system?


Aidan or Ray can likely comment since you are using ReaxFF.
200 atm with a std of 750, for a target of 0 does not sound that
bad. I am guessing the volume change for 200 to 0 is quite
small for a solid.


I am baffled by your question. Your result:

229.4 atm (std. dev 747.4 atm)

is statistically indistinguishable from the target pressure of 0 atm. Why
do you say this this not equilibrated? Were you expecting 0 atm +/- 0?
That would only be possible at T=0. I suppose a valid question would be:
is <delP^2> ~ (700 atm)^2 a reasonable level of fluctuation for this
system. The answer depends on things like the volume, the temperature, and
the bulk modulus. It is not a LAMMPS question.