diff --git a/src/atom.cpp b/src/atom.cpp index 8cfb2fd65b..7cc2115625 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1036,6 +1036,8 @@ void Atom::set_mass(const char *str) mass[itype] = mass_one; mass_setflag[itype] = 1; + + if (mass[itype] <= 0.0) error->all("Invalid mass value"); } /* ---------------------------------------------------------------------- @@ -1050,6 +1052,8 @@ void Atom::set_mass(int itype, double value) mass[itype] = value; mass_setflag[itype] = 1; + + if (mass[itype] <= 0.0) error->all("Invalid mass value"); } /* ---------------------------------------------------------------------- @@ -1068,6 +1072,8 @@ void Atom::set_mass(int narg, char **arg) for (int itype = lo; itype <= hi; itype++) { mass[itype] = atof(arg[1]); mass_setflag[itype] = 1; + + if (mass[itype] <= 0.0) error->all("Invalid mass value"); } } @@ -1116,6 +1122,9 @@ void Atom::set_shape(const char *str) shape[itype][1] = 0.5*b; shape[itype][2] = 0.5*c; shape_setflag[itype] = 1; + + if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || shape[itype][2] < 0.0) + error->all("Invalid shape value"); } /* ---------------------------------------------------------------------- @@ -1139,6 +1148,10 @@ void Atom::set_shape(int narg, char **arg) shape[itype][1] = 0.5*atof(arg[2]); shape[itype][2] = 0.5*atof(arg[3]); shape_setflag[itype] = 1; + + if (shape[itype][0] < 0.0 || shape[itype][1] < 0.0 || + shape[itype][2] < 0.0) + error->all("Invalid shape value"); } } @@ -1176,13 +1189,15 @@ void Atom::set_dipole(const char *str) { if (dipole == NULL) error->all("Cannot set dipole for this atom style"); - int i; + int itype; double dipole_one; - int n = sscanf(str,"%d %lg",&i,&dipole_one); - if (n != 2) error->all("Invalid shape line in data file"); + int n = sscanf(str,"%d %lg",&itype,&dipole_one); + if (n != 2) error->all("Invalid dipole line in data file"); - dipole[i] = dipole_one; - dipole_setflag[i] = 1; + dipole[itype] = dipole_one; + dipole_setflag[itype] = 1; + + if (dipole[itype] < 0.0) error->all("Invalid dipole value"); } /* ---------------------------------------------------------------------- @@ -1201,6 +1216,8 @@ void Atom::set_dipole(int narg, char **arg) for (int itype = lo; itype <= hi; itype++) { dipole[itype] = atof(arg[1]); dipole_setflag[itype] = 1; + + if (dipole[itype] < 0.0) error->all("Invalid dipole value"); } } diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 4c5aa6d9be..98b0e1591f 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -66,7 +66,7 @@ void ComputeERotateSphere::init() for (int i = 1; i <= atom->ntypes; i++) { if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) error->all("Compute erotate/sphere requires spherical particle shapes"); - inertia[i] = 0.25*shape[i][0]*shape[i][0] * mass[i]; + inertia[i] = shape[i][0]*shape[i][0] * mass[i]; } } } diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 0eb00ef45f..cbd5f1bd08 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -98,7 +98,7 @@ void ComputeTempSphere::init() for (int i = 1; i <= atom->ntypes; i++) { if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2]) error->all("Compute temp/sphere requires spherical particle shapes"); - inertia[i] = INERTIA * 0.25*shape[i][0]*shape[i][0] * mass[i]; + inertia[i] = INERTIA * shape[i][0]*shape[i][0] * mass[i]; } } } diff --git a/src/verlet.cpp b/src/verlet.cpp index 881f5e499f..d271a83cdf 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -221,26 +221,70 @@ void Verlet::force_clear() { int i; - // clear global force array - // nall includes ghosts only if either newton flag is set + // clear force on all particles + // if either newton flag is set, also include ghosts - int nall; - if (force->newton) nall = atom->nlocal + atom->nghost; - else nall = atom->nlocal; + if (neighbor->includegroup == 0) { + int nall; + if (force->newton) nall = atom->nlocal + atom->nghost; + else nall = atom->nlocal; - double **f = atom->f; - for (i = 0; i < nall; i++) { - f[i][0] = 0.0; - f[i][1] = 0.0; - f[i][2] = 0.0; - } - - if (torqueflag) { - double **torque = atom->torque; + double **f = atom->f; for (i = 0; i < nall; i++) { - torque[i][0] = 0.0; - torque[i][1] = 0.0; - torque[i][2] = 0.0; + f[i][0] = 0.0; + f[i][1] = 0.0; + f[i][2] = 0.0; + } + + if (torqueflag) { + double **torque = atom->torque; + for (i = 0; i < nall; i++) { + torque[i][0] = 0.0; + torque[i][1] = 0.0; + torque[i][2] = 0.0; + } + } + + // neighbor includegroup flag is set + // clear force only on initial nfirst particles + // if either newton flag is set, also include ghosts + + } else { + int nall = atom->nfirst; + + double **f = atom->f; + for (i = 0; i < nall; i++) { + f[i][0] = 0.0; + f[i][1] = 0.0; + f[i][2] = 0.0; + } + + if (torqueflag) { + double **torque = atom->torque; + for (i = 0; i < nall; i++) { + torque[i][0] = 0.0; + torque[i][1] = 0.0; + torque[i][2] = 0.0; + } + } + + if (force->newton) { + nall = atom->nlocal + atom->nghost; + + for (i = atom->nlocal; i < nall; i++) { + f[i][0] = 0.0; + f[i][1] = 0.0; + f[i][2] = 0.0; + } + + if (torqueflag) { + double **torque = atom->torque; + for (i = atom->nlocal; i < nall; i++) { + torque[i][0] = 0.0; + torque[i][1] = 0.0; + torque[i][2] = 0.0; + } + } } } }