fix_rigid_small.cpp
8bee4ef5966alammm-devel

fix_rigid_small.cpp
View Options

			/* ----------------------------------------------------------------------
			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.
			------------------------------------------------------------------------- */

			#include "math.h"
			#include "stdio.h"
			#include "stdlib.h"
			#include "string.h"
			#include "fix_rigid_small.h"
			#include "math_extra.h"
			#include "atom.h"
			#include "atom_vec_ellipsoid.h"
			#include "atom_vec_line.h"
			#include "atom_vec_tri.h"
			#include "molecule.h"
			#include "domain.h"
			#include "update.h"
			#include "respa.h"
			#include "modify.h"
			#include "group.h"
			#include "comm.h"
			#include "force.h"
			#include "output.h"
			#include "random_mars.h"
			#include "math_const.h"
			#include "memory.h"
			#include "error.h"

			#include <map>

			using namespace LAMMPS_NS;
			using namespace FixConst;
			using namespace MathConst;

			// allocate space for static class variable

			FixRigidSmall *FixRigidSmall::frsptr;

			#define MAXLINE 256
			#define CHUNK 1024
			#define ATTRIBUTE_PERBODY 11

			#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.4 // 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) :
			Fix(lmp, narg, arg)
			{
			int i;

			scalar_flag = 1;
			extscalar = 0;
			global_freq = 1;
			time_integrate = 1;
			rigid_flag = 1;
			virial_flag = 1;
			create_attribute = 1;

			MPI_Comm_rank(world,&me);
			MPI_Comm_size(world,&nprocs);

			// perform initial allocation of atom-based arrays
			// register with Atom class

			extended = orientflag = dorientflag = 0;
			bodyown = NULL;
			bodytag = NULL;
			atom2body = NULL;
			displace = NULL;
			eflags = NULL;
			orient = NULL;
			dorient = NULL;
			grow_arrays(atom->nmax);
			atom->add_callback(0);

			// 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");

			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 #

			int *mask = atom->mask;
			tagint *molecule = atom->molecule;
			int nlocal = atom->nlocal;

			maxmol = -1;
			for (i = 0; i < nlocal; i++)
			if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);

			tagint itmp;
			MPI_Allreduce(&maxmol,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
			maxmol = itmp;

			// parse optional args

			int seed;
			langflag = 0;
			infile = NULL;
			onemol = NULL;

			tstat_flag = 0;
			pstat_flag = 0;
			allremap = 1;
			id_dilate = NULL;
			t_chain = 10;
			t_iter = 1;
			t_order = 3;
			p_chain = 10;

			pcouple = NONE;
			pstyle = ANISO;

			for (int i = 0; i < 3; i++) {
			p_start[i] = p_stop[i] = p_period[i] = 0.0;
			p_flag[i] = 0;
			}

			int iarg = 4;
			while (iarg < narg) {
			if (strcmp(arg[iarg],"langevin") == 0) {
			if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/small") != 0)
			error->all(FLERR,"Illegal fix rigid/small command");
			langflag = 1;
			t_start = force->numeric(FLERR,arg[iarg+1]);
			t_stop = force->numeric(FLERR,arg[iarg+2]);
			t_period = force->numeric(FLERR,arg[iarg+3]);
			seed = force->inumeric(FLERR,arg[iarg+4]);
			if (t_period <= 0.0)
			error->all(FLERR,"Fix rigid/small langevin period must be > 0.0");
			if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
			iarg += 5;
			} else if (strcmp(arg[iarg],"infile") == 0) {
			if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			delete [] infile;
			int n = strlen(arg[iarg+1]) + 1;
			infile = new char[n];
			strcpy(infile,arg[iarg+1]);
			restart_file = 1;
			iarg += 2;
			} else if (strcmp(arg[iarg],"mol") == 0) {
			if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			int imol = atom->find_molecule(arg[iarg+1]);
			if (imol == -1)
			error->all(FLERR,"Molecule template ID for "
			"fix rigid/small does not exist");
			if (atom->molecules[imol]->nset > 1 && comm->me == 0)
			error->warning(FLERR,"Molecule template for "
			"fix rigid/small has multiple molecules");
			onemol = atom->molecules[imol];
			iarg += 2;

			} else if (strcmp(arg[iarg],"temp") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/nvt/small") != 0 &&
			strcmp(style,"rigid/npt/small") != 0)
			error->all(FLERR,"Illegal fix rigid command");
			tstat_flag = 1;
			t_start = force->numeric(FLERR,arg[iarg+1]);
			t_stop = force->numeric(FLERR,arg[iarg+2]);
			t_period = force->numeric(FLERR,arg[iarg+3]);
			iarg += 4;

			} else if (strcmp(arg[iarg],"iso") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/npt/small") != 0 &&
			strcmp(style,"rigid/nph/small") != 0)
			error->all(FLERR,"Illegal fix rigid/small command");
			pcouple = XYZ;
			p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
			p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
			p_period[0] = p_period[1] = p_period[2] =
			force->numeric(FLERR,arg[iarg+3]);
			p_flag[0] = p_flag[1] = p_flag[2] = 1;
			if (domain->dimension == 2) {
			p_start[2] = p_stop[2] = p_period[2] = 0.0;
			p_flag[2] = 0;
			}
			iarg += 4;

			} else if (strcmp(arg[iarg],"aniso") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/npt/small") != 0 &&
			strcmp(style,"rigid/nph/small") != 0)
			error->all(FLERR,"Illegal fix rigid/small command");
			p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
			p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
			p_period[0] = p_period[1] = p_period[2] =
			force->numeric(FLERR,arg[iarg+3]);
			p_flag[0] = p_flag[1] = p_flag[2] = 1;
			if (domain->dimension == 2) {
			p_start[2] = p_stop[2] = p_period[2] = 0.0;
			p_flag[2] = 0;
			}
			iarg += 4;

			} else if (strcmp(arg[iarg],"x") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			p_start[0] = force->numeric(FLERR,arg[iarg+1]);
			p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
			p_period[0] = force->numeric(FLERR,arg[iarg+3]);
			p_flag[0] = 1;
			iarg += 4;

			} else if (strcmp(arg[iarg],"y") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			p_start[1] = force->numeric(FLERR,arg[iarg+1]);
			p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
			p_period[1] = force->numeric(FLERR,arg[iarg+3]);
			p_flag[1] = 1;
			iarg += 4;

			} else if (strcmp(arg[iarg],"z") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			p_start[2] = force->numeric(FLERR,arg[iarg+1]);
			p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
			p_period[2] = force->numeric(FLERR,arg[iarg+3]);
			p_flag[2] = 1;
			iarg += 4;

			} else if (strcmp(arg[iarg],"couple") == 0) {
			if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
			else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
			else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
			else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
			else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
			else error->all(FLERR,"Illegal fix rigid/small command");
			iarg += 2;

			} else if (strcmp(arg[iarg],"dilate") == 0) {
			if (iarg+2 > narg)
			error->all(FLERR,"Illegal fix rigid/small nvt/npt/nph command");
			if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
			else {
			allremap = 0;
			delete [] id_dilate;
			int n = strlen(arg[iarg+1]) + 1;
			id_dilate = new char[n];
			strcpy(id_dilate,arg[iarg+1]);
			int idilate = group->find(id_dilate);
			if (idilate == -1)
			error->all(FLERR,"Fix rigid/small nvt/npt/nph dilate group ID "
			"does not exist");
			}
			iarg += 2;

			} else if (strcmp(arg[iarg],"tparam") == 0) {
			if (iarg+4 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/nvt/small") != 0 &&
			strcmp(style,"rigid/npt/small") != 0)
			error->all(FLERR,"Illegal fix rigid/small command");
			t_chain = force->numeric(FLERR,arg[iarg+1]);
			t_iter = force->numeric(FLERR,arg[iarg+2]);
			t_order = force->numeric(FLERR,arg[iarg+3]);
			iarg += 4;

			} else if (strcmp(arg[iarg],"pchain") == 0) {
			if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
			if (strcmp(style,"rigid/npt/small") != 0 &&
			strcmp(style,"rigid/nph/small") != 0)
			error->all(FLERR,"Illegal fix rigid/small command");
			p_chain = force->numeric(FLERR,arg[iarg+1]);
			iarg += 2;


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

			// error check and further setup for Molecule template

			if (onemol) {
			if (onemol->xflag == 0)
			error->all(FLERR,"Fix rigid/small molecule must have coordinates");
			if (onemol->typeflag == 0)
			error->all(FLERR,"Fix rigid/small molecule must have atom types");

			// fix rigid/small uses center, masstotal, COM, inertia of molecule

			onemol->compute_center();
			onemol->compute_mass();
			onemol->compute_com();
			onemol->compute_inertia();
			}

			// set pstat_flag

			pstat_flag = 0;
			for (int i = 0; i < 3; i++)
			if (p_flag[i]) pstat_flag = 1;

			if (pcouple == XYZ \|\| (domain->dimension == 2 && pcouple == XY)) pstyle = ISO;
			else pstyle = ANISO;

			// create rigid bodies based on molecule ID
			// sets bodytag for owned atoms
			// body attributes are computed later by setup_bodies()

			create_bodies();

			// set nlocal_body and allocate bodies I own

			tagint *tag = atom->tag;

			nlocal_body = nghost_body = 0;
			for (i = 0; i < nlocal; i++)
			if (bodytag[i] == tag[i]) nlocal_body++;

			nmax_body = 0;
			while (nmax_body < nlocal_body) nmax_body += DELTA_BODY;
			body = (Body ) memory->smalloc(nmax_bodysizeof(Body),
			"rigid/small:body");

			// set bodyown for owned atoms

			nlocal_body = 0;
			for (i = 0; i < nlocal; i++)
			if (bodytag[i] == tag[i]) {
			body[nlocal_body].ilocal = i;
			bodyown[i] = nlocal_body++;
			} else bodyown[i] = -1;


			// bodysize = sizeof(Body) in doubles

			bodysize = sizeof(Body)/sizeof(double);
			if (bodysize*sizeof(double) != sizeof(Body)) bodysize++;

			// set max comm sizes needed by this fix

			comm_forward = 1 + bodysize;
			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

			avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
			avec_line = (AtomVecLine *) atom->style_match("line");
			avec_tri = (AtomVecTri *) atom->style_match("tri");

			// print statistics

			int one = 0;
			bigint atomone = 0;
			for (int i = 0; i < nlocal; i++) {
			if (bodyown[i] >= 0) one++;
			if (bodytag[i] > 0) atomone++;
			}
			MPI_Allreduce(&one,&nbody,1,MPI_INT,MPI_SUM,world);
			bigint atomall;
			MPI_Allreduce(&atomone,&atomall,1,MPI_LMP_BIGINT,MPI_SUM,world);

			if (me == 0) {
			if (screen) {
			fprintf(screen,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
			nbody,atomall);
			fprintf(screen," %g = max distance from body owner to body atom\n",
			maxextent);
			}
			if (logfile) {
			fprintf(logfile,"%d rigid bodies with " BIGINT_FORMAT " atoms\n",
			nbody,atomall);
			fprintf(logfile," %g = max distance from body owner to body atom\n",
			maxextent);
			}
			}

			// initialize Marsaglia RNG with processor-unique seed

			maxlang = 0;
			langextra = NULL;
			random = NULL;
			if (langflag) random = new RanMars(lmp,seed + comm->me);

			// mass vector for granular pair styles

			mass_body = NULL;
			nmax_mass = 0;

			// firstflag = 1 triggers one-time initialization of rigid body attributes

			firstflag = 1;
			}

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

			FixRigidSmall::~FixRigidSmall()
			{
			// unregister callbacks to this fix from Atom class

			atom->delete_callback(id,0);

			// delete locally stored arrays

			memory->sfree(body);

			memory->destroy(bodyown);
			memory->destroy(bodytag);
			memory->destroy(atom2body);
			memory->destroy(displace);
			memory->destroy(eflags);
			memory->destroy(orient);
			memory->destroy(dorient);

			delete random;
			delete [] infile;

			memory->destroy(langextra);
			memory->destroy(mass_body);
			}

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

			int FixRigidSmall::setmask()
			{
			int mask = 0;
			mask \|= INITIAL_INTEGRATE;
			mask \|= FINAL_INTEGRATE;
			if (langflag) mask \|= POST_FORCE;
			mask \|= PRE_NEIGHBOR;
			mask \|= INITIAL_INTEGRATE_RESPA;
			mask \|= FINAL_INTEGRATE_RESPA;
			return mask;
			}

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

			void FixRigidSmall::init()
			{
			int i;

			triclinic = domain->triclinic;

			// warn if more than one rigid fix

			int count = 0;
			for (i = 0; i < modify->nfix; i++)
			if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
			if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid");

			// error if npt,nph fix comes before rigid fix

			for (i = 0; i < modify->nfix; i++) {
			if (strcmp(modify->fix[i]->style,"npt") == 0) break;
			if (strcmp(modify->fix[i]->style,"nph") == 0) break;
			}
			if (i < modify->nfix) {
			for (int j = i; j < modify->nfix; j++)
			if (strcmp(modify->fix[j]->style,"rigid") == 0)
			error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
			}

			// timestep info

			dtv = update->dt;
			dtf = 0.5 * update->dt * force->ftm2v;
			dtq = 0.5 * update->dt;

			if (strstr(update->integrate_style,"respa"))
			step_respa = ((Respa *) update->integrate)->step;
			}

			/* ----------------------------------------------------------------------
			one-time initialization of rigid body attributes via local comm
			extended flags, mass, COM, inertia tensor, displacement of each atom
			performed after init() b/c requires communication stencil
			has been setup by comm->borders()
			------------------------------------------------------------------------- */

			void FixRigidSmall::setup_pre_neighbor()
			{
			if (firstflag) {
			firstflag = 0;
			setup_bodies_static();
			setup_bodies_dynamic();
			} else pre_neighbor();
			}

			/* ----------------------------------------------------------------------
			compute initial fcm and torque on bodies, also initial virial
			reset all particle velocities to be consistent with vcm and omega
			------------------------------------------------------------------------- */

			void FixRigidSmall::setup(int vflag)
			{
			int i,n,ibody;

			//check(1);

			// sum fcm, torque across all rigid bodies
			// fcm = force on COM
			// torque = torque around COM

			double **x = atom->x;
			double **f = atom->f;
			imageint *image = atom->image;
			int nlocal = atom->nlocal;

			double xcm,fcm,*tcm;
			double dx,dy,dz;
			double unwrap[3];

			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			fcm = body[ibody].fcm;
			fcm[0] = fcm[1] = fcm[2] = 0.0;
			tcm = body[ibody].torque;
			tcm[0] = tcm[1] = tcm[2] = 0.0;
			}

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			fcm = b->fcm;
			fcm[0] += f[i][0];
			fcm[1] += f[i][1];
			fcm[2] += f[i][2];

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			dx = unwrap[0] - xcm[0];
			dy = unwrap[1] - xcm[1];
			dz = unwrap[2] - xcm[2];

			tcm = b->torque;
			tcm[0] += dy * f[i][2] - dz * f[i][1];
			tcm[1] += dz * f[i][0] - dx * f[i][2];
			tcm[2] += dx * f[i][1] - dy * f[i][0];
			}

			// extended particles add their rotation/torque to angmom/torque of body

			if (extended) {
			double **torque = atom->torque;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];
			if (eflags[i] & TORQUE) {
			tcm = b->torque;
			tcm[0] += torque[i][0];
			tcm[1] += torque[i][1];
			tcm[2] += torque[i][2];
			}
			}
			}

			// reverse communicate fcm, torque of all bodies

			commflag = FORCE_TORQUE;
			comm->reverse_comm_variable_fix(this);

			// virial setup before call to set_v

			if (vflag) v_setup(vflag);
			else evflag = 0;

			// compute and forward communicate vcm and omega of all bodies

			for (ibody = 0; ibody < nlocal_body; ibody++) {
			Body *b = &body[ibody];
			MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
			b->ez_space,b->inertia,b->omega);
			}

			commflag = FINAL;
			comm->forward_comm_variable_fix(this);

			// set velocity/rotation of atoms in rigid bodues

			set_v();

			// guesstimate virial as 2x the set_v contribution

			if (vflag_global)
			for (n = 0; n < 6; n++) virial[n] *= 2.0;
			if (vflag_atom) {
			for (i = 0; i < nlocal; i++)
			for (n = 0; n < 6; n++)
			vatom[i][n] *= 2.0;
			}
			}

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

			void FixRigidSmall::initial_integrate(int vflag)
			{
			double dtfm;

			//check(2);

			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			Body *b = &body[ibody];

			// update vcm by 1/2 step

			dtfm = dtf / b->mass;
			b->vcm[0] += dtfm * b->fcm[0];
			b->vcm[1] += dtfm * b->fcm[1];
			b->vcm[2] += dtfm * b->fcm[2];

			// update xcm by full step

			b->xcm[0] += dtv * b->vcm[0];
			b->xcm[1] += dtv * b->vcm[1];
			b->xcm[2] += dtv * b->vcm[2];

			// update angular momentum by 1/2 step

			b->angmom[0] += dtf * b->torque[0];
			b->angmom[1] += dtf * b->torque[1];
			b->angmom[2] += dtf * b->torque[2];

			// compute omega at 1/2 step from angmom at 1/2 step and current q
			// update quaternion a full step via Richardson iteration
			// returns new normalized quaternion, also updated omega at 1/2 step
			// update ex,ey,ez to reflect new quaternion

			MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
			b->ez_space,b->inertia,b->omega);
			MathExtra::richardson(b->quat,b->angmom,b->omega,b->inertia,dtq);
			MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
			}

			// virial setup before call to set_xv

			if (vflag) v_setup(vflag);
			else evflag = 0;

			// forward communicate updated info of all bodies

			commflag = INITIAL;
			comm->forward_comm_variable_fix(this);

			// set coords/orient and velocity/rotation of atoms in rigid bodies

			set_xv();
			}

			/* ----------------------------------------------------------------------
			apply Langevin thermostat to all 6 DOF of rigid bodies I own
			unlike fix langevin, this stores extra force in extra arrays,
			which are added in when final_integrate() calculates a new fcm/torque
			------------------------------------------------------------------------- */

			void FixRigidSmall::post_force(int vflag)
			{
			double gamma1,gamma2;

			// grow langextra if needed

			if (nlocal_body > maxlang) {
			memory->destroy(langextra);
			maxlang = nlocal_body + nghost_body;
			memory->create(langextra,maxlang,6,"rigid/small:langextra");
			}

			double delta = update->ntimestep - update->beginstep;
			delta /= update->endstep - update->beginstep;
			double t_target = t_start + delta * (t_stop-t_start);
			double tsqrt = sqrt(t_target);

			double boltz = force->boltz;
			double dt = update->dt;
			double mvv2e = force->mvv2e;
			double ftm2v = force->ftm2v;

			double vcm,omega,*inertia;

			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			vcm = body[ibody].vcm;
			omega = body[ibody].omega;
			inertia = body[ibody].inertia;

			gamma1 = -body[ibody].mass / t_period / ftm2v;
			gamma2 = sqrt(body[ibody].mass) * tsqrt *
			sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
			langextra[ibody][0] = gamma1vcm[0] + gamma2(random->uniform()-0.5);
			langextra[ibody][1] = gamma1vcm[1] + gamma2(random->uniform()-0.5);
			langextra[ibody][2] = gamma1vcm[2] + gamma2(random->uniform()-0.5);

			gamma1 = -1.0 / t_period / ftm2v;
			gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
			langextra[ibody][3] = inertia[0]gamma1omega[0] +
			sqrt(inertia[0])gamma2(random->uniform()-0.5);
			langextra[ibody][4] = inertia[1]gamma1omega[1] +
			sqrt(inertia[1])gamma2(random->uniform()-0.5);
			langextra[ibody][5] = inertia[2]gamma1omega[2] +
			sqrt(inertia[2])gamma2(random->uniform()-0.5);
			}
			}

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

			void FixRigidSmall::final_integrate()
			{
			int i,ibody;
			double dtfm;

			//check(3);

			// sum over atoms to get force and torque on rigid body

			imageint *image = atom->image;
			double **x = atom->x;
			double **f = atom->f;
			int nlocal = atom->nlocal;

			double dx,dy,dz;
			double unwrap[3];
			double xcm,fcm,*tcm;

			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			fcm = body[ibody].fcm;
			fcm[0] = fcm[1] = fcm[2] = 0.0;
			tcm = body[ibody].torque;
			tcm[0] = tcm[1] = tcm[2] = 0.0;
			}

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			fcm = b->fcm;
			fcm[0] += f[i][0];
			fcm[1] += f[i][1];
			fcm[2] += f[i][2];

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			dx = unwrap[0] - xcm[0];
			dy = unwrap[1] - xcm[1];
			dz = unwrap[2] - xcm[2];

			tcm = b->torque;
			tcm[0] += dyf[i][2] - dzf[i][1];
			tcm[1] += dzf[i][0] - dxf[i][2];
			tcm[2] += dxf[i][1] - dyf[i][0];
			}

			// extended particles add their torque to torque of body

			if (extended) {
			double **torque = atom->torque;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;

			if (eflags[i] & TORQUE) {
			tcm = body[atom2body[i]].torque;
			tcm[0] += torque[i][0];
			tcm[1] += torque[i][1];
			tcm[2] += torque[i][2];
			}
			}
			}

			// reverse communicate fcm, torque of all bodies

			commflag = FORCE_TORQUE;
			comm->reverse_comm_variable_fix(this);

			// include Langevin thermostat forces and torques

			if (langflag) {
			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			fcm = body[ibody].fcm;
			fcm[0] += langextra[ibody][0];
			fcm[1] += langextra[ibody][1];
			fcm[2] += langextra[ibody][2];
			tcm = body[ibody].torque;
			tcm[0] += langextra[ibody][3];
			tcm[1] += langextra[ibody][4];
			tcm[2] += langextra[ibody][5];
			}
			}

			// update vcm and angmom, recompute omega

			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			Body *b = &body[ibody];

			// update vcm by 1/2 step

			dtfm = dtf / b->mass;
			b->vcm[0] += dtfm * b->fcm[0];
			b->vcm[1] += dtfm * b->fcm[1];
			b->vcm[2] += dtfm * b->fcm[2];

			// update angular momentum by 1/2 step

			b->angmom[0] += dtf * b->torque[0];
			b->angmom[1] += dtf * b->torque[1];
			b->angmom[2] += dtf * b->torque[2];

			MathExtra::angmom_to_omega(b->angmom,b->ex_space,b->ey_space,
			b->ez_space,b->inertia,b->omega);
			}

			// forward communicate updated info of all bodies

			commflag = FINAL;
			comm->forward_comm_variable_fix(this);

			// set velocity/rotation of atoms in rigid bodies
			// virial is already setup from initial_integrate

			set_v();
			}

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

			void FixRigidSmall::initial_integrate_respa(int vflag, int ilevel, int iloop)
			{
			dtv = step_respa[ilevel];
			dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
			dtq = 0.5 * step_respa[ilevel];

			if (ilevel == 0) initial_integrate(vflag);
			else final_integrate();
			}

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

			void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
			{
			dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
			final_integrate();
			}

			/* ----------------------------------------------------------------------
			atom to processor assignments have changed,
			so acquire ghost bodies and setup atom2body
			remap xcm of each rigid body back into periodic simulation box
			done during pre_neighbor so will be after call to pbc()
			and after fix_deform::pre_exchange() may have flipped box
			use domain->remap() in case xcm is far away from box
			due to 1st definition of rigid body or due to box flip
			if don't do this, then atoms of a body which drifts far away
			from a triclinic box will be remapped back into box
			with huge displacements when the box tilt changes via set_x()
			adjust image flag of body and image flags of all atoms in body
			------------------------------------------------------------------------- */

			void FixRigidSmall::pre_neighbor()
			{
			// remap xcm and image flags of each body as needed

			imageint original,oldimage,newimage;

			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			Body *b = &body[ibody];

			original = b->image;
			domain->remap(b->xcm,b->image);

			if (original == b->image) b->remapflag[3] = 0;
			else {
			oldimage = original & IMGMASK;
			newimage = b->image & IMGMASK;
			b->remapflag[0] = newimage - oldimage;
			oldimage = (original >> IMGBITS) & IMGMASK;
			newimage = (b->image >> IMGBITS) & IMGMASK;
			b->remapflag[1] = newimage - oldimage;
			oldimage = original >> IMG2BITS;
			newimage = b->image >> IMG2BITS;
			b->remapflag[2] = newimage - oldimage;
			b->remapflag[3] = 1;
			}
			}

			// acquire ghost bodies via forward comm
			// also gets new remapflags needed for adjusting atom image flags
			// reset atom2body for owned atoms
			// forward comm sets atom2body for ghost atoms

			nghost_body = 0;
			commflag = FULL_BODY;
			comm->forward_comm_variable_fix(this);
			reset_atom2body();
			//check(4);

			// adjust image flags of any atom in a rigid body whose xcm was remapped

			imageint *image = atom->image;
			int nlocal = atom->nlocal;

			imageint idim,otherdims;

			for (int i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (b->remapflag[3] == 0) continue;

			if (b->remapflag[0]) {
			idim = image[i] & IMGMASK;
			otherdims = image[i] ^ idim;
			idim -= b->remapflag[0];
			idim &= IMGMASK;
			image[i] = otherdims \| idim;
			}
			if (b->remapflag[1]) {
			idim = (image[i] >> IMGBITS) & IMGMASK;
			otherdims = image[i] ^ (idim << IMGBITS);
			idim -= b->remapflag[1];
			idim &= IMGMASK;
			image[i] = otherdims \| (idim << IMGBITS);
			}
			if (b->remapflag[2]) {
			idim = image[i] >> IMG2BITS;
			otherdims = image[i] ^ (idim << IMG2BITS);
			idim -= b->remapflag[2];
			idim &= IMGMASK;
			image[i] = otherdims \| (idim << IMG2BITS);
			}
			}
			}

			/* ----------------------------------------------------------------------
			count # of DOF removed by rigid bodies for atoms in igroup
			return total count of DOF
			------------------------------------------------------------------------- */

			int FixRigidSmall::dof(int tgroup)
			{
			int i,j;

			// cannot count DOF correctly unless setup_bodies_static() has been called

			if (firstflag) {
			if (comm->me == 0)
			error->warning(FLERR,"Cannot count rigid body degrees-of-freedom "
			"before bodies are fully initialized");
			return 0;
			}

			int tgroupbit = group->bitmask[tgroup];

			// counts = 3 values per rigid body I own
			// 0 = # of point particles in rigid body and in temperature group
			// 1 = # of finite-size particles in rigid body and in temperature group
			// 2 = # of particles in rigid body, disregarding temperature group

			memory->create(counts,nlocal_body+nghost_body,3,"rigid/small:counts");
			for (int i = 0; i < nlocal_body+nghost_body; i++)
			counts[i][0] = counts[i][1] = counts[i][2] = 0;

			// tally counts from my owned atoms
			// 0 = # of point particles in rigid body and in temperature group
			// 1 = # of finite-size particles in rigid body and in temperature group
			// 2 = # of particles in rigid body, disregarding temperature group

			int *mask = atom->mask;
			int nlocal = atom->nlocal;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			j = atom2body[i];
			counts[j][2]++;
			if (mask[i] & tgroupbit) {
			if (extended && eflags[i]) counts[j][1]++;
			else counts[j][0]++;
			}
			}

			commflag = DOF;
			comm->reverse_comm_variable_fix(this);

			// nall = count0 = # of point particles in each rigid body
			// mall = count1 = # of finite-size particles in each rigid body
			// warn if nall+mall != nrigid for any body included in temperature group

			int flag = 0;
			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			if (counts[ibody][0]+counts[ibody][1] > 0 &&
			counts[ibody][0]+counts[ibody][1] != counts[ibody][2]) flag = 1;
			}
			int flagall;
			MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
			if (flagall && me == 0)
			error->warning(FLERR,"Computing temperature of portions of rigid bodies");

			// remove appropriate DOFs for each rigid body wholly in temperature group
			// N = # of point particles in body
			// M = # of finite-size particles in body
			// 3d body has 3N + 6M dof to start with
			// 2d body has 2N + 3M dof to start with
			// 3d point-particle body with all non-zero I should have 6 dof, remove 3N-6
			// 3d point-particle body (linear) with a 0 I should have 5 dof, remove 3N-5
			// 2d point-particle body should have 3 dof, remove 2N-3
			// 3d body with any finite-size M should have 6 dof, remove (3N+6M) - 6
			// 2d body with any finite-size M should have 3 dof, remove (2N+3M) - 3

			double *inertia;

			int n = 0;
			if (domain->dimension == 3) {
			for (int ibody = 0; ibody < nlocal_body; ibody++) {
			if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2]) {
			n += 3counts[ibody][0] + 6counts[ibody][1] - 6;
			inertia = body[ibody].inertia;
			if (inertia[0] == 0.0 \|\| inertia[1] == 0.0 \|\| inertia[2] == 0.0) n++;
			}
			}
			} else if (domain->dimension == 2) {
			for (int ibody = 0; ibody < nlocal_body; ibody++)
			if (counts[ibody][0]+counts[ibody][1] == counts[ibody][2])
			n += 2counts[ibody][0] + 3counts[ibody][1] - 3;
			}

			memory->destroy(counts);

			int nall;
			MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
			return nall;
			}

			/* ----------------------------------------------------------------------
			adjust xcm of each rigid body due to box deformation
			called by various fixes that change box size/shape
			flag = 0/1 means map from box to lamda coords or vice versa
			------------------------------------------------------------------------- */

			void FixRigidSmall::deform(int flag)
			{
			if (flag == 0)
			for (int ibody = 0; ibody < nlocal_body; ibody++)
			domain->x2lamda(body[ibody].xcm,body[ibody].xcm);
			else
			for (int ibody = 0; ibody < nlocal_body; ibody++)
			domain->lamda2x(body[ibody].xcm,body[ibody].xcm);
			}

			/* ----------------------------------------------------------------------
			set space-frame coords and velocity of each atom in each rigid body
			set orientation and rotation of extended particles
			x = Q displace + Xcm, mapped back to periodic box
			v = Vcm + (W cross (x - Xcm))
			------------------------------------------------------------------------- */

			void FixRigidSmall::set_xv()
			{
			int xbox,ybox,zbox;
			double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
			double ione[3],exone[3],eyone[3],ezone[3],vr[6],p[3][3];

			double xprd = domain->xprd;
			double yprd = domain->yprd;
			double zprd = domain->zprd;
			double xy = domain->xy;
			double xz = domain->xz;
			double yz = domain->yz;

			imageint *image = atom->image;
			double **x = atom->x;
			double **v = atom->v;
			double **f = atom->f;
			double *rmass = atom->rmass;
			double *mass = atom->mass;
			int *type = atom->type;
			int nlocal = atom->nlocal;

			// set x and v of each atom

			for (int i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			xbox = (image[i] & IMGMASK) - IMGMAX;
			ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
			zbox = (image[i] >> IMG2BITS) - IMGMAX;

			// save old positions and velocities for virial

			if (evflag) {
			if (triclinic == 0) {
			x0 = x[i][0] + xbox*xprd;
			x1 = x[i][1] + ybox*yprd;
			x2 = x[i][2] + zbox*zprd;
			} else {
			x0 = x[i][0] + xboxxprd + yboxxy + zbox*xz;
			x1 = x[i][1] + yboxyprd + zboxyz;
			x2 = x[i][2] + zbox*zprd;
			}
			v0 = v[i][0];
			v1 = v[i][1];
			v2 = v[i][2];
			}

			// x = displacement from center-of-mass, based on body orientation
			// v = vcm + omega around center-of-mass

			MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],x[i]);

			v[i][0] = b->omega[1]x[i][2] - b->omega[2]x[i][1] + b->vcm[0];
			v[i][1] = b->omega[2]x[i][0] - b->omega[0]x[i][2] + b->vcm[1];
			v[i][2] = b->omega[0]x[i][1] - b->omega[1]x[i][0] + b->vcm[2];

			// add center of mass to displacement
			// map back into periodic box via xbox,ybox,zbox
			// for triclinic, add in box tilt factors as well

			if (triclinic == 0) {
			x[i][0] += b->xcm[0] - xbox*xprd;
			x[i][1] += b->xcm[1] - ybox*yprd;
			x[i][2] += b->xcm[2] - zbox*zprd;
			} else {
			x[i][0] += b->xcm[0] - xboxxprd - yboxxy - zbox*xz;
			x[i][1] += b->xcm[1] - yboxyprd - zboxyz;
			x[i][2] += b->xcm[2] - zbox*zprd;
			}

			// virial = unwrapped coords dotted into body constraint force
			// body constraint force = implied force due to v change minus f external
			// assume f does not include forces internal to body
			// 1/2 factor b/c final_integrate contributes other half
			// assume per-atom contribution is due to constraint force on that atom

			if (evflag) {
			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];
			fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
			fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
			fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];

			vr[0] = 0.5x0fc0;
			vr[1] = 0.5x1fc1;
			vr[2] = 0.5x2fc2;
			vr[3] = 0.5x0fc1;
			vr[4] = 0.5x0fc2;
			vr[5] = 0.5x1fc2;

			v_tally(1,&i,1.0,vr);
			}
			}

			// set orientation, omega, angmom of each extended particle

			if (extended) {
			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 = atom->omega;
			double **angmom = 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 (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (eflags[i] & SPHERE) {
			omega[i][0] = b->omega[0];
			omega[i][1] = b->omega[1];
			omega[i][2] = b->omega[2];
			} else if (eflags[i] & ELLIPSOID) {
			shape = ebonus[ellipsoid[i]].shape;
			quatatom = ebonus[ellipsoid[i]].quat;
			MathExtra::quatquat(b->quat,orient[i],quatatom);
			MathExtra::qnormalize(quatatom);
			ione[0] = EINERTIArmass[i] (shape[1]shape[1] + shape[2]shape[2]);
			ione[1] = EINERTIArmass[i] (shape[0]shape[0] + shape[2]shape[2]);
			ione[2] = EINERTIArmass[i] (shape[0]shape[0] + shape[1]shape[1]);
			MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
			MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,angmom[i]);
			} else if (eflags[i] & LINE) {
			if (b->quat[3] >= 0.0) theta_body = 2.0*acos(b->quat[0]);
			else theta_body = -2.0*acos(b->quat[0]);
			theta = orient[i][0] + theta_body;
			while (theta <= MINUSPI) theta += TWOPI;
			while (theta > MY_PI) theta -= TWOPI;
			lbonus[line[i]].theta = theta;
			omega[i][0] = b->omega[0];
			omega[i][1] = b->omega[1];
			omega[i][2] = b->omega[2];
			} else if (eflags[i] & TRIANGLE) {
			inertiaatom = tbonus[tri[i]].inertia;
			quatatom = tbonus[tri[i]].quat;
			MathExtra::quatquat(b->quat,orient[i],quatatom);
			MathExtra::qnormalize(quatatom);
			MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
			MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,
			inertiaatom,angmom[i]);
			}
			if (eflags[i] & DIPOLE) {
			MathExtra::quat_to_mat(b->quat,p);
			MathExtra::matvec(p,dorient[i],mu[i]);
			MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
			}
			}
			}
			}

			/* ----------------------------------------------------------------------
			set space-frame velocity of each atom in a rigid body
			set omega and angmom of extended particles
			v = Vcm + (W cross (x - Xcm))
			------------------------------------------------------------------------- */

			void FixRigidSmall::set_v()
			{
			int xbox,ybox,zbox;
			double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
			double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6];

			double xprd = domain->xprd;
			double yprd = domain->yprd;
			double zprd = domain->zprd;
			double xy = domain->xy;
			double xz = domain->xz;
			double yz = domain->yz;

			double **x = atom->x;
			double **v = atom->v;
			double **f = atom->f;
			double *rmass = atom->rmass;
			double *mass = atom->mass;
			int *type = atom->type;
			imageint *image = atom->image;
			int nlocal = atom->nlocal;

			// set v of each atom

			for (int i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			MathExtra::matvec(b->ex_space,b->ey_space,b->ez_space,displace[i],delta);

			// save old velocities for virial

			if (evflag) {
			v0 = v[i][0];
			v1 = v[i][1];
			v2 = v[i][2];
			}

			v[i][0] = b->omega[1]delta[2] - b->omega[2]delta[1] + b->vcm[0];
			v[i][1] = b->omega[2]delta[0] - b->omega[0]delta[2] + b->vcm[1];
			v[i][2] = b->omega[0]delta[1] - b->omega[1]delta[0] + b->vcm[2];

			// virial = unwrapped coords dotted into body constraint force
			// body constraint force = implied force due to v change minus f external
			// assume f does not include forces internal to body
			// 1/2 factor b/c initial_integrate contributes other half
			// assume per-atom contribution is due to constraint force on that atom

			if (evflag) {
			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];
			fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
			fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
			fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];

			xbox = (image[i] & IMGMASK) - IMGMAX;
			ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
			zbox = (image[i] >> IMG2BITS) - IMGMAX;

			if (triclinic == 0) {
			x0 = x[i][0] + xbox*xprd;
			x1 = x[i][1] + ybox*yprd;
			x2 = x[i][2] + zbox*zprd;
			} else {
			x0 = x[i][0] + xboxxprd + yboxxy + zbox*xz;
			x1 = x[i][1] + yboxyprd + zboxyz;
			x2 = x[i][2] + zbox*zprd;
			}

			vr[0] = 0.5x0fc0;
			vr[1] = 0.5x1fc1;
			vr[2] = 0.5x2fc2;
			vr[3] = 0.5x0fc1;
			vr[4] = 0.5x0fc2;
			vr[5] = 0.5x1fc2;

			v_tally(1,&i,1.0,vr);
			}
			}

			// set omega, angmom of each extended particle

			if (extended) {
			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 = atom->omega;
			double **angmom = atom->angmom;
			int *ellipsoid = atom->ellipsoid;
			int *tri = atom->tri;

			for (int i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (eflags[i] & SPHERE) {
			omega[i][0] = b->omega[0];
			omega[i][1] = b->omega[1];
			omega[i][2] = b->omega[2];
			} else if (eflags[i] & ELLIPSOID) {
			shape = ebonus[ellipsoid[i]].shape;
			quatatom = ebonus[ellipsoid[i]].quat;
			ione[0] = EINERTIArmass[i] (shape[1]shape[1] + shape[2]shape[2]);
			ione[1] = EINERTIArmass[i] (shape[0]shape[0] + shape[2]shape[2]);
			ione[2] = EINERTIArmass[i] (shape[0]shape[0] + shape[1]shape[1]);
			MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
			MathExtra::omega_to_angmom(b->omega,exone,eyone,ezone,ione,
			angmom[i]);
			} else if (eflags[i] & LINE) {
			omega[i][0] = b->omega[0];
			omega[i][1] = b->omega[1];
			omega[i][2] = b->omega[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(b->omega,exone,eyone,ezone,
			inertiaatom,angmom[i]);
			}
			}
			}
			}

			/* ----------------------------------------------------------------------
			one-time identification of which atoms are in which rigid bodies
			set bodytag for all owned atoms
			------------------------------------------------------------------------- */

			void FixRigidSmall::create_bodies()
			{
			int i,m,n;
			double unwrap[3];

			// error check on image flags of atoms in rigid bodies

			imageint *image = atom->image;
			int *mask = atom->mask;
			int nlocal = atom->nlocal;

			int *periodicity = domain->periodicity;
			int xbox,ybox,zbox;

			int flag = 0;
			for (i = 0; i < nlocal; i++) {
			if (!(mask[i] & groupbit)) continue;
			xbox = (image[i] & IMGMASK) - IMGMAX;
			ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
			zbox = (image[i] >> IMG2BITS) - IMGMAX;
			if ((xbox && !periodicity[0]) \|\| (ybox && !periodicity[1]) \|\|
			(zbox && !periodicity[2])) flag = 1;
			}

			int flagall;
			MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
			if (flagall) error->all(FLERR,"Fix rigid/small atom has non-zero image flag "
			"in a non-periodic dimension");

			// allocate buffer for passing messages around ring of procs
			// percount = max number of values to put in buffer for each of ncount

			int ncount = 0;
			for (i = 0; i < nlocal; i++)
			if (mask[i] & groupbit) ncount++;

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

			// create map hash for storing unique molecule IDs of my atoms
			// key = molecule ID
			// value = index into per-body data structure
			// n = # of entries in hash

			hash = new std::map<tagint,int>();
			hash->clear();

			// setup hash
			// key = body ID
			// 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++;
			}

			// bbox = bounding box of each rigid body my atoms are part of

			memory->create(bbox,n,6,"rigid/small:bbox");

			for (i = 0; i < n; i++) {
			bbox[i][0] = bbox[i][2] = bbox[i][4] = BIG;
			bbox[i][1] = bbox[i][3] = bbox[i][5] = -BIG;
			}

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

			double **x = atom->x;

			m = 0;
			for (i = 0; i < nlocal; i++) {
			if (!(mask[i] & groupbit)) continue;
			domain->unmap(x[i],image[i],unwrap);
			buf[m++] = molecule[i];
			buf[m++] = unwrap[0];
			buf[m++] = unwrap[1];
			buf[m++] = unwrap[2];
			}

			// pass buffer around ring of procs
			// func = update bbox with atom coords from every proc
			// when done, have full bbox for every rigid body my atoms are part of

			frsptr = this;
			comm->ring(m,sizeof(double),buf,1,ring_bbox,NULL);

			// ctr = center pt of each rigid body my atoms are part of

			memory->create(ctr,n,6,"rigid/small:bbox");

			for (i = 0; i < n; i++) {
			ctr[i][0] = 0.5 * (bbox[i][0] + bbox[i][1]);
			ctr[i][1] = 0.5 * (bbox[i][2] + bbox[i][3]);
			ctr[i][2] = 0.5 * (bbox[i][4] + bbox[i][5]);
			}

			// idclose = ID of atom in body closest to center pt (smaller ID if tied)
			// rsqclose = distance squared from idclose to center pt

			memory->create(idclose,n,"rigid/small:idclose");
			memory->create(rsqclose,n,"rigid/small:rsqclose");

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

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

			tagint *tag = atom->tag;

			m = 0;
			for (i = 0; i < nlocal; i++) {
			if (!(mask[i] & groupbit)) continue;
			domain->unmap(x[i],image[i],unwrap);
			buf[m++] = molecule[i];
			buf[m++] = ubuf(tag[i]).d;
			buf[m++] = unwrap[0];
			buf[m++] = unwrap[1];
			buf[m++] = unwrap[2];
			}

			// pass buffer around ring of procs
			// func = update idclose,rsqclose with atom IDs from every proc
			// when done, have idclose for every rigid body my atoms are part of

			frsptr = this;
			comm->ring(m,sizeof(double),buf,2,ring_nearest,NULL);

			// set bodytag of all owned atoms, based on idclose
			// find max value of rsqclose across all procs

			double rsqmax = 0.0;
			for (i = 0; i < nlocal; i++) {
			bodytag[i] = 0;
			if (!(mask[i] & groupbit)) continue;
			m = hash->find(molecule[i])->second;
			bodytag[i] = idclose[m];
			rsqmax = MAX(rsqmax,rsqclose[m]);
			}

			// pack my atoms into buffer as bodytag of owning atom, unwrapped coords

			m = 0;
			for (i = 0; i < nlocal; i++) {
			if (!(mask[i] & groupbit)) continue;
			domain->unmap(x[i],image[i],unwrap);
			buf[m++] = ubuf(bodytag[i]).d;
			buf[m++] = unwrap[0];
			buf[m++] = unwrap[1];
			buf[m++] = unwrap[2];
			}

			// pass buffer around ring of procs
			// func = update rsqfar for atoms belonging to bodies I own
			// when done, have rsqfar for all atoms in bodies I own

			rsqfar = 0.0;
			frsptr = this;
			comm->ring(m,sizeof(double),buf,3,ring_farthest,NULL);

			// find maxextent of rsqfar across all procs
			// if defined, include molecule->maxextent

			MPI_Allreduce(&rsqfar,&maxextent,1,MPI_DOUBLE,MPI_MAX,world);
			maxextent = sqrt(maxextent);
			if (onemol) maxextent = MAX(maxextent,onemol->maxextent);

			// clean up

			delete hash;
			memory->destroy(buf);
			memory->destroy(bbox);
			memory->destroy(ctr);
			memory->destroy(idclose);
			memory->destroy(rsqclose);
			}

			/* ----------------------------------------------------------------------
			process rigid body atoms from another proc
			update bounding box for rigid bodies my atoms are part of
			------------------------------------------------------------------------- */

			void FixRigidSmall::ring_bbox(int n, char *cbuf)
			{
			std::map<tagint,int> *hash = frsptr->hash;
			double **bbox = frsptr->bbox;

			double buf = (double ) cbuf;
			int ndatums = n/4;

			int j,imol;
			double *x;

			int m = 0;
			for (int i = 0; i < ndatums; i++, m += 4) {
			imol = static_cast<int> (buf[m]);
			if (hash->find(imol) != hash->end()) {
			j = hash->find(imol)->second;
			x = &buf[m+1];
			bbox[j][0] = MIN(bbox[j][0],x[0]);
			bbox[j][1] = MAX(bbox[j][1],x[0]);
			bbox[j][2] = MIN(bbox[j][2],x[1]);
			bbox[j][3] = MAX(bbox[j][3],x[1]);
			bbox[j][4] = MIN(bbox[j][4],x[2]);
			bbox[j][5] = MAX(bbox[j][5],x[2]);
			}
			}
			}

			/* ----------------------------------------------------------------------
			process rigid body atoms from another proc
			update nearest atom to body center for rigid bodies my atoms are part of
			------------------------------------------------------------------------- */

			void FixRigidSmall::ring_nearest(int n, char *cbuf)
			{
			std::map<tagint,int> *hash = frsptr->hash;
			double **ctr = frsptr->ctr;
			tagint *idclose = frsptr->idclose;
			double *rsqclose = frsptr->rsqclose;

			double buf = (double ) cbuf;
			int ndatums = n/5;

			int j,imol;
			tagint tag;
			double delx,dely,delz,rsq;
			double *x;

			int m = 0;
			for (int i = 0; i < ndatums; i++, m += 5) {
			imol = static_cast<int> (buf[m]);
			if (hash->find(imol) != hash->end()) {
			j = hash->find(imol)->second;
			tag = (tagint) ubuf(buf[m+1]).i;
			x = &buf[m+2];
			delx = x[0] - ctr[j][0];
			dely = x[1] - ctr[j][1];
			delz = x[2] - ctr[j][2];
			rsq = delxdelx + delydely + delz*delz;
			if (rsq <= rsqclose[j]) {
			if (rsq == rsqclose[j] && tag > idclose[j]) continue;
			idclose[j] = tag;
			rsqclose[j] = rsq;
			}
			}
			}
			}

			/* ----------------------------------------------------------------------
			process rigid body atoms from another proc
			update rsqfar = distance from owning atom to other atom
			------------------------------------------------------------------------- */

			void FixRigidSmall::ring_farthest(int n, char *cbuf)
			{
			double **x = frsptr->atom->x;
			imageint *image = frsptr->atom->image;
			int nlocal = frsptr->atom->nlocal;

			double buf = (double ) cbuf;
			int ndatums = n/4;

			int iowner;
			tagint tag;
			double delx,dely,delz,rsq;
			double *xx;
			double unwrap[3];

			int m = 0;
			for (int i = 0; i < ndatums; i++, m += 4) {
			tag = (tagint) ubuf(buf[m]).i;
			iowner = frsptr->atom->map(tag);
			if (iowner < 0 \|\| iowner >= nlocal) continue;
			frsptr->domain->unmap(x[iowner],image[iowner],unwrap);
			xx = &buf[m+1];
			delx = xx[0] - unwrap[0];
			dely = xx[1] - unwrap[1];
			delz = xx[2] - unwrap[2];
			rsq = delxdelx + delydely + delz*delz;
			frsptr->rsqfar = MAX(frsptr->rsqfar,rsq);
			}
			}

			/* ----------------------------------------------------------------------
			one-time initialization of rigid body attributes
			extended flags, mass, center-of-mass
			Cartesian and diagonalized inertia tensor
			read per-body attributes from infile if specified
			------------------------------------------------------------------------- */

			void FixRigidSmall::setup_bodies_static()
			{
			int i,ibody;

			// extended = 1 if any particle in a rigid body is finite size
			// or has a dipole moment

			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->line_flag \|\|
			atom->tri_flag \|\| atom->mu_flag) {
			int flag = 0;
			for (i = 0; i < nlocal; i++) {
			if (bodytag[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;
			}

			MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
			}

			// extended = 1 if using molecule template with finite-size particles

			if (onemol && onemol->radiusflag) extended = 1;

			// grow extended arrays and set extended flags for each particle
			// 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;
			grow_arrays(atom->nmax);

			for (i = 0; i < nlocal; i++) {
			eflags[i] = 0;
			if (bodytag[i] == 0) continue;

			// set to POINT or SPHERE or ELLIPSOID or LINE

			if (radius && radius[i] > 0.0) {
			eflags[i] \|= SPHERE;
			eflags[i] \|= OMEGA;
			eflags[i] \|= TORQUE;
			} else if (ellipsoid && ellipsoid[i] >= 0) {
			eflags[i] \|= ELLIPSOID;
			eflags[i] \|= ANGMOM;
			eflags[i] \|= TORQUE;
			} 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] \|= DIPOLE;
			}
			}

			// acquire ghost bodies via forward comm
			// set atom2body for ghost atoms via forward comm
			// set atom2body for other owned atoms via reset_atom2body()

			nghost_body = 0;
			commflag = FULL_BODY;
			comm->forward_comm_variable_fix(this);
			reset_atom2body();

			// compute mass & center-of-mass of each rigid body

			double **x = atom->x;
			imageint *image = atom->image;

			double *xcm;

			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			xcm = body[ibody].xcm;
			xcm[0] = xcm[1] = xcm[2] = 0.0;
			body[ibody].mass = 0.0;
			}

			double unwrap[3];
			double massone;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			xcm[0] += unwrap[0] * massone;
			xcm[1] += unwrap[1] * massone;
			xcm[2] += unwrap[2] * massone;
			b->mass += massone;
			}

			// reverse communicate xcm, mass of all bodies

			commflag = XCM_MASS;
			comm->reverse_comm_variable_fix(this);

			for (ibody = 0; ibody < nlocal_body; ibody++) {
			xcm = body[ibody].xcm;
			xcm[0] /= body[ibody].mass;
			xcm[1] /= body[ibody].mass;
			xcm[2] /= body[ibody].mass;
			}

			// overwrite masstotal and center-of-mass with file values
			// inbody[i] = 0/1 if Ith rigid body is initialized by file

			int *inbody;
			if (infile) {
			memory->create(inbody,nlocal_body,"rigid/small:inbody");
			for (ibody = 0; ibody < nlocal_body; ibody++) inbody[ibody] = 0;
			readfile(0,NULL,inbody);
			}

			// set image flags for each rigid body to default values
			// then remap the xcm of each body back into simulation box if needed

			for (ibody = 0; ibody < nlocal_body; ibody++)
			body[ibody].image = ((imageint) IMGMAX << IMG2BITS) \|
			((imageint) IMGMAX << IMGBITS) \| IMGMAX;

			pre_neighbor();

			// compute 6 moments of inertia of each body in Cartesian reference frame
			// dx,dy,dz = coords relative to center-of-mass
			// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector

			memory->create(itensor,nlocal_body+nghost_body,6,"rigid/small:itensor");
			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
			for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;

			double dx,dy,dz;
			double *inertia;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			dx = unwrap[0] - xcm[0];
			dy = unwrap[1] - xcm[1];
			dz = unwrap[2] - xcm[2];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			inertia = itensor[atom2body[i]];
			inertia[0] += massone * (dydy + dzdz);
			inertia[1] += massone * (dxdx + dzdz);
			inertia[2] += massone * (dxdx + dydy);
			inertia[3] -= massone * dy*dz;
			inertia[4] -= massone * dx*dz;
			inertia[5] -= massone * dx*dy;
			}

			// extended particles may contribute extra terms to moments of inertia

			if (extended) {
			double ivec[6];
			double shape,quatatom,*inertiaatom;
			double length,theta;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			inertia = itensor[atom2body[i]];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			if (eflags[i] & SPHERE) {
			inertia[0] += SINERTIAmassone radius[i]*radius[i];
			inertia[1] += SINERTIAmassone radius[i]*radius[i];
			inertia[2] += SINERTIAmassone 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);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[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);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[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);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[5] += ivec[5];
			}
			}
			}

			// reverse communicate inertia tensor of all bodies

			commflag = ITENSOR;
			comm->reverse_comm_variable_fix(this);

			// overwrite Cartesian inertia tensor with file values

			if (infile) readfile(1,itensor,inbody);

			// diagonalize inertia tensor for each body via Jacobi rotations
			// inertia = 3 eigenvalues = principal moments of inertia
			// evectors and exzy_space = 3 evectors = principal axes of rigid body

			int ierror;
			double cross[3];
			double tensor[3][3],evectors[3][3];
			double ex,ey,*ez;

			for (ibody = 0; ibody < nlocal_body; ibody++) {
			tensor[0][0] = itensor[ibody][0];
			tensor[1][1] = itensor[ibody][1];
			tensor[2][2] = itensor[ibody][2];
			tensor[1][2] = tensor[2][1] = itensor[ibody][3];
			tensor[0][2] = tensor[2][0] = itensor[ibody][4];
			tensor[0][1] = tensor[1][0] = itensor[ibody][5];

			inertia = body[ibody].inertia;
			ierror = MathExtra::jacobi(tensor,inertia,evectors);
			if (ierror) error->all(FLERR,
			"Insufficient Jacobi rotations for rigid body");

			ex = body[ibody].ex_space;
			ex[0] = evectors[0][0];
			ex[1] = evectors[1][0];
			ex[2] = evectors[2][0];
			ey = body[ibody].ey_space;
			ey[0] = evectors[0][1];
			ey[1] = evectors[1][1];
			ey[2] = evectors[2][1];
			ez = body[ibody].ez_space;
			ez[0] = evectors[0][2];
			ez[1] = evectors[1][2];
			ez[2] = evectors[2][2];

			// if any principal moment < scaled EPSILON, set to 0.0

			double max;
			max = MAX(inertia[0],inertia[1]);
			max = MAX(max,inertia[2]);

			if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
			if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
			if (inertia[2] < EPSILON*max) inertia[2] = 0.0;

			// enforce 3 evectors as a right-handed coordinate system
			// flip 3rd vector if needed

			MathExtra::cross3(ex,ey,cross);
			if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);

			// create initial quaternion

			MathExtra::exyz_to_q(ex,ey,ez,body[ibody].quat);
			}

			// forward communicate updated info of all bodies

			commflag = INITIAL;
			comm->forward_comm_variable_fix(this);

			// displace = initial atom coords in basis of principal axes
			// set displace = 0.0 for atoms not in any rigid body
			// for extended particles, set their orientation wrt to rigid body

			double qc[4],delta[3];
			double *quatatom;
			double theta_body;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) {
			displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
			continue;
			}

			Body *b = &body[atom2body[i]];

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			delta[0] = unwrap[0] - xcm[0];
			delta[1] = unwrap[1] - xcm[1];
			delta[2] = unwrap[2] - xcm[2];
			MathExtra::transpose_matvec(b->ex_space,b->ey_space,b->ez_space,
			delta,displace[i]);

			if (extended) {
			if (eflags[i] & ELLIPSOID) {
			quatatom = ebonus[ellipsoid[i]].quat;
			MathExtra::qconjugate(b->quat,qc);
			MathExtra::quatquat(qc,quatatom,orient[i]);
			MathExtra::qnormalize(orient[i]);
			} else if (eflags[i] & LINE) {
			if (b->quat[3] >= 0.0) 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;
			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(b->quat,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(b->ex_space,b->ey_space,b->ez_space,
			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;
			}
			}

			// test for valid principal moments & axes
			// recompute moments of inertia around new axes
			// 3 diagonal moments should equal principal moments
			// 3 off-diagonal moments should be 0.0
			// extended particles may contribute extra terms to moments of inertia

			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
			for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			inertia = itensor[atom2body[i]];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			inertia[0] += massone *
			(displace[i][1]displace[i][1] + displace[i][2]displace[i][2]);
			inertia[1] += massone *
			(displace[i][0]displace[i][0] + displace[i][2]displace[i][2]);
			inertia[2] += massone *
			(displace[i][0]displace[i][0] + displace[i][1]displace[i][1]);
			inertia[3] -= massone * displace[i][1]*displace[i][2];
			inertia[4] -= massone * displace[i][0]*displace[i][2];
			inertia[5] -= massone * displace[i][0]*displace[i][1];
			}

			if (extended) {
			double ivec[6];
			double shape,inertiaatom;
			double length;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			inertia = itensor[atom2body[i]];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			if (eflags[i] & SPHERE) {
			inertia[0] += SINERTIAmassone radius[i]*radius[i];
			inertia[1] += SINERTIAmassone radius[i]*radius[i];
			inertia[2] += SINERTIAmassone radius[i]*radius[i];
			} else if (eflags[i] & ELLIPSOID) {
			shape = ebonus[ellipsoid[i]].shape;
			MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[5] += ivec[5];
			} else if (eflags[i] & LINE) {
			length = lbonus[line[i]].length;
			MathExtra::inertia_line(length,orient[i][0],massone,ivec);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[5] += ivec[5];
			} else if (eflags[i] & TRIANGLE) {
			inertiaatom = tbonus[tri[i]].inertia;
			MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec);
			inertia[0] += ivec[0];
			inertia[1] += ivec[1];
			inertia[2] += ivec[2];
			inertia[3] += ivec[3];
			inertia[4] += ivec[4];
			inertia[5] += ivec[5];
			}
			}
			}

			// reverse communicate inertia tensor of all bodies

			commflag = ITENSOR;
			comm->reverse_comm_variable_fix(this);

			// error check that re-computed momemts of inertia match diagonalized ones
			// do not do test for bodies with params read from infile

			double norm;
			for (ibody = 0; ibody < nlocal_body; ibody++) {
			if (infile && inbody[ibody]) continue;
			inertia = body[ibody].inertia;

			if (inertia[0] == 0.0) {
			if (fabs(itensor[ibody][0]) > TOLERANCE)
			error->all(FLERR,"Fix rigid: Bad principal moments");
			} else {
			if (fabs((itensor[ibody][0]-inertia[0])/inertia[0]) >
			TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
			}
			if (inertia[1] == 0.0) {
			if (fabs(itensor[ibody][1]) > TOLERANCE)
			error->all(FLERR,"Fix rigid: Bad principal moments");
			} else {
			if (fabs((itensor[ibody][1]-inertia[1])/inertia[1]) >
			TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
			}
			if (inertia[2] == 0.0) {
			if (fabs(itensor[ibody][2]) > TOLERANCE)
			error->all(FLERR,"Fix rigid: Bad principal moments");
			} else {
			if (fabs((itensor[ibody][2]-inertia[2])/inertia[2]) >
			TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
			}
			norm = (inertia[0] + inertia[1] + inertia[2]) / 3.0;
			if (fabs(itensor[ibody][3]/norm) > TOLERANCE \|\|
			fabs(itensor[ibody][4]/norm) > TOLERANCE \|\|
			fabs(itensor[ibody][5]/norm) > TOLERANCE)
			error->all(FLERR,"Fix rigid: Bad principal moments");
			}

			// clean up

			memory->destroy(itensor);
			if (infile) memory->destroy(inbody);
			}

			/* ----------------------------------------------------------------------
			one-time initialization of dynamic rigid body attributes
			Vcm and angmom, computed explicitly from constituent particles
			even if wrong for overlapping particles, is OK,
			since is just setting initial time=0 Vcm and angmom of the body
			which can be estimated value
			------------------------------------------------------------------------- */

			void FixRigidSmall::setup_bodies_dynamic()
			{
			int i,ibody;
			double massone,radone;

			// sum vcm, angmom across all rigid bodies
			// vcm = velocity of COM
			// angmom = angular momentum around COM

			double **x = atom->x;
			double **v = atom->v;
			double *rmass = atom->rmass;
			double *mass = atom->mass;
			int *type = atom->type;
			imageint *image = atom->image;
			int nlocal = atom->nlocal;

			double xcm,vcm,*acm;
			double dx,dy,dz;
			double unwrap[3];

			for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			vcm = body[ibody].vcm;
			vcm[0] = vcm[1] = vcm[2] = 0.0;
			acm = body[ibody].angmom;
			acm[0] = acm[1] = acm[2] = 0.0;
			}

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (rmass) massone = rmass[i];
			else massone = mass[type[i]];

			vcm = b->vcm;
			vcm[0] += v[i][0] * massone;
			vcm[1] += v[i][1] * massone;
			vcm[2] += v[i][2] * massone;

			domain->unmap(x[i],image[i],unwrap);
			xcm = b->xcm;
			dx = unwrap[0] - xcm[0];
			dy = unwrap[1] - xcm[1];
			dz = unwrap[2] - xcm[2];

			acm = b->angmom;
			acm[0] += dy * massonev[i][2] - dz massone*v[i][1];
			acm[1] += dz * massonev[i][0] - dx massone*v[i][2];
			acm[2] += dx * massonev[i][1] - dy massone*v[i][0];
			}

			// extended particles add their rotation to angmom of body

			if (extended) {
			AtomVecLine::Bonus *lbonus;
			if (avec_line) lbonus = avec_line->bonus;
			double **omega = atom->omega;
			double **angmom = atom->angmom;
			double *radius = atom->radius;
			int *line = atom->line;

			for (i = 0; i < nlocal; i++) {
			if (atom2body[i] < 0) continue;
			Body *b = &body[atom2body[i]];

			if (eflags[i] & OMEGA) {
			if (eflags[i] & SPHERE) {
			radone = radius[i];
			acm = b->angmom;
			acm[0] += SINERTIArmass[i] radoneradone omega[i][0];
			acm[1] += SINERTIArmass[i] radoneradone omega[i][1];
			acm[2] += SINERTIArmass[i] radoneradone omega[i][2];
			} else if (eflags[i] & LINE) {
			radone = lbonus[line[i]].length;
			b->angmom[2] += LINERTIArmass[i] radoneradone omega[i][2];
			}
			}
			if (eflags[i] & ANGMOM) {
			acm = b->angmom;
			acm[0] += angmom[i][0];
			acm[1] += angmom[i][1];
			acm[2] += angmom[i][2];
			}
			}
			}

			// reverse communicate vcm, angmom of all bodies

			commflag = VCM_ANGMOM;
			comm->reverse_comm_variable_fix(this);

			// normalize velocity of COM

			for (ibody = 0; ibody < nlocal_body; ibody++) {
			vcm = body[ibody].vcm;
			vcm[0] /= body[ibody].mass;
			vcm[1] /= body[ibody].mass;
			vcm[2] /= body[ibody].mass;
			}
			}

			/* ----------------------------------------------------------------------
			read per rigid body info from user-provided file
			which = 0 to read total mass and center-of-mass
			which = 1 to read 6 moments of inertia, store in array
			flag inbody = 0 for local bodies whose info is read from file
			nlines = # of lines of rigid body info, 0 is OK
			one line = rigid-ID mass xcm ycm zcm ixx iyy izz ixy ixz iyz
			and rigid-ID = mol-ID for fix rigid/small
			------------------------------------------------------------------------- */

			void FixRigidSmall::readfile(int which, double *array, int inbody)
			{
			int i,j,m,nchunk,eofflag,nlines;
			tagint id;
			FILE *fp;
			char eof,start,next,buf;
			char line[MAXLINE];

			// create local hash with key/value pairs
			// key = mol ID of bodies my atoms own
			// value = index into local body array

			int nlocal = atom->nlocal;

			hash = new std::map<tagint,int>();
			for (i = 0; i < nlocal; i++)
			if (bodyown[i] >= 0) (*hash)[atom->molecule[i]] = bodyown[i];

			// open file and read header

			if (me == 0) {
			fp = fopen(infile,"r");
			if (fp == NULL) {
			char str[128];
			sprintf(str,"Cannot open fix rigid/small infile %s",infile);
			error->one(FLERR,str);
			}

			while (1) {
			eof = fgets(line,MAXLINE,fp);
			if (eof == NULL)
			error->one(FLERR,"Unexpected end of fix rigid/small file");
			start = &line[strspn(line," \t\n\v\f\r")];
			if (start != '\0' && start != '#') break;
			}

			sscanf(line,"%d",&nlines);
			}

			MPI_Bcast(&nlines,1,MPI_INT,0,world);

			char buffer = new char[CHUNKMAXLINE];
			char *values = new char[ATTRIBUTE_PERBODY];

			int nread = 0;
			while (nread < nlines) {
			nchunk = MIN(nlines-nread,CHUNK);
			eofflag = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
			if (eofflag) error->all(FLERR,"Unexpected end of fix rigid/small file");

			buf = buffer;
			next = strchr(buf,'\n');
			*next = '\0';
			int nwords = atom->count_words(buf);
			*next = '\n';

			if (nwords != ATTRIBUTE_PERBODY)
			error->all(FLERR,"Incorrect rigid body format in fix rigid/small file");

			// loop over lines of rigid body attributes
			// tokenize the line into values
			// id = rigid body ID = mol-ID
			// for which = 0, store mass/com in vec/array
			// for which = 1, store interia tensor array, invert 3,4,5 values to Voigt

			for (int i = 0; i < nchunk; i++) {
			next = strchr(buf,'\n');

			values[0] = strtok(buf," \t\n\r\f");
			for (j = 1; j < nwords; j++)
			values[j] = strtok(NULL," \t\n\r\f");

			id = ATOTAGINT(values[0]);
			if (id <= 0 \|\| id > maxmol)
			error->all(FLERR,"Invalid rigid body ID in fix rigid/small file");
			if (hash->find(id) == hash->end()) {
			buf = next + 1;
			continue;
			}
			m = (*hash)[id];
			inbody[m] = 1;

			if (which == 0) {
			body[m].mass = atof(values[1]);
			body[m].xcm[0] = atof(values[2]);
			body[m].xcm[1] = atof(values[3]);
			body[m].xcm[2] = atof(values[4]);
			} else {
			array[m][0] = atof(values[5]);
			array[m][1] = atof(values[6]);
			array[m][2] = atof(values[7]);
			array[m][3] = atof(values[10]);
			array[m][4] = atof(values[9]);
			array[m][5] = atof(values[8]);
			}

			buf = next + 1;
			}

			nread += nchunk;
			}

			if (me == 0) fclose(fp);

			delete [] buffer;
			delete [] values;
			delete hash;
			}

			/* ----------------------------------------------------------------------
			write out restart info for mass, COM, inertia tensor to file
			identical format to infile option, so info can be read in when restarting
			each proc contributes info for rigid bodies it owns
			------------------------------------------------------------------------- */

			void FixRigidSmall::write_restart_file(char *file)
			{
			FILE *fp;

			// do not write file if bodies have not yet been intialized

			if (firstflag) return;

			// proc 0 opens file and writes header

			if (me == 0) {
			char outfile[128];
			sprintf(outfile,"%s.rigid",file);
			fp = fopen(outfile,"w");
			if (fp == NULL) {
			char str[128];
			sprintf(str,"Cannot open fix rigid restart file %s",outfile);
			error->one(FLERR,str);
			}

			fprintf(fp,"# fix rigid mass, COM, inertia tensor info for "
			"%d bodies on timestep " BIGINT_FORMAT "\n\n",
			nbody,update->ntimestep);
			fprintf(fp,"%d\n",nbody);
			}

			// communication buffer for all my rigid body info
			// max_size = largest buffer needed by any proc

			int ncol = 11;
			int sendrow = nlocal_body;
			int maxrow;
			MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);

			double **buf;
			if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"rigid/small:buf");
			else memory->create(buf,MAX(1,sendrow),ncol,"rigid/small:buf");

			// pack my rigid body info into buf
			// compute I tensor against xyz axes from diagonalized I and current quat
			// Ispace = P Idiag P_transpose
			// P is stored column-wise in exyz_space

			double p[3][3],pdiag[3][3],ispace[3][3];

			for (int i = 0; i < nlocal_body; i++) {
			MathExtra::col2mat(body[i].ex_space,body[i].ey_space,body[i].ez_space,p);
			MathExtra::times3_diag(p,body[i].inertia,pdiag);
			MathExtra::times3_transpose(pdiag,p,ispace);

			buf[i][0] = atom->molecule[body[i].ilocal];
			buf[i][1] = body[i].mass;
			buf[i][2] = body[i].xcm[0];
			buf[i][3] = body[i].xcm[1];
			buf[i][4] = body[i].xcm[2];
			buf[i][5] = ispace[0][0];
			buf[i][6] = ispace[1][1];
			buf[i][7] = ispace[2][2];
			buf[i][8] = ispace[0][1];
			buf[i][9] = ispace[0][2];
			buf[i][10] = ispace[1][2];
			}

			// write one chunk of rigid body info per proc to file
			// proc 0 pings each proc, receives its chunk, writes to file
			// all other procs wait for ping, send their chunk to proc 0

			int tmp,recvrow;
			MPI_Status status;
			MPI_Request request;

			if (me == 0) {
			for (int iproc = 0; iproc < nprocs; iproc++) {
			if (iproc) {
			MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
			MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
			MPI_Wait(&request,&status);
			MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
			recvrow /= ncol;
			} else recvrow = sendrow;

			for (int i = 0; i < recvrow; i++)
			fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e "
			"%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n",
			static_cast<int> (buf[i][0]),buf[i][1],
			buf[i][2],buf[i][3],buf[i][4],
			buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],buf[i][10]);
			}

			} else {
			MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
			MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
			}

			// clean up and close file

			memory->destroy(buf);
			if (me == 0) fclose(fp);
			}

			/* ----------------------------------------------------------------------
			allocate local atom-based arrays
			------------------------------------------------------------------------- */

			void FixRigidSmall::grow_arrays(int nmax)
			{
			memory->grow(bodyown,nmax,"rigid/small:bodyown");
			memory->grow(bodytag,nmax,"rigid/small:bodytag");
			memory->grow(atom2body,nmax,"rigid/small:atom2body");
			memory->grow(displace,nmax,3,"rigid/small:displace");
			if (extended) {
			memory->grow(eflags,nmax,"rigid/small:eflags");
			if (orientflag) memory->grow(orient,nmax,orientflag,"rigid/small:orient");
			if (dorientflag) memory->grow(dorient,nmax,3,"rigid/small:dorient");
			}
			}

			/* ----------------------------------------------------------------------
			copy values within local atom-based arrays
			------------------------------------------------------------------------- */

			void FixRigidSmall::copy_arrays(int i, int j, int delflag)
			{
			bodytag[j] = bodytag[i];
			displace[j][0] = displace[i][0];
			displace[j][1] = displace[i][1];
			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 deleting atom J via delflag and J owns a body, then delete it

			if (delflag && bodyown[j] >= 0) {
			bodyown[body[nlocal_body-1].ilocal] = bodyown[j];
			memcpy(&body[bodyown[j]],&body[nlocal_body-1],sizeof(Body));
			nlocal_body--;
			}

			// if atom I owns a body, reset I's body.ilocal to loc J
			// do NOT do this if self-copy (I=J) since I's body is already deleted

			if (bodyown[i] >= 0 && i != j) body[bodyown[i]].ilocal = j;
			bodyown[j] = bodyown[i];
			}

			/* ----------------------------------------------------------------------
			initialize one atom's array values, called when atom is created
			------------------------------------------------------------------------- */

			void FixRigidSmall::set_arrays(int i)
			{
			bodyown[i] = -1;
			bodytag[i] = 0;
			atom2body[i] = -1;
			displace[i][0] = 0.0;
			displace[i][1] = 0.0;
			displace[i][2] = 0.0;
			}

			/* ----------------------------------------------------------------------
			initialize a molecule inserted by another fix, e.g. deposit or pour
			called when molecule is created
			nlocalprev = # of atoms on this proc before molecule inserted
			tagprev = atom ID previous to new atoms in the molecule
			xgeom = geometric center of new molecule
			vcm = COM velocity of new molecule
			quat = rotation of new molecule (around geometric center)
			relative to template in Molecule class
			------------------------------------------------------------------------- */

			void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev,
			double xgeom, double vcm, double *quat)
			{
			int m;
			double ctr2com[3],ctr2com_rotate[3];
			double rotmat[3][3];

			int nlocal = atom->nlocal;
			if (nlocalprev == nlocal) return;

			tagint *tag = atom->tag;

			for (int i = nlocalprev; i < nlocal; i++) {
			bodytag[i] = tagprev + onemol->comatom;
			if (tag[i]-tagprev == onemol->comatom) bodyown[i] = nlocal_body;

			m = tag[i] - tagprev-1;
			displace[i][0] = onemol->dxbody[m][0];
			displace[i][1] = onemol->dxbody[m][1];
			displace[i][2] = onemol->dxbody[m][2];

			eflags[i] = 0;
			if (onemol->radiusflag) {
			eflags[i] \|= SPHERE;
			eflags[i] \|= OMEGA;
			eflags[i] \|= TORQUE;
			}

			if (bodyown[i] >= 0) {
			if (nlocal_body == nmax_body) grow_body();
			Body *b = &body[nlocal_body];
			b->mass = onemol->masstotal;

			// new COM = Q (onemol->xcm - onemol->center) + xgeom
			// Q = rotation matrix associated with quat

			MathExtra::quat_to_mat(quat,rotmat);
			MathExtra::sub3(onemol->com,onemol->center,ctr2com);
			MathExtra::matvec(rotmat,ctr2com,ctr2com_rotate);
			MathExtra::add3(ctr2com_rotate,xgeom,b->xcm);

			b->vcm[0] = vcm[0];
			b->vcm[1] = vcm[1];
			b->vcm[2] = vcm[2];
			b->inertia[0] = onemol->inertia[0];
			b->inertia[1] = onemol->inertia[1];
			b->inertia[2] = onemol->inertia[2];

			// final quat is product of insertion quat and original quat
			// true even if insertion rotation was not around COM

			MathExtra::quatquat(quat,onemol->quat,b->quat);
			MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);

			b->angmom[0] = b->angmom[1] = b->angmom[2] = 0.0;
			b->omega[0] = b->omega[1] = b->omega[2] = 0.0;
			b->image = ((imageint) IMGMAX << IMG2BITS) \|
			((imageint) IMGMAX << IMGBITS) \| IMGMAX;
			b->ilocal = i;
			nlocal_body++;
			}
			}

			// increment total # of rigid bodies

			nbody++;
			}

			/* ----------------------------------------------------------------------
			pack values in local atom-based arrays for exchange with another proc
			------------------------------------------------------------------------- */

			int FixRigidSmall::pack_exchange(int i, double *buf)
			{
			buf[0] = ubuf(bodytag[i]).d;
			buf[1] = displace[i][0];
			buf[2] = displace[i][1];
			buf[3] = displace[i][2];

			// extended attribute info

			int m = 4;
			if (extended) {
			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];
			}
			}

			// atom not in a rigid body

			if (!bodytag[i]) return m;

			// atom does not own its rigid body

			if (bodyown[i] < 0) {
			buf[m++] = 0;
			return m;
			}

			// body info for atom that owns a rigid body

			buf[m++] = 1;
			memcpy(&buf[m],&body[bodyown[i]],sizeof(Body));
			m += bodysize;
			return m;
			}

			/* ----------------------------------------------------------------------
			unpack values in local atom-based arrays from exchange with another proc
			------------------------------------------------------------------------- */

			int FixRigidSmall::unpack_exchange(int nlocal, double *buf)
			{
			bodytag[nlocal] = (tagint) ubuf(buf[0]).i;
			displace[nlocal][0] = buf[1];
			displace[nlocal][1] = buf[2];
			displace[nlocal][2] = buf[3];

			// extended attribute info

			int m = 4;
			if (extended) {
			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++];
			}
			}

			// atom not in a rigid body

			if (!bodytag[nlocal]) {
			bodyown[nlocal] = -1;
			return m;
			}

			// atom does not own its rigid body

			bodyown[nlocal] = static_cast<int> (buf[m++]);
			if (bodyown[nlocal] == 0) {
			bodyown[nlocal] = -1;
			return m;
			}

			// body info for atom that owns a rigid body

			if (nlocal_body == nmax_body) grow_body();
			memcpy(&body[nlocal_body],&buf[m],sizeof(Body));
			m += bodysize;
			body[nlocal_body].ilocal = nlocal;
			bodyown[nlocal] = nlocal_body++;

			return m;
			}

			/* ----------------------------------------------------------------------
			only pack body info if own or ghost atom owns the body
			for FULL_BODY, send 0/1 flag with every atom
			------------------------------------------------------------------------- */

			int FixRigidSmall::pack_comm(int n, int list, double buf,
			int pbc_flag, int *pbc)
			{
			int i,j;
			double xcm,vcm,quat,omega,ex_space,ey_space,ez_space,conjqm;

			int m = 0;

			if (commflag == INITIAL) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			xcm = body[bodyown[j]].xcm;
			buf[m++] = xcm[0];
			buf[m++] = xcm[1];
			buf[m++] = xcm[2];
			vcm = body[bodyown[j]].vcm;
			buf[m++] = vcm[0];
			buf[m++] = vcm[1];
			buf[m++] = vcm[2];
			quat = body[bodyown[j]].quat;
			buf[m++] = quat[0];
			buf[m++] = quat[1];
			buf[m++] = quat[2];
			buf[m++] = quat[3];
			omega = body[bodyown[j]].omega;
			buf[m++] = omega[0];
			buf[m++] = omega[1];
			buf[m++] = omega[2];
			ex_space = body[bodyown[j]].ex_space;
			buf[m++] = ex_space[0];
			buf[m++] = ex_space[1];
			buf[m++] = ex_space[2];
			ey_space = body[bodyown[j]].ey_space;
			buf[m++] = ey_space[0];
			buf[m++] = ey_space[1];
			buf[m++] = ey_space[2];
			ez_space = body[bodyown[j]].ez_space;
			buf[m++] = ez_space[0];
			buf[m++] = ez_space[1];
			buf[m++] = ez_space[2];
			conjqm = body[bodyown[j]].conjqm;
			buf[m++] = conjqm[0];
			buf[m++] = conjqm[1];
			buf[m++] = conjqm[2];
			buf[m++] = conjqm[3];
			}

			} else if (commflag == FINAL) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			vcm = body[bodyown[j]].vcm;
			buf[m++] = vcm[0];
			buf[m++] = vcm[1];
			buf[m++] = vcm[2];
			omega = body[bodyown[j]].omega;
			buf[m++] = omega[0];
			buf[m++] = omega[1];
			buf[m++] = omega[2];
			conjqm = body[bodyown[j]].conjqm;
			buf[m++] = conjqm[0];
			buf[m++] = conjqm[1];
			buf[m++] = conjqm[2];
			buf[m++] = conjqm[3];
			}

			} else if (commflag == FULL_BODY) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) buf[m++] = 0;
			else {
			buf[m++] = 1;
			memcpy(&buf[m],&body[bodyown[j]],sizeof(Body));
			m += bodysize;
			}
			}
			}

			return m;
			}

			/* ----------------------------------------------------------------------
			only ghost atoms are looped over
			for FULL_BODY, store a new ghost body if this atom owns it
			for other commflag values, only unpack body info if atom owns it
			------------------------------------------------------------------------- */

			void FixRigidSmall::unpack_comm(int n, int first, double *buf)
			{
			int i,j,last;
			double xcm,vcm,quat,omega,ex_space,ey_space,ez_space,conjqm;

			int m = 0;
			last = first + n;

			if (commflag == INITIAL) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			xcm = body[bodyown[i]].xcm;
			xcm[0] = buf[m++];
			xcm[1] = buf[m++];
			xcm[2] = buf[m++];
			vcm = body[bodyown[i]].vcm;
			vcm[0] = buf[m++];
			vcm[1] = buf[m++];
			vcm[2] = buf[m++];
			quat = body[bodyown[i]].quat;
			quat[0] = buf[m++];
			quat[1] = buf[m++];
			quat[2] = buf[m++];
			quat[3] = buf[m++];
			omega = body[bodyown[i]].omega;
			omega[0] = buf[m++];
			omega[1] = buf[m++];
			omega[2] = buf[m++];
			ex_space = body[bodyown[i]].ex_space;
			ex_space[0] = buf[m++];
			ex_space[1] = buf[m++];
			ex_space[2] = buf[m++];
			ey_space = body[bodyown[i]].ey_space;
			ey_space[0] = buf[m++];
			ey_space[1] = buf[m++];
			ey_space[2] = buf[m++];
			ez_space = body[bodyown[i]].ez_space;
			ez_space[0] = buf[m++];
			ez_space[1] = buf[m++];
			ez_space[2] = buf[m++];
			conjqm = body[bodyown[i]].conjqm;
			conjqm[0] = buf[m++];
			conjqm[1] = buf[m++];
			conjqm[2] = buf[m++];
			conjqm[3] = buf[m++];
			}

			} else if (commflag == FINAL) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			vcm = body[bodyown[i]].vcm;
			vcm[0] = buf[m++];
			vcm[1] = buf[m++];
			vcm[2] = buf[m++];
			omega = body[bodyown[i]].omega;
			omega[0] = buf[m++];
			omega[1] = buf[m++];
			omega[2] = buf[m++];
			conjqm = body[bodyown[i]].conjqm;
			conjqm[0] = buf[m++];
			conjqm[1] = buf[m++];
			conjqm[2] = buf[m++];
			conjqm[3] = buf[m++];
			}

			} else if (commflag == FULL_BODY) {
			for (i = first; i < last; i++) {
			bodyown[i] = static_cast<int> (buf[m++]);
			if (bodyown[i] == 0) bodyown[i] = -1;
			else {
			j = nlocal_body + nghost_body;
			if (j == nmax_body) grow_body();
			memcpy(&body[j],&buf[m],sizeof(Body));
			m += bodysize;
			body[j].ilocal = i;
			bodyown[i] = j;
			nghost_body++;
			}
			}
			}
			}

			/* ----------------------------------------------------------------------
			only ghost atoms are looped over
			only pack body info if atom owns it
			------------------------------------------------------------------------- */

			int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
			{
			int i,j,m,last;
			double fcm,torque,vcm,angmom,*xcm;

			m = 0;
			last = first + n;

			if (commflag == FORCE_TORQUE) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			fcm = body[bodyown[i]].fcm;
			buf[m++] = fcm[0];
			buf[m++] = fcm[1];
			buf[m++] = fcm[2];
			torque = body[bodyown[i]].torque;
			buf[m++] = torque[0];
			buf[m++] = torque[1];
			buf[m++] = torque[2];
			}

			} else if (commflag == VCM_ANGMOM) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			vcm = body[bodyown[i]].vcm;
			buf[m++] = vcm[0];
			buf[m++] = vcm[1];
			buf[m++] = vcm[2];
			angmom = body[bodyown[i]].angmom;
			buf[m++] = angmom[0];
			buf[m++] = angmom[1];
			buf[m++] = angmom[2];
			}

			} else if (commflag == XCM_MASS) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			xcm = body[bodyown[i]].xcm;
			buf[m++] = xcm[0];
			buf[m++] = xcm[1];
			buf[m++] = xcm[2];
			buf[m++] = body[bodyown[i]].mass;
			}

			} else if (commflag == ITENSOR) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			j = bodyown[i];
			buf[m++] = itensor[j][0];
			buf[m++] = itensor[j][1];
			buf[m++] = itensor[j][2];
			buf[m++] = itensor[j][3];
			buf[m++] = itensor[j][4];
			buf[m++] = itensor[j][5];
			}

			} else if (commflag == DOF) {
			for (i = first; i < last; i++) {
			if (bodyown[i] < 0) continue;
			j = bodyown[i];
			buf[m++] = counts[j][0];
			buf[m++] = counts[j][1];
			buf[m++] = counts[j][2];
			}
			}

			return m;
			}

			/* ----------------------------------------------------------------------
			only unpack body info if own or ghost atom owns the body
			------------------------------------------------------------------------- */

			void FixRigidSmall::unpack_reverse_comm(int n, int list, double buf)
			{
			int i,j,k;
			double fcm,torque,vcm,angmom,*xcm;

			int m = 0;

			if (commflag == FORCE_TORQUE) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			fcm = body[bodyown[j]].fcm;
			fcm[0] += buf[m++];
			fcm[1] += buf[m++];
			fcm[2] += buf[m++];
			torque = body[bodyown[j]].torque;
			torque[0] += buf[m++];
			torque[1] += buf[m++];
			torque[2] += buf[m++];
			}

			} else if (commflag == VCM_ANGMOM) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			vcm = body[bodyown[j]].vcm;
			vcm[0] += buf[m++];
			vcm[1] += buf[m++];
			vcm[2] += buf[m++];
			angmom = body[bodyown[j]].angmom;
			angmom[0] += buf[m++];
			angmom[1] += buf[m++];
			angmom[2] += buf[m++];
			}

			} else if (commflag == XCM_MASS) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			xcm = body[bodyown[j]].xcm;
			xcm[0] += buf[m++];
			xcm[1] += buf[m++];
			xcm[2] += buf[m++];
			body[bodyown[j]].mass += buf[m++];
			}

			} else if (commflag == ITENSOR) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			k = bodyown[j];
			itensor[k][0] += buf[m++];
			itensor[k][1] += buf[m++];
			itensor[k][2] += buf[m++];
			itensor[k][3] += buf[m++];
			itensor[k][4] += buf[m++];
			itensor[k][5] += buf[m++];
			}

			} else if (commflag == DOF) {
			for (i = 0; i < n; i++) {
			j = list[i];
			if (bodyown[j] < 0) continue;
			k = bodyown[j];
			counts[k][0] += static_cast<int> (buf[m++]);
			counts[k][1] += static_cast<int> (buf[m++]);
			counts[k][2] += static_cast<int> (buf[m++]);
			}
			}
			}

			/* ----------------------------------------------------------------------
			grow body data structure
			------------------------------------------------------------------------- */

			void FixRigidSmall::grow_body()
			{
			nmax_body += DELTA_BODY;
			body = (Body ) memory->srealloc(body,nmax_bodysizeof(Body),
			"rigid/small:body");
			}

			/* ----------------------------------------------------------------------
			reset atom2body for all owned atoms
			do this via bodyown of atom that owns the body the owned atom is in
			atom2body values can point to original body or any image of the body
			------------------------------------------------------------------------- */

			void FixRigidSmall::reset_atom2body()
			{
			int iowner;

			// iowner = index of atom that owns the body that atom I is in

			int nlocal = atom->nlocal;

			for (int i = 0; i < nlocal; i++) {
			atom2body[i] = -1;
			if (bodytag[i]) {
			iowner = atom->map(bodytag[i]);
			if (iowner == -1) {
			char str[128];
			sprintf(str,
			"Rigid body atoms " TAGINT_FORMAT " " TAGINT_FORMAT
			" missing on proc %d at step " BIGINT_FORMAT,
			atom->tag[i],bodytag[i],comm->me,update->ntimestep);
			error->one(FLERR,str);

			}
			atom2body[i] = bodyown[iowner];
			}
			}
			}

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

			void FixRigidSmall::reset_dt()
			{
			dtv = update->dt;
			dtf = 0.5 * update->dt * force->ftm2v;
			dtq = 0.5 * update->dt;
			}

			/* ----------------------------------------------------------------------
			zero linear momentum of each rigid body
			set Vcm to 0.0, then reset velocities of particles via set_v()
			------------------------------------------------------------------------- */

			void FixRigidSmall::zero_momentum()
			{
			double *vcm;
			for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			vcm = body[ibody].vcm;
			vcm[0] = vcm[1] = vcm[2] = 0.0;
			}

			// forward communicate of vcm to all ghost copies

			commflag = FINAL;
			comm->forward_comm_variable_fix(this);

			// set velocity of atoms in rigid bodues

			evflag = 0;
			set_v();
			}

			/* ----------------------------------------------------------------------
			zero angular momentum of each rigid body
			set angmom/omega to 0.0, then reset velocities of particles via set_v()
			------------------------------------------------------------------------- */

			void FixRigidSmall::zero_rotation()
			{
			double angmom,omega;
			for (int ibody = 0; ibody < nlocal_body+nghost_body; ibody++) {
			angmom = body[ibody].angmom;
			angmom[0] = angmom[1] = angmom[2] = 0.0;
			omega = body[ibody].omega;
			omega[0] = omega[1] = omega[2] = 0.0;
			}

			// forward communicate of omega to all ghost copies

			commflag = FINAL;
			comm->forward_comm_variable_fix(this);

			// set velocity of atoms in rigid bodues

			evflag = 0;
			set_v();
			}

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

			void FixRigidSmall::extract(const char str, int &dim)
			{
			if (strcmp(str,"body") == 0) {
			dim = 1;
			return atom2body;
			}

			if (strcmp(str,"onemol") == 0) {
			dim = 0;
			return onemol;
			}

			// return vector of rigid body masses, for owned+ghost bodies
			// used by granular pair styles, indexed by atom2body

			if (strcmp(str,"masstotal") == 0) {
			dim = 1;

			if (nmax_mass < nmax_body) {
			memory->destroy(mass_body);
			nmax_mass = nmax_body;
			memory->create(mass_body,nmax_mass,"rigid:mass_body");
			}

			int n = nlocal_body + nghost_body;
			for (int i = 0; i < n; i++)
			mass_body[i] = body[i].mass;

			return mass_body;
			}

			return NULL;
			}

			/* ----------------------------------------------------------------------
			return translational KE for all rigid bodies
			KE = 1/2 M Vcm^2
			sum local body results across procs
			------------------------------------------------------------------------- */

			double FixRigidSmall::extract_ke()
			{
			double *vcm;

			double ke = 0.0;
			for (int i = 0; i < nlocal_body; i++) {
			vcm = body[i].vcm;
			ke += body[i].mass * (vcm[0]vcm[0] + vcm[1]vcm[1] + vcm[2]*vcm[2]);
			}

			double keall;
			MPI_Allreduce(&ke,&keall,1,MPI_DOUBLE,MPI_SUM,world);

			return 0.5*keall;
			}

			/* ----------------------------------------------------------------------
			return rotational KE for all rigid bodies
			Erotational = 1/2 I wbody^2
			------------------------------------------------------------------------- */

			double FixRigidSmall::extract_erotational()
			{
			double wbody[3],rot[3][3];
			double *inertia;

			double erotate = 0.0;
			for (int i = 0; i < nlocal_body; i++) {

			// for Iw^2 rotational term, need wbody = angular velocity in body frame
			// not omega = angular velocity in space frame

			inertia = body[i].inertia;
			MathExtra::quat_to_mat(body[i].quat,rot);
			MathExtra::transpose_matvec(rot,body[i].angmom,wbody);
			if (inertia[0] == 0.0) wbody[0] = 0.0;
			else wbody[0] /= inertia[0];
			if (inertia[1] == 0.0) wbody[1] = 0.0;
			else wbody[1] /= inertia[1];
			if (inertia[2] == 0.0) wbody[2] = 0.0;
			else wbody[2] /= inertia[2];

			erotate += inertia[0]wbody[0]wbody[0] + inertia[1]wbody[1]wbody[1] +
			inertia[2]wbody[2]wbody[2];
			}

			double erotateall;
			MPI_Allreduce(&erotate,&erotateall,1,MPI_DOUBLE,MPI_SUM,world);

			return 0.5*erotateall;
			}

			/* ----------------------------------------------------------------------
			return temperature of collection of rigid bodies
			non-active DOF are removed by fflag/tflag and in tfactor
			------------------------------------------------------------------------- */

			double FixRigidSmall::compute_scalar()
			{
			double wbody[3],rot[3][3];

			double vcm,inertia;

			double t = 0.0;

			for (int i = 0; i < nlocal_body; i++) {
			vcm = body[i].vcm;
			t += body[i].mass * (vcm[0]vcm[0] + vcm[1]vcm[1] + vcm[2]*vcm[2]);

			// for Iw^2 rotational term, need wbody = angular velocity in body frame
			// not omega = angular velocity in space frame

			inertia = body[i].inertia;
			MathExtra::quat_to_mat(body[i].quat,rot);
			MathExtra::transpose_matvec(rot,body[i].angmom,wbody);
			if (inertia[0] == 0.0) wbody[0] = 0.0;
			else wbody[0] /= inertia[0];
			if (inertia[1] == 0.0) wbody[1] = 0.0;
			else wbody[1] /= inertia[1];
			if (inertia[2] == 0.0) wbody[2] = 0.0;
			else wbody[2] /= inertia[2];

			t += inertia[0]wbody[0]wbody[0] + inertia[1]wbody[1]wbody[1] +
			inertia[2]wbody[2]wbody[2];
			}

			double tall;
			MPI_Allreduce(&t,&tall,1,MPI_DOUBLE,MPI_SUM,world);

			double tfactor = force->mvv2e / (6.0nbody force->boltz);
			tall *= tfactor;
			return tall;
			}

			/* ----------------------------------------------------------------------
			memory usage of local atom-based arrays
			------------------------------------------------------------------------- */

			double FixRigidSmall::memory_usage()
			{
			int nmax = atom->nmax;
			double bytes = 2 * nmax * sizeof(int);
			bytes += nmax3 sizeof(double);
			bytes += maxvatom6 sizeof(double); // vatom
			if (extended) {
			bytes += nmax * sizeof(int);
			if (orientflag) bytes = nmaxorientflag sizeof(double);
			if (dorientflag) bytes = nmax3 sizeof(double);
			}
			bytes += nmax_body * sizeof(Body);
			return bytes;
			}

			/* ----------------------------------------------------------------------
			debug method for sanity checking of atom/body data pointers
			------------------------------------------------------------------------- */

			/*
			void FixRigidSmall::check(int flag)
			{
			for (int i = 0; i < atom->nlocal; i++) {
			if (bodyown[i] >= 0) {
			if (bodytag[i] != atom->tag[i]) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD AAA");
			}
			if (bodyown[i] < 0 \|\| bodyown[i] >= nlocal_body) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD BBB");
			}
			if (atom2body[i] != bodyown[i]) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD CCC");
			}
			if (body[bodyown[i]].ilocal != i) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD DDD");
			}
			}
			}

			for (int i = 0; i < atom->nlocal; i++) {
			if (bodyown[i] < 0 && bodytag[i] > 0) {
			if (atom2body[i] < 0 \|\| atom2body[i] >= nlocal_body+nghost_body) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD EEE");
			}
			if (bodytag[i] != atom->tag[body[atom2body[i]].ilocal]) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD FFF");
			}
			}
			}

			for (int i = atom->nlocal; i < atom->nlocal + atom->nghost; i++) {
			if (bodyown[i] >= 0) {
			if (bodyown[i] < nlocal_body \|\|
			bodyown[i] >= nlocal_body+nghost_body) {
			printf("Values %d %d: %d %d %d\n",
			i,atom->tag[i],bodyown[i],nlocal_body,nghost_body);
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD GGG");
			}
			if (body[bodyown[i]].ilocal != i) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD HHH");
			}
			}
			}

			for (int i = 0; i < nlocal_body; i++) {
			if (body[i].ilocal < 0 \|\| body[i].ilocal >= atom->nlocal) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD III");
			}
			if (bodytag[body[i].ilocal] != atom->tag[body[i].ilocal] \|\|
			bodyown[body[i].ilocal] != i) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD JJJ");
			}
			}

			for (int i = nlocal_body; i < nlocal_body + nghost_body; i++) {
			if (body[i].ilocal < atom->nlocal \|\|
			body[i].ilocal >= atom->nlocal + atom->nghost) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD KKK");
			}
			if (bodyown[body[i].ilocal] != i) {
			printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
			errorx->one(FLERR,"BAD LLL");
			}
			}
			}
			*/

lammps/src/RIGID/fix_rigid_small.cpp8bee4ef5966alammm-devel

fix_rigid_small.cppView Options

lammps/src/RIGID/fix_rigid_small.cpp
8bee4ef5966alammm-devel

fix_rigid_small.cpp
View Options