From 1dec0d587bbeda8654fcb9f2e2fca369106da61b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 9 Oct 2007 17:20:33 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@987
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/compute.h | 4 +-
src/fix.cpp | 3 +-
src/fix.h | 8 ++-
src/fix_ave_spatial.cpp | 3 +-
src/fix_ave_time.cpp | 136 ++++++++++++++++++++++++++-------------
src/fix_ave_time.h | 5 +-
src/fix_indent.cpp | 83 +++++++++++++++++-------
src/fix_indent.h | 7 +-
src/fix_nph.cpp | 9 +--
src/fix_nph.h | 2 +-
src/fix_npt.cpp | 9 +--
src/fix_npt.h | 2 +-
src/fix_nvt.cpp | 10 +--
src/fix_nvt.h | 2 +-
src/fix_orient_fcc.cpp | 32 ++++-----
src/fix_orient_fcc.h | 5 +-
src/fix_set_force.cpp | 23 ++++---
src/fix_set_force.h | 5 +-
src/fix_temp_rescale.cpp | 9 ++-
src/fix_temp_rescale.h | 2 +-
src/fix_wall_lj126.cpp | 55 +++++++++++-----
src/fix_wall_lj126.h | 6 +-
src/fix_wall_lj93.cpp | 57 +++++++++++-----
src/fix_wall_lj93.h | 6 +-
src/modify.cpp | 4 +-
src/thermo.cpp | 11 +++-
src/variable.cpp | 37 ++++++++++-
27 files changed, 360 insertions(+), 175 deletions(-)
diff --git a/src/compute.h b/src/compute.h
index 9bcfe3a90e..559bd6ed1c 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -30,9 +30,9 @@ class Compute : protected Pointers {
int scalar_flag; // 0/1 if compute_scalar() function exists
int vector_flag; // 0/1 if compute_vector() function exists
- int peratom_flag; // 0/1 if compute_peratom() function exists
int size_vector; // N = size of global vector
- int size_peratom; // 0 = just scalar_atom, N = size of vector_atom
+ int peratom_flag; // 0/1 if compute_peratom() function exists
+ int size_peratom; // 0 = scalar_atom, N = size of vector_atom
int extensive; // 0/1 if scalar,vector are intensive/extensive values
int tempflag; // 1 if Compute can be used as temperature
diff --git a/src/fix.cpp b/src/fix.cpp
index ae48cf4764..df6d7acb29 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -42,7 +42,8 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
rigid_flag = 0;
virial_flag = 0;
no_change_box = 0;
- peratom_flag = 0;
+
+ scalar_flag = vector_flag = peratom_flag = 0;
comm_forward = comm_reverse = 0;
diff --git a/src/fix.h b/src/fix.h
index f0a2330a41..e77f06aac7 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -35,6 +35,11 @@ class Fix : protected Pointers {
int virial_flag; // 1 if Fix contributes to virial, 0 if not
int no_change_box; // 1 if cannot swap ortho <-> triclinic
+ int scalar_flag; // 0/1 if compute_scalar() function exists
+ int vector_flag; // 0/1 if compute_vector() function exists
+ int size_vector; // N = size of global vector
+ int scalar_vector_freq; // frequency compute s/v data is available at
+
int peratom_flag; // 0/1 if per-atom data is stored
int size_peratom; // 0 = scalar_atom, N = size of vector_atom
double *scalar_atom; // computed per-atom scalar
@@ -93,7 +98,8 @@ class Fix : protected Pointers {
virtual int pack_reverse_comm(int, int, double *) {return 0;}
virtual void unpack_reverse_comm(int, int *, double *) {}
- virtual double thermo(int) {return 0.0;}
+ virtual double compute_scalar() {return 0.0;}
+ virtual double compute_vector(int) {return 0.0;}
virtual int dof(int) {return 0;}
virtual void deform(int) {}
diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp
index f9697c053b..22d4a15a5a 100644
--- a/src/fix_ave_spatial.cpp
+++ b/src/fix_ave_spatial.cpp
@@ -260,10 +260,9 @@ void FixAveSpatial::init()
if (ifix < 0)
error->all("Fix ID for fix ave/spatial does not exist");
fix = modify->fix[ifix];
- if (nevery % modify->fix[ifix]->peratom_freq)
+ if (nevery % fix->peratom_freq)
error->all("Fix ave/spatial and fix not computed at compatible times");
}
-
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp
index ac8f7b1adf..d487ec6f2d 100644
--- a/src/fix_ave_time.cpp
+++ b/src/fix_ave_time.cpp
@@ -26,29 +26,35 @@
using namespace LAMMPS_NS;
+enum{COMPUTE,FIX};
+
/* ---------------------------------------------------------------------- */
FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
- if (narg != 9) error->all("Illegal fix ave/time command");
+ if (narg != 10) error->all("Illegal fix ave/time command");
nevery = atoi(arg[3]);
nrepeat = atoi(arg[4]);
nfreq = atoi(arg[5]);
- int n = strlen(arg[6]) + 1;
- id_compute = new char[n];
- strcpy(id_compute,arg[6]);
+ if (strcmp(arg[6],"compute") == 0) which = COMPUTE;
+ else if (strcmp(arg[6],"fix") == 0) which = FIX;
+ else error->all("Illegal fix ave/time command");
+
+ int n = strlen(arg[7]) + 1;
+ id = new char[n];
+ strcpy(id,arg[7]);
- int flag = atoi(arg[7]);
+ int flag = atoi(arg[8]);
MPI_Comm_rank(world,&me);
if (me == 0) {
- fp = fopen(arg[8],"w");
+ fp = fopen(arg[9],"w");
if (fp == NULL) {
char str[128];
- sprintf(str,"Cannot open fix ave/time file %s",arg[8]);
+ sprintf(str,"Cannot open fix ave/time file %s",arg[9]);
error->one(str);
}
}
@@ -59,31 +65,52 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
error->all("Illegal fix ave/time command");
- int icompute = modify->find_compute(id_compute);
- if (icompute < 0) error->all("Compute ID for fix ave/time does not exist");
+ int icompute,ifix;
+ if (which == COMPUTE) {
+ icompute = modify->find_compute(id);
+ if (icompute < 0) error->all("Compute ID for fix ave/time does not exist");
+ } else {
+ ifix = modify->find_fix(id);
+ if (ifix < 0) error->all("Fix ID for fix ave/time does not exist");
+ }
if (flag < 0 || flag > 2) error->all("Illegal fix ave/time command");
sflag = vflag = 0;
if (flag == 0 || flag == 2) sflag = 1;
if (flag == 1 || flag == 2) vflag = 1;
- if (sflag && modify->compute[icompute]->scalar_flag == 0)
- error->all("Fix ave/time compute does not calculate a scalar");
- if (vflag && modify->compute[icompute]->vector_flag == 0)
- error->all("Fix ave/time compute does not calculate a vector");
+ if (which == COMPUTE) {
+ if (sflag && modify->compute[icompute]->scalar_flag == 0)
+ error->all("Fix ave/time compute does not calculate a scalar");
+ if (vflag && modify->compute[icompute]->vector_flag == 0)
+ error->all("Fix ave/time compute does not calculate a vector");
+ } else {
+ if (sflag && modify->fix[ifix]->scalar_flag == 0)
+ error->all("Fix ave/time fix does not calculate a scalar");
+ if (vflag && modify->fix[ifix]->vector_flag == 0)
+ error->all("Fix ave/time fix does not calculate a vector");
+ }
- if (modify->compute[icompute]->pressflag) pressure_every = nevery;
+ if (which == COMPUTE &&
+ modify->compute[icompute]->pressflag) pressure_every = nevery;
// setup list of computes to call, including pre-computes
- ncompute = 1 + modify->compute[icompute]->npre;
- compute = new Compute*[ncompute];
+ compute = NULL;
+ if (which == COMPUTE) {
+ ncompute = 1 + modify->compute[icompute]->npre;
+ compute = new Compute*[ncompute];
+ } else ncompute = 0;
// print header into file
if (me == 0) {
- fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n",
- id,group->names[modify->compute[icompute]->igroup],id_compute);
+ if (which == COMPUTE)
+ fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n",
+ id,group->names[modify->compute[icompute]->igroup],id);
+ else
+ fprintf(fp,"Time-averaged data for fix %s, group %s, and fix %s\n",
+ id,group->names[modify->fix[ifix]->igroup],id);
if (sflag and !vflag)
fprintf(fp,"TimeStep Value\n");
else if (!sflag and vflag)
@@ -94,7 +121,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
vector = NULL;
if (vflag) {
- size_vector = modify->compute[icompute]->size_vector;
+ if (which == COMPUTE) size_vector = modify->compute[icompute]->size_vector;
+ else size_vector = modify->fix[ifix]->size_vector;
vector = new double[size_vector];
}
@@ -111,7 +139,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
FixAveTime::~FixAveTime()
{
- delete [] id_compute;
+ delete [] id;
if (me == 0) fclose(fp);
delete [] compute;
delete [] vector;
@@ -132,20 +160,34 @@ void FixAveTime::init()
{
// set ptrs to one or more computes called each end-of-step
- int icompute = modify->find_compute(id_compute);
- if (icompute < 0)
- error->all("Compute ID for fix ave/time does not exist");
+ if (which == COMPUTE) {
+ int icompute = modify->find_compute(id);
+ if (icompute < 0)
+ error->all("Compute ID for fix ave/time does not exist");
- ncompute = 0;
- if (modify->compute[icompute]->npre)
- for (int i = 0; i < modify->compute[icompute]->npre; i++) {
- int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]);
- if (ic < 0)
- error->all("Precompute ID for fix ave/time does not exist");
- compute[ncompute++] = modify->compute[ic];
- }
-
- compute[ncompute++] = modify->compute[icompute];
+ ncompute = 0;
+ if (modify->compute[icompute]->npre)
+ for (int i = 0; i < modify->compute[icompute]->npre; i++) {
+ int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]);
+ if (ic < 0)
+ error->all("Precompute ID for fix ave/time does not exist");
+ compute[ncompute++] = modify->compute[ic];
+ }
+
+ compute[ncompute++] = modify->compute[icompute];
+ }
+
+ // set ptr to fix ID
+ // check that fix frequency is acceptable
+
+ if (which == FIX) {
+ int ifix = modify->find_fix(id);
+ if (ifix < 0)
+ error->all("Fix ID for fix ave/time does not exist");
+ fix = modify->fix[ifix];
+ if (nevery % fix->scalar_vector_freq)
+ error->all("Fix ave/time and fix not computed at compatible times");
+ }
}
/* ---------------------------------------------------------------------- */
@@ -166,17 +208,25 @@ void FixAveTime::end_of_step()
for (i = 0; i < size_vector; i++) vector[i] = 0.0;
}
- // accumulate results of compute to local copy
+ // accumulate results of compute or fix to local copy
- if (sflag) {
- double value;
- for (i = 0; i < ncompute; i++) value = compute[i]->compute_scalar();
- scalar += value;
- }
- if (vflag) {
- for (i = 0; i < ncompute; i++) compute[i]->compute_vector();
- double *cvector = compute[ncompute-1]->vector;
- for (i = 0; i < size_vector; i++) vector[i] += cvector[i];
+ if (which == COMPUTE) {
+ if (sflag) {
+ double value;
+ for (i = 0; i < ncompute; i++) value = compute[i]->compute_scalar();
+ scalar += value;
+ }
+ if (vflag) {
+ for (i = 0; i < ncompute; i++) compute[i]->compute_vector();
+ double *cvector = compute[ncompute-1]->vector;
+ for (i = 0; i < size_vector; i++) vector[i] += cvector[i];
+ }
+
+ } else {
+ if (sflag) scalar += fix->compute_scalar();
+ if (vflag)
+ for (i = 0; i < size_vector; i++)
+ vector[i] += fix->compute_vector(i);
}
irepeat++;
diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h
index 01fb897c48..af96d59da7 100644
--- a/src/fix_ave_time.h
+++ b/src/fix_ave_time.h
@@ -29,8 +29,8 @@ class FixAveTime : public Fix {
private:
int me;
- int nrepeat,nfreq,nvalid,irepeat;
- char *id_compute;
+ int nrepeat,nfreq,nvalid,irepeat,which,ifix;
+ char *id;
FILE *fp;
int sflag,vflag;
@@ -38,6 +38,7 @@ class FixAveTime : public Fix {
double scalar,*vector;
int ncompute;
class Compute **compute;
+ class Fix *fix;
};
}
diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp
index c4244393da..9f15aad46d 100644
--- a/src/fix_indent.cpp
+++ b/src/fix_indent.cpp
@@ -37,6 +37,12 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix indent command");
+
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 3;
+ scalar_vector_freq = 1;
+
k = atof(arg[3]);
// set defaults
@@ -122,7 +128,6 @@ void FixIndent::init()
void FixIndent::setup()
{
- eflag_enable = 1;
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(1);
else {
@@ -130,14 +135,12 @@ void FixIndent::setup()
post_force_respa(1,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
- eflag_enable = 0;
}
/* ---------------------------------------------------------------------- */
void FixIndent::min_setup()
{
- eflag_enable = 1;
post_force(1);
}
@@ -145,10 +148,6 @@ void FixIndent::min_setup()
void FixIndent::post_force(int vflag)
{
- bool eflag = false;
- if (eflag_enable) eflag = true;
- else if (output->next_thermo == update->ntimestep) eflag = true;
-
// set current r0
// for minimization, always set to r0_stop
@@ -160,8 +159,10 @@ void FixIndent::post_force(int vflag)
r0 = r0_start + delta * (r0_stop-r0_start);
}
- double eng;
- if (eflag) eng = 0.0;
+ // indenter values, 0 = energy, 1-3 = force components
+
+ indenter_flag = 0;
+ indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
// spherical indenter
@@ -179,7 +180,7 @@ void FixIndent::post_force(int vflag)
int *mask = atom->mask;
int nlocal = atom->nlocal;
- double delx,dely,delz,r,dr,fmag;
+ double delx,dely,delz,r,dr,fmag,fx,fy,fz;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -190,10 +191,16 @@ void FixIndent::post_force(int vflag)
dr = r - r0;
if (dr >= 0.0) continue;
fmag = k*dr*dr;
- f[i][0] += delx*fmag/r;
- f[i][1] += dely*fmag/r;
- f[i][2] += delz*fmag/r;
- if (eflag) eng -= k3 * dr*dr*dr;
+ fx = delx*fmag/r;
+ fy = dely*fmag/r;
+ fz = delz*fmag/r;
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ indenter[0] -= k3 * dr*dr*dr;
+ indenter[1] -= fx;
+ indenter[2] -= fy;
+ indenter[3] -= fz;
}
// cylindrical indenter
@@ -220,7 +227,7 @@ void FixIndent::post_force(int vflag)
int *mask = atom->mask;
int nlocal = atom->nlocal;
- double delx,dely,delz,r,dr,fmag;
+ double delx,dely,delz,r,dr,fmag,fx,fy,fz;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -241,14 +248,18 @@ void FixIndent::post_force(int vflag)
dr = r - r0;
if (dr >= 0.0) continue;
fmag = k*dr*dr;
- f[i][0] += delx*fmag/r;
- f[i][1] += dely*fmag/r;
- f[i][2] += delz*fmag/r;
- if (eflag) eng -= k3 * dr*dr*dr;
+ fx = delx*fmag/r;
+ fy = dely*fmag/r;
+ fz = delz*fmag/r;
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ indenter[0] -= k3 * dr*dr*dr;
+ indenter[1] -= fx;
+ indenter[2] -= fy;
+ indenter[3] -= fz;
}
}
-
- if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
@@ -265,12 +276,34 @@ void FixIndent::min_post_force(int vflag)
post_force(vflag);
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ energy of indenter interaction
+------------------------------------------------------------------------- */
+
+double FixIndent::compute_scalar()
+{
+ // only sum across procs one time
-double FixIndent::thermo(int n)
+ if (indenter_flag == 0) {
+ MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world);
+ indenter_flag = 1;
+ }
+ return indenter_all[0];
+}
+
+/* ----------------------------------------------------------------------
+ components of force on indenter
+------------------------------------------------------------------------- */
+
+double FixIndent::compute_vector(int n)
{
- if (n == 0) return etotal;
- else return 0.0;
+ // only sum across procs one time
+
+ if (indenter_flag == 0) {
+ MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world);
+ indenter_flag = 1;
+ }
+ return indenter_all[n];
}
/* ----------------------------------------------------------------------
diff --git a/src/fix_indent.h b/src/fix_indent.h
index f6e53912ce..f64b5bf62e 100644
--- a/src/fix_indent.h
+++ b/src/fix_indent.h
@@ -28,12 +28,15 @@ class FixIndent : public Fix {
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
- double thermo(int);
+ double compute_scalar();
+ double compute_vector(int);
private:
int istyle,scaleflag,radflag,thermo_flag,eflag_enable;
- double k,k3,eng,etotal;
+ double k,k3;
double x0,y0,z0,r0_stop,r0_start;
+ int indenter_flag;
+ double indenter[4],indenter_all[4];
int cdim;
double c1,c2;
double vx,vy,vz;
diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp
index d8b2076d88..d3ff39ed12 100644
--- a/src/fix_nph.cpp
+++ b/src/fix_nph.cpp
@@ -46,6 +46,8 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
restart_global = 1;
pressure_every = 1;
box_change = 1;
+ scalar_flag = 1;
+ scalar_vector_freq = 1;
double p_period[3];
if (strcmp(arg[3],"xyz") == 0) {
@@ -306,7 +308,7 @@ void FixNPH::init()
void FixNPH::setup()
{
- p_target[0] = p_start[0]; // needed by thermo() method
+ p_target[0] = p_start[0]; // used by compute_scalar()
p_target[1] = p_start[1];
p_target[2] = p_start[2];
@@ -767,7 +769,7 @@ int FixNPH::modify_param(int narg, char **arg)
/* ---------------------------------------------------------------------- */
-double FixNPH::thermo(int n)
+double FixNPH::compute_scalar()
{
double volume;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
@@ -781,6 +783,5 @@ double FixNPH::thermo(int n)
energy += 0.5*nkt*omega_dot[i]*omega_dot[i] /
(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
- if (n == 0) return energy;
- else return 0.0;
+ return energy;
}
diff --git a/src/fix_nph.h b/src/fix_nph.h
index 60b79846fb..09ec8d540d 100644
--- a/src/fix_nph.h
+++ b/src/fix_nph.h
@@ -29,7 +29,7 @@ class FixNPH : public Fix {
void final_integrate();
void initial_integrate_respa(int,int);
void final_integrate_respa(int);
- double thermo(int);
+ double compute_scalar();
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp
index f57df6fc25..a0dd52261b 100644
--- a/src/fix_npt.cpp
+++ b/src/fix_npt.cpp
@@ -45,6 +45,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
restart_global = 1;
pressure_every = 1;
box_change = 1;
+ scalar_flag = 1;
+ scalar_vector_freq = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
@@ -307,7 +309,7 @@ void FixNPT::init()
void FixNPT::setup()
{
- t_target = t_start; // used by thermo()
+ t_target = t_start; // used by compute_scalar()
p_target[0] = p_start[0];
p_target[1] = p_start[1];
p_target[2] = p_start[2];
@@ -796,7 +798,7 @@ int FixNPT::modify_param(int narg, char **arg)
/* ---------------------------------------------------------------------- */
-double FixNPT::thermo(int n)
+double FixNPT::compute_scalar()
{
double ke = temperature->dof * boltz * t_target;
double keplus = atom->natoms * boltz * t_target;
@@ -812,6 +814,5 @@ double FixNPT::thermo(int n)
energy += 0.5*keplus*omega_dot[i]*omega_dot[i] /
(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
- if (n == 0) return energy;
- else return 0.0;
+ return energy;
}
diff --git a/src/fix_npt.h b/src/fix_npt.h
index 0a32120ba4..687756cddc 100644
--- a/src/fix_npt.h
+++ b/src/fix_npt.h
@@ -29,7 +29,7 @@ class FixNPT : public Fix {
virtual void final_integrate();
void initial_integrate_respa(int, int);
void final_integrate_respa(int);
- double thermo(int);
+ double compute_scalar();
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp
index 7a56cbd23b..859080eb17 100644
--- a/src/fix_nvt.cpp
+++ b/src/fix_nvt.cpp
@@ -39,6 +39,8 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
if (narg < 6) error->all("Illegal fix nvt command");
restart_global = 1;
+ scalar_flag = 1;
+ scalar_vector_freq = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
@@ -129,7 +131,7 @@ void FixNVT::init()
void FixNVT::setup()
{
- t_target = t_start; // needed by thermo() method
+ t_target = t_start; // used by compute_scalar()
t_current = temperature->compute_scalar();
}
@@ -377,11 +379,9 @@ void FixNVT::reset_target(double t_new)
/* ---------------------------------------------------------------------- */
-double FixNVT::thermo(int n)
+double FixNVT::compute_scalar()
{
double ke = temperature->dof * force->boltz * t_target;
double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq));
-
- if (n == 0) return energy;
- else return 0.0;
+ return energy;
}
diff --git a/src/fix_nvt.h b/src/fix_nvt.h
index 7cd50c68a6..12818afd1c 100644
--- a/src/fix_nvt.h
+++ b/src/fix_nvt.h
@@ -29,7 +29,7 @@ class FixNVT : public Fix {
virtual void final_integrate();
virtual void initial_integrate_respa(int,int);
void final_integrate_respa(int);
- double thermo(int);
+ double compute_scalar();
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp
index 29ce6fb0ff..a50474bf23 100644
--- a/src/fix_orient_fcc.cpp
+++ b/src/fix_orient_fcc.cpp
@@ -46,6 +46,9 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
if (narg != 11) error->all("Illegal fix orient/fcc command");
+ scalar_flag = 1;
+ scalar_vector_freq = 1;
+
nstats = atoi(arg[3]);
direction_of_motion = atoi(arg[4]);
a = atof(arg[5]);
@@ -194,7 +197,6 @@ void FixOrientFCC::init_list(int id, NeighList *ptr)
void FixOrientFCC::setup()
{
- eflag_enable = 1;
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(1);
else {
@@ -202,7 +204,6 @@ void FixOrientFCC::setup()
post_force_respa(1,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
- eflag_enable = 0;
}
/* ---------------------------------------------------------------------- */
@@ -211,16 +212,12 @@ void FixOrientFCC::post_force(int vflag)
{
int i,j,k,ii,jj,inum,jnum,m,n,nn,nsort,id_self;
int *ilist,*jlist,*numneigh,**firstneigh;
- double edelta,added_energy,omega;
+ double edelta,omega;
double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other;
double dxi[3];
double *dxiptr;
bool found_myself;
- bool eflag = false;
- if (eflag_enable) eflag = true;
- else if (output->next_thermo == update->ntimestep) eflag = true;
-
// set local ptrs
double **x = atom->x;
@@ -246,7 +243,7 @@ void FixOrientFCC::post_force(int vflag)
// loop over owned atoms and build Nbr data structure of neighbors
// use full neighbor list
- if (eflag) added_energy = 0.0;
+ added_energy = 0.0;
int count = 0;
int mincount = BIG;
int maxcount = 0;
@@ -319,16 +316,16 @@ void FixOrientFCC::post_force(int vflag)
if (xi_total < xi0) {
nbr[i].duxi = 0.0;
- if (eflag) edelta = 0.0;
+ edelta = 0.0;
} else if (xi_total > xi1) {
nbr[i].duxi = 0.0;
- if (eflag) edelta = Vxi;
+ edelta = Vxi;
} else {
omega = (0.5*PI)*(xi_total-xi0) / (xi1-xi0);
nbr[i].duxi = PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0));
- if (eflag) edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
+ edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
}
- if (eflag) added_energy += edelta;
+ added_energy += edelta;
}
if (maxcount) delete [] sort;
@@ -391,11 +388,6 @@ void FixOrientFCC::post_force(int vflag)
}
}
- // sum energy across procs
-
- if (eflag)
- MPI_Allreduce(&added_energy,&total_added_e,1,MPI_DOUBLE,MPI_SUM,world);
-
// print statistics every nstats timesteps
if (nstats && update->ntimestep % nstats == 0) {
@@ -433,10 +425,10 @@ void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop)
/* ---------------------------------------------------------------------- */
-double FixOrientFCC::thermo(int n)
+double FixOrientFCC::compute_scalar()
{
- if (n == 0) return total_added_e;
- else return 0.0;
+ MPI_Allreduce(&added_energy,&total_added_e,1,MPI_DOUBLE,MPI_SUM,world);
+ return total_added_e;
}
/* ---------------------------------------------------------------------- */
diff --git a/src/fix_orient_fcc.h b/src/fix_orient_fcc.h
index 1f19a2c0b3..3293ba6700 100644
--- a/src/fix_orient_fcc.h
+++ b/src/fix_orient_fcc.h
@@ -45,7 +45,7 @@ class FixOrientFCC : public Fix {
void setup();
void post_force(int);
void post_force_respa(int, int, int);
- double thermo(int);
+ double compute_scalar();
int pack_comm(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
double memory_usage();
@@ -64,9 +64,8 @@ class FixOrientFCC : public Fix {
char *xifilename, *chifilename; // file names for 2 crystal orientations
bool use_xismooth;
- int eflag_enable;
double Rxi[12][3],Rchi[12][3],half_xi_chi_vec[2][6][3];
- double xiid,xi0,xi1,xicutoffsq,cutsq,total_added_e;
+ double xiid,xi0,xi1,xicutoffsq,cutsq,added_energy,total_added_e;
int half_fcc_nn,nmax;
Nbr *nbr;
diff --git a/src/fix_set_force.cpp b/src/fix_set_force.cpp
index 5ef89e9069..6474e095c3 100644
--- a/src/fix_set_force.cpp
+++ b/src/fix_set_force.cpp
@@ -28,6 +28,10 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 6) error->all("Illegal fix setforce command");
+ vector_flag = 1;
+ size_vector = 3;
+ scalar_vector_freq = 1;
+
flagx = flagy = flagz = 1;
if (strcmp(arg[3],"NULL") == 0) flagx = 0;
else xvalue = atof(arg[3]);
@@ -86,6 +90,7 @@ void FixSetForce::post_force(int vflag)
int nlocal = atom->nlocal;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
+ force_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -111,6 +116,7 @@ void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
int nlocal = atom->nlocal;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
+ force_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -132,15 +138,16 @@ void FixSetForce::min_post_force(int vflag)
}
/* ----------------------------------------------------------------------
- allow setforce values to be printed with thermo output
- n = 1,2,3 = components of total force on fix group before reset
+ return components of total force on fix group before reset
------------------------------------------------------------------------- */
-double FixSetForce::thermo(int n)
+double FixSetForce::compute_vector(int n)
{
- if (n >= 1 && n <= 3) {
- double ftotal;
- MPI_Allreduce(&foriginal[n-1],&ftotal,1,MPI_DOUBLE,MPI_SUM,world);
- return ftotal;
- } else return 0.0;
+ // only sum across procs one time
+
+ if (force_flag == 0) {
+ MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world);
+ force_flag = 1;
+ }
+ return foriginal_all[n];
}
diff --git a/src/fix_set_force.h b/src/fix_set_force.h
index c6252fb6d5..15ea1a3b70 100644
--- a/src/fix_set_force.h
+++ b/src/fix_set_force.h
@@ -28,12 +28,13 @@ class FixSetForce : public Fix {
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
- double thermo(int);
+ double compute_vector(int);
private:
int flagx,flagy,flagz;
double xvalue,yvalue,zvalue;
- double foriginal[3];
+ double foriginal[3],foriginal_all[3];
+ int force_flag;
int nlevels_respa;
};
diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp
index 7d5f200030..bd645be70c 100644
--- a/src/fix_temp_rescale.cpp
+++ b/src/fix_temp_rescale.cpp
@@ -36,9 +36,13 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 8) error->all("Illegal fix temp/rescale command");
+
nevery = atoi(arg[3]);
if (nevery <= 0) error->all("Illegal fix temp/rescale command");
+ scalar_flag = 1;
+ scalar_vector_freq = nevery;
+
t_start = atof(arg[4]);
t_end = atof(arg[5]);
t_window = atof(arg[6]);
@@ -223,8 +227,7 @@ int FixTempRescale::modify_param(int narg, char **arg)
/* ---------------------------------------------------------------------- */
-double FixTempRescale::thermo(int n)
+double FixTempRescale::compute_scalar()
{
- if (n == 0) return energy;
- else return 0.0;
+ return energy;
}
diff --git a/src/fix_temp_rescale.h b/src/fix_temp_rescale.h
index 53d6db4188..11137bbbc4 100644
--- a/src/fix_temp_rescale.h
+++ b/src/fix_temp_rescale.h
@@ -25,7 +25,7 @@ class FixTempRescale : public Fix {
int setmask();
void init();
void end_of_step();
- double thermo(int);
+ double compute_scalar();
int modify_param(int, char **);
private:
diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp
index a0ff286373..5b43044287 100644
--- a/src/fix_wall_lj126.cpp
+++ b/src/fix_wall_lj126.cpp
@@ -35,6 +35,11 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 8) error->all("Illegal fix wall/lj126 command");
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 3;
+ scalar_vector_freq = 1;
+
if (strcmp(arg[3],"xlo") == 0) {
dim = 0;
side = -1;
@@ -101,7 +106,6 @@ void FixWallLJ126::init()
void FixWallLJ126::setup()
{
- eflag_enable = 1;
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(1);
else {
@@ -109,14 +113,12 @@ void FixWallLJ126::setup()
post_force_respa(1,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
- eflag_enable = 0;
}
/* ---------------------------------------------------------------------- */
void FixWallLJ126::min_setup()
{
- eflag_enable = 1;
post_force(1);
}
@@ -124,17 +126,14 @@ void FixWallLJ126::min_setup()
void FixWallLJ126::post_force(int vflag)
{
- bool eflag = false;
- if (eflag_enable) eflag = true;
- else if (output->next_thermo == update->ntimestep) eflag = true;
-
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
- double delta,rinv,r2inv,r6inv,eng;
- if (eflag) eng = 0.0;
+ double delta,rinv,r2inv,r6inv,fwall;
+ wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+ wall_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -145,11 +144,11 @@ void FixWallLJ126::post_force(int vflag)
rinv = 1.0/delta;
r2inv = rinv*rinv;
r6inv = r2inv*r2inv*r2inv;
- f[i][dim] -= r6inv*(coeff1*r6inv - coeff2) * side;
- if (eflag) eng += r6inv*(coeff3*r6inv - coeff4) - offset;
+ fwall = r6inv*(coeff1*r6inv - coeff2) * side;
+ f[i][dim] -= fwall;
+ wall[0] += r6inv*(coeff3*r6inv - coeff4) - offset;
+ wall[dim] += fwall;
}
-
- if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
@@ -166,10 +165,32 @@ void FixWallLJ126::min_post_force(int vflag)
post_force(vflag);
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ energy of wall interaction
+------------------------------------------------------------------------- */
+
+double FixWallLJ126::compute_scalar()
+{
+ // only sum across procs one time
-double FixWallLJ126::thermo(int n)
+ if (wall_flag == 0) {
+ MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ wall_flag = 1;
+ }
+ return wall_all[0];
+}
+
+/* ----------------------------------------------------------------------
+ components of force on wall
+------------------------------------------------------------------------- */
+
+double FixWallLJ126::compute_vector(int n)
{
- if (n == 0) return etotal;
- else return 0.0;
+ // only sum across procs one time
+
+ if (wall_flag == 0) {
+ MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ wall_flag = 1;
+ }
+ return wall_all[n];
}
diff --git a/src/fix_wall_lj126.h b/src/fix_wall_lj126.h
index e5087f2d86..2118af4f43 100644
--- a/src/fix_wall_lj126.h
+++ b/src/fix_wall_lj126.h
@@ -28,13 +28,15 @@ class FixWallLJ126 : public Fix {
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
- double thermo(int);
+ double compute_scalar();
+ double compute_vector(int);
private:
int dim,side,thermo_flag,eflag_enable;
double coord,epsilon,sigma,cutoff;
double coeff1,coeff2,coeff3,coeff4,offset;
- double etotal;
+ double wall[4],wall_all[4];
+ int wall_flag;
int nlevels_respa;
};
diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp
index 7d40af983d..99d020ebb6 100644
--- a/src/fix_wall_lj93.cpp
+++ b/src/fix_wall_lj93.cpp
@@ -31,6 +31,11 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 8) error->all("Illegal fix wall/lj93 command");
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 3;
+ scalar_vector_freq = 1;
+
if (strcmp(arg[3],"xlo") == 0) {
dim = 0;
side = -1;
@@ -98,7 +103,6 @@ void FixWallLJ93::init()
void FixWallLJ93::setup()
{
- eflag_enable = 1;
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(1);
else {
@@ -106,14 +110,12 @@ void FixWallLJ93::setup()
post_force_respa(1,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
- eflag_enable = 0;
}
/* ---------------------------------------------------------------------- */
void FixWallLJ93::min_setup()
{
- eflag_enable = 1;
post_force(1);
}
@@ -121,17 +123,14 @@ void FixWallLJ93::min_setup()
void FixWallLJ93::post_force(int vflag)
{
- bool eflag = false;
- if (eflag_enable) eflag = true;
- else if (output->next_thermo == update->ntimestep) eflag = true;
-
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
- double delta,rinv,r2inv,r4inv,r10inv,eng;
- if (eflag) eng = 0.0;
+ double delta,rinv,r2inv,r4inv,r10inv,fwall;
+ wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+ wall_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@@ -143,11 +142,11 @@ void FixWallLJ93::post_force(int vflag)
r2inv = rinv*rinv;
r4inv = r2inv*r2inv;
r10inv = r4inv*r4inv*r2inv;
- f[i][dim] -= (coeff1*r10inv - coeff2*r4inv) * side;
- if (eflag) eng += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+ fwall = (coeff1*r10inv - coeff2*r4inv) * side;
+ f[i][dim] -= fwall;
+ wall[0] += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+ wall[dim] += fwall;
}
-
- if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
@@ -164,10 +163,34 @@ void FixWallLJ93::min_post_force(int vflag)
post_force(vflag);
}
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ energy of wall interaction
+------------------------------------------------------------------------- */
+
+double FixWallLJ93::compute_scalar()
+{
+ // only sum across procs one time
+
+ if (wall_flag == 0) {
+ MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ wall_flag = 1;
+ }
+ return wall_all[0];
+}
+
+/* ----------------------------------------------------------------------
+ components of force on wall
+------------------------------------------------------------------------- */
-double FixWallLJ93::thermo(int n)
+double FixWallLJ93::compute_vector(int n)
{
- if (n == 0) return etotal;
- else return 0.0;
+ // only sum across procs one time
+
+ if (wall_flag == 0) {
+ MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ wall_flag = 1;
+ }
+ return wall_all[n];
}
+
+
diff --git a/src/fix_wall_lj93.h b/src/fix_wall_lj93.h
index ae846ba6d0..9a6d2cc6f5 100644
--- a/src/fix_wall_lj93.h
+++ b/src/fix_wall_lj93.h
@@ -29,13 +29,15 @@ class FixWallLJ93 : public Fix {
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
- double thermo(int);
+ double compute_scalar();
+ double compute_vector(int);
private:
int dim,side,thermo_flag,eflag_enable;
double coord,epsilon,sigma,cutoff;
double coeff1,coeff2,coeff3,coeff4,offset;
- double etotal;
+ double wall[4],wall_all[4];
+ int wall_flag;
int nlevels_respa;
};
diff --git a/src/modify.cpp b/src/modify.cpp
index 3cb076ee7f..e0fb599bb8 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -235,14 +235,14 @@ void Modify::end_of_step()
/* ----------------------------------------------------------------------
thermo energy call only for relevant fixes
called by Thermo class
- arg to Fix thermo() is 0, so fix will return its energy contribution
+ compute_scalar() is fix call to return energy
------------------------------------------------------------------------- */
double Modify::thermo_energy()
{
double energy = 0.0;
for (int i = 0; i < n_thermo_energy; i++)
- energy += fix[list_thermo_energy[i]]->thermo(0);
+ energy += fix[list_thermo_energy[i]]->compute_scalar();
return energy;
}
diff --git a/src/thermo.cpp b/src/thermo.cpp
index 2304a79a14..b6c2bb9fd4 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -793,12 +793,20 @@ void Thermo::parse_fields(char *str)
for (int ic = 0; ic < modify->compute[n]->npre; ic++)
int tmp = add_compute(modify->compute[n]->id_pre[ic],
arg_object[nfield]);
+
field2object[nfield] = add_compute(id,arg_object[nfield]);
addfield(copy,&Thermo::compute_compute,FLOAT);
} else if (word[0] == 'f') {
n = modify->find_fix(id);
if (n < 0) error->all("Could not find thermo custom fix ID");
+ if (arg_object[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
+ error->all("Thermo fix ID does not compute scalar info");
+ if (arg_object[nfield] > 0 && modify->fix[n]->vector_flag == 0)
+ error->all("Thermo fix ID does not compute vector info");
+ if (arg_object[nfield] > 0 &&
+ arg_object[nfield] > modify->fix[n]->size_vector)
+ error->all("Thermo fix ID vector is not large enough");
field2object[nfield] = add_fix(id);
addfield(copy,&Thermo::compute_fix,FLOAT);
@@ -1027,7 +1035,8 @@ void Thermo::compute_compute()
void Thermo::compute_fix()
{
int index = field2object[ifield];
- dvalue = fixes[index]->thermo(arg_object[ifield]);
+ if (arg_object[ifield] == 0) dvalue = fixes[index]->compute_scalar();
+ else dvalue = fixes[index]->compute_vector(arg_object[ifield]-1);
if (normflag) dvalue /= natoms;
}
diff --git a/src/variable.cpp b/src/variable.cpp
index cb7311647f..55563edaed 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -24,6 +24,7 @@
#include "domain.h"
#include "modify.h"
#include "compute.h"
+#include "fix.h"
#include "output.h"
#include "thermo.h"
#include "memory.h"
@@ -407,6 +408,7 @@ void Variable::copy(int narg, char **from, char **to)
group function = mass(group), xcm(group,x), ...
atom vector = x[123], y[3], vx[34], ...
compute vector = c_mytemp[0], c_thermo_press[3], ...
+ fix vector = f_indent[0], f_setforce[3], ...
numbers start with a digit or "." or "-" (no parens or brackets)
keywords start with a lowercase letter (no parens or brackets)
functions contain ()
@@ -414,7 +416,7 @@ void Variable::copy(int narg, char **from, char **to)
vectors contain []
single arg must be integer
for atom vectors, arg is global ID of atom
- for compute vectors, 0 is the scalar value, 1-N are vector values
+ for compute/fix vectors, 0 is the scalar value, 1-N are vector values
see lists of valid functions & vectors below
return answer = value of string
------------------------------------------------------------------------- */
@@ -662,6 +664,8 @@ double Variable::evaluate(char *str, Tree *tree)
// index = everything between [] evaluated as integer
// if vector name starts with "c_", trailing chars are compute ID
// check if compute ID exists, invoke it with index as arg
+ // if vector name starts with "cf_", trailing chars are fix ID
+ // check if fix ID exists, call its scalar or vector fn with index as arg
// else is atom vector
// find which proc owns atom via atom->map()
// grab atom-based value with index as global atom ID
@@ -696,7 +700,7 @@ double Variable::evaluate(char *str, Tree *tree)
delete [] id;
modify->compute[icompute]->init();
- // call compute() if index = 0, else compute_vector()
+ // call compute_scalar() if index = 0, else compute_vector()
// make pre-call to Compute object's pre-compute(s) if defined
int index = atoi(arg);
@@ -713,7 +717,7 @@ double Variable::evaluate(char *str, Tree *tree)
answer = modify->compute[icompute]->compute_scalar();
} else if (index > 0) {
if (modify->compute[icompute]->vector_flag == 0)
- error->all("Variable compute ID does not compute scalar info");
+ error->all("Variable compute ID does not compute vector info");
if (index > modify->compute[icompute]->size_vector)
error->all("Variable compute ID vector is not large enough");
if (modify->compute[icompute]->npre)
@@ -727,6 +731,33 @@ double Variable::evaluate(char *str, Tree *tree)
answer = modify->compute[icompute]->vector[index-1];
} else error->all("Invalid compute ID index in variable");
+ } else if (strncmp(vector,"f_",2) == 0) {
+ n = strlen(vector) - 2 + 1;
+ char *id = new char[n];
+ strcpy(id,&vector[2]);
+
+ int ifix;
+ for (ifix = 0; ifix < modify->nfix; ifix++)
+ if (strcmp(id,modify->fix[ifix]->id) == 0) break;
+ if (ifix == modify->nfix)
+ error->all("Invalid fix ID in variable");
+ delete [] id;
+
+ // call compute_scalar() if index = 0, else compute_vector()
+
+ int index = atoi(arg);
+ if (index == 0) {
+ if (modify->fix[ifix]->scalar_flag == 0)
+ error->all("Variable fix ID does not compute scalar info");
+ answer = modify->fix[ifix]->compute_scalar();
+ } else if (index > 0) {
+ if (modify->fix[ifix]->vector_flag == 0)
+ error->all("Variable fix ID does not compute vector info");
+ if (index > modify->fix[ifix]->size_vector)
+ error->all("Variable fix ID vector is not large enough");
+ answer = modify->fix[ifix]->compute_vector(index-1);
+ } else error->all("Invalid fix ID index in variable");
+
} else if (tree) {
if (strlen(arg)) error->all("Invalid atom vector in variable");
--
GitLab