How to Acheive density independent

Dear lammps users,

I am simulating methane system.
I want to achieve slush state of methane. That is 50% solid and 50% liquid state.
The continuum state pressure = 11.7 kpas, Temp = 90.7 K, and slush density observed is 481 kg/m3.
Now, when i consider a simulation box, 50 Angstrom length each side, with 752 methane molecule particle, and use npt fix with same pressure and temperature for continuum study. But the density observed is 0.25 kg/m3.
How can i ensure my system achieve slush state with this density? Should I change the number of particle or simulation box size?

The melting point would be a (T,P) combination in which the two phases are equally stable, so in fact there is no driving force for the transformation (i.e., phase transition) to occur (and also you should be able to have all possible % partitions into solid/liquid without getting out of the given (T,P) condition). I never simulated coexistance of phases, but my best guess would be to first do a simulation of the solid (at a thermodynamic state at which this is the stable phase) and then evolve it up to the melting point of choice (hoping that no phase transition would occur). Then you do the same thing with the liquid in another simulation. If you save a final configuration for these two independent systems, you could then build something that resembles a configuration for the 50/50 system (obviously you would need to use the same amount of each phase to build this new configuration). (PS: note that the density of the liquid and solid will be different). Since there is no driving force, your system should somewhat maintain this composition apart from (hopefully) soft changes during an equilibration.

Also, you need to see if your force field is actually able to reproduce the phase diagram (or at least this specific melting point), otherwise it will not work out. The force field would also need to be able to predict the density of the two phases at your target thermodynamic state if you are specifically worried about reproducing the density of the overall system.

I dont see how you would get to a system of such specific composition by starting from one phase only.

If you want to have a specific density, you don’t use fix npt, but fix nvt. Then your density is an input parameter.

If you want to simulate coexistence at the molecular scale you have to follow the advice of @ceciliaalvares and set up a two phase system directly by from your input. Phase transitions have a hysteresis and at the atomic scale that one is huge. You can have a solid way beyond the melting point and vice versa since both melting and freezing require nucleation.

However, there is a simpler to achieve this. Define two groups, solid and liquid, where you designate the atoms that you want in either phase starting from a bulk solid. First you equilibrate the whole system to the melting point without melting. Now you do a run with fix nvt only applied to the liquid group at a temperature way above the melting point until your system is melted and then slowly cool it down to the melting point. Then you re-quilibrate the solid only by using fix nvt only applied to the solid group. You may need to alternate between the two one or two more times until either side is properly equilibrated to the melting point and then you can go on to simulate the whole system with nvt.

Please note that the relevant temperature and density is determined by the potential you are using, not the experimental values. Thus it is advisable to first set up your system so that you have two halves. When running the whole system you can determine the melting point as the temperature where neither the liquid nor the solid phase is growing or shrinking. Otherwise you need to adjust the temperature accordingly until you have found the suitable one for your model.

There is lots of published literature on this and also plenty of discussion here in the forum archives.

Thank you so much.

Thanks a lot

I still facing the problem to achieve pressure, density and temperature at desired value. I created a system with particular density, say rho-1. I run with NPT first, so pressure and temperature I attained. But after equilibrium reached, i shifted to NVT. At the same time, my system density rho-1 changed to the NPT equilibrium density, rho-2. If i use only NVT, i wont get required pressure.
Can you please suggest me a solution?

You didn’t pay sufficient attention to what I told you.

If you don’t get the desired pressure at the given density, then the cause is with the force field. If you run with fix npt, you VERY OBVIOUSLY will not get the desired density since fix npt will change it to obtain the desired pressure. So, if the pressure is not matched with fix nvt using fix npt will not cure it. That is just logical thinking!

How far off is the pressure? Are you taking the time average of it?

The pressure is far different. But let me check with the force field once again. I consider OPLS force field.

That is a very unscientific statement. Why can’t you just specify some numbers?
You also didn’t answer my second question.

Liquids and solids are not very compressible, so what may seem large to you, may in real life not be a big deal considering that pressure on atomic scale simulations can fluctuate a lot.

Please also note that none of this is a LAMMPS problem, but rather an issue of understanding simulations.