Simulating Water Vapor

LAMMPS LAMMPS Beginners
1

Hello, I’m a beginner and still learning about LAMMPS and currently studying about water. Usually, I just read old mailing question from LAMMPS website and not confident to ask in the forum, but currently, it’s 404 not found. I wanted to know, is it possible to simulate water vaporization with LAMMPS? are there any example of it?
I can’t find material or paper about it, mostly I found about condensation or vapor-liquid system. I have tried to “modify” SPCE water model script I found but I got an error.

or maybe any other advice? A little background, I tried to simulate water adsorption on silica surrounded with argon gas for my undergraduate study. My supervisor told me to simulate the individual components so I also could learn how to use LAMMPS, I have simulated silica, water SPCE and argon, but the water is not a gas. My other supervisor advised me to use fix deposit instead to insert argon and water at the silica surface than make water vapor simulation.

2

To simulate water vapor you have to simulate the correct conditions: you need to set and maintain suitable values for temperature and density.

You need to keep in mind that if you use a specific water potential, that those do not follow exactly the density/temperature/pressure relations that you have from experimental water. As a beginner it is also recommended to stay away from phase transitions. on the other hand, you don’t want to go near or beyond the critical point.

it is impossible to give any specific advice on

Axel

3

Is it enough to model water vapor to increase its temperature to convert it from a fluid phase to vapor?

4

The answer to this question is the same as before:

5

How to adjust the density value? Is it possible to explain more?

6

Density is essentially an input parameter.
You define the volume of the system, you create the atoms/molecules. Thus you control the density of the entire system. Whether you have a vapor or drops of liquid or grains of solid then is determined by the phase diagram. Getting the actual density is a simple exercise in math and thermodynamics.

Same if you read your system from a data file.

7

Dear Axel and all,

I am also trying to simulate the vapor phase of water, starting with 1283 molecules of mW water model in a (20,20,20) nm^3 box, using GCMC and NVT fix commands. While running the simulation, the temperature and pressure increases sharply, with 2 order of magnitude. I have gone through the previous threads and try to fix this issue using 0.1 fs timestep but the issue is unsolved.
I would like to share my input and output files, if required. I would be grateful to you for the same.

Thanks,
Moid

8

A few questions:

  • what happened to the system when the temperature and pressure changed so dramatically when visualizing the trajectory? Is there some corresponding change visible?
  • what happens when you run without fix gcmc? What is the purpose of using fix gcmc in this case?
  • when you change the time step, what happens to the dramatic change event? Is it after the same number of steps or the same time nor neither?
  • what is the exact temperature you are thermostatting to? Is it sufficiently beyond the boiling point of your water model?
9

I have an additional question:

10

Dear Axel,

Thank you very much for your prompt response.
I am observing the following thing while this change occurs in temperature and pressure

  1. The water molecules, which forms droplet (collection of few molecules) get frozen however the vapor molecules move randomly but get arrested once associate with those droplets.

  2. When I am running without fix gcmc, the simulation goes well with NVT simulation

  3. When I am using 0.1 fs time step, it occurs around 16500 step however around 10500 step in the case of 1 time step.

  4. I am running the simulation at 300K which is below the boiling point. I am interested to find out the equilibrium chemical potential at 300K for liquid vapor coexistence.

I just wanted to mention that when we use pressure command in fix gcmc (please see the script below), the chemical potential would be ignored therefore i have run the another simulation without pressure command but the issue has not been resolved.

Dear Simon,

Yes i have used the compute_modify dynamic/dof. Please let me know if you find someting wrong in the script below.

Input script:

variable vol equal 8000000

units real
atom_style full
boundary p p p

region 1 block 0 200 0 200 0 200
read_data water_mW.data

neighbor 2 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 54654

pair_style sw
pair_coeff * * system.in.sw mW

group 1 region 1

mass 1 18.015
fix 1 all momentum 1 linear 1 1 1

timestep 1

min_style cg
minimize 1e-20 1e-20 10000 10000

min_style sd
minimize 1e-40 1e-40 10000 10000

reset_timestep 0

print “================================================”
print “GCMD dynamics for 2000ps”
print “================================================”

fix gcmc1 1 gcmc 1000 100 0 1 99999 300 -10.4 0.5 region 1 pressure 0.1 fugacity_coeff 0.95 group 1 grouptype 1 1 tfac_insert 1

compute t1 1 temp/region 1
compute_modify t1 dynamic/dof yes

variable 1 atom “type==1”
group 3 dynamic all region 1 every 1

compute p1 3 stress/atom t1
compute p1rd3 3 reduce sum c_p1[1] c_p1[2] c_p1[3]
variable press1 equal (-(c_p1rd3[1]+c_p1rd3[2]+c_p1rd3[3])/(3*(vol/3)))*0.101325

fix p1 3 ave/time 1 100 100 v_press1 ave running start 20000 title1 “Average Press[MPa]” file avep0.1_t0.1fs.txt

variable p1ave equal f_p1
variable 3 equal count(3)

thermo_style custom step temp press pe ke density atoms v_3 c_t1 v_press1 v_p1ave
thermo 1000
thermo_modify lost ignore

dump 1 all atom 100000 GCMD_only.dump

run 2000000

unfix gcmc1
undump 1
unfix p1

print “==================================================================”
print “NVT simulation with GCMD”
print “=================================================================”

reset_timestep 0
timestep 0.1

fix gcmc1 1 gcmc 1 20 0 1 99999 300 -10.4 0.5 region 1 pressure 0.1 fugacity_coeff 0.95 group 1 grouptype 1 1 tfac_insert 1
fix 5 1 nvt temp 300.0 300.0 100.0

compute t3 1 temp/region 1
compute_modify t3 dynamic/dof yes

compute p3 3 stress/atom t3
compute p3rd3 3 reduce sum c_p3[1] c_p3[2] c_p3[3]
variable press3 equal (-(c_p3rd3[1]+c_p3rd3[2]+c_p3rd3[3])/(3*(vol/3)))*0.101325

fix p5 3 ave/time 1 100 100 v_press3 ave running start 30000 title1 “Average Press[MPa]” file avep0.1nvt_0.1fs.txt
variable p5ave equal f_p5

thermo_style custom step temp press pe ke density atoms v_3 c_t3 v_press3
thermo 100
thermo_modify lost ignore

dump 3 all atom 1000 GCMD_NVT.dump
run 1000000

unfix 5

Thanks,
Moid

11

There really isn’t much that I can advise beyond some general remarks here. The major open question is whether the kind of system you have set up is actually suitable for the intended purpose. That I cannot answer from remote and even if I had the means to reproduce the simulation, I don’t have the time to go into this since it is nowhere near my interests which are focused on code development and not on applications.

Some observations:

This seems to suggest that you have an oversaturated vapor and the fix gcmc insertions eventually provide the nucleation for condensation. That will release large amounts of kinetic energy and your system may be too small to absorb that. Of course, this is mostly speculation, but something worth investigating.

That is consistent with the hyposthesis above.

I am confused here. This sounds like you didn’t read the documentation for fix gcmc well enough. Specifying “pressure” (of the reservoir) is just an alternate way to specific the chemical potential and thus those must be exclusive. So you just need to use one of the two.

Your input has thermo_modify lost ignore, that is a very bad idea. That way you will not be made aware of you simulation getting into a bad state. Since you have full periodic boundaries, no atoms should be lost.

12

You could try to scan the pressure range by running either a series of fix NVT simulations at different box sizes or fix NPT simulations at different pressures first and observing what happens.

13

Hello,

Can I see your water data file? I think the problem can be seen from there.

Nur.