Preform a Monte Carlo simulation with fixed molecular number and SHAKE constraints

Hello,

I am going to study ice nucleation by Monte Carlo method. The system is in a canonical (NVT) ensemble with fixed number of water molecules. Also, the SHAKE algorithm is applied to each water molecules.

I am grateful that LAMMPS has developed the MC package for Monte Carlo simulations. But I am not sure which commands (“fix tfmc” or “fix gcmc”) are more appropriate for the problem I am studying.

fix tfmc:
The “fix tfmc” command is not compatible with “fix shake”. The water model I used is TIP4P/2005 or TIP4P/Ice in which water molecules are rigid. If I have to use “fix tfmc”, I need to use a flexible water model with elastic constants of bonds and angles. The flexible water model is more time consuming and not sure to be applicable for TIP4P/2005 or TIP4P/Ice.

fix gcmc:
The “fix gcmc” command is used to grand canonical Monte Carlo simulation with exchange of atoms or molecules. It is compatible to “fix shake”. But I am not sure if I could perform a canonical simulations (fixed number of molecules) by “fix gcmc”.

I am thinking about some solutions as blew:

  1. Set the average number of GCMC exchanges “X” as zero
    That makes “random_int_fraction” be always smaller than or equal to “nmcmoves”. All the “ncycles” are used for molecular translations or rotations.
  2. Set the chemical potential of imaginary reservoir “mu” as zero
  3. Set the probability of GCMC exchange as zero directly in “fix_gcmc.cpp”

I am not sure which solutions works and better. It would be appreciated if you could provide any suggestions.

Best,
Jason

Hello,

I am going to study ice nucleation by Monte Carlo method. The system is in a
canonical (NVT) ensemble with fixed number of water molecules. Also, the
SHAKE algorithm is applied to each water molecules.

I am grateful that LAMMPS has developed the MC package for Monte Carlo
simulations. But I am not sure which commands ("fix tfmc" or "fix gcmc") are
more appropriate for the problem I am studying.

fix tfmc:
The "fix tfmc" command is not compatible with "fix shake". The water model I
used is TIP4P/2005 or TIP4P/Ice in which water molecules are rigid. If I
have to use "fix tfmc", I need to use a flexible water model with elastic
constants of bonds and angles. The flexible water model is more time
consuming and not sure to be applicable for TIP4P/2005 or TIP4P/Ice.

fix gcmc:
The "fix gcmc" command is used to grand canonical Monte Carlo simulation
with exchange of atoms or molecules. It is compatible to "fix shake". But I
am not sure if I could perform a canonical simulations (fixed number of
molecules) by "fix gcmc".

I am thinking about some solutions as blew:
1. Set the average number of GCMC exchanges "X" as zero
That makes "random_int_fraction" be always smaller than or equal to
"nmcmoves". All the "ncycles" are used for molecular translations or
rotations.
2. Set the chemical potential of imaginary reservoir "mu" as zero
3. Set the probability of GCMC exchange as zero directly in "fix_gcmc.cpp"

I am not sure which solutions works and better. It would be appreciated if
you could provide any suggestions.

neither is good. you should use a real MC code.
LAMMPS is still an MD code and the MC package just offers a little "MC
flavoring" to address specific problems.

axel.

Aidan can answer whether what you want is possible and/or efficient with fix gcmc.

Axel is also correct, that it may be that what you want to do will

be faster in a dedicated MC code.

Steve

"study ice nucleation using Monte Carlo." Formally, you can't do
that. Nucleation is a kinetic phenomenon. All kinetic effects in
Monte Carlo simulations are artifical, controlled by arbitrary
parameters like the average displacement per move. That said, it
might be possible to learn something from a Monte Carlo simulation
with physically-inspired kinetics, which is what time-stamped
force-bias Monte Carlo gives you (fix tfmc). I suggest you also look
at the papers of Brown and Yamada who simulated ice nucleation using
straight MD. You could also try the PRD method in LAMMPS.

Dear Axel and Steve,

Thank you so much for your suggestions. Very helpful.

Best,
Jason

Dear Aidan,

Thank you for your suggestions and insights on the modeling of ice nucleation.

In my study, a super-cooled water droplet is placed initially on the substrate of diverse texture, geometry, and surface chemistry. What I focus on is not the growth rate of ice nucleation, but the nucleation path or the final configuration of crystallized ice on different substrates.

The growth rate of ice crystal is reported to be less than 0.1 nm/ns for typical low-energy planes. Also, ice crystal nuclei have minimum size, which limits us to scale down the size of the original droplet. So it is very expensive in computational cost for straight MD of ice nucleation modeling with TIP4P models. I tried droplet of different sizes (2nm, 3nm, 4nm, 5nm), but they don’t start to freeze at 220K after 100ns simulations. Some literature of straight MD reported a droplet of 512 water molecules required over 300ns to be completely frozen. That is why I am interested in a Monte Carlo approach.

I have found a paper by Brown and Yamada as below. Is it the paper you were talking about?

Brown, W. M., & Yamada, M. (2013). Implementing molecular dynamics on hybrid high performance computers—Three-body potentials. Computer Physics Communications, 184(12), 2785-2793.

In Brown and Yamada’s paper, a coarse-grained water model (mW) is used. It reminds us of another paper using mW model for ice nucleation. I will try this coarse-grained model instead of current TIP4P/Ice model. But I am still not very sure the mW model can be used for ice nucleation, at least accepted by reviewers.

I will also try PRD method for acceleration. Thanks again.

Best,
Jason