ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7146 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2011-10-20 14:56:49 +00:00
parent 83b55b844c
commit ca8dfd3d7c
11 changed files with 1688 additions and 568 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1164
src/SRD/fix_srd.cpp
View File

File diff suppressed because it is too large Load Diff

52
src/SRD/fix_srd.h
View File

@ -43,9 +43,9 @@ class FixSRD : public Fix {
int me,nprocs;
int bigexist,biggroup,biggroupbit;
int collidestyle,lamdaflag,overlap,insideflag,exactflag,maxbounceallow;
int cubicflag,shiftuser,shiftseed,shiftflag,streamflag;
int cubicflag,shiftuser,shiftseed,shiftflag,tstat;
double gridsrd,gridsearch,lamda,radfactor,cubictol;
int triclinic,change_size,change_shape;
int triclinic,change_size,change_shape,deformflag;
double dt_big,dt_srd;
double mass_big,mass_srd;
@ -64,6 +64,8 @@ class FixSRD : public Fix {
double walltrigger;
class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line;
class AtomVecTri *avec_tri;
// for orthogonal box, these are in box units
// for triclinic box, these are in lamda units
@ -92,18 +94,21 @@ class FixSRD : public Fix {
struct Big {
int index; // local index of particle/wall
int type; // SPHERE or ELLIPSOID or WALL
int type; // SPHERE or ELLIPSOID or LINE or TRI or WALL
double radius,radsq; // radius of sphere
double aradsqinv; // 3 ellipsoid radii
double bradsqinv;
double cradsqinv;
double length; // length of line segment
double normbody[3]; // normal of tri in body-frame
double cutbinsq; // add big to bin if within this distance
double omega[3]; // current omega for sphere or ellipsoid
double ex[3],ey[3],ez[3]; // current orientation vecs for ellipsoid
double omega[3]; // current omega for sphere/ellipsoid/tri/line
double ex[3],ey[3],ez[3]; // current orientation vecs for ellipsoid/tri
double norm[3]; // current unit normal of tri in space-frame
double theta; // current orientation of line
};
Big *biglist; // list of info for each owned & ghost big and wall
int any_ellipsoids; // 1 if any big particles are ellipsoids
int torqueflag; // 1 if any big particle is torqued
// current size of particle-based arrays
@ -123,7 +128,7 @@ class FixSRD : public Fix {
int owner; // 1 if I am owner of this bin, 0 if not
int n; // # of SRD particles in bin
double xctr[3]; // center point of bin, only used for triclinic
double vave[3]; // sum of v components for SRD particles in bin
double vsum[3]; // sum of v components for SRD particles in bin
double random; // random value if I am owner
};
@ -167,6 +172,17 @@ class FixSRD : public Fix {
int maxstencil; // max # of bins stencil array can hold
int **stencil; // list of 3d bin offsets a big particle can overlap
// persistent data for line/tri collision calculations
double tfraction,theta0,theta1;
double xs0[3],xs1[3],xsc[3];
double xb0[3],xb1[3],xbc[3];
double nbc[3];
// shared data for triangle collision calculations
// private functions
void reset_velocities();
void vbin_comm(int);
void vbin_pack(BinAve *, int, int *, double *);
@ -177,6 +193,8 @@ class FixSRD : public Fix {
int inside_sphere(double *, double *, Big *);
int inside_ellipsoid(double *, double *, Big *);
int inside_line(double *, double *, double *, double *, Big *, double);
int inside_tri(double *, double *, double *, double *, Big *, double);
int inside_wall(double *, int);
double collision_sphere_exact(double *, double *, double *, double *,
@ -187,18 +205,19 @@ class FixSRD : public Fix {
Big *, double *, double *, double *);
void collision_ellipsoid_inexact(double *, double *,
Big *, double *, double *, double *);
double collision_line_exact(double *, double *, double *, double *,
Big *, double, double *, double *, double *);
double collision_tri_exact(double *, double *, double *, double *,
Big *, double, double *, double *, double *);
double collision_wall_exact(double *, int, double *,
double *, double *, double *);
void collision_wall_inexact(double *, int, double *, double *, double *);
void slip_sphere(double *, double *, double *, double *);
void slip_ellipsoid(double *, double *, double *, Big *,
double *, double *, double *);
void slip(double *, double *, double *, Big *,
double *, double *, double *);
void slip_wall(double *, int, double *, double *);
void noslip(double *, double *, double *, Big *,
void noslip(double *, double *, double *, Big *, int,
double *, double *, double *);
void noslip_wall(double *, int, double *, double *, double *);
void force_torque(double *, double *, double *,
double *, double *, double *);
@ -217,11 +236,12 @@ class FixSRD : public Fix {
double point_bin_distance(double *, int, int, int);
double bin_bin_distance(int, int, int);
void exyz_from_q(double *, double *, double *, double *);
void omega_from_mq(double *, double *, double *, double *,
double, double *, double *);
void velocity_stats(int);
double newton_raphson(double, double);
void lineside(double, double &, double &);
void triside(double, double &, double &);
double distance(int, int);
void print_collision(int, int, int, double, double,
double *, double *, double *, int);

12
src/atom.cpp
View File

@ -75,7 +75,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
radius = rmass = NULL;
vfrac = s0 = NULL;
x0 = NULL;
ellipsoid = NULL;
ellipsoid = line = tri = NULL;
spin = NULL;
eradius = ervel = erforce = NULL;
cs = csforce = vforce = ervelforce = NULL;
@ -106,7 +106,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// initialize atom style and array existence flags
// customize by adding new flag
sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
sphere_flag = ellipsoid_flag = line_flag = tri_flag = 0;
peri_flag = electron_flag = 0;
wavepacket_flag = sph_flag = 0;
molecule_flag = q_flag = mu_flag = 0;
@ -185,6 +186,8 @@ Atom::~Atom()
memory->destroy(s0);
memory->destroy(x0);
memory->destroy(ellipsoid);
memory->destroy(line);
memory->destroy(tri);
memory->destroy(spin);
memory->destroy(eradius);
memory->destroy(ervel);
@ -259,8 +262,9 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
// may have been set by old avec
// customize by adding new flag
sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
sphere_flag = ellipsoid_flag = line_flag = tri_flag = 0;
peri_flag = electron_flag = 0;
molecule_flag = q_flag = mu_flag = 0;
rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;

4
src/atom.h
View File

@ -51,7 +51,7 @@ class Atom : protected Pointers {
double **omega,**angmom,**torque;
double *radius,*rmass,*vfrac,*s0;
double **x0;
int *ellipsoid;
int *ellipsoid,*line,*tri;
int *spin;
double *eradius,*ervel,*erforce,*ervelforce;
double *cs,*csforce,*vforce;
@ -84,7 +84,7 @@ class Atom : protected Pointers {
// atom style and per-atom array existence flags
// customize by adding new flag
int sphere_flag,ellipsoid_flag,peri_flag,electron_flag;
int sphere_flag,ellipsoid_flag,line_flag,tri_flag,peri_flag,electron_flag;
int wavepacket_flag,sph_flag;
int molecule_flag,q_flag,mu_flag;

428
src/compute_property_atom.cpp
View File

@ -11,10 +11,14 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "compute_property_atom.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
@ -173,39 +177,39 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
pack_choice[i] = &ComputePropertyAtom::pack_angmomz;
} else if (strcmp(arg[iarg],"shapex") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapex;
} else if (strcmp(arg[iarg],"shapey") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapey;
} else if (strcmp(arg[iarg],"shapez") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapez;
} else if (strcmp(arg[iarg],"quatw") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatw;
} else if (strcmp(arg[iarg],"quati") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quati;
} else if (strcmp(arg[iarg],"quatj") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatj;
} else if (strcmp(arg[iarg],"quatk") == 0) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatk;
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
@ -244,6 +248,83 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_erforce;
} else if (strcmp(arg[iarg],"end1x") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end1x;
} else if (strcmp(arg[iarg],"end1y") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end1y;
} else if (strcmp(arg[iarg],"end1z") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end1z;
} else if (strcmp(arg[iarg],"end2x") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end2x;
} else if (strcmp(arg[iarg],"end2y") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end2y;
} else if (strcmp(arg[iarg],"end2z") == 0) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_end2z;
} else if (strcmp(arg[iarg],"corner1x") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner1x;
} else if (strcmp(arg[iarg],"corner1y") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner1y;
} else if (strcmp(arg[iarg],"corner1z") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner1z;
} else if (strcmp(arg[iarg],"corner2x") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner2x;
} else if (strcmp(arg[iarg],"corner2y") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner2y;
} else if (strcmp(arg[iarg],"corner2z") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner2z;
} else if (strcmp(arg[iarg],"corner3x") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner3x;
} else if (strcmp(arg[iarg],"corner3y") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner3y;
} else if (strcmp(arg[iarg],"corner3z") == 0) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri) error->all(FLERR,"Compute property/atom for "
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_corner3z;
} else error->all(FLERR,"Invalid keyword in compute property/atom command");
}
@ -994,7 +1075,7 @@ void ComputePropertyAtom::pack_angmomz(int n)
void ComputePropertyAtom::pack_shapex(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1011,7 +1092,7 @@ void ComputePropertyAtom::pack_shapex(int n)
void ComputePropertyAtom::pack_shapey(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1028,7 +1109,7 @@ void ComputePropertyAtom::pack_shapey(int n)
void ComputePropertyAtom::pack_shapez(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1045,7 +1126,7 @@ void ComputePropertyAtom::pack_shapez(int n)
void ComputePropertyAtom::pack_quatw(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1062,7 +1143,7 @@ void ComputePropertyAtom::pack_quatw(int n)
void ComputePropertyAtom::pack_quati(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1079,7 +1160,7 @@ void ComputePropertyAtom::pack_quati(int n)
void ComputePropertyAtom::pack_quatj(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1096,7 +1177,7 @@ void ComputePropertyAtom::pack_quatj(int n)
void ComputePropertyAtom::pack_quatk(int n)
{
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -1213,3 +1294,294 @@ void ComputePropertyAtom::pack_erforce(int n)
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end1x(int n)
{
AtomVecLine::Bonus *bonus = avec_line->bonus;
int *line = atom->line;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0)
buf[n] = x[i][0] - 0.5*bonus[line[i]].length*cos(bonus[line[i]].theta);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end1y(int n)
{
AtomVecLine::Bonus *bonus = avec_line->bonus;
int *line = atom->line;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0)
buf[n] = x[i][1] - 0.5*bonus[line[i]].length*sin(bonus[line[i]].theta);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end1z(int n)
{
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = x[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end2x(int n)
{
AtomVecLine::Bonus *bonus = avec_line->bonus;
int *line = atom->line;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0)
buf[n] = x[i][0] + 0.5*bonus[line[i]].length*cos(bonus[line[i]].theta);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end2y(int n)
{
AtomVecLine::Bonus *bonus = avec_line->bonus;
int *line = atom->line;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0)
buf[n] = x[i][1] + 0.5*bonus[line[i]].length*sin(bonus[line[i]].theta);
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_end2z(int n)
{
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = x[i][2];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner1x(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c1,c);
buf[n] = x[i][0] + c[0];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner1y(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c1,c);
buf[n] = x[i][1] + c[1];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner1z(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c1,c);
buf[n] = x[i][2] + c[2];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner2x(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c2,c);
buf[n] = x[i][0] + c[0];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner2y(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c2,c);
buf[n] = x[i][1] + c[1];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner2z(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c2,c);
buf[n] = x[i][2] + c[2];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner3x(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c3,c);
buf[n] = x[i][0] + c[0];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner3y(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c3,c);
buf[n] = x[i][1] + c[1];
} else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyAtom::pack_corner3z(int n)
{
AtomVecTri::Bonus *bonus = avec_tri->bonus;
int *tri = atom->tri;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double p[3][3],c[3];
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && tri[i] >= 0) {
MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
MathExtra::matvec(p,bonus[tri[i]].c3,c);
buf[n] = x[i][2] + c[2];
} else buf[n] = 0.0;
n += nvalues;
}
}

19
src/compute_property_atom.h
View File

@ -38,7 +38,9 @@ class ComputePropertyAtom : public Compute {
double *vector;
double **array;
double *buf;
class AtomVecEllipsoid *avec;
class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line;
class AtomVecTri *avec_tri;
typedef void (ComputePropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
@ -101,6 +103,21 @@ class ComputePropertyAtom : public Compute {
void pack_eradius(int);
void pack_ervel(int);
void pack_erforce(int);
void pack_end1x(int);
void pack_end1y(int);
void pack_end1z(int);
void pack_end2x(int);
void pack_end2y(int);
void pack_end2z(int);
void pack_corner1x(int);
void pack_corner1y(int);
void pack_corner1z(int);
void pack_corner2x(int);
void pack_corner2y(int);
void pack_corner2z(int);
void pack_corner3x(int);
void pack_corner3y(int);
void pack_corner3z(int);
};
}

325
src/fix_rigid.cpp
View File

@ -19,6 +19,8 @@
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
@ -28,14 +30,20 @@
#include "random_mars.h"
#include "force.h"
#include "output.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
/* ---------------------------------------------------------------------- */
FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
@ -56,12 +64,12 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays
// register with Atom class
extended = dorientflag = qorientflag = 0;
extended = orientflag = dorientflag = 0;
body = NULL;
displace = NULL;
eflags = NULL;
orient = NULL;
dorient = NULL;
qorient = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
@ -278,7 +286,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else error->all(FLERR,"Illegal fix rigid command");
if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0))
error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation");
error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation");
int count = 0;
for (int m = mlo; m <= mhi; m++) {
@ -384,18 +392,24 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// bitmasks for properties of extended particles
INERTIA_POINT = 1;
INERTIA_SPHERE = 2;
INERTIA_ELLIPSOID = 4;
ORIENT_DIPOLE = 8;
ORIENT_QUAT = 16;
OMEGA = 32;
ANGMOM = 64;
TORQUE = 128;
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
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec_line = (AtomVecLine *) atom->style_match("line");
avec_tri = (AtomVecTri *) atom->style_match("tri");
// print statistics
@ -423,8 +437,8 @@ FixRigid::~FixRigid()
memory->destroy(body);
memory->destroy(displace);
memory->destroy(eflags);
memory->destroy(orient);
memory->destroy(dorient);
memory->destroy(qorient);
// delete nbody-length arrays
@ -469,7 +483,7 @@ int FixRigid::setmask()
void FixRigid::init()
{
int i,ibody;
int i,itype,ibody;
triclinic = domain->triclinic;
@ -504,24 +518,33 @@ void FixRigid::init()
// extended = 1 if any particle in a rigid body is finite size
// or has a dipole moment
extended = dorientflag = qorientflag = 0;
extended = orientflag = dorientflag = 0;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **mu = atom->mu;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *type = atom->type;
int nlocal = atom->nlocal;
if (atom->radius_flag || atom->ellipsoid_flag || atom->mu_flag) {
if (atom->radius_flag || atom->ellipsoid_flag || atom->line_flag ||
atom->tri_flag || atom->mu_flag) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
if (line && line[i] >= 0) flag = 1;
if (tri && tri[i] >= 0) flag = 1;
if (mu && mu[i][3] > 0.0) flag = 1;
}
@ -529,35 +552,45 @@ void FixRigid::init()
}
// grow extended arrays and set extended flags for each particle
// qorientflag = 1 if any particle stores quat orientation
// orientflag = 4 if any particle stores ellipsoid or tri orientation
// orientflag = 1 if any particle stores line orientation
// dorientflag = 1 if any particle stores dipole orientation
if (extended) {
if (atom->ellipsoid_flag) orientflag = 4;
if (atom->line_flag) orientflag = 1;
if (atom->tri_flag) orientflag = 4;
if (atom->mu_flag) dorientflag = 1;
if (atom->ellipsoid_flag) qorientflag = 1;
grow_arrays(atom->nmax);
for (i = 0; i < nlocal; i++) {
eflags[i] = 0;
if (body[i] < 0) continue;
// set INERTIA to POINT or SPHERE or ELLIPSOID
// set to POINT or SPHERE or ELLIPSOID or LINE
if (radius && radius[i] > 0.0) {
eflags[i] |= INERTIA_SPHERE;
eflags[i] |= SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
} else if (ellipsoid && ellipsoid[i] >= 0) {
eflags[i] |= INERTIA_ELLIPSOID;
eflags[i] |= ORIENT_QUAT;
eflags[i] |= ELLIPSOID;
eflags[i] |= ANGMOM;
eflags[i] |= TORQUE;
} else eflags[i] |= INERTIA_POINT;
} else if (line && line[i] >= 0) {
eflags[i] |= LINE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
} else if (tri && tri[i] >= 0) {
eflags[i] |= TRIANGLE;
eflags[i] |= ANGMOM;
eflags[i] |= TORQUE;
} else eflags[i] |= POINT;
// set DIPOLE if atom->mu and mu[3] > 0.0
if (atom->mu_flag && mu[i][3] > 0.0)
eflags[i] |= ORIENT_DIPOLE;
eflags[i] |= DIPOLE;
}
}
@ -592,8 +625,8 @@ void FixRigid::init()
if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) ||
(zbox && !periodicity[2]))
error->one(FLERR,"Fix rigid atom has non-zero image flag "
"in a non-periodic dimension");
error->one(FLERR,"Fix rigid atom has non-zero image flag "
"in a non-periodic dimension");
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
@ -629,7 +662,7 @@ void FixRigid::init()
// dx,dy,dz = coords relative to center-of-mass
// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
double dx,dy,dz;
double dx,dy,dz,rad;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@ -671,18 +704,19 @@ void FixRigid::init()
if (extended) {
double ivec[6];
double *shape,*quatatom;
double *shape,*quatatom,*inertiaatom;
double length,theta;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
massone = rmass[i];
if (eflags[i] & INERTIA_SPHERE) {
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
} else if (eflags[i] & INERTIA_ELLIPSOID) {
if (eflags[i] & SPHERE) {
sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i];
sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i];
sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
@ -692,6 +726,26 @@ void FixRigid::init()
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
} else if (eflags[i] & LINE) {
length = lbonus[line[i]].length;
theta = lbonus[line[i]].theta;
MathExtra::inertia_line(length,theta,massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::inertia_triangle(inertiaatom,quatatom,massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
}
}
}
@ -715,7 +769,8 @@ void FixRigid::init()
tensor[0][1] = tensor[1][0] = all[ibody][5];
ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
if (ierror) error->all(FLERR,"Insufficient Jacobi rotations for rigid body");
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for rigid body");
ex_space[ibody][0] = evectors[0][0];
ex_space[ibody][1] = evectors[1][0];
@ -756,6 +811,7 @@ void FixRigid::init()
double qc[4],delta[3];
double *quatatom;
double theta_body;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) {
@ -786,20 +842,34 @@ void FixRigid::init()
ez_space[ibody],delta,displace[i]);
if (extended) {
if (eflags[i] & ORIENT_DIPOLE) {
if (eflags[i] & ELLIPSOID) {
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::qconjugate(quat[ibody],qc);
MathExtra::quatquat(qc,quatatom,orient[i]);
MathExtra::qnormalize(orient[i]);
} else if (eflags[i] & LINE) {
if (quat[ibody][3] >= 0.0) 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;
while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI;
while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI;
if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (eflags[i] & TRIANGLE) {
quatatom = tbonus[tri[i]].quat;
MathExtra::qconjugate(quat[ibody],qc);
MathExtra::quatquat(qc,quatatom,orient[i]);
MathExtra::qnormalize(orient[i]);
} else if (orientflag == 4) {
orient[i][0] = orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
} else if (orientflag == 1)
orient[i][0] = 0.0;
if (eflags[i] & DIPOLE) {
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],
ez_space[ibody],mu[i],dorient[i]);
MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]);
} else if (dorientflag)
dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
if (eflags[i] & ORIENT_QUAT) {
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::qconjugate(quat[ibody],qc);
MathExtra::quatquat(qc,quatatom,qorient[i]);
MathExtra::qnormalize(qorient[i]);
} else if (qorientflag)
qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
}
}
@ -831,20 +901,39 @@ void FixRigid::init()
if (extended) {
double ivec[6];
double *shape;
double *shape,*inertiaatom;
double length;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
massone = rmass[i];
if (eflags[i] & INERTIA_SPHERE) {
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
} else if (eflags[i] & INERTIA_ELLIPSOID) {
if (eflags[i] & SPHERE) {
sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i];
sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i];
sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i];
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
} else if (eflags[i] & LINE) {
length = lbonus[line[i]].length;
MathExtra::inertia_line(length,orient[i][0],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
@ -894,7 +983,8 @@ void FixRigid::init()
ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
}
tfactor = force->mvv2e / (ndof * force->boltz);
if (ndof > 0.0) tfactor = force->mvv2e / (ndof * force->boltz);
else tfactor = 0.0;
}
/* ---------------------------------------------------------------------- */
@ -996,24 +1086,29 @@ void FixRigid::setup(int vflag)
// extended particles add their rotation/torque to angmom/torque of body
if (extended) {
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **torque_one = atom->torque;
double *radius = atom->radius;
int *line = atom->line;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & OMEGA) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
radone = radius[i];
sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0];
sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1];
sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2];
if (eflags[i] & SPHERE) {
radone = radius[i];
sum[ibody][0] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0];
sum[ibody][1] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1];
sum[ibody][2] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2];
} else if (eflags[i] & LINE) {
radone = lbonus[line[i]].length;
sum[ibody][2] += LINERTIA*rmass[i] * radone*radone * omega_one[i][2];
}
}
if (eflags[i] & ANGMOM) {
sum[ibody][0] += angmom_one[i][0];
sum[ibody][1] += angmom_one[i][1];
@ -1516,7 +1611,7 @@ void FixRigid::deform(int flag)
void FixRigid::set_xv()
{
int ibody;
int ibody,itype;
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
@ -1622,43 +1717,64 @@ void FixRigid::set_xv()
// set orientation, omega, angmom of each extended particle
if (extended) {
double *shape,*quatatom;
double theta_body,theta;
double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecLine::Bonus *lbonus;
if (avec_line) lbonus = avec_line->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **mu = atom->mu;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & ORIENT_DIPOLE) {
MathExtra::quat_to_mat(quat[ibody],p);
MathExtra::matvec(p,dorient[i],mu[i]);
MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}
if (eflags[i] & ORIENT_QUAT) {
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::quatquat(quat[ibody],qorient[i],quatatom);
MathExtra::qnormalize(quatatom);
}
if (eflags[i] & OMEGA) {
if (eflags[i] & SPHERE) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::quatquat(quat[ibody],orient[i],quatatom);
MathExtra::qnormalize(quatatom);
ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
} else if (eflags[i] & LINE) {
if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
else theta_body = -2.0*acos(quat[ibody][0]);
theta = orient[i][0] + theta_body;
while (theta <= MINUSPI) theta += TWOPI;
while (theta > MY_PI) theta -= TWOPI;
lbonus[line[i]].theta = theta;
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::quatquat(quat[ibody],orient[i],quatatom);
MathExtra::qnormalize(quatatom);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,
inertiaatom,angmom_one[i]);
}
if (eflags[i] & DIPOLE) {
MathExtra::quat_to_mat(quat[ibody],p);
MathExtra::matvec(p,dorient[i],mu[i]);
MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}
}
}
@ -1672,8 +1788,9 @@ void FixRigid::set_xv()
void FixRigid::set_v()
{
int ibody;
int ibody,itype;
int xbox,ybox,zbox;
double dx,dy,dz;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6];
@ -1761,32 +1878,44 @@ void FixRigid::set_v()
// set omega, angmom of each extended particle
if (extended) {
double *shape,*quatatom;
double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
AtomVecTri::Bonus *tbonus;
if (avec_tri) tbonus = avec_tri->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
int *ellipsoid = atom->ellipsoid;
int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & OMEGA) {
if (eflags[i] & SPHERE) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
}
if (eflags[i] & ANGMOM) {
} else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
} else if (eflags[i] & LINE) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
} else if (eflags[i] & TRIANGLE) {
inertiaatom = tbonus[tri[i]].inertia;
quatatom = tbonus[tri[i]].quat;
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,
inertiaatom,angmom_one[i]);
}
}
}
@ -1804,8 +1933,8 @@ double FixRigid::memory_usage()
bytes += maxvatom*6 * sizeof(double);
if (extended) {
bytes += nmax * sizeof(int);
if (orientflag) bytes = nmax*orientflag * sizeof(double);
if (dorientflag) bytes = nmax*3 * sizeof(double);
if (qorientflag) bytes = nmax*4 * sizeof(double);
}
return bytes;
}
@ -1820,8 +1949,8 @@ void FixRigid::grow_arrays(int nmax)
memory->grow(displace,nmax,3,"rigid:displace");
if (extended) {
memory->grow(eflags,nmax,"rigid:eflags");
if (orientflag) memory->grow(orient,nmax,orientflag,"rigid:orient");
if (dorientflag) memory->grow(dorient,nmax,3,"rigid:dorient");
if (qorientflag) memory->grow(qorient,nmax,4,"rigid:qorient");
}
}
@ -1837,17 +1966,13 @@ void FixRigid::copy_arrays(int i, int j)
displace[j][2] = displace[i][2];
if (extended) {
eflags[j] = eflags[i];
for (int k = 0; k < orientflag; k++)
orient[j][k] = orient[i][k];
if (dorientflag) {
dorient[j][0] = dorient[i][0];
dorient[j][1] = dorient[i][1];
dorient[j][2] = dorient[i][2];
}
if (qorientflag) {
qorient[j][0] = qorient[i][0];
qorient[j][1] = qorient[i][1];
qorient[j][2] = qorient[i][2];
qorient[j][3] = qorient[i][3];
}
}
}
@ -1877,17 +2002,13 @@ int FixRigid::pack_exchange(int i, double *buf)
int m = 4;
buf[m++] = eflags[i];
for (int j = 0; j < orientflag; j++)
buf[m++] = orient[i][j];
if (dorientflag) {
buf[m++] = dorient[i][0];
buf[m++] = dorient[i][1];
buf[m++] = dorient[i][2];
}
if (qorientflag) {
buf[m++] = qorient[i][0];
buf[m++] = qorient[i][1];
buf[m++] = qorient[i][2];
buf[m++] = qorient[i][3];
}
return m;
}
@ -1905,17 +2026,13 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
int m = 4;
eflags[nlocal] = static_cast<int> (buf[m++]);
for (int j = 0; j < orientflag; j++)
orient[nlocal][j] = buf[m++];
if (dorientflag) {
dorient[nlocal][0] = buf[m++];
dorient[nlocal][1] = buf[m++];
dorient[nlocal][2] = buf[m++];
}
if (qorientflag) {
qorient[nlocal][0] = buf[m++];
qorient[nlocal][1] = buf[m++];
qorient[nlocal][2] = buf[m++];
qorient[nlocal][3] = buf[m++];
}
return m;
}

11
src/fix_rigid.h
View File

@ -56,6 +56,7 @@ class FixRigid : public Fix {
double dtv,dtf,dtq;
double *step_respa;
int triclinic;
double MINUSPI,TWOPI;
int nbody; // # of rigid bodies
int *nrigid; // # of atoms in each rigid body
@ -82,11 +83,11 @@ class FixRigid : public Fix {
int **remapflag; // PBC remap flags for each rigid body
int extended; // 1 if any particles have extended attributes
int orientflag; // 1 if particles store spatial orientation
int dorientflag; // 1 if particles store dipole orientation
int qorientflag; // 1 if particles store quat orientation
int *eflags; // flags for extended particles
double **qorient; // rotation state of ext particle wrt rigid body
double **orient; // orientation vector of particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body
double tfactor; // scale factor on temperature of rigid bodies
@ -104,10 +105,10 @@ class FixRigid : public Fix {
class RanMars *random;
class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line;
class AtomVecTri *avec_tri;
// bitmasks for eflags
int INERTIA_POINT,INERTIA_SPHERE,INERTIA_ELLIPSOID;
int ORIENT_DIPOLE,ORIENT_QUAT;
int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE;
void no_squish_rotate(int, double *, double *, double *, double);

99
src/read_data.cpp
View File

@ -21,6 +21,8 @@
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "comm.h"
#include "update.h"
#include "force.h"
@ -42,7 +44,10 @@ using namespace LAMMPS_NS;
#define DELTA 4 // must be 2 or larger
// customize for new sections
#define NSECTIONS 21 // change when add to header::section_keywords
#define NSECTIONS 23 // change when add to header::section_keywords
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
/* ---------------------------------------------------------------------- */
@ -60,6 +65,10 @@ ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
nellipsoids = 0;
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
nlines = 0;
avec_line = (AtomVecLine *) atom->style_match("line");
ntris = 0;
avec_tri = (AtomVecTri *) atom->style_match("tri");
}
/* ---------------------------------------------------------------------- */
@ -151,7 +160,17 @@ void ReadData::command(int narg, char **arg)
if (!avec_ellipsoid)
error->all(FLERR,"Invalid data file section: Ellipsoids");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Ellipsoids");
ellipsoids();
bonus(nellipsoids,(AtomVec *) avec_ellipsoid,"ellipsoids");
} else if (strcmp(keyword,"Lines") == 0) {
if (!avec_line)
error->all(FLERR,"Invalid data file section: Lines");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Lines");
bonus(nlines,(AtomVec *) avec_line,"lines");
} else if (strcmp(keyword,"Triangles") == 0) {
if (!avec_tri)
error->all(FLERR,"Invalid data file section: Triangles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles");
bonus(ntris,(AtomVec *) avec_tri,"triangles");
} else if (strcmp(keyword,"Bonds") == 0) {
if (atom->avec->bonds_allow == 0)
@ -222,25 +241,31 @@ void ReadData::command(int narg, char **arg)
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before MiddleBondTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
dihedralcoeffs(1);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before EndBondTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
dihedralcoeffs(2);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before AngleTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
dihedralcoeffs(3);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before AngleAngleTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
dihedralcoeffs(4);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
@ -274,7 +299,8 @@ void ReadData::command(int narg, char **arg)
// error if natoms > 0 yet no atoms were read
if (atom->natoms > 0 && atomflag == 0) error->all(FLERR,"No atoms in data file");
if (atom->natoms > 0 && atomflag == 0)
error->all(FLERR,"No atoms in data file");
// create bond topology now that system is defined
@ -303,7 +329,7 @@ void ReadData::header(int flag)
// customize for new sections
char *section_keywords[NSECTIONS] =
{"Atoms","Velocities","Ellipsoids",
{"Atoms","Velocities","Ellipsoids","Lines","Triangles",
"Bonds","Angles","Dihedrals","Impropers",
"Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
@ -373,6 +399,14 @@ void ReadData::header(int flag)
if (!avec_ellipsoid)
error->all(FLERR,"No ellipsoids allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nellipsoids);
} else if (strstr(line,"lines")) {
if (!avec_line)
error->all(FLERR,"No lines allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nlines);
} else if (strstr(line,"triangles")) {
if (!avec_tri)
error->all(FLERR,"No triangles allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&ntris);
}
else if (strstr(line,"xlo xhi"))
@ -475,7 +509,8 @@ void ReadData::atoms()
if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms);
}
if (natoms != atom->natoms) error->all(FLERR,"Did not assign all atoms correctly");
if (natoms != atom->natoms)
error->all(FLERR,"Did not assign all atoms correctly");
// if any atom ID < 0, error
// if all atom IDs = 0, tag_enable = 0
@ -568,11 +603,11 @@ void ReadData::velocities()
}
/* ----------------------------------------------------------------------
read all ellipsoids
read all bonus data
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
void ReadData::ellipsoids()
void ReadData::bonus(bigint nbonus, AtomVec *ptr, char *type)
{
int i,m,nchunk;
@ -585,7 +620,7 @@ void ReadData::ellipsoids()
}
bigint nread = 0;
bigint natoms = nellipsoids;
bigint natoms = nbonus;
while (nread < natoms) {
if (natoms-nread > CHUNK) nchunk = CHUNK;
@ -603,7 +638,7 @@ void ReadData::ellipsoids()
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
atom->data_bonus(nchunk,buffer,avec_ellipsoid);
atom->data_bonus(nchunk,buffer,ptr);
nread += nchunk;
}
@ -613,8 +648,8 @@ void ReadData::ellipsoids()
}
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " ellipsoids\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " ellipsoids\n",natoms);
if (screen) fprintf(screen," " BIGINT_FORMAT " %s\n",natoms,type);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " %s\n",natoms,type);
}
}
@ -660,7 +695,8 @@ void ReadData::bonds()
if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor);
}
if (sum != factor*atom->nbonds) error->all(FLERR,"Bonds assigned incorrectly");
if (sum != factor*atom->nbonds)
error->all(FLERR,"Bonds assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
@ -705,7 +741,8 @@ void ReadData::angles()
if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor);
}
if (sum != factor*atom->nangles) error->all(FLERR,"Angles assigned incorrectly");
if (sum != factor*atom->nangles)
error->all(FLERR,"Angles assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
@ -1010,6 +1047,8 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
int natoms = static_cast<int> (atom->natoms);
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
int ellipsoid_flag = 0;
int line_flag = 0;
int tri_flag = 0;
// customize for new sections
// allocate topology counting vector
@ -1031,6 +1070,14 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
error->all(FLERR,"Invalid data file section: Ellipsoids");
ellipsoid_flag = 1;
skip_lines(nellipsoids);
} else if (strcmp(keyword,"Lines") == 0) {
if (!avec_line) error->all(FLERR,"Invalid data file section: Lines");
line_flag = 1;
skip_lines(nlines);
} else if (strcmp(keyword,"Triangles") == 0) {
if (!avec_tri) error->all(FLERR,"Invalid data file section: Triangles");
tri_flag = 1;
skip_lines(ntris);
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
@ -1077,25 +1124,31 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before MiddleBondTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before EndBondTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before AngleTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before AngleAngleTorsion Coeffs");
error->all(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
@ -1252,6 +1305,10 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
if (nellipsoids && !ellipsoid_flag)
error->one(FLERR,"Needed bonus data not in data file");
if (nlines && !line_flag)
error->one(FLERR,"Needed bonus data not in data file");
if (ntris && !tri_flag)
error->one(FLERR,"Needed bonus data not in data file");
}
/* ----------------------------------------------------------------------

6
src/read_data.h
View File

@ -40,6 +40,10 @@ class ReadData : protected Pointers {
bigint nellipsoids;
class AtomVecEllipsoid *avec_ellipsoid;
bigint nlines;
class AtomVecLine *avec_line;
bigint ntris;
class AtomVecTri *avec_tri;
void open(char *);
void scan(int &, int &, int &, int &);
@ -51,7 +55,7 @@ class ReadData : protected Pointers {
void atoms();
void velocities();
void ellipsoids();
void bonus(bigint, class AtomVec *, char *);
void bonds();
void angles();

136
src/set.cpp
View File

@ -19,6 +19,8 @@
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "domain.h"
#include "region.h"
#include "group.h"
@ -35,8 +37,8 @@ using namespace LAMMPS_NS;
using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,ANGMOM,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
MESO_E,MESO_CV,MESO_RHO};
@ -150,6 +152,22 @@ void Set::command(int narg, char **arg)
}
set(SHAPE);
iarg += 4;
} else if (strcmp(arg[iarg],"length") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
dvalue = atof(arg[iarg+1]);
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (dvalue < 0.0) error->all(FLERR,"Invalid length in set command");
set(LENGTH);
iarg += 2;
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
dvalue = atof(arg[iarg+1]);
if (!atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (dvalue < 0.0) error->all(FLERR,"Invalid length in set command");
set(TRI);
iarg += 2;
} else if (strcmp(arg[iarg],"dipole") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
xvalue = atof(arg[iarg+1]);
@ -177,19 +195,36 @@ void Set::command(int narg, char **arg)
yvalue = atof(arg[iarg+2]);
zvalue = atof(arg[iarg+3]);
wvalue = atof(arg[iarg+4]);
if (!atom->ellipsoid_flag)
if (!atom->ellipsoid_flag && !atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(QUAT);
iarg += 5;
} else if (strcmp(arg[iarg],"quat/random") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = atoi(arg[iarg+1]);
if (!atom->ellipsoid_flag)
if (!atom->ellipsoid_flag && !atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
setrandom(QUAT_RANDOM);
iarg += 2;
} else if (strcmp(arg[iarg],"theta") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
dvalue = atof(arg[iarg+1]);
dvalue *= MY_PI/180.0;
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(THETA);
iarg += 2;
} else if (strcmp(arg[iarg],"angmom") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
xvalue = atof(arg[iarg+1]);
yvalue = atof(arg[iarg+2]);
zvalue = atof(arg[iarg+3]);
if (!atom->ellipsoid_flag && !atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(ANGMOM);
iarg += 4;
} else if (strcmp(arg[iarg],"diameter") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
dvalue = atof(arg[iarg+1]);
@ -385,6 +420,8 @@ void Set::set(int keyword)
{
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
selection(atom->nlocal);
@ -405,23 +442,49 @@ void Set::set(int keyword)
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
// set shape
// set shape of ellipsoidal particle
else if (keyword == SHAPE)
avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue);
// set length of line particle
else if (keyword == LENGTH)
avec_line->set_length(i,dvalue);
// set corners of tri particle
else if (keyword == TRI)
avec_tri->set_equilateral(i,dvalue);
// set rmass via density
// if radius > 0.0, treat as sphere
// if shape > 0.0, treat as ellipsoid
// if length > 0.0, treat as line
// if area > 0.0, treat as tri
// else set rmass to density directly
else if (keyword == DENSITY) {
if (atom->radius_flag && atom->radius[i] > 0.0)
atom->rmass[i] = 4.0*MY_PI*THIRD *
atom->rmass[i] = 4.0*MY_PI/3.0 *
atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
atom->rmass[i] = 4.0*MY_PI*THIRD * shape[0]*shape[1]*shape[2] * dvalue;
atom->rmass[i] = 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
} else if (atom->line_flag && atom->line[i] >= 0) {
double length = avec_line->bonus[atom->line[i]].length;
atom->rmass[i] = length * dvalue;
} else if (atom->tri_flag && atom->tri[i] >= 0) {
double *c1 = avec_tri->bonus[atom->tri[i]].c1;
double *c2 = avec_tri->bonus[atom->tri[i]].c2;
double *c3 = avec_tri->bonus[atom->tri[i]].c3;
double c2mc1[2],c3mc1[3];
MathExtra::sub3(c2,c1,c2mc1);
MathExtra::sub3(c3,c1,c3mc1);
double norm[3];
MathExtra::cross3(c2mc1,c3mc1,norm);
double area = 0.5 * MathExtra::len3(norm);
atom->rmass[i] = area * dvalue;
} else atom->rmass[i] = dvalue;
// reset any or all of 3 image flags
@ -446,13 +509,17 @@ void Set::set(int keyword)
mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] +
mu[i][2]*mu[i][2]);
// set quaternion orientation
// set quaternion orientation of ellipsoid or tri particle
} else if (keyword == QUAT) {
if (atom->ellipsoid[i] < 0)
error->one(FLERR,
"Cannot set quaternion for atom that is not an ellipsoid");
double *quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
double *quat;
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
double theta2 = MY_PI2 * wvalue/180.0;
double sintheta2 = sin(theta2);
quat[0] = cos(theta2);
@ -460,7 +527,22 @@ void Set::set(int keyword)
quat[2] = yvalue * sintheta2;
quat[3] = zvalue * sintheta2;
MathExtra::qnormalize(quat);
// set theta of line particle
} else if (keyword == THETA) {
if (atom->line[i] < 0)
error->one(FLERR,"Cannot set theta for atom that is not a line");
avec_line->bonus[atom->line[i]].theta = dvalue;
// set angmom of ellipsoidal or tri particle
} else if (keyword == ANGMOM) {
atom->angmom[i][0] = xvalue;
atom->angmom[i][1] = yvalue;
atom->angmom[i][2] = zvalue;
}
count++;
}
}
@ -475,6 +557,11 @@ void Set::setrandom(int keyword)
{
int i;
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
selection(atom->nlocal);
RanPark *random = new RanPark(lmp,1);
double **x = atom->x;
@ -535,13 +622,10 @@ void Set::setrandom(int keyword)
}
// set quaternions to random orientations in 3d or 2d
// no need to normalize quats since creations algorithms already do
} else if (keyword == QUAT_RANDOM) {
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *tri = atom->tri;
int nlocal = atom->nlocal;
double *quat;
@ -549,10 +633,13 @@ void Set::setrandom(int keyword)
double s,t1,t2,theta1,theta2;
for (i = 0; i < nlocal; i++)
if (select[i]) {
if (ellipsoid[i] < 0)
error->one(FLERR,"Cannot set quaternion for atom "
"that is not an ellipsoid");
quat = bonus[ellipsoid[i]].quat;
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]);
s = random->uniform();
t1 = sqrt(1.0-s);
@ -570,10 +657,11 @@ void Set::setrandom(int keyword)
double theta2;
for (i = 0; i < nlocal; i++)
if (select[i]) {
if (ellipsoid[i] < 0)
error->one(FLERR,"Cannot set quaternion for atom "
"that is not an ellipsoid");
quat = bonus[ellipsoid[i]].quat;
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]);
theta2 = MY_PI*random->uniform();
quat[0] = cos(theta2);