Trouble preserving density

Hello LAMMPS users,

I intend to calculate some of the thermo-mechanical properties of bulk alumina. First, the box needs to be relaxed. For relaxation, I usually employ a “fix npt” at a specific pressure and temperature. However, I came to realize that the box undergoes a decrease in density. The density starts at about 3.90, which is close to the density of bulk alumina, but after 50,000 steps I get the value of about 2.90, which is way below the desired density. I tried a “fix nvt” right before the “fix npt”, however, this also couldn’t maintain the density.

Could you please help me with that? What is the best to relax a box at a specific temperature and pressure?

Thank you all for your consideration.

It is difficult to make an assessment without knowing more details.

It is not very likely that the choice of time integration algorithm has a significant impact.
It is much more likely that your force field and/or starting geometry is unsuitable.

Thank you for your response.
Actually I am using the ReaxFF reactive force field to evaluate the thermal conductivity of bulk alumina.
For this, the box has dimensions of 30x30x500A.
Could it be due to timestep?

I already gave you my assessment based on the available information. It has not changed since you have not provided significant additional information.

Oh, sorry.
The 3D simulation box has dimensions of 30x30x500A.
I am trying to simulate bulk alumina at the temperature of 0K and the pressure of 1atm using the ReaxFF potential, and I use the timestep of 0.5fs. Using the data file, this is the code I use:

#---------- Simulation

thermo 1000
timestep 0.5

pair_style reax/c NULL safezone 2.9 mincap 5000
pair_coeff * * ffield.reax.AlO Al O
fix reaxfix all qeq/reax 1 0.0 10.0 1e-6 param.qeq

velocity all create 300 456783 mom yes rot yes dist gaussian

thermo_style custom step density temp etotal press pxx pyy pzz pxy pxz pyz vol lx ly lz

dump 1 all xyz 5000 dump.xyz
dump_modify 1 element Al O

fix NPT_fix all npt temp 300 300 5 aniso 1 1 12.5 drag 2
restart 10000 restart.*
run 50000

This is incomplete and thus does not allow to recheck/repeat your runs and neither does it show energies and pressures from your output. Also, how can I know that your force field parameters are parameterized for your type of system?

Here, I uploaded the input file along with the potential file I use for my simulations.
I really appreciate it if you can help with this issue.
input.in (3.2 KB)
ffield.reax.AlO (28.0 KB)

How do you know that this potential file is suitable? From the description it could just as well be meant for use with hydrocarbons and/or carbohydrates or similar interacting with metallic aluminium.

Thank you so much for your consideration.
I think, the main paper, proposing it, employed it to model oxide-coated aluminum nanoparticles. Also, our group utilized this function for many years and derived good results. Did you find anything else in the code or in the log file?

If people in your group have done simulations with this potential and on similar systems, you should ask them for assistance, not here. They most certainly have more experience than anybody here.

Sure. Thank you again for your time.

What stands out is that you are starting with a density of 3.9 but also have very high pressure.
So you should observe the geometry during those 1000 steps. Also you should do the same 1000 steps with fix nvt and compare. If the force field is suitable, then the geometry may be not correct. And that would become visible by a significant change in geometry. At 300K there should be none of that.

1 Like

Actually, you’re completely right. But what is the usual way you use to relax a box at a temperature and pressure? Do you use a single npt ensemble? Or maybe a series of different ensembles? I’m just curious to know what is the best way based on your own experience?

Your are too fixated on the command(s) to do the box relaxation.

The cause of your problems, however, must be a different one. If your initial geometry is bad, then no choice of integrators or simulation settings can salvage this. So you first need to confirm that your geometry is correct. If that is not the case, then it must be due to the force field parameters (despite your claims to the contrary) or settings related to that. You are looking at an extremely hard compound, so its compressibility will be small and thus energies and forces and pressure rise quickly if the input geometry is not perfect (and perfect means perfect for the force field parameters, not what you get from experiment). If I was in your place, I would research the literature for alternate force fields for the same compound and do experiments with those. Will they reproduce the expected density closer than what you are currently using? That would be a strong hint toward problems with the parameters you currently use. Or will they show the same behavior as the current one? That would be a strong hint toward issues with the geometry/box etc.

I am so very grateful for your help.

There are not many functions that can model the Al/O interactions. Other than ReaxFF, there is Streitz-Mintmire’s EAM+CTI potential function which can be used. However, it doesn’t support OPENMP, and KSPACE package must be employed, which is another trouble as I’m performing the simulations on Windows.

Thank you again for your time.

How is using windows a problem?
You can always use MPI parallelism and it is usually more efficient than OpenMP.

Can I use this KSPACE package on windows?
Or, are you saying that MPI supports “pair_style coul/streitz command”?
On the website I read that " The coul/cut/global , coul/long , coul/msm , coul/streitz , and tip4p/long styles are part of the KSPACE package. They are only enabled if LAMMPS was built with that package. See the Build package doc page for more info."

Since I’m new to LAMMPS, I don’t really know how to access this KSPACE package. Can you please give me instructions on how to use it?

That depends on how you obtained your executable. The info config command will tell you which packages are included.

If your LAMMPS executable has been compiled with MPI enabled, then yes. All of LAMMPS supports MPI. And if certain combinations of options cannot be used you with more than one MPI task, then you will get a suitable error.

1 Like

Since you said you are doing research in a group that has experience with LAMMPS you should ask to get local tutoring from your more experienced group members. Trying to teach you over the net is extremely inefficient (and don’t have the time to do it. You are stretching my patience to the limit already).