lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4494 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2010-08-17 21:58:16 +00:00
parent 315c8187dd
commit 4529d7187c
10 changed files with 996 additions and 254 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
src/compute_reduce.cpp
View File

@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
enum{SUM,MINN,MAXX,AVE};
enum{X,V,F,COMPUTE,FIX,VARIABLE};
enum{GLOBAL,PERATOM,LOCAL};
enum{PERATOM,LOCAL};
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
@ -181,20 +181,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all("Compute ID for compute reduce does not exist");
if (modify->compute[icompute]->vector_flag ||
modify->compute[icompute]->array_flag) {
flavor[i] = GLOBAL;
if (argindex[i] == 0 &&
modify->compute[icompute]->vector_flag == 0)
error->all("Compute reduce compute does not "
"calculate a global vector");
if (argindex[i] && modify->compute[icompute]->array_flag == 0)
error->all("Compute reduce compute does not "
"calculate a global array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_array_cols)
error->all("Compute reduce compute array is accessed out-of-range");
} else if (modify->compute[icompute]->peratom_flag) {
if (modify->compute[icompute]->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
@ -218,26 +205,13 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_local_cols)
error->all("Compute reduce compute array is accessed out-of-range");
}
} else error->all("Compute reduce compute calculates global values");
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all("Fix ID for compute reduce does not exist");
if (modify->fix[ifix]->vector_flag ||
modify->fix[ifix]->array_flag) {
flavor[i] = GLOBAL;
if (argindex[i] == 0 &&
modify->fix[ifix]->vector_flag == 0)
error->all("Compute reduce fix does not "
"calculate a global vector");
if (argindex[i] && modify->fix[ifix]->array_flag == 0)
error->all("Compute reduce fix does not "
"calculate a global array");
if (argindex[i] &&
argindex[i] > modify->fix[ifix]->size_array_cols)
error->all("Compute reduce fix array is accessed out-of-range");
} else if (modify->fix[ifix]->peratom_flag) {
if (modify->fix[ifix]->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 &&
modify->fix[ifix]->size_peratom_cols != 0)
@ -261,7 +235,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
if (argindex[i] &&
argindex[i] > modify->fix[ifix]->size_local_cols)
error->all("Compute reduce fix array is accessed out-of-range");
}
} else error->all("Compute reduce fix calculates global values");
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
@ -354,25 +328,13 @@ double ComputeReduce::compute_scalar()
double one = compute_one(0,-1);
if (mode == SUM) {
if (flavor[0] == GLOBAL)
scalar = one;
else
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
} else if (mode == MINN) {
if (flavor[0] == GLOBAL)
scalar = one;
else
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
} else if (mode == MAXX) {
if (flavor[0] == GLOBAL)
scalar = one;
else
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
} else if (mode == AVE) {
if (flavor[0] == GLOBAL)
scalar = one;
else
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
scalar /= count(0);
}
@ -392,34 +354,22 @@ void ComputeReduce::compute_vector()
}
if (mode == SUM) {
for (int m = 0; m < nvalues; m++) {
if (flavor[m] == GLOBAL)
vector[m] = onevec[m];
else
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
}
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
} else if (mode == MINN) {
if (!replace) {
for (int m = 0; m < nvalues; m++) {
if (flavor[m] == GLOBAL)
vector[m] = onevec[m];
else
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MIN,world);
}
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MIN,world);
} else {
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
if (flavor[m] == GLOBAL) {
vector[m] = onevec[m];
owner[m] = me;
} else {
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
@ -431,25 +381,17 @@ void ComputeReduce::compute_vector()
} else if (mode == MAXX) {
if (!replace) {
for (int m = 0; m < nvalues; m++) {
if (flavor[m] == GLOBAL)
vector[m] = onevec[m];
else
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MAX,world);
}
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MAX,world);
} else {
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
if (flavor[m] == GLOBAL) {
vector[m] = onevec[m];
owner[m] = me;
} else {
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
@ -461,10 +403,7 @@ void ComputeReduce::compute_vector()
} else if (mode == AVE) {
for (int m = 0; m < nvalues; m++) {
if (flavor[m] == GLOBAL)
vector[m] = onevec[m];
else
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
vector[m] /= count(m);
}
}
@ -524,33 +463,7 @@ double ComputeReduce::compute_one(int m, int flag)
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[vidx];
if (flavor[m] == GLOBAL) {
if (aidx == 0) {
if (!(compute->invoked_flag & INVOKED_VECTOR)) {
compute->compute_vector();
compute->invoked_flag |= INVOKED_VECTOR;
}
double *comp_vec = compute->vector;
int n = compute->size_vector;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,comp_vec[i],i);
else one = comp_vec[flag];
} else {
if (!(compute->invoked_flag & INVOKED_ARRAY)) {
compute->compute_array();
compute->invoked_flag |= INVOKED_ARRAY;
}
double **carray = compute->array;
int n = compute->size_array_rows;
int aidx_1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,carray[i][aidx_1],i);
else one = carray[flag][aidx_1];
}
} else if (flavor[m] == PERATOM) {
if (flavor[m] == PERATOM) {
if (!(compute->invoked_flag & INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= INVOKED_PERATOM;
@ -604,22 +517,7 @@ double ComputeReduce::compute_one(int m, int flag)
error->all("Fix used in compute reduce not computed at compatible time");
Fix *fix = modify->fix[vidx];
if (flavor[m] == GLOBAL) {
if (aidx == 0) {
int n = fix->size_vector;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,fix->compute_vector(i),i);
else one = fix->compute_vector(flag);
} else {
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < nlocal; i++)
combine(one,fix->compute_array(i,aidxm1),i);
else one = fix->compute_array(flag,aidxm1);
}
} else if (flavor[m] == PERATOM) {
if (flavor[m] == PERATOM) {
if (aidx == 0) {
double *fix_vector = fix->vector_atom;
int n = nlocal;
@ -686,10 +584,7 @@ double ComputeReduce::count(int m)
return group->count(igroup);
else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[vidx];
if (flavor[m] == GLOBAL) {
if (aidx == 0) return(1.0*compute->size_vector);
else return(1.0*compute->size_array_rows);
} else if (flavor[m] == PERATOM) {
if (flavor[m] == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
double ncount = compute->size_local_rows;
@ -699,10 +594,7 @@ double ComputeReduce::count(int m)
}
} else if (which[m] == FIX) {
Fix *fix = modify->fix[vidx];
if (flavor[m] == GLOBAL) {
if (aidx == 0) return(1.0*fix->size_vector);
else return(1.0*fix->size_array_rows);
} else if (flavor[m] == PERATOM) {
if (flavor[m] == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
double ncount = fix->size_local_rows;

6
src/fix_ave_atom.cpp
View File

@ -118,9 +118,9 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0) error->all("Illegal fix ave/atom command");
if (peratom_freq < nevery || peratom_freq % nevery ||
(nrepeat-1)*nevery >= peratom_freq)
if (nevery <= 0 || nrepeat <= 0 || peratom_freq <= 0)
error->all("Illegal fix ave/atom command");
if (peratom_freq % nevery || (nrepeat-1)*nevery >= peratom_freq)
error->all("Illegal fix ave/atom command");
for (int i = 0; i < nvalues; i++) {

574
src/fix_ave_correlate.cpp Normal file
View File

@ -0,0 +1,574 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Benoit Leblanc, (Materials Design)
Reese Jones (Sandia)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "fix_ave_correlate.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{COMPUTE,FIX,VARIABLE};
enum{ONE,RUNNING};
enum{AUTO,UPPER,LOWER,AUTOUPPER,AUTOLOWER,FULL};
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
/* ---------------------------------------------------------------------- */
FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
Fix (lmp, narg, arg)
{
if (narg < 7) error->all ("Illegal fix ave/correlate command");
MPI_Comm_rank(world,&me);
nevery = atoi(arg[3]);
nrepeat = atoi(arg[4]);
nfreq = atoi(arg[5]);
global_freq = nfreq;
time_depend = 1;
// parse values until one isn't recognized
which = new int[narg-6];
argindex = new int[narg-6];
ids = new char*[narg-6];
value2index = new int[narg-6];
nvalues = 0;
int iarg = 6;
while (iarg < narg) {
if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"v_",2) == 0) {
if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all("Illegal fix ave/correlate command");
argindex[nvalues] = atoi(ptr+1);
*ptr = '\0';
} else argindex[nvalues] = 0;
n = strlen(suffix) + 1;
ids[nvalues] = new char[n];
strcpy(ids[nvalues],suffix);
delete [] suffix;
nvalues++;
iarg++;
} else break;
}
// optional args
type = AUTO;
ave = ONE;
startstep = 0;
prefactor = 1.0;
fp = NULL;
char *title1 = NULL;
char *title2 = NULL;
char *title3 = NULL;
while (iarg < narg) {
if (strcmp(arg[iarg],"type") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
if (strcmp(arg[iarg+1],"auto") == 0) type = AUTO;
else if (strcmp(arg[iarg+1],"upper") == 0) type = UPPER;
else if (strcmp(arg[iarg+1],"lower") == 0) type = LOWER;
else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER;
else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER;
else if (strcmp(arg[iarg+1],"full") == 0) type = FULL;
else error->all("Illegal fix ave/correlate command");
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
else error->all("Illegal fix ave/correlate command");
iarg += 2;
} else if (strcmp(arg[iarg],"start") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
startstep = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"prefactor") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
prefactor = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
if (me == 0) {
fp = fopen(arg[iarg+1],"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix ave/correlate file %s",arg[iarg+1]);
error->one(str);
}
}
iarg += 2;
} else if (strcmp(arg[iarg],"title1") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
delete [] title1;
int n = strlen(arg[iarg+1]) + 1;
title1 = new char[n];
strcpy(title1,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title2") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
delete [] title2;
int n = strlen(arg[iarg+1]) + 1;
title2 = new char[n];
strcpy(title2,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title3") == 0) {
if (iarg+2 > narg) error->all("Illegal fix ave/correlate command");
delete [] title3;
int n = strlen(arg[iarg+1]) + 1;
title3 = new char[n];
strcpy(title3,arg[iarg+1]);
iarg += 2;
} else error->all("Illegal fix ave/correlate command");
}
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all("Illegal fix ave/correlate command");
if (nfreq % nevery || (nrepeat-1)*nevery > nfreq)
error->all("Illegal fix ave/correlate command");
for (int i = 0; i < nvalues; i++) {
if (which[i] == COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all ("Compute ID for fix ave/correlate does not exist");
if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0)
error->all ("Fix ave/correlate compute does not calculate a scalar");
if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
error->all ("Fix ave/correlate compute does not calculate a vector");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
error->all ("Fix ave/correlate compute vector "
"is accessed out-of-range");
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all ("Fix ID for fix ave/correlate does not exist");
if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0)
error->all ("Fix ave/correlate fix does not calculate a scalar");
if (argindex[i] && modify->fix[ifix]->vector_flag == 0)
error->all ("Fix ave/correlate fix does not calculate a vector");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
error->all ("Fix ave/correlate fix vector is accessed out-of-range");
if (nevery % modify->fix[ifix]->global_freq)
error->all("Fix for fix ave/correlate "
"not computed at compatible time");
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all ("Variable name for fix ave/correlate does not exist");
if (input->variable->equalstyle(ivariable) == 0)
error->all ("Fix ave/correlate variable is not equal-style variable");
}
}
// npair = # of correlation pairs to calculate
if (type == AUTO) npair = nvalues;
if (type == UPPER || type == LOWER) npair = nvalues*(nvalues-1)/2;
if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues*(nvalues+1)/2;
if (type == FULL) npair = nvalues*nvalues;
// print file comment lines
if (fp && me == 0) {
if (title1) fprintf(fp,"%s\n",title1);
else fprintf(fp,"# Time-correlated data for fix %s\n",id);
if (title2) fprintf(fp,"%s\n",title2);
else fprintf(fp,"# Timestep Number-of-time-windows\n");
if (title3) fprintf(fp,"%s\n",title3);
else {
fprintf(fp,"# Index TimeDelta Ncount");
if (type == AUTO)
for (int i = 0; i < nvalues; i++)
fprintf(fp," %s*%s",arg[6+i],arg[6+i]);
else if (type == UPPER)
for (int i = 0; i < nvalues; i++)
for (int j = i+1; j < nvalues; j++)
fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
else if (type == LOWER)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < i-1; j++)
fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
else if (type == AUTOUPPER)
for (int i = 0; i < nvalues; i++)
for (int j = i; j < nvalues; j++)
fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
else if (type == AUTOLOWER)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < i; j++)
fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
else if (type == FULL)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < nvalues; j++)
fprintf(fp," %s*%s",arg[6+i],arg[6+j]);
fprintf(fp,"\n");
}
}
delete [] title1;
delete [] title2;
delete [] title3;
// allocate and initialize memory for averaging
// set count and corr to zero since they accumulate
// also set save versions to zero in case accessed via compute_array()
values = memory->create_2d_double_array(nrepeat,nvalues,
"ave/correlate:values");
count = (int *) memory->smalloc(nrepeat*sizeof(int),
"ave/correlate:count");
save_count = (int *) memory->smalloc(nrepeat*sizeof(int),
"ave/correlate:save_count");
corr = memory->create_2d_double_array(nrepeat,npair,
"ave/correlate:corr");
save_corr = memory->create_2d_double_array(nrepeat,npair,
"ave/correlate:save_corr");
int i,j;
for (i = 0; i < nrepeat; i++) {
save_count[i] = count[i] = 0;
for (j = 0; j < npair; j++)
save_corr[i][j] = corr[i][j] = 0.0;
}
lastindex = -1;
firstindex = 0;
nsample = 0;
// this fix produces a global array
array_flag = 1;
size_array_rows = nrepeat;
size_array_cols = npair+2;
extarray = 0;
// nvalid = next step on which end_of_step does something
// this step if multiple of nevery, else next multiple
// startstep is lower bound
nvalid = update->ntimestep;
if (startstep > nvalid) nvalid = startstep;
if (nvalid % nevery) nvalid = (nvalid/nevery)*nevery + nevery;
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can set timestep for ones actually invoked
modify->addstep_compute_all(nvalid);
}
/* ---------------------------------------------------------------------- */
FixAveCorrelate::~FixAveCorrelate()
{
delete [] which;
delete [] argindex;
delete [] value2index;
for (int i = 0; i < nvalues; i++) delete [] ids[i];
delete [] ids;
memory->destroy_2d_double_array(values);
memory->sfree(count);
memory->sfree(save_count);
memory->destroy_2d_double_array(corr);
memory->destroy_2d_double_array(save_corr);
if (fp && me == 0) fclose(fp);
}
/* ---------------------------------------------------------------------- */
int FixAveCorrelate::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixAveCorrelate::init()
{
// set current indices for all computes,fixes,variables
for (int i = 0; i < nvalues; i++) {
if (which[i] == COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all("Compute ID for fix ave/correlate does not exist");
value2index[i] = icompute;
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all("Fix ID for fix ave/correlate does not exist");
value2index[i] = ifix;
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all("Variable name for fix ave/correlate does not exist");
value2index[i] = ivariable;
}
}
}
/* ----------------------------------------------------------------------
only does something if nvalid = current timestep
------------------------------------------------------------------------- */
void FixAveCorrelate::setup(int vflag)
{
end_of_step();
}
/* ---------------------------------------------------------------------- */
void FixAveCorrelate::end_of_step()
{
int i,j,k,m;
double scalar;
// skip if not step which requires doing something
int ntimestep = update->ntimestep;
if (ntimestep != nvalid) return;
// 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 to file
if (fp && me == 0) {
fprintf(fp,"%d %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++)
fprintf(fp," %g",prefactor*corr[i][j]/count[i]);
else
for (j = 0; j < npair; j++)
fprintf(fp," 0.0");
fprintf(fp,"\n");
}
fflush(fp);
}
// zero accumulation if requested
// recalculate Cij(0)
if (ave == ONE) {
for (i = 0; i < nrepeat; i++) {
count[i] = 0;
for (j = 0; j < npair; j++)
corr[i][j] = 0.0;
}
nsample = 1;
accumulate();
}
}
/* ----------------------------------------------------------------------
accumulate correlation data using more recently added values
------------------------------------------------------------------------- */
void FixAveCorrelate::accumulate()
{
int i,j,k,m,n,ipair;
for (k = 0; k < nsample; k++) count[k]++;
if (type == AUTO) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++) {
corr[k][ipair++] += values[m][i]*values[n][i];
}
m--;
if (m < 0) m = nrepeat-1;
}
} else if (type == UPPER) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = i+1; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
} else if (type == LOWER) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j < i-1; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
} else if (type == AUTOUPPER) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = i; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
} else if (type == AUTOLOWER) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j < i; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
} else if (type == FULL) {
m = n = lastindex;
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
}
}
/* ----------------------------------------------------------------------
return I,J array value
------------------------------------------------------------------------- */
double FixAveCorrelate::compute_array(int i, int j)
{
if (j == 0) return 1.0*i*nevery;
else if (j == 1) return 1.0*save_count[i];
else if (save_count[i]) return save_corr[i][j-2];
return 0.0;
}

66
src/fix_ave_correlate.h Normal file
View File

@ -0,0 +1,66 @@
/* ----------------------------------------------------------------------
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/correlate,FixAveCorrelate)
#else
#ifndef LMP_FIX_AVE_CORRELATE_H
#define LMP_FIX_AVE_CORRELATE_H
#include "stdio.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixAveCorrelate : public Fix {
public:
FixAveCorrelate(class LAMMPS *, int, char **);
~FixAveCorrelate();
int setmask();
void init();
void setup(int);
void end_of_step();
double compute_array(int,int);
private:
int me,nvalues;
int nrepeat,nfreq,nvalid;
int *which,*argindex,*value2index;
char **ids;
FILE *fp;
int type,ave,startstep;
double prefactor;
char *title1,*title2,*title3;
int firstindex; // index in values ring of earliest time sample
int lastindex; // index in values ring of latest time sample
int nsample; // number of time samples in values ring
int npair; // number of correlation pairs to calculate
int *count;
double **values,**corr;
int *save_count; // saved values at Nfreq for output via compute_array()
double **save_corr;
void accumulate();
};
}
#endif
#endif

6
src/fix_ave_histo.cpp
View File

@ -238,10 +238,10 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
// kind = inputs are all global, or all per-atom, or all local
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0) error->all("Illegal fix ave/histo command");
if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all("Illegal fix ave/histo command");
if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
error->all("Illegal fix ave/histo command");
if (lo >= hi) error->all("Illegal fix ave/histo command");
if (nbins <= 0) error->all("Illegal fix ave/histo command");

6
src/fix_ave_spatial.cpp
View File

@ -225,10 +225,10 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
// setup and error check
if (nevery <= 0) error->all("Illegal fix ave/spatial command");
if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all("Illegal fix ave/spatial command");
if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
error->all("Illegal fix ave/spatial command");
if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
for (int i = 0; i < nvalues; i++) {

7
src/fix_ave_time.cpp
View File

@ -165,8 +165,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0) error->all("Illegal fix ave/time command");
if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all("Illegal fix ave/time command");
if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
error->all("Illegal fix ave/time command");
for (int i = 0; i < nvalues; i++) {
@ -284,7 +285,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
delete [] title2;
delete [] title3;
// allocate and initialize memory for averaging
// allocate memory for averaging
vector = vector_total = NULL;
vector_list = NULL;

2
src/math_extra.h
View File

@ -83,7 +83,7 @@ namespace MathExtra {
inline void multiply_shape_shape(const double *one, const double *two,
double *ans);
};
}
/* ----------------------------------------------------------------------
normalize a vector

406
src/variable.cpp
View File

@ -35,6 +35,9 @@ using namespace LAMMPS_NS;
#define VARDELTA 4
#define MAXLEVEL 4
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
@ -42,14 +45,17 @@ enum{ARG,OP};
enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
EQ,NE,LT,LE,GT,GE,AND,OR,
SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,
CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,
CEIL,FLOOR,ROUND,
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
enum{SUM,XMIN,XMAX,AVE,TRAP};
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
#define INVOKED_PERATOM 8
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
@ -758,10 +764,10 @@ double Variable::evaluate(char *str, Tree **tree)
} else if (nbracket == 2 && compute->array_flag) {
if (index1 > compute->size_array_rows)
error->all("Variable formula compute vector "
error->all("Variable formula compute array "
"is accessed out-of-range");
if (index2 > compute->size_array_cols)
error->all("Variable formula compute vector "
error->all("Variable formula compute array "
"is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_array != update->ntimestep)
@ -804,7 +810,7 @@ double Variable::evaluate(char *str, Tree **tree)
compute->size_peratom_cols > 0) {
if (index2 > compute->size_peratom_cols)
error->all("Variable formula compute vector "
error->all("Variable formula compute array "
"is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
@ -850,7 +856,7 @@ double Variable::evaluate(char *str, Tree **tree)
if (tree == NULL)
error->all("Per-atom compute in equal-style variable formula");
if (index1 > compute->size_peratom_cols)
error->all("Variable formula compute vector "
error->all("Variable formula compute array "
"is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
@ -946,9 +952,9 @@ double Variable::evaluate(char *str, Tree **tree)
} else if (nbracket == 2 && fix->array_flag) {
if (index1 > fix->size_array_rows)
error->all("Variable formula fix vector is accessed out-of-range");
error->all("Variable formula fix array is accessed out-of-range");
if (index2 > fix->size_array_cols)
error->all("Variable formula fix vector is accessed out-of-range");
error->all("Variable formula fix array is accessed out-of-range");
if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
error->all("Fix in variable not computed at compatible time");
@ -979,7 +985,7 @@ double Variable::evaluate(char *str, Tree **tree)
fix->size_peratom_cols > 0) {
if (index2 > fix->size_peratom_cols)
error->all("Variable formula fix vector is accessed out-of-range");
error->all("Variable formula fix array is accessed out-of-range");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
@ -1014,7 +1020,7 @@ double Variable::evaluate(char *str, Tree **tree)
if (tree == NULL)
error->all("Per-atom fix in equal-style variable formula");
if (index1 > fix->size_peratom_cols)
error->all("Variable formula fix vector is accessed out-of-range");
error->all("Variable formula fix array is accessed out-of-range");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
@ -1098,13 +1104,13 @@ double Variable::evaluate(char *str, Tree **tree)
delete [] id;
// ----------------
// math/group function or atom value/vector or thermo keyword
// math/group/special function or atom value/vector or thermo keyword
// ----------------
} else {
// ----------------
// math or group function
// math or group or special function
// ----------------
if (str[i] == '(') {
@ -1116,7 +1122,10 @@ double Variable::evaluate(char *str, Tree **tree)
treestack,ntreestack,argstack,nargstack));
else if (group_function(word,contents,tree,
treestack,ntreestack,argstack,nargstack));
else error->all("Invalid math or group function in variable formula");
else if (special_function(word,contents,tree,
treestack,ntreestack,argstack,nargstack));
else error->all("Invalid math/group/special function "
"in variable formula");
delete [] contents;
// ----------------
@ -1311,7 +1320,7 @@ double Variable::evaluate(char *str, Tree **tree)
process an evaulation tree
customize by adding a math function:
sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan()
ceil(),floor(),round(),ramp(),stagger(),logfreq()
ceil(),floor(),round()
---------------------------------------------------------------------- */
double Variable::eval_tree(Tree *tree, int i)
@ -1424,43 +1433,6 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == ROUND)
return MYROUND(eval_tree(tree->left,i));
if (tree->type == RAMP) {
if (update->whichflag == 0)
error->all("Cannot use ramp in variable formula between runs");
double arg1 = eval_tree(tree->left,i);
double arg2 = eval_tree(tree->right,i);
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
return arg1 + delta*(arg2-arg1);
}
if (tree->type == STAGGER) {
int arg1 = static_cast<int> (eval_tree(tree->left,i));
int arg2 = static_cast<int> (eval_tree(tree->right,i));
if (arg1 <= 0 || arg2 <= 0 || arg1 <= arg2)
error->all("Invalid math function in variable formula");
int lower = update->ntimestep/arg1 * arg1;
int delta = update->ntimestep - lower;
if (delta < arg2) return 1.0*(lower+arg2);
else return 1.0*(lower+arg1);
}
if (tree->type == LOGFREQ) {
int arg1 = static_cast<int> (eval_tree(tree->left,i));
int arg2 = static_cast<int> (eval_tree(tree->middle,i));
int arg3 = static_cast<int> (eval_tree(tree->right,i));
if (arg1 <= 0 || arg2 <= 0 || arg3 <= 0 || arg2 >= arg3)
error->all("Invalid math function in variable formula");
if (update->ntimestep < arg1) return 1.0*arg1;
else {
int lower = arg1;
while (update->ntimestep >= arg3*lower) lower *= arg3;
int multiple = update->ntimestep/lower;
if (multiple < arg2) return 1.0*(multiple+1)*lower;
else return 1.0*(lower*arg3);
}
}
return 0.0;
}
@ -1518,19 +1490,19 @@ int Variable::int_between_brackets(char *&ptr)
while (*ptr && *ptr != ']') {
if (!isdigit(*ptr))
error->all("Non digit character between brackets in input command");
error->all("Non digit character between brackets in variable");
ptr++;
}
if (*ptr != ']') error->all("Mismatched brackets in input command");
if (ptr == start) error->all("Empty brackets in input command");
if (*ptr != ']') error->all("Mismatched brackets in variable");
if (ptr == start) error->all("Empty brackets in variable");
*ptr = '\0';
int index = atoi(start);
*ptr = ']';
if (index == 0)
error->all("Index between input command brackets must be positive");
error->all("Index between variable brackets must be positive");
return index;
}
@ -1542,7 +1514,7 @@ int Variable::int_between_brackets(char *&ptr)
return 0 if not a match, 1 if successfully processed
customize by adding a math function in 2 places:
sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan()
ceil(),floor(),round(),ramp(),stagger(),logfreq()
ceil(),floor(),round()
------------------------------------------------------------------------- */
int Variable::math_function(char *word, char *contents, Tree **tree,
@ -1557,8 +1529,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
strcmp(word,"tan") && strcmp(word,"asin") &&
strcmp(word,"acos") && strcmp(word,"atan") &&
strcmp(word,"ceil") && strcmp(word,"floor") &&
strcmp(word,"round") && strcmp(word,"ramp") &&
strcmp(word,"stagger") && strcmp(word,"logfreq"))
strcmp(word,"round"))
return 0;
// parse contents for arg1,arg2,arg3 separated by commas
@ -1709,50 +1680,6 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
if (narg != 1) error->all("Invalid math function in variable formula");
if (tree) newtree->type = ROUND;
else argstack[nargstack++] = MYROUND(value1);
} else if (strcmp(word,"ramp") == 0) {
if (narg != 2) error->all("Invalid math function in variable formula");
if (update->whichflag == 0)
error->all("Cannot use ramp in variable formula between runs");
if (tree) newtree->type = RAMP;
else {
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
argstack[nargstack++] = value1 + delta*(value2-value1);
}
} else if (strcmp(word,"stagger") == 0) {
if (narg != 2) error->all("Invalid math function in variable formula");
if (tree) newtree->type = STAGGER;
else {
int ivalue1 = static_cast<int> (value1);
int ivalue2 = static_cast<int> (value2);
if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
error->all("Invalid math function in variable formula");
int lower = update->ntimestep/ivalue1 * ivalue1;
int delta = update->ntimestep - lower;
if (delta < value2) argstack[nargstack++] = lower+ivalue2;
else argstack[nargstack++] = lower+ivalue1;
}
} else if (strcmp(word,"logfreq") == 0) {
if (narg != 3) error->all("Invalid math function in variable formula");
if (tree) newtree->type = LOGFREQ;
else {
int ivalue1 = static_cast<int> (value1);
int ivalue2 = static_cast<int> (value2);
int ivalue3 = static_cast<int> (value3);
if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
error->all("Invalid math function in variable formula");
if (update->ntimestep < ivalue1) argstack[nargstack++] = ivalue1;
else {
int lower = ivalue1;
while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
int multiple = update->ntimestep/lower;
if (multiple < ivalue2) argstack[nargstack++] = (multiple+1)*lower;
else argstack[nargstack++] = lower*ivalue3;
}
}
}
delete [] arg1;
@ -2007,6 +1934,236 @@ int Variable::region_function(char *id)
return iregion;
}
/* ----------------------------------------------------------------------
process a special function in formula
push result onto tree or arg stack
word = special function
contents = str between parentheses with one,two,three args
return 0 if not a match, 1 if successfully processed
customize by adding a special function in 2 places
ramp(x,y),stagger(x,y),logfreq(x,y,z),sum(x),min(x),max(x),
ave(x),trap(x)
------------------------------------------------------------------------- */
int Variable::special_function(char *word, char *contents, Tree **tree,
Tree **treestack, int &ntreestack,
double *argstack, int &nargstack)
{
// word not a match to any special function
if (strcmp(word,"ramp") && strcmp(word,"stagger") &&
strcmp(word,"logfreq") && strcmp(word,"sum") &&
strcmp(word,"min") && strcmp(word,"max") &&
strcmp(word,"ave") && strcmp(word,"trap"))
return 0;
// parse contents for arg1,arg2,arg3 separated by commas
// ptr1,ptr2 = location of 1st and 2nd comma, NULL if none
char *arg1,*arg2,*arg3;
char *ptr1,*ptr2;
ptr1 = strchr(contents,',');
if (ptr1) {
*ptr1 = '\0';
ptr2 = strchr(ptr1+1,',');
if (ptr2) *ptr2 = '\0';
} else ptr2 = NULL;
int n = strlen(contents) + 1;
arg1 = new char[n];
strcpy(arg1,contents);
int narg = 1;
if (ptr1) {
n = strlen(ptr1+1) + 1;
arg2 = new char[n];
strcpy(arg2,ptr1+1);
narg = 2;
} else arg2 = NULL;
if (ptr2) {
n = strlen(ptr2+1) + 1;
arg3 = new char[n];
strcpy(arg3,ptr2+1);
narg = 3;
} else arg3 = NULL;
// match word to special function
double value;
if (strcmp(word,"ramp") == 0) {
if (narg != 2) error->all("Invalid special function in variable formula");
if (update->whichflag == 0)
error->all("Cannot use ramp in variable formula between runs");
double value1 = numeric(arg1);
double value2 = numeric(arg2);
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
value = value1 + delta*(value2-value1);
} else if (strcmp(word,"stagger") == 0) {
if (narg != 2) error->all("Invalid special function in variable formula");
int ivalue1 = inumeric(arg1);
int ivalue2 = inumeric(arg2);
if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
error->all("Invalid special function in variable formula");
int lower = update->ntimestep/ivalue1 * ivalue1;
int delta = update->ntimestep - lower;
if (delta < ivalue2) value = lower+ivalue2;
else value = lower+ivalue1;
} else if (strcmp(word,"logfreq") == 0) {
if (narg != 3) error->all("Invalid special function in variable formula");
int ivalue1 = inumeric(arg1);
int ivalue2 = inumeric(arg2);
int ivalue3 = inumeric(arg3);
if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
error->all("Invalid special function in variable formula");
if (update->ntimestep < ivalue1) value = ivalue1;
else {
int lower = ivalue1;
while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
int multiple = update->ntimestep/lower;
if (multiple < ivalue2) value = (multiple+1)*lower;
else value = lower*ivalue3;
}
} else {
if (narg != 1) error->all("Invalid special function in variable formula");
int method;
if (strcmp(word,"sum") == 0) method = SUM;
else if (strcmp(word,"min") == 0) method = XMIN;
else if (strcmp(word,"max") == 0) method = XMAX;
else if (strcmp(word,"ave") == 0) method = AVE;
else if (strcmp(word,"trap") == 0) method = TRAP;
Compute *compute = NULL;
Fix *fix = NULL;
int index,nvec,nstride;
if (strstr(arg1,"c_") == arg1) {
ptr1 = strchr(arg1,'[');
if (ptr1) {
ptr2 = ptr1;
index = int_between_brackets(ptr2);
*ptr1 = '\0';
} else index = 0;
int icompute = modify->find_compute(&arg1[2]);
if (icompute < 0) error->all("Invalid compute ID in variable formula");
compute = modify->compute[icompute];
if (index == 0 && compute->vector_flag) {
if (update->whichflag == 0) {
if (compute->invoked_vector != update->ntimestep)
error->all("Compute used in variable between runs is not current");
} else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
compute->compute_vector();
compute->invoked_flag |= INVOKED_VECTOR;
}
nvec = compute->size_vector;
nstride = 1;
} else if (index && compute->array_flag) {
if (index > compute->size_array_cols)
error->all("Variable formula compute array "
"is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_array != update->ntimestep)
error->all("Compute used in variable between runs is not current");
} else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
compute->compute_array();
compute->invoked_flag |= INVOKED_ARRAY;
}
nvec = compute->size_array_rows;
nstride = compute->size_array_cols;
} else error->all("Mismatched compute in variable formula");
} else if (strstr(arg1,"f_") == arg1) {
ptr1 = strchr(arg1,'[');
if (ptr1) {
ptr2 = ptr1;
index = int_between_brackets(ptr2);
*ptr1 = '\0';
} else index = 0;
int ifix = modify->find_fix(&arg1[2]);
if (ifix < 0) error->all("Invalid fix ID in variable formula");
fix = modify->fix[ifix];
if (index == 0 && fix->vector_flag) {
if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
error->all("Fix in variable not computed at compatible time");
nvec = fix->size_vector;
nstride = 1;
} else if (index && fix->array_flag) {
if (index > fix->size_array_cols)
error->all("Variable formula fix array is accessed out-of-range");
if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
error->all("Fix in variable not computed at compatible time");
nvec = fix->size_array_rows;
nstride = fix->size_array_cols;
} else error->all("Mismatched fix in variable formula");
} else error->all("Invalid special function in variable formula");
value = 0.0;
if (method == XMIN) value = BIG;
if (method == XMAX) value = -BIG;
if (compute) {
double *vec;
if (index) vec = &compute->array[0][index];
else vec = compute->vector;
int j = 0;
for (int i = 0; i < nvec; i++) {
if (method == SUM) value += vec[j];
else if (method == XMIN) value = MIN(value,vec[j]);
else if (method == XMAX) value = MAX(value,vec[j]);
else if (method == AVE) value += vec[j];
else if (method == TRAP) {
if (i > 0 && i < nvec-1) value += vec[j];
else value += 0.5*vec[j];
}
j += nstride;
}
}
if (fix) {
double one;
for (int i = 0; i < nvec; i++) {
if (index) one = fix->compute_array(i,index-1);
else one = fix->compute_vector(i);
if (method == SUM) value += one;
else if (method == XMIN) value = MIN(value,one);
else if (method == XMAX) value = MAX(value,one);
else if (method == AVE) value += one;
else if (method == TRAP) {
if (i > 1 && i < nvec) value += one;
else value += 0.5*one;
}
}
}
if (method == AVE) value /= nvec;
}
delete [] arg1;
delete [] arg2;
delete [] arg3;
// save value in tree or on argstack
if (tree) {
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value;
newtree->left = newtree->middle = newtree->right = NULL;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value;
return 1;
}
/* ----------------------------------------------------------------------
extract a global value from a per-atom quantity in a formula
flag = 0 -> word is an atom vector
@ -2127,3 +2284,50 @@ void Variable::atom_vector(char *word, Tree **tree,
else if (strcmp(word,"fy") == 0) newtree->array = &atom->f[0][1];
else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
}
/* ----------------------------------------------------------------------
read a floating point value from a string
generate an error if not a legitimate floating point value
------------------------------------------------------------------------- */
double Variable::numeric(char *str)
{
int n = strlen(str);
for (int i = 0; i < n; i++) {
if (isdigit(str[i])) continue;
if (str[i] == '-' || str[i] == '+' || str[i] == '.') continue;
if (str[i] == 'e' || str[i] == 'E') continue;
error->all("Expected floating point parameter in variable definition");
}
return atof(str);
}
/* ----------------------------------------------------------------------
read an integer value from a string
generate an error if not a legitimate integer value
------------------------------------------------------------------------- */
int Variable::inumeric(char *str)
{
int n = strlen(str);
for (int i = 0; i < n; i++) {
if (isdigit(str[i]) || str[i] == '-' || str[i] == '+') continue;
error->all("Expected integer parameter in variable definition");
}
return atoi(str);
}
/* ----------------------------------------------------------------------
debug routine for printing formula tree recursively
------------------------------------------------------------------------- */
void Variable::print_tree(Tree *tree, int level)
{
printf("TREE %d: %d %g\n",level,tree->type,tree->value);
if (tree->left) print_tree(tree->left,level+1);
if (tree->middle) print_tree(tree->middle,level+1);
if (tree->right) print_tree(tree->right,level+1);
return;
}

5
src/variable.h
View File

@ -63,10 +63,15 @@ class Variable : protected Pointers {
int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);
int group_function(char *, char *, Tree **, Tree **, int &, double *, int &);
int region_function(char *);
int special_function(char *, char *, Tree **, Tree **,
int &, double *, int &);
void peratom2global(int, char *, double *, int, int,
Tree **, Tree **, int &, double *, int &);
int is_atom_vector(char *);
void atom_vector(char *, Tree **, Tree **, int &);
double numeric(char *);
int inumeric(char *);
void print_tree(Tree *, int);
};
}