Density Profile for Tersoff Parametrization

I am having some difficulty acquiring an accurate density profile for the Tersoff parametrization for quartz. I first minimized the structure, then ran an NVT ensemble to equilibrate the system at 300K. The density after minimization was around 2.42 g/cc. Once I run the NPT ensemble, I see an increase in density. Now the minimization takes place at 0K. The NVT takes place at 300K but dimensions of the box cannot change. Running NPT (where the density can change), I would expect density to decrease due to the temperature increase; however, it increases to about 2.46 g/cc. Is what I am seeing a shortcoming of these parameters or could it be an error in my script?

I feel like what I did was fairly straightforward and am concerned with how the density profile looks. I attached a graph. At timestep 20,000 (x-axis), you can see the density increase. This is where the switch from NVT to NPT takes place.

#Equilibration for quartz using the Tersoff Potential
units metal
atom_style full

#Data format is read in below
read_data …/…/quartz.data

####################Potential#############################

pair_style tersoff
pair_style tersoff
pair_coeff * * …/…/…/…/potentials/SiO.tersoff O Si

Tersoff_density.PNG

I am having some difficulty acquiring an accurate density profile for the
Tersoff parametrization for quartz. I first minimized the structure, then
ran an NVT ensemble to equilibrate the system at 300K. The density after
minimization was around 2.42 g/cc. Once I run the NPT ensemble, I see an
increase in density. Now the minimization takes place at 0K. The NVT takes
place at 300K but dimensions of the box cannot change. Running NPT (where
the density can change), I would expect density to decrease due to the
temperature increase; however, it increases to about 2.46 g/cc. Is what I am
seeing a shortcoming of these parameters or could it be an error in my
script?

what makes you think this is bad?

the information that i am missing in your description is the density
that was reported in the reference where your potential was taken
from.

axel.

The reference said that the density was 2.42 g/cc at 0K and 0atm. This is indeed what I found. There was no reference to what happens at a higher temperature. Increasing temperature should decrease the density. This does seem to be the case when I heat it up from 300K to 5000K. I do not understand this initial jump in density at a higher temperature and it seems unphysical to me.

Ben

The reference said that the density was 2.42 g/cc at 0K and 0atm. This is
indeed what I found. There was no reference to what happens at a higher
temperature. Increasing temperature should decrease the density. This does
seem to be the case when I heat it up from 300K to 5000K. I do not
understand this initial jump in density at a higher temperature and it seems
unphysical to me.

well, you should carefully compare the structures and also
cool/minimize the 300K system down to 0K to see, if you get back the
original data.

axel.

I tried running it back down to 0.0001 K (since I cannot use 0K with NPT). The density actually increases up to about 2.5 g/cc when temperature is cooled from 300K to almost 0K. This does seem physical to me. That density is increases as temperature decreases. I do not understand the why the initial structure minimizes to 2.42 g/cc.

Metastable geometry due to symmetry?

Axel, I am not sure what that means. Is there a reference you can point me to?

Ben

This is crystallography 101 stuff. Do as I suggested and carefully compare the geometries.

I can see right off the bat that the geometries are different - But I am not entirely sure what I should be looking for. Are you suggesting that my initial configuration was in a meta-stable state (or at a local minimum) only due to its symmetry and that heating and cooling provides a state that is no longer meta-stable?

Ben

After some work, I think I have found some solutions but request verification. I was minimizing the structure to 0 bars and 0K to start off that returned a density of 2.42 g/cc. I heated it up using NVT ensemble to 300K. From there, I noticed that when I switched to NPT ensemble at 300K and 0 bars, the density was increasing, which did not make any sense at first. So I knew something unphysical was taking place. Then I realized that the pressure was the reason why density was increasing. I printed out values of the pressure from the NVT ensemble and noticed that it converged to -3000 bars. If I input this pressure into the NPT ensemble instead of 0 bars, the density did decrease as it should. My concern is the input parameter of -3000 bars. This is on the order of -300 MPa which seems astronomical.

Another thing noted is the fact that when I heat up using NVT and then cool down using NPT, the density is not the same as in the minimized structure. It is actually about 2.5 g/cc, slightly higher than the original 2.42 g/cc. Axel suggested that this was due to metastable geometry due to symmetry.

I am not sure whether proceeding with an input pressure or -3000 bars makes sense physically. Is this quantity meaningful? Or would it make more sense to get out of the meta-stable geometry by heating and cooling and proceeding in that manner?

Best Regards,

Ben