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