From 64b10074dc6a8eba50005055100e2344f670e85a Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 21 Feb 2007 00:15:29 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@322 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/KSPACE/ewald.h | 3 +- src/compute.h | 6 +- src/fix.cpp | 2 + src/fix.h | 6 +- src/fix_nph.cpp | 1 + src/fix_npt.cpp | 1 + src/modify.cpp | 51 +++++- src/modify.h | 7 +- src/respa.cpp | 98 ++++++++--- src/respa.h | 15 +- src/style.h | 8 + src/variable.cpp | 427 ++++++++++++++++++++++++++++----------------- src/variable.h | 18 +- src/verlet.h | 14 +- 14 files changed, 455 insertions(+), 202 deletions(-) diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h index e6f5a1a19f..c953e749b0 100644 --- a/src/KSPACE/ewald.h +++ b/src/KSPACE/ewald.h @@ -50,7 +50,8 @@ class Ewald : public KSpace { void slabcorr(int); }; -#endif } +#endif + diff --git a/src/compute.h b/src/compute.h index 49cb3903c1..f0d9222e04 100644 --- a/src/compute.h +++ b/src/compute.h @@ -23,15 +23,15 @@ class Compute : protected Pointers { char *id,*style; int igroup,groupbit; - double scalar; // computed scalar - double *vector; // computed vector + double scalar; // computed global scalar + double *vector; // computed global vector double *scalar_atom; // computed per-atom scalar double **vector_atom; // computed per-atom vector int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int peratom_flag; // 0/1 if compute_peratom() function exists - int size_vector; // N = size of vector + int size_vector; // N = size of global vector int size_peratom; // 0 = just scalar_atom, N = size of vector_atom int extensive; // 0/1 if scalar,vector are intensive/extensive values diff --git a/src/fix.cpp b/src/fix.cpp index 2fb690f1ce..9b0b2cc3c1 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -37,6 +37,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) restart_peratom = 0; force_reneighbor = 0; thermo_energy = 0; + pressure_every = 0; comm_forward = comm_reverse = 0; neigh_half_once = neigh_half_every = 0; @@ -55,6 +56,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) POST_FORCE_RESPA = 256; FINAL_INTEGRATE_RESPA = 512; MIN_POST_FORCE = 1024; + MIN_ENERGY = 2048; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix.h b/src/fix.h index 564eefcc9d..176cc11ec8 100644 --- a/src/fix.h +++ b/src/fix.h @@ -29,6 +29,7 @@ class Fix : protected Pointers { int next_reneighbor; // next timestep to force a reneighboring int thermo_energy; // 1 if ThEng enabled via fix_modify, 0 if not int nevery; // how often to call an end_of_step fix + int pressure_every; // how often fix needs virial computed int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) @@ -42,7 +43,7 @@ class Fix : protected Pointers { int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR; // mask settings int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY; int INITIAL_INTEGRATE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA; - int MIN_POST_FORCE; + int MIN_POST_FORCE,MIN_ENERGY; Fix(class LAMMPS *, int, char **); virtual ~Fix(); @@ -76,6 +77,9 @@ class Fix : protected Pointers { virtual void final_integrate_respa(int) {} virtual void min_post_force(int) {} + virtual double min_energy(double *, double *) {return 0.0;} + virtual int min_dof() {return 0;} + virtual void min_xinitial(double *) {} virtual int pack_comm(int, int *, double *, int *) {return 0;} virtual void unpack_comm(int, int, double *) {} diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index ce149ec859..a760575075 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -44,6 +44,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : if (narg < 4) error->all("Illegal fix nph command"); restart_global = 1; + pressure_every = 1; double p_period[3]; if (strcmp(arg[3],"xyz") == 0) { diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 68e3e213ae..cb73f8e296 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -44,6 +44,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : if (narg < 7) error->all("Illegal fix npt command"); restart_global = 1; + pressure_every = 1; t_start = atof(arg[3]); t_stop = atof(arg[4]); diff --git a/src/modify.cpp b/src/modify.cpp index cc8415db5b..5f7e568cf2 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -47,6 +47,7 @@ using namespace LAMMPS_NS; #define POST_FORCE_RESPA 256 #define FINAL_INTEGRATE_RESPA 512 #define MIN_POST_FORCE 1024 +#define MIN_ENERGY 2048 /* ---------------------------------------------------------------------- */ @@ -57,7 +58,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) n_pre_exchange = n_pre_neighbor = 0; n_post_force = n_final_integrate = n_end_of_step = n_thermo_energy = 0; n_initial_integrate_respa = n_post_force_respa = n_final_integrate_respa = 0; - n_min_post_force = 0; + n_min_post_force = n_min_energy = 0; fix = NULL; fmask = NULL; @@ -67,7 +68,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) list_thermo_energy = NULL; list_initial_integrate_respa = list_post_force_respa = NULL; list_final_integrate_respa = NULL; - list_min_post_force = NULL; + list_min_post_force = list_min_energy = NULL; end_of_step_every = NULL; @@ -103,6 +104,7 @@ Modify::~Modify() delete [] list_post_force_respa; delete [] list_final_integrate_respa; delete [] list_min_post_force; + delete [] list_min_energy; delete [] end_of_step_every; @@ -149,6 +151,7 @@ void Modify::init() n_final_integrate_respa,list_final_integrate_respa); list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force); + list_init(MIN_ENERGY,n_min_energy,list_min_energy); // init each compute @@ -283,6 +286,50 @@ void Modify::min_post_force(int vflag) fix[list_min_post_force[i]]->min_post_force(vflag); } +/* ---------------------------------------------------------------------- + minimizer energy, force evaluation only for relevant fixes + return energy and forces on extra degrees of freedom +------------------------------------------------------------------------- */ + +double Modify::min_energy(double *xextra, double *fextra) +{ + int ifix,index; + + index = 0; + double energy_extra = 0.0; + for (int i = 0; i < n_min_energy; i++) { + ifix = list_min_energy[i]; + energy_extra += fix[ifix]->min_energy(&xextra[index],&fextra[index]); + index += fix[ifix]->min_dof(); + } + return energy_extra; +} + +/* ---------------------------------------------------------------------- + minimizer extra degrees of freedom from relevant fixes +------------------------------------------------------------------------- */ + +int Modify::min_dof() +{ + int ndof = 0; + for (int i = 0; i < n_min_energy; i++) + ndof += fix[list_min_energy[i]]->min_dof(); + return ndof; +} + +/* ---------------------------------------------------------------------- + minimizer initial xextra values only from relevant fixes +------------------------------------------------------------------------- */ + +void Modify::min_xinitial(double *xextra) +{ + int index = 0; + for (int i = 0; i < n_min_energy; i++) { + fix[list_min_energy[i]]->min_xinitial(&xextra[index]); + index += fix[list_min_energy[i]]->min_dof(); + } +} + /* ---------------------------------------------------------------------- add a new fix or replace one with same ID ------------------------------------------------------------------------- */ diff --git a/src/modify.h b/src/modify.h index 2a6ae0f064..e49bef7f74 100644 --- a/src/modify.h +++ b/src/modify.h @@ -25,7 +25,7 @@ class Modify : protected Pointers { int n_initial_integrate,n_pre_decide,n_pre_exchange,n_pre_neighbor; int n_post_force,n_final_integrate,n_end_of_step,n_thermo_energy; int n_initial_integrate_respa,n_post_force_respa,n_final_integrate_respa; - int n_min_post_force; + int n_min_post_force,n_min_energy; int nfix_restart_peratom; class Fix **fix; // list of fixes @@ -52,6 +52,9 @@ class Modify : protected Pointers { void final_integrate_respa(int); void min_post_force(int); + double min_energy(double *, double *); + int min_dof(); + void min_xinitial(double *); void add_fix(int, char **); void modify_fix(int, char **); @@ -77,7 +80,7 @@ class Modify : protected Pointers { int *list_thermo_energy; int *list_initial_integrate_respa,*list_post_force_respa; int *list_final_integrate_respa; - int *list_min_post_force; + int *list_min_post_force,*list_min_energy; int *end_of_step_every; diff --git a/src/respa.cpp b/src/respa.cpp index 995a0f8a13..becb5949ce 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -40,6 +40,9 @@ using namespace LAMMPS_NS; +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + /* ---------------------------------------------------------------------- */ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) @@ -273,25 +276,30 @@ void Respa::init() if (force->pair && force->pair->respa_enable == 0) error->all("Pair style does not support rRESPA inner/middle/outer"); - // set flags for how virial should be computed when needed - // pressure_flag is 1 if NPT,NPH - // virial_every is how virial should be computed every timestep - // 0 = not computed, 1 = computed explicity by pair - // virial_thermo is how virial should be computed on thermo timesteps - // 1 = computed explicity by pair - // unlike Verlet, virial is never computed implicitly + // setup virial computations for timestepping + // virial_style = 1 (explicit) since never computed implicitly like Verlet + // virial_every = 1 if computed every timestep (NPT,NPH) + // fix arrays store info on fixes that need virial computed occasionally - int pressure_flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1; - if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1; + virial_style = 1; + + virial_every = 0; + nfix_virial = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->pressure_every == 1) virial_every = 1; + else if (modify->fix[i]->pressure_every > 1) nfix_virial++; + + if (nfix_virial) { + delete [] fix_virial_every; + delete [] next_fix_virial; + fix_virial_every = new int[nfix_virial]; + next_fix_virial = new int[nfix_virial]; + nfix_virial = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->pressure_every > 1) + fix_virial_every[nfix_virial++] = modify->fix[i]->pressure_every; } - if (pressure_flag) virial_every = 1; - else virial_every = 0; - - virial_thermo = 1; - // step[] = timestep for each level step[nlevels-1] = update->dt; @@ -338,8 +346,8 @@ void Respa::setup() // compute all forces - int eflag = 1; - int vflag = virial_thermo; + eflag = 1; + vflag = virial_style; for (int ilevel = 0; ilevel < nlevels; ilevel++) { force_clear(newton[ilevel]); @@ -370,6 +378,21 @@ void Respa::setup() modify->setup(); sum_flevel_f(); output->setup(1); + + // setup virial computations for timestepping + + int ntimestep = update->ntimestep; + next_virial = 0; + if (virial_every) next_virial = ntimestep + 1; + else { + for (int ivirial = 0; ivirial < nfix_virial; ivirial++) { + next_fix_virial[ivirial] = + (ntimestep/fix_virial_every[ivirial])*fix_virial_every[ivirial] + + fix_virial_every[ivirial]; + if (ivirial) next_virial = MIN(next_virial,next_fix_virial[ivirial]); + else next_virial = next_fix_virial[0]; + } + } } /* ---------------------------------------------------------------------- @@ -378,22 +401,28 @@ void Respa::setup() void Respa::iterate(int n) { + int ntimestep; + for (int i = 0; i < n; i++) { - update->ntimestep++; + ntimestep = ++update->ntimestep; - eflag = 0; - vflag = virial_every; - if (output->next_thermo == update->ntimestep) { - eflag = 1; - vflag = virial_thermo; - } + // eflag/vflag = 0/1 for energy/virial computation + + if (ntimestep == output->next_thermo) eflag = 1; + else eflag = 0; + + if (ntimestep == output->next_thermo || ntimestep == next_virial) { + vflag = virial_style; + if (virial_every) next_virial++; + else next_virial = fix_virial(ntimestep); + } else vflag = 0; recurse(nlevels-1); if (modify->n_end_of_step) modify->end_of_step(); - if (output->next == update->ntimestep) { + if (ntimestep == output->next) { timer->stamp(); sum_flevel_f(); output->write(update->ntimestep); @@ -582,3 +611,20 @@ void Respa::sum_flevel_f() } } } + +/* ---------------------------------------------------------------------- + return next timestep virial should be computed + based on one or more fixes that need virial computed periodically +------------------------------------------------------------------------- */ + +int Respa::fix_virial(int ntimestep) +{ + int next; + for (int ivirial = 0; ivirial < nfix_virial; ivirial++) { + if (ntimestep == next_fix_virial[ivirial]) + next_fix_virial[ivirial] += fix_virial_every[ivirial]; + if (ivirial) next = MIN(next,next_fix_virial[ivirial]); + else next = next_fix_virial[0]; + } + return next; +} diff --git a/src/respa.h b/src/respa.h index d11a16ab4b..e592f77fd8 100644 --- a/src/respa.h +++ b/src/respa.h @@ -44,16 +44,21 @@ class Respa : public Integrate { void copy_flevel_f(int); private: - int virial_every; // what vflag should be on every timestep (0,1) - int virial_thermo; // what vflag should be on thermo steps (1) + int eflag,vflag; // flags for energy/virial computation + int virial_style; // compute virial explicitly (not implicit) + int virial_every; // 1 if virial computed every step + int next_virial; // next timestep to compute virial + int nfix_virial; // # of fixes that need virial occasionally + int *fix_virial_every; // frequency they require it + int *next_fix_virial; // next timestep they need it - int *newton; // newton flag at each level - int eflag,vflag; // flags for energy/virial computation - class FixRespa *fix_respa; // Fix to store the force level array + int *newton; // newton flag at each level + class FixRespa *fix_respa; // Fix to store the force level array void recurse(int); void force_clear(int); void sum_flevel_f(); + int fix_virial(int); }; } diff --git a/src/style.h b/src/style.h index d97fcf7f02..8c1e094c08 100644 --- a/src/style.h +++ b/src/style.h @@ -84,6 +84,7 @@ CommandStyle(write_restart,WriteRestart) #include "compute_temp_partial.h" #include "compute_temp_ramp.h" #include "compute_temp_region.h" +#include "compute_variable_atom.h" #endif #ifdef ComputeClass @@ -99,6 +100,7 @@ ComputeStyle(temp,ComputeTemp) ComputeStyle(temp/partial,ComputeTempPartial) ComputeStyle(temp/ramp,ComputeTempRamp) ComputeStyle(temp/region,ComputeTempRegion) +ComputeStyle(variable/atom,ComputeVariableAtom) #endif #ifdef DihedralInclude @@ -124,6 +126,9 @@ DumpStyle(xyz,DumpXYZ) #ifdef FixInclude #include "fix_add_force.h" #include "fix_ave_force.h" +#include "fix_ave_spatial.h" +#include "fix_ave_time.h" + //#include "fix_box_relax.h" #include "fix_com.h" #include "fix_drag.h" #include "fix_deposit.h" @@ -168,6 +173,9 @@ DumpStyle(xyz,DumpXYZ) #ifdef FixClass FixStyle(addforce,FixAddForce) FixStyle(aveforce,FixAveForce) +FixStyle(ave/spatial,FixAveSpatial) +FixStyle(ave/time,FixAveTime) + //FixStyle(box/relax,FixBoxRelax) FixStyle(com,FixCOM) FixStyle(drag,FixDrag) FixStyle(deposit,FixDeposit) diff --git a/src/variable.cpp b/src/variable.cpp index 96f3ee5ecf..dd687eb6d4 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -33,7 +33,8 @@ using namespace LAMMPS_NS; #define VARDELTA 4 -enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP}; +enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP,ATOM}; +enum{VALUE,ATOMARRAY,TYPEARRAY,ADD,SUB,MULT,DIV,NEG,POW,EXP,LN,SQRT}; /* ---------------------------------------------------------------------- */ @@ -73,9 +74,9 @@ void Variable::set(int narg, char **arg) { if (narg < 3) error->all("Illegal variable command"); - // if var already exists, just skip (except EQUAL vars) + // if var already exists, just skip - if (find(arg[0]) >= 0 && strcmp(arg[1],"equal") != 0) return; + if (find(arg[0]) >= 0) return; // make space for new variable @@ -90,6 +91,11 @@ void Variable::set(int narg, char **arg) memory->srealloc(data,maxvar*sizeof(char **),"var:data"); } + // set name of variable + + names[nvar] = new char[strlen(arg[0]) + 1]; + strcpy(names[nvar],arg[0]); + // INDEX // num = listed args, index = 1st value, data = copied args @@ -111,16 +117,10 @@ void Variable::set(int narg, char **arg) for (int i = 0; i < num[nvar]; i++) data[nvar][i] = NULL; // EQUAL - // remove pre-existing var if also style EQUAL (allows it to be reset) // num = 2, index = 1st value // data = 2 values, 1st is string to eval, 2nd is filled on retrieval } else if (strcmp(arg[1],"equal") == 0) { - if (find(arg[0]) >= 0) { - if (style[find(arg[0])] != EQUAL) - error->all("Cannot redefine variable as a different style"); - remove(find(arg[0])); - } style[nvar] = EQUAL; num[nvar] = 2; index[nvar] = 0; @@ -187,17 +187,24 @@ void Variable::set(int narg, char **arg) arg[0],index[nvar]+1,universe->iworld); } + // ATOM + // num = 1, index = 1st value + // data = 1 value, string to eval + + } else if (strcmp(arg[1],"atom") == 0) { + style[nvar] = ATOM; + num[nvar] = 1; + index[nvar] = 0; + data[nvar] = new char*[num[nvar]]; + copy(1,&arg[2],data[nvar]); + } else error->all("Illegal variable command"); - // set variable name (after possible EQUAL remove) - - names[nvar] = new char[strlen(arg[0]) + 1]; - strcpy(names[nvar],arg[0]); nvar++; } /* ---------------------------------------------------------------------- - single-value INDEX variable created by command-line argument + single-value EQUAL variable created by command-line argument ------------------------------------------------------------------------- */ void Variable::set(char *name, char *value) @@ -205,25 +212,12 @@ void Variable::set(char *name, char *value) int ivar = find(name); if (ivar >= 0) error->all("Command-line variable already exists"); - if (nvar == maxvar) { - maxvar += VARDELTA; - names = (char **) - memory->srealloc(names,maxvar*sizeof(char *),"var:names"); - style = (int *) memory->srealloc(style,maxvar*sizeof(int),"var:style"); - num = (int *) memory->srealloc(num,maxvar*sizeof(int),"var:num"); - index = (int *) memory->srealloc(index,maxvar*sizeof(int),"var:index"); - data = (char ***) - memory->srealloc(data,maxvar*sizeof(char **),"var:data"); - } - - names[nvar] = new char[strlen(name) + 1]; - strcpy(names[nvar],name); - style[nvar] = INDEX; - num[nvar] = 1; - index[nvar] = 0; - data[nvar] = new char*[num[nvar]]; - copy(1,&value,data[nvar]); - nvar++; + char **newarg = new char*[3]; + newarg[0] = name; + newarg[1] = "equal"; + newarg[2] = value; + set(3,newarg); + delete [] newarg; } /* ---------------------------------------------------------------------- @@ -247,10 +241,10 @@ int Variable::next(int narg, char **arg) error->all("All variables in next command must be same style"); } - // check for invalid styles EQUAL or WORLD + // invalid styles EQUAL or WORLD or ATOM int istyle = style[find(arg[0])]; - if (istyle == EQUAL || istyle == WORLD) + if (istyle == EQUAL || istyle == WORLD || istyle == ATOM) error->all("Invalid variable style with next command"); // increment all variables in list @@ -336,14 +330,16 @@ char *Variable::retrieve(char *name) delete [] value; str = data[ivar][0]; } else if (style[ivar] == EQUAL) { - char *value = evaluate(data[ivar][0]); - int n = strlen(value) + 1; + char result[32]; + double answer = evaluate(data[ivar][0],NULL); + sprintf(result,"%.10g",answer); + int n = strlen(result) + 1; if (data[ivar][1]) delete [] data[ivar][1]; data[ivar][1] = new char[n]; - strcpy(data[ivar][1],value); - delete [] value; + strcpy(data[ivar][1],result); str = data[ivar][1]; - } + } else if (style[ivar] == ATOM) return NULL; + return str; } @@ -404,34 +400,34 @@ void Variable::copy(int narg, char **from, char **to) group function = mass(group), xcm(group,x), ... atom vector = x[123], y[3], vx[34], ... compute vector = c_mytemp[0], c_thermo_press[3], ... + numbers start with a digit or "." or "-" (no parens or brackets) + keywords start with a lowercase letter (no parens or brackets) functions contain () can have 1 or 2 args, each of which can be a "string" vectors contain [] single arg must be integer - for atom vectors, it is global ID of atom + for atom vectors, arg is global ID of atom for compute vectors, 0 is the scalar value, 1-N are vector values - keywords start with a lowercase letter (no parens or brackets) - numbers start with a digit or "." or "-" (no parens or brackets) - see lists of valid functions, vectors, keywords below - when string is evaluated, put result in a newly allocated string - return the address of result string (will be freed by caller) + see lists of valid functions & vectors below + return answer = value of string ------------------------------------------------------------------------- */ -char *Variable::evaluate(char *str) +double Variable::evaluate(char *str, Tree *tree) { - // allocate a new string for the eventual result - - char *result = new char[32]; - double answer; + double answer = 0.0; + if (tree) { + tree->type = VALUE; + tree->left = tree->right = NULL; + } // string is a "function" since contains () // grab one or two args, separated by "," // evaulate args recursively - // then evaluate math or group function + // if tree = NULL, evaluate math or group function + // else store as leaf in tree if (strchr(str,'(')) { - if (str[strlen(str)-1] != ')') - error->all("Cannot evaluate variable equal command"); + if (str[strlen(str)-1] != ')') error->all("Cannot evaluate variable"); char *ptr = strchr(str,'('); int n = ptr - str; @@ -442,8 +438,7 @@ char *Variable::evaluate(char *str) char *comma = ++ptr; int level = 0; while (1) { - if (*comma == '\0') - error->all("Cannot evaluate variable equal command"); + if (*comma == '\0') error->all("Cannot evaluate variable"); else if (*comma == ',' && level == 0) break; else if (*comma == ')' && level == 0) break; else if (*comma == '(') level++; @@ -467,106 +462,132 @@ char *Variable::evaluate(char *str) } else arg2 = NULL; double value1,value2; - char *strarg1 = NULL; - char *strarg2 = NULL; // customize by adding math function to this list and to if statement // add(x,y),sub(x,y),mult(x,y),div(x,y),neg(x), // pow(x,y),exp(x),ln(x),sqrt(x) if (strcmp(func,"add") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - strarg2 = evaluate(arg2); - value1 = atof(strarg1); - value2 = atof(strarg2); - answer = value1 + value2; + if (!arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = ADD; + tree->left = new Tree(); + tree->right = new Tree(); + value1 = evaluate(arg1,tree->left); + value2 = evaluate(arg2,tree->right); + } else answer = evaluate(arg1,NULL) + evaluate(arg2,NULL); } else if (strcmp(func,"sub") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - strarg2 = evaluate(arg2); - value1 = atof(strarg1); - value2 = atof(strarg2); - answer = value1 - value2; + if (!arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = SUB; + tree->left = new Tree(); + tree->right = new Tree(); + value1 = evaluate(arg1,tree->left); + value2 = evaluate(arg2,tree->right); + } else answer = evaluate(arg1,NULL) - evaluate(arg2,NULL); } else if (strcmp(func,"mult") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - strarg2 = evaluate(arg2); - value1 = atof(strarg1); - value2 = atof(strarg2); - answer = value1 * value2; + if (!arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = MULT; + tree->left = new Tree(); + tree->right = new Tree(); + value1 = evaluate(arg1,tree->left); + value2 = evaluate(arg2,tree->right); + } else answer = evaluate(arg1,NULL) * evaluate(arg2,NULL); } else if (strcmp(func,"div") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - strarg2 = evaluate(arg2); - value1 = atof(strarg1); - value2 = atof(strarg2); - if (value2 == 0.0) - error->all("Cannot evaluate variable equal command"); - answer = value1 / value2; + if (!arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = DIV; + tree->left = new Tree(); + tree->right = new Tree(); + value1 = evaluate(arg1,tree->left); + value2 = evaluate(arg2,tree->right); + } else { + value1 = evaluate(arg1,NULL); + value2 = evaluate(arg2,NULL); + if (value2 == 0.0) error->all("Cannot evaluate variable"); + answer = value1 / value2; + } } else if (strcmp(func,"neg") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - value1 = atof(strarg1); - answer = -value1; + if (arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = NEG; + tree->left = new Tree(); + value1 = evaluate(arg1,tree->left); + } else answer = -evaluate(arg1,NULL); } else if (strcmp(func,"pow") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - strarg2 = evaluate(arg2); - value1 = atof(strarg1); - value2 = atof(strarg2); - if (value2 == 0.0) - error->all("Cannot evaluate variable equal command"); - answer = pow(value1,value2); + if (!arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = POW; + tree->left = new Tree(); + tree->right = new Tree(); + value1 = evaluate(arg1,tree->left); + value2 = evaluate(arg2,tree->right); + } else { + value1 = evaluate(arg1,NULL); + value2 = evaluate(arg2,NULL); + if (value2 == 0.0) error->all("Cannot evaluate variable"); + answer = pow(value1,value2); + } } else if (strcmp(func,"exp") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - value1 = atof(strarg1); - answer = exp(value1); + if (arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = EXP; + tree->left = new Tree(); + value1 = evaluate(arg1,tree->left); + } else answer = exp(evaluate(arg1,NULL)); } else if (strcmp(func,"ln") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - value1 = atof(strarg1); - if (value1 == 0.0) - error->all("Cannot evaluate variable equal command"); - answer = log(value1); + if (arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = LN; + tree->left = new Tree(); + value1 = evaluate(arg1,tree->left); + } else { + value1 = evaluate(arg1,NULL); + if (value1 == 0.0) error->all("Cannot evaluate variable"); + answer = log(value1); + } } else if (strcmp(func,"sqrt") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); - strarg1 = evaluate(arg1); - value1 = atof(strarg1); - if (value1 == 0.0) - error->all("Cannot evaluate variable equal command"); - answer = sqrt(value1); + if (arg2) error->all("Cannot evaluate variable"); + if (tree) { + tree->type = SQRT; + tree->left = new Tree(); + value1 = evaluate(arg1,tree->left); + } else { + value1 = evaluate(arg1,NULL); + if (value1 == 0.0) error->all("Cannot evaluate variable"); + answer = sqrt(value1); + } // customize by adding group function to this list and to if statement // mass(group),charge(group),xcm(group,dim),vcm(group,dim), // bound(group,xmin),gyration(group) } else if (strcmp(func,"mass") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); + if (arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); atom->check_mass(); answer = group->mass(igroup); } else if (strcmp(func,"charge") == 0) { - if (arg2) error->all("Cannot evaluate variable equal command"); + if (arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); answer = group->charge(igroup); } else if (strcmp(func,"xcm") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); + if (!arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); atom->check_mass(); double masstotal = group->mass(igroup); double xcm[3]; @@ -574,12 +595,12 @@ char *Variable::evaluate(char *str) if (strcmp(arg2,"x") == 0) answer = xcm[0]; else if (strcmp(arg2,"y") == 0) answer = xcm[1]; else if (strcmp(arg2,"z") == 0) answer = xcm[2]; - else error->all("Cannot evaluate variable equal command"); + else error->all("Cannot evaluate variable"); } else if (strcmp(func,"vcm") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); + if (!arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); atom->check_mass(); double masstotal = group->mass(igroup); double vcm[3]; @@ -587,12 +608,12 @@ char *Variable::evaluate(char *str) if (strcmp(arg2,"x") == 0) answer = vcm[0]; else if (strcmp(arg2,"y") == 0) answer = vcm[1]; else if (strcmp(arg2,"z") == 0) answer = vcm[2]; - else error->all("Cannot evaluate variable equal command"); + else error->all("Cannot evaluate variable"); } else if (strcmp(func,"bound") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); + if (!arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); double minmax[6]; group->bounds(igroup,minmax); if (strcmp(arg2,"xmin") == 0) answer = minmax[0]; @@ -601,28 +622,23 @@ char *Variable::evaluate(char *str) else if (strcmp(arg2,"ymax") == 0) answer = minmax[3]; else if (strcmp(arg2,"zmin") == 0) answer = minmax[4]; else if (strcmp(arg2,"zmax") == 0) answer = minmax[5]; - else error->all("Cannot evaluate variable equal command"); + else error->all("Cannot evaluate variable"); } else if (strcmp(func,"gyration") == 0) { - if (!arg2) error->all("Cannot evaluate variable equal command"); + if (!arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1); - if (igroup == -1) error->all("Variable equal group ID does not exist"); + if (igroup == -1) error->all("Variable group ID does not exist"); atom->check_mass(); double masstotal = group->mass(igroup); double xcm[3]; group->xcm(igroup,masstotal,xcm); answer = group->gyration(igroup,masstotal,xcm); - } else error->all("Invalid math/group function in variable equal command"); + } else error->all("Invalid math/group function in variable"); delete [] func; delete [] arg1; - delete [] strarg1; - if (arg2) { - delete [] arg2; - delete [] strarg2; - } - sprintf(result,"%.10g",answer); + delete [] arg2; // string is a "vector" since contains [] // index = everything between [] evaluated as integer @@ -633,8 +649,7 @@ char *Variable::evaluate(char *str) // grab atom-based value with index as global atom ID } else if (strchr(str,'[')) { - if (str[strlen(str)-1] != ']') - error->all("Cannot evaluate variable equal command"); + if (str[strlen(str)-1] != ']') error->all("Cannot evaluate variable"); char *ptr = strchr(str,'['); int n = ptr - str; @@ -659,7 +674,7 @@ char *Variable::evaluate(char *str) for (icompute = 0; icompute < modify->ncompute; icompute++) if (strcmp(id,modify->compute[icompute]->id) == 0) break; if (icompute == modify->ncompute) - error->all("Invalid compute ID in variable equal"); + error->all("Invalid compute ID in variable"); delete [] id; modify->compute[icompute]->init(); @@ -672,7 +687,7 @@ char *Variable::evaluate(char *str) error->all("Variable compute ID does not compute scalar info"); if (modify->compute[icompute]->id_pre) { int ipre = modify->find_compute(modify->compute[icompute]->id_pre); - if (ipre < 0) error->all("Could not pre-compute in variable equal"); + if (ipre < 0) error->all("Could not pre-compute in variable"); answer = modify->compute[ipre]->compute_scalar(); } answer = modify->compute[icompute]->compute_scalar(); @@ -683,28 +698,62 @@ char *Variable::evaluate(char *str) error->all("Variable compute ID vector is not large enough"); if (modify->compute[icompute]->id_pre) { int ipre = modify->find_compute(modify->compute[icompute]->id_pre); - if (ipre < 0) error->all("Could not pre-compute in variable equal"); + if (ipre < 0) error->all("Could not pre-compute in variable"); modify->compute[ipre]->compute_vector(); } modify->compute[icompute]->compute_vector(); answer = modify->compute[icompute]->vector[index-1]; - } else error->all("Invalid compute ID index in variable equal"); + } else error->all("Invalid compute ID index in variable"); + + } else if (tree) { + + if (strlen(arg)) error->all("Invalid atom vector argument in variable"); + + // customize by adding atom vector to this list and to if statement + // mass,x,y,z,vx,vy,vz,fx,fy,fz + + tree->type = ATOMARRAY; + tree->nstride = 3; + + if (strcmp(vector,"mass") == 0) { + if (atom->mass) { + tree->type = TYPEARRAY; + tree->array = atom->mass; + } else { + tree->nstride = 1; + tree->array = atom->rmass; + } + } + else if (strcmp(vector,"x") == 0) tree->array = &atom->x[0][0]; + else if (strcmp(vector,"y") == 0) tree->array = &atom->x[0][1]; + else if (strcmp(vector,"z") == 0) tree->array = &atom->x[0][2]; + else if (strcmp(vector,"vx") == 0) tree->array = &atom->v[0][0]; + else if (strcmp(vector,"vy") == 0) tree->array = &atom->v[0][1]; + else if (strcmp(vector,"vz") == 0) tree->array = &atom->v[0][2]; + else if (strcmp(vector,"fx") == 0) tree->array = &atom->f[0][0]; + else if (strcmp(vector,"fy") == 0) tree->array = &atom->f[0][1]; + else if (strcmp(vector,"fz") == 0) tree->array = &atom->f[0][2]; + + else error->all("Invalid atom vector in variable"); } else { if (atom->map_style == 0) - error->all("Cannot use atom vectors in variable " - "unless atom map exists"); + error->all("Cannot use atom vector in variable unless atom map exists"); int index = atom->map(atoi(arg)); // customize by adding atom vector to this list and to if statement - // x,y,z,vx,vy,vz,fx,fy,fz + // mass,x,y,z,vx,vy,vz,fx,fy,fz double mine; if (index >= 0 && index < atom->nlocal) { - if (strcmp(vector,"x") == 0) mine = atom->x[index][0]; + if (strcmp(vector,"mass") == 0) { + if (atom->mass) mine = atom->mass[atom->type[index]]; + else mine = atom->rmass[index]; + } + else if (strcmp(vector,"x") == 0) mine = atom->x[index][0]; else if (strcmp(vector,"y") == 0) mine = atom->x[index][1]; else if (strcmp(vector,"z") == 0) mine = atom->x[index][2]; else if (strcmp(vector,"vx") == 0) mine = atom->v[index][0]; @@ -714,7 +763,7 @@ char *Variable::evaluate(char *str) else if (strcmp(vector,"fy") == 0) mine = atom->f[index][1]; else if (strcmp(vector,"fz") == 0) mine = atom->f[index][2]; - else error->one("Invalid atom vector in variable equal command"); + else error->one("Invalid atom vector in variable"); } else mine = 0.0; @@ -723,7 +772,6 @@ char *Variable::evaluate(char *str) delete [] vector; delete [] arg; - sprintf(result,"%.10g",answer); // string is "keyword" since starts with lowercase letter // if keyword starts with "v_", trailing chars are variable name @@ -731,8 +779,6 @@ char *Variable::evaluate(char *str) // else is thermo keyword // evaluate it via evaluate_keyword() - // compute appropriate value via thermo - } else if (islower(str[0])) { if (strncmp(str,"v_",2) == 0) { @@ -741,31 +787,98 @@ char *Variable::evaluate(char *str) strcpy(id,&str[2]); char *v = retrieve(id); - if (v == NULL) - error->all("Invalid variable name in variable equal command"); + if (v == NULL) error->all("Invalid variable name in variable"); delete [] id; answer = atof(v); } else { int flag = output->thermo->evaluate_keyword(str,&answer); - if (flag) error->all("Invalid thermo keyword in variable equal command"); + if (flag) error->all("Invalid thermo keyword in variable"); } - sprintf(result,"%.10g",answer); - // string is a number since starts with digit or "." or "-" // just copy to result } else if (isdigit(str[0]) || str[0] == '.' || str[0] == '-') { - strcpy(result,str); + answer = atof(str); // string is an error - } else error->all("Cannot evaluate variable equal command"); + } else error->all("Cannot evaluate variable"); - // return newly allocated string + // store answer in tree and return it - return result; + if (tree) tree->value = answer; + return answer; +} + +/* ---------------------------------------------------------------------- */ + +void Variable::build_parse_tree(int ivar) +{ + if (style[ivar] != ATOM) error->all(""); + ptree = new Tree(); + double tmp = evaluate(data[ivar][0],ptree); +} + +/* ---------------------------------------------------------------------- */ + +void Variable::evaluate_parse_tree(int igroup, double *result) +{ + int groupbit = group->bitmask[igroup]; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] && groupbit) result[i] = eval_tree(ptree,i); + else result[i] = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +void Variable::free_parse_tree() +{ + free_tree(ptree); +} + +/* ---------------------------------------------------------------------- */ + +double Variable::eval_tree(Tree *tree, int i) +{ + if (tree->type == VALUE) return tree->value; + if (tree->type == ATOMARRAY) return tree->array[i*tree->nstride]; + if (tree->type == TYPEARRAY) return tree->array[atom->type[i]]; + + if (tree->type == ADD) + return eval_tree(tree->left,i) + eval_tree(tree->right,i); + if (tree->type == SUB) + return eval_tree(tree->left,i) - eval_tree(tree->right,i); + if (tree->type == MULT) + return eval_tree(tree->left,i) * eval_tree(tree->right,i); + if (tree->type == DIV) + return eval_tree(tree->left,i) / eval_tree(tree->right,i); + if (tree->type == NEG) + return -eval_tree(tree->left,i); + if (tree->type == POW) + return pow(eval_tree(tree->left,i),eval_tree(tree->right,i)); + if (tree->type == EXP) + return exp(eval_tree(tree->left,i)); + if (tree->type == LN) + return log(eval_tree(tree->left,i)); + if (tree->type == SQRT) + return sqrt(eval_tree(tree->left,i)); + + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +void Variable::free_tree(Tree *tree) +{ + if (tree->left) free_tree(tree->left); + if (tree->right) free_tree(tree->right); + delete tree; } diff --git a/src/variable.h b/src/variable.h index 860f93200d..c1b5f2d919 100644 --- a/src/variable.h +++ b/src/variable.h @@ -28,6 +28,10 @@ class Variable : protected Pointers { int find(char *); char *retrieve(char *); + void build_parse_tree(int); + void evaluate_parse_tree(int, double *); + void free_parse_tree(); + private: int me; int nvar; // # of defined variables @@ -38,9 +42,21 @@ class Variable : protected Pointers { int *index; // next available value for each variable char ***data; // str value of each variable's values + struct Tree { + double value; + double *array; + int nstride; + int type; + Tree *left,*right; + }; + + Tree *ptree; // parse tree for an ATOM variable + void copy(int, char **, char **); - char *evaluate(char *); + double evaluate(char *, Tree *); void remove(int); + double eval_tree(Tree *, int); + void free_tree(Tree *); }; } diff --git a/src/verlet.h b/src/verlet.h index ef9fbb0f84..1aae82c7f6 100644 --- a/src/verlet.h +++ b/src/verlet.h @@ -21,19 +21,25 @@ namespace LAMMPS_NS { class Verlet : public Integrate { public: Verlet(class LAMMPS *, int, char **); + ~Verlet(); void init(); void setup(); void iterate(int); private: - int virial_every; // what vflag should be on every timestep (0,1,2) - int virial_thermo; // what vflag should be on thermo steps (1,2) - int pairflag,torqueflag,granflag; + int virial_style; // compute virial explicitly or implicitly + int virial_every; // 1 if virial computed every step + int next_virial; // next timestep to compute virial + int nfix_virial; // # of fixes that need virial occasionally + int *fix_virial_every; // frequency they require it + int *next_fix_virial; // next timestep they need it - int maxpair; // copies of Update quantities + int pairflag,torqueflag,granflag; // arrays to zero out every step + int maxpair; // local copies of Update quantities double **f_pair; void force_clear(int); + int fix_virial(int); }; }