fix gcmc with shake keyword

Hej all,

is the shake keyword available in fix gcmc yet? Because the LAMMPS version from 21Jan15 results in an error if I use the shake keyword in fix gcmc. I searched for the keyword in fix_gcmc.cpp and found some handling for this but apparently when it comes to setting the shakeflag in FixGCMC::options(...) this is simply not implemented. So I'm asking myself if this is already finished and, however, not correctly or a feature still under construction. Aside from this the keyword is already mentioned in the official LAMMPS documentation.

Best regards
Robert

You can use fix GCMC with shake, but there are some restrictions. First, you have to use the mol keyword in your fix GCMC command. You also have to specify the relevant shake fix that you're using in your fix GCMC command, which it sound like you are doing. If this still doesn't work for you, please let us know which error message you're seeing and send me your input.

Also, I'll note that we're working on a new version of fix GCMC that is more general than the current version. If you'd like to help test this newer version, please let me know.

Paul

Hej Paul,

thanks for helping me working this out! Of course I would like to test the new version. As I said, there is no possibility to activate the shake feature in fix gcmc as the flag is not set in the fix_gcmc.cpp routine. I compared this to the method in fix_deposit.cpp and it seems that there are a few things missing in fix_gcmc.cpp.

Snippet from fix_gcmc.cpp:

void FixGCMC::options(int narg, char **arg)
{
  if (narg < 0) error->all(FLERR,"Illegal fix gcmc command");

  // defaults

  mode = ATOM;
  max_rotation_angle = 10*MY_PI/180;
  regionflag = 0;
  iregion = -1;
  region_volume = 0;
  max_region_attempts = 1000;
  rotation_group = 0;
  rotation_groupbit = 0;
  rotation_inversegroupbit = 0;
  pressure_flag = false;
  pressure = 0.0;
  fugacity_coeff = 1.0;
  shakeflag = 0;
  idshake = NULL;
  
  int iarg = 0;
  while (iarg < narg) {
  if (strcmp(arg[iarg],"mol") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
      imol = atom->find_molecule(arg[iarg+1]);
      if (imol == -1)
        error->all(FLERR,"Molecule template ID for fix gcmc does not exist");
      if (atom->molecules[imol]->nset > 1 && comm->me == 0)
        error->warning(FLERR,"Molecule template for "
                       "fix gcmc has multiple molecules");
      mode = MOLECULE;
      onemols = &atom->molecules[imol];
      nmol = onemols[0]->nset;
      iarg += 2;
    } else if (strcmp(arg[iarg],"region") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
      iregion = domain->find_region(arg[iarg+1]);
      if (iregion == -1)
        error->all(FLERR,"Region ID for fix gcmc does not exist");
      int n = strlen(arg[iarg+1]) + 1;
      idregion = new char[n];
      strcpy(idregion,arg[iarg+1]);
      regionflag = 1;
      iarg += 2;
    } else if (strcmp(arg[iarg],"maxangle") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
      max_rotation_angle = force->numeric(FLERR,arg[iarg+1]);
      max_rotation_angle *= MY_PI/180;
      iarg += 2;
    } else if (strcmp(arg[iarg],"pressure") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
      pressure = force->numeric(FLERR,arg[iarg+1]);
      pressure_flag = true;
      iarg += 2;
    } else if (strcmp(arg[iarg],"fugacity_coeff") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
      fugacity_coeff = force->numeric(FLERR,arg[iarg+1]);
      iarg += 2;
    } else error->all(FLERR,"Illegal fix gcmc command");
  }
}

Apparently there is no setting of the shakeflag. Furthermore there are some infos missing in line 944-945:

    if (shakeflag)
      fixshake->set_molecule(nlocalprev,maxtag_all,imol,NULL,NULL,NULL);

I think NULL is not correct for letting shake know how the molecule is build up. I used com_coord for the first NULL (as mentioned in fix_deposit.cpp) and created a new array for vnew with the random velocity of the new molecule. However, I'm stucking at the problem of the quat word corresponding somehow to a rotation matrix...

Best regards,

Robert

Ok, I'll look into it. I'll also send you the newer (beta) version separately.

Paul