Recently, I have been using GCMC to simulate the adsorption of water molecules on silicon surfaces. But we encountered some problems.
Firstly, the pressure in the water molecule region is significantly higher than the pressure set in the fix gcmc command. Secondly, under my operating conditions, water molecules in the gas phase region agglomerate, resulting in lower pressure and higher density of water molecules in the gas phase. (The temperature of the system is 500K)
这是我的in文件
units metal
dimension 3
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
region box block -53.76513 49.9240 -53.21274 53.21274 -214.866 214.866 units box
create_box 6 box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
molecule sder H2O.txt
lattice diamond 5.431 origin 0 0 0 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region 1 block -53.76513 49.9240 -53.21274 53.21274 205.45923 214.866 units box
create_atoms 3 region 1
mass 3 28.09
group wall1 type 3
lattice diamond 5.431 origin 0 0 0 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region 2 block -53.76513 49.9240 -53.21274 53.21274 186.64569 205.45923 units box
create_atoms 4 region 2
mass 4 28.09
group cu1 type 4
region 3 block -48 48 -53 53 -176 176 units box
create_atoms 0 random 2000 33345 3 mol sder 2367 units box
mass 1 15.9994
mass 2 1.008
group ar type 1 2
lattice diamond 5.431 origin 0 0 0 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region 4 block -53.76513 49.9240 -53.21274 53.21274 -214.866 -205.45923 units box
create_atoms 5 region 4
mass 5 28.09
group wall2 type 5
lattice diamond 5.431 origin 0 0 0 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region 5 block -53.76513 49.9240 -53.21274 53.21274 -205.45922 -186.64569 units box
create_atoms 6 region 5
mass 6 28.09
group cu2 type 6
kspace_style pppm 1.0e-6
kspace_modify order 6
fix 1111 wall1 nve/noforce
fix 1112 wall2 nve/noforce
pair_style hybrid lj/cut 12.0 sw lj/cut/coul/long 10
pair_coeff 1 1 lj/cut/coul/long 0.006739 3.166
pair_coeff 2 2 lj/cut/coul/long 0.0000 0.0000
pair_coeff 2 3 lj/cut 0.0000 0.00000
pair_coeff 2 4 lj/cut 0.0000 0.00000
pair_coeff 2 5 lj/cut 0.0000 0.00000
pair_coeff 2 6 lj/cut 0.0000 0.00000
pair_coeff 1 2 lj/cut/coul/long 0.00000 0.00000
pair_coeff 1 3 lj/cut 0.035646 2.6305
pair_coeff 1 4 lj/cut 0.035646 2.6305
pair_coeff 1 5 lj/cut 0.035646 2.6305
pair_coeff 1 6 lj/cut 0.035646 2.6305
pair_coeff 3 4 lj/cut 2.095 2.168
pair_coeff 3 5 lj/cut 2.095 2.168
pair_coeff 3 6 lj/cut 2.095 2.168
pair_coeff 4 5 lj/cut 2.095 2.168
pair_coeff 4 6 lj/cut 2.095 2.168
pair_coeff 5 6 lj/cut 2.095 2.168
pair_coeff * * sw Si.sw NULL NULL Si Si Si Si
bond_coeff 1 100.00 1
angle_coeff 1 300.00 109.47
group cu union cu1 cu2
group wall union wall1 wall2
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
region GCMC block -45 45 -45 45 -150 150 units box
velocity ar create 500 4928495 rot yes dist gaussian
velocity cu create 500 49995 rot yes dist gaussian
velocity wall set 0.0 0.0 0.0 units box
################ minimize
min_style cg
minimize 1.0e-4 1.0e-6 10000 10000
#output
timestep 0.001
dump 1 all xyz 10000 dump8500.xyz
fix 2 wall setforce 0.0 0.0 0.0
write_data ab.data
compute 1 all pe/atom
compute 2 ar ke/atom
compute 22 cu1 ke/atom
compute 3 ar cluster/atom 5.0
compute peratom ar stress/atom NULL
compute p ar reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
fix 33 ar ave/atom 10 20 1000 c_peratom[1] c_peratom[2] c_peratom[3]
compute hk ar chunk/atom bin/1d z lower 0.5
fix 23 ar ave/chunk 10000 1 10000 hk f_33[*] file press8500.profile
variable press equal -(c_p[1]+c_p[2]+c_p[3])
compute pe all reduce sum c_1
dump myDump all atom 10000 dump8500.atom
thermo 10000
thermo_style custom step temp pe etotal press v_press c_pe
thermo_modify lost ignore flush yes
variable tempar atom c_2
variable tempcu1 atom c_22
compute h ar chunk/atom bin/1d z lower 2
fix 3 ar ave/chunk 10000 1 10000 h density/mass file dens8500.profile
fix 4 ar ave/chunk 1000 10 10000 h v_tempar file temp8500.profile
compute cluster ar cluster/atom 5.0
compute cc1 ar chunk/atom c_cluster compress yes
compute size ar property/chunk cc1 count
fix 5 ar ave/histo 10000 1 10000 0 30 30 c_size mode vector ave running beyond ignore &
file tmp8500.histo
compute hhh cu1 chunk/atom bin/3d x lower 1 y lower 1 z lower 1
fix 7 cu1 ave/chunk 1000 10 10000 hhh v_tempcu1 file temp8500.profile
dump 2 all custom 10000 dump8500.lammpstrj id type mol x y z vx vy vz fx fy fz
dump 3 all custom 10000 dump8500.mtForce.* id type x y z c_3 c_1 &
c_peratom[1] c_peratom[2] c_peratom[3] fx fy fz
dump_modify 3 pad 5
fix 22 cu nvt temp 500.0 500.0 0.1
fix 22222 ar nvt temp 500.0 500.0 0.1
fix 10 ar shake 0.0001 20 0 b 1 a 1 mol sder
compute_modify thermo_temp dynamic yes
compute_modify 22_temp dynamic yes
compute_modify 22222_temp dynamic yes
run_style verlet
run 200000
unfix 22222
variable tfac equal 5.0/3.0
#(3 trans + 2 rot)/(3 trans)
fix 55322 ar gcmc 40000 1000 0 0 45321 500 -3 0.5 mol sder shake 10 region GCMC pressure 4.0 group ar full_energy tfac_insert ${tfac}
fix myar ar nvt temp 500.0 500.0 0.1
compute_modify myar_temp dynamic yes
thermo 10000
thermo_style custom step temp pe etotal press v_press c_pe
thermo_modify lost ignore flush yes
run 1000000
unfix 55322
unfix myar
fix my2 ar nvt temp 500.0 500.0 0.1
compute_modify my2_temp dynamic yes
run_style verlet
run 500000
Hello,
I don’t really understand what pressure you are referring to. Also, if you have interfaces and solid material in your system, do you really expect the pressure in the (GCMC) reservoir and the pressure in the water near the solid to be the same? The chemical potentials should be the same at equilibrium, but the pressure?
Also you said that there is an agglomeration in the “gas phase”. What about testing whether you see the same behavior in a box containing only water vapor. What is expected for pure water system is quite well understood and documented, and it is easy to make sure that the GCMC simulation returns the correct behavior in this particular case.
Simon
I’m very sorry, I couldn’t describe my problem very well. Due to the inability to add attachments, you can check the link I provided.
https://pan.baidu.com/s/1eRCYgHPB2no2b9UWojUbkA?pwd=1234
Here are some explanations of my program. The working fluids I used were water and silicon, and the simulated box sizes were x (-53,53), y (-53,53), and z (-215，215). The upper and lower regions are silicon substrates, and the middle region is water molecules.You can check the image called snapshot in my link.
I added water molecules using GCMC in the box z (-150，150) region. The pressure I am referring to is the pressure in the X (EDGE, EDGE), Y (EDGE, EDGE), Z (-150150) region, which is the gas phase pressure.
I calculated the pressure and density in the gas phase region Z (-150，150). You can view the pressure density image in the link. It can be observed that as the gas phase pressure increases, the aggregation of water molecules causes my gas phase pressure to be lower than the standard pressure in the NIST library.
Out of curiosity, what is your units system, I don’t see it being defined in your input, and units lj is the default (edit: the lj units system would be a very bad choice given the values of your parameters).
Thank you for taking the time to reply to my question. I am sorry about that，I missed the unit while copying the in file
the units I choose is metal, The water molecule model is spc/e
OK.
Well in any case if you want to test the robustness of an approach, you have to simplify and compare your results with known experiments/theoretical values. In this case, I would definitely make pure vapor simulations. Once you are 100% sure of your GCMC approach in vapor, you can come back to this more complex system.
There seems to be a misunderstanding here. The pressure parameter is a different way to enter the chemical potential (the more common way to define the probability of a GCMC run to insert or remove atoms for the given setup). Please see the discussion of the pressure keyword in the fix gcmc documentation. It is the pressure in a fictitious gas reservoir that provides the atoms to insert and does not refer to the gas phase part of your system.
Thank you for taking the time to reply to my question.
I know that the pressure in the gas reservoir and the pressure in the GCMC area are not the same value, and what confuses me is, is there any relationship between the two? The pressure of the gas reservoir I specified is 4.0 bars, and the pressure of the system I want is about 0-2.6 MPa. The pressure of the gas reservoir is much lower than the pressure of the system. Why can the gas reservoir still add molecules to the system so quickly.
If I want to keep the final pressure of my system at 2.0MPa, how should I specify the pressure in the fix gcmc command?
Thank you for your reply. I will go test it now
Whether fix gcmc inserts atoms depends on the change in energy caused by the insertion. Like with any Monte Carlo method a lowering of energy will usually result in an accepted move and a move raising the energy should be accepted with a probability determined by a Boltzman factor. Similar for MC rotations and translations. How the chemical potential plays into that is something that you should learn from a text book on the GCMC method and is beyond the scope of this forum. Typically there will be a “critical” value of “mu”; if you stray from that just a small bit, the system will either be depleted or filled rather quickly.