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