These are additional compute styles in USER packages, which can be diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index feb2faf9f6..3775f339a8 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -641,6 +641,7 @@ KOKKOS, o = USER-OMP, t = OPT. "erotate/sphere"_compute_erotate_sphere.html, "erotate/sphere/atom"_compute_erotate_sphere_atom.html, "event/displace"_compute_event_displace.html, +"force/molecule"_compute_force_molecule.html, "group/group"_compute_group_group.html, "gyration"_compute_gyration.html, "gyration/molecule"_compute_gyration_molecule.html, @@ -681,6 +682,7 @@ KOKKOS, o = USER-OMP, t = OPT. "temp/sphere"_compute_temp_sphere.html, "ti"_compute_ti.html, "vacf"_compute_vacf.html, +"vcm/molecule"_compute_vcm_molecule.html, "voronoi/atom"_compute_voronoi_atom.html :tb(c=6,ea=c) These are additional compute styles in USER packages, which can be diff --git a/doc/compute_force_molecule.html b/doc/compute_force_molecule.html new file mode 100644 index 0000000000..0ef54b3feb --- /dev/null +++ b/doc/compute_force_molecule.html @@ -0,0 +1,63 @@ + +
Syntax: +
+compute ID group-ID force/molecule ++
Examples: +
+compute 1 fluid force/molecule ++
Description: +
+Define a computation that calculates the total force on individual +molecules, as a sum over the forces on atoms in the molecule. The +x,y,z components of the force on each molecule are computed. +
+The force on a particular molecule is only computed if one or more of +its atoms are in the specified group. Normally all atoms in the +molecule should be in the group, however this is not required. LAMMPS +will warn you if this is not the case. Only atoms in the group +contribute to the force on the molecule. +
+The ordering of per-molecule quantities produced by this compute is +consistent with the ordering produced by other compute commands that +generate per-molecule datums. Conceptually, them molecule IDs will be +in ascending order for any molecule with one or more of its atoms in +the specified group. +
+Output info: +
+This compute calculates a global array where the number of rows = +Nmolecules and the number of columns = 3 for the fx,fy,fz components +of force on each molecule. These values can be accessed by any +command that uses global array values from a compute as input. See +Section_howto 15 for an overview of +LAMMPS output options. +
+The array values are "intensive". The array values will be in +distance units. +
+Restrictions: none +
+Related commands: +
+ +Default: none +
+ diff --git a/doc/compute_force_molecule.txt b/doc/compute_force_molecule.txt new file mode 100644 index 0000000000..5021788158 --- /dev/null +++ b/doc/compute_force_molecule.txt @@ -0,0 +1,58 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute force/molecule command :h3 + +[Syntax:] + +compute ID group-ID force/molecule :pre + +ID, group-ID are documented in "compute"_compute.html command +force/molecule = style name of this compute command :ul + +[Examples:] + +compute 1 fluid force/molecule :pre + +[Description:] + +Define a computation that calculates the total force on individual +molecules, as a sum over the forces on atoms in the molecule. The +x,y,z components of the force on each molecule are computed. + +The force on a particular molecule is only computed if one or more of +its atoms are in the specified group. Normally all atoms in the +molecule should be in the group, however this is not required. LAMMPS +will warn you if this is not the case. Only atoms in the group +contribute to the force on the molecule. + +The ordering of per-molecule quantities produced by this compute is +consistent with the ordering produced by other compute commands that +generate per-molecule datums. Conceptually, them molecule IDs will be +in ascending order for any molecule with one or more of its atoms in +the specified group. + +[Output info:] + +This compute calculates a global array where the number of rows = +Nmolecules and the number of columns = 3 for the fx,fy,fz components +of force on each molecule. These values can be accessed by any +command that uses global array values from a compute as input. See +"Section_howto 15"_Section_howto.html#howto_15 for an overview of +LAMMPS output options. + +The array values are "intensive". The array values will be in +distance "units"_units.html. + +[Restrictions:] none + +[Related commands:] + +"compute vcm/molecule"_compute_vcm_molecule.html + +[Default:] none diff --git a/doc/compute_vcm_molecule.html b/doc/compute_vcm_molecule.html new file mode 100644 index 0000000000..0bc881b1d7 --- /dev/null +++ b/doc/compute_vcm_molecule.html @@ -0,0 +1,65 @@ + +Syntax: +
+compute ID group-ID vcm/molecule ++
Examples: +
+compute 1 fluid vcm/molecule ++
Description: +
+Define a computation that calculates the center-of-mass velocity of +individual molecules. The x,y,z components of the velocity of each +molecule are computed. This is calcualted by summing mass*velocity +for each atom in the molecule and dividing the sum by the total mass +of the molecule. +
+The velocity of a particular molecule is only computed if one or more +of its atoms are in the specified group. Normally all atoms in the +molecule should be in the group, however this is not required. LAMMPS +will warn you if this is not the case. Only atoms in the group +contribute to the velocity calculation for the molecule. +
+The ordering of per-molecule quantities produced by this compute is +consistent with the ordering produced by other compute commands that +generate per-molecule datums. Conceptually, them molecule IDs will be +in ascending order for any molecule with one or more of its atoms in +the specified group. +
+Output info: +
+This compute calculates a global array where the number of rows = +Nmolecules and the number of columns = 3 for the vx,vy,vz components +of the center-of-mass velocity of each molecule. These values can be +accessed by any command that uses global array values from a compute +as input. See Section_howto 15 for an +overview of LAMMPS output options. +
+The array values are "intensive". The array values will be in +distance units. +
+Restrictions: none +
+Related commands: +
+ +Default: none +
+ diff --git a/doc/compute_vcm_molecule.txt b/doc/compute_vcm_molecule.txt new file mode 100644 index 0000000000..24a8fc5031 --- /dev/null +++ b/doc/compute_vcm_molecule.txt @@ -0,0 +1,60 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute vcm/molecule command :h3 + +[Syntax:] + +compute ID group-ID vcm/molecule :pre + +ID, group-ID are documented in "compute"_compute.html command +vcm/molecule = style name of this compute command :ul + +[Examples:] + +compute 1 fluid vcm/molecule :pre + +[Description:] + +Define a computation that calculates the center-of-mass velocity of +individual molecules. The x,y,z components of the velocity of each +molecule are computed. This is calcualted by summing mass*velocity +for each atom in the molecule and dividing the sum by the total mass +of the molecule. + +The velocity of a particular molecule is only computed if one or more +of its atoms are in the specified group. Normally all atoms in the +molecule should be in the group, however this is not required. LAMMPS +will warn you if this is not the case. Only atoms in the group +contribute to the velocity calculation for the molecule. + +The ordering of per-molecule quantities produced by this compute is +consistent with the ordering produced by other compute commands that +generate per-molecule datums. Conceptually, them molecule IDs will be +in ascending order for any molecule with one or more of its atoms in +the specified group. + +[Output info:] + +This compute calculates a global array where the number of rows = +Nmolecules and the number of columns = 3 for the vx,vy,vz components +of the center-of-mass velocity of each molecule. These values can be +accessed by any command that uses global array values from a compute +as input. See "Section_howto 15"_Section_howto.html#howto_15 for an +overview of LAMMPS output options. + +The array values are "intensive". The array values will be in +distance "units"_units.html. + +[Restrictions:] none + +[Related commands:] + +"compute force/molecule"_compute_force_molecule.html + +[Default:] none diff --git a/doc/lattice.html b/doc/lattice.html index e3dd79db8a..2e75e2157d 100644 --- a/doc/lattice.html +++ b/doc/lattice.html @@ -82,6 +82,15 @@ are the edge vectors of the unit cell. This is the nomenclature for unit cell they determine does not have to be a "primitive cell" of minimum volume. +Note that the lattice command can be used multiple times in an input +script. Each time it is invoked, the lattice attributes are +re-defined and are used for all subsequent commands (that use lattice +attributes). For example, a sequence of lattice, +region, and create_atoms commands +can be repeated multiple times to build a poly-crystalline model with +different geometric regions populated with atoms in different lattice +orientations. +
A lattice of style none does not define a unit cell and basis set, @@ -168,61 +177,73 @@ mapping it into the simulation box. The dim argument is one of the the crystallographic direction in the lattice that you want to orient along that axis, specified as integers. E.g. "orient x 2 1 0" means the x-axis in the simulation box will be the [210] lattice -direction. The 3 lattice directions you specify must be mutually +direction, and similarly for y and z. The 3 lattice directions you +specify do not have to be unit vectors, but they must be mutually orthogonal and obey the right-hand rule, i.e. (X cross Y) points in -the Z direction. Note that this description is really only valid for -orthogonal lattices. If you are using the more general lattice style -custom with non-orthogonal a1,a2,a3 vectors, then think of the 3 -orient options as creating a 3x3 rotation matrix which is applied to -a1,a2,a3 to rotate the original unit cell to a new orientation in the -simulation box. +the Z direction. +
+IMPORTANT NOTE: The preceding paragraph describing lattice directions +is only valid for orthogonal cubic unit cells (or square in 2d). If +you are using a hcp or hex lattice or the more general lattice +style custom with non-orthogonal a1,a2,a3 vectors, then you should +think of the 3 orient vectors as creating a 3x3 rotation matrix +which is applied to a1,a2,a3 to rotate the original unit cell to a new +orientation in the simulation box.
Several LAMMPS commands have the option to use distance units that are -inferred from "lattice spacing" in the x,y,z box directions. E.g. the -region command can create a block of size 10x20x20, -where 10 means 10 lattice spacings in the x direction. +inferred from "lattice spacings" in the x,y,z box directions. +E.g. the region command can create a block of size +10x20x20, where 10 means 10 lattice spacings in the x direction.
-The spacing option sets the 3 lattice spacings directly. All must -be non-zero (use 1.0 for dz in a 2d simulation). The specified values -are multiplied by the multiplicative factor described above that is -associated with the scale factor. Thus a spacing of 1.0 means one -unit cell independent of the scale factor. This option can be useful -if the spacings LAMMPS computes are inconvenient to use in subsequent -commands, which can be the case for non-orthogonal or rotated -lattices. +
IMPORTANT NOTE: Though they are called lattice spacings, all the +commands that have a "units lattice" option, simply use the 3 values +as scale factors on the distance units defined by the +units command. Thus if you do not like the lattice +spacings computed by LAMMPS (e.g. for a non-orthogonal or rotated unit +cell), you can define the 3 values to be whatever you wish, via the +spacing option.
If the spacing option is not specified, the lattice spacings are computed by LAMMPS in the following way. A unit cell of the lattice -is mapped into the simulation box (scaled, shifted, rotated), so that -it now has (perhaps) a modified size and orientation. The lattice -spacing in X is defined as the difference between the min/max extent -of the x coordinates of the 8 corner points of the modified unit cell. -Similarly, the Y and Z lattice spacings are defined as the difference -in the min/max of the y and z coordinates. +is mapped into the simulation box (scaled and rotated), so that it now +has (perhaps) a modified size and orientation. The lattice spacing in +X is defined as the difference between the min/max extent of the x +coordinates of the 8 corner points of the modified unit cell (4 in +2d). Similarly, the Y and Z lattice spacings are defined as the +difference in the min/max of the y and z coordinates.
-Note that if the unit cell is orthogonal with axis-aligned edges (not -rotated via the orient keyword), then the lattice spacings in each +
Note that if the unit cell is orthogonal with axis-aligned edges (no +rotation via the orient keyword), then the lattice spacings in each dimension are simply the scale factor (described above) multiplied by the length of a1,a2,a3. Thus a hex style lattice with a scale factor of 3.0 Angstroms, would have a lattice spacing of 3.0 in x and 3*sqrt(3.0) in y.
IMPORTANT NOTE: For non-orthogonal unit cells and/or when a rotation -is applied via the orient keyword, then the lattice spacings may be -less intuitive. In particular, in these cases, there is no guarantee -that the lattice spacing is an integer multiple of the periodicity of -the lattice in that direction. Thus, if you create an orthogonal -periodic simulation box whose size in a dimension is a multiple of the -lattice spacing, and then fill it with atoms via the -create_atoms command, you will NOT necessarily -create a periodic system. I.e. atoms may overlap incorrectly at the -faces of the simulation box. +is applied via the orient keyword, then the lattice spacings +computed by LAMMPS are typically less intuitive. In particular, in +these cases, there is no guarantee that a particular lattice spacing +is an integer multiple of the periodicity of the lattice in that +direction. Thus, if you create an orthogonal periodic simulation box +whose size in a dimension is a multiple of the lattice spacing, and +then fill it with atoms via the create_atoms +command, you will NOT necessarily create a periodic system. +I.e. atoms may overlap incorrectly at the faces of the simulation box.
-Regardless of these issues, the values of the lattice spacings LAMMPS -calculates are printed out, so their effect in commands that use the -spacings should be decipherable. +
The spacing option sets the 3 lattice spacings directly. All must +be non-zero (use 1.0 for dz in a 2d simulation). The specified values +are multiplied by the multiplicative factor described above that is +associated with the scale factor. Thus a spacing of 1.0 means one +unit cell edge length independent of the scale factor. As mentioned +above, this option can be useful if the spacings LAMMPS computes are +inconvenient to use in subsequent commands, which can be the case for +non-orthogonal or rotated lattices. +
+Note that whenever the lattice command is used, the values of the +lattice spacings LAMMPS calculates are printed out. Thus their effect +in commands that use the spacings should be decipherable.
The stride(x,y,z) function uses the current timestep to generate a new timestep. X,y >= 0 and z > 0 and x <= y are required. The generated -timesteps increase in increments of z, from x to y, I.e. it generates +timesteps increase in increments of z, from x to y, i.e. it generates the sequece x,x+z,x+2z,...,y. If y-x is not a multiple of z, then similar to the way a for loop operates, the last value will be one that does not exceed y. For any current timestep, the next timestep -in the sequence is returned. Thus if stagger(1000,2000,100) is used +in the sequence is returned. Thus if stride(1000,2000,100) is used in a variable by the dump_modify every command, it will generate the sequence of output timesteps:
1000,1100,1200, ... ,1900,2000+
The stride2(x,y,z,a,b,c) function is similar to the stride() function +except it generates two sets of strided timesteps, one at a coarser +level and one at a finer level. Thus it is useful for debugging, +e.g. to produce output every timestep at the point in simulation when +a problem occurs. X,y >= 0 and z > 0 and x <= y are required, as are +a,b >= 0 and c > 0 and a < b. Also, a >= x and b <= y are required so +that the second stride is inside the first. The generated timesteps +increase in increments of z, starting at x, until a is reached. At +that point the timestep increases in increments of c, from a to b, +then after b, increments by z are resumed until y is reached. For any +current timestep, the next timestep in the sequence is returned. Thus +if stride(1000,2000,100,1350,1360,1) is used in a variable by the +dump_modify every command, it will generate the +sequence of output timesteps: +
+1000,1100,1200,1300,1350,1351,1352, ... 1359,1360,1400,1500, ... ,2000 +
The vdisplace(x,y) function takes 2 arguments: x = value0 and y =
velocity, and uses the elapsed time to change the value by a linear
displacement due to the applied velocity over the course of a run,
diff --git a/doc/variable.txt b/doc/variable.txt
index c2805fed8a..868c0d70b3 100644
--- a/doc/variable.txt
+++ b/doc/variable.txt
@@ -49,7 +49,7 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
math functions = sqrt(x), exp(x), ln(x), log(x), abs(x),
sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x),
random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x)
- ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
+ ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
group functions = count(group), mass(group), charge(group),
xcm(group,dim), vcm(group,dim), fcm(group,dim),
bound(group,dir), gyration(group), ke(group),
@@ -371,7 +371,7 @@ Constant: PI
Thermo keywords: vol, pe, ebond, etc
Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x%y,
Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x%y, x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x || y, !x
-Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
+Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
Group functions: count(ID), mass(ID), charge(ID), xcm(ID,dim), \
vcm(ID,dim), fcm(ID,dim), bound(ID,dir), \
gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), \
@@ -542,16 +542,33 @@ output timesteps:
The stride(x,y,z) function uses the current timestep to generate a new
timestep. X,y >= 0 and z > 0 and x <= y are required. The generated
-timesteps increase in increments of z, from x to y, I.e. it generates
+timesteps increase in increments of z, from x to y, i.e. it generates
the sequece x,x+z,x+2z,...,y. If y-x is not a multiple of z, then
similar to the way a for loop operates, the last value will be one
that does not exceed y. For any current timestep, the next timestep
-in the sequence is returned. Thus if stagger(1000,2000,100) is used
+in the sequence is returned. Thus if stride(1000,2000,100) is used
in a variable by the "dump_modify every"_dump_modify.html command, it
will generate the sequence of output timesteps:
1000,1100,1200, ... ,1900,2000 :pre
+The stride2(x,y,z,a,b,c) function is similar to the stride() function
+except it generates two sets of strided timesteps, one at a coarser
+level and one at a finer level. Thus it is useful for debugging,
+e.g. to produce output every timestep at the point in simulation when
+a problem occurs. X,y >= 0 and z > 0 and x <= y are required, as are
+a,b >= 0 and c > 0 and a < b. Also, a >= x and b <= y are required so
+that the second stride is inside the first. The generated timesteps
+increase in increments of z, starting at x, until a is reached. At
+that point the timestep increases in increments of c, from a to b,
+then after b, increments by z are resumed until y is reached. For any
+current timestep, the next timestep in the sequence is returned. Thus
+if stride(1000,2000,100,1350,1360,1) is used in a variable by the
+"dump_modify every"_dump_modify.html command, it will generate the
+sequence of output timesteps:
+
+1000,1100,1200,1300,1350,1351,1352, ... 1359,1360,1400,1500, ... ,2000 :pre
+
The vdisplace(x,y) function takes 2 arguments: x = value0 and y =
velocity, and uses the elapsed time to change the value by a linear
displacement due to the applied velocity over the course of a run,
diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h
index 9128d8e16f..68176fd668 100644
--- a/src/QEQ/fix_qeq.h
+++ b/src/QEQ/fix_qeq.h
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp
index 848ca501c3..6099560245 100644
--- a/src/QEQ/fix_qeq_dynamic.cpp
+++ b/src/QEQ/fix_qeq_dynamic.cpp
@@ -169,6 +169,7 @@ void FixQEqDynamic::pre_force(int vflag)
}
if (force->kspace) force->kspace->setup();
+
}
/* ---------------------------------------------------------------------- */
@@ -213,6 +214,7 @@ double FixQEqDynamic::compute_eneg()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
+ j &= NEIGHMASK;
delr[0] = x[i][0] - x[j][0];
delr[1] = x[i][1] - x[j][1];
diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp
index 5f25332c44..28499b5f8d 100644
--- a/src/QEQ/fix_qeq_point.cpp
+++ b/src/QEQ/fix_qeq_point.cpp
@@ -28,8 +28,8 @@
#include "neigh_request.h"
#include "update.h"
#include "force.h"
-#include "kspace.h"
#include "group.h"
+#include "kspace.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
@@ -53,8 +53,8 @@ void FixQEqPoint::init()
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
- neighbor->requests[irequest]->half = 1;
- neighbor->requests[irequest]->full = 0;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding");
@@ -146,6 +146,7 @@ void FixQEqPoint::compute_H()
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
+ j &= NEIGHMASK;
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp
index a2c8611af7..d4b061575c 100644
--- a/src/QEQ/fix_qeq_shielded.cpp
+++ b/src/QEQ/fix_qeq_shielded.cpp
@@ -28,8 +28,8 @@
#include "neigh_request.h"
#include "update.h"
#include "force.h"
-#include "kspace.h"
#include "group.h"
+#include "kspace.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
@@ -190,6 +190,7 @@ void FixQEqShielded::compute_H()
for( jj = 0; jj < jnum; jj++ ) {
j = jlist[jj];
+ j &= NEIGHMASK;
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];
diff --git a/src/compute_force_molecule.cpp b/src/compute_force_molecule.cpp
new file mode 100644
index 0000000000..d05b162e3f
--- /dev/null
+++ b/src/compute_force_molecule.cpp
@@ -0,0 +1,106 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_force_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeForceMolecule::ComputeForceMolecule(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal compute force/molecule command");
+
+ if (atom->molecular == 0)
+ error->all(FLERR,"Compute force/molecule requires molecular atom style");
+
+ array_flag = 1;
+ size_array_cols = 3;
+ extarray = 0;
+
+ // setup molecule-based data
+
+ nmolecules = molecules_in_group(idlo,idhi);
+ size_array_rows = nmolecules;
+
+ memory->create(force,nmolecules,3,"force/molecule:force");
+ memory->create(forceall,nmolecules,3,"force/molecule:forceall");
+ array = forceall;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeForceMolecule::~ComputeForceMolecule()
+{
+ memory->destroy(force);
+ memory->destroy(forceall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeForceMolecule::init()
+{
+ int ntmp = molecules_in_group(idlo,idhi);
+ if (ntmp != nmolecules)
+ error->all(FLERR,"Molecule count changed in compute force/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeForceMolecule::compute_array()
+{
+ tagint imol;
+ double massone;
+ double unwrap[3];
+
+ invoked_array = update->ntimestep;
+
+ for (int i = 0; i < nmolecules; i++)
+ force[i][0] = force[i][1] = force[i][2] = 0.0;
+
+ double **f = atom->f;
+ int *mask = atom->mask;
+ tagint *molecule = atom->molecule;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ imol = molecule[i];
+ if (molmap) imol = molmap[imol-idlo];
+ else imol--;
+ force[imol][0] += f[i][0];
+ force[imol][1] += f[i][1];
+ force[imol][2] += f[i][2];
+ }
+
+ MPI_Allreduce(&force[0][0],&forceall[0][0],3*nmolecules,
+ MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeForceMolecule::memory_usage()
+{
+ double bytes = 0;
+ if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
+ bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+ return bytes;
+}
diff --git a/src/compute_force_molecule.h b/src/compute_force_molecule.h
new file mode 100644
index 0000000000..ec90e3c9fe
--- /dev/null
+++ b/src/compute_force_molecule.h
@@ -0,0 +1,63 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(force/molecule,ComputeForceMolecule)
+
+#else
+
+#ifndef LMP_COMPUTE_FORCE_MOLECULE_H
+#define LMP_COMPUTE_FORCE_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeForceMolecule : public Compute {
+ public:
+ ComputeForceMolecule(class LAMMPS *, int, char **);
+ ~ComputeForceMolecule();
+ void init();
+ void compute_array();
+ double memory_usage();
+
+ private:
+ int nmolecules;
+ tagint idlo,idhi;
+
+ double **force,**forceall;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute force/molecule requires molecular atom style
+
+Self-explanatory.
+
+E: Molecule count changed in compute force/molecule
+
+Number of molecules must remain constant over time.
+
+*/
diff --git a/src/compute_vcm_molecule.cpp b/src/compute_vcm_molecule.cpp
new file mode 100644
index 0000000000..ab86d2ee0a
--- /dev/null
+++ b/src/compute_vcm_molecule.cpp
@@ -0,0 +1,147 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_vcm_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVCMMolecule::ComputeVCMMolecule(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal compute vcm/molecule command");
+
+ if (atom->molecular == 0)
+ error->all(FLERR,"Compute vcm/molecule requires molecular atom style");
+
+ array_flag = 1;
+ size_array_cols = 3;
+ extarray = 0;
+
+ // setup molecule-based data
+
+ nmolecules = molecules_in_group(idlo,idhi);
+ size_array_rows = nmolecules;
+
+ memory->create(massproc,nmolecules,"vcm/molecule:massproc");
+ memory->create(masstotal,nmolecules,"vcm/molecule:masstotal");
+ memory->create(vcm,nmolecules,3,"vcm/molecule:vcm");
+ memory->create(vcmall,nmolecules,3,"vcm/molecule:vcmall");
+ array = vcmall;
+
+ // compute masstotal for each molecule
+
+ int *mask = atom->mask;
+ tagint *molecule = atom->molecule;
+ int *type = atom->type;
+ double *mass = atom->mass;
+ double *rmass = atom->rmass;
+ int nlocal = atom->nlocal;
+
+ tagint imol;
+ double massone;
+
+ for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+ imol = molecule[i];
+ if (molmap) imol = molmap[imol-idlo];
+ else imol--;
+ massproc[imol] += massone;
+ }
+
+ MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVCMMolecule::~ComputeVCMMolecule()
+{
+ memory->destroy(massproc);
+ memory->destroy(masstotal);
+ memory->destroy(vcm);
+ memory->destroy(vcmall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVCMMolecule::init()
+{
+ int ntmp = molecules_in_group(idlo,idhi);
+ if (ntmp != nmolecules)
+ error->all(FLERR,"Molecule count changed in compute vcm/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVCMMolecule::compute_array()
+{
+ tagint imol;
+ double massone;
+ double unwrap[3];
+
+ invoked_array = update->ntimestep;
+
+ for (int i = 0; i < nmolecules; i++)
+ vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
+
+ double **v = atom->v;
+ int *mask = atom->mask;
+ tagint *molecule = atom->molecule;
+ int *type = atom->type;
+ imageint *image = atom->image;
+ double *mass = atom->mass;
+ double *rmass = atom->rmass;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+ imol = molecule[i];
+ if (molmap) imol = molmap[imol-idlo];
+ else imol--;
+ vcm[imol][0] += v[i][0] * massone;
+ vcm[imol][1] += v[i][1] * massone;
+ vcm[imol][2] += v[i][2] * massone;
+ }
+
+ MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nmolecules,
+ MPI_DOUBLE,MPI_SUM,world);
+ for (int i = 0; i < nmolecules; i++) {
+ vcmall[i][0] /= masstotal[i];
+ vcmall[i][1] /= masstotal[i];
+ vcmall[i][2] /= masstotal[i];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeVCMMolecule::memory_usage()
+{
+ double bytes = (bigint) nmolecules * 2 * sizeof(double);
+ if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
+ bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+ return bytes;
+}
diff --git a/src/compute_vcm_molecule.h b/src/compute_vcm_molecule.h
new file mode 100644
index 0000000000..b5bce320ca
--- /dev/null
+++ b/src/compute_vcm_molecule.h
@@ -0,0 +1,64 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(vcm/molecule,ComputeVCMMolecule)
+
+#else
+
+#ifndef LMP_COMPUTE_VCM_MOLECULE_H
+#define LMP_COMPUTE_VCM_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeVCMMolecule : public Compute {
+ public:
+ ComputeVCMMolecule(class LAMMPS *, int, char **);
+ ~ComputeVCMMolecule();
+ void init();
+ void compute_array();
+ double memory_usage();
+
+ private:
+ int nmolecules;
+ tagint idlo,idhi;
+
+ double *massproc,*masstotal;
+ double **vcm,**vcmall;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute vcm/molecule requires molecular atom style
+
+Self-explanatory.
+
+E: Molecule count changed in compute vcm/molecule
+
+Number of molecules must remain constant over time.
+
+*/
diff --git a/src/variable.cpp b/src/variable.cpp
index f371512128..f3cd960ea8 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -45,6 +45,7 @@ using namespace MathConst;
#define MAXLINE 256
#define CHUNK 1024
#define VALUELENGTH 64
+#define MAXFUNCARG 6
#define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
@@ -59,7 +60,7 @@ enum{ARG,OP};
enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
NOT,EQ,NE,LT,LE,GT,GE,AND,OR,
SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
- RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE,
+ RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE,STRIDE2,
VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY};
@@ -973,7 +974,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = atof(number);
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->extra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = atof(number);
@@ -1057,7 +1059,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1082,7 +1085,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1110,7 +1114,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1179,7 +1184,8 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->array = compute->vector_atom;
newtree->nstride = 1;
newtree->selfalloc = 0;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
// c_ID[i] = vector from per-atom array
@@ -1210,7 +1216,8 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->array = NULL;
newtree->nstride = compute->size_peratom_cols;
newtree->selfalloc = 0;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else error->all(FLERR,"Mismatched compute in variable formula");
@@ -1265,7 +1272,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1284,7 +1292,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1306,7 +1315,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1360,7 +1370,8 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->array = fix->vector_atom;
newtree->nstride = 1;
newtree->selfalloc = 0;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
// f_ID[i] = vector from per-atom array
@@ -1385,7 +1396,8 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->array = NULL;
newtree->nstride = fix->size_peratom_cols;
newtree->selfalloc = 0;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else error->all(FLERR,"Mismatched fix in variable formula");
@@ -1430,7 +1442,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = atof(var);
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = atof(var);
@@ -1458,7 +1471,8 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->array = reader[ivar]->fix->vstore;
newtree->nstride = 1;
newtree->selfalloc = 0;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
// v_name[N] = scalar from atom-style variable
@@ -1548,7 +1562,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
@@ -1568,7 +1583,8 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = VALUE;
newtree->value = value1;
- newtree->left = newtree->middle = newtree->right = NULL;
+ newtree->first = newtree->second = NULL;
+ newtree->nextra = 0;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
}
@@ -1645,12 +1661,13 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = opprevious;
if (opprevious == UNARY) {
- newtree->left = treestack[--ntreestack];
- newtree->middle = newtree->right = NULL;
+ newtree->first = treestack[--ntreestack];
+ newtree->second = NULL;
+ newtree->nextra = 0;
} else {
- newtree->right = treestack[--ntreestack];
- newtree->middle = NULL;
- newtree->left = treestack[--ntreestack];
+ newtree->second = treestack[--ntreestack];
+ newtree->first = treestack[--ntreestack];
+ newtree->nextra = 0;
}
treestack[ntreestack++] = newtree;
@@ -1762,36 +1779,36 @@ double Variable::collapse_tree(Tree *tree)
if (tree->type == BIGINTARRAY) return 0.0;
if (tree->type == ADD) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = arg1 + arg2;
return tree->value;
}
if (tree->type == SUBTRACT) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = arg1 - arg2;
return tree->value;
}
if (tree->type == MULTIPLY) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = arg1 * arg2;
return tree->value;
}
if (tree->type == DIVIDE) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg2 == 0.0) error->one(FLERR,"Divide by 0 in variable formula");
tree->value = arg1 / arg2;
@@ -1799,9 +1816,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == MODULO) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg2 == 0.0) error->one(FLERR,"Modulo 0 in variable formula");
tree->value = fmod(arg1,arg2);
@@ -1809,9 +1826,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == CARAT) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg2 == 0.0) error->one(FLERR,"Power by 0 in variable formula");
tree->value = pow(arg1,arg2);
@@ -1819,16 +1836,16 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == UNARY) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = -arg1;
return tree->value;
}
if (tree->type == NOT) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 == 0.0) tree->value = 1.0;
else tree->value = 0.0;
@@ -1836,9 +1853,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == EQ) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 == arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1846,9 +1863,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == NE) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 != arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1856,9 +1873,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == LT) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 < arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1866,9 +1883,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == LE) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 <= arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1876,9 +1893,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == GT) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 > arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1886,9 +1903,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == GE) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 >= arg2) tree->value = 1.0;
else tree->value = 0.0;
@@ -1896,9 +1913,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == AND) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 != 0.0 && arg2 != 0.0) tree->value = 1.0;
else tree->value = 0.0;
@@ -1906,9 +1923,9 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == OR) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 != 0.0 || arg2 != 0.0) tree->value = 1.0;
else tree->value = 0.0;
@@ -1916,8 +1933,8 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == SQRT) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 < 0.0)
error->one(FLERR,"Sqrt of negative value in variable formula");
@@ -1926,16 +1943,16 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == EXP) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = exp(arg1);
return tree->value;
}
if (tree->type == LN) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 <= 0.0)
error->one(FLERR,"Log of zero/negative value in variable formula");
@@ -1944,8 +1961,8 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == LOG) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 <= 0.0)
error->one(FLERR,"Log of zero/negative value in variable formula");
@@ -1954,40 +1971,40 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == ABS) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = fabs(arg1);
return tree->value;
}
if (tree->type == SIN) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = sin(arg1);
return tree->value;
}
if (tree->type == COS) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = cos(arg1);
return tree->value;
}
if (tree->type == TAN) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = tan(arg1);
return tree->value;
}
if (tree->type == ASIN) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 < -1.0 || arg1 > 1.0)
error->one(FLERR,"Arcsin of invalid value in variable formula");
@@ -1996,8 +2013,8 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == ACOS) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
if (arg1 < -1.0 || arg1 > 1.0)
error->one(FLERR,"Arccos of invalid value in variable formula");
@@ -2006,17 +2023,17 @@ double Variable::collapse_tree(Tree *tree)
}
if (tree->type == ATAN) {
- arg1 = collapse_tree(tree->left);
- if (tree->left->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ if (tree->first->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = atan(arg1);
return tree->value;
}
if (tree->type == ATAN2) {
- arg1 = collapse_tree(tree->left);
- arg2 = collapse_tree(tree->right);
- if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
+ arg1 = collapse_tree(tree->first);
+ arg2 = collapse_tree(tree->second);
+ if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
tree->type = VALUE;
tree->value = atan2(arg1,arg2);
return tree->value;
@@ -2025,10 +2042,10 @@ double Variable::collapse_tree(Tree *tree)
// random() or normal() do not become a single collapsed value
if (tree->type == RANDOM) {
- collapse_tree(tree->left);
- collapse_tree(tree->middle);
+ collapse_tree(tree->first);
+ collapse_tree(tree->second);
if (randomatom == NULL) {
- int seed = static_cast