Hello,
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?
Thanks