You have to set the temperature and pressures values at which you want to obtain the properties of liquid methane. Try to first relax the liquid by using fix nve + berendsen termostat.

What kind of system are you simulating? Interaction potentials, structure..?

Please attache the file which contains the structure and those which includes all LAMMPS commands.

My system is 3375 methane molecules. I used periodic boundary conditions, real unites and full atom style in my input file. Other commands in my input file are:

pair_style lj/cut/coul/cut 12.0

bond_style harmonic

angle_style harmonic

pair_modify mix geometric tail yes

special_bonds lj/coul 0.0 0.0 0.5

pair_coeff   1 1  0.066 3.5      # CT

pair_coeff   2 2  0.03 2.5      # HC

bond_coeff   1   340.0 1.09     # CT-HC

angle_coeff  1    33.0 107.8    # HC-CT-HC

velocity all create 93.7 53244 dist gaussian mom no rot no

fix 1 all nve

run 1000

minimize 1.0e-6 1.0e-8 10000 100000

timestep 0.25

reset_timestep 0

neighbor 1.5 bin

neigh_modify every 10 delay 20 check yes

thermo 400

thermo_style multi

unfix 1

fix 2 all npt temp 93.7 93.7 50 iso 1.45 1.45 10000

run 4000000

The mean pressure in 4000000 step is typically -2.5 to 9 atm (in separate simulations), but, in a single run, the pressure varies between -3200 to 1500 atm. Is it my fault or these types of fluctuations are in the nature of MD simulations?

In deed, the pressure is one of the most fluctuating properties in MD simulations. But your values seems to be very high. Please send me your system config file via personal message or attache it to an answer. I will have a detailed look at it.

I sent you my data file and 3 figures (pressure, density and temperature) during simulation  for 3375 methane molecule. I made a cubic box of CH4 molecules which the methane molecules were in the edges of a cube like NaCl crystallographic structure. This is my starting configuration of the liquid methane. Temperature was kept constant in 93.7K and the pressure at 1.45 atm. I hope that the methane could be liquid in this temperature and pressure.

I would be highly appreciated if you could tell me if I am wrong or not?

Dear Mohammad,

first a general device due to such questions in RG:

You can make the life of the tutors much easier if you give straight a way the full information about your problem/system. Which means in your case: attach _all_ necessary files. Your copy/paste action of the config file neglects the header, so I needed to work out the atom_style and units first before I could start a reproduction run.

About your issue:

As mentioned already, the pressure is in general a properties which fluctuates much stronger, than for example the temperature in general in a MD run. The fluctuations depend on the compressibility. Therefore the fluctuations are stronger for a liquid, then for a vapor phase.

You mentioned earlier that your pressure fluct. from  -3200 to +1500 atm, which is only the case for the very first steps. I would guess your preparation run (minimize) is to short.

The expected density value for 93K and 1.45atm is 0.44849 g/cm3, as taken from NIST. As the density vs time plot shows, the system is not equilibrated at the end, so you haven't reached the final density at all. I don't now which force field you are using (reference?) so I don't know if it is capable to model methane in this T/p region.

Best,

Frank

Dear Frank

Thank you for your advices. I sent you my data file in previous post. Here is my complete input file.

Best Wishes

Mohammad

If you check the density profile, your system is not arrived the equilibrium. Firstly, you need longer run.

What do you mean by "separate simulation", "single run"?

Dear Baris

I repeat my simulation several times and the pressure varies between -2.5 to 9 atm.

If your system is in equilibrium, the average pressure should be equal to that you set in your input file. In addition, the kinetic term contributes to the total pressure of your system. If you run your simulation at 1 K (via fix npt you can't set 0 K as indicated in the manual) you'll have less fluctuations.

Dear Baris

Tanke you or yor suggestion; but, I want to simulate the liquid methane. Is it valid to set the temperature at 1K? If I do this, can I get the correct properties of liquid methane?

You have to set the temperature and pressures values at which you want to obtain the properties of liquid methane. Try to first relax the liquid by using fix nve + berendsen termostat.

Dear Baris

I used your suggestions in my simulation (3375 methane molecules in cubic box at 93.7 K and 1.45 atm (1ns)) and the results are:

Mean Temperature In Simulation Is (K):                    9.368719e+001

Mean Pressure In Simulation Is (atm):                     -3.082390e+000

Mean Total Energy In Simulation Is (Kcal/mol):       1.999199e+003

Mean Volume In Simulation Is (A^3):                         1.815718e+005

Mean Density In Simulation Is (g/cc):                         4.952243e-001

The figures of Density, Pressure and Temperature are attached to this letter. Do you think that this simulation is correct? Can I use the output of this simulation for calculation of Radial Distribution Function and other thermodynamic properties of liquid methane?

fix NVE + langevin thermostat (or/and barostat) are not physical ensembles. After relaxing the structure, you should use a real ensemble simulation, fix NVT or fix NPT.

Dear Baris

My fix commands are:

fix 1 all nve

fix 2 all langevin 93.7 93.7 50.0 6243

fix 3 all npt temp 93.7 93.7 50.0 iso 1.54 1.54 10000

Is the previous posted results suitable for calculation of radial distribution function?

The above combination of fix's is not a real ensemble. It's just for relaxing the structure. After relaxing, you have to use 'real' ensembles, such as fix nvt, fix npt for the production of data from which you can calculate RDFs. 

Dear Baris

Thanks very much. I will try your suggestions.

Best Wishes

Hi, 

The pressure does fluctuate a lot but given a required period of time it will average out to the desired value, (iso 1.45 1.45 100000) I think your damping factor is quite high, try setting it to a lower value like around 100.