Unverified Commit 4d4ae93e authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add missing group names in fix widom. refactor group definition in fix widom and gcmc

parent 93ed07f4
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -75,7 +75,7 @@ eim: NaCl using the EIM potential
ellipse:  ellipsoidal particles in spherical solvent, 2d system
flow:     Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
gcmc:     Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
gcmc:     Grand Canonical MC with fix gcmc, Widom insertion with fix widom
gjf:      use of fix langevin Gronbech-Jensen/Farago option
granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics
+10 −29
Original line number Diff line number Diff line
@@ -493,11 +493,7 @@ void FixGCMC::init()
    }
  }

  if (full_flag) {
    char *id_pe = (char *) "thermo_pe";
    int ipe = modify->find_compute(id_pe);
    c_pe = modify->compute[ipe];
  }
  if (full_flag) c_pe = modify->compute[modify->find_compute("thermo_pe")];

  int *type = atom->type;

@@ -580,18 +576,12 @@ void FixGCMC::init()
  // skip if already exists from previous init()

  if (full_flag && !exclusion_group_bit) {
    char **group_arg = new char*[4];

    // create unique group name for atoms to be excluded

    int len = strlen(id) + 30;
    group_arg[0] = new char[len];
    sprintf(group_arg[0],"FixGCMC:gcmc_exclusion_group:%s",id);
    group_arg[1] = (char *) "subtract";
    group_arg[2] = (char *) "all";
    group_arg[3] = (char *) "all";
    group->assign(4,group_arg);
    exclusion_group = group->find(group_arg[0]);
    auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id;
    group->assign(group_id + " subtract all all");
    exclusion_group = group->find(group_id);
    if (exclusion_group == -1)
      error->all(FLERR,"Could not find fix gcmc exclusion group ID");
    exclusion_group_bit = group->bitmask[exclusion_group];
@@ -603,34 +593,25 @@ void FixGCMC::init()
    char **arg = new char*[narg];;
    arg[0] = (char *) "exclude";
    arg[1] = (char *) "group";
    arg[2] = group_arg[0];
    arg[2] = (char *) group_id.c_str();
    arg[3] = (char *) "all";
    neighbor->modify_params(narg,arg);
    delete [] group_arg[0];
    delete [] group_arg;
    delete [] arg;
  }

  // create a new group for temporary use with selected molecules

  if (exchmode == EXCHMOL || movemode == MOVEMOL) {
    char **group_arg = new char*[3];

    // create unique group name for atoms to be rotated
    int len = strlen(id) + 30;
    group_arg[0] = new char[len];
    sprintf(group_arg[0],"FixGCMC:rotation_gas_atoms:%s",id);
    group_arg[1] = (char *) "molecule";
    char digits[12];
    sprintf(digits,"%d",-1);
    group_arg[2] = digits;
    group->assign(3,group_arg);
    molecule_group = group->find(group_arg[0]);

    auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id;
    group->assign(group_id + " molecule -1");
    molecule_group = group->find(group_id);
    if (molecule_group == -1)
      error->all(FLERR,"Could not find fix gcmc rotation group ID");
    molecule_group_bit = group->bitmask[molecule_group];
    molecule_group_inversebit = molecule_group_bit ^ ~0;
    delete [] group_arg[0];
    delete [] group_arg;
  }

  // get all of the needed molecule data if exchanging
+11 −37
Original line number Diff line number Diff line
@@ -298,13 +298,7 @@ void FixWidom::init()
    }
  }

  if (full_flag) {
    char *id_pe = (char *) "thermo_pe";
    int ipe = modify->find_compute(id_pe);
    c_pe = modify->compute[ipe];
  }

  int *type = atom->type;
  if (full_flag) c_pe = modify->compute[modify->find_compute("thermo_pe")];

  if (exchmode == EXCHATOM) {
    if (nwidom_type <= 0 || nwidom_type > atom->ntypes)
@@ -340,19 +334,14 @@ void FixWidom::init()
  // skip if already exists from previous init()

  if (full_flag && !exclusion_group_bit) {
    char **group_arg = new char*[4];

    // create unique group name for atoms to be excluded

    int len = strlen(id) + 30;
    group_arg[0] = new char[len];
    group_arg[1] = (char *) "subtract";
    group_arg[2] = (char *) "all";
    group_arg[3] = (char *) "all";
    group->assign(4,group_arg);
    exclusion_group = group->find(group_arg[0]);
    auto group_id = std::string("FixWidom:widom_exclusion_group:") + id;
    group->assign(group_id + " subtract all all");
    exclusion_group = group->find(group_id);
    if (exclusion_group == -1)
      error->all(FLERR,"Could not find fix Widom exclusion group ID");
      error->all(FLERR,"Could not find fix widom exclusion group ID");
    exclusion_group_bit = group->bitmask[exclusion_group];

    // neighbor list exclusion setup
@@ -362,33 +351,23 @@ void FixWidom::init()
    char **arg = new char*[narg];;
    arg[0] = (char *) "exclude";
    arg[1] = (char *) "group";
    arg[2] = group_arg[0];
    arg[2] = (char *) group_id.c_str();
    arg[3] = (char *) "all";
    neighbor->modify_params(narg,arg);
    delete [] group_arg[0];
    delete [] group_arg;
    delete [] arg;
  }

  // create a new group for temporary use with selected molecules

  if (exchmode == EXCHMOL) {
    char **group_arg = new char*[3];
    // create unique group name for atoms to be rotated
    int len = strlen(id) + 30;
    group_arg[0] = new char[len];
    group_arg[1] = (char *) "molecule";
    char digits[12];
    sprintf(digits,"%d",-1);
    group_arg[2] = digits;
    group->assign(3,group_arg);
    molecule_group = group->find(group_arg[0]);

    auto group_id = std::string("FixWidom:rotation_gas_atoms:") + id;
    group->assign(group_id + " molecule -1");
    molecule_group = group->find(group_id);
    if (molecule_group == -1)
      error->all(FLERR,"Could not find fix Widom rotation group ID");
      error->all(FLERR,"Could not find fix widom rotation group ID");
    molecule_group_bit = group->bitmask[molecule_group];
    molecule_group_inversebit = molecule_group_bit ^ ~0;
    delete [] group_arg[0];
    delete [] group_arg;
  }

  // get all of the needed molecule data if exchanging
@@ -836,9 +815,6 @@ void FixWidom::attempt_molecule_insertion_full()
  MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);

  for (int imove = 0; imove < ninsertions; imove++) {

    int nlocalprev = atom->nlocal;

    double com_coord[3];
    if (regionflag) {
      int region_attempt = 0;
@@ -1051,8 +1027,6 @@ double FixWidom::molecule_energy(tagint gas_molecule_id)

double FixWidom::energy_full()
{
  int imolecule;

  if (triclinic) domain->x2lamda(atom->nlocal);
  domain->pbc();
  comm->exchange();