Very high negative pressure in MD simulation of liquid quenching

Question

I am trying to prepare turbostratic carbon structure from an initial random carbon sample with a density of 1.8gcm$^{-3}$. The initial sample image is attached here:

The problem is the high negative pressure (up to -100,000 bars) during the quenching process. You can see the pressure distribution as the blue curve of the following image:

I understand pressure can be negative in MD simulations. However, this is a very high negative pressure in my opinion, and something is going wrong here I'm guessing. Moreover, it creates problems later on when a barostat is applied to equilibrate sample at room pressure (e.g. stress distribution in atoms changes completely after barostating, due to large adjustment of atom positions corresponding to large pressure adjustment). I'd really appreciate any advice on getting rid of such high negative pressure. I have used "Bussi1" thermostat along with NVE integrator to control temperature and update atom positions and velocity respectively. below is my script

# ------------------------------------------
# ----- Initialization----------------------
# ------------------------------------------

units metal
atom_style atomic
boundary p p p
newton on
dimension 3

# reading initial/restart sample -----------
read_data MDSample.dat


# ------------------------------------------
# -------- Simulation setting --------------
# ------------------------------------------

# pair potential----------------------------
pair_style tersoff
pair_coeff * * SiCGe.tersoff C

# time step---------------------------------
#reset_timestep 0
timestep 0.00005

# minimization -----------------------------
min_style cg
minimize 1.0e-20 1.0e-20 1000 10000

# set intial temperature -------------------
velocity all create 100 1201

# compute virial stress of each atoms ------
compute atomstress all stress/atom NULL virial

# compute ke of each atoms -----------------
compute atomke all ke/atom

# compute pe of each atoms -----------------
compute atompe all pe/atom

# compute global msd ----------------------- 
compute mymsd all msd

# macroscopic/ensemble average output variables --------------------------------
thermo_style custom step dt time atoms temp press pe ke etotal evdwl vol density c_mymsd[*]
thermo 1000


# -----------------------------------------
# ---- Melting simulation(8000K) (5ps)-----
# -----------------------------------------

# set ensemble -----------------------------
fix 1 all nve
fix 2 all temp/csvr 8000.0 8000.0 0.1 54324

# print restart file -----------------------
restart 100000 MeltSample.restart

# print sample file for visualization -----
dump sample all custom 100000 MeltSample.txt id x y z c_atomstress[*] c_atomke c_atompe

run 100000 #--------------------------------
unfix 1
unfix 2


# -----------------------------------------
# ---- Quenching simulation(300K) (5ps)----
# -----------------------------------------

# set ensemble -----------------------------
fix 3 all nve
fix 4 all temp/csvr 300.0 300.0 0.1 54324

# print restart file -----------------------
restart 200000 Quench300K.restart

# print sample file for visualization -----
dump sample1 all custom 200000 Quench300K.txt id x y z c_atomstress[*] c_atomke c_atompe

run 100000 #--------------------------------
unfix 3
unfix 4


# -----------------------------------------
# ---Annealing simulation (5000K)(50ps)---
# -----------------------------------------

# set ensemble -----------------------------
fix 5 all nve
fix 6 all temp/csvr 5000.0 5000.0 0.1 54324

# print restart file -----------------------
restart 1200000 Ann5000K.restart

# print sample file for visualization -----
dump sample2 all custom 1200000 Ann5000K.txt id x y z c_atomstress[*] c_atomke c_atompe

run 1000000 #--------------------------------
unfix 5
unfix 6


#############################################
print "All done"
#############################################

Answer 1

That kind of pressure feels reasonable.

I don't think there's any cause for concern based on the data you have shown. Consider the following back-of-the-envelope calculation:

• The thermal expansion coefficient of carbon is (according to WebElements) roughly 7.1 ppm/K. That's the linear (length) expansion coefficient, the volumetric expansion coefficient is approximately three times larger. Let's take $\alpha_V$ = 20 ppm/K as a rough guess for your material.
• Assuming that the thermal expansion coefficient is constant across the entire temperature range (which is a really crude assumption), cooling your material down by 7700 K will make its volume shrink by 15.4%.
• The bulk modulus of carbon is (again according to WebElements) K = 33 GPa. Since $K = - V \frac{dP}{dV}$, let's assume it's again constant during the whole quenching procedure and thus $\Delta P = - K \frac{\Delta V}{V}$. This gives us $\Delta P = -33 \cdot 0.154 = -5\,\mathrm{GPa}$, or -50000 bar.

Given the (in)accuracy of this estimate, it is fairly OK that your observed pressure change is just a factor of 2 off. You could make a better estimate using a real equation of state for carbon, but note also that the empirical force field will on its own likely have a limited accuracy in reproducing the true EOS.