C++ Compute::compute_vector方法代码示例

void *lammps_extract_compute(void *ptr, char *id, int style, int type)
{
  LAMMPS *lmp = (LAMMPS *) ptr;

  int icompute = lmp->modify->find_compute(id);
  if (icompute < 0) return NULL;
  Compute *compute = lmp->modify->compute[icompute];

  if (style == 0) {
    if (type == 0) {
      if (!compute->scalar_flag) return NULL;
      if (compute->invoked_scalar != lmp->update->ntimestep)
	compute->compute_scalar();
      return (void *) &compute->scalar;
    }
    if (type == 1) {
      if (!compute->vector_flag) return NULL;
      if (compute->invoked_vector != lmp->update->ntimestep)
	compute->compute_vector();
      return (void *) compute->vector;
    }
    if (type == 2) {
      if (!compute->array_flag) return NULL;
      if (compute->invoked_array != lmp->update->ntimestep)
	compute->compute_array();
      return (void *) compute->array;
    }
  }

  if (style == 1) {
    if (!compute->peratom_flag) return NULL;
    if (type == 1) {
      if (compute->invoked_peratom != lmp->update->ntimestep)
	compute->compute_peratom();
      return (void *) compute->vector_atom;
    }
    if (type == 2) {
      if (compute->invoked_peratom != lmp->update->ntimestep)
	compute->compute_peratom();
      return (void *) compute->array_atom;
    }
  }

  if (style == 2) {
    if (!compute->local_flag) return NULL;
    if (type == 1) {
      if (compute->invoked_local != lmp->update->ntimestep)
	compute->compute_local();
      return (void *) compute->vector_local;
    }
    if (type == 2) {
      if (compute->invoked_local != lmp->update->ntimestep)
	compute->compute_local();
      return (void *) compute->array_local;
    }
  }

  return NULL;
}

void ComputeSlice::extract_one(int m, double *vec, int stride)
{
  int i,j;

  // invoke the appropriate compute if needed

  if (which[m] == COMPUTE) {
    Compute *compute = modify->compute[value2index[m]];

    if (argindex[m] == 0) {
      if (!(compute->invoked_flag & INVOKED_VECTOR)) {
	compute->compute_vector();
	compute->invoked_flag |= INVOKED_VECTOR;
      }
      double *cvector = compute->vector;
      j = 0;
      for (i = nstart; i < nstop; i += nskip) {
	vec[j] = cvector[i-1];
	j += stride;
      }
      
    } else {
      if (!(compute->invoked_flag & INVOKED_ARRAY)) {
	compute->compute_array();
	compute->invoked_flag |= INVOKED_ARRAY;
      }
      double **carray = compute->array;
      int icol = argindex[m]-1;
      j = 0;
      for (i = nstart; i < nstop; i += nskip) {
	vec[j] = carray[i-1][icol];
	j += stride;
      }
    }

  // access fix fields, check if fix frequency is a match
    
  } else if (which[m] == FIX) {
    if (update->ntimestep % modify->fix[value2index[m]]->global_freq)
      error->all(FLERR,"Fix used in compute slice not computed at compatible time");
    Fix *fix = modify->fix[value2index[m]];

    if (argindex[m] == 0) {
      j = 0;
      for (i = nstart; i < nstop; i += nskip) {
	vec[j] = fix->compute_vector(i-1);
	j += stride;
      }
    } else {
      int icol = argindex[m]-1;
      j = 0;
      for (i = nstart; i < nstop; i += nskip) {
	vec[j] = fix->compute_array(i-1,icol);
	j += stride;
      }
    }
  }
}

void FixAveHisto::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;
  }

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (i = 0; i < nvalues; i++) {
    m = value2index[i];
    j = argindex[i];

    // atom attributes

    if (which[i] == X)
      bin_atoms(&atom->x[0][j],3);
    else if (which[i] == V)
      bin_atoms(&atom->v[0][j],3);
    else if (which[i] == F)
      bin_atoms(&atom->f[0][j],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;
          }
          bin_one(compute->scalar);
        } else {
          if (!(compute->invoked_flag & INVOKED_VECTOR)) {
            compute->compute_vector();
            compute->invoked_flag |= INVOKED_VECTOR;
          }
          bin_one(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;
          }
          bin_vector(compute->size_vector,compute->vector,1);
        } else {
          if (!(compute->invoked_flag & INVOKED_ARRAY)) {
            compute->compute_array();
            compute->invoked_flag |= INVOKED_ARRAY;
          }
          if (compute->array)
            bin_vector(compute->size_array_rows,&compute->array[0][j-1],
                       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)
          bin_atoms(compute->vector_atom,1);
        else if (compute->array_atom)
          bin_atoms(&compute->array_atom[0][j-1],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)
          bin_vector(compute->size_local_rows,compute->vector_local,1);
        else if (compute->array_local)
          bin_vector(compute->size_local_rows,&compute->array_local[0][j-1],
                     compute->size_local_cols);
      }

      // access fix fields, guaranteed to be ready

//.........这里部分代码省略.........

void FixAveTime::invoke_vector(bigint ntimestep)
{
  int i,j,m;
  
  // zero if first step

  if (irepeat == 0)
    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) array[i][j] = 0.0;
  
  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add
  
  modify->clearstep_compute();
  
  for (j = 0; j < nvalues; j++) {
    m = value2index[j];
    
    // invoke compute if not previously invoked
    
    if (which[j] == COMPUTE) {
      Compute *compute = modify->compute[m];
      
      if (argindex[j] == 0) {
	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
	  compute->compute_vector();
	  compute->invoked_flag |= INVOKED_VECTOR;
	}
	double *cvector = compute->vector;
	for (i = 0; i < nrows; i++)
	  column[i] = cvector[i];
	
      } else {
	if (!(compute->invoked_flag & INVOKED_ARRAY)) {
	  compute->compute_array();
	  compute->invoked_flag |= INVOKED_ARRAY;
	}
	double **carray = compute->array;
	int icol = argindex[j]-1;
	for (i = 0; i < nrows; i++)
	  column[i] = carray[i][icol];
      }
      
    // access fix fields, guaranteed to be ready
      
    } else if (which[j] == FIX) {
      Fix *fix = modify->fix[m];
      if (argindex[j] == 0)
	for (i = 0; i < nrows; i++)
	  column[i] = fix->compute_vector(i);
      else {
	int icol = argindex[j]-1;
	for (i = 0; i < nrows; i++)
	  column[i] = fix->compute_array(i,icol);
      }
    }
    
    // add columns of values to array or just set directly if offcol is set
    
    if (offcol[j]) {
      for (i = 0; i < nrows; i++)
	array[i][j] = column[i];
    } else {
      for (i = 0; i < nrows; i++)
	array[i][j] += column[i];
    }
  }
  
  // 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);
  
  // average the final result for the Nfreq timestep
  
  double repeat = nrepeat;
  for (i = 0; i < nrows; i++)
    for (j = 0; j < nvalues; j++)
      if (offcol[j] == 0) array[i][j] /= repeat;
  
  // 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 (ntimestep >= startstep) {
    if (ave == ONE) {
      for (i = 0; i < nrows; i++)
	for (j = 0; j < nvalues; j++) array_total[i][j] = array[i][j];
      norm = 1;
      
    } else if (ave == RUNNING) {
//.........这里部分代码省略.........

void FixAveTime::invoke_scalar(bigint ntimestep)
{
  int i,m;
  double scalar;

  // zero if first step

  if (irepeat == 0)
    for (i = 0; i < nvalues; i++) vector[i] = 0.0;

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (i = 0; i < nvalues; i++) {
    m = value2index[i];
    
    // invoke compute if not previously invoked
    
    if (which[i] == COMPUTE) {
      Compute *compute = modify->compute[m];
      
      if (argindex[i] == 0) {
	if (!(compute->invoked_flag & INVOKED_SCALAR)) {
	  compute->compute_scalar();
	  compute->invoked_flag |= INVOKED_SCALAR;
	}
	scalar = compute->scalar;
      } else {
	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
	  compute->compute_vector();
	  compute->invoked_flag |= INVOKED_VECTOR;
	}
	scalar = compute->vector[argindex[i]-1];
      }
      
    // access fix fields, guaranteed to be ready
      
    } else if (which[i] == FIX) {
      if (argindex[i] == 0) 
	scalar = modify->fix[m]->compute_scalar();
      else
	scalar = modify->fix[m]->compute_vector(argindex[i]-1);
      
    // evaluate equal-style variable
      
    } else if (which[i] == VARIABLE)
      scalar = input->variable->compute_equal(m);
    
    // add value to vector or just set directly if offcol is set
    
    if (offcol[i]) vector[i] = scalar;
    else vector[i] += scalar;
  }

  // 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);

  // average the final result for the Nfreq timestep

  double repeat = nrepeat;
  for (i = 0; i < nvalues; i++)
    if (offcol[i] == 0) vector[i] /= repeat;

  // 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 (ntimestep >= startstep) {
    if (ave == ONE) {
      for (i = 0; i < nvalues; i++) vector_total[i] = vector[i];
      norm = 1;

    } else if (ave == RUNNING) {
      for (i = 0; i < nvalues; i++) vector_total[i] += vector[i];
      norm++;
      
    } else if (ave == WINDOW) {
      for (i = 0; i < nvalues; i++) {
	vector_total[i] += vector[i];
	if (window_limit) vector_total[i] -= vector_list[iwindow][i];
	vector_list[iwindow][i] = vector[i];
      }
      
      iwindow++;
      if (iwindow == nwindow) {
	iwindow = 0;
	window_limit = 1;
//.........这里部分代码省略.........

void FixAveTime::invoke_vector(bigint ntimestep)
{
  int i,j,m;

  // first sample within single Nfreq epoch
  // zero out arrays that accumulate over many samples, but not across epochs
  // invoke setup_chunks() to determine current nchunk
  //   re-allocate per-chunk arrays if needed
  // invoke lock() in two cases:
  //   if nrepeat > 1: so nchunk cannot change until Nfreq epoch is over,
  //     will be unlocked on last repeat of this Nfreq
  //   if ave = RUNNING/WINDOW and not yet locked:
  //     set forever, will be unlocked in fix destructor
  // wrap setup_chunks in clearstep/addstep b/c it may invoke computes
  //   both nevery and nfreq are future steps,
  //   since call below to cchunk->ichunk()
  //     does not re-invoke internal cchunk compute on this same step

  if (irepeat == 0) {
    if (any_variable_length) {
      modify->clearstep_compute();
      int nrows_new = column_length(1);
      modify->addstep_compute(ntimestep+nevery);
      modify->addstep_compute(ntimestep+nfreq);

      if (all_variable_length && nrows_new != nrows) {
        nrows = nrows_new;
        memory->destroy(column);
        memory->create(column,nrows,"ave/time:column");
        allocate_arrays();
      }

      bigint ntimestep = update->ntimestep;
      int lockforever_flag = 0;
      for (i = 0; i < nvalues; i++) {
        if (!varlen[i]) continue;
        if (nrepeat > 1 && ave == ONE) {
          Compute *compute = modify->compute[value2index[i]];
          compute->lock(this,ntimestep,ntimestep+(nrepeat-1)*nevery);
        } else if ((ave == RUNNING || ave == WINDOW) && !lockforever) {
          Compute *compute = modify->compute[value2index[i]];
          compute->lock(this,update->ntimestep,-1);
          lockforever_flag = 1;
        }
      }
      if (lockforever_flag) lockforever = 1;
    }

    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) array[i][j] = 0.0;
  }

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (j = 0; j < nvalues; j++) {
    m = value2index[j];

    // invoke compute if not previously invoked

    if (which[j] == COMPUTE) {
      Compute *compute = modify->compute[m];

      if (argindex[j] == 0) {
        if (!(compute->invoked_flag & INVOKED_VECTOR)) {
          compute->compute_vector();
          compute->invoked_flag |= INVOKED_VECTOR;
        }
        double *cvector = compute->vector;
        for (i = 0; i < nrows; i++)
          column[i] = cvector[i];

      } else {
        if (!(compute->invoked_flag & INVOKED_ARRAY)) {
          compute->compute_array();
          compute->invoked_flag |= INVOKED_ARRAY;
        }
        double **carray = compute->array;
        int icol = argindex[j]-1;
        for (i = 0; i < nrows; i++)
          column[i] = carray[i][icol];
      }

    // access fix fields, guaranteed to be ready

    } else if (which[j] == FIX) {
      Fix *fix = modify->fix[m];
      if (argindex[j] == 0)
        for (i = 0; i < nrows; i++)
          column[i] = fix->compute_vector(i);
      else {
        int icol = argindex[j]-1;
        for (i = 0; i < nrows; i++)
          column[i] = fix->compute_array(i,icol);
      }
    }

    // add columns of values to array or just set directly if offcol is set
//.........这里部分代码省略.........

void FixAveTime::invoke_scalar(bigint ntimestep)
{
  int i,m;
  double scalar;

  // zero if first sample within single Nfreq epoch
  // NOTE: doc this
  // are not checking for returned length, just initialize it
  // check for exceeding length is done below

  if (irepeat == 0) {
    if (any_variable_length) {
      modify->clearstep_compute();
      column_length(1);
      modify->addstep_compute(ntimestep+nevery);
      modify->addstep_compute(ntimestep+nfreq);
    }
    for (i = 0; i < nvalues; i++) vector[i] = 0.0;
  }

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (i = 0; i < nvalues; i++) {
    m = value2index[i];

    // invoke compute if not previously invoked

    if (which[i] == COMPUTE) {
      Compute *compute = modify->compute[m];

      if (argindex[i] == 0) {
        if (!(compute->invoked_flag & INVOKED_SCALAR)) {
          compute->compute_scalar();
          compute->invoked_flag |= INVOKED_SCALAR;
        }
        scalar = compute->scalar;
      } else {
        if (!(compute->invoked_flag & INVOKED_VECTOR)) {
          compute->compute_vector();
          compute->invoked_flag |= INVOKED_VECTOR;
        }

        // insure no out-of-range access to variable-length compute vector

        if (varlen[i] && compute->size_vector < argindex[i]) scalar = 0.0;
        else scalar = compute->vector[argindex[i]-1];
      }

    // access fix fields, guaranteed to be ready

    } else if (which[i] == FIX) {
      if (argindex[i] == 0)
        scalar = modify->fix[m]->compute_scalar();
      else
        scalar = modify->fix[m]->compute_vector(argindex[i]-1);

    // evaluate equal-style variable

    } else if (which[i] == VARIABLE)
      scalar = input->variable->compute_equal(m);

    // add value to vector or just set directly if offcol is set

    if (offcol[i]) vector[i] = scalar;
    else vector[i] += scalar;
  }

  // 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);

  // average the final result for the Nfreq timestep

  double repeat = nrepeat;
  for (i = 0; i < nvalues; i++)
    if (offcol[i] == 0) vector[i] /= repeat;

  // 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) {
    for (i = 0; i < nvalues; i++) vector_total[i] = vector[i];
    norm = 1;

  } else if (ave == RUNNING) {
    for (i = 0; i < nvalues; i++) vector_total[i] += vector[i];
//.........这里部分代码省略.........

void ComputeChunkSpreadAtom::compute_peratom()
{
  invoked_peratom = update->ntimestep;

  // grow local vector_atom or array_atom if necessary

  if (atom->nmax > nmax) {
    if (nvalues == 1) {
      memory->destroy(vector_atom);
      nmax = atom->nmax;
      memory->create(vector_atom,nmax,"chunk/spread/atom:vector_atom");
    } else {
      memory->destroy(array_atom);
      nmax = atom->nmax;
      memory->create(array_atom,nmax,nvalues,"chunk/spread/atom:array_atom");
    }
  }

  // compute chunk/atom assigns atoms to chunk IDs
  // extract ichunk index vector from compute
  // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms

  int nchunk = cchunk->setup_chunks();
  cchunk->compute_ichunk();
  int *ichunk = cchunk->ichunk;

  // loop over values, access compute or fix
  // loop over atoms, use chunk ID of each atom to store value from compute/fix

  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  int i,m,n,index,nstride;
  double *ptr;

  for (m = 0; m < nvalues; m++) {
    n = value2index[m];

    // copy compute/fix values into vector_atom or array_atom
    // nstride between values for each atom

    if (nvalues == 1) {
      ptr = vector_atom;
      nstride = 1;
    } else {
      ptr = &array_atom[0][m];
      nstride = nvalues;
    }

    // invoke compute if not previously invoked

    if (which[m] == COMPUTE) {
      Compute *compute = modify->compute[n];

      if (argindex[m] == 0) {
        if (!(compute->invoked_flag & INVOKED_VECTOR)) {
          compute->compute_vector();
          compute->invoked_flag |= INVOKED_VECTOR;
        }
	double *cvector = compute->vector;
	for (i = 0; i < nlocal; i++, ptr += nstride) {
	  *ptr = 0.0;
	  if (!(mask[i] & groupbit)) continue;
	  index = ichunk[i]-1;
	  if (index < 0 || index >= nchunk) continue;
	  *ptr = cvector[index];
	}

      } else {
        if (!(compute->invoked_flag & INVOKED_ARRAY)) {
          compute->compute_array();
          compute->invoked_flag |= INVOKED_ARRAY;
        }
        int icol = argindex[m]-1;
        double **carray = compute->array;
	for (i = 0; i < nlocal; i++, ptr += nstride) {
	  *ptr = 0.0;
	  if (!(mask[i] & groupbit)) continue;
	  index = ichunk[i]-1;
	  if (index < 0 || index >= nchunk) continue;
	  *ptr = carray[index][icol];
	}
      }

    // access fix data, check if fix frequency is a match
    // are assuming the fix global vector/array is per-chunk data
    // check if index exceeds fix output length/rows

    } else if (which[m] == FIX) {
      Fix *fix = modify->fix[n];
      if (update->ntimestep % fix->global_freq)
        error->all(FLERR,"Fix used in compute chunk/spread/atom not "
                   "computed at compatible time");

      if (argindex[m] == 0) {
        int nfix = fix->size_vector;
        for (i = 0; i < nlocal; i++, ptr += nstride) {
          *ptr = 0.0;
          if (!(mask[i] & groupbit)) continue;
          index = ichunk[i]-1;
//.........这里部分代码省略.........

void FixAveCorrelate::end_of_step()
{
  int i,j,m;
  double scalar;

  // 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/correlate");
  if (ntimestep != nvalid) return;
  nvalid_last = nvalid;

  // accumulate results of computes,fixes,variables to origin
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  // lastindex = index in values ring of latest time sample

  lastindex++;
  if (lastindex == nrepeat) lastindex = 0;

  for (i = 0; i < nvalues; i++) {
    m = value2index[i];

    // invoke compute if not previously invoked

    if (which[i] == COMPUTE) {
      Compute *compute = modify->compute[m];

      if (argindex[i] == 0) {
        if (!(compute->invoked_flag & INVOKED_SCALAR)) {
          compute->compute_scalar();
          compute->invoked_flag |= INVOKED_SCALAR;
        }
        scalar = compute->scalar;
      } else {
        if (!(compute->invoked_flag & INVOKED_VECTOR)) {
          compute->compute_vector();
          compute->invoked_flag |= INVOKED_VECTOR;
        }
        scalar = compute->vector[argindex[i]-1];
      }

    // access fix fields, guaranteed to be ready

    } else if (which[i] == FIX) {
      if (argindex[i] == 0)
        scalar = modify->fix[m]->compute_scalar();
      else
        scalar = modify->fix[m]->compute_vector(argindex[i]-1);

    // evaluate equal-style variable

    } else if (which[i] == VARIABLE)
      scalar = input->variable->compute_equal(m);

    values[lastindex][i] = scalar;
  }

  // fistindex = index in values ring of earliest time sample
  // nsample = number of time samples in values ring

  if (nsample < nrepeat) nsample++;
  else {
    firstindex++;
    if (firstindex == nrepeat) firstindex = 0;
  }

  nvalid += nevery;
  modify->addstep_compute(nvalid);

  // calculate all Cij() enabled by latest values

  accumulate();
  if (ntimestep % nfreq) return;

  // save results in save_count and save_corr

  for (i = 0; i < nrepeat; i++) {
    save_count[i] = count[i];
    if (count[i])
      for (j = 0; j < npair; j++)
        save_corr[i][j] = prefactor*corr[i][j]/count[i];
    else
      for (j = 0; j < npair; j++)
        save_corr[i][j] = 0.0;
  }

  // output result to file

  if (fp && me == 0) {
    if (overwrite) fseek(fp,filepos,SEEK_SET);
    fprintf(fp,BIGINT_FORMAT " %d\n",ntimestep,nrepeat);
    for (i = 0; i < nrepeat; i++) {
      fprintf(fp,"%d %d %d",i+1,i*nevery,count[i]);
      if (count[i])
        for (j = 0; j < npair; j++)
//.........这里部分代码省略.........

void FixSAEDVTK::invoke_vector(bigint ntimestep)
{
  // zero if first step
  int icompute = modify->find_compute(ids);
  if (icompute < 0)
    error->all(FLERR,"Compute ID for fix saed/vtk does not exist");

  if (irepeat == 0)
    for (int i = 0; i < nrows; i++)
       vector[i] = 0.0;

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  // invoke compute if not previously invoked

  Compute *compute = modify->compute[icompute];

  if (!(compute->invoked_flag & INVOKED_VECTOR)) {
    compute->compute_vector();
    compute->invoked_flag |= INVOKED_VECTOR;
  }

  double *vector = compute->vector;

  // 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);

  // average the final result for the Nfreq timestep

  double repeat = nrepeat;
  for ( int i = 0; i < nrows; i++)
    vector[i] /= repeat;

  // 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) {
    for (int i = 0; i < nrows; i++) vector_total[i] = vector[i];
    norm = 1;

  } else if (ave == RUNNING) {
    for (int i = 0; i < nrows; i++) vector_total[i] += vector[i];
    norm++;

  } else if (ave == WINDOW) {
    for (int i = 0; i < nrows; i++) {
      vector_total[i] += vector[i];
      if (window_limit) vector_total[i] -= vector_list[iwindow][i];
        vector_list[iwindow][i] = vector[i];
    }

    iwindow++;
    if (iwindow == nwindow) {
      iwindow = 0;
      window_limit = 1;
    }
    if (window_limit) norm = nwindow;
    else norm = iwindow;
  }

  // output result to file

  if (fp && me == 0) {
    if (nOutput > 0) {
      fclose(fp);

      char nName [128];
      sprintf(nName,"%s.%d.vtk",filename,nOutput);
      fp = fopen(nName,"w");

      if (fp == NULL) {
        char str[128];
        sprintf(str,"Cannot open fix saed/vtk file %s",nName);
        error->one(FLERR,str);
      }
    }

    fprintf(fp,"# vtk DataFile Version 3.0 c_%s\n",ids);
    fprintf(fp,"Image data set\n");
    fprintf(fp,"ASCII\n");
    fprintf(fp,"DATASET STRUCTURED_POINTS\n");
    fprintf(fp,"DIMENSIONS %d %d %d\n", Dim[0],  Dim[1], Dim[2]);
    fprintf(fp,"ASPECT_RATIO %g %g %g\n", dK[0], dK[1], dK[2]);
    fprintf(fp,"ORIGIN %g %g %g\n", Knmin[0] * dK[0],  Knmin[1] * dK[1], Knmin[2] * dK[2]);
    fprintf(fp,"POINT_DATA %d\n",  Dim[0] *  Dim[1] * Dim[2] );
//.........这里部分代码省略.........

void FixAveTimeSAED::invoke_vector(bigint ntimestep)
{

  int i,j,m;
  // zero if first step

  if (irepeat == 0)
    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) array[i][j] = 0.0;

  // accumulate results of computes,fixes,variables to local copy
  // compute/fix/variable may invoke computes so wrap with clear/add

  modify->clearstep_compute();

  for (j = 0; j < nvalues; j++) {
    m = value2index[j];

    // invoke compute if not previously invoked

    if (which[j] == COMPUTE) {
      Compute *compute = modify->compute[m];

      if (argindex[j] == 0) {
        if (!(compute->invoked_flag & INVOKED_VECTOR)) {
          compute->compute_vector();
          compute->invoked_flag |= INVOKED_VECTOR;
        }
        double *cvector = compute->vector;
        for (i = 0; i < nrows; i++)
          column[i] = cvector[i];

      } else {
        if (!(compute->invoked_flag & INVOKED_ARRAY)) {
          compute->compute_array();
          compute->invoked_flag |= INVOKED_ARRAY;
        }
        double **carray = compute->array;
        int icol = argindex[j]-1;
        for (i = 0; i < nrows; i++)
          column[i] = carray[i][icol];
      }
    }

    // add columns of values to array or just set directly if offcol is set

    if (offcol[j]) {
      for (i = 0; i < nrows; i++)
        array[i][j] = column[i];
    } else {
      for (i = 0; i < nrows; i++)
        array[i][j] += column[i];
    }
  }

  // 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);

  // average the final result for the Nfreq timestep

  double repeat = nrepeat;
  for (i = 0; i < nrows; i++)
    for (j = 0; j < nvalues; j++)
      if (offcol[j] == 0) array[i][j] /= repeat;

  // 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) {
    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) array_total[i][j] = array[i][j];
    norm = 1;

  } else if (ave == RUNNING) {
    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) array_total[i][j] += array[i][j];
    norm++;

  } else if (ave == WINDOW) {
    for (i = 0; i < nrows; i++)
      for (j = 0; j < nvalues; j++) {
        array_total[i][j] += array[i][j];
        if (window_limit) array_total[i][j] -= array_list[iwindow][i][j];
        array_list[iwindow][i][j] = array[i][j];
      }

    iwindow++;
    if (iwindow == nwindow) {
//.........这里部分代码省略.........