Merge pull request #1670 from akohlmey/fix-rigid-nh-no-extended

Consolidate constants and enumerators across rigid fix styles
This commit is contained in:
Axel Kohlmeyer
2019-09-13 08:04:16 -04:00
committed by GitHub
10 changed files with 82 additions and 111 deletions

View File

@ -34,25 +34,12 @@
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
using namespace RigidConst;
enum{SINGLE,MOLECULE,GROUP};
enum{NONE,XYZ,XY,YZ,XZ};
enum{ISO,ANISO,TRICLINIC};
#define MAXLINE 1024
#define CHUNK 1024
#define ATTRIBUTE_PERBODY 20
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -605,21 +592,6 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
for (ibody = 0; ibody < nbody; ibody++) for (ibody = 0; ibody < nbody; ibody++)
if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body"); if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body");
// bitmasks for properties of extended particles
POINT = 1;
SPHERE = 2;
ELLIPSOID = 4;
LINE = 8;
TRIANGLE = 16;
DIPOLE = 32;
OMEGA = 64;
ANGMOM = 128;
TORQUE = 256;
MINUSPI = -MY_PI;
TWOPI = 2.0*MY_PI;
// wait to setup bodies until first init() using current atom properties // wait to setup bodies until first init() using current atom properties
setupflag = 0; setupflag = 0;
@ -1472,8 +1444,8 @@ void FixRigid::set_xv()
if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
else theta_body = -2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]);
theta = orient[i][0] + theta_body; theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI; while (theta <= -MY_PI) theta += MY_2PI;
while (theta > MY_PI) theta -= TWOPI; while (theta > MY_PI) theta -= MY_2PI;
lbonus[line[i]].theta = theta; lbonus[line[i]].theta = theta;
omega_one[i][0] = omega[ibody][0]; omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1]; omega_one[i][1] = omega[ibody][1];
@ -2018,8 +1990,8 @@ void FixRigid::setup_bodies_static()
if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
else theta_body = -2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]);
orient[i][0] = lbonus[line[i]].theta - theta_body; orient[i][0] = lbonus[line[i]].theta - theta_body;
while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI; while (orient[i][0] <= -MY_PI) orient[i][0] += MY_2PI;
while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI; while (orient[i][0] > MY_PI) orient[i][0] -= MY_2PI;
if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0; if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (eflags[i] & TRIANGLE) { } else if (eflags[i] & TRIANGLE) {
quatatom = tbonus[tri[i]].quat; quatatom = tbonus[tri[i]].quat;

View File

@ -65,7 +65,6 @@ class FixRigid : public Fix {
double dtv,dtf,dtq; double dtv,dtf,dtq;
double *step_respa; double *step_respa;
int triclinic; int triclinic;
double MINUSPI,TWOPI;
char *inpfile; // file to read rigid body attributes from char *inpfile; // file to read rigid body attributes from
int rstyle; // SINGLE,MOLECULE,GROUP int rstyle; // SINGLE,MOLECULE,GROUP
@ -137,9 +136,6 @@ class FixRigid : public Fix {
class AtomVecLine *avec_line; class AtomVecLine *avec_line;
class AtomVecTri *avec_tri; class AtomVecTri *avec_tri;
int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE;
void image_shift(); void image_shift();
void set_xv(); void set_xv();
void set_v(); void set_v();

View File

@ -33,14 +33,11 @@
#include "kspace.h" #include "kspace.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace RigidConst;
enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
#define EPSILON 1.0e-7
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

View File

@ -34,17 +34,12 @@
#include "kspace.h" #include "kspace.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathExtra; using namespace MathExtra;
using namespace RigidConst;
enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
#define EPSILON 1.0e-7
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

View File

@ -38,34 +38,17 @@
#include "hashlittle.h" #include "hashlittle.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "rigid_const.h"
#include <map> #include <map>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
using namespace RigidConst;
#define RVOUS 1 // 0 for irregular, 1 for all2all #define RVOUS 1 // 0 for irregular, 1 for all2all
#define MAXLINE 1024
#define CHUNK 1024
#define ATTRIBUTE_PERBODY 20
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define BIG 1.0e20
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
#define DELTA_BODY 10000
enum{NONE,XYZ,XY,YZ,XZ}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
@ -455,21 +438,6 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 1 + bodysize; comm_forward = 1 + bodysize;
comm_reverse = 6; comm_reverse = 6;
// bitmasks for properties of extended particles
POINT = 1;
SPHERE = 2;
ELLIPSOID = 4;
LINE = 8;
TRIANGLE = 16;
DIPOLE = 32;
OMEGA = 64;
ANGMOM = 128;
TORQUE = 256;
MINUSPI = -MY_PI;
TWOPI = 2.0*MY_PI;
// atom style pointers to particles that store extra info // atom style pointers to particles that store extra info
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
@ -1384,8 +1352,8 @@ void FixRigidSmall::set_xv()
if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]); if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
else theta_body = -2.0*acos(b->quat[0]); else theta_body = -2.0*acos(b->quat[0]);
theta = orient[i][0] + theta_body; theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI; while (theta <= -MY_PI) theta += MY_2PI;
while (theta > MY_PI) theta -= TWOPI; while (theta > MY_PI) theta -= MY_2PI;
lbonus[line[i]].theta = theta; lbonus[line[i]].theta = theta;
omega[i][0] = b->omega[0]; omega[i][0] = b->omega[0];
omega[i][1] = b->omega[1]; omega[i][1] = b->omega[1];
@ -2155,8 +2123,8 @@ void FixRigidSmall::setup_bodies_static()
if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]); if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
else theta_body = -2.0*acos(b->quat[0]); else theta_body = -2.0*acos(b->quat[0]);
orient[i][0] = lbonus[line[i]].theta - theta_body; orient[i][0] = lbonus[line[i]].theta - theta_body;
while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI; while (orient[i][0] <= -MY_PI) orient[i][0] += MY_2PI;
while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI; while (orient[i][0] > MY_PI) orient[i][0] -= MY_2PI;
if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0; if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (eflags[i] & TRIANGLE) { } else if (eflags[i] & TRIANGLE) {
quatatom = tbonus[tri[i]].quat; quatatom = tbonus[tri[i]].quat;

View File

@ -72,7 +72,6 @@ class FixRigidSmall : public Fix {
double dtv,dtf,dtq; double dtv,dtf,dtq;
double *step_respa; double *step_respa;
int triclinic; int triclinic;
double MINUSPI,TWOPI;
char *inpfile; // file to read rigid body attributes from char *inpfile; // file to read rigid body attributes from
int setupflag; // 1 if body properties are setup, else 0 int setupflag; // 1 if body properties are setup, else 0
@ -129,9 +128,6 @@ class FixRigidSmall : public Fix {
int dorientflag; // 1 if particles store dipole orientation int dorientflag; // 1 if particles store dipole orientation
int reinitflag; // 1 if re-initialize rigid bodies between runs int reinitflag; // 1 if re-initialize rigid bodies between runs
int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE;
class AtomVecEllipsoid *avec_ellipsoid; class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line; class AtomVecLine *avec_line;
class AtomVecTri *avec_tri; class AtomVecTri *avec_tri;

54
src/RIGID/rigid_const.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_RIGID_CONST_H
#define LMP_RIGID_CONST_H
namespace LAMMPS_NS {
namespace RigidConst {
enum{SINGLE,MOLECULE,GROUP};
enum{NONE,XYZ,XY,YZ,XZ};
enum{ISO,ANISO,TRICLINIC};
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
enum {POINT = 1<<0,
SPHERE = 1<<1,
ELLIPSOID = 1<<2,
LINE = 1<<3,
TRIANGLE = 1<<4,
DIPOLE = 1<<5,
OMEGA = 1<<6,
ANGMOM = 1<<7,
TORQUE = 1<<8
};
static const double TOLERANCE = 1.0e-6;
static const double EPSILON = 1.0e-7;
static const double BIG = 1.0e20;
// moment of inertia prefactor for sphere
static const double SINERTIA = 0.4;
// moment of inertia prefactor for ellipsoid
static const double EINERTIA = 0.2;
// moment of inertia prefactor for line segment
static const double LINERTIA = 1.0/12.0;
static const int MAXLINE = 1024;
static const int CHUNK = 1024;
static const int DELTA_BODY = 10000;
static const int ATTRIBUTE_PERBODY = 20;
}
}
#endif

View File

@ -38,15 +38,12 @@
#include "math_extra.h" #include "math_extra.h"
#include "math_const.h" #include "math_const.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
using namespace RigidConst;
enum{SINGLE,MOLECULE,GROUP}; // same as in FixRigid
enum{ISO,ANISO,TRICLINIC}; // same as in FixRigid
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid
typedef struct { double x,y,z; } dbl3_t; typedef struct { double x,y,z; } dbl3_t;
@ -770,8 +767,8 @@ void FixRigidNHOMP::set_xv_thr()
if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
else theta_body = -2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]);
theta = orient[i][0] + theta_body; theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI; while (theta <= -MY_PI) theta += MY_2PI;
while (theta > MY_PI) theta -= TWOPI; while (theta > MY_PI) theta -= MY_2PI;
lbonus[line[i]].theta = theta; lbonus[line[i]].theta = theta;
omega_one[i][0] = omega[ibody][0]; omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1]; omega_one[i][1] = omega[ibody][1];

View File

@ -33,14 +33,12 @@
#include "math_extra.h" #include "math_extra.h"
#include "math_const.h" #include "math_const.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
using namespace RigidConst;
enum{SINGLE,MOLECULE,GROUP}; // same as in FixRigid
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
typedef struct { double x,y,z; } dbl3_t; typedef struct { double x,y,z; } dbl3_t;
@ -488,8 +486,8 @@ void FixRigidOMP::set_xv_thr()
if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]); if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
else theta_body = -2.0*acos(quat[ibody][0]); else theta_body = -2.0*acos(quat[ibody][0]);
theta = orient[i][0] + theta_body; theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI; while (theta <= -MY_PI) theta += MY_2PI;
while (theta > MY_PI) theta -= TWOPI; while (theta > MY_PI) theta -= MY_2PI;
lbonus[line[i]].theta = theta; lbonus[line[i]].theta = theta;
omega_one[i][0] = omega[ibody][0]; omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1]; omega_one[i][1] = omega[ibody][1];

View File

@ -31,14 +31,12 @@
#include "math_extra.h" #include "math_extra.h"
#include "math_const.h" #include "math_const.h"
#include "rigid_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
using namespace RigidConst;
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};
typedef struct { double x,y,z; } dbl3_t; typedef struct { double x,y,z; } dbl3_t;
@ -427,8 +425,8 @@ void FixRigidSmallOMP::set_xv_thr()
if (b.quat[3] >= 0.0) theta_body = 2.0*acos(b.quat[0]); if (b.quat[3] >= 0.0) theta_body = 2.0*acos(b.quat[0]);
else theta_body = -2.0*acos(b.quat[0]); else theta_body = -2.0*acos(b.quat[0]);
theta = orient[i][0] + theta_body; theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI; while (theta <= -MY_PI) theta += MY_2PI;
while (theta > MY_PI) theta -= TWOPI; while (theta > MY_PI) theta -= MY_2PI;
lbonus[line[i]].theta = theta; lbonus[line[i]].theta = theta;
omega[i][0] = b.omega[0]; omega[i][0] = b.omega[0];
omega[i][1] = b.omega[1]; omega[i][1] = b.omega[1];