GCMC + NVT Simulations with SHAKE

Hello LAMMPS users,

I am trying to run a gcmc simulation for SPC/E water to obtain water adsorption isotherms in carbon nanotube.
At first, I put CNT in the center of simulation box.
Then, I inserted 10 water molecules randomly into the simulation box.
The configuration seemed to be fine after minimization, so I ran a GCMC + NVT simulation with fix shake.

Here are some lines that use fix shake and shake keyword in fix gcmc, combined with fix nvt to simulate the grand canonical ensemble.

fix mynvt all nvt temp \{temp\} {temp} 100
fix wshake h2o shake 0.0001 50 1 b 2 a 2 mol h2omol
fix mygcmc h2o gcmc 100 100 0 0 54341 \{temp\} {mu} \{disp\} mol h2omol tfac\_insert {tfac} group h2o shake wshake

However, LAMMPS stopped with below error on step 3.

ERROR on proc 0: Out of range atoms - cannot compute PPPM

So, I checked the movie and SHAKE stats.
Then, I noticed that 2 H atoms in a water molecule overlapped although I used fix shake command for water molecule.
Below is the output for SHAKE stats.

SHAKE stats (type/ave/delta) on step 3
  1 1.26653 0.0600891 18 :bond
  1 0.8368 0.488548 :angle

I think fix shake command doesn't work.

Instead, when the shake keyword is excluded from the fix gcmc:

fix mygcmc h2o gcmc 100 100 0 0 54341 \{temp\} {mu} \{disp\} mol h2omol&                  tfac\_insert {tfac} group h2o ### shake wshake

the LAMMPS run normally.

How can I solve this problem?
I attach my input script.
I'll be pleased to attach more input scripts if needed.
Any help would be very much appreciated.

Best Regards,

Yuki

in.GCMC.CNT (2.23 KB)

what is your LAMMPS version?

Dear Mr.Axel,

Thank you for your reply.
My LAMMPS version was 16 Mar 2018.
Now, I have updated to the latest version (22 Aug 2018) and run the same simulation again.
Although, the situation does not change.

Best regards,
Yuki

there is not much to suggest here. fix gcmc is tricky and LAMMPS is
not really a code designed for MC applications.
one thing you could try is to use pair styles from the CORESHELL
package instead of the regular styles. those have a small modification
to the coulomb interactions, that avoids the singularity for
overlapping atoms. not sure if that helps, but since it is easy to do,
it is worth a try.

axel.

Dear Mr. Axel,

I have changed pair styles (lj/cut/coul/long -> lj/cut/coul/long/cs) and run the same simulation again.
However, the situation does not change.
Anyway, thank you for your advice.

Best regards,
Yuki

Dear Mr. Axel,

I have changed pair styles (lj/cut/coul/long -> lj/cut/coul/long/cs) and run the same simulation again.
However, the situation does not change.

ok. there are two more options, that you could try: 1) you could try
using fix rigid instead of fix shake, 2) you could try a flexible
water potential

please keep in mind, that fix gcmc was initially conceived to study
noble gas diffusion in zeolithes. it has been expanded in scope
multiple times and what you are trying is the most complex combination
of options (molecules vs. atoms, long ranging coulomb vs. short-range
LJ interactions, rigid vs. flexible molecules).

so trying to avoid any of those complications might help. but
ultimately, it may just be that LAMMPS isn't a good tool for your
study.

axel.

Dear Axel,

Thank you for your reply.
I have decided to use fix rigid and the simulation runs successfully now.

OK. I will keep it in mind.
Thank you for your advises.

Best regards,
Yuki