# Very high negative pressure in MD simulation of liquid quenching

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

# ------------------------------------------
# -------- 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"
#############################################


• off hand, 60 picoseconds is not long, but I will admit, it is usually long enough to equilibrate mechanical properties like pressure. What is your desired temperature you want to run at for the production stage? Commented Jan 22, 2022 at 15:00
• The desired temperature ranges between 1K and 300K for the production stage.
– Abd
Commented Jan 22, 2022 at 23:56
• During cooling down the annealed sample to 1K or 300K, the pressure goes around -100,000 bars. To readjust this pressure to 1 bar changes atoms repositions in a dramatic (mean-square-deviation up to 25 Angstrom in 5 ps) manner along with density change from 1.8 to 1.96 gcm**-3. I would like to avoid such a large variation, that appears only due to barostat, which is the last step of the sample preparation. Or at least I need to know that this is normal and the final sample can be used for investigation. Thanks!
– Abd
Commented Jan 23, 2022 at 0:33
• Try using the Parrinelo-Rahman barostat instead and see how that works Commented Jan 23, 2022 at 8:36
• I use lammps for simulation. The npt ensemble in lammps uses Parrinelo-Rahman strain energy description, and this is the one I've used for basrostating. Not sure if you are talking about a different barostat. There's another ensemble in lammps names "npt/cauchy" uses the Parrinelo-Rahman strain energy and changes simulation box size. Haven't used this one yet. After your comment, I think this is what I'm going to try. Thank you.
– Abd
Commented Jan 23, 2022 at 10:50

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.

• Thanks for taking the time! this is really helpful and I agree with you. I'm moving forward with the results I got and applying barostat before production.
– Abd
Commented Jan 31, 2022 at 23:42