diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index 3dd2bbb28..48b443f04 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -1,401 +1,401 @@ /* ---------------------------------------------------------------------- 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 "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"); if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; else tempbias = 1; // init and setup bias compute because // this compute's setup()->dof_compute() may be called first tbias->init(); tbias->setup(); } } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::dof_compute() { adjust_dof_fix(); // 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; tbias->dof_remove_pre(); 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.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); 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(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); 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/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp index 73edbd131..807a132e9 100644 --- a/src/CORESHELL/compute_temp_cs.cpp +++ b/src/CORESHELL/compute_temp_cs.cpp @@ -1,499 +1,499 @@ /* ---------------------------------------------------------------------- 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: Hendrik Heenen (Technical University of Munich) (hendrik.heenen at mytum.com) ------------------------------------------------------------------------- */ #include "mpi.h" #include "stdlib.h" #include "string.h" #include "math.h" #include "compute_temp_cs.h" #include "atom.h" #include "atom_vec.h" #include "domain.h" #include "update.h" #include "force.h" #include "group.h" #include "modify.h" #include "fix.h" #include "fix_store.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 5) error->all(FLERR,"Illegal compute temp/cs command"); if (atom->avec->bonds_allow == 0) error->all(FLERR,"Compute temp/cs used when bonds are not allowed"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; extarray = 0; // find and define groupbits for core and shell groups cgroup = group->find(arg[3]); if (cgroup == -1) error->all(FLERR,"Cannot find specified group ID for core particles"); groupbit_c = group->bitmask[cgroup]; sgroup = group->find(arg[4]); if (sgroup == -1) error->all(FLERR,"Cannot find specified group ID for shell particles"); groupbit_s = group->bitmask[sgroup]; // create a new fix STORE style // id = compute-ID + COMPUTE_STORE, fix group = compute group int n = strlen(id) + strlen("_COMPUTE_STORE") + 1; id_fix = new char[n]; strcpy(id_fix,id); strcat(id_fix,"_COMPUTE_STORE"); char **newarg = new char*[5]; newarg[0] = id_fix; newarg[1] = group->names[igroup]; newarg[2] = (char *) "STORE"; newarg[3] = (char *) "0"; newarg[4] = (char *) "1"; modify->add_fix(5,newarg); fix = (FixStore *) modify->fix[modify->nfix-1]; delete [] newarg; // set fix store values = 0 for now // fill them in via setup() once Comm::borders() has been called // skip if resetting from restart file if (fix->restart_reset) { fix->restart_reset = 0; firstflag = 0; } else { double *partner = fix->vstore; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) partner[i] = ubuf(0).d; firstflag = 1; } // allocate memory vector = new double[6]; maxatom = 0; vint = NULL; // set comm size needed by this Compute comm_reverse = 1; } /* ---------------------------------------------------------------------- */ ComputeTempCS::~ComputeTempCS() { // check nfix in case all fixes have already been deleted if (modify->nfix) modify->delete_fix(id_fix); delete [] id_fix; delete [] vector; memory->destroy(vint); } /* ---------------------------------------------------------------------- */ void ComputeTempCS::init() { if (comm->ghost_velocity == 0) error->all(FLERR,"Compute temp/cs requires ghost atoms store velocity"); } /* ---------------------------------------------------------------------- */ void ComputeTempCS::setup() { if (firstflag) { firstflag = 0; // insure # of core atoms = # of shell atoms int ncores = group->count(cgroup); nshells = group->count(sgroup); if (ncores != nshells) error->all(FLERR,"Number of core atoms != number of shell atoms"); // for each C/S pair: // set partner IDs of both atoms if this atom stores bond between them // will set partner IDs for ghost atoms if needed by another proc // nall loop insures all ghost atom partner IDs are set before reverse comm int *num_bond = atom->num_bond; tagint **bond_atom = atom->bond_atom; tagint *tag = atom->tag; int *mask = atom->mask; int nlocal = atom->nlocal; double *partner = fix->vstore; tagint partnerID; int nall = nlocal + atom->nghost; for (int i = nlocal; i < nall; i++) partner[i] = ubuf(0).d; int i,j,m,match; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit_c || mask[i] & groupbit_s) { for (m = 0; m < num_bond[i]; m++) { partnerID = bond_atom[i][m]; j = atom->map(partnerID); if (j == -1) error->one(FLERR,"Core/shell partner atom not found"); match = 0; if (mask[i] & groupbit_c && mask[j] & groupbit_s) match = 1; if (mask[i] & groupbit_s && mask[j] & groupbit_c) match = 1; if (match) { partner[i] = ubuf(partnerID).d; partner[j] = ubuf(tag[i]).d; } } } } // reverse comm to acquire unknown partner IDs from ghost atoms // only needed if newton_bond = on if (force->newton_bond) comm->reverse_comm_compute(this); // check that all C/S partners were found int flag = 0; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit_c || mask[i] & groupbit_s) { partnerID = (tagint) ubuf(partner[i]).i; if (partnerID == 0) flag = 1; } } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) error->all(FLERR,"Core/shell partners were not all found"); } // calculate DOF for temperature dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempCS::dof_compute() { adjust_dof_fix(); int nper = domain->dimension; double natoms = group->count(igroup); dof = nper * natoms; dof -= nper * nshells; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempCS::compute_scalar() { double vthermal[3]; invoked_scalar = update->ntimestep; vcm_pairs(); // calculate thermal scalar in respect to atom velocities as center-of-mass // velocities of its according core/shell pairs double **v = atom->v; int *mask = atom->mask; int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++){ if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vint[i][0]; vthermal[1] = v[i][1] - vint[i][1]; vthermal[2] = v[i][2] - vint[i][2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempCS::compute_vector() { double vthermal[3]; invoked_vector = update->ntimestep; vcm_pairs(); // calculate thermal vector in respect to atom velocities as center-of-mass // velocities of its according C/S pairs double **v = atom->v; int *mask = atom->mask; int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; int nlocal = atom->nlocal; double massone; double t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++){ if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vint[i][0]; vthermal[1] = v[i][1] - vint[i][1]; vthermal[2] = v[i][2] - vint[i][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- */ void ComputeTempCS::vcm_pairs() { int i,j; double massone,masstwo; double vcm[3]; // reallocate vint if necessary int nlocal = atom->nlocal; if (nlocal > maxatom) { memory->destroy(vint); maxatom = atom->nmax; memory->create(vint,maxatom,3,"temp/cs:vint"); } // vcm = COM velocity of each CS pair // vint = internal velocity of each C/S atom, used as bias double **v = atom->v; int *mask = atom->mask; int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; double *partner = fix->vstore; tagint partnerID; for (i = 0; i < nlocal; i++) { if ((mask[i] & groupbit) && (mask[i] & groupbit_c || mask[i] & groupbit_s)) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; vcm[0] = v[i][0]*massone; vcm[1] = v[i][1]*massone; vcm[2] = v[i][2]*massone; partnerID = (tagint) ubuf(partner[i]).i; j = atom->map(partnerID); if (j == -1) error->one(FLERR,"Core/shell partner atom not found"); if (rmass) masstwo = rmass[j]; else masstwo = mass[type[j]]; vcm[0] += v[j][0]*masstwo; vcm[1] += v[j][1]*masstwo; vcm[2] += v[j][2]*masstwo; vcm[0] /= (massone + masstwo); vcm[1] /= (massone + masstwo); vcm[2] /= (massone + masstwo); vint[i][0] = v[i][0] - vcm[0]; vint[i][1] = v[i][1] - vcm[1]; vint[i][2] = v[i][2] - vcm[2]; } else vint[i][0] = vint[i][1] = vint[i][2] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity thermal velocity in this case is COM velocity of C/S pair ------------------------------------------------------------------------- */ void ComputeTempCS::remove_bias(int i, double *v) { v[0] -= vint[i][0]; v[1] -= vint[i][1]; v[2] -= vint[i][2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity thermal velocity in this case is COM velocity of C/S pair ------------------------------------------------------------------------- */ void ComputeTempCS::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] -= vint[i][0]; v[i][1] -= vint[i][1]; v[i][2] -= vint[i][2]; } } /* ---------------------------------------------------------------------- reset thermal velocity of all atoms to be consistent with bias called from velocity command after it creates thermal velocities this resets each atom's velocity to COM velocity of C/S pair ------------------------------------------------------------------------- */ void ComputeTempCS::reapply_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; // recalculate current COM velocities vcm_pairs(); // zero vint after using ti so that Velocity call to restore_bias_all() // will not further alter the velocities within a C/S pair for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] -= vint[i][0]; v[i][1] -= vint[i][1]; v[i][2] -= vint[i][2]; vint[i][0] = 0.0; vint[i][1] = 0.0; vint[i][2] = 0.0; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempCS::restore_bias(int i, double *v) { v[0] += vint[i][0]; v[1] += vint[i][1]; v[2] += vint[i][2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempCS::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vint[i][0]; v[i][1] += vint[i][1]; v[i][2] += vint[i][2]; } } /* ---------------------------------------------------------------------- */ int ComputeTempCS::pack_reverse_comm(int n, int first, double *buf) { int i,m,last; double *partner = fix->vstore; m = 0; last = first + n; for (i = first; i < last; i++) buf[m++] = partner[i]; return m; } /* ---------------------------------------------------------------------- */ void ComputeTempCS::unpack_reverse_comm(int n, int *list, double *buf) { int i,j,m; double *partner = fix->vstore; tagint partnerID; m = 0; for (i = 0; i < n; i++) { j = list[i]; partnerID = (tagint) ubuf(buf[m++]).i; if (partnerID) partner[j] = ubuf(partnerID).d; } } /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ double ComputeTempCS::memory_usage() { double bytes = (bigint) maxatom * 3 * sizeof(double); return bytes; } diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 8b5697a81..e4cdaaca1 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -1,3463 +1,3471 @@ /* ---------------------------------------------------------------------- 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 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; dof_flag = 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; xcmimage = 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; onemols = 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"); onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; 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 (onemols) { for (int i = 0; i < nmol; i++) { if (onemols[i]->xflag == 0) error->all(FLERR,"Fix rigid/small molecule must have coordinates"); if (onemols[i]->typeflag == 0) error->all(FLERR,"Fix rigid/small molecule must have atom types"); // fix rigid/small uses center, masstotal, COM, inertia of molecule onemols[i]->compute_center(); onemols[i]->compute_mass(); onemols[i]->compute_com(); onemols[i]->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_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 (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; staticflag = 0; } /* ---------------------------------------------------------------------- */ 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(xcmimage); 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; } /* ---------------------------------------------------------------------- setup static/dynamic properties of rigid bodies, using current atom info allows resetting of atom properties like mass between runs - only do static initialization once if using an infile + only do static initialization must be done once + do not do again if using an infile or mol keyword + this is b/c the infile and Molecule class define the COM, etc + so should not reset those values cannot do this until now, b/c requires comm->setup() to have setup stencil invoke pre_neighbor() to insure body xcmimage flags are reset needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run setup_bodies_static() invokes pre_neighbor itself ------------------------------------------------------------------------- */ void FixRigidSmall::setup_pre_neighbor() { - if (!staticflag || !infile) setup_bodies_static(); + if (!staticflag || (!infile && !onemols)) setup_bodies_static(); else pre_neighbor(); - setup_bodies_dynamic(); + + //if (!staticflag || (!infile)) setup_bodies_static(); + //else pre_neighbor(); + + //setup_bodies_dynamic(); + if (!onemols) setup_bodies_dynamic(); } /* ---------------------------------------------------------------------- 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; 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],xcmimage[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_fix(this,6); // 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_fix(this,10); // 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_fix(this,26); // 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; //check(3); // sum over atoms to get force and torque on rigid body 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],xcmimage[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_fix(this,6); // 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_fix(this,10); // 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(); } /* ---------------------------------------------------------------------- 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 first-time definition of rigid body in setup_bodies_static() or due to box flip also adjust imagebody = rigid body image flags, due to xcm remap - also adjust rgdimage flags of all atoms in bodies for two effects + then communicate bodies so other procs will know of changes to body xcm + then adjust xcmimage flags of all atoms in bodies via image_shift() + for two effects (1) change in true image flags due to pbc() call during exchange (2) change in imagebody due to xcm remap - rgdimage flags are always -1,0,-1 so that body can be unwrapped + xcmimage flags are always -1,0,-1 so that body can be unwrapped around in-box xcm and stay close to simulation box - if don't do this, then a body could end up very far away - when unwrapped by true image flags, and set_xv() will - compute huge displacements every step to reset coords of + if just inferred unwrapped from atom image flags, + then a body could end up very far away + when unwrapped by true image flags + then set_xv() will compute huge displacements every step to reset coords of all the body atoms to be back inside the box, ditto for triclinic box flip note: so just want to avoid that numeric probem? ------------------------------------------------------------------------- */ void FixRigidSmall::pre_neighbor() { - // 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 + for (int ibody = 0; ibody < nlocal_body; ibody++) { + Body *b = &body[ibody]; + domain->remap(b->xcm,b->image); + } nghost_body = 0; commflag = FULL_BODY; comm->forward_comm_fix(this); reset_atom2body(); //check(4); - for (int ibody = 0; ibody < nlocal_body; ibody++) { - Body *b = &body[ibody]; - domain->remap(b->xcm,b->image); - } image_shift(); } /* ---------------------------------------------------------------------- reset body xcmimage flags of atoms in bodies xcmimage flags are relative to xcm so that body can be unwrapped xcmimage = true image flag - imagebody flag ------------------------------------------------------------------------- */ void FixRigidSmall::image_shift() { imageint tdim,bdim,xdim[3]; imageint *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { if (atom2body[i] < 0) continue; Body *b = &body[atom2body[i]]; tdim = image[i] & IMGMASK; bdim = b->image & IMGMASK; xdim[0] = IMGMAX + tdim - bdim; tdim = (image[i] >> IMGBITS) & IMGMASK; bdim = (b->image >> IMGBITS) & IMGMASK; xdim[1] = IMGMAX + tdim - bdim; tdim = image[i] >> IMG2BITS; bdim = b->image >> IMG2BITS; xdim[2] = IMGMAX + tdim - bdim; xcmimage[i] = (xdim[2] << IMG2BITS) | (xdim[1] << IMGBITS) | xdim[0]; } } /* ---------------------------------------------------------------------- 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 (!staticflag) { 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_fix(this,3); // 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 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; 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 = (xcmimage[i] & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (xcmimage[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 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; 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 = (xcmimage[i] & IMGMASK) - IMGMAX; ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (xcmimage[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 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(); 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 (onemols) { for (int i = 0; i < nmol; i++) maxextent = MAX(maxextent,onemols[i]->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 *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::ring_nearest(int n, char *cbuf) { std::map *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 (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 = delx*delx + dely*dely + 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 = delx*delx + dely*dely + delz*delz; frsptr->rsqfar = MAX(frsptr->rsqfar,rsq); } } /* ---------------------------------------------------------------------- initialization of rigid body attributes called at setup, so body/atom properties can be changed between runs - unless reading from infile, in which case only called once + unless reading from infile or defined by Molecule class, + in which case only called once sets extended flags, masstotal, center-of-mass sets Cartesian and diagonalized inertia tensor sets body image flags, but only on first call ------------------------------------------------------------------------- */ 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 // require all molecules in template to have consistent radiusflag if (onemols) { int radiusflag = onemols[0]->radiusflag; for (i = 1; i < nmol; i++) { if (onemols[i]->radiusflag != radiusflag) error->all(FLERR,"Inconsistent use of finite-size particles " "by molecule template molecules"); } if (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; } } // first-time setting of body xcmimage flags = true image flags if (!staticflag) { imageint *image = atom->image; for (i = 0; i < nlocal; i++) if (bodytag[i] >= 0) xcmimage[i] = image[i]; else xcmimage[i] = 0; } // 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_fix(this); reset_atom2body(); // compute mass & center-of-mass of each rigid body double **x = atom->x; 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],xcmimage[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_fix(this,4); 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); } // one-time set of rigid body image flags to default values // staticflag insures this is only done once, not on successive runs // then remap the xcm of each body back into simulation box // and reset body xcmimage flags via pre_neighbor() if (!staticflag) { 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],xcmimage[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 * (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 = itensor[atom2body[i]]; 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]; } } } // reverse communicate inertia tensor of all bodies commflag = ITENSOR; comm->reverse_comm_fix(this,6); // 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_fix(this,26); // 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],xcmimage[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] += 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_fix(this,6); // 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); // static properties have now been initialized once // used to prevent re-initialization which would re-read infile staticflag = 1; } /* ---------------------------------------------------------------------- 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 + 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; 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],xcmimage[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]; } // 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] += 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]; } } } // reverse communicate vcm, angmom of all bodies commflag = VCM_ANGMOM; comm->reverse_comm_fix(this,6); // 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(); 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[CHUNK*MAXLINE]; 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 inertia 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 (!staticflag) 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; if (me == 0) { MPI_Status status; MPI_Request request; 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 (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,MPI_STATUS_IGNORE); 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(xcmimage,nmax,"rigid/small:xcmimage"); 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"); } // check for regrow of vatom // must be done whether per-atom virial is accumulated on this step or not // b/c this is only time grow_array() may be called // need to regrow b/c vatom is calculated before and after atom migration if (nmax > maxvatom) { maxvatom = atom->nmax; memory->grow(vatom,maxvatom,6,"fix:vatom"); } } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ void FixRigidSmall::copy_arrays(int i, int j, int delflag) { bodytag[j] = bodytag[i]; xcmimage[j] = xcmimage[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]; } } // must also copy vatom if per-atom virial calculated on this timestep // since vatom is calculated before and after atom migration if (vflag_atom) for (int k = 0; k < 6; k++) vatom[j][k] = vatom[i][k]; // 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; xcmimage[i] = 0; displace[i][0] = 0.0; displace[i][1] = 0.0; displace[i][2] = 0.0; // must also zero vatom if per-atom virial calculated on this timestep // since vatom is calculated before and after atom migration if (vflag_atom) for (int k = 0; k < 6; k++) vatom[i][k] = 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, int imol, 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 + onemols[imol]->comatom; if (tag[i]-tagprev == onemols[imol]->comatom) bodyown[i] = nlocal_body; m = tag[i] - tagprev-1; displace[i][0] = onemols[imol]->dxbody[m][0]; displace[i][1] = onemols[imol]->dxbody[m][1]; displace[i][2] = onemols[imol]->dxbody[m][2]; if (extended) { eflags[i] = 0; if (onemols[imol]->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 = onemols[imol]->masstotal; // new COM = Q (onemols[imol]->xcm - onemols[imol]->center) + xgeom // Q = rotation matrix associated with quat MathExtra::quat_to_mat(quat,rotmat); MathExtra::sub3(onemols[imol]->com,onemols[imol]->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] = onemols[imol]->inertia[0]; b->inertia[1] = onemols[imol]->inertia[1]; b->inertia[2] = onemols[imol]->inertia[2]; // final quat is product of insertion quat and original quat // true even if insertion rotation was not around COM MathExtra::quatquat(quat,onemols[imol]->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] = ubuf(xcmimage[i]).d; buf[2] = displace[i][0]; buf[3] = displace[i][1]; buf[4] = displace[i][2]; // extended attribute info int m = 5; 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; // must also pack vatom if per-atom virial calculated on this timestep // since vatom is calculated before and after atom migration if (vflag_atom) for (int k = 0; k < 6; k++) buf[m++] = vatom[i][k]; // 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; xcmimage[nlocal] = (imageint) ubuf(buf[1]).i; displace[nlocal][0] = buf[2]; displace[nlocal][1] = buf[3]; displace[nlocal][2] = buf[4]; // extended attribute info int m = 5; 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]) { bodyown[nlocal] = -1; return m; } // must also unpack vatom if per-atom virial calculated on this timestep // since vatom is calculated before and after atom migration if (vflag_atom) for (int k = 0; k < 6; k++) vatom[nlocal][k] = buf[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_forward_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_forward_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 (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 (buf[m++]); counts[k][1] += static_cast (buf[m++]); counts[k][2] += static_cast (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 " 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_fix(this,10); // 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_fix(this,10); // 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 onemols; } // 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.0*nbody * force->boltz); tall *= tfactor; return tall; } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixRigidSmall::memory_usage() { int nmax = atom->nmax; double bytes = nmax*2 * sizeof(int); bytes += nmax * sizeof(imageint); 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; } /* ---------------------------------------------------------------------- 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"); } } } */ diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 06449da27..5dacdbf9a 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -1,319 +1,319 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(rigid/small,FixRigidSmall) #else #ifndef LMP_FIX_RIGID_SMALL_H #define LMP_FIX_RIGID_SMALL_H #include "fix.h" // replace this later #include namespace LAMMPS_NS { class FixRigidSmall : public Fix { public: // static variable for ring communication callback to access class data static FixRigidSmall *frsptr; FixRigidSmall(class LAMMPS *, int, char **); virtual ~FixRigidSmall(); virtual int setmask(); virtual void init(); virtual void setup(int); virtual void initial_integrate(int); void post_force(int); virtual void final_integrate(); void initial_integrate_respa(int, int, int); void final_integrate_respa(int, int); void write_restart_file(char *); void grow_arrays(int); void copy_arrays(int, int, int); void set_arrays(int); void set_molecule(int, tagint, int, double *, double *, double *); int pack_exchange(int, double *); int unpack_exchange(int, double *); int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); int pack_reverse_comm(int, int, double *); void unpack_reverse_comm(int, int *, double *); void setup_pre_neighbor(); void pre_neighbor(); int dof(int); void deform(int); void reset_dt(); void zero_momentum(); void zero_rotation(); void *extract(const char*, int &); double extract_ke(); double extract_erotational(); double compute_scalar(); double memory_usage(); protected: int me,nprocs; double dtv,dtf,dtq; double *step_respa; int triclinic; double MINUSPI,TWOPI; char *infile; // file to read rigid body attributes from int staticflag; // 1 if static body properties are setup, else 0 int commflag; // various modes of forward/reverse comm int nbody; // total # of rigid bodies tagint maxmol; // max mol-ID double maxextent; // furthest distance from body owner to body atom struct Body { double mass; // total mass of body double xcm[3]; // COM position double vcm[3]; // COM velocity double fcm[3]; // force on COM double torque[3]; // torque around COM double quat[4]; // quaternion for orientation of body double inertia[3]; // 3 principal components of inertia double ex_space[3]; // principal axes in space coords double ey_space[3]; double ez_space[3]; double angmom[3]; // space-frame angular momentum of body double omega[3]; // space-frame omega of body double conjqm[4]; // conjugate quaternion momentum imageint image; // image flags of xcm int remapflag[4]; // PBC remap flags int ilocal; // index of owning atom }; Body *body; // list of rigid bodies, owned and ghost int nlocal_body; // # of owned rigid bodies int nghost_body; // # of ghost rigid bodies int nmax_body; // max # of bodies that body can hold int bodysize; // sizeof(Body) in doubles // per-atom quantities // only defined for owned atoms, except bodyown for own+ghost int *bodyown; // index of body if atom owns a body, -1 if not tagint *bodytag; // ID of body this atom is in, 0 if none // ID = tag of atom that owns body int *atom2body; // index of owned/ghost body this atom is in, -1 if not // can point to original or any image of the body imageint *xcmimage; // internal image flags for atoms in rigid bodies // set relative to in-box xcm of each body double **displace; // displacement of each atom in body coords int *eflags; // flags for extended particles double **orient; // orientation vector of particle wrt rigid body double **dorient; // orientation of dipole mu wrt rigid body int extended; // 1 if any particles have extended attributes int orientflag; // 1 if particles store spatial orientation int dorientflag; // 1 if particles store dipole orientation int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags int OMEGA,ANGMOM,TORQUE; class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; class AtomVecTri *avec_tri; // temporary per-body storage int **counts; // counts of atom types in bodies double **itensor; // 6 space-frame components of inertia tensor // mass per body, accessed by granular pair styles double *mass_body; int nmax_mass; // Langevin thermostatting int langflag; // 0/1 = no/yes Langevin thermostat double t_start,t_stop,t_period; // thermostat params double **langextra; // Langevin thermostat forces and torques int maxlang; // max size of langextra class RanMars *random; // RNG int tstat_flag,pstat_flag; // 0/1 = no/yes thermostat/barostat int t_chain,t_iter,t_order; double p_start[3],p_stop[3]; double p_period[3],p_freq[3]; int p_flag[3]; int pcouple,pstyle; int p_chain; int allremap; // remap all atoms int dilate_group_bit; // mask for dilation group char *id_dilate; // group name to dilate double p_current[3],p_target[3]; // molecules added on-the-fly as rigid bodies class Molecule **onemols; int nmol; // class data used by ring communication callbacks std::map *hash; double **bbox; double **ctr; tagint *idclose; double *rsqclose; double rsqfar; void image_shift(); void set_xv(); void set_v(); void create_bodies(); void setup_bodies_static(); void setup_bodies_dynamic(); void readfile(int, double **, int *); void grow_body(); void reset_atom2body(); // callback functions for ring communication static void ring_bbox(int, char *); static void ring_nearest(int, char *); static void ring_farthest(int, char *); // debug //void check(int); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Fix rigid/small requires atom attribute molecule Self-explanatory. E: Fix rigid/small requires an atom map, see atom_modify Self-explanatory. E: Fix rigid/small langevin period must be > 0.0 Self-explanatory. E: Molecule template ID for fix rigid/small does not exist Self-explanatory. E: Fix rigid/small nvt/npt/nph dilate group ID does not exist Self-explanatory. E: Fix rigid/small molecule must have coordinates The defined molecule does not specify coordinates. E: Fix rigid/small molecule must have atom types The defined molecule does not specify atom types. W: More than one fix rigid It is not efficient to use fix rigid more than once. E: Rigid fix must come before NPT/NPH fix NPT/NPH fix must be defined in input script after all rigid fixes, else the rigid fix contribution to the pressure virial is incorrect. -W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized +W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized h This means the temperature associated with the rigid bodies may be incorrect on this timestep. W: Computing temperature of portions of rigid bodies The group defined by the temperature compute does not encompass all the atoms in one or more rigid bodies, so the change in degrees-of-freedom for the atoms in those partial rigid bodies will not be accounted for. E: Fix rigid/small atom has non-zero image flag in a non-periodic dimension Image flags for non-periodic dimensions should not be set. E: Inconsistent use of finite-size particles by molecule template molecules Not all of the molecules define a radius for their constituent particles. E: Insufficient Jacobi rotations for rigid body Eigensolve for rigid body was not sufficiently accurate. E: Fix rigid: Bad principal moments The principal moments of inertia computed for a rigid body are not within the required tolerances. E: Cannot open fix rigid/small infile %s The specified file cannot be opened. Check that the path and name are correct. E: Unexpected end of fix rigid/small file A read operation from the file failed. E: Incorrect rigid body format in fix rigid/small file The number of fields per line is not what expected. E: Invalid rigid body ID in fix rigid/small file The ID does not match the number of an existing ID of rigid bodies that are defined by the fix rigid/small command. E: Cannot open fix rigid restart file %s The specified file cannot be opened. Check that the path and name are correct. E: Rigid body atoms %d %d missing on proc %d at step %ld This means that an atom cannot find the atom that owns the rigid body it is part of, or vice versa. The solution is to use the communicate cutoff command to insure ghost atoms are acquired from far enough away to encompass the max distance printed when the fix rigid/small command was invoked. */ diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index 7472b4d23..92ebdb298 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -1,327 +1,327 @@ /* ---------------------------------------------------------------------- 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: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "stdlib.h" #include "compute_temp_deform_eff.h" #include "domain.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "group.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp /* ---------------------------------------------------------------------- */ ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/deform/eff command"); if (!atom->electron_flag) error->all(FLERR,"Compute temp/deform/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempDeformEff::~ComputeTempDeformEff() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::init() { // check fix deform remap settings int i; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP && comm->me == 0) error->warning(FLERR,"Using compute temp/deform/eff with inconsistent " "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) error->warning(FLERR, "Using compute temp/deform/eff with no fix deform defined"); } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; // just include nuclear dof int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (abs(spin[i]) == 1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempDeformEff::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo // vthermal = thermal velocity = v - vstream double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (mass) { t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2])* mass[type[i]]; if (abs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::compute_vector() { double lamda[3],vstream[3],vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; if (abs(spin[i])==1) { t[0] += mefactor * massone * ervel[i]*ervel[i]; t[1] += mefactor * massone * ervel[i]*ervel[i]; t[2] += mefactor * massone * ervel[i]*ervel[i]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeformEff::remove_bias(int i, double *v) { double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; domain->x2lamda(atom->x[i],lamda); vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity NOTE: only removes translational velocity bias from electrons ------------------------------------------------------------------------- */ void ComputeTempDeformEff::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/deform/eff:vbiasall"); } double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(atom->x[i],lamda); vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[i][0] -= vbiasall[i][0]; v[i][1] -= vbiasall[i][1]; v[i][2] -= vbiasall[i][2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeformEff::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeformEff::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempDeformEff::memory_usage() { double bytes = maxbias * sizeof(double); return bytes; } diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 9f1bafdd4..9b24e0a70 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -1,167 +1,167 @@ /* ---------------------------------------------------------------------- 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: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "stdlib.h" #include "compute_temp_eff.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Compute temp/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[size_vector]; } /* ---------------------------------------------------------------------- */ ComputeTempEff::~ComputeTempEff() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempEff::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempEff::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (abs(spin[i])==1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempEff::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double t = 0.0; if (mass) { 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]) * mass[type[i]]; if (abs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempEff::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = mass[type[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]; if (abs(spin[i])==1) { t[0] += mefactor*massone*ervel[i]*ervel[i]; t[1] += mefactor*massone*ervel[i]*ervel[i]; t[2] += mefactor*massone*ervel[i]*ervel[i]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index 4ad653b2f..5a472549a 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -1,137 +1,137 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "string.h" #include "compute_temp.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "comm.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTemp::~ComputeTemp() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTemp::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTemp::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTemp::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { 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]; } else { 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]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTemp::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[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]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp index 2d36bdf0b..afbe02129 100644 --- a/src/compute_temp_chunk.cpp +++ b/src/compute_temp_chunk.cpp @@ -1,854 +1,856 @@ /* ---------------------------------------------------------------------- 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 "string.h" #include "compute_temp_chunk.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{TEMP,KECOM,INTERNAL}; /* ---------------------------------------------------------------------- */ ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 4) error->all(FLERR,"Illegal compute temp/chunk command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; // ID of compute chunk/atom int n = strlen(arg[3]) + 1; idchunk = new char[n]; strcpy(idchunk,arg[3]); biasflag = 0; init(); // optional per-chunk values nvalues = narg-4; which = new int[nvalues]; nvalues = 0; int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"temp") == 0) which[nvalues] = TEMP; else if (strcmp(arg[iarg],"kecom") == 0) which[nvalues] = KECOM; else if (strcmp(arg[iarg],"internal") == 0) which[nvalues] = INTERNAL; else break; iarg++; nvalues++; } // optional args comflag = 0; biasflag = 0; id_bias = NULL; adof = domain->dimension; cdof = 0.0; while (iarg < narg) { if (strcmp(arg[iarg],"com") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); if (strcmp(arg[iarg+1],"yes") == 0) comflag = 1; else if (strcmp(arg[iarg+1],"no") == 0) comflag = 0; else error->all(FLERR,"Illegal compute temp/chunk command"); iarg += 2; } else if (strcmp(arg[iarg],"bias") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); biasflag = 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],"adof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); adof = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"cdof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); cdof = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else error->all(FLERR,"Illegal compute temp/chunk command"); } // error check on bias compute if (biasflag) { 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"); } // this compute only calculates a bias, if comflag is set // won't be two biases since comflag and biasflag cannot both be set if (comflag && biasflag) error->all(FLERR,"Cannot use both com and bias with compute temp/chunk"); if (comflag) tempbias = 1; // vector data vector = new double[6]; // chunk-based data nchunk = 1; maxchunk = 0; sum = sumall = NULL; count = countall = NULL; massproc = masstotal = NULL; vcm = vcmall = NULL; array = NULL; if (nvalues) { array_flag = 1; size_array_cols = nvalues; size_array_rows = 0; size_array_rows_variable = 1; extarray = 0; } allocate(); comstep = -1; } /* ---------------------------------------------------------------------- */ ComputeTempChunk::~ComputeTempChunk() { delete [] idchunk; delete [] which; delete [] id_bias; delete [] vector; memory->destroy(sum); memory->destroy(sumall); memory->destroy(count); memory->destroy(countall); memory->destroy(array); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(vcm); memory->destroy(vcmall); } /* ---------------------------------------------------------------------- */ void ComputeTempChunk::init() { int icompute = modify->find_compute(idchunk); if (icompute < 0) error->all(FLERR,"Chunk/atom compute does not exist for " "compute temp/chunk"); cchunk = (ComputeChunkAtom *) modify->compute[icompute]; if (strcmp(cchunk->style,"chunk/atom") != 0) error->all(FLERR,"Compute temp/chunk does not use chunk/atom compute"); if (biasflag) { 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]; } } /* ---------------------------------------------------------------------- */ double ComputeTempChunk::compute_scalar() { int i,index; invoked_scalar = update->ntimestep; // calculate chunk assignments, // since only atoms in chunks contribute to global temperature // compute chunk/atom assigns atoms to chunk IDs // extract ichunk index vector from compute // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms nchunk = cchunk->setup_chunks(); cchunk->compute_ichunk(); int *ichunk = cchunk->ichunk; // remove velocity bias if (biasflag) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } // calculate COM velocity for each chunk // won't be invoked with bias also removed = 2 biases if (comflag && comstep != update->ntimestep) vcm_compute(); // calculate global temperature, optionally removing COM velocity double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; int mycount = 0; if (!comflag) { if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; mycount++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; mycount++; } } } else { double vx,vy,vz; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; t += (vx*vx + vy*vy + vz*vz) * rmass[i]; mycount++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; t += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; mycount++; } } } // restore velocity bias if (biasflag) tbias->restore_bias_all(); // final temperature MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); double rcount = mycount; double allcount; MPI_Allreduce(&rcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world); double dof = nchunk*cdof + adof*allcount; + if (dof < 0.0 && allcount > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); double tfactor = 0.0; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempChunk::compute_vector() { int i,index; invoked_vector = update->ntimestep; // calculate chunk assignments, // since only atoms in chunks contribute to global temperature // compute chunk/atom assigns atoms to chunk IDs // extract ichunk index vector from compute // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms nchunk = cchunk->setup_chunks(); cchunk->compute_ichunk(); int *ichunk = cchunk->ichunk; // remove velocity bias if (biasflag) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } // calculate COM velocity for each chunk // won't be invoked with bias also removed = 2 biases if (comflag && comstep != update->ntimestep) vcm_compute(); // calculate KE tensor, optionally removing COM velocity double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; if (!comflag) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; if (rmass) massone = rmass[i]; else massone = mass[type[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]; } } else { double vx,vy,vz; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; t[0] += massone * vx*vx; t[1] += massone * vy*vy; t[2] += massone * vz*vz; t[3] += massone * vx*vy; t[4] += massone * vx*vz; t[5] += massone * vy*vz; } } // restore velocity bias if (biasflag) tbias->restore_bias_all(); // final KE MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- */ void ComputeTempChunk::compute_array() { invoked_array = update->ntimestep; // compute chunk/atom assigns atoms to chunk IDs // extract ichunk index vector from compute // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms nchunk = cchunk->setup_chunks(); cchunk->compute_ichunk(); if (nchunk > maxchunk) allocate(); size_array_rows = nchunk; // remove velocity bias if (biasflag) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } // calculate COM velocity for each chunk whether comflag set or not // needed by some values even if comflag not set // important to do this after velocity bias is removed // otherwise per-chunk values that use both v and vcm will be inconsistent if (comstep != update->ntimestep) vcm_compute(); // compute each value for (int i = 0; i < nvalues; i++) { if (which[i] == TEMP) temperature(i); else if (which[i] == KECOM) kecom(i); else if (which[i] == INTERNAL) internal(i); } // restore velocity bias if (biasflag) tbias->restore_bias_all(); } /* ---------------------------------------------------------------------- calculate velocity of COM for each chunk ------------------------------------------------------------------------- */ void ComputeTempChunk::vcm_compute() { int i,index; double massone; // avoid re-computing VCM more than once per step comstep = update->ntimestep; int *ichunk = cchunk->ichunk; for (int i = 0; i < nchunk; i++) { vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; massproc[i] = 0.0; } double **v = atom->v; int *mask = atom->mask; int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; vcm[index][0] += v[i][0] * massone; vcm[index][1] += v[i][1] * massone; vcm[index][2] += v[i][2] * massone; massproc[index] += massone; } MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < nchunk; i++) { vcmall[i][0] /= masstotal[i]; vcmall[i][1] /= masstotal[i]; vcmall[i][2] /= masstotal[i]; } } /* ---------------------------------------------------------------------- temperature of each chunk ------------------------------------------------------------------------- */ void ComputeTempChunk::temperature(int icol) { int i,index; int *ichunk = cchunk->ichunk; // zero local per-chunk values for (int i = 0; i < nchunk; i++) { count[i] = 0; sum[i] = 0.0; } // per-chunk temperature, option for removing COM velocity double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; if (!comflag) { if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; sum[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; count[index]++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; sum[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; count[index]++; } } } else { double vx,vy,vz; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; count[index]++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; count[index]++; } } } // sum across procs MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(count,countall,nchunk,MPI_INT,MPI_SUM,world); // normalize temperatures by per-chunk DOF double dof,tfactor; double mvv2e = force->mvv2e; double boltz = force->boltz; for (int i = 0; i < nchunk; i++) { dof = cdof + adof*countall[i]; if (dof > 0.0) tfactor = mvv2e / (dof * boltz); else tfactor = 0.0; array[i][icol] = tfactor * sumall[i]; } } /* ---------------------------------------------------------------------- KE of entire chunk moving at VCM ------------------------------------------------------------------------- */ void ComputeTempChunk::kecom(int icol) { int index; int *ichunk = cchunk->ichunk; // zero local per-chunk values for (int i = 0; i < nchunk; i++) sum[i] = 0.0; // per-chunk COM KE double *mass = atom->mass; double *rmass = atom->rmass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; double vx,vy,vz; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = vcmall[index][0]; vy = vcmall[index][1]; vz = vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = vcmall[index][0]; vy = vcmall[index][1]; vz = vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; } } // sum across procs MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); double mvv2e = force->mvv2e; for (int i = 0; i < nchunk; i++) array[i][icol] = 0.5 * mvv2e * sumall[i]; } /* ---------------------------------------------------------------------- internal KE of each chunk around its VCM computed using per-atom velocities with chunk VCM subtracted off ------------------------------------------------------------------------- */ void ComputeTempChunk::internal(int icol) { int index; int *ichunk = cchunk->ichunk; // zero local per-chunk values for (int i = 0; i < nchunk; i++) sum[i] = 0.0; // per-chunk internal KE double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; double vx,vy,vz; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; } } // sum across procs MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); double mvv2e = force->mvv2e; for (int i = 0; i < nchunk; i++) array[i][icol] = 0.5 * mvv2e * sumall[i]; } /* ---------------------------------------------------------------------- bias methods: called by thermostats ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempChunk::remove_bias(int i, double *v) { int index = cchunk->ichunk[i]; if (index < 0) return; v[0] -= vcmall[index][0]; v[1] -= vcmall[index][1]; v[2] -= vcmall[index][2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempChunk::remove_bias_all() { int index; int *ichunk = cchunk->ichunk; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]; if (index < 0) continue; v[i][0] -= vbias[0]; v[i][1] -= vbias[1]; v[i][2] -= vbias[2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempChunk::restore_bias(int i, double *v) { int index = cchunk->ichunk[i]; if (index < 0) return; v[0] += vcmall[index][0]; v[1] += vcmall[index][1]; v[2] += vcmall[index][2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempChunk::restore_bias_all() { int index; int *ichunk = cchunk->ichunk; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]; if (index < 0) continue; v[i][0] += vbias[0]; v[i][1] += vbias[1]; v[i][2] += vbias[2]; } } /* ---------------------------------------------------------------------- lock methods: called by fix ave/time these methods insure vector/array size is locked for Nfreq epoch by passing lock info along to compute chunk/atom ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- increment lock counter ------------------------------------------------------------------------- */ void ComputeTempChunk::lock_enable() { cchunk->lockcount++; } /* ---------------------------------------------------------------------- decrement lock counter in compute chunk/atom, it if still exists ------------------------------------------------------------------------- */ void ComputeTempChunk::lock_disable() { int icompute = modify->find_compute(idchunk); if (icompute >= 0) { cchunk = (ComputeChunkAtom *) modify->compute[icompute]; cchunk->lockcount--; } } /* ---------------------------------------------------------------------- calculate and return # of chunks = length of vector/array ------------------------------------------------------------------------- */ int ComputeTempChunk::lock_length() { nchunk = cchunk->setup_chunks(); return nchunk; } /* ---------------------------------------------------------------------- set the lock from startstep to stopstep ------------------------------------------------------------------------- */ void ComputeTempChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) { cchunk->lock(fixptr,startstep,stopstep); } /* ---------------------------------------------------------------------- unset the lock ------------------------------------------------------------------------- */ void ComputeTempChunk::unlock(Fix *fixptr) { cchunk->unlock(fixptr); } /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeTempChunk::allocate() { memory->destroy(sum); memory->destroy(sumall); memory->destroy(count); memory->destroy(countall); memory->destroy(array); maxchunk = nchunk; memory->create(sum,maxchunk,"temp/chunk:sum"); memory->create(sumall,maxchunk,"temp/chunk:sumall"); memory->create(count,maxchunk,"temp/chunk:count"); memory->create(countall,maxchunk,"temp/chunk:countall"); memory->create(array,maxchunk,nvalues,"temp/chunk:array"); if (comflag || nvalues) { memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(vcm); memory->destroy(vcmall); memory->create(massproc,maxchunk,"vcm/chunk:massproc"); memory->create(masstotal,maxchunk,"vcm/chunk:masstotal"); memory->create(vcm,maxchunk,3,"vcm/chunk:vcm"); memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall"); } } /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ double ComputeTempChunk::memory_usage() { double bytes = (bigint) maxchunk * 2 * sizeof(double); bytes = (bigint) maxchunk * 2 * sizeof(int); bytes = (bigint) maxchunk * nvalues * sizeof(double); if (comflag || nvalues) { bytes += (bigint) maxchunk * 2 * sizeof(double); bytes += (bigint) maxchunk * 2*3 * sizeof(double); } return bytes; } diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index 8e0db3595..e9f31dbfe 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -1,221 +1,220 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_com.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "domain.h" #include "lattice.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempCOM::~ComputeTempCOM() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::init() { masstotal = group->mass(igroup); } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); - int nper = domain->dimension; - dof = nper * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempCOM::compute_scalar() { double vthermal[3]; invoked_scalar = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vbias[0]; vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::compute_vector() { int i; double vthermal[3]; invoked_vector = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vbias[0]; vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } 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 ComputeTempCOM::remove_bias(int i, double *v) { v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempCOM::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] -= vbias[0]; v[i][1] -= vbias[1]; v[i][2] -= vbias[2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempCOM::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempCOM::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbias[0]; v[i][1] += vbias[1]; v[i][2] += vbias[2]; } } diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index e608ba3ce..7513f6fb1 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -1,294 +1,294 @@ /* ---------------------------------------------------------------------- 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: Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "string.h" #include "compute_temp_deform.h" #include "domain.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "group.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp /* ---------------------------------------------------------------------- */ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/deform command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempDeform::~ComputeTempDeform() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::init() { int i; // check fix deform remap settings for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP && comm->me == 0) error->warning(FLERR,"Using compute temp/deform with inconsistent " "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) error->warning(FLERR, "Using compute temp/deform with no fix deform defined"); } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempDeform::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo // vthermal = thermal velocity = v - vstream double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::compute_vector() { double lamda[3],vstream[3],vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeform::remove_bias(int i, double *v) { double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; domain->x2lamda(atom->x[i],lamda); vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeform::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/deform:vbiasall"); } double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(atom->x[i],lamda); vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[i][0] -= vbiasall[i][0]; v[i][1] -= vbiasall[i][1]; v[i][2] -= vbiasall[i][2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeform::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeform::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempDeform::memory_usage() { double bytes = 3*maxbias * sizeof(double); return bytes; } diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 2b13dd092..9f0ff3f31 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -1,564 +1,563 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_profile.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "fix.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{TENSOR,BIN}; /* ---------------------------------------------------------------------- */ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 7) error->all(FLERR,"Illegal compute temp/profile command"); scalar_flag = 1; extscalar = 0; tempflag = 1; tempbias = 1; xflag = force->inumeric(FLERR,arg[3]); yflag = force->inumeric(FLERR,arg[4]); zflag = force->inumeric(FLERR,arg[5]); if (zflag && domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot use vz for 2d systemx"); ncount = 0; ivx = ivy = ivz = 0; if (xflag) ivx = ncount++; if (yflag) ivy = ncount++; if (zflag) ivz = ncount++; ncount += 2; nbinx = nbiny = nbinz = 1; int iarg = 6; if (strcmp(arg[iarg],"x") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbinx = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"y") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbiny = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"z") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinz = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"xy") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbiny = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"yz") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbiny = force->inumeric(FLERR,arg[iarg+1]); nbinz = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"xz") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbinz = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"xyz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbiny = force->inumeric(FLERR,arg[iarg+2]); nbinz = force->inumeric(FLERR,arg[iarg+3]); iarg += 4; } else error->all(FLERR,"Illegal compute temp/profile command"); // optional keywords outflag = TENSOR; while (iarg < narg) { if (strcmp(arg[iarg],"out") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (strcmp(arg[iarg+1],"tensor") == 0) outflag = TENSOR; else if (strcmp(arg[iarg+1],"bin") == 0) outflag = BIN; else error->all(FLERR,"Illegal compute temp/profile command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/profile command"); } // setup nbins = nbinx*nbiny*nbinz; if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command"); memory->create(vbin,nbins,ncount,"temp/profile:vbin"); memory->create(binave,nbins,ncount,"temp/profile:binave"); if (outflag == TENSOR) { vector_flag = 1; size_vector = 6; extvector = 1; vector = new double[size_vector]; } else { array_flag = 1; size_array_rows = nbins; size_array_cols = 2; extarray = 0; memory->create(tbin,nbins,"temp/profile:tbin"); memory->create(tbinall,nbins,"temp/profile:tbinall"); memory->create(array,nbins,2,"temp/profile:array"); } maxatom = 0; bin = NULL; } /* ---------------------------------------------------------------------- */ ComputeTempProfile::~ComputeTempProfile() { memory->destroy(vbin); memory->destroy(binave); memory->destroy(bin); if (outflag == TENSOR) delete [] vector; else { memory->destroy(tbin); memory->destroy(tbinall); memory->destroy(array); } } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::init() { dof_compute(); // ptrs to domain data box_change = domain->box_change; triclinic = domain->triclinic; periodicity = domain->periodicity; if (triclinic) { boxlo = domain->boxlo_lamda; boxhi = domain->boxhi_lamda; prd = domain->prd_lamda; } else { boxlo = domain->boxlo; boxhi = domain->boxhi; prd = domain->prd; } if (!box_change) bin_setup(); } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); - int nper = domain->dimension; - dof = nper * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempProfile::compute_scalar() { int ibin; double vthermal[3]; invoked_scalar = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::compute_vector() { int i,ibin; double vthermal[3]; invoked_vector = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::compute_array() { int i,ibin; double vthermal[3]; invoked_array = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nbins; i++) tbin[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) tbin[ibin] += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else tbin[ibin] += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(tbin,tbinall,nbins,MPI_DOUBLE,MPI_SUM,world); int nper = domain->dimension; for (i = 0; i < nbins; i++) { array[i][0] = binave[i][ncount-1]; if (array[i][0] > 0.0) { dof = nper*array[i][0] - extra_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; array[i][1] = tfactor*tbinall[i]; } else array[i][1] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempProfile::remove_bias(int i, double *v) { int ibin = bin[i]; if (xflag) v[0] -= binave[ibin][ivx]; if (yflag) v[1] -= binave[ibin][ivy]; if (zflag) v[2] -= binave[ibin][ivz]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempProfile::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; int ibin; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) v[i][0] -= binave[ibin][ivx]; if (yflag) v[i][1] -= binave[ibin][ivy]; if (zflag) v[i][2] -= binave[ibin][ivz]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempProfile::restore_bias(int i, double *v) { int ibin = bin[i]; if (xflag) v[0] += binave[ibin][ivx]; if (yflag) v[1] += binave[ibin][ivy]; if (zflag) v[2] += binave[ibin][ivz]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempProfile::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; int ibin; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) v[i][0] += binave[ibin][ivx]; if (yflag) v[i][1] += binave[ibin][ivy]; if (zflag) v[i][2] += binave[ibin][ivz]; } } /* ---------------------------------------------------------------------- compute average COM velocity in each bin ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_average() { int i,j,ibin; if (box_change) bin_setup(); bin_assign(); // clear bins, including particle mass and count for (i = 0; i < nbins; i++) for (j = 0; j < ncount; j++) vbin[i][j] = 0.0; // sum each particle's mass-weighted velocity, mass, count to appropriate bin double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; int nc2 = ncount-2; int nc1 = ncount-1; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vbin[ibin][ivx] += rmass[i]*v[i][0]; if (yflag) vbin[ibin][ivy] += rmass[i]*v[i][1]; if (zflag) vbin[ibin][ivz] += rmass[i]*v[i][2]; vbin[ibin][nc2] += rmass[i]; vbin[ibin][nc1] += 1.0; } } else { double onemass; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; onemass = mass[type[i]]; if (xflag) vbin[ibin][ivx] += onemass*v[i][0]; if (yflag) vbin[ibin][ivy] += onemass*v[i][1]; if (zflag) vbin[ibin][ivz] += onemass*v[i][2]; vbin[ibin][nc2] += onemass; vbin[ibin][nc1] += 1.0; } } // sum bins across processors MPI_Allreduce(vbin[0],binave[0],nbins*ncount,MPI_DOUBLE,MPI_SUM,world); // compute ave COM velocity in each bin, checking for no particles for (i = 0; i < nbins; i++) if (binave[i][nc1] > 0.0) for (j = 0; j < nc2; j++) binave[i][j] /= binave[i][nc2]; } /* ---------------------------------------------------------------------- set bin sizes, redo if box size changes ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_setup() { invdelta[0] = nbinx / prd[0]; invdelta[1] = nbiny / prd[1]; invdelta[2] = nbinz / prd[2]; } /* ---------------------------------------------------------------------- assign all atoms to bins ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_assign() { // reallocate bin array if necessary if (atom->nlocal > maxatom) { maxatom = atom->nmax; memory->destroy(bin); memory->create(bin,maxatom,"temp/profile:bin"); } // assign each atom to a bin, accounting for PBC // if triclinic, do this in lamda space double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; int ibinx,ibiny,ibinz; double coord; if (triclinic) domain->x2lamda(nlocal); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (nbinx > 1) { coord = x[i][0]; if (periodicity[0]) { if (coord < boxlo[0]) coord += prd[0]; if (coord >= boxhi[0]) coord -= prd[0]; } ibinx = static_cast ((coord - boxlo[0]) * invdelta[0]); ibinx = MAX(ibinx,0); ibinx = MIN(ibinx,nbinx-1); } else ibinx = 0; if (nbiny > 1) { coord = x[i][1]; if (periodicity[1]) { if (coord < boxlo[1]) coord += prd[1]; if (coord >= boxhi[1]) coord -= prd[1]; } ibiny = static_cast ((coord - boxlo[1]) * invdelta[1]); ibiny = MAX(ibiny,0); ibiny = MIN(ibiny,nbiny-1); } else ibiny = 0; if (nbinz > 1) { coord = x[i][2]; if (periodicity[2]) { if (coord < boxlo[2]) coord += prd[2]; if (coord >= boxhi[2]) coord -= prd[2]; } ibinz = static_cast ((coord - boxlo[2]) * invdelta[2]); ibinz = MAX(ibinz,0); ibinz = MIN(ibinz,nbinz-1); } else ibinz = 0; bin[i] = nbinx*nbiny*ibinz + nbinx*ibiny + ibinx; } if (triclinic) domain->lamda2x(nlocal); } /* ---------------------------------------------------------------------- */ double ComputeTempProfile::memory_usage() { double bytes = maxatom * sizeof(int); bytes += nbins*ncount * sizeof(double); return bytes; } diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 76b47a39a..0c38ed811 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -1,297 +1,296 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_ramp.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "fix.h" #include "domain.h" #include "lattice.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 9) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; // parse optional args scaleflag = 1; int iarg = 9; while (iarg < narg) { if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/ramp command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all(FLERR,"Illegal compute temp/ramp command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/ramp command"); } // setup scaling if (scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; // read standard args and apply scaling if (strcmp(arg[3],"vx") == 0) v_dim = 0; else if (strcmp(arg[3],"vy") == 0) v_dim = 1; else if (strcmp(arg[3],"vz") == 0) v_dim = 2; else error->all(FLERR,"Illegal compute temp/ramp command"); if (v_dim == 0) { v_lo = xscale*force->numeric(FLERR,arg[4]); v_hi = xscale*force->numeric(FLERR,arg[5]); } else if (v_dim == 1) { v_lo = yscale*force->numeric(FLERR,arg[4]); v_hi = yscale*force->numeric(FLERR,arg[5]); } else if (v_dim == 2) { v_lo = zscale*force->numeric(FLERR,arg[4]); v_hi = zscale*force->numeric(FLERR,arg[5]); } if (strcmp(arg[6],"x") == 0) coord_dim = 0; else if (strcmp(arg[6],"y") == 0) coord_dim = 1; else if (strcmp(arg[6],"z") == 0) coord_dim = 2; else error->all(FLERR,"Illegal compute temp/ramp command"); if (coord_dim == 0) { coord_lo = xscale*force->numeric(FLERR,arg[7]); coord_hi = xscale*force->numeric(FLERR,arg[8]); } else if (coord_dim == 1) { coord_lo = yscale*force->numeric(FLERR,arg[7]); coord_hi = yscale*force->numeric(FLERR,arg[8]); } else if (coord_dim == 2) { coord_lo = zscale*force->numeric(FLERR,arg[7]); coord_hi = zscale*force->numeric(FLERR,arg[8]); } maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempRamp::~ComputeTempRamp() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::dof_compute() { adjust_dof_fix(); double natoms = group->count(igroup); - int nper = domain->dimension; - dof = nper * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempRamp::compute_scalar() { double fraction,vramp,vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vramp = v_lo + fraction*(v_hi - v_lo); vthermal[0] = v[i][0]; vthermal[1] = v[i][1]; vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::compute_vector() { int i; double fraction,vramp,vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vramp = v_lo + fraction*(v_hi - v_lo); vthermal[0] = v[i][0]; vthermal[1] = v[i][1]; vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } 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 ComputeTempRamp::remove_bias(int i, double *v) { double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vbias[v_dim] = v_lo + fraction*(v_hi - v_lo); v[v_dim] -= vbias[v_dim]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRamp::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/ramp:vbiasall"); } double fraction; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vbiasall[i][v_dim] = v_lo + fraction*(v_hi - v_lo); v[i][v_dim] -= vbiasall[i][v_dim]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempRamp::restore_bias(int i, double *v) { v[v_dim] += vbias[v_dim]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempRamp::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][v_dim] += vbiasall[i][v_dim]; } /* ---------------------------------------------------------------------- */ double ComputeTempRamp::memory_usage() { double bytes = 3*maxbias * sizeof(double); return bytes; } diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 70f13d40b..050a98793 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -1,339 +1,337 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "string.h" #include "compute_temp_sphere.h" #include "atom.h" #include "atom_vec.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "comm.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; enum{ROTATE,ALL}; #define INERTIA 0.4 // moment of inertia prefactor for sphere /* ---------------------------------------------------------------------- */ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 3) error->all(FLERR,"Illegal compute temp/sphere 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/sphere 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/sphere 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/sphere command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/sphere command"); } vector = new double[6]; // error checks if (!atom->sphere_flag) error->all(FLERR,"Compute temp/sphere requires atom style sphere"); } /* ---------------------------------------------------------------------- */ ComputeTempSphere::~ComputeTempSphere() { delete [] id_bias; delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::init() { 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"); if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; else tempbias = 1; // init and setup bias compute because // this compute's setup()->dof_compute() may be called first tbias->init(); tbias->setup(); } } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::setup() { dynamic = 0; if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::dof_compute() { int count,count_all; adjust_dof_fix(); + double natoms = group->count(igroup); // 6 or 3 dof for extended/point particles for 3d // 3 or 2 dof for extended/point particles 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 *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; count = 0; if (domain->dimension == 3) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (radius[i] == 0.0) { if (mode == ALL) count += 3; } else { if (mode == ALL) count += 6; else count += 3; } } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (radius[i] == 0.0) { if (mode == ALL) count += 2; } else { if (mode == ALL) count += 3; else count += 1; } } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof = count_all; // additional adjustments to dof if (tempbias == 1) { - if (mode == ALL) { - double natoms = group->count(igroup); - dof -= tbias->dof_remove(-1) * natoms; - } + if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms; } else if (tempbias == 2) { int *mask = atom->mask; int nlocal = atom->nlocal; tbias->dof_remove_pre(); count = 0; if (domain->dimension == 3) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (tbias->dof_remove(i)) { if (radius[i] == 0.0) { if (mode == ALL) count += 3; } else { if (mode == ALL) count += 6; else count += 3; } } } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (tbias->dof_remove(i)) { if (radius[i] == 0.0) { if (mode == ALL) count += 2; } else { if (mode == ALL) count += 3; else count += 1; } } } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof -= count_all; } dof -= extra_dof + fix_dof; + if (dof < 0.0 && natoms > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempSphere::compute_scalar() { invoked_scalar = update->ntimestep; if (tempbias) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } double **v = atom->v; double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; // point particles will not contribute rotation due to radius = 0 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]; t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i]; } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic || tempbias == 2) dof_compute(); - if (tfactor == 0.0 && atom->natoms != 0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::compute_vector() { invoked_vector = update->ntimestep; if (tempbias) { if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } double **v = atom->v; double **omega = atom->omega; double *rmass = atom->rmass; double *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; // point particles will not contribute rotation due to radius = 0 double massone,inertiaone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; if (mode == ALL) { for (int 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]; inertiaone = INERTIA*rmass[i]*radius[i]*radius[i]; t[0] += inertiaone * omega[i][0]*omega[i][0]; t[1] += inertiaone * omega[i][1]*omega[i][1]; t[2] += inertiaone * omega[i][2]*omega[i][2]; t[3] += inertiaone * omega[i][0]*omega[i][1]; t[4] += inertiaone * omega[i][0]*omega[i][2]; t[5] += inertiaone * omega[i][1]*omega[i][2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { inertiaone = INERTIA*rmass[i]*radius[i]*radius[i]; t[0] += inertiaone * omega[i][0]*omega[i][0]; t[1] += inertiaone * omega[i][1]*omega[i][1]; t[2] += inertiaone * omega[i][2]*omega[i][2]; t[3] += inertiaone * omega[i][0]*omega[i][1]; t[4] += inertiaone * omega[i][0]*omega[i][2]; t[5] += inertiaone * omega[i][1]*omega[i][2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempSphere::remove_bias(int i, double *v) { tbias->remove_bias(i,v); } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempSphere::restore_bias(int i, double *v) { tbias->restore_bias(i,v); }