We first equilibrate the systems for 10 ns at a constant pressure P = 1 atm and a constant temperature T = 600 K. The Nóse-Hoover method is employed for both the barostat and the thermostat with a coupling constant τ P = 1 ps and τ T = 0.1 ps, respectively. Subsequently, we quench the system with a cooling rate 1 K/ps to low temperature. Further equilibration with 10 ns runs under the NPT ensemble is taken to ensure that the density of the system converges to a plateau value. The production run of the system at each temperature is carried out under the NVT ensemble for calculations of various observables. Simulation observables, such as radial distribution functions and pressure, are calculated from the production runs.