From 60cd71620a54202dfb5e6be06dd9c8618c2bb176 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 22 Jul 2015 15:20:12 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13673 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_ave_histo_weight.cpp | 562 +++++++++++++++++++++++++++++++++++ src/fix_ave_histo_weight.h | 193 ++++++++++++ 2 files changed, 755 insertions(+) create mode 100644 src/fix_ave_histo_weight.cpp create mode 100644 src/fix_ave_histo_weight.h diff --git a/src/fix_ave_histo_weight.cpp b/src/fix_ave_histo_weight.cpp new file mode 100644 index 0000000000..722e1678b3 --- /dev/null +++ b/src/fix_ave_histo_weight.cpp @@ -0,0 +1,562 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "string.h" +#include "unistd.h" +#include "fix_ave_histo_weight.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "group.h" +#include "input.h" +#include "variable.h" +#include "memory.h" +#include "error.h" +#include "force.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{X,V,F,COMPUTE,FIX,VARIABLE}; +enum{ONE,RUNNING}; +enum{SCALAR,VECTOR,WINDOW}; +enum{GLOBAL,PERATOM,LOCAL}; +enum{IGNORE,END,EXTRA}; +enum{SINGLE,VALUE}; + +#define INVOKED_SCALAR 1 +#define INVOKED_VECTOR 2 +#define INVOKED_ARRAY 4 +#define INVOKED_PERATOM 8 +#define INVOKED_LOCAL 16 + +#define BIG 1.0e20 + +/* ---------------------------------------------------------------------- */ + +FixAveHistoWeight::FixAveHistoWeight(LAMMPS *lmp, int narg, char **arg) : + FixAveHisto(lmp, narg, arg) +{ + // nvalues = 2 required for histo/weight + + if (nvalues != 2) error->all(FLERR,"Illegal fix ave/histo/weight command"); + + // check that length of 2 values is the same + + int size[2]; + + for (int i = 0; i < nvalues; i++) { + if (which[i] == X || which[i] == V || which[i] == F) { + size[i] = atom->nlocal; + } else if (which[i] == COMPUTE && kind == GLOBAL && mode == SCALAR) { + int icompute = modify->find_compute(ids[i]); + size[i] = modify->compute[icompute]->size_vector; + } else if (which[i] == COMPUTE && kind == GLOBAL && mode == VECTOR) { + int icompute = modify->find_compute(ids[i]); + size[i] = modify->compute[icompute]->size_array_rows; + } else if (which[i] == COMPUTE && kind == PERATOM) { + size[i] = atom->nlocal; + } else if (which[i] == COMPUTE && kind == LOCAL) { + int icompute = modify->find_compute(ids[i]); + size[i] = modify->compute[icompute]->size_local_rows; + } else if (which[i] == FIX && kind == GLOBAL && mode == SCALAR) { + int ifix = modify->find_fix(ids[i]); + size[i] = modify->fix[ifix]->size_vector; + } else if (which[i] == FIX && kind == GLOBAL && mode == VECTOR) { + int ifix = modify->find_fix(ids[i]); + size[i]= modify->fix[ifix]->size_array_rows; + } else if (which[i] == FIX && kind == PERATOM) { + size[i] = atom->nlocal; + } else if (which[i] == FIX && kind == LOCAL) { + int ifix = modify->find_fix(ids[i]); + size[i] = modify->fix[ifix]->size_local_rows; + } else if (which[i] == VARIABLE && kind == PERATOM) { + size[i] = atom->nlocal; + } + } + + if (size[0] != size[1]) + error->all(FLERR,"Fix ave/histo/weight value and weight vector " + "lengths do not match"); +} + +/* ---------------------------------------------------------------------- */ + +void FixAveHistoWeight::end_of_step() +{ + int i,j,m; + + // skip if not step which requires doing something + // error check if timestep was reset in an invalid manner + + bigint ntimestep = update->ntimestep; + if (ntimestep < nvalid_last || ntimestep > nvalid) + error->all(FLERR,"Invalid timestep reset for fix ave/histo"); + if (ntimestep != nvalid) return; + nvalid_last = nvalid; + + // zero if first step + + if (irepeat == 0) { + stats[0] = stats[1] = 0.0; + stats[2] = BIG; + stats[3] = -BIG; + for (i = 0; i < nbins; i++) bin[i] = 0.0; + } + + // first calculate weight factors, then bin single value + // accumulate results of computes,fixes,variables to local copy + // compute/fix/variable may invoke computes so wrap with clear/add + + modify->clearstep_compute(); + + // calcualte weight factors which are 2nd value (i = 1) + + double weight = 0.0; + double *weights = NULL; + int stride = 0; + i = 1; + + m = value2index[i]; + j = argindex[i]; + + // atom attributes + + if (which[i] == X) { + weights = &atom->x[0][j]; + stride = 3; + } else if (which[i] == V){ + weights = &atom->v[0][j]; + stride = 3; + bin_atoms(&atom->v[0][j],3); + } else if (which[i] == F) { + weights = &atom->f[0][j]; + stride = 3; + } + + // invoke compute if not previously invoked + + if (which[i] == COMPUTE) { + Compute *compute = modify->compute[m]; + if (kind == GLOBAL && mode == SCALAR) { + if (j == 0) { + if (!(compute->invoked_flag & INVOKED_SCALAR)) { + compute->compute_scalar(); + compute->invoked_flag |= INVOKED_SCALAR; + } + weight = compute->scalar; + } else { + if (!(compute->invoked_flag & INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= INVOKED_VECTOR; + } + weight = compute->vector[j-1]; + } + } else if (kind == GLOBAL && mode == VECTOR) { + if (j == 0) { + if (!(compute->invoked_flag & INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= INVOKED_VECTOR; + } + weights = compute->vector; + stride = 1; + } else { + if (!(compute->invoked_flag & INVOKED_ARRAY)) { + compute->compute_array(); + compute->invoked_flag |= INVOKED_ARRAY; + } + if (compute->array) weights = &compute->array[0][j-1]; + stride = compute->size_array_cols; + } + } else if (kind == PERATOM) { + if (!(compute->invoked_flag & INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= INVOKED_PERATOM; + } + if (j == 0) { + weights = compute->vector_atom; + stride = 1; + } else if (compute->array_atom) { + weights = &compute->array_atom[0][j-1]; + stride = compute->size_peratom_cols; + } + } else if (kind == LOCAL) { + if (!(compute->invoked_flag & INVOKED_LOCAL)) { + compute->compute_local(); + compute->invoked_flag |= INVOKED_LOCAL; + } + if (j == 0) { + weights = compute->vector_local; + stride = 1; + } else if (compute->array_local) { + weights = &compute->array_local[0][j-1]; + stride = compute->size_local_cols; + } + } + + // access fix fields, guaranteed to be ready + + } else if (which[i] == FIX) { + + Fix *fix = modify->fix[m]; + + if (kind == GLOBAL && mode == SCALAR) { + if (j == 0) weight = fix->compute_scalar(); + else weight = fix->compute_vector(j-1); + + } else if (kind == GLOBAL && mode == VECTOR) { + + error->all(FLERR,"Illegal fix ave/spatial command"); + + if (j == 0) { + int n = fix->size_vector; + for (i = 0; i < n; i++) weights[n] = fix->compute_vector(i); + } else { + int n = fix->size_vector; + for (i = 0; i < n; i++) weights[n] = fix->compute_array(i,j-1); + } + + } else if (kind == PERATOM) { + if (j == 0) { + weights = fix->vector_atom; + stride = 1; + } else if (fix->array_atom) { + weights = fix->array_atom[j-1]; + stride = fix->size_peratom_cols; + } + } else if (kind == LOCAL) { + if (j == 0) { + weights = fix->vector_local; + stride = 1; + } else if (fix->array_local) { + weights = &fix->array_local[0][j-1]; + stride = fix->size_local_cols; + } + } + + // evaluate equal-style variable + + } else if (which[i] == VARIABLE && kind == GLOBAL) { + weight = input->variable->compute_equal(m); + + } else if (which[i] == VARIABLE && kind == PERATOM) { + if (atom->nlocal > maxatom) { + memory->destroy(vector); + maxatom = atom->nmax; + memory->create(vector,maxatom,"ave/histo/weights:vector"); + } + input->variable->compute_atom(m,igroup,vector,1,0); + weights = vector; + stride = 1; + } + + // bin values using weights, values are 1st value (i = 0) + + i = 0; + m = value2index[i]; + j = argindex[i]; + + // atom attributes + + if (which[i] == X) + bin_atoms_weights(&atom->x[0][j],3,weights,stride); + else if (which[i] == V) + bin_atoms_weights(&atom->v[0][j],3,weights,stride); + else if (which[i] == F) + bin_atoms_weights(&atom->f[0][j],3,weights,stride); + + // invoke compute if not previously invoked + + if (which[i] == COMPUTE) { + Compute *compute = modify->compute[m]; + if (kind == GLOBAL && mode == SCALAR) { + if (j == 0) { + if (!(compute->invoked_flag & INVOKED_SCALAR)) { + compute->compute_scalar(); + compute->invoked_flag |= INVOKED_SCALAR; + } + bin_one_weights(compute->scalar,weight); + } else { + if (!(compute->invoked_flag & INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= INVOKED_VECTOR; + } + bin_one_weights(compute->vector[j-1],weight); + } + } else if (kind == GLOBAL && mode == VECTOR) { + if (j == 0) { + if (!(compute->invoked_flag & INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= INVOKED_VECTOR; + } + bin_vector_weights(compute->size_vector,compute->vector,1, + weights,stride); + } else { + if (!(compute->invoked_flag & INVOKED_ARRAY)) { + compute->compute_array(); + compute->invoked_flag |= INVOKED_ARRAY; + } + if (compute->array) + bin_vector_weights(compute->size_array_rows,&compute->array[0][j-1], + compute->size_array_cols,weights,stride); + } + + } else if (kind == PERATOM) { + if (!(compute->invoked_flag & INVOKED_PERATOM)) { + compute->compute_peratom(); + compute->invoked_flag |= INVOKED_PERATOM; + } + if (j == 0) + bin_atoms_weights(compute->vector_atom,1,weights, stride); + else if (compute->array_atom) + bin_atoms_weights(&compute->array_atom[0][j-1], + compute->size_peratom_cols,weights,stride); + + } else if (kind == LOCAL) { + if (!(compute->invoked_flag & INVOKED_LOCAL)) { + compute->compute_local(); + compute->invoked_flag |= INVOKED_LOCAL; + } + if (j == 0) + bin_vector_weights(compute->size_local_rows, + compute->vector_local,1,weights,stride); + else if (compute->array_local) + bin_vector_weights(compute->size_local_rows, + &compute->array_local[0][j-1], + compute->size_local_cols,weights,stride); + } + + // access fix fields, guaranteed to be ready + + } else if (which[i] == FIX) { + + Fix *fix = modify->fix[m]; + + if (kind == GLOBAL && mode == SCALAR) { + if (j == 0) bin_one_weights(fix->compute_scalar(),weight); + else bin_one_weights(fix->compute_vector(j-1),weight); + + } else if (kind == GLOBAL && mode == VECTOR) { + if (j == 0) { + int n = fix->size_vector; + for (i = 0; i < n; i++) + bin_one_weights(fix->compute_vector(i),weights[i*stride]); + } else { + int n = fix->size_vector; + for (i = 0; i < n; i++) + bin_one_weights(fix->compute_array(i,j-1),weights[i*stride]); + } + + } else if (kind == PERATOM) { + if (j == 0) + bin_atoms_weights(fix->vector_atom,1,weights,stride); + else if (fix->array_atom) + bin_atoms_weights(fix->array_atom[j-1],fix->size_peratom_cols, + weights,stride); + + + } else if (kind == LOCAL) { + if (j == 0) bin_vector_weights(fix->size_local_rows,fix->vector_local,1, + weights,stride); + else if (fix->array_local) + bin_vector_weights(fix->size_local_rows,&fix->array_local[0][j-1], + fix->size_local_cols,weights,stride); + } + + // evaluate equal-style variable + + } else if (which[i] == VARIABLE && kind == GLOBAL) { + bin_one_weights(input->variable->compute_equal(m),weight); + + } else if (which[i] == VARIABLE && kind == PERATOM) { + if (atom->nlocal > maxatom) { + memory->destroy(vector); + maxatom = atom->nmax; + memory->create(vector,maxatom,"ave/histo/weights:vector"); + } + input->variable->compute_atom(m,igroup,vector,1,0); + bin_atoms_weights(vector,1,weights,stride); + } + + // code beyond this point is identical to FixAveHisto + + // done if irepeat < nrepeat + // else reset irepeat and nvalid + + irepeat++; + if (irepeat < nrepeat) { + nvalid += nevery; + modify->addstep_compute(nvalid); + return; + } + + irepeat = 0; + nvalid = ntimestep + nfreq - (nrepeat-1)*nevery; + modify->addstep_compute(nvalid); + + // merge histogram stats across procs if necessary + + if (kind == PERATOM || kind == LOCAL) { + MPI_Allreduce(stats,stats_all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&stats[2],&stats_all[2],1,MPI_DOUBLE,MPI_MIN,world); + MPI_Allreduce(&stats[3],&stats_all[3],1,MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(bin,bin_all,nbins,MPI_DOUBLE,MPI_SUM,world); + + stats[0] = stats_all[0]; + stats[1] = stats_all[1]; + stats[2] = stats_all[2]; + stats[3] = stats_all[3]; + for (i = 0; i < nbins; i++) bin[i] = bin_all[i]; + } + + // if ave = ONE, only single Nfreq timestep value is needed + // if ave = RUNNING, combine with all previous Nfreq timestep values + // if ave = WINDOW, combine with nwindow most recent Nfreq timestep values + + if (ave == ONE) { + stats_total[0] = stats[0]; + stats_total[1] = stats[1]; + stats_total[2] = stats[2]; + stats_total[3] = stats[3]; + for (i = 0; i < nbins; i++) bin_total[i] = bin[i]; + + } else if (ave == RUNNING) { + stats_total[0] += stats[0]; + stats_total[1] += stats[1]; + stats_total[2] = MIN(stats_total[2],stats[2]); + stats_total[3] = MAX(stats_total[3],stats[3]); + for (i = 0; i < nbins; i++) bin_total[i] += bin[i]; + + } else if (ave == WINDOW) { + stats_total[0] += stats[0]; + if (window_limit) stats_total[0] -= stats_list[iwindow][0]; + stats_list[iwindow][0] = stats[0]; + stats_total[1] += stats[1]; + if (window_limit) stats_total[1] -= stats_list[iwindow][1]; + stats_list[iwindow][1] = stats[1]; + + if (window_limit) m = nwindow; + else m = iwindow+1; + + stats_list[iwindow][2] = stats[2]; + stats_total[2] = stats_list[0][2]; + for (i = 1; i < m; i++) + stats_total[2] = MIN(stats_total[2],stats_list[i][2]); + stats_list[iwindow][3] = stats[3]; + stats_total[3] = stats_list[0][3]; + for (i = 1; i < m; i++) + stats_total[3] = MAX(stats_total[3],stats_list[i][3]); + + for (i = 0; i < nbins; i++) { + bin_total[i] += bin[i]; + if (window_limit) bin_total[i] -= bin_list[iwindow][i]; + bin_list[iwindow][i] = bin[i]; + } + + iwindow++; + if (iwindow == nwindow) { + iwindow = 0; + window_limit = 1; + } + } + + // output result to file + + if (fp && me == 0) { + if (overwrite) fseek(fp,filepos,SEEK_SET); + fprintf(fp,BIGINT_FORMAT " %d %g %g %g %g\n",ntimestep,nbins, + stats_total[0],stats_total[1],stats_total[2],stats_total[3]); + if (stats_total[0] != 0.0) + for (i = 0; i < nbins; i++) + fprintf(fp,"%d %g %g %g\n", + i+1,coord[i],bin_total[i],bin_total[i]/stats_total[0]); + else + for (i = 0; i < nbins; i++) + fprintf(fp,"%d %g %g %g\n",i+1,coord[i],0.0,0.0); + fflush(fp); + if (overwrite) { + long fileend = ftell(fp); + ftruncate(fileno(fp),fileend); + } + } +} + +/* ---------------------------------------------------------------------- + bin a single value with weight) +------------------------------------------------------------------------- */ + +void FixAveHistoWeight::bin_one_weights(double value, double weight) +{ + stats[2] = MIN(stats[2],value); + stats[3] = MAX(stats[3],value); + + if (value < lo) { + if (beyond == IGNORE) { + stats[1] += weight; + return; + } else bin[0] += weight; + } else if (value > hi) { + if (beyond == IGNORE) { + stats[1] += weight; + return; + } else bin[nbins-1] += weight; + } else { + int ibin = static_cast ((value-lo)*bininv); + ibin = MIN(ibin,nbins-1); + if (beyond == EXTRA) ibin++; + bin[ibin] += weight; + } + + stats[0] += weight; +} + +/* ---------------------------------------------------------------------- + bin a vector of values with weights + values and weights each have a stride +------------------------------------------------------------------------- */ + +void FixAveHistoWeight::bin_vector_weights(int n, double *values, + int stride, double *weights, + int stridewt) +{ + int m = 0; + int m2 = 0; + for (int i = 0; i < n; i++) { + bin_one_weights(values[m],weights[m2]); + m += stride; + m2 += stridewt; + } +} + +/* ---------------------------------------------------------------------- + bin a per-atom vector of values with weights + values and weights each have a stride + only bin if atom is in group +------------------------------------------------------------------------- */ + +void FixAveHistoWeight::bin_atoms_weights(double *values, int stride, + double *weights, int stridewt) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int m = 0; + int m2 = 0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) bin_one_weights(values[m],weights[m2]); + m += stride; + m2 += stridewt; + } +} diff --git a/src/fix_ave_histo_weight.h b/src/fix_ave_histo_weight.h new file mode 100644 index 0000000000..99db7333d8 --- /dev/null +++ b/src/fix_ave_histo_weight.h @@ -0,0 +1,193 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(ave/histo/weight,FixAveHistoWeight) + +#else + +#ifndef LMP_FIX_AVE_HISTO_WEIGHT_H +#define LMP_FIX_AVE_HISTO_WEIGHT_H + +#include "stdio.h" +#include "fix_ave_histo.h" + +namespace LAMMPS_NS { + +class FixAveHistoWeight : public FixAveHisto { + public: + FixAveHistoWeight(class LAMMPS *, int, char **); + ~FixAveHistoWeight() {} + void end_of_step(); + + private: + void bin_one_weights(double, double); + void bin_vector_weights(int, double *, int, double *, int); + void bin_atoms_weights(double *, int, double *, int); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +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. + +E: Compute ID for fix ave/histo does not exist + +Self-explanatory. + +E: Fix ID for fix ave/histo does not exist + +Self-explanatory. + +E: Fix ave/histo input is invalid compute + +Self-explanatory. + +E: Fix ave/histo input is invalid fix + +Self-explanatory. + +E: Fix ave/histo input is invalid variable + +Self-explanatory. + +E: Fix ave/histo inputs are not all global, peratom, or local + +All inputs in a single fix ave/histo command must be of the +same style. + +E: Fix ave/histo cannot input per-atom values in scalar mode + +Self-explanatory. + +E: Fix ave/histo cannot input local values in scalar mode + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a global scalar + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a global vector + +Self-explanatory. + +E: Fix ave/histo compute vector is accessed out-of-range + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a global array + +Self-explanatory. + +E: Fix ave/histo compute array is accessed out-of-range + +Self-explanatory. + +E: Fix ave/histo compute does not calculate per-atom values + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a per-atom vector + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a per-atom array + +Self-explanatory. + +E: Fix ave/histo compute does not calculate local values + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a local vector + +Self-explanatory. + +E: Fix ave/histo compute does not calculate a local array + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a global scalar + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a global vector + +Self-explanatory. + +E: Fix ave/histo fix vector is accessed out-of-range + +Self-explanatory. + +E: Fix for fix ave/histo not computed at compatible time + +Fixes generate their values on specific timesteps. Fix ave/histo is +requesting a value on a non-allowed timestep. + +E: Fix ave/histo fix does not calculate a global array + +Self-explanatory. + +E: Fix ave/histo fix array is accessed out-of-range + +Self-explanatory. + +E: Fix ave/histo fix does not calculate per-atom values + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a per-atom vector + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a per-atom array + +Self-explanatory. + +E: Fix ave/histo fix does not calculate local values + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a local vector + +Self-explanatory. + +E: Fix ave/histo fix does not calculate a local array + +Self-explanatory. + +E: Variable name for fix ave/histo does not exist + +Self-explanatory. + +E: Invalid timestep reset for fix ave/histo + +Resetting the timestep has invalidated the sequence of timesteps this +fix needs to process. + +E: Cannot open fix ave/histo file %s + +The specified file cannot be opened. Check that the path and name are +correct. + +*/