Commit 961b9763 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add/correct code to define rigid bodies with atomstyle/atomfile variables and... 

add/correct code to define rigid bodies with atomstyle/atomfile variables and via custom per-atom properties
parent ac6434e4

src/RIGID/fix_rigid.cpp

+9 −4
Original line number Diff line number Diff line
@@ -137,9 +137,10 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : 
    int custom_flag = strcmp(arg[3],"custom") == 0;

    if (custom_flag) {

      if (narg < 5) error->all(FLERR,"Illegal fix rigid command");



      // determine whether atom-style variable or atom property is used.

      if (strstr(arg[4],"i_") == arg[4]) {

        int is_double;

        int is_double=0;

        int custom_index = atom->find_custom(arg[4]+2,is_double);

        if (custom_index == -1)

          error->all(FLERR,"Fix rigid custom requires previously defined property/atom");
@@ -153,13 +154,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : 
        MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world);

        molecule = new tagint[nlocal];

        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit) molecule[i] = (tagint)(value[i] - minval + 1);

          if (mask[i] & groupbit)

            molecule[i] = (tagint)(value[i] - minval + 1);

          else

            molecule[i] = 0;



      } else if (strstr(arg[4],"v_") == arg[4]) {

        int ivariable = input->variable->find(arg[4]+2);

        if (ivariable < 0)

          error->all(FLERR,"Variable name for fix rigid custom does not exist");

        if (input->variable->atomstyle(ivariable) == 0)

          error->all(FLERR,"Fix rigid custom variable is not atom-style variable");

          error->all(FLERR,"Fix rigid custom variable is no atom-style variable");

        double *value = new double[nlocal];

        input->variable->compute_atom(ivariable,0,value,1,0);

        int minval = INT_MAX;
@@ -169,7 +174,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : 
        MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world);

        molecule = new tagint[nlocal];

        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit) molecule[i] = (tagint)((int)value[i] - minval + 1);

          if (mask[i] & groupbit) molecule[i] = (tagint)((tagint)value[i] - minval + 1);

        delete[] value;

      } else error->all(FLERR,"Unsupported fix rigid custom property");

    } else {

src/RIGID/fix_rigid_small.cpp

+80 −30
Original line number Diff line number Diff line
@@ -29,7 +29,9 @@ 
#include "group.h"

#include "comm.h"

#include "force.h"

#include "input.h"

#include "output.h"

#include "variable.h"

#include "random_mars.h"

#include "math_const.h"

#include "memory.h"
@@ -67,8 +69,9 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 
  infile(NULL), body(NULL), bodyown(NULL), bodytag(NULL), atom2body(NULL),

  xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL),

  avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL),

  itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL), id_dilate(NULL), 

  onemols(NULL), hash(NULL), bbox(NULL), ctr(NULL), idclose(NULL), rsqclose(NULL)

  itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL),

  id_dilate(NULL), onemols(NULL), hash(NULL), bbox(NULL), ctr(NULL),

  idclose(NULL), rsqclose(NULL)

{

  int i;
@@ -89,7 +92,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 
  // perform initial allocation of atom-based arrays

  // register with Atom class



  extended = orientflag = dorientflag = 0;

  extended = orientflag = dorientflag = customflag = 0;

  bodyown = NULL;

  bodytag = NULL;

  atom2body = NULL;
@@ -103,24 +106,70 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 


  // parse args for rigid body specification



  if (narg < 4) error->all(FLERR,"Illegal fix rigid/small command");

  if (strcmp(arg[3],"molecule") != 0)

    error->all(FLERR,"Illegal fix rigid/small command");

  int *mask = atom->mask;

  tagint *bodyid = NULL;

  int nlocal = atom->nlocal;



  if (narg < 4) error->all(FLERR,"Illegal fix rigid/small command");

  if (strcmp(arg[3],"molecule") == 0) {

    if (atom->molecule_flag == 0)

      error->all(FLERR,"Fix rigid/small requires atom attribute molecule");

  if (atom->map_style == 0)

    error->all(FLERR,"Fix rigid/small requires an atom map, see atom_modify");



  // maxmol = largest molecule #

  } else if (strcmp(arg[3],"custom") == 0) {

    if (narg < 5) error->all(FLERR,"Illegal fix rigid/small command");

      bodyid = new tagint[nlocal];

      customflag = 1;



      // determine whether atom-style variable or atom property is used.

      if (strstr(arg[4],"i_") == arg[4]) {

        int is_double=0;

        int custom_index = atom->find_custom(arg[4]+2,is_double);

        if (custom_index == -1)

          error->all(FLERR,"Fix rigid/small custom requires previously defined property/atom");

        else if (is_double)

          error->all(FLERR,"Fix rigid/small custom requires integer-valued property/atom");



        int minval = INT_MAX;

        int *value = atom->ivector[custom_index];

        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit) minval = MIN(minval,value[i]);

        int vmin = minval;

        MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world);



  int *mask = atom->mask;

  tagint *molecule = atom->molecule;

  int nlocal = atom->nlocal;

        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit) 

            bodyid[i] = (tagint)(value[i] - minval + 1);

          else bodyid[i] = 0;



      } else if (strstr(arg[4],"v_") == arg[4]) {

        int ivariable = input->variable->find(arg[4]+2);

        if (ivariable < 0)

          error->all(FLERR,"Variable name for fix rigid/small custom does not exist");

        if (input->variable->atomstyle(ivariable) == 0)

          error->all(FLERR,"Fix rigid/small custom variable is no atom-style variable");

        double *value = new double[nlocal];

        input->variable->compute_atom(ivariable,0,value,1,0);

        int minval = INT_MAX;

        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit) minval = MIN(minval,(int)value[i]);

        int vmin = minval;

        MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world);



        for (i = 0; i < nlocal; i++)

          if (mask[i] & groupbit)

            bodyid[i] = (tagint)((tagint)value[i] - minval + 1);

          else bodyid[0] = 0;

        delete[] value;

      } else error->all(FLERR,"Unsupported fix rigid custom property");

  } else error->all(FLERR,"Illegal fix rigid/small command");



  if (atom->map_style == 0)

    error->all(FLERR,"Fix rigid/small requires an atom map, see atom_modify");



  // maxmol = largest bodyid #

  maxmol = -1;

  for (i = 0; i < nlocal; i++)

    if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);

    if (mask[i] & groupbit) maxmol = MAX(maxmol,bodyid[i]);



  tagint itmp;

  MPI_Allreduce(&maxmol,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
@@ -155,6 +204,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 
  }



  int iarg = 4;

  if (customflag) ++iarg;



  while (iarg < narg) {

    if (strcmp(arg[iarg],"langevin") == 0) {

      if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid/small command");
@@ -344,11 +395,12 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 
  if (pcouple == XYZ || (domain->dimension == 2 && pcouple == XY)) pstyle = ISO;

  else pstyle = ANISO;



  // create rigid bodies based on molecule ID

  // create rigid bodies based on molecule or custom ID

  // sets bodytag for owned atoms

  // body attributes are computed later by setup_bodies()



  create_bodies();

  create_bodies(bodyid);

  if (customflag) delete [] bodyid;



  // set nlocal_body and allocate bodies I own
@@ -1424,7 +1476,7 @@ void FixRigidSmall::set_v() 
   set bodytag for all owned atoms

------------------------------------------------------------------------- */



void FixRigidSmall::create_bodies()

void FixRigidSmall::create_bodies(tagint *bodyid)

{

  int i,m,n;

  double unwrap[3];
@@ -1464,8 +1516,8 @@ void FixRigidSmall::create_bodies() 
  double *buf;

  memory->create(buf,ncount*percount,"rigid/small:buf");



  // create map hash for storing unique molecule IDs of my atoms

  // key = molecule ID

  // create map hash for storing unique body IDs of my atoms

  // key = body ID

  // value = index into per-body data structure

  // n = # of entries in hash
@@ -1477,12 +1529,10 @@ void FixRigidSmall::create_bodies() 
  // value = index into N-length data structure

  // n = count of unique bodies my atoms are part of



  tagint *molecule = atom->molecule;



  n = 0;

  for (i = 0; i < nlocal; i++) {

    if (!(mask[i] & groupbit)) continue;

    if (hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = n++;

    if (hash->find(bodyid[i]) == hash->end()) (*hash)[bodyid[i]] = n++;

  }



  // bbox = bounding box of each rigid body my atoms are part of
@@ -1494,7 +1544,7 @@ void FixRigidSmall::create_bodies() 
    bbox[i][1] = bbox[i][3] = bbox[i][5] = -BIG;

  }



  // pack my atoms into buffer as molecule ID, unwrapped coords

  // pack my atoms into buffer as body ID, unwrapped coords



  double **x = atom->x;
@@ -1502,7 +1552,7 @@ void FixRigidSmall::create_bodies() 
  for (i = 0; i < nlocal; i++) {

    if (!(mask[i] & groupbit)) continue;

    domain->unmap(x[i],image[i],unwrap);

    buf[m++] = molecule[i];

    buf[m++] = bodyid[i];

    buf[m++] = unwrap[0];

    buf[m++] = unwrap[1];

    buf[m++] = unwrap[2];
@@ -1542,7 +1592,7 @@ void FixRigidSmall::create_bodies() 


  for (i = 0; i < n; i++) rsqclose[i] = BIG;



  // pack my atoms into buffer as molecule ID, atom ID, unwrapped coords

  // pack my atoms into buffer as body ID, atom ID, unwrapped coords



  tagint *tag = atom->tag;
@@ -1550,7 +1600,7 @@ void FixRigidSmall::create_bodies() 
  for (i = 0; i < nlocal; i++) {

    if (!(mask[i] & groupbit)) continue;

    domain->unmap(x[i],image[i],unwrap);

    buf[m++] = molecule[i];

    buf[m++] = bodyid[i];

    buf[m++] = ubuf(tag[i]).d;

    buf[m++] = unwrap[0];

    buf[m++] = unwrap[1];
@@ -1570,7 +1620,7 @@ void FixRigidSmall::create_bodies() 
  for (i = 0; i < nlocal; i++) {

    bodytag[i] = 0;

    if (!(mask[i] & groupbit)) continue;

    m = hash->find(molecule[i])->second;

    m = hash->find(bodyid[i])->second;

    bodytag[i] = idclose[m];

    rsqmax = MAX(rsqmax,rsqclose[m]);

  }

src/RIGID/fix_rigid_small.h

+2 −1
Original line number Diff line number Diff line
@@ -79,6 +79,7 @@ class FixRigidSmall : public Fix { 
  char *infile;             // file to read rigid body attributes from

  int setupflag;            // 1 if body properties are setup, else 0

  int commflag;             // various modes of forward/reverse comm

  int customflag;           // 1 if custom property/variable define bodies

  int nbody;                // total # of rigid bodies

  int nlinear;              // total # of linear rigid bodies

  tagint maxmol;            // max mol-ID
@@ -187,7 +188,7 @@ class FixRigidSmall : public Fix { 
  void image_shift();

  void set_xv();

  void set_v();

  void create_bodies();

  void create_bodies(tagint *);

  void setup_bodies_static();

  void setup_bodies_dynamic();

  void readfile(int, double **, int *);

CRA Git | Maintained and supported by SUSTech CRA and CCSE