Temperature-control and pressure-control in fix-gcmc

LAMMPS LAMMPS Mailing List Mirror
1

Dear lammps developers and users,

I have problems with temperature-control and pressure-control in fix-gcmc.

I started with argon, the mono-atomic system, for which nvt, npt, nvt+gcmc have been tested. For these three types of simulations, temperature T, pressure P, and density rho, agree with each other very well. Here the main commands in nvt+gcmc, related to T and P, are as follows:

timestep 1.0
fix fixmomentum all momentum 100 linear 1 1 1
compute mdtemp all temp
compute_modify mdtemp dynamic yes
fix mdnvt all nvt temp 140.46 140.46 100.0
fix_modify mdnvt temp mdtemp

fix fixgcmc all gcmc 10 20 10 1 666 140.46 0.0 0.5 pressure 50.0
compute_modify thermo_temp dynamic yes

When I turn to CO2, the molecular system, the results from nvt and npt still agree with each other (also agree with the experimental values).

As for nvt+gcmc, the results seem not correct although I use the same commands above to control T and P, sometimes T and P go extremely high.

To me, the problems may comes from the fluctuation of the number of molecules in the simulation box, and thus temperature and pressure and even some potentials were not calculated correctly (but why it is not a problem for mono-atomic system?). After some tests, I use the following method:

timestep 0.5
compute mdtemp all temp
compute_modify mdtemp dynamic yes

compute_modify thermo_pe dynamic yes
compute_modify thermo_press dynamic yes
compute_modify thermo_temp dynamic yes

fix tc all temp/rescale 1 280 280 10 1
compute_modify tc_temp dynamic yes

fix fixmomentum all momentum 100 linear 1 1 1

fix fixgcmc all gcmc 10 1 1 0 779877 280 -0.0 0.1 maxangle 30 mol CO2 pressure 5.0 tfac_insert 0.5
fix mdnvt all nvt temp 280 280 100
fix_modify mdnvt temp mdtemp

In this way, the temperature fluctuate between 260-300 (set to 280 K), but the pressure has a average about 15 (set to 5 bar), the density fluctuates around 0.015 (the experimental value is 0.02 g/cm^3). Note that I used tail correction and “pair_style lj/cut/coul/cut 10.0”

It seems, I am very close to the correct results but there are still some parameters to be improved (although I have done many tests).

Could you please give me some help with this? Thanks!

Yongbiao

2

Dear lammps developers and users,

[...]

It seems, I am very close to the correct results but there are still some
parameters to be improved (although I have done many tests).

you are not close at all. no simulation, that uses temp/rescale, let
alone uses *two* thermostats for the same set of atoms can be anywhere
close to be correct (in the sense of producing statistical mechanical
meaningful results).

what seems to be the major issue here is a major lack of background
knowledge and realization of some rather fundamental principles.

first off, you have to realize that fix gcmc was not conceived to do
the kind of simulation that you are trying to do (inserting CO2 into
water). it was designed to study distribution and diffusion of gases
in porous media.

with GCMC simulations, you have to realize that insertions are always
more problematic than deletions. obviously so, since a deletion cannot
cause any overlaps, which insertions can. the lower the density of a
medium, the weaker the interactions and the smaller the object to
insert, the easier it is to have successful insertions by just
inserting particles at a random location with a random orientation.
your argon system is a bit more problematic than dilute gases, but -
as you noticed - still manageable. If you have a condensed system,
you'll have to do increasingly more insertion attempts, because you
need to fine a suitable "hole" that can easily accept your inserted
particle. these kinds of "holes" form through thermal fluctuations and
the size of those is limited by the system size. also, "searching" for
a suitable "hole" by just using random coordinates, is not a very
effective way to go about that.

however, water molecules have in comparison quite strong interactions,
thus it is not as likely to find a suitable hole. it is less
compressible than argon, thus forcing in a molecule will increase the
potential energy much more. also, you have additional complications by
having to break hydrogen bonds, for example and finally (and most
importantly), CO2 is large and oblong, so you not only have to find a
suitably sized hole, it also has to have a specific shape and you also
have to identify correspondingly the best orientation to insert the
molecule. common sense tells us, that that is going to be extremely
unlikely to find from just random inserts at random orientations, so
if you set up your simulation to actually generate insertions within a
reasonable time, they are very likely inducing a significant amount of
overlap and thus very high potential energy. the result of which is
what you are trying to "fix" by messing with simulation parameters,
but ultimately, what you are up against, is an intrinsic problem of
the kind of insertion algorithm you are applying and consequently, fix
gcmc in LAMMPS is not the right tool to set up such a simulation. you
need to use a *proper* monte carlo code for this.

in other words, i don't see any immediate solution to solve your
problem through features that currently exist in LAMMPS. you are
trying to do something with fix gcmc for what it is not designed to
do. GI-GO.

axel.

3

Yongbiao,

Axel is too pessimistic. It is possible to simulate CO2 in water using fix gcmc, although it is true that a pure Monte Carlo code (such as TOWHEE) might be more efficient. In any case, the steps you are taking are good. First you simulated a dense CO2 vapor using NVT and NPT MD to establish a reliable benchmark. fix gcmc should be able to reproduce the same thermodynamic state as you measured in NVT and NPT. However, there are a couple of points to remember:

  1. In fix gcmc the pressure keyword is the pressure in a fictitious reservoir of non-interacting CO2 molecules (P = rho k T), not the actual pressure of real CO2 at the same density.
  2. The large fluctuations may be caused by insufficient number of molecules (fluctuations ~ 1/sqrt(N)). Using special fixes like fix momentum and temp/rescale is risky, because they can hide underlying problems.
  3. You could try using fix gcmc with no MD.

Aidan