diff --git a/src/fix.cpp b/src/fix.cpp index e6e8a292b..80fa00f4b 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -1,318 +1,393 @@ /* ---------------------------------------------------------------------- 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 #include #include "fix.h" #include "atom.h" #include "group.h" #include "force.h" #include "comm.h" #include "atom_masks.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; // allocate space for static class instance variable and initialize it int Fix::instance_total = 0; /* ---------------------------------------------------------------------- */ -Fix::Fix(LAMMPS *lmp, int narg, char **arg) : +Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp), id(NULL), style(NULL), extlist(NULL), vector_atom(NULL), array_atom(NULL), vector_local(NULL), array_local(NULL), eatom(NULL), vatom(NULL) { instance_me = instance_total++; // fix ID, group, and style // ID must be all alphanumeric chars or underscores int n = strlen(arg[0]) + 1; id = new char[n]; strcpy(id,arg[0]); for (int i = 0; i < n-1; i++) if (!isalnum(id[i]) && id[i] != '_') error->all(FLERR,"Fix ID must be alphanumeric or underscore characters"); igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find fix group ID"); groupbit = group->bitmask[igroup]; n = strlen(arg[2]) + 1; style = new char[n]; strcpy(style,arg[2]); restart_global = restart_peratom = restart_file = 0; force_reneighbor = 0; box_change_size = box_change_shape = box_change_domain = 0; thermo_energy = 0; + thermo_virial = 0; rigid_flag = 0; peatom_flag = 0; virial_flag = 0; no_change_box = 0; time_integrate = 0; time_depend = 0; create_attribute = 0; restart_pbc = 0; wd_header = wd_section = 0; dynamic_group_allow = 0; dynamic = 0; dof_flag = 0; special_alter_flag = 0; enforce2d_flag = 0; respa_level_support = 0; respa_level = -1; scalar_flag = vector_flag = array_flag = 0; peratom_flag = local_flag = 0; size_vector_variable = size_array_rows_variable = 0; comm_forward = comm_reverse = comm_border = 0; restart_reset = 0; // reasonable defaults // however, each fix that uses these values should explicitly set them nevery = 1; global_freq = 1; // per-atom virial // set vflag_atom = 0 b/c some fixes grow vatom in grow_arrays() // which may occur outside of timestepping maxeatom = maxvatom = 0; vflag_atom = 0; // KOKKOS per-fix data masks execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; kokkosable = 0; copymode = 0; } /* ---------------------------------------------------------------------- */ Fix::~Fix() { if (copymode) return; delete [] id; delete [] style; memory->destroy(eatom); memory->destroy(vatom); } /* ---------------------------------------------------------------------- process params common to all fixes here if unknown param, call modify_param specific to the fix ------------------------------------------------------------------------- */ void Fix::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal fix_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"dynamic/dof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (strcmp(arg[iarg+1],"no") == 0) dynamic = 0; else if (strcmp(arg[iarg+1],"yes") == 0) dynamic = 1; else error->all(FLERR,"Illegal fix_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"energy") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0; - else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1; - else error->all(FLERR,"Illegal fix_modify command"); + else if (strcmp(arg[iarg+1],"yes") == 0) { + if (!(THERMO_ENERGY & setmask())) + error->all(FLERR,"Illegal fix_modify command"); + thermo_energy = 1; + } else error->all(FLERR,"Illegal fix_modify command"); + iarg += 2; + } else if (strcmp(arg[iarg],"virial") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); + if (strcmp(arg[iarg+1],"no") == 0) thermo_virial = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) { + if (virial_flag == 0) + error->all(FLERR,"Illegal fix_modify command"); + thermo_virial = 1; + } else error->all(FLERR,"Illegal fix_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"respa") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (!respa_level_support) error->all(FLERR,"Illegal fix_modify command"); int lvl = force->inumeric(FLERR,arg[iarg+1]); if (lvl < 0) error->all(FLERR,"Illegal fix_modify command"); respa_level = lvl-1; iarg += 2; } else { int n = modify_param(narg-iarg,&arg[iarg]); if (n == 0) error->all(FLERR,"Illegal fix_modify command"); iarg += n; } } } /* ---------------------------------------------------------------------- setup for energy, virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) fixes call this if they use ev_tally() ------------------------------------------------------------------------- */ void Fix::ev_setup(int eflag, int vflag) { int i,n; evflag = 1; eflag_either = eflag; eflag_global = eflag % 2; eflag_atom = eflag / 2; vflag_either = vflag; vflag_global = vflag % 4; vflag_atom = vflag / 4; // reallocate per-atom arrays if necessary if (eflag_atom && atom->nlocal > maxeatom) { maxeatom = atom->nmax; memory->destroy(eatom); memory->create(eatom,maxeatom,"fix:eatom"); } if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,maxvatom,6,"fix:vatom"); } // zero accumulators // no global energy variable to zero (unlike pair,bond,angle,etc) // fixes tally it individually via fix_modify energy yes and compute_scalar() if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (eflag_atom) { n = atom->nlocal; for (i = 0; i < n; i++) eatom[i] = 0.0; } if (vflag_atom) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; vatom[i][2] = 0.0; vatom[i][3] = 0.0; vatom[i][4] = 0.0; vatom[i][5] = 0.0; } } } /* ---------------------------------------------------------------------- - setup for virial computation - see integrate::ev_set() for values of vflag (0-6) - fixes call this if use v_tally() + if thermo_virial is on: + setup for virial computation + see integrate::ev_set() for values of vflag (0-6) + fixes call this if use v_tally() + else: set evflag=0 ------------------------------------------------------------------------- */ void Fix::v_setup(int vflag) { int i,n; + if (!thermo_virial) { + evflag = 0; + return; + } + evflag = 1; vflag_global = vflag % 4; vflag_atom = vflag / 4; // reallocate per-atom array if necessary if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,maxvatom,6,"fix:vatom"); } // zero accumulators if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (vflag_atom) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; vatom[i][2] = 0.0; vatom[i][3] = 0.0; vatom[i][4] = 0.0; vatom[i][5] = 0.0; } } } /* ---------------------------------------------------------------------- tally per-atom energy and global/per-atom virial into accumulators n = # of local owned atoms involved, with local indices in list eng = total energy for the interaction involving total atoms v = total virial for the interaction involving total atoms increment per-atom energy of each atom in list by 1/total fraction v_tally tallies virial this method can be used when fix computes energy/forces in post_force() e.g. fix cmap: compute energy and virial only on owned atoms whether newton_bond is on or off other procs will tally left-over fractions for atoms they own ------------------------------------------------------------------------- */ void Fix::ev_tally(int n, int *list, double total, double eng, double *v) { if (eflag_atom) { double fraction = eng/total; for (int i = 0; i < n; i++) eatom[list[i]] += fraction; } v_tally(n,list,total,v); } /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators n = # of local owned atoms involved, with local indices in list v = total virial for the interaction involving total atoms increment global virial by n/total fraction increment per-atom virial of each atom in list by 1/total fraction this method can be used when fix computes forces in post_force() e.g. fix shake, fix rigid: compute virial only on owned atoms whether newton_bond is on or off other procs will tally left-over fractions for atoms they own ------------------------------------------------------------------------- */ void Fix::v_tally(int n, int *list, double total, double *v) { int m; if (vflag_global) { double fraction = n/total; virial[0] += fraction*v[0]; virial[1] += fraction*v[1]; virial[2] += fraction*v[2]; virial[3] += fraction*v[3]; virial[4] += fraction*v[4]; virial[5] += fraction*v[5]; } if (vflag_atom) { double fraction = 1.0/total; for (int i = 0; i < n; i++) { m = list[i]; vatom[m][0] += fraction*v[0]; vatom[m][1] += fraction*v[1]; vatom[m][2] += fraction*v[2]; vatom[m][3] += fraction*v[3]; vatom[m][4] += fraction*v[4]; vatom[m][5] += fraction*v[5]; } } } + +/* ---------------------------------------------------------------------- + tally virial into global and per-atom accumulators + i = local index of atom + v = total virial for the interaction + increment global virial by v + increment per-atom virial by v + this method can be used when fix computes forces in post_force() + and the force depends on a distance to some external object + e.g. fix wall/lj93: compute virial only on owned atoms +------------------------------------------------------------------------- */ + +void Fix::v_tally(int i, double *v) +{ + + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + vatom[i][0] += v[0]; + vatom[i][1] += v[1]; + vatom[i][2] += v[2]; + vatom[i][3] += v[3]; + vatom[i][4] += v[4]; + vatom[i][5] += v[5]; + } +} + +/* ---------------------------------------------------------------------- + tally virial component into global and per-atom accumulators + n = index of virial component (0-5) + i = local index of atom + vn = nth component of virial for the interaction + increment nth component of global virial by vn + increment nth component of per-atom virial by vn + this method can be used when fix computes forces in post_force() + and the force depends on a distance to some external object + e.g. fix wall/lj93: compute virial only on owned atoms +------------------------------------------------------------------------- */ + +void Fix::v_tally(int n, int i, double vn) +{ + + if (vflag_global) + virial[n] += vn; + + if (vflag_atom) + vatom[i][n] += vn; +} diff --git a/src/fix.h b/src/fix.h index d91937848..3f3289530 100644 --- a/src/fix.h +++ b/src/fix.h @@ -1,285 +1,288 @@ /* -*- c++ -*- ---------------------------------------------------------- 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 LMP_FIX_H #define LMP_FIX_H #include "pointers.h" namespace LAMMPS_NS { class Fix : protected Pointers { public: static int instance_total; // # of Fix classes ever instantiated char *id,*style; int igroup,groupbit; int restart_global; // 1 if Fix saves global state, 0 if not int restart_peratom; // 1 if Fix saves peratom state, 0 if not int restart_file; // 1 if Fix writes own restart file, 0 if not int force_reneighbor; // 1 if Fix forces reneighboring, 0 if not int box_change_size; // 1 if Fix changes box size, 0 if not int box_change_shape; // 1 if Fix changes box shape, 0 if not int box_change_domain; // 1 if Fix changes proc sub-domains, 0 if not bigint next_reneighbor; // next timestep to force a reneighboring int thermo_energy; // 1 if fix_modify enabled ThEng, 0 if not + int thermo_virial; // 1 if fix_modify enabled ThVir, 0 if not int nevery; // how often to call an end_of_step fix int rigid_flag; // 1 if Fix integrates rigid bodies, 0 if not int peatom_flag; // 1 if Fix contributes per-atom eng, 0 if not int virial_flag; // 1 if Fix contributes to virial, 0 if not int no_change_box; // 1 if cannot swap ortho <-> triclinic int time_integrate; // 1 if fix performs time integration, 0 if no int time_depend; // 1 if requires continuous timestepping int create_attribute; // 1 if fix stores attributes that need // setting when a new atom is created int restart_pbc; // 1 if fix moves atoms (except integrate) // so write_restart must remap to PBC int wd_header; // # of header values fix writes to data file int wd_section; // # of sections fix writes to data file int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 int dof_flag; // 1 if has dof() method (not min_dof()) int special_alter_flag; // 1 if has special_alter() meth for spec lists int enforce2d_flag; // 1 if has enforce2d method int respa_level_support; // 1 if fix supports fix_modify respa int respa_level; // which respa level to apply fix (1-Nrespa) int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int array_flag; // 0/1 if compute_array() function exists int size_vector; // length of global vector int size_array_rows; // rows in global array int size_array_cols; // columns in global array int size_vector_variable; // 1 if vec length is unknown in advance int size_array_rows_variable; // 1 if array rows is unknown in advance int global_freq; // frequency s/v data is available at int peratom_flag; // 0/1 if per-atom data is stored int size_peratom_cols; // 0 = vector, N = columns in peratom array int peratom_freq; // frequency per-atom data is available at int local_flag; // 0/1 if local data is stored int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array int local_freq; // frequency local data is available at int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component int extarray; // 0/1 if global array is intensive/extensive double *vector_atom; // computed per-atom vector double **array_atom; // computed per-atom array double *vector_local; // computed local vector double **array_local; // computed local array int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int comm_border; // size of border communication (0 if none) double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial int restart_reset; // 1 if restart just re-initialized fix // KOKKOS host/device flag and data masks int kokkosable; // 1 if Kokkos fix ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; Fix(class LAMMPS *, int, char **); virtual ~Fix(); void modify_params(int, char **); virtual int setmask() = 0; virtual void post_constructor() {} virtual void init() {} virtual void init_list(int, class NeighList *) {} virtual void setup(int) {} virtual void setup_pre_exchange() {} virtual void setup_pre_neighbor() {} virtual void setup_pre_force(int) {} virtual void setup_pre_reverse(int, int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} virtual void pre_force(int) {} virtual void pre_reverse(int,int) {} virtual void post_force(int) {} virtual void final_integrate() {} virtual void end_of_step() {} virtual void post_run() {} virtual void write_restart(FILE *) {} virtual void write_restart_file(char *) {} virtual void restart(char *) {} virtual void grow_arrays(int) {} virtual void copy_arrays(int, int, int) {} virtual void set_arrays(int) {} virtual void update_arrays(int, int) {} virtual void set_molecule(int, tagint, int, double *, double *, double *) {} virtual void clear_bonus() {} virtual int pack_border(int, int *, double *) {return 0;} virtual int unpack_border(int, int, double *) {return 0;} virtual int pack_exchange(int, double *) {return 0;} virtual int unpack_exchange(int, double *) {return 0;} virtual int pack_restart(int, double *) {return 0;} virtual void unpack_restart(int, int) {} virtual int size_restart(int) {return 0;} virtual int maxsize_restart() {return 0;} virtual void setup_pre_force_respa(int, int) {} virtual void initial_integrate_respa(int, int, int) {} virtual void post_integrate_respa(int, int) {} virtual void pre_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} virtual void final_integrate_respa(int, int) {} virtual void min_pre_exchange() {} virtual void min_pre_neighbor() {} virtual void min_pre_force(int) {} virtual void min_pre_reverse(int,int) {} virtual void min_post_force(int) {} virtual double min_energy(double *) {return 0.0;} virtual void min_store() {} virtual void min_clearstore() {} virtual void min_pushstore() {} virtual void min_popstore() {} virtual int min_reset_ref() {return 0;} virtual void min_step(double, double *) {} virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm_size(int, int) {return 0;} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual double compute_scalar() {return 0.0;} virtual double compute_vector(int) {return 0.0;} virtual double compute_array(int,int) {return 0.0;} virtual int dof(int) {return 0;} virtual void deform(int) {} virtual void reset_target(double) {} virtual void reset_dt() {} virtual void enforce2d() {} virtual void read_data_header(char *) {} virtual void read_data_section(char *, int, char *, tagint) {} virtual bigint read_data_skip_lines(char *) {return 0;} virtual void write_data_header(FILE *, int) {} virtual void write_data_section_size(int, int &, int &) {} virtual void write_data_section_pack(int, double **) {} virtual void write_data_section_keyword(int, FILE *) {} virtual void write_data_section(int, FILE *, int, double **, int) {} virtual void zero_momentum() {} virtual void zero_rotation() {} virtual void rebuild_special() {} virtual int image(int *&, double **&) {return 0;} virtual int modify_param(int, char **) {return 0;} virtual void *extract(const char *, int &) {return NULL;} virtual double memory_usage() {return 0.0;} protected: int instance_me; // which Fix class instantiation I am int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; int copymode; // if set, do not deallocate during destruction // required when classes are used as functors by Kokkos int dynamic; // recount atoms for temperature computes void ev_setup(int, int); void ev_tally(int, int *, double, double, double *); void v_setup(int); void v_tally(int, int *, double, double *); + void v_tally(int, double *); + void v_tally(int, int, double); // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation union ubuf { double d; int64_t i; ubuf(double arg) : d(arg) {} ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; }; namespace FixConst { static const int INITIAL_INTEGRATE = 1<<0; static const int POST_INTEGRATE = 1<<1; static const int PRE_EXCHANGE = 1<<2; static const int PRE_NEIGHBOR = 1<<3; static const int PRE_FORCE = 1<<4; static const int PRE_REVERSE = 1<<5; static const int POST_FORCE = 1<<6; static const int FINAL_INTEGRATE = 1<<7; static const int END_OF_STEP = 1<<8; static const int POST_RUN = 1<<9; static const int THERMO_ENERGY = 1<<10; static const int INITIAL_INTEGRATE_RESPA = 1<<11; static const int POST_INTEGRATE_RESPA = 1<<12; static const int PRE_FORCE_RESPA = 1<<13; static const int POST_FORCE_RESPA = 1<<14; static const int FINAL_INTEGRATE_RESPA = 1<<15; static const int MIN_PRE_EXCHANGE = 1<<16; static const int MIN_PRE_NEIGHBOR = 1<<17; static const int MIN_PRE_FORCE = 1<<18; static const int MIN_PRE_REVERSE = 1<<19; static const int MIN_POST_FORCE = 1<<20; static const int MIN_ENERGY = 1<<21; static const int FIX_CONST_LAST = 1<<22; } } #endif /* ERROR/WARNING messages: E: Fix ID must be alphanumeric or underscore characters Self-explanatory. E: Could not find fix group ID A group ID used in the fix command does not exist. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. */