diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index edcff5dc6..82ec19437 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -1,390 +1,391 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Mike Brown (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "string.h" #include "compute_temp_asphere.h" #include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{ROTATE,ALL}; #define INERTIA 0.2 // moment of inertia prefactor for ellipsoid /* ---------------------------------------------------------------------- */ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 3) error->all(FLERR,"Illegal compute temp/asphere command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 0; id_bias = NULL; mode = ALL; int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"bias") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command"); tempbias = 1; int n = strlen(arg[iarg+1]) + 1; id_bias = new char[n]; strcpy(id_bias,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"dof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command"); if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE; else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL; else error->all(FLERR,"Illegal compute temp/asphere command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/asphere command"); } vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempAsphere::~ComputeTempAsphere() { delete [] id_bias; delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::init() { // error check avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); if (!avec) error->all(FLERR,"Compute temp/asphere requires atom style ellipsoid"); // check that all particles are finite-size, no point particles allowed int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (ellipsoid[i] < 0) error->one(FLERR,"Compute temp/asphere requires extended particles"); if (tempbias) { int i = modify->find_compute(id_bias); if (i < 0) error->all(FLERR,"Could not find compute ID for temperature bias"); tbias = modify->compute[i]; if (tbias->tempflag == 0) error->all(FLERR,"Bias compute does not calculate temperature"); if (tbias->tempbias == 0) error->all(FLERR,"Bias compute does not calculate a velocity bias"); if (tbias->igroup != igroup) error->all(FLERR,"Bias compute group does not match compute group"); tbias->init(); + tbias->setup(); if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; else tempbias = 1; } } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::setup() { fix_dof = 0; for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::dof_compute() { // 6 dof for 3d, 3 dof for 2d // which dof are included also depends on mode // assume full rotation of extended particles // user should correct this via compute_modify if needed double natoms = group->count(igroup); int nper; if (domain->dimension == 3) { if (mode == ALL) nper = 6; else nper = 3; } else { if (mode == ALL) nper = 3; else nper = 1; } dof = nper*natoms; // additional adjustments to dof if (tempbias == 1) { if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms; } else if (tempbias == 2) { int *mask = atom->mask; int nlocal = atom->nlocal; int count = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (tbias->dof_remove(i)) count++; int count_all; MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof -= nper*count_all; } dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempAsphere::compute_scalar() { invoked_scalar = update->ntimestep; if (tempbias) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } AtomVecEllipsoid::Bonus *bonus = avec->bonus; double **v = atom->v; double **angmom = atom->angmom; double *rmass = atom->rmass; int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; double *shape,*quat; double wbody[3],inertia[3]; double rot[3][3]; // sum translational and rotational energy for each particle // no point particles since divide by inertia double t = 0.0; if (mode == ALL) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; t += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; t += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic || tempbias == 2) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::compute_vector() { int i; invoked_vector = update->ntimestep; if (tempbias) { if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } AtomVecEllipsoid::Bonus *bonus = avec->bonus; double **v = atom->v; double **angmom = atom->angmom; double *rmass = atom->rmass; int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; double *shape,*quat; double wbody[3],inertia[3],t[6]; double rot[3][3]; double massone; // sum translational and rotational energy for each particle // no point particles since divide by inertia for (i = 0; i < 6; i++) t[i] = 0.0; if (mode == ALL) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = rmass[i]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; // rotational kinetic energy t[0] += inertia[0]*wbody[0]*wbody[0]; t[1] += inertia[1]*wbody[1]*wbody[1]; t[2] += inertia[2]*wbody[2]*wbody[2]; t[3] += inertia[0]*wbody[0]*wbody[1]; t[4] += inertia[1]*wbody[0]*wbody[2]; t[5] += inertia[2]*wbody[1]*wbody[2]; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; massone = rmass[i]; inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; // rotational kinetic energy t[0] += inertia[0]*wbody[0]*wbody[0]; t[1] += inertia[1]*wbody[1]*wbody[1]; t[2] += inertia[2]*wbody[2]*wbody[2]; t[3] += inertia[0]*wbody[0]*wbody[1]; t[4] += inertia[1]*wbody[0]*wbody[2]; t[5] += inertia[2]*wbody[1]*wbody[2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempAsphere::remove_bias(int i, double *v) { if (tbias) tbias->remove_bias(i,v); } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempAsphere::restore_bias(int i, double *v) { if (tbias) tbias->restore_bias(i,v); } diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index daf58cd4f..261eeef81 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -1,2550 +1,2550 @@ /* ---------------------------------------------------------------------- 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 "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" 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 1000 +#define DELTA_BODY 10000 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,ibody; scalar_flag = 1; extscalar = 0; 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"); // parse optional args int seed; langflag = 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 = atof(arg[iarg+1]); t_stop = atof(arg[iarg+2]); t_period = atof(arg[iarg+3]); seed = atoi(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 error->all(FLERR,"Illegal fix rigid/small command"); } // create rigid bodies based on molecule ID // sets bodytag for owned atoms // body attributes are computed later by setup_bodies() if (atom->molecule_flag == 0) error->all(FLERR,"Fix rigid/small requires atom attribute molecule"); create_bodies(); // set nlocal_body and allocate bodies I own int *tag = atom->tag; int nlocal = atom->nlocal; 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_body*sizeof(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 (nbody == 0) error->all(FLERR,"No rigid bodies defined"); if (me == 0) { if (screen) fprintf(screen, "%d rigid bodies with " BIGINT_FORMAT " atoms\n", nbody,atomall); if (logfile) fprintf(logfile, "%d rigid bodies with " BIGINT_FORMAT " atoms\n", nbody,atomall); } // initialize Marsaglia RNG with processor-unique seed maxlang = 0; langextra = NULL; random = NULL; if (langflag) random = new RanMars(lmp,seed + comm->me); // 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; memory->destroy(langextra); } /* ---------------------------------------------------------------------- */ 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,ibody; 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; } /* ---------------------------------------------------------------------- */ void FixRigidSmall::setup(int vflag) { int i,n,ibody; double massone,radone; // one-time initialization of rigid body attributes via local comm // extended flags, mass, COM, inertia tensor, displacement of each atom if (firstflag) { firstflag = 0; setup_bodies(); } //check(1); // sum vcm, fcm, angmom, torque across all rigid bodies // vcm = velocity of COM // fcm = force on COM // angmom = angular momentum around COM // torque = torque around COM 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 *image = atom->image; int nlocal = atom->nlocal; double *xcm,*vcm,*fcm,*acm,*tcm; 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; fcm = body[ibody].fcm; fcm[0] = fcm[1] = fcm[2] = 0.0; acm = body[ibody].angmom; acm[0] = acm[1] = acm[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]]; 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; 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]; acm = b->angmom; acm[0] += dy * massone*v[i][2] - dz * massone*v[i][1]; acm[1] += dz * massone*v[i][0] - dx * massone*v[i][2]; acm[2] += dx * massone*v[i][1] - dy * massone*v[i][0]; 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) { AtomVecLine::Bonus *lbonus; if (avec_line) lbonus = avec_line->bonus; double **omega = atom->omega; double **angmom = atom->angmom; double **torque = atom->torque; 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] += SINERTIA*rmass[i] * radone*radone * omega[i][0]; acm[1] += SINERTIA*rmass[i] * radone*radone * omega[i][1]; acm[2] += SINERTIA*rmass[i] * radone*radone * omega[i][2]; } else if (eflags[i] & LINE) { radone = lbonus[line[i]].length; b->angmom[2] += LINERTIA*rmass[i] * radone*radone * 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]; } if (eflags[i] & TORQUE) { tcm = b->torque; tcm[0] += torque[i][0]; tcm[1] += torque[i][1]; tcm[2] += torque[i][2]; } } } // reverse communicate vcm, fcm, angmom, torque of all bodies commflag = VCM_ANGMOM; comm->reverse_comm_variable_fix(this); commflag = FORCE_TORQUE; 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; } // virial setup before call to set_v if (vflag) v_setup(vflag); else evflag = 0; // compute and forward communicate updated info 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] = gamma1*vcm[0] + gamma2*(random->uniform()-0.5); langextra[ibody][1] = gamma1*vcm[1] + gamma2*(random->uniform()-0.5); langextra[ibody][2] = gamma1*vcm[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]*gamma1*omega[0] + sqrt(inertia[0])*gamma2*(random->uniform()-0.5); langextra[ibody][4] = inertia[1]*gamma1*omega[1] + sqrt(inertia[1])*gamma2*(random->uniform()-0.5); langextra[ibody][5] = inertia[2]*gamma1*omega[2] + sqrt(inertia[2])*gamma2*(random->uniform()-0.5); } } /* ---------------------------------------------------------------------- */ void FixRigidSmall::final_integrate() { int i,ibody; double dtfm,xy,xz,yz; //check(3); // sum over atoms to get force and torque on rigid body tagint *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] += 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 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[atom2body[i]].fcm; fcm[0] += langextra[ibody][0]; fcm[1] += langextra[ibody][1]; fcm[2] += langextra[ibody][2]; tcm = body[atom2body[i]].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() { // do nothing if bodies are not yet initialized // i.e. pre_neighbor() is being called from Verlet::setup() if (firstflag) return; // remap xcm and image flags of each body as needed tagint 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 tagint *image = atom->image; int nlocal = atom->nlocal; int ibody; tagint 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; int tgroupbit = group->bitmask[tgroup]; // if firstflag, just return 0 // is being called by first run init() // local comm stencil is not setup and rigid body inertia is not calculated // will be called again by temperature compute setup() if (firstflag) return 0; // 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 += 3*counts[ibody][0] + 6*counts[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 += 2*counts[ibody][0] + 3*counts[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 ibody,itype; 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; tagint *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] + xbox*xprd + ybox*xy + zbox*xz; x1 = x[i][1] + ybox*yprd + zbox*yz; 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] - xbox*xprd - ybox*xy - zbox*xz; x[i][1] += b->xcm[1] - ybox*yprd - zbox*yz; 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.5*x0*fc0; vr[1] = 0.5*x1*fc1; vr[2] = 0.5*x2*fc2; vr[3] = 0.5*x0*fc1; vr[4] = 0.5*x0*fc2; vr[5] = 0.5*x1*fc2; 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] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]); ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]); ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]); MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(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 ibody,itype; int xbox,ybox,zbox; double dx,dy,dz; 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; tagint *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] + xbox*xprd + ybox*xy + zbox*xz; x1 = x[i][1] + ybox*yprd + zbox*yz; x2 = x[i][2] + zbox*zprd; } vr[0] = 0.5*x0*fc0; vr[1] = 0.5*x1*fc1; vr[2] = 0.5*x2*fc2; vr[3] = 0.5*x0*fc1; vr[4] = 0.5*x0*fc2; vr[5] = 0.5*x1*fc2; 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] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]); ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]); ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]); MathExtra::q_to_exyz(quatatom,exone,eyone,ezone); MathExtra::omega_to_angmom(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 int *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(); hash->clear(); // setup hash // key = body ID // value = index into N-length data structure // n = count of unique bodies my atoms are part of int *molecule = atom->molecule; n = 0; for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; if (hash->find(molecule[i]) == hash->end()) { hash->insert(std::pair (molecule[i],n)); 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,bbox_update,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) // rsq = 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 int *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++] = tag[i]; buf[m++] = unwrap[0]; buf[m++] = unwrap[1]; buf[m++] = unwrap[2]; } // pass buffer around ring of procs // func = 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,nearest_update,NULL); // set bodytag of all owned atoms, based on idclose for (i = 0; i < nlocal; i++) { bodytag[i] = 0; if (!(mask[i] & groupbit)) continue; m = hash->find(molecule[i])->second; bodytag[i] = idclose[m]; } // 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::bbox_update(int n, char *cbuf) { std::map *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 (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::nearest_update(int n, char *cbuf) { std::map *hash = frsptr->hash; double **ctr = frsptr->ctr; int *idclose = frsptr->idclose; double *rsqclose = frsptr->rsqclose; double *buf = (double *) cbuf; int ndatums = n/5; int j,imol,tag; double delx,dely,delz,rsq; double *x; int m = 0; for (int i = 0; i < ndatums; i++, m += 5) { imol = static_cast (buf[m]); if (hash->find(imol) != hash->end()) { j = hash->find(imol)->second; tag = static_cast (buf[m+1]); x = &buf[m+2]; delx = x[0] - ctr[j][0]; dely = x[1] - ctr[j][1]; delz = x[2] - ctr[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq <= rsqclose[j]) { if (rsq == rsqclose[j] && tag > idclose[j]) continue; idclose[j] = tag; rsqclose[j] = rsq; } } } } /* ---------------------------------------------------------------------- one-time initialization of rigid body attributes extended flags, mass, center-of-mass Cartesian and diagonalized inertia tensor ------------------------------------------------------------------------- */ void FixRigidSmall::setup_bodies() { int i,itype,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 (atom2body[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); } // 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 (atom2body[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 // error if image flag is not 0 in a non-periodic dim double **x = atom->x; tagint *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; } // 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 = ((tagint) IMGMAX << IMG2BITS) | ((tagint) 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 double dx,dy,dz,rad; double *inertia; for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) for (i = 0; i < 6; i++) body[ibody].itensor[i] = 0.0; 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 = b->itensor; inertia[0] += massone * (dy*dy + dz*dz); inertia[1] += massone * (dx*dx + dz*dz); inertia[2] += massone * (dx*dx + dy*dy); 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 = body[atom2body[i]].itensor; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if (eflags[i] & SPHERE) { inertia[0] += SINERTIA*massone * radius[i]*radius[i]; inertia[1] += SINERTIA*massone * radius[i]*radius[i]; inertia[2] += SINERTIA*massone * radius[i]*radius[i]; } else if (eflags[i] & ELLIPSOID) { shape = ebonus[ellipsoid[i]].shape; quatatom = ebonus[ellipsoid[i]].quat; MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec); 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]; } } } for (ibody = 0; ibody < nlocal_body; ibody++) { inertia = body[atom2body[i]].itensor; } // reverse communicate inertia tensor of all bodies commflag = ITENSOR; comm->reverse_comm_variable_fix(this); // 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 *itensor,*ex,*ey,*ez; for (ibody = 0; ibody < nlocal_body; ibody++) { itensor = body[ibody].itensor; tensor[0][0] = itensor[0]; tensor[1][1] = itensor[1]; tensor[2][2] = itensor[2]; tensor[1][2] = tensor[2][1] = itensor[3]; tensor[0][2] = tensor[2][0] = itensor[4]; tensor[0][1] = tensor[1][0] = itensor[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++) body[ibody].itensor[i] = 0.0; for (i = 0; i < nlocal; i++) { if (atom2body[i] < 0) continue; inertia = body[atom2body[i]].itensor; 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 = body[atom2body[i]].itensor; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; if (eflags[i] & SPHERE) { inertia[0] += SINERTIA*massone * radius[i]*radius[i]; inertia[1] += SINERTIA*massone * radius[i]*radius[i]; inertia[2] += SINERTIA*massone * 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 double *inew; double norm; for (ibody = 0; ibody < nlocal_body; ibody++) { inertia = body[ibody].inertia; itensor = body[ibody].itensor; if (inertia[0] == 0.0) { if (fabs(itensor[0]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((itensor[0]-inertia[0])/inertia[0]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } if (inertia[1] == 0.0) { if (fabs(itensor[1]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((itensor[1]-inertia[1])/inertia[1]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } if (inertia[2] == 0.0) { if (fabs(itensor[2]) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } else { if (fabs((itensor[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[3]/norm) > TOLERANCE || fabs(itensor[4]/norm) > TOLERANCE || fabs(itensor[5]/norm) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixRigidSmall::memory_usage() { int nmax = atom->nmax; double bytes = 2 * nmax * sizeof(int); bytes += nmax*3 * sizeof(double); bytes += maxvatom*6 * sizeof(double); // vatom if (extended) { bytes += nmax * sizeof(int); if (orientflag) bytes = nmax*orientflag * sizeof(double); if (dorientflag) bytes = nmax*3 * sizeof(double); } bytes += nmax_body * sizeof(Body); return bytes; } /* ---------------------------------------------------------------------- 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; } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ int FixRigidSmall::pack_exchange(int i, double *buf) { buf[0] = bodytag[i]; 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] = static_cast (buf[0]); 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 (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]) return m; // atom does not own its rigid body bodyown[nlocal] = static_cast (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; int nlocal = atom->nlocal; 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]; } } 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]; } } 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,flag; double *xcm,*vcm,*quat,*omega,*ex_space,*ey_space,*ez_space; 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++]; } } 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++]; } } else if (commflag == FULL_BODY) { for (i = first; i < last; i++) { bodyown[i] = static_cast (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,*itensor; 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; itensor = body[bodyown[i]].itensor; buf[m++] = itensor[0]; buf[m++] = itensor[1]; buf[m++] = itensor[2]; buf[m++] = itensor[3]; buf[m++] = itensor[4]; buf[m++] = itensor[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,*itensor; int *tag = atom->tag; int nlocal = atom->nlocal; 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; itensor = body[bodyown[j]].itensor; itensor[0] += buf[m++]; itensor[1] += buf[m++]; itensor[2] += buf[m++]; itensor[3] += buf[m++]; itensor[4] += buf[m++]; itensor[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] += buf[m++]; counts[k][1] += buf[m++]; counts[k][2] += buf[m++]; } } } /* ---------------------------------------------------------------------- grow body data structure ------------------------------------------------------------------------- */ void FixRigidSmall::grow_body() { nmax_body += DELTA_BODY; body = (Body *) memory->srealloc(body,nmax_body*sizeof(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 %d %d 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; } /* ---------------------------------------------------------------------- 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); error->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); error->one(FLERR,"BAD BBB"); } if (atom2body[i] != bodyown[i]) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); error->one(FLERR,"BAD CCC"); } if (body[bodyown[i]].ilocal != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); error->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); error->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); error->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); error->one(FLERR,"BAD GGG"); } if (body[bodyown[i]].ilocal != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); error->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); error->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); error->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); error->one(FLERR,"BAD KKK"); } if (bodyown[body[i].ilocal] != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); error->one(FLERR,"BAD LLL"); } } */ }