Dear lammps users,
I’m trying to use fix gcmc for some simulations.
First of all I though to perform some preliminary test with a very simple system composed by just one atom type that can be inserted with the fix gcmc.
I share below the testing script
dimension 3
boundary p p p
units real
atom_style atomic
lattice sc 1
region large_box block -100 100 -100 100 -100 100
create_box 1 large_box
create_atoms 1 random 100 1234 large_box
pair_style lj/cut 20
pair_modify shift yes
pair_coeff * * 10 10
mass 1 1
minimize 1e-6 1e-6 1000 1000
write_data initial_topology.data
velocity all create 300 87287 mom yes rot yes dist gaussian
variable n equal count(all)
neighbor 2 bin
neigh_modify delay 0 every 1
timestep 1
fix 1 all gcmc 1 1 0 1 12345 300 -1 0 region large_box
fix 4 all print 1 "$(step) $(v_n)" file mu_-1_ni_100.txt title "step molecules" screen no
fix 5 all nvt temp 300 300 1
thermo_style custom step v_n temp press etotal
#thermo_modify lost warn
thermo 1
run 100000
Given this I have two questions :
The first is about the fix gcmc command
For what I undestand when I use it I fix also the chemical potential of the imaginary reservoir (in this case -1 kcal/mol) and lammps use this value to compare with the chemical potential inside the region that I gave him.
I tried to run the script above in two cases one time starting from 100 beads and another from 200 (both time of the same type). I was expecting that the simulations would reach the same number of bead after a certain time (reaching like an equilibrium number of bead). But the two simulations reaches very different number of beads. This can be also due to the fact that I should have waited for more time.
The second question is about the chemical potential of the species inside the specified region.
In Lammps is there a way to compute easily the chemica potential inside the region ?
Many thanks if someone can help me I’m using lammps-23Jun2022.
Hope that the questions are clear an thanks for your time.
Best regards
Daniele
Please note the following remark from the fix gcmc documentation:
If used with fix nvt, the temperature of the imaginary reservoir, T, should be set to be equivalent to the target temperature used in fix nvt. Otherwise, the imaginary reservoir will not be in thermal equilibrium with the simulation cell. Also, it is important that the temperature used by fix nvt be dynamic/dof, which can be achieved as follows:
compute mdtemp mdatoms temp
compute_modify mdtemp dynamic/dof yes
fix mdnvt mdatoms nvt temp 300.0 300.0 10.0
fix_modify mdnvt temp mdtemp
The chemical potential is an input parameter for fix gcmc. Please also note this comment from the documentation:
As an alternative to specifying mu directly, the ideal gas reservoir can be defined by its pressure P using the pressure keyword, in which case the user-specified chemical potential is ignored. The user may also specify the fugacity coefficient using the fugacity_coeff keyword, which defaults to unity.
To determine the current chemical potential is a different calculation. You can only determine the excess chemical potential via the fix widom command — LAMMPS documentation. So you have to compute the ideal gas contribution separately.
Sorry for my late reply !!
Many thanks for all your nice advice and your time.
Thanks
I’m sorry to reissue this topic.
I bet that it is a real trivial question but if i try as axel said and as de documentation said with
compute mdtemp mdatoms temp
compute_modify mdtemp dynamic/dof yes
fix mdnvt mdatoms nvt temp 300.0 300.0 10.0
fix_modify mdnvt temp mdtemp
I still don’t reach the same final number of beads starting from different number.
Am I getting something wrong on how fix gcmc should work ?
I’ll attach my code
dimension 3
boundary p p p
units real
atom_style atomic
lattice sc 1
region large_box block -50 50 -50 50 -50 50
create_box 1 large_box
create_atoms 1 random 900 1234 large_box
pair_style lj/cut 20
pair_modify shift yes
pair_coeff * * 10 10
mass 1 1
minimize 1e-6 1e-6 1000 1000
write_data initial_topology.data
velocity all create 300 87287 mom yes rot yes dist gaussian
variable n equal count(all)
neighbor 2 bin
neigh_modify delay 0 every 1
timestep 1
fix 1 all gcmc 1 100 0 1 12345 300 -1 0 region large_box
fix 4 all print 1 "$(step) $(v_n)" file mu_-1_ni_1000.txt title "step molecules" screen no
compute mdtemp all temp
compute_modify mdtemp dynamic/dof yes
fix 5 all nvt temp 300 300 10
fix_modify 5 temp mdtemp
thermo_style custom step v_n temp press etotal
#thermo_modify lost warn
thermo 1
run 100000
Again thanks for your time and sorry for having so many questions.
Best Daniele
Where do these parameters come from? They are very unusual for real units.
For Argon (see examples/UNITS) in real units you have:
pair_style lj/cut 8.76
pair_modify shift yes
pair_coeff * * 0.2349 3.504
mass 1 39.95
With an over 40 times stronger pairwise interaction and 40 times lighter atoms, there is not going to happen much in the MD part and thus it does not seem unreasonable that you won’t reach equilibrium with the number of steps you are using.
thank u for the suggestion.
It does make sense and now seems to work !!
Thanks
You have your problem solved, but you left my curiosity unsatisfied.
Why didn’t you answer my question?
I’m sorry u right I completely forgot.
Those were just random parameters that I put to try to simulate a LJ fluid like.
But when u wrote to me I realized that I should’ve think much more about the importance of those parameters.
I tried with the parameter that u suggest.
But to be sure about the results I also tried with the parameter suggested in Frenkel and Smit, Understanding Molecular Simulation, Academic Press, London, 2002.
To have a reference. The parameters that I took was the ones in the case study 9 at page 133 for lj units.
Using those and setting either the pressure or the chimical potential I was able to recover the density in the figure 5.8.
I’m sorry for my negligence in choosing those parameters I should’ve think more about the phisical meaning.
Many thanks for your time though.
Best
Daniele