From 8f02ff3b29a0811d75591ef8b8b89e81f58b4dd7 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 6 Sep 2010 16:19:47 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4690 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MANYBODY/pair_comb.cpp | 4 +- src/PERI/fix_peri_neigh.cpp | 2 - src/PERI/pair_peri_lps.cpp | 75 ++++++++++++++++++++++--------------- src/compute_temp_region.cpp | 2 +- src/fix_adapt.cpp | 4 +- src/fix_box_relax.cpp | 4 +- src/fix_indent.cpp | 2 +- src/fix_nh.cpp | 2 +- src/fix_nh_sphere.cpp | 2 +- src/fix_nve_sphere.cpp | 2 +- src/fix_press_berendsen.cpp | 18 ++++----- src/fix_rigid.cpp | 2 +- src/fix_wall_reflect.cpp | 12 +++--- 13 files changed, 71 insertions(+), 60 deletions(-) diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index 0dbcc648d6..83eebfae7e 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -691,7 +691,7 @@ void PairComb::read_file(char *file) params[nparams].powermint = int(params[nparams].powerm); - // parameter set sanity checks + // parameter sanity checks if (params[nparams].lam11 < 0.0 || params[nparams].lam12 < 0.0 || params[nparams].c < 0.0 || params[nparams].d < 0.0 || @@ -723,7 +723,7 @@ void PairComb::read_file(char *file) params[nparams].lam12 < params[nparams].lam22 || params[nparams].biga1< params[nparams].bigb1 || params[nparams].biga2< params[nparams].bigb2) - error->all("Improper COMB pair parameters"); + error->all("Illegal COMB parameter"); nparams++; } diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index 3744baa7f0..14dbf8cbba 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -281,8 +281,6 @@ void FixPeriNeigh::setup(int vflag) else if (pairlps != NULL) // call the PairPeriLPS influence function wvolume[i] += pairlps->influence_function(delx0,dely0,delz0) * rsq0 * vfrac[j] * vfrac_scale; - else - error->all("Unknown peridynamic pair style in FixPeriNeigh."); } } diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp index 69aa232254..1a3173134d 100644 --- a/src/PERI/pair_peri_lps.cpp +++ b/src/PERI/pair_peri_lps.cpp @@ -59,7 +59,8 @@ PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp) // set comm size needed by this Pair // comm_reverse not needed - comm_forward = 1; // For passing dilatation (theta) + + comm_forward = 1; // for passing dilatation (theta) } /* ---------------------------------------------------------------------- */ @@ -287,10 +288,14 @@ void PairPeriLPS::compute(int eflag, int vflag) (1.0 + ((delta - half_lc)/(2*half_lc) ) ); else vfrac_scale = 1.0; - omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); - omega_minus = influence_function(delx0,dely0,delz0); - rk = ( (3.0 * bulkmodulus[itype][itype]) - (5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ; - rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_minus = influence_function(delx0,dely0,delz0); + rk = ( (3.0 * bulkmodulus[itype][itype]) - + (5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale * + ( (omega_plus * theta[i] / wvolume[i]) + + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj]; + rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * + ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; if (r > 0.0) fbond = -(rk/r); else fbond = 0.0; @@ -302,7 +307,9 @@ void PairPeriLPS::compute(int eflag, int vflag) // since I-J is double counted, set newton off & use 1/2 factor and I,I double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0); - if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * omega_plus * ( deviatoric_extension * deviatoric_extension ) * vfrac[j] * vfrac_scale; + if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus*(deviatoric_extension * deviatoric_extension) * + vfrac[j] * vfrac_scale; if (evflag) ev_tally(i,i,nlocal,0, 0.5*evdwl,0.0,0.5*fbond,delx,dely,delz); @@ -346,8 +353,8 @@ void PairPeriLPS::allocate() setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); - bulkmodulus = memory->create_2d_double_array(n+1,n+1,"pair:bulkmodulus"); - shearmodulus = memory->create_2d_double_array(n+1,n+1,"pair:shearmodulus"); + bulkmodulus = memory->create_2d_double_array(n+1,n+1,"pair:bulkmodulus"); + shearmodulus = memory->create_2d_double_array(n+1,n+1,"pair:shearmodulus"); s00 = memory->create_2d_double_array(n+1,n+1,"pair:s00"); alpha = memory->create_2d_double_array(n+1,n+1,"pair:alpha"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); @@ -514,8 +521,8 @@ void PairPeriLPS::read_restart(FILE *fp) /* ---------------------------------------------------------------------- */ -double PairPeriLPS::single(int i, int j, int itype, int jtype, double rsq, - double factor_coul, double factor_lj, +double PairPeriLPS::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { double delx0,dely0,delz0,rsq0; @@ -547,10 +554,10 @@ double PairPeriLPS::single(int i, int j, int itype, int jtype, double rsq, if (r < d_ij) { dr = r - d_ij; - // kshort resembles the short-range force constant of the bond-based theory in 3-D ! - kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / - ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype] ); - rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype])); + // kshort resembles short-range force constant of bond-based theory in 3d + kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / + ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]); + rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype])); if (r > 0.0) fforce += -(rk/r); energy += 0.5*rk*dr; } @@ -587,13 +594,19 @@ double PairPeriLPS::single(int i, int j, int itype, int jtype, double rsq, (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc))); else vfrac_scale = 1.0; - omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); omega_minus = influence_function(delx0,dely0,delz0); - rk = (3.0* bulkmodulus[itype][itype] -5.0 * shearmodulus[itype][itype] )* vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + ( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj] ; - rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; + rk = (3.0* bulkmodulus[itype][itype] -5.0 * shearmodulus[itype][itype]) * + vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + + (omega_minus * theta[j] / wvolume[j])) * + r0[i][jj]; + rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * + ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; if (r > 0.0) fforce += -(rk/r); - energy += 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * omega_plus * ( dr - theta[i]* r0[i][jj] / 3.0 ) * ( dr - theta[i]* r0[i][jj] / 3.0 ) * vfrac[j] * vfrac_scale; + energy += 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus * ( dr - theta[i]* r0[i][jj] / 3.0 ) * + ( dr - theta[i]* r0[i][jj] / 3.0 ) * vfrac[j] * vfrac_scale; } } @@ -612,7 +625,7 @@ double PairPeriLPS::memory_usage() } /* ---------------------------------------------------------------------- - Influence function definition + influence function definition ------------------------------------------------------------------------- */ double PairPeriLPS::influence_function(double xi_x, double xi_y, double xi_z) @@ -620,10 +633,9 @@ double PairPeriLPS::influence_function(double xi_x, double xi_y, double xi_z) double r = sqrt(xi_x*xi_x + xi_y*xi_y + xi_z*xi_z); double omega; - if (fabs(r) < 2.2204e-016) error->one("pair_peri_lps - influence_function: Divide by 0"); - + if (fabs(r) < 2.2204e-016) + error->one("Divide by 0 in influence function of pair peri/lps"); omega = 1.0/r; - return omega; } @@ -655,7 +667,7 @@ void PairPeriLPS::compute_dilatation() int nonperiodic = domain->nonperiodic; // compute the dilatation theta - // loop over my particles + for (i = 0; i < nlocal; i++) { xtmp = x[i][0]; ytmp = x[i][1]; @@ -667,7 +679,6 @@ void PairPeriLPS::compute_dilatation() theta[i] = 0.0; itype = type[i]; - // loop over partners of particle i for (jj = 0; jj < jnum; jj++) { // if bond already broken, skip this partner @@ -704,23 +715,25 @@ void PairPeriLPS::compute_dilatation() (1.0 + ((delta - half_lc)/(2*half_lc) ) ); else vfrac_scale = 1.0; - theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * vfrac[j] * vfrac_scale; + theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * + vfrac[j] * vfrac_scale; - } // end loop over all neighbors jj + } + + // if wvolume[i] is zero, then particle i has no bonds + // therefore, the dilatation is set to - // If wvolume[i] is zero, then particle i has no bonds. Therefore, the dilatation is set to 0. if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i]; else theta[i] = 0; - - } // end loop over all particles i - + } } /* ---------------------------------------------------------------------- communication routines ---------------------------------------------------------------------- */ -int PairPeriLPS::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +int PairPeriLPS::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) { int i,j,m; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index f1d43f83e5..f8a7adec76 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -68,7 +68,7 @@ void ComputeTempRegion::init() iregion = domain->find_region(idregion); if (iregion == -1) - error->all("Region ID for temp reduce/region does not exist"); + error->all("Region ID for compute temp reduce/region does not exist"); dof = 0.0; } diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 83ffe62cec..d874a1e173 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -164,14 +164,14 @@ void FixAdapt::init() if (strcmp(param[m],"diameter") == 0) { awhich[m] = DIAMETER; if (!atom->radius_flag) - error->all("Fix adapt requires atom attribute radius"); + error->all("Fix adapt requires atom attribute diameter"); } else error->all("Fix adapt atom attribute is not recognized"); } ivar[m] = input->variable->find(var[m]); if (ivar[m] < 0) error->all("Variable name for fix adapt does not exist"); if (!input->variable->equalstyle(ivar[m])) - error->all("Variable for fix adapt is not equal style"); + error->all("Variable for fix adapt is invalid style"); } } diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp index 99e3e44f16..f10e87f312 100644 --- a/src/fix_box_relax.cpp +++ b/src/fix_box_relax.cpp @@ -219,7 +219,7 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : if (pcouple == XYZ && dimension == 3 && (p_target[0] != p_target[1] || p_target[0] != p_target[2])) - error->all("Invalid fix nvt/npt/nph pressure settings"); + error->all("Invalid fix box/relax pressure settings"); if (pcouple == XYZ && dimension == 2 && p_target[0] != p_target[1]) error->all("Invalid fix box/relax pressure settings"); if (pcouple == XY && p_target[0] != p_target[1]) @@ -611,7 +611,7 @@ void FixBoxRelax::remap() domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0-ctr)*ds[i]/s0[i]; domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0-ctr)*ds[i]/s0[i]; if (domain->boxlo[i] >= domain->boxhi[i]) - error->all("fix box/relax generated negative dimension"); + error->all("Fix box/relax generated negative box length"); } if (pstyle == TRICLINIC) { diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 55f5f9cf86..8bd4ac2a2c 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -116,7 +116,7 @@ void FixIndent::init() xvar = input->variable->find(xstr); if (xvar < 0) error->all("Variable name for fix indent does not exist"); if (!input->variable->equalstyle(xvar)) - error->all("Variable for fix indent is not equal style"); + error->all("Variable for fix indent is invalid style"); } if (ystr) { yvar = input->variable->find(ystr); diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 35803f4fbe..da9e470ec5 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -96,7 +96,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) t_stop = atof(arg[iarg+2]); t_period = atof(arg[iarg+3]); if (t_start < 0.0 || t_stop <= 0.0) - error->all("Target T for fix nvt/npt/nph cannot be 0.0"); + error->all("Target temperature for fix nvt/npt/nph cannot be 0.0"); iarg += 4; } else if (strcmp(arg[iarg],"iso") == 0) { diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp index 956fbf8c58..6dbf9566e5 100644 --- a/src/fix_nh_sphere.cpp +++ b/src/fix_nh_sphere.cpp @@ -36,7 +36,7 @@ FixNHSphere::FixNHSphere(LAMMPS *lmp, int narg, char **arg) : "atom attributes omega, torque"); if (!atom->radius_flag && !atom->avec->shape_type) error->all("Fix nvt/nph/npt sphere requires " - "atom attribute radius or shape"); + "atom attribute diameter or shape"); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index 7a8f66934a..7c21b63804 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -56,7 +56,7 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) : if (!atom->omega_flag || !atom->torque_flag) error->all("Fix nve/sphere requires atom attributes omega, torque"); if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Fix nve/sphere requires atom attribute radius or shape"); + error->all("Fix nve/sphere requires atom attribute diameter or shape"); if (extra == DIPOLE && !atom->mu_flag) error->all("Fix nve/sphere requires atom attribute mu"); } diff --git a/src/fix_press_berendsen.cpp b/src/fix_press_berendsen.cpp index a4c05d96c8..2da7bde7b7 100644 --- a/src/fix_press_berendsen.cpp +++ b/src/fix_press_berendsen.cpp @@ -116,7 +116,7 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : p_flag[2] = 1; iarg += 4; if (dimension == 2) - error->all("Invalid fix press/berendsen command for a 2d simulation"); + error->all("Invalid fix press/berendsen for a 2d simulation"); } else if (strcmp(arg[iarg],"couple") == 0) { if (iarg+2 > narg) error->all("Illegal fix press/berendsen command"); @@ -148,20 +148,20 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : // error checks if (dimension == 2 && p_flag[2]) - error->all("Invalid fix press/berendsen command for a 2d simulation"); + error->all("Invalid fix press/berendsen for a 2d simulation"); if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) - error->all("Invalid fix press/berendsen command for a 2d simulation"); + error->all("Invalid fix press/berendsen for a 2d simulation"); if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) - error->all("Invalid fix press/berendsen command pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (p_flag[0] && domain->xperiodic == 0) error->all("Cannot use fix press/berendsen on a non-periodic dimension"); @@ -178,7 +178,7 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : if (pcouple == XYZ && dimension == 2 && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) - error->all("Invalid fix press_berendsen pressure settings"); + error->all("Invalid fix press/berendsen pressure settings"); if (pcouple == XY && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 051cdb46cd..92c610b5e9 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -1332,7 +1332,7 @@ void FixRigid::no_squish_rotate(int k, double *p, double *q, kq[1] = q[2]; kp[1] = p[2]; kq[2] = -q[1]; kp[2] = -p[1]; kq[3] = q[0]; kp[3] = p[0]; - } else error->all("Invalid no-squish rotate pivot value"); + } // obtain phi, cosines and sines diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index 26d35980e6..ba9aea7b4a 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -210,42 +210,42 @@ void FixWallReflect::init() if (xlovar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(xlovar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } if (xhistr) { xhivar = input->variable->find(xhistr); if (xhivar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(xhivar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } if (ylostr) { ylovar = input->variable->find(ylostr); if (ylovar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(ylovar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } if (yhistr) { yhivar = input->variable->find(yhistr); if (yhivar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(yhivar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } if (zlostr) { zlovar = input->variable->find(zlostr); if (zlovar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(zlovar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } if (zhistr) { zhivar = input->variable->find(zhistr); if (zhivar < 0) error->all("Variable name for fix wall/relect does not exist"); if (!input->variable->equalstyle(zhivar)) - error->all("Variable for fix wall/reflect is not equal style"); + error->all("Variable for fix wall/reflect is invalid style"); } int nrigid = 0;