GCMC + NVT Simulations

Hi all,

I am trying to perform water adsorption in silicate minerals using GCMC module in LAMMPS. My surface contains 688 atoms,10 atom types with 1 bond type and no angles. Therefore, in my lammps data file, i have defined 11 and 12 for water O and H respectively with bond type as 2 and angle type as 1. I have attached the data file for reference.

I will be using rigid SPC water model defined in the molecule file “spc.txt” which contains the template for single water molecule (attached). Initially, i was not able to run GCMC simulations but i think i have passed that part now. However, I have some questions on it.

  1. If i run the GCMC simulations with only the chemical potential (mu) defined and not the pressure, i get low adsorption of H2O molecules. On the other hand, if i define the pressure explicitly thereby ignoring the mu value,i get more H2O in to the system. Does it mean that i had given a low mu value than the one calculated by LAMMPS using the defined pressure. If i am correct the value of

  2. This is also related to 1st question. The temperature fluctuations for simulations performed with defined chemical potential is close to the one i have specified in the input file. But for the simulations, where the mu is computed based on the given pressure values, the temperature is much less the temperature of the reservoir. For instance, it is 50K less than the temperature of the reservoir bath. Is it related to the system getting cold after insertion move. Should i have to use the tfac_insert keyword to over come this issue.

  3. Last, as mentioned earlier, i want to use rigid water model. So i have defined the Shake Flags and related sections in the molecule file. However, when i use fix gcmc command with shake keyword i get an error “Illegal fix gcmc command”. The fix ID for the shake is 1 so i define the fix gcmc as follows

" fix 2 wat gcmc 1 10 10 0 3456543 348.15 -11.07 0.5 mol 99 pressure 88.83 shake 1".

I have attached the complete input files, data and molecule file for your kind attention. Please let me know if you find the attachements broken as i faced similar issues last time.

Thanks in advance,

Regards,

Naresh

spc.txt (591 Bytes)

in.lammps (1.1 KB)

smec.lammps14 (50 KB)

Hi all,
        I am trying to perform water adsorption in silicate minerals using
GCMC module in LAMMPS. My surface contains 688 atoms,10 atom types with 1
bond type and no angles. Therefore, in my lammps data file, i have defined
11 and 12 for water O and H respectively with bond type as 2 and angle type
as 1. I have attached the data file for reference.

I will be using rigid SPC water model defined in the molecule file
"spc.txt" which contains the template for single water molecule (attached).
Initially, i was not able to run GCMC simulations but i think i have passed
that part now. However, I have some questions on it.

1) If i run the GCMC simulations with only the chemical potential (mu)
defined and not the pressure, i get low adsorption of H2O molecules. On the
other hand, if i define the pressure explicitly thereby ignoring the mu
value,i get more H2O in to the system. Does it mean that i had given a low
mu value than the one calculated by LAMMPS using the defined pressure. If i
am correct the value of

You can set mu and/or pressure to whatever you like. In each case, the
system will equilibrate to a state with a chemical potential matching what
you specified. The idea that you get more adsorption when using pressure
than when using mu makes as much sense as saying that Fahrenheit
temperatures are hotter than Celsius temperatures.

2) This is also related to 1st question. The temperature fluctuations for
simulations performed with defined chemical potential is close to the one i
have specified in the input file. But for the simulations, where the mu is
computed based on the given pressure values, the temperature is much less
the temperature of the reservoir. For instance, it is 50K less than the
temperature of the reservoir bath. Is it related to the system getting cold
after insertion move. Should i have to use the tfac_insert keyword to over
come this issue.

I we assume that your molecule template is in a configuration that is a lot
lower in energy than the average energy of the molecule at the temperature
of interest, then yes, adding a lot of molecules will cool the system down.
This can be counteracted using the tfac_insert keyword.

3) Last, as mentioned earlier, i want to use rigid water model. So i have
defined the Shake Flags and related sections in the molecule file. However,
when i use fix gcmc command with shake keyword i get an error "Illegal fix
gcmc command". The fix ID for the shake is 1 so i define the fix gcmc as
follows

" fix 2 wat gcmc 1 10 10 0 3456543 348.15 -11.07 0.5 mol 99 pressure 88.83
shake 1".

You must be running the wrong input file or something, because your script
works fine for me and it contains the same line that you show above:

      fix 2 wat gcmc 1 10 10 0 3456543 348.15 -11.07 0.5 mol 99 pressure
88.83 shake 1

Aidan

it could also be due to using an older version of LAMMPS that is
missing support for shake and molecule templates.

axel.

Hi,

Thanks for the reply. Yes. I was using 2014 version and i updated to the recent one. It is running now but still ends in error.
ERROR on proc 0: Shake determinant = 0.0 (…/fix_shake.cpp:1998). I get this error at around 200 fs. Please let me know if you think anything not correct in the molecule file.

Also, i have one other question. As mentioned earlier, i am looking to get the adsorption of H2O in silicate minerals at a interlayer distance of 12.0. However, i see that the layers are expanding because i did not keep them as rigid.

Since, the interlayer water content significantly depends on the relative displacement of one layer over the other, i kept the layers flexible. I think this leads to the reason for expansion or contraction during minimization process.

Therefore, i am thinking to allow the silicate layer flexible along lateral (X and Y directions) but keeping the Z coordinate fixed (so that i can keep the interlayer distance fixed).

I looked in to the manual and found the keyword “enforce2d”. When used, i end up getting the error message saying “ERROR: Cannot use fix enforce2d with 3d simulation”.

There are fixs like “rigid” and “oneway”. But i did could not a way to enforce the movement of group of particles along xy and not z directions.

Please give me your suggestions.

Thanks alot,

Naresh

“ERROR on proc 0: Shake determinant = 0.0” This indicates that the SHAKE matrix is singular. This should never happen for a physical simulation. Probably something else bad is happening long before you see this error message. Check to see that all your velocities and forces are in the normal range. If not, figure out why that is.

If you read the documentation, you would have seen right away that enforce2d has nothing to do with what you are looking for. There are many ways to stop things from moving in LAMMPS, see e.g fix spring and setforce commands.

Aidan

Hi,
   Thanks for the reply. Yes. I was using 2014 version and i updated to the
recent one. It is running now but still ends in error.
   ERROR on proc 0: Shake determinant = 0.0 (../fix_shake.cpp:1998). I get
this error at around 200 fs. Please let me know if you think anything not
correct in the molecule file.

this is more likely a consequence of bad force field parameters or a
bad geometry.
please visualize your trajectory. you will likely see one (or more)
atoms moving at very high velocity.

   Also, i have one other question. As mentioned earlier, i am looking to
get the adsorption of H2O in silicate minerals at a interlayer distance of
12.0. However, i see that the layers are expanding because i did not keep
them as rigid.

please note, that "rigid" is not the correct term for what you are
describing, but "immobile".

   Since, the interlayer water content significantly depends on the relative
displacement of one layer over the other, i kept the layers flexible. I
think this leads to the reason for expansion or contraction during
minimization process.
   Therefore, i am thinking to allow the silicate layer flexible along
lateral (X and Y directions) but keeping the Z coordinate fixed (so that i
can keep the interlayer distance fixed).

   I looked in to the manual and found the keyword "enforce2d". When used, i
end up getting the error message saying "ERROR: Cannot use fix enforce2d
with 3d simulation".
  There are fixs like "rigid" and "oneway". But i did could not a way to
enforce the movement of group of particles along xy and not z directions.

you can try the following:

velocity layergroup set NULL NULL 0.0
fix f1 layergroup setforce NULL NULL 0.0

if a particle has no z-direction velocity and not z-direction force,
it will not move in that direction.

axel.