diff --git a/doc/src/compute_angle_local.txt b/doc/src/compute_angle_local.txt
index 3a321965ef71135104bdf12f8dc13111dd515564..8acaec94d5331fb8ea8bd4898e4642333d77cf69 100644
--- a/doc/src/compute_angle_local.txt
+++ b/doc/src/compute_angle_local.txt
@@ -10,20 +10,27 @@ compute angle/local command :h3
[Syntax:]
-compute ID group-ID angle/local value1 value2 ... :pre
+compute ID group-ID angle/local value1 value2 ... keyword args ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
angle/local = style name of this compute command :l
one or more values may be appended :l
-value = {theta} or {eng} :l
+value = {theta} or {eng} or {v_name} :l
{theta} = tabulate angles
- {eng} = tabulate angle energies :pre
+ {eng} = tabulate angle energies
+ {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+ {set} args = theta name
+ theta = only currently allowed arg
+ name = name of variable to set with theta :pre
:ule
[Examples:]
compute 1 all angle/local theta
-compute 1 all angle/local eng theta :pre
+compute 1 all angle/local eng theta
+compute 1 all angle/local theta v_cos set theta t :pre
[Description:]
@@ -36,6 +43,47 @@ The value {theta} is the angle for the 3 atoms in the interaction.
The value {eng} is the interaction energy for the angle.
+The value {v_name} can be used together with the {set} keyword to
+compute a user-specified function of the angle theta. The {name}
+specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the angle theta. This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything. It must be an
+internal-style variable, because this command resets its value
+directly. The {set} keyword is used to identify the name of this
+other variable associated with theta.
+
+Note that the value of theta for each angle which stored in the
+internal variable is in radians, not degrees.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the cosine and cosine^2 of every angle in the system
+and output the statistics in various ways:
+
+variable t internal 0.0
+variable cos equal cos(v_t)
+variable cossq equal cos(v_t)*cos(v_t) :pre
+
+compute 1 all property/local aatom1 aatom2 aatom3 atype
+compute 2 all angle/local eng theta v_cos v_cossq set theta t
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the energy, angle,
+cosine(angle), cosine^2(angle) for every angle in the system. The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output. And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the cosine(angle) values and write them to a
+file.
+
+:line
+
The local data stored by this command is generated by looping over all
the atoms owned on a processor and their angles. An angle will only
be included if all 3 atoms in the angle are in the specified compute
@@ -65,12 +113,12 @@ dump 1 all local 1000 tmp.dump index c_1\[1\] c_1\[2\] c_1\[3\] c_1\[4\] c_2\[1\
[Output info:]
This compute calculates a local vector or local array depending on the
-number of keywords. The length of the vector or number of rows in the
-array is the number of angles. If a single keyword is specified, a
-local vector is produced. If two or more keywords are specified, a
+number of values. The length of the vector or number of rows in the
+array is the number of angles. If a single value is specified, a
+local vector is produced. If two or more values are specified, a
local array is produced where the number of columns = the number of
-keywords. The vector or array can be accessed by any command that
-uses local values from a compute as input. See the "Howto
+values. The vector or array can be accessed by any command that uses
+local values from a compute as input. See the "Howto
output"_Howto_output.html doc page for an overview of LAMMPS output
options.
diff --git a/doc/src/compute_bond_local.txt b/doc/src/compute_bond_local.txt
index c3dc1cc4afe92ae09118e8f56c0e678a380cee71..4afd1aec402e550e5e9f18d13cf597ba09b0104e 100644
--- a/doc/src/compute_bond_local.txt
+++ b/doc/src/compute_bond_local.txt
@@ -10,12 +10,12 @@ compute bond/local command :h3
[Syntax:]
-compute ID group-ID bond/local value1 value2 ... :pre
+compute ID group-ID bond/local value1 value2 ... keyword args ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
bond/local = style name of this compute command :l
one or more values may be appended :l
-value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {omega} or {velvib} :l
+value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {omega} or {velvib} or {v_name} :l
{dist} = bond distance
{engpot} = bond potential energy
{force} = bond force :pre
@@ -23,13 +23,22 @@ value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {
{engrot} = bond kinetic energy of rotation
{engtrans} = bond kinetic energy of translation
{omega} = magnitude of bond angular velocity
- {velvib} = vibrational velocity along the bond length :pre
+ {velvib} = vibrational velocity along the bond length
+ {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+ {set} args = dist name
+ dist = only currently allowed arg
+ name = name of variable to set with distance (dist) :pre
+:ule
+
:ule
[Examples:]
compute 1 all bond/local engpot
compute 1 all bond/local dist engpot force :pre
+compute 1 all angle/local dist v_distsq set dist d :pre
[Description:]
@@ -38,6 +47,10 @@ interactions. The number of datums generated, aggregated across all
processors, equals the number of bonds in the system, modified by the
group parameter as explained below.
+All these properties are computed for the pair of atoms in a bond,
+whether the 2 atoms represent a simple diatomic molecule, or are part
+of some larger molecule.
+
The value {dist} is the current length of the bond.
The value {engpot} is the potential energy for the bond,
@@ -79,9 +92,41 @@ two atoms in the bond towards each other. A negative value means the
2 atoms are moving toward each other; a positive value means they are
moving apart.
-Note that all these properties are computed for the pair of atoms in a
-bond, whether the 2 atoms represent a simple diatomic molecule, or are
-part of some larger molecule.
+The value {v_name} can be used together with the {set} keyword to
+compute a user-specified function of the bond distance. The {name}
+specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the bond distance. This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything. It must be an
+internal-style variable, because this command resets its value
+directly. The {set} keyword is used to identify the name of this
+other variable associated with theta.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the distance^2 of every bond in the system and
+output the statistics in various ways:
+
+variable d internal 0.0
+variable dsq equal v_d*v_d :pre
+
+compute 1 all property/local batom1 batom2 btype
+compute 2 all bond/local engpot dist v_dsq set dist d
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the energy, distance,
+distance^2 for every bond in the system. The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output. And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the distance^2 values and write them to a file.
+
+:line
The local data stored by this command is generated by looping over all
the atoms owned on a processor and their bonds. A bond will only be
@@ -111,12 +156,12 @@ dump 1 all local 1000 tmp.dump index c_1\[*\] c_2\[*\] :pre
[Output info:]
This compute calculates a local vector or local array depending on the
-number of keywords. The length of the vector or number of rows in the
-array is the number of bonds. If a single keyword is specified, a
-local vector is produced. If two or more keywords are specified, a
-local array is produced where the number of columns = the number of
-keywords. The vector or array can be accessed by any command that
-uses local values from a compute as input. See the "Howto
+number of values. The length of the vector or number of rows in the
+array is the number of bonds. If a single value is specified, a local
+vector is produced. If two or more values are specified, a local
+array is produced where the number of columns = the number of values.
+The vector or array can be accessed by any command that uses local
+values from a compute as input. See the "Howto
output"_Howto_output.html doc page for an overview of LAMMPS output
options.
diff --git a/doc/src/compute_dihedral_local.txt b/doc/src/compute_dihedral_local.txt
index 77812699d35fc826c95a344ddc7bf0207caf2ec4..951b360f19516e36433ed4a10a1966c83d10bcfb 100644
--- a/doc/src/compute_dihedral_local.txt
+++ b/doc/src/compute_dihedral_local.txt
@@ -10,18 +10,25 @@ compute dihedral/local command :h3
[Syntax:]
-compute ID group-ID dihedral/local value1 value2 ... :pre
+compute ID group-ID dihedral/local value1 value2 ... keyword args ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
dihedral/local = style name of this compute command :l
one or more values may be appended :l
-value = {phi} :l
- {phi} = tabulate dihedral angles :pre
+value = {phi} or {v_name} :l
+ {phi} = tabulate dihedral angles
+ {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+ {set} args = phi name
+ phi = only currently allowed arg
+ name = name of variable to set with phi :pre
:ule
[Examples:]
compute 1 all dihedral/local phi :pre
+compute 1 all dihedral/local phi v_cos set phi p :pre
[Description:]
@@ -33,6 +40,47 @@ by the group parameter as explained below.
The value {phi} is the dihedral angle, as defined in the diagram on
the "dihedral_style"_dihedral_style.html doc page.
+The value {v_name} can be used together with the {set} keyword to
+compute a user-specified function of the dihedral angle phi. The
+{name} specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the angle phi. This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything. It must be an
+internal-style variable, because this command resets its value
+directly. The {set} keyword is used to identify the name of this
+other variable associated with phi.
+
+Note that the value of phi for each angle which stored in the internal
+variable is in radians, not degrees.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the cosine and cosine^2 of every dihedral angle in
+the system and output the statistics in various ways:
+
+variable p internal 0.0
+variable cos equal cos(v_p)
+variable cossq equal cos(v_p)*cos(v_p) :pre
+
+compute 1 all property/local datom1 datom2 datom3 datom4 dtype
+compute 2 all dihedral/local phi v_cos v_cossq set phi p
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the angle,
+cosine(angle), cosine^2(angle) for every dihedral in the system. The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output. And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the cosine(angle) values and write them to a
+file.
+
+:line
+
The local data stored by this command is generated by looping over all
the atoms owned on a processor and their dihedrals. A dihedral will
only be included if all 4 atoms in the dihedral are in the specified
@@ -57,12 +105,12 @@ dump 1 all local 1000 tmp.dump index c_1\[1\] c_1\[2\] c_1\[3\] c_1\[4\] c_1\[5\
[Output info:]
This compute calculates a local vector or local array depending on the
-number of keywords. The length of the vector or number of rows in the
-array is the number of dihedrals. If a single keyword is specified, a
-local vector is produced. If two or more keywords are specified, a
+number of values. The length of the vector or number of rows in the
+array is the number of dihedrals. If a single value is specified, a
+local vector is produced. If two or more values are specified, a
local array is produced where the number of columns = the number of
-keywords. The vector or array can be accessed by any command that
-uses local values from a compute as input. See the "Howto
+values. The vector or array can be accessed by any command that uses
+local values from a compute as input. See the "Howto
output"_Howto_output.html doc page for an overview of LAMMPS output
options.
diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp
index 7137077158eb33e1dfc2abbea7378426639e0f65..95faad54ad6e898cc89fe24663a002fb563b5353 100644
--- a/src/compute_angle_local.cpp
+++ b/src/compute_angle_local.cpp
@@ -21,6 +21,8 @@
#include "domain.h"
#include "force.h"
#include "angle.h"
+#include "input.h"
+#include "variable.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
@@ -30,11 +32,13 @@ using namespace MathConst;
#define DELTA 10000
+enum{THETA,ENG,VARIABLE};
+
/* ---------------------------------------------------------------------- */
ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- vlocal(NULL), alocal(NULL)
+ bstyle(NULL), vstr(NULL), vvar(NULL), tstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute angle/local command");
@@ -42,19 +46,82 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Compute angle/local used when angles are not allowed");
local_flag = 1;
- nvalues = narg - 3;
- if (nvalues == 1) size_local_cols = 0;
- else size_local_cols = nvalues;
- tflag = eflag = -1;
+ // style args
+
+ nvalues = narg - 3;
+ bstyle = new int[nvalues];
+ vstr = new char*[nvalues];
+ vvar = new int[nvalues];
+
nvalues = 0;
+ tflag = 0;
+ nvar = 0;
+
+ int iarg;
+ for (iarg = 3; iarg < narg; iarg++) {
+ if (strcmp(arg[iarg],"theta") == 0) {
+ bstyle[nvalues++] = THETA;
+ tflag = 1;
+ } else if (strcmp(arg[iarg],"eng") == 0) {
+ bstyle[nvalues++] = ENG;
+ } else if (strncmp(arg[iarg],"v_",2) == 0) {
+ bstyle[nvalues++] = VARIABLE;
+ int n = strlen(arg[iarg]);
+ vstr[nvar] = new char[n];
+ strcpy(vstr[nvar],&arg[iarg][2]);
+ nvar++;
+ } else break;
+ }
- for (int iarg = 3; iarg < narg; iarg++) {
- if (strcmp(arg[iarg],"theta") == 0) tflag = nvalues++;
- else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
- else error->all(FLERR,"Invalid keyword in compute angle/local command");
+ // optional args
+
+ setflag = 0;
+ tstr = NULL;
+
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"set") == 0) {
+ setflag = 1;
+ if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command");
+ if (strcmp(arg[iarg+1],"theta") == 0) {
+ delete [] tstr;
+ int n = strlen(arg[iarg+2]) + 1;
+ tstr = new char[n];
+ strcpy(tstr,arg[iarg+2]);
+ tflag = 1;
+ } else error->all(FLERR,"Illegal compute angle/local command");
+ iarg += 3;
+ } else error->all(FLERR,"Illegal compute angle/local command");
}
+ // error check
+
+ if (nvar) {
+ if (!setflag)
+ error->all(FLERR,"Compute angle/local variable requires a set variable");
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,"Variable name for copute angle/local does not exist");
+ if (!input->variable->equalstyle(vvar[i]))
+ error->all(FLERR,"Variable for compute angle/local is invalid style");
+ }
+
+ if (tstr) {
+ tvar = input->variable->find(tstr);
+ if (tvar < 0)
+ error->all(FLERR,"Variable name for compute angle/local does not exist");
+ if (!input->variable->internalstyle(tvar))
+ error->all(FLERR,"Variable for compute angle/local is invalid style");
+ }
+ } else if (setflag)
+ error->all(FLERR,"Compute angle/local set with no variable");
+
+ // initialize output
+
+ if (nvalues == 1) size_local_cols = 0;
+ else size_local_cols = nvalues;
+
nmax = 0;
vlocal = NULL;
alocal = NULL;
@@ -64,6 +131,13 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeAngleLocal::~ComputeAngleLocal()
{
+ delete [] bstyle;
+ for (int i = 0; i < nvar; i++) delete [] vstr[i];
+ delete [] vstr;
+ delete [] vvar;
+
+ delete [] tstr;
+
memory->destroy(vlocal);
memory->destroy(alocal);
}
@@ -75,6 +149,20 @@ void ComputeAngleLocal::init()
if (force->angle == NULL)
error->all(FLERR,"No angle style is defined for compute angle/local");
+ if (nvar) {
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,"Variable name for compute angle/local does not exist");
+ }
+
+ if (tstr) {
+ tvar = input->variable->find(tstr);
+ if (tvar < 0)
+ error->all(FLERR,"Variable name for compute angle/local does not exist");
+ }
+ }
+
// do initial memory allocation so that memory_usage() is correct
ncount = compute_angles(0);
@@ -109,11 +197,11 @@ void ComputeAngleLocal::compute_local()
int ComputeAngleLocal::compute_angles(int flag)
{
- int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype;
+ int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar;
tagint tagprev;
double delx1,dely1,delz1,delx2,dely2,delz2;
- double rsq1,rsq2,r1,r2,c;
- double *tbuf,*ebuf;
+ double rsq1,rsq2,r1,r2,c,theta;
+ double *ptr;
double **x = atom->x;
tagint *tag = atom->tag;
@@ -131,17 +219,7 @@ int ComputeAngleLocal::compute_angles(int flag)
int nlocal = atom->nlocal;
int molecular = atom->molecular;
- if (flag) {
- if (nvalues == 1) {
- if (tflag >= 0) tbuf = vlocal;
- if (eflag >= 0) ebuf = vlocal;
- } else {
- if (tflag >= 0 && alocal) tbuf = &alocal[0][tflag];
- else tbuf = NULL;
- if (eflag >= 0 && alocal) ebuf = &alocal[0][eflag];
- else ebuf = NULL;
- }
- }
+ // loop over all atoms and their angles
Angle *angle = force->angle;
@@ -175,39 +253,62 @@ int ComputeAngleLocal::compute_angles(int flag)
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atype == 0) continue;
- if (flag) {
- if (tflag >= 0) {
- delx1 = x[atom1][0] - x[atom2][0];
- dely1 = x[atom1][1] - x[atom2][1];
- delz1 = x[atom1][2] - x[atom2][2];
- domain->minimum_image(delx1,dely1,delz1);
-
- rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
- r1 = sqrt(rsq1);
-
- delx2 = x[atom3][0] - x[atom2][0];
- dely2 = x[atom3][1] - x[atom2][1];
- delz2 = x[atom3][2] - x[atom2][2];
- domain->minimum_image(delx2,dely2,delz2);
-
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- r2 = sqrt(rsq2);
+ if (!flag) {
+ m++;
+ continue;
+ }
- // c = cosine of angle
+ // theta needed by one or more outputs
+
+ if (tflag) {
+ delx1 = x[atom1][0] - x[atom2][0];
+ dely1 = x[atom1][1] - x[atom2][1];
+ delz1 = x[atom1][2] - x[atom2][2];
+ domain->minimum_image(delx1,dely1,delz1);
+
+ rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
+ r1 = sqrt(rsq1);
+
+ delx2 = x[atom3][0] - x[atom2][0];
+ dely2 = x[atom3][1] - x[atom2][1];
+ delz2 = x[atom3][2] - x[atom2][2];
+ domain->minimum_image(delx2,dely2,delz2);
+
+ rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ r2 = sqrt(rsq2);
+
+ // c = cosine of angle
+ // theta = angle in radians
+
+ c = delx1*delx2 + dely1*dely2 + delz1*delz2;
+ c /= r1*r2;
+ if (c > 1.0) c = 1.0;
+ if (c < -1.0) c = -1.0;
+ theta = acos(c);
+ }
- c = delx1*delx2 + dely1*dely2 + delz1*delz2;
- c /= r1*r2;
- if (c > 1.0) c = 1.0;
- if (c < -1.0) c = -1.0;
- tbuf[n] = 180.0*acos(c)/MY_PI;
- }
+ if (nvalues == 1) ptr = &vlocal[m];
+ else ptr = alocal[m];
+
+ if (nvar) {
+ ivar = 0;
+ if (tstr) input->variable->internal_set(tvar,theta);
+ }
- if (eflag >= 0) {
- if (atype > 0)
- ebuf[n] = angle->single(atype,atom1,atom2,atom3);
- else ebuf[n] = 0.0;
+ for (n = 0; n < nvalues; n++) {
+ switch (bstyle[n]) {
+ case THETA:
+ ptr[n] = 180.0*theta/MY_PI;
+ break;
+ case ENG:
+ if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3);
+ else ptr[n] = 0.0;
+ break;
+ case VARIABLE:
+ ptr[n] = input->variable->compute_equal(vvar[ivar]);
+ ivar++;
+ break;
}
- n += nvalues;
}
m++;
diff --git a/src/compute_angle_local.h b/src/compute_angle_local.h
index 6f59d7f097c74d4db8b47ea59dcc404f40e846ad..a95e23206d18d4244a51c804a5f3dcfffe3e3381 100644
--- a/src/compute_angle_local.h
+++ b/src/compute_angle_local.h
@@ -33,8 +33,12 @@ class ComputeAngleLocal : public Compute {
double memory_usage();
private:
- int nvalues,tflag,eflag;
- int ncount;
+ int nvalues,nvar,ncount,setflag,tflag;
+
+ int tvar;
+ int *bstyle,*vvar;
+ char *tstr;
+ char **vstr;
int nmax;
double *vlocal;
diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp
index 56abb9b07113722f9fdf39e6cc6368c89334e352..6a179cf1b4484e3949339c21d139bff52f33f7aa 100644
--- a/src/compute_bond_local.cpp
+++ b/src/compute_bond_local.cpp
@@ -21,8 +21,10 @@
#include "domain.h"
#include "force.h"
#include "bond.h"
-#include "math_extra.h"
#include "comm.h"
+#include "input.h"
+#include "variable.h"
+#include "math_extra.h"
#include "memory.h"
#include "error.h"
@@ -31,13 +33,13 @@ using namespace LAMMPS_NS;
#define DELTA 10000
#define EPSILON 1.0e-12
-enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE};
+enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,VARIABLE};
/* ---------------------------------------------------------------------- */
ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- bstyle(NULL), vlocal(NULL), alocal(NULL)
+ bstyle(NULL), vstr(NULL), vvar(NULL), dstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute bond/local command");
@@ -47,14 +49,18 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
local_flag = 1;
comm_forward = 3;
- nvalues = narg - 3;
- if (nvalues == 1) size_local_cols = 0;
- else size_local_cols = nvalues;
+ // style args
+ nvalues = narg - 3;
bstyle = new int[nvalues];
-
+ vstr = new char*[nvalues];
+ vvar = new int[nvalues];
+
nvalues = 0;
- for (int iarg = 3; iarg < narg; iarg++) {
+ nvar = 0;
+
+ int iarg;
+ for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
@@ -63,9 +69,58 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS;
else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA;
else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB;
- else error->all(FLERR,"Invalid keyword in compute bond/local command");
+ else if (strncmp(arg[iarg],"v_",2) == 0) {
+ bstyle[nvalues++] = VARIABLE;
+ int n = strlen(arg[iarg]);
+ vstr[nvar] = new char[n];
+ strcpy(vstr[nvar],&arg[iarg][2]);
+ nvar++;
+ } else break;
+ }
+
+ // optional args
+
+ setflag = 0;
+ dstr = NULL;
+
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"set") == 0) {
+ setflag = 1;
+ if (iarg+3 > narg) error->all(FLERR,"Illegal compute bond/local command");
+ if (strcmp(arg[iarg+1],"dist") == 0) {
+ delete [] dstr;
+ int n = strlen(arg[iarg+2]) + 1;
+ dstr = new char[n];
+ strcpy(dstr,arg[iarg+2]);
+ } else error->all(FLERR,"Illegal compute bond/local command");
+ iarg += 3;
+ } else error->all(FLERR,"Illegal compute bond/local command");
}
+ // error check
+
+ if (nvar) {
+ if (!setflag)
+ error->all(FLERR,"Compute bond/local variable requires a set variable");
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,"Variable name for copute bond/local does not exist");
+ if (!input->variable->equalstyle(vvar[i]))
+ error->all(FLERR,"Variable for compute bond/local is invalid style");
+ }
+
+ if (dstr) {
+ dvar = input->variable->find(dstr);
+ if (dvar < 0)
+ error->all(FLERR,"Variable name for compute bond/local does not exist");
+ if (!input->variable->internalstyle(dvar))
+ error->all(FLERR,"Variable for compute bond/local is invalid style");
+ }
+ } else if (setflag)
+ error->all(FLERR,"Compute bond/local set with no variable");
+
+
// set singleflag if need to call bond->single()
// set velflag if compute any quantities based on velocities
@@ -77,6 +132,11 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
}
+ // initialize output
+
+ if (nvalues == 1) size_local_cols = 0;
+ else size_local_cols = nvalues;
+
nmax = 0;
vlocal = NULL;
alocal = NULL;
@@ -86,9 +146,15 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeBondLocal::~ComputeBondLocal()
{
+ delete [] bstyle;
+ for (int i = 0; i < nvar; i++) delete [] vstr[i];
+ delete [] vstr;
+ delete [] vvar;
+
+ delete [] dstr;
+
memory->destroy(vlocal);
memory->destroy(alocal);
- delete [] bstyle;
}
/* ---------------------------------------------------------------------- */
@@ -98,6 +164,20 @@ void ComputeBondLocal::init()
if (force->bond == NULL)
error->all(FLERR,"No bond style is defined for compute bond/local");
+ if (nvar) {
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,"Variable name for compute bond/local does not exist");
+ }
+
+ if (dstr) {
+ dvar = input->variable->find(dstr);
+ if (dvar < 0)
+ error->all(FLERR,"Variable name for compute bond/local does not exist");
+ }
+ }
+
// set ghostvelflag if need to acquire ghost atom velocities
if (velflag && !comm->ghost_velocity) ghostvelflag = 1;
@@ -140,7 +220,7 @@ void ComputeBondLocal::compute_local()
int ComputeBondLocal::compute_bonds(int flag)
{
- int i,m,n,nb,atom1,atom2,imol,iatom,btype;
+ int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar;
tagint tagprev;
double dx,dy,dz,rsq;
double mass1,mass2,masstotal,invmasstotal;
@@ -297,6 +377,11 @@ int ComputeBondLocal::compute_bonds(int flag)
if (nvalues == 1) ptr = &vlocal[m];
else ptr = alocal[m];
+ if (nvar) {
+ ivar = 0;
+ if (dstr) input->variable->internal_set(dvar,sqrt(rsq));
+ }
+
for (n = 0; n < nvalues; n++) {
switch (bstyle[n]) {
case DIST:
@@ -323,6 +408,10 @@ int ComputeBondLocal::compute_bonds(int flag)
case VELVIB:
ptr[n] = vvib;
break;
+ case VARIABLE:
+ ptr[n] = input->variable->compute_equal(vvar[ivar]);
+ ivar++;
+ break;
}
}
}
diff --git a/src/compute_bond_local.h b/src/compute_bond_local.h
index b3e99fc61b9ac12e4c2bd6ce97d4355d5ddeedfd..17111a941c36d8b7020e965e385ad05db9d55316 100644
--- a/src/compute_bond_local.h
+++ b/src/compute_bond_local.h
@@ -35,10 +35,13 @@ class ComputeBondLocal : public Compute {
double memory_usage();
private:
- int nvalues;
- int ncount;
- int *bstyle;
+ int nvalues,nvar,ncount,setflag;
+
int singleflag,velflag,ghostvelflag,initflag;
+ int dvar;
+ int *bstyle,*vvar;
+ char *dstr;
+ char **vstr;
int nmax;
double *vlocal;
diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp
index 42d1476ad25fb368b08b36bf99195f3a0750235b..7444630090f231aadf4aabdb2c6fa66e00c2e989 100644
--- a/src/compute_dihedral_local.cpp
+++ b/src/compute_dihedral_local.cpp
@@ -21,6 +21,9 @@
#include "domain.h"
#include "force.h"
#include "dihedral.h"
+#include "input.h"
+#include "variable.h"
+
#include "math_const.h"
#include "memory.h"
#include "error.h"
@@ -31,11 +34,13 @@ using namespace MathConst;
#define DELTA 10000
#define SMALL 0.001
+enum{PHI,VARIABLE};
+
/* ---------------------------------------------------------------------- */
ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- vlocal(NULL), alocal(NULL)
+ bstyle(NULL), vstr(NULL), vvar(NULL), pstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command");
@@ -44,18 +49,80 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
"Compute dihedral/local used when dihedrals are not allowed");
local_flag = 1;
- nvalues = narg - 3;
- if (nvalues == 1) size_local_cols = 0;
- else size_local_cols = nvalues;
- pflag = -1;
+ // style args
+
+ nvalues = narg - 3;
+ bstyle = new int[nvalues];
+ vstr = new char*[nvalues];
+ vvar = new int[nvalues];
+
nvalues = 0;
+ nvar = 0;
+
+ int iarg;
+ for (iarg = 3; iarg < narg; iarg++) {
+ if (strcmp(arg[iarg],"phi") == 0) {
+ bstyle[nvalues++] = PHI;
+ } else if (strncmp(arg[iarg],"v_",2) == 0) {
+ bstyle[nvalues++] = VARIABLE;
+ int n = strlen(arg[iarg]);
+ vstr[nvar] = new char[n];
+ strcpy(vstr[nvar],&arg[iarg][2]);
+ nvar++;
+ } else break;
+ }
- for (int iarg = 3; iarg < narg; iarg++) {
- if (strcmp(arg[iarg],"phi") == 0) pflag = nvalues++;
- else error->all(FLERR,"Invalid keyword in compute dihedral/local command");
+ // optional args
+
+ setflag = 0;
+ pstr = NULL;
+
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"set") == 0) {
+ setflag = 1;
+ if (iarg+3 > narg)
+ error->all(FLERR,"Illegal compute dihedral/local command");
+ if (strcmp(arg[iarg+1],"phi") == 0) {
+ delete [] pstr;
+ int n = strlen(arg[iarg+2]) + 1;
+ pstr = new char[n];
+ strcpy(pstr,arg[iarg+2]);
+ } else error->all(FLERR,"Illegal compute dihedral/local command");
+ iarg += 3;
+ } else error->all(FLERR,"Illegal compute dihedral/local command");
}
+ // error check
+
+ if (nvar) {
+ if (!setflag)
+ error->all(FLERR,"Compute dihedral/local variable requires a set variable");
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,
+ "Variable name for copute dihedral/local does not exist");
+ if (!input->variable->equalstyle(vvar[i]))
+ error->all(FLERR,"Variable for compute dihedral/local is invalid style");
+ }
+
+ if (pstr) {
+ pvar = input->variable->find(pstr);
+ if (pvar < 0)
+ error->all(FLERR,
+ "Variable name for compute dihedral/local does not exist");
+ if (!input->variable->internalstyle(pvar))
+ error->all(FLERR,"Variable for compute dihedral/local is invalid style");
+ }
+ } else if (setflag)
+ error->all(FLERR,"Compute dihedral/local set with no variable");
+
+ // initialize output
+
+ if (nvalues == 1) size_local_cols = 0;
+ else size_local_cols = nvalues;
+
nmax = 0;
vlocal = NULL;
alocal = NULL;
@@ -65,6 +132,13 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeDihedralLocal::~ComputeDihedralLocal()
{
+ delete [] bstyle;
+ for (int i = 0; i < nvar; i++) delete [] vstr[i];
+ delete [] vstr;
+ delete [] vvar;
+
+ delete [] pstr;
+
memory->destroy(vlocal);
memory->destroy(alocal);
}
@@ -76,6 +150,22 @@ void ComputeDihedralLocal::init()
if (force->dihedral == NULL)
error->all(FLERR,"No dihedral style is defined for compute dihedral/local");
+ if (nvar) {
+ for (int i = 0; i < nvar; i++) {
+ vvar[i] = input->variable->find(vstr[i]);
+ if (vvar[i] < 0)
+ error->all(FLERR,
+ "Variable name for compute dihedral/local does not exist");
+ }
+
+ if (pstr) {
+ pvar = input->variable->find(pstr);
+ if (pvar < 0)
+ error->all(FLERR,
+ "Variable name for compute dihedral/local does not exist");
+ }
+ }
+
// do initial memory allocation so that memory_usage() is correct
ncount = compute_dihedrals(0);
@@ -107,12 +197,12 @@ void ComputeDihedralLocal::compute_local()
int ComputeDihedralLocal::compute_dihedrals(int flag)
{
- int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom;
+ int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar;
tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv;
- double s,c;
- double *pbuf;
+ double s,c,phi;
+ double *ptr;
double **x = atom->x;
tagint *tag = atom->tag;
@@ -130,14 +220,7 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
int nlocal = atom->nlocal;
int molecular = atom->molecular;
- if (flag) {
- if (nvalues == 1) {
- if (pflag >= 0) pbuf = vlocal;
- } else {
- if (pflag >= 0 && alocal) pbuf = &alocal[0][pflag];
- else pbuf = NULL;
- }
- }
+ // loop over all atoms and their dihedrals
m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) {
@@ -169,56 +252,75 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
- if (flag) {
-
- // phi calculation from dihedral style harmonic
-
- if (pflag >= 0) {
- vb1x = x[atom1][0] - x[atom2][0];
- vb1y = x[atom1][1] - x[atom2][1];
- vb1z = x[atom1][2] - x[atom2][2];
- domain->minimum_image(vb1x,vb1y,vb1z);
-
- vb2x = x[atom3][0] - x[atom2][0];
- vb2y = x[atom3][1] - x[atom2][1];
- vb2z = x[atom3][2] - x[atom2][2];
- domain->minimum_image(vb2x,vb2y,vb2z);
-
- vb2xm = -vb2x;
- vb2ym = -vb2y;
- vb2zm = -vb2z;
- domain->minimum_image(vb2xm,vb2ym,vb2zm);
-
- vb3x = x[atom4][0] - x[atom3][0];
- vb3y = x[atom4][1] - x[atom3][1];
- vb3z = x[atom4][2] - x[atom3][2];
- domain->minimum_image(vb3x,vb3y,vb3z);
-
- ax = vb1y*vb2zm - vb1z*vb2ym;
- ay = vb1z*vb2xm - vb1x*vb2zm;
- az = vb1x*vb2ym - vb1y*vb2xm;
- bx = vb3y*vb2zm - vb3z*vb2ym;
- by = vb3z*vb2xm - vb3x*vb2zm;
- bz = vb3x*vb2ym - vb3y*vb2xm;
-
- rasq = ax*ax + ay*ay + az*az;
- rbsq = bx*bx + by*by + bz*bz;
- rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
- rg = sqrt(rgsq);
-
- ra2inv = rb2inv = 0.0;
- if (rasq > 0) ra2inv = 1.0/rasq;
- if (rbsq > 0) rb2inv = 1.0/rbsq;
- rabinv = sqrt(ra2inv*rb2inv);
-
- c = (ax*bx + ay*by + az*bz)*rabinv;
- s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
-
- if (c > 1.0) c = 1.0;
- if (c < -1.0) c = -1.0;
- pbuf[n] = 180.0*atan2(s,c)/MY_PI;
+ if (!flag) {
+ m++;
+ continue;
+ }
+
+ // phi calculation from dihedral style harmonic
+
+ vb1x = x[atom1][0] - x[atom2][0];
+ vb1y = x[atom1][1] - x[atom2][1];
+ vb1z = x[atom1][2] - x[atom2][2];
+ domain->minimum_image(vb1x,vb1y,vb1z);
+
+ vb2x = x[atom3][0] - x[atom2][0];
+ vb2y = x[atom3][1] - x[atom2][1];
+ vb2z = x[atom3][2] - x[atom2][2];
+ domain->minimum_image(vb2x,vb2y,vb2z);
+
+ vb2xm = -vb2x;
+ vb2ym = -vb2y;
+ vb2zm = -vb2z;
+ domain->minimum_image(vb2xm,vb2ym,vb2zm);
+
+ vb3x = x[atom4][0] - x[atom3][0];
+ vb3y = x[atom4][1] - x[atom3][1];
+ vb3z = x[atom4][2] - x[atom3][2];
+ domain->minimum_image(vb3x,vb3y,vb3z);
+
+ ax = vb1y*vb2zm - vb1z*vb2ym;
+ ay = vb1z*vb2xm - vb1x*vb2zm;
+ az = vb1x*vb2ym - vb1y*vb2xm;
+ bx = vb3y*vb2zm - vb3z*vb2ym;
+ by = vb3z*vb2xm - vb3x*vb2zm;
+ bz = vb3x*vb2ym - vb3y*vb2xm;
+
+ rasq = ax*ax + ay*ay + az*az;
+ rbsq = bx*bx + by*by + bz*bz;
+ rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
+ rg = sqrt(rgsq);
+
+ ra2inv = rb2inv = 0.0;
+ if (rasq > 0) ra2inv = 1.0/rasq;
+ if (rbsq > 0) rb2inv = 1.0/rbsq;
+ rabinv = sqrt(ra2inv*rb2inv);
+
+ c = (ax*bx + ay*by + az*bz)*rabinv;
+ s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
+
+ if (c > 1.0) c = 1.0;
+ if (c < -1.0) c = -1.0;
+ phi = atan2(s,c);
+
+ if (nvalues == 1) ptr = &vlocal[m];
+ else ptr = alocal[m];
+
+ if (nvar) {
+ ivar = 0;
+ if (pstr) input->variable->internal_set(pvar,phi);
+ }
+
+ for (n = 0; n < nvalues; n++) {
+ switch (bstyle[n]) {
+ case PHI:
+ ptr[n] = 180.0*phi/MY_PI;
+ break;
+ case VARIABLE:
+ ptr[n] = input->variable->compute_equal(vvar[ivar]);
+ ivar++;
+ break;
}
- n += nvalues;
}
m++;
diff --git a/src/compute_dihedral_local.h b/src/compute_dihedral_local.h
index c077e52dbd3470b75672d70773a0706d0e3b8498..d5f6a641833435a4f97de69f650df94e357aae63 100644
--- a/src/compute_dihedral_local.h
+++ b/src/compute_dihedral_local.h
@@ -33,8 +33,12 @@ class ComputeDihedralLocal : public Compute {
double memory_usage();
private:
- int nvalues,pflag;
- int ncount;
+ int nvalues,nvar,ncount,setflag,tflag;
+
+ int pvar;
+ int *bstyle,*vvar;
+ char *pstr;
+ char **vstr;
int nmax;
double *vlocal;