Unverified Commit 92b508f1 authored by Richard Berger's avatar Richard Berger Committed by GitHub
Browse files

Merge pull request #1097 from lammps/localvars

add variable option to compute bond/angle/dihedral local
parents a33f45f1 bcecc038
Loading
Loading
Loading
Loading
+57 −9
Original line number Diff line number Diff line
@@ -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.

+57 −12
Original line number Diff line number Diff line
@@ -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.

+56 −8
Original line number Diff line number Diff line
@@ -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.

+153 −52
Original line number Diff line number Diff line
@@ -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;

  // style args

  nvalues = narg - 3;
  if (nvalues == 1) size_local_cols = 0;
  else size_local_cols = nvalues;
  bstyle = new int[nvalues];
  vstr = new char*[nvalues];
  vvar = new int[nvalues];
  
  tflag = eflag = -1;
  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;
  }

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

  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");
  // 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,8 +253,14 @@ int ComputeAngleLocal::compute_angles(int flag)
      if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
      if (atype == 0) continue;

      if (flag) {
        if (tflag >= 0) {
      if (!flag) {
        m++;
        continue;
      }

      // 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];
@@ -194,20 +278,37 @@ int ComputeAngleLocal::compute_angles(int flag)
	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;
          tbuf[n] = 180.0*acos(c)/MY_PI;
	theta = acos(c);
      }

      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++;
+6 −2
Original line number Diff line number Diff line
@@ -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;
Loading