Need help for NPT problem

LAMMPS LAMMPS Beginners
1

Dear experts,
I am student studying glass molecular dynamic simulation using Lammps-29Aug2024.

I am currently facing issues with the initial energy stabilization and a drastic increase in the simulation box size after the NPT process, which results in an excessively low density. I initially set the box size based on the experimentally obtained density value.

I am using the following input file for my simulation:

Ewald Summation for Coulomb

kspace_style ewald 1.0e-4

minimize 1.0e-6 1.0e-8 50000 100000

Simulation settings

timestep 1.0
thermo 100
dump 1 all xyz 1000 dump.xyz

velocity all create 300.0 12345 mom no rot no
velocity all scale 300.0

fix 1 all nvt temp 300.0 300.0 100.0
run 40000
unfix 1

Initial Equilibration (6000K, NVT)

fix equil1 all nvt temp 300.0 6000.0 1000.0
run 100000
unfix equil1

Melt Equilibration (5000K, NVT)

fix melt all nvt temp 5000.0 5000.0 100.0
run 40000
unfix melt

Quenching (5000K → 300K, NVT)

fix cool all nvt temp 5000.0 300.0 2500.0
run 500000
unfix cool

Final Equilibration (300K, NPT)

fix final_eq all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run 30000
unfix final_eq

fix nvt_eq all nvt temp 300.0 300.0 100.0
run 40000
unfix nvt_eq

When running this simulation, I observe that as the NVT stage ends and the NPT stage begins, the pressure value is initially high. As the system attempts to reach 1.0 atm, the simulation box expands significantly, leading to a very low density.

To resolve this issue, I am wondering if I need to further stabilize the initial energy. Additionally, if there are any problems with my input file, I would appreciate your guidance.

2

Are you absolutely positive that you are inputting the correct potential parameters? This problem could come from a mistake there. In this case your atoms would be experiencing an unrealistically high repulsive force that justifies the high pressure and thus the increase in volume to get to 1.0 atm.

3

Why do you use “ewald” and not “pppm”?
Also, a convergence of 1.0e-4 results in (mostly) converged forces, but not converged stresses.
Try 1.0e-6.

Why the different thermostat time constants?
Why not using a dissipative thermostat, e.g. fix nve + fix langevin?

Why such a short simulation time for this part?
What is the average temperature of the system when you start fix npt?

You need to get this from a local tutor or adviser that is training you in doing meaningful MD simulations and can look over your head and discuss general details about best practices in performance MD simulations. An online forum does not work well as a teacher.

4

Thank you for your guidance. I will double check on the potential parameters.

5

@akohlmey First of all, thank you for your detailed guidance. It is very helpful since I am studying by my self. I will study and fix each sentence one by one.
I am practicing glass simulation by modifying https://doi.org/10.1111/jace.16082 this paper.
Since I was trying to verifiy the success of building the structure model caused the short simulation time.

6

That is a very bad idea. To learn MD properly, you need proper in-person guidance (at least for the initial time) from an experienced practitioner. Otherwise, you are going to make lots of trivial mistakes and not all of them will manifest in LAMMPS stopping or crashing. An online forum cannot replace that kind of tutoring and there is no proper textbook that contains all the little practical pieces of advice that you can only get from a practitioner. Doing MD simulations is as much a craft as it is a science.

Please don’t just take my word for it, but read from a person that has been there: How to calculate the amount of absorbed gas - #5 by srtee

7

Thank you for your wise advice. My co-workers at my institution have experience in MD simulations but lack knowledge of LAMMPS programming. I will find a proper teacher to learn from in person. Thank you so much!