GCMC lammps-version after lammps-30Jan15

After my first trials to reproduce a \mu vs \rho from
dx.doi.org/10.1021/jp302428b (Fig. 10 left), I'd like
to ask some questions. Below are some only outlining
remarks, I'd be glad to communicate further.

1) the said paper's simulations have been made with LAMMPS,
the corresponding author did, to my knowledge, provide
contributions to the GCMC development in LAMMPS. In the
paper, it is said: "all simulations were made with long
range (pppm) electrostatics".
    I started out w/pre lammps-30Jan15 (serial) version and
had immediate success, by using a SPC model with kspace
(pair_style lj/cut/coul/long 10.0, kspace_style pppm 1.0e-5),
shake on the water molecules, starting with a (large enough)
box containing 1 molecule and providing molecules to
be gcmc-inserted from a template, I was able to trace
this curve surprisingly exact.

2) after trying lammps-30Jan15, this approach would no
longer work. kspace-solver will now bail out immediately with
"ERROR: Cannot use kspace solver on system with no charge (../kspace.cpp:282)"
    OK, it has been warned about kspace in the GCMC documentation
anyway, but what now? To my big surprise, when experimenting w/
pair_style lj/cut/coul/cut 10.0 in connection with "dielectric 78"
(it's a diluted system, \rho between 8x10^-4 and 5x10^-2 kg/m³)
the values were, again, almost accurately hit. This means,
the kspace_style pppm 1.0e-5 in the old versions had the effect
of "dielectric 78" in non-kspace systems (rc = 10A, both program
versions). I'm not sure how to handle this further, therefore
my inquiry.

3) any attempts to combine fix gcmc w/fic nvt on water molecules
failed on "error: missing bond xx to yy" when invoking nvt after
gcmc. Outputting coordinates plus atom id's reveals, id's are
constantly increasing during the gcmc steps, therefore, they
are not continous when entering the nvt steps. I'm not sure
this is the reason but my attempts to check trough fix_gcmc.cpp
(maybe atom id's can be re-sorted dynamically) were not fruitful.
Would such an (hybrid mc/md) approach be generally possible
at this point?

4) I'd be glad to contribute in the solution of this problem
(after eventually understanding the intricacies).

Regards & Thanks


1) I hadn't seen that paper, but one of the authors (Sabine Leroch, CC'd) has posted on this list in the past about GCMC. She may wish to comment further, but my understanding is that she used her own independent GCMC implementation to perform these calculations. The standard version of fix GCMC distributed from the LAMMPS website prior to Jan 30, 2015 was not written to work with kspace, and the documentation made that clear. It may have produced results, but those results would not have included the long-range electrostatics (kspace) contributions. The new version as of Jan 30th now allows use of kspace, and the current documentation gives details. The long-range electrostatics (kspace) contributions are included.

2) The error message you're seeing is because at some point in your simulation, you have zero point charges in your system. There are several ways around this, but the quick work-around would be to make sure you never have zero point charges in your system. The better fix, which we'll need to do, is enable LAMMPS kspace solvers to work on systems with no charge. Your experimentation has likely demonstrated fortuitous agreement between setups that all neglect or dampen the long-range electrostatics.

3) You should be able to use the current version of fix GCMC with fix NVT. Please make sure you're following the instructions in the fix GCMC documentation in this regard, especially with respect to the compute_modify command. If you still see errors, please send me an example problem and I'll take a look.

4) Thanks for your post and observations.