From 5d66783527634495925ddff422044627f896a8a1 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 1 Dec 2010 00:01:42 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5343 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/atom.cpp | 9 ++- src/atom.h | 3 + src/compute_erotate_sphere.cpp | 14 +++- src/compute_property_atom.cpp | 131 +++++++++++++++++++++++++++++++++ src/compute_property_atom.h | 6 ++ src/input.cpp | 27 +++++-- 6 files changed, 181 insertions(+), 9 deletions(-) diff --git a/src/atom.cpp b/src/atom.cpp index 16f21748f0..f9145e2ebc 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -76,6 +76,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) spin = NULL; eradius = ervel = erforce = NULL; + length = theta = NULL; + maxspecial = 1; nspecial = NULL; special = NULL; @@ -102,6 +104,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) quat_flag = omega_flag = angmom_flag = torque_flag = 0; radius_flag = density_flag = rmass_flag = vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; + line_flag = 0; // ntype-length arrays @@ -185,6 +188,9 @@ Atom::~Atom() memory->sfree(ervel); memory->sfree(erforce); + memory->sfree(length); + memory->sfree(theta); + memory->sfree(molecule); memory->destroy_2d_int_array(nspecial); @@ -262,7 +268,8 @@ void Atom::create_avec(const char *style, int narg, char **arg) quat_flag = omega_flag = angmom_flag = torque_flag = 0; radius_flag = density_flag = rmass_flag = vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; - + line_flag = 0; + avec = new_avec(style,narg,arg); int n = strlen(style) + 1; atom_style = new char[n]; diff --git a/src/atom.h b/src/atom.h index b12b337062..41706d21a2 100644 --- a/src/atom.h +++ b/src/atom.h @@ -55,6 +55,8 @@ class Atom : protected Pointers { int *spin; double *eradius,*ervel,*erforce; + double *length,*theta; + int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs int **special; // IDs of 1-2,1-3,1-4 neighs of each atom int maxspecial; // special[nlocal][maxspecial] @@ -84,6 +86,7 @@ class Atom : protected Pointers { int quat_flag,omega_flag,angmom_flag,torque_flag; int radius_flag,density_flag,rmass_flag,vfrac_flag; int spin_flag,eradius_flag,ervel_flag,erforce_flag; + int line_flag; // extra peratom info in restart file destined for fix & diag diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 1558196cce..27e121bcdc 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -23,7 +23,8 @@ using namespace LAMMPS_NS; -#define INERTIA 0.4 // moment of inertia for sphere +//#define INERTIA 0.4 // moment of inertia for sphere +#define INERTIA (1.0/12.0) // moment of inertia for sphere /* ---------------------------------------------------------------------- */ @@ -37,11 +38,13 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) : // error checks + /* if (!atom->omega_flag) error->all("Compute erotate/sphere requires atom attribute omega"); if (!atom->radius_flag && !atom->avec->shape_type) error->all("Compute erotate/sphere requires atom attribute " "radius or shape"); + */ } /* ---------------------------------------------------------------------- */ @@ -53,6 +56,7 @@ void ComputeERotateSphere::init() // if shape used, check that all particles are spherical // point particles are allowed + /* if (atom->radius == NULL) { double **shape = atom->shape; int *type = atom->type; @@ -68,6 +72,7 @@ void ComputeERotateSphere::init() "spherical particle shapes"); } } + */ pfactor = 0.5 * force->mvv2e * INERTIA; } @@ -112,6 +117,13 @@ double ComputeERotateSphere::compute_scalar() } else { if (rmass) { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + + omega[i][2]*omega[i][2]) * + atom->length[i]*atom->length[i]*rmass[i]; + } + } else if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { itype = type[i]; diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index 4447d9f06f..926e5d2e6f 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "math.h" #include "string.h" #include "compute_property_atom.h" #include "atom.h" @@ -218,6 +219,37 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : "atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_erforce; + } else if (strcmp(arg[iarg],"end1x") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end1x; + } else if (strcmp(arg[iarg],"end1y") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end1y; + } else if (strcmp(arg[iarg],"end1z") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end1z; + } else if (strcmp(arg[iarg],"end2x") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end2x; + } else if (strcmp(arg[iarg],"end2y") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end2y; + } else if (strcmp(arg[iarg],"end2z") == 0) { + if (!atom->line_flag) + error->all("Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_end2z; + } else error->all("Invalid keyword in compute property/atom command"); } @@ -1100,3 +1132,102 @@ void ComputePropertyAtom::pack_erforce(int n) n += nvalues; } } + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end1x(int n) +{ + double **x = atom->x; + double *length = atom->length; + double *theta = atom->theta; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][0] - 0.5*length[i]*cos(theta[i]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end1y(int n) +{ + double **x = atom->x; + double *length = atom->length; + double *theta = atom->theta; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][1] - 0.5*length[i]*sin(theta[i]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end1z(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end2x(int n) +{ + double **x = atom->x; + double *length = atom->length; + double *theta = atom->theta; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][0] + 0.5*length[i]*cos(theta[i]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end2y(int n) +{ + double **x = atom->x; + double *length = atom->length; + double *theta = atom->theta; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][1] + 0.5*length[i]*sin(theta[i]); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputePropertyAtom::pack_end2z(int n) +{ + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = x[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h index 1d8cd07f23..196a22b735 100644 --- a/src/compute_property_atom.h +++ b/src/compute_property_atom.h @@ -94,6 +94,12 @@ class ComputePropertyAtom : public Compute { void pack_eradius(int); void pack_ervel(int); void pack_erforce(int); + void pack_end1x(int); + void pack_end1y(int); + void pack_end1z(int); + void pack_end2x(int); + void pack_end2y(int); + void pack_end2z(int); }; } diff --git a/src/input.cpp b/src/input.cpp index 6fd8e99f66..46eb4c57d5 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -333,13 +333,14 @@ void Input::parse() } /* ---------------------------------------------------------------------- - substitute for $ variables in str + substitute for $ variables in str and return it + str assumed to be long enough to hold expanded version print updated string if flag is set and not searching for label ------------------------------------------------------------------------- */ void Input::substitute(char *str, int flag) { - // use work[] as scratch space to expand str + // use work[] as scratch space to expand str, then copy back to str // do not replace $ inside single/double quotes // var = pts at variable name, ended by NULL // if $ is followed by '{', trailing '}' becomes NULL @@ -528,12 +529,15 @@ void Input::ifthenelse() // substitute for variables in Boolean expression for "if" // in case expression was enclosed in quotes + // must substitute on copy of arg else will step on subsequent args - substitute(arg[0],0); + char *scopy = new char[MAXLINE]; + strcpy(scopy,arg[0]); + substitute(scopy,0); // evaluate Boolean expression for "if" - double btest = variable->evaluate_boolean(arg[0]); + double btest = variable->evaluate_boolean(scopy); // bound "then" commands @@ -546,6 +550,8 @@ void Input::ifthenelse() iarg++; int last = iarg-1; + printf("FIRSTLAST %d %d\n",first,last); + // execute "then" commands // make copies of all arg string commands // required because re-parsing a command via one() will wipe out args @@ -569,23 +575,29 @@ void Input::ifthenelse() for (int i = 0; i < ncommands; i++) delete [] commands[i]; delete [] commands; + delete [] scopy; return; } // done if no "elif" or "else" - if (iarg == narg) return; + if (iarg == narg) { + delete [] scopy; + return; + } // check "elif" or "else" until find commands to execute // substitute for variables and evaluate Boolean expression for "elif" + // must substitute on copy of arg else will step on subsequent args // bound and execute "elif" or "else" commands while (1) { if (iarg+2 > narg) error->all("Illegal if command"); if (strcmp(arg[iarg],"elif") == 0) { - substitute(arg[iarg+1],0); - btest = variable->evaluate_boolean(arg[iarg+1]); + strcpy(scopy,arg[iarg+1]); + substitute(scopy,0); + btest = variable->evaluate_boolean(scopy); first = iarg+2; } else { btest = 1.0; @@ -622,6 +634,7 @@ void Input::ifthenelse() for (int i = 0; i < ncommands; i++) delete [] commands[i]; delete [] commands; + delete [] scopy; return; }