From 16ef076c318d1d2965b9a4d31f2065f468fb777e Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 29 Apr 2009 16:54:06 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2780 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute.h | 4 + src/compute_pressure.cpp | 4 +- src/compute_temp_partial.cpp | 47 ++-- src/compute_temp_partial.h | 3 - src/compute_temp_profile.cpp | 494 +++++++++++++++++++++++++++++++++++ src/compute_temp_profile.h | 59 +++++ src/compute_temp_ramp.cpp | 16 +- src/compute_temp_ramp.h | 3 - src/compute_temp_region.h | 3 - src/group.cpp | 46 ++-- src/style.h | 2 + 11 files changed, 627 insertions(+), 54 deletions(-) create mode 100644 src/compute_temp_profile.cpp create mode 100644 src/compute_temp_profile.h diff --git a/src/compute.h b/src/compute.h index 233e61f9ea..cc4a24fb17 100644 --- a/src/compute.h +++ b/src/compute.h @@ -96,6 +96,10 @@ class Compute : protected Pointers { int extra_dof; // extra DOF for temperature computes int dynamic; // recount atoms for temperature computes int thermoflag; // 1 if include fix PE for PE computes + + double vbias[3]; // stored velocity bias for one atom + double **vbiasall; // stored velocity bias for all atoms + int maxbias; // size of vbiasall array }; } diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp index 21f8b68499..2e9d4157ab 100644 --- a/src/compute_pressure.cpp +++ b/src/compute_pressure.cpp @@ -110,8 +110,8 @@ void ComputePressure::init() nktv2p = force->nktv2p; dimension = domain->dimension; - // set temperature used by pressure - // must be done in init() since user can change it via modify command + // set temperature compute, must be done in init() + // fixes could have changed or compute_modify could have changed it int icompute = modify->find_compute(id_pre); if (icompute < 0) error->all("Could not find compute pressure temp ID"); diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 088aca39ab..faae430bee 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -33,10 +33,6 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : { if (narg != 6) error->all("Illegal compute temp/partial command"); - xflag = atoi(arg[3]); - yflag = atoi(arg[4]); - zflag = atoi(arg[5]); - scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; @@ -44,6 +40,10 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : tempflag = 1; tempbias = 1; + xflag = atoi(arg[3]); + yflag = atoi(arg[4]); + zflag = atoi(arg[5]); + maxbias = 0; vbiasall = NULL; vector = new double[6]; @@ -193,21 +193,27 @@ void ComputeTempPartial::remove_bias_all() "compute/temp:vbiasall"); } - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (!xflag) { + if (!xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][0] = v[i][0]; v[i][0] = 0.0; } - if (!yflag) { + } + if (!yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][1] = v[i][1]; v[i][1] = 0.0; } - if (!zflag) { + } + if (!zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { vbiasall[i][2] = v[i][2]; v[i][2] = 0.0; } - } + } } /* ---------------------------------------------------------------------- @@ -233,12 +239,21 @@ void ComputeTempPartial::restore_bias_all() int *mask = atom->mask; int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (!xflag) v[i][0] += vbiasall[i][0]; - if (!yflag) v[i][1] += vbiasall[i][1]; - if (!zflag) v[i][2] += vbiasall[i][2]; - } + if (!xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][0] += vbiasall[i][0]; + } + if (!yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][1] += vbiasall[i][1]; + } + if (!zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][2] += vbiasall[i][2]; + } } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_temp_partial.h b/src/compute_temp_partial.h index 0505265ed8..9d414f2221 100644 --- a/src/compute_temp_partial.h +++ b/src/compute_temp_partial.h @@ -37,9 +37,6 @@ class ComputeTempPartial : public Compute { int xflag,yflag,zflag; int fix_dof; double tfactor; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array void dof_compute(); }; diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp new file mode 100644 index 0000000000..8733fe6b29 --- /dev/null +++ b/src/compute_temp_profile.cpp @@ -0,0 +1,494 @@ +/* ---------------------------------------------------------------------- + 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 "modify.h" +#include "fix.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 7) error->all("Illegal compute temp/profile command"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = 1; + tempflag = 1; + tempbias = 1; + + xflag = atoi(arg[3]); + yflag = atoi(arg[4]); + zflag = atoi(arg[5]); + + nbinx = nbiny = nbinz = 1; + + if (strcmp(arg[6],"x") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + nbinx = atoi(arg[7]); + ivx = 1; + ncount = 1; + } else if (strcmp(arg[6],"y") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + nbiny = atoi(arg[7]); + ivy = 0; + ncount = 1; + } else if (strcmp(arg[6],"z") == 0) { + if (narg != 8) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinz = atoi(arg[7]); + ivz = 0; + ncount = 1; + } else if (strcmp(arg[6],"xy") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + nbinx = atoi(arg[7]); + nbiny = atoi(arg[8]); + ivx = 0; + ivy = 1; + ncount = 2; + } else if (strcmp(arg[6],"yz") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbiny = atoi(arg[7]); + nbinz = atoi(arg[8]); + ivy = 0; + ivz = 1; + ncount = 2; + } else if (strcmp(arg[6],"xz") == 0) { + if (narg != 9) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinx = atoi(arg[7]); + nbinz = atoi(arg[8]); + ivx = 0; + ivz = 1; + ncount = 2; + } else if (strcmp(arg[6],"xyz") == 0) { + if (narg != 10) error->all("Illegal compute temp/profile command"); + if (domain->dimension == 2) + error->all("Compute temp/profile command cannot bin z for 2d systems"); + nbinx = atoi(arg[7]); + nbiny = atoi(arg[8]); + nbinz = atoi(arg[9]); + ivx = 0; + ivy = 1; + ivz = 2; + ncount = 3; + } else error->all("Illegal compute temp/profile command"); + + nbins = nbinx*nbiny*nbinz; + if (nbins <= 0) error->all("Illegal compute temp/profile command"); + + vbin = memory->create_2d_double_array(nbins,ncount+1,"temp/profile:vbin"); + binave = memory->create_2d_double_array(nbins,ncount+1, + "temp/profile:binave"); + + maxatom = 0; + bin = NULL; + + maxbias = 0; + vbiasall = NULL; + vector = new double[6]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempProfile::~ComputeTempProfile() +{ + memory->destroy_2d_double_array(vbin); + memory->destroy_2d_double_array(binave); + memory->sfree(bin); + memory->destroy_2d_double_array(vbiasall); + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempProfile::init() +{ + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + 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::dof_compute() +{ + double natoms = group->count(igroup); + int nper = domain->dimension; + dof = nper * natoms; + dof -= extra_dof + fix_dof; + 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 **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) { + 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(); + scalar *= tfactor; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempProfile::compute_vector() +{ + int i,ibin; + double vthermal[3]; + + invoked_vector = update->ntimestep; + + bin_average(); + + 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) { + 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; +} + +/* ---------------------------------------------------------------------- + remove velocity bias from atom I to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempProfile::remove_bias(int i, double *v) +{ + int ibin = bin[i]; + if (xflag) { + vbias[0] = binave[ibin][ivx]; + v[0] -= vbias[0]; + } + if (yflag) { + vbias[1] = binave[ibin][ivy]; + v[1] -= vbias[1]; + } + if (zflag) { + vbias[2] = binave[ibin][ivz]; + v[2] -= vbias[2]; + } +} + +/* ---------------------------------------------------------------------- + remove velocity bias from all atoms to leave thermal velocity +------------------------------------------------------------------------- */ + +void ComputeTempProfile::remove_bias_all() +{ + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (nlocal > maxbias) { + memory->destroy_2d_double_array(vbiasall); + maxbias = atom->nmax; + vbiasall = memory->create_2d_double_array(maxbias,3, + "compute/temp:vbiasall"); + } + + int ibin; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) { + vbiasall[i][0] = binave[ibin][ivx]; + v[i][0] -= vbiasall[i][0]; + } + if (yflag) { + vbiasall[i][1] = binave[ibin][ivy]; + v[i][1] -= vbiasall[i][1]; + } + if (zflag) { + vbiasall[i][2] = binave[ibin][ivz]; + 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 ComputeTempProfile::restore_bias(int i, double *v) +{ + if (xflag) v[0] += vbias[0]; + if (yflag) v[1] += vbias[1]; + if (zflag) 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 ComputeTempProfile::restore_bias_all() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (xflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][0] += vbiasall[i][0]; + } + if (yflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][1] += vbiasall[i][1]; + } + if (zflag) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + v[i][2] += vbiasall[i][2]; + } +} + +/* ---------------------------------------------------------------------- + compute average velocity in each bin +------------------------------------------------------------------------- */ + +void ComputeTempProfile::bin_average() +{ + int i,j,ibin; + + if (box_change) bin_setup(); + bin_assign(); + + // clear bins, including particle count + + int nc1 = ncount + 1; + for (i = 0; i < nbins; i++) + for (j = 0; j < nc1; j++) + vbin[i][j] = 0.0; + + // sum each particle's velocity to appropriate bin + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + ibin = bin[i]; + if (xflag) vbin[ibin][ivx] += v[i][0]; + if (yflag) vbin[ibin][ivy] += v[i][1]; + if (zflag) vbin[ibin][ivz] += v[i][2]; + vbin[ibin][ncount] += 1.0; + } + + // sum bins across processors + + MPI_Allreduce(vbin[0],binave[0],nbins*ncount,MPI_DOUBLE,MPI_SUM,world); + + // compute ave velocity in each bin, checking for no particles + + for (i = 0; i < nbins; i++) + if (binave[i][ncount] > 0.0) + for (j = 0; j < ncount; j++) + binave[i][j] /= binave[i][ncount]; +} + +/* ---------------------------------------------------------------------- + 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->sfree(bin); + bin = (int *) memory->smalloc(maxatom*sizeof(int),"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 = maxbias * sizeof(double); + bytes += maxatom * sizeof(int); + bytes += nbins*ncount * sizeof(double); + return bytes; +} diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h new file mode 100644 index 0000000000..9194e90c48 --- /dev/null +++ b/src/compute_temp_profile.h @@ -0,0 +1,59 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_TEMP_PROFILE_H +#define COMPUTE_TEMP_PROFILE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTempProfile : public Compute { + public: + ComputeTempProfile(class LAMMPS *, int, char **); + ~ComputeTempProfile(); + void init(); + double compute_scalar(); + void compute_vector(); + + void remove_bias(int, double *); + void remove_bias_all(); + void restore_bias(int, double *); + void restore_bias_all(); + double memory_usage(); + + private: + int xflag,yflag,zflag,ncount; + int nbinx,nbiny,nbinz,nbins; + int ivx,ivy,ivz; + int fix_dof; + double tfactor; + + int box_change,triclinic; + int *periodicity; + double *boxlo,*boxhi,*prd; + double invdelta[3]; + + int maxatom; + int *bin; + double **vbin,**binave; + + void dof_compute(); + void bin_average(); + void bin_setup(); + void bin_assign(); +}; + +} + +#endif diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 5e13d3e823..897b8b8744 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -38,6 +38,13 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : { if (narg < 9) error->all("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; @@ -99,15 +106,6 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : coord_hi = zscale*atof(arg[8]); } - // settings - - scalar_flag = vector_flag = 1; - size_vector = 6; - extscalar = 0; - extvector = 1; - tempflag = 1; - tempbias = 1; - maxbias = 0; vbiasall = NULL; vector = new double[6]; diff --git a/src/compute_temp_ramp.h b/src/compute_temp_ramp.h index 986ee09a34..ae1148836a 100644 --- a/src/compute_temp_ramp.h +++ b/src/compute_temp_ramp.h @@ -39,9 +39,6 @@ class ComputeTempRamp : public Compute { double v_lo,v_hi; int scaleflag,fix_dof; double tfactor,xscale,yscale,zscale; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array void dof_compute(); }; diff --git a/src/compute_temp_region.h b/src/compute_temp_region.h index 7dc5087a64..43ac880997 100644 --- a/src/compute_temp_region.h +++ b/src/compute_temp_region.h @@ -35,9 +35,6 @@ class ComputeTempRegion : public Compute { private: int iregion; - double vbias[3]; // stored velocity bias for one atom - double **vbiasall; // stored velocity bias for all atoms - int maxbias; // size of vbiasall array }; } diff --git a/src/group.cpp b/src/group.cpp index d75de55e9d..2a7a395960 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -774,9 +774,11 @@ void Group::xcm(int igroup, double masstotal, double *cm) } MPI_Allreduce(cmone,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -833,9 +835,11 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion) } MPI_Allreduce(cmone,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -877,9 +881,11 @@ void Group::vcm(int igroup, double masstotal, double *cm) } MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -923,9 +929,11 @@ void Group::vcm(int igroup, double masstotal, double *cm, int iregion) } MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world); - cm[0] /= masstotal; - cm[1] /= masstotal; - cm[2] /= masstotal; + if (masstotal > 0.0) { + cm[0] /= masstotal; + cm[1] /= masstotal; + cm[2] /= masstotal; + } } /* ---------------------------------------------------------------------- @@ -1092,7 +1100,8 @@ double Group::gyration(int igroup, double masstotal, double *cm) double rg_all; MPI_Allreduce(&rg,&rg_all,1,MPI_DOUBLE,MPI_SUM,world); - return sqrt(rg_all/masstotal); + if (masstotal > 0.0) return sqrt(rg_all/masstotal); + return 0.0; } /* ---------------------------------------------------------------------- @@ -1136,7 +1145,8 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion) double rg_all; MPI_Allreduce(&rg,&rg_all,1,MPI_DOUBLE,MPI_SUM,world); - return sqrt(rg_all/masstotal); + if (masstotal > 0.0) return sqrt(rg_all/masstotal); + return 0.0; } /* ---------------------------------------------------------------------- @@ -1264,10 +1274,10 @@ void Group::omega(double *angmom, double inertia[3][3], double *w) inertia[0][1]*inertia[1][0]*inertia[2][2] - inertia[2][0]*inertia[1][1]*inertia[0][2]; - int i,j; - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - inverse[i][j] /= determinant; + if (determinant > 0.0) + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + inverse[i][j] /= determinant; w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] + inverse[0][2]*angmom[2]; diff --git a/src/style.h b/src/style.h index 4af24afe8b..051622abf5 100644 --- a/src/style.h +++ b/src/style.h @@ -94,6 +94,7 @@ CommandStyle(write_restart,WriteRestart) #include "compute_temp_com.h" #include "compute_temp_deform.h" #include "compute_temp_partial.h" +#include "compute_temp_profile.h" #include "compute_temp_ramp.h" #include "compute_temp_region.h" #include "compute_temp_sphere.h" @@ -118,6 +119,7 @@ ComputeStyle(temp,ComputeTemp) ComputeStyle(temp/com,ComputeTempCOM) ComputeStyle(temp/deform,ComputeTempDeform) ComputeStyle(temp/partial,ComputeTempPartial) +ComputeStyle(temp/profile,ComputeTempProfile) ComputeStyle(temp/ramp,ComputeTempRamp) ComputeStyle(temp/region,ComputeTempRegion) ComputeStyle(temp/sphere,ComputeTempSphere)