Merge pull request lammps#4473 from akohlmey/collected-small-changes

Collected small changes and fixes

Author: akohlmey

doc/src/Commands_compute.rst

@@ -178,6 +178,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`ti <compute_ti>`
    * :doc:`torque/chunk <compute_torque_chunk>`
    * :doc:`vacf <compute_vacf>`
+   * :doc:`vacf/chunk <compute_vacf_chunk>`
    * :doc:`vcm/chunk <compute_vcm_chunk>`
    * :doc:`viscosity/cos <compute_viscosity_cos>`
    * :doc:`voronoi/atom <compute_voronoi_atom>`

doc/src/compute.rst

@@ -356,6 +356,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`ti <compute_ti>` - thermodynamic integration free energy values
 * :doc:`torque/chunk <compute_torque_chunk>` - torque applied on each chunk
 * :doc:`vacf <compute_vacf>` - velocity auto-correlation function of group of atoms
+* :doc:`vacf/chunk <compute_vacf_chunk>` - velocity auto-correlation for the center of mass velocities of chunks of atoms
 * :doc:`vcm/chunk <compute_vcm_chunk>` - velocity of center-of-mass for each chunk
 * :doc:`viscosity/cos <compute_viscosity_cos>` - velocity profile under cosine-shaped acceleration
 * :doc:`voronoi/atom <compute_voronoi_atom>` - Voronoi volume and neighbors for each atom

doc/src/compute_msd.rst

@@ -116,7 +116,9 @@ Compute *msd* cannot be used with a dynamic group.
 Related commands
 """"""""""""""""
 
-:doc:`compute msd/nongauss <compute_msd_nongauss>`, :doc:`compute displace_atom <compute_displace_atom>`, :doc:`fix store/state <fix_store_state>`, :doc:`compute msd/chunk <compute_msd_chunk>`
+:doc:`compute msd/nongauss <compute_msd_nongauss>`,
+:doc:`compute displace_atom <compute_displace_atom>`, :doc:`fix store/state <fix_store_state>`,
+:doc:`compute msd/chunk <compute_msd_chunk>`
 
 Default
 """""""

doc/src/compute_msd_chunk.rst

@@ -131,7 +131,7 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`compute msd <compute_msd>`
+:doc:`compute msd <compute_msd>`, :doc:`compute vacf/chunk <compute_vacf_chunk>`
 
 Default
 """""""

doc/src/compute_vacf.rst

@@ -76,7 +76,7 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`compute msd <compute_msd>`
+:doc:`compute msd <compute_msd>`, :doc:`compute vacf/chunk <compute_vacf_chunk>`
 
 Default
 """""""

doc/src/compute_vacf_chunk.rst

@@ -0,0 +1,124 @@
+.. index:: compute vacf/chunk
+
+compute vacf/chunk command
+==========================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID vacf/chunk chunkID
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* vacf/chunk = style name of this compute command
+* chunkID = ID of :doc:`compute chunk/atom <compute_chunk_atom>` command
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all vacf/chunk molchunk
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Define a computation that calculates the velocity auto-correlation
+function (VACF) for multiple chunks of atoms.
+
+In LAMMPS, chunks are collections of atoms defined by a :doc:`compute
+chunk/atom <compute_chunk_atom>` command, which assigns each atom to a
+single chunk (or no chunk).  The ID for this command is specified as
+chunkID.  For example, a single chunk could be the atoms in a molecule
+or atoms in a spatial bin.  See the :doc:`compute chunk/atom
+<compute_chunk_atom>` and :doc:`Howto chunk <Howto_chunk>` doc pages for
+details of how chunks can be defined and examples of how they can be
+used to measure properties of a system.
+
+Four quantities are calculated by this compute for each chunk.  The
+first 3 quantities are the product of the initial center of mass
+velocity (VCM) for each chunk in *x*, *y*, and *z* direction with the
+current center of mass velocity in the same direction.  The fourth
+component is the total VACF, i.e. the sum of the three components.
+
+Note that only atoms in the specified group contribute to the
+calculation.  The :doc:`compute chunk/atom <compute_chunk_atom>` command
+defines its own group; atoms will have a chunk ID = 0 if they are not in
+that group, signifying they are not assigned to a chunk, and will thus
+also not contribute to this calculation.  You can specify the "all"
+group for this command if you simply want to include atoms with non-zero
+chunk IDs.
+
+The integral of the VACF versus time is proportional to the diffusion
+coefficient of the diffusing chunks.
+
+.. note::
+
+   The number of chunks *Nchunk* calculated by the
+   :doc:`compute chunk/atom <compute_chunk_atom>` command must remain constant
+   each time this compute is invoked, so that the dot product for each chunk
+   from its original position can be computed consistently.  If *Nchunk*
+   does not remain constant, an error will be generated.  If needed, you
+   can enforce a constant *Nchunk* by using the *nchunk once* or *ids once*
+   options when specifying the :doc:`compute chunk/atom <compute_chunk_atom>`
+   command.
+
+.. note::
+
+   This compute stores the original center-of-mass velocities of each
+   chunk.  When a VACF is calculated on a later timestep, it is assumed
+   that the same atoms are assigned to the same chunk ID.  However
+   LAMMPS has no simple way to ensure this is the case, though you can
+   use the *ids once* option when specifying the :doc:`compute
+   chunk/atom <compute_chunk_atom>` command.  Note that if this is not
+   the case, the VACF calculation does not have a sensible meaning.
+
+.. note::
+
+   If you want the quantities calculated by this compute to be
+   continuous when running from a :doc:`restart file <read_restart>`, then
+   you should use the same ID for this compute, as in the original run.
+   This is so that the fix this compute creates to store per-chunk
+   quantities will also have the same ID, and thus be initialized
+   correctly with chunk reference positions from the restart file.
+
+The simplest way to output the results of the compute vacf/chunk
+calculation to a file is to use the :doc:`fix ave/time <fix_ave_time>`
+command, for example:
+
+.. code-block:: LAMMPS
+
+   compute cc1 all chunk/atom molecule
+   compute myChunk all vacf/chunk cc1
+   fix 1 all ave/time 100 1 100 c_myChunk[*] file tmp.out mode vector
+
+Output info
+"""""""""""
+
+This compute calculates a global array where the number of rows = the
+number of chunks *Nchunk* as calculated by the specified :doc:`compute
+chunk/atom <compute_chunk_atom>` command.  The number of columns = 4 for
+the *x*, *y*, *z*, component and the total VACF.  These values can be
+accessed by any command that uses global array values from a compute as
+input.  See the :doc:`Howto output <Howto_output>` page for an overview
+of LAMMPS output options.
+
+The array values are "intensive".  The array values will be in
+distance\ :math:`^2` divided by time\ :math:`^2` :doc:`units <units>`.
+
+Restrictions
+""""""""""""
+ none
+
+Related commands
+""""""""""""""""
+
+:doc:`compute vacf <compute_vacf>`, :doc:`compute msd/chunk <compute_msd_chunk>`
+
+Default
+"""""""
+
+none

python/lammps/core.py

@@ -17,7 +17,6 @@
 from __future__ import print_function
 
 import os
-import sys
 from ctypes import CDLL, POINTER, RTLD_GLOBAL, CFUNCTYPE, py_object, byref, cast, sizeof, \
   create_string_buffer, c_int, c_int32, c_int64, c_double, c_void_p, c_char_p, c_char,    \
   pythonapi

src/BOCS/compute_pressure_bocs.cpp

@@ -43,7 +43,7 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
   vptr(nullptr), id_temp(nullptr)
 {
   if (narg < 4) utils::missing_cmd_args(FLERR,"compute pressure/bocs", error);
-  if (igroup) error->all(FLERR,"Compute pressure/bocs must use group all");
+  if (igroup) error->all(FLERR, 1, "Compute pressure/bocs must use group all");
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
@@ -64,9 +64,9 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
 
     temperature = modify->get_compute_by_id(id_temp);
     if (!temperature)
-      error->all(FLERR,"Could not find compute pressure/bocs temperature compute {}", id_temp);
+      error->all(FLERR, 3, "Could not find compute pressure/bocs temperature compute {}", id_temp);
     if (temperature->tempflag == 0)
-      error->all(FLERR,"Compute pressure/bocs temperature compute {} does not compute "
+      error->all(FLERR, 3, "Compute pressure/bocs temperature compute {} does not compute "
                  "temperature", id_temp);
   }
 
@@ -104,8 +104,8 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
   // error check
 
   if (keflag && id_temp == nullptr)
-    error->all(FLERR,"Compute pressure/bocs requires temperature ID "
-               "to include kinetic energy");
+    error->all(FLERR, 3,
+               "Compute pressure/bocs requires temperature ID to include kinetic energy");
 
   vector = new double[size_vector];
   nvirial = 0;
@@ -119,9 +119,9 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
 
 ComputePressureBocs::~ComputePressureBocs()
 {
-  delete [] id_temp;
-  delete [] vector;
-  delete [] vptr;
+  delete[] id_temp;
+  delete[] vector;
+  delete[] vptr;
   if (phi_coeff) free(phi_coeff);
 }
 
@@ -139,7 +139,8 @@ void ComputePressureBocs::init()
   if (keflag) {
     temperature = modify->get_compute_by_id(id_temp);
     if (!temperature)
-      error->all(FLERR,"Could not find compute pressure/bocs temperature compute {}", id_temp);
+      error->all(FLERR, Error::NOLASTLINE,
+                 "Could not find compute pressure/bocs temperature compute {}", id_temp);
   }
 
   // detect contributions to virial
@@ -214,10 +215,10 @@ double ComputePressureBocs::find_index(double * grid, double value)
 
   if (value >= grid[i] && value <= (grid[i] + spacing)) { return i; }
 
-  error->all(FLERR,"find_index could not find value in grid for value: {}", value);
-  for (int i = 0; i < gridsize; ++i)
-  {
-    fprintf(stderr, "grid %d: %f\n",i,grid[i]);
+  error->all(FLERR, Error::NOLASTLINE,
+             "find_index could not find value in grid for value: {}", value);
+  for (int i = 0; i < gridsize; ++i) {
+    fprintf(stderr, "grid %d: %f\n", i, grid[i]);
   }
 
   exit(1);
@@ -237,7 +238,7 @@ double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,
     return grid[1][i] + (deltax) * ( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );
   else if (basis_type == BASIS_CUBIC_SPLINE)
     return grid[1][i] + (grid[2][i] * deltax) + (grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));
-  else error->all(FLERR,"bad spline type passed to get_cg_p_corr()\n");
+  else error->all(FLERR, Error::NOLASTLINE, "bad spline type passed to get_cg_p_corr()");
   return 0.0;
 }
 
@@ -251,7 +252,7 @@ void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,
                                        double sent_vavg)
 {
   if (basis_type == BASIS_ANALYTIC) p_basis_type = BASIS_ANALYTIC;
-  else error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
+  else error->all(FLERR, Error::NOLASTLINE, "Incorrect basis type passed to ComputePressureBocs");
 
   p_match_flag = 1;
 
@@ -276,7 +277,7 @@ void ComputePressureBocs::send_cg_info(int basis_type,
   else if (basis_type == BASIS_CUBIC_SPLINE) { p_basis_type = BASIS_CUBIC_SPLINE; }
   else
   {
-    error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
+    error->all(FLERR, Error::NOLASTLINE, "Incorrect basis type passed to ComputePressureBocs");
   }
   splines = in_splines;
   spline_length = gridsize;
@@ -293,7 +294,7 @@ double ComputePressureBocs::compute_scalar()
 {
   invoked_scalar = update->ntimestep;
   if (update->vflag_global != invoked_scalar)
-    error->all(FLERR,"Virial was not tallied on needed timestep");
+    error->all(FLERR, Error::NOLASTLINE, "Virial was not tallied on needed timestep");
 
   // invoke temperature if it hasn't been already
 
@@ -328,9 +329,8 @@ double ComputePressureBocs::compute_scalar()
       scalar = (virial[0] + virial[1] + virial[2]) / 3.0 *
                inv_volume * nktv2p + (correction);
   } else {
-    if (p_match_flag)
-    {
-      error->all(FLERR,"Pressure matching not implemented in 2-d.\n");
+    if (p_match_flag) {
+      error->all(FLERR, Error::NOLASTLINE, "Pressure matching not implemented in 2-d.");
       exit(1);
     } // The rest of this can probably be deleted.
     inv_volume = 1.0 / (domain->xprd * domain->yprd);
@@ -354,10 +354,10 @@ void ComputePressureBocs::compute_vector()
 {
   invoked_vector = update->ntimestep;
   if (update->vflag_global != invoked_vector)
-    error->all(FLERR,"Virial was not tallied on needed timestep");
+    error->all(FLERR, Error::NOLASTLINE, "Virial was not tallied on needed timestep");
 
   if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag)
-    error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' for "
+    error->all(FLERR, Error::NOLASTLINE, "Must use 'kspace_modify pressure/scalar no' for "
                "tensor components with kspace_style msm");
 
   // invoke temperature if it hasn't been already

src/EXTRA-COMMAND/region2vmd.cpp

@@ -765,85 +765,6 @@ void Region2VMD::write_region(FILE *fp, Region *region)
                    ellipsoid->c);
     }
 
-  } else if (regstyle == "plane") {
-    const auto plane = dynamic_cast<RegPlane *>(region);
-    if (!plane) {
-      error->one(FLERR, Error::NOLASTLINE, "Region {} is not of style 'plane'", region->id);
-    } else {
-      double v1[3], v2[3], v3[3], v4[3], v5[3], v6[3];
-      MathExtra::copy3(plane->normal, v1);
-      MathExtra::norm3(v1);
-      v2[0] = plane->xp + dx;
-      v2[1] = plane->yp + dy;
-      v2[2] = plane->zp + dz;
-      // determine the largest component of the plane normal and use that to determine the triangle edges
-      int dim1, dim2, dim3;
-      if (fabs(v1[0]) > fabs(v1[1])) {      // x > y
-        if (fabs(v1[0]) > fabs(v1[2])) {    // x > z => x is largest
-          dim3 = 0;
-          dim1 = 1;
-          dim2 = 2;
-        } else {    // z is largest
-          dim3 = 2;
-          dim1 = 0;
-          dim2 = 1;
-        }
-      } else {                              // x < y
-        if (fabs(v1[1]) > fabs(v1[2])) {    // y > z => y is largest
-          dim3 = 1;
-          dim1 = 0;
-          dim2 = 2;
-        } else {    // z is largest
-          dim3 = 2;
-          dim1 = 0;
-          dim2 = 1;
-        }
-      }
-
-      MathExtra::scale3(2.0, domain->boxlo, v3);
-      v3[dim3] = 0.0;
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v5);
-      v3[dim1] = domain->boxhi[dim1];
-      v3[dim2] = domain->boxlo[dim2];
-      v3[dim3] = 0.0;
-      MathExtra::scale3(2.0, v3);
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v6);
-      utils::print(fp, "draw triangle {{{} {} {}}} {{{} {} {}}} {{{} {} {}}}\n", v2[dim1], v2[dim2],
-                   v2[dim3], v5[dim1], v5[dim2], v5[dim3], v6[dim1], v6[dim2], v6[dim3]);
-      v3[dim1] = domain->boxlo[dim1];
-      v3[dim2] = domain->boxhi[dim2];
-      v3[dim3] = 0.0;
-      MathExtra::scale3(2.0, v3);
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v6);
-      utils::print(fp, "draw triangle {{{} {} {}}} {{{} {} {}}} {{{} {} {}}}\n", v2[dim1], v2[dim2],
-                   v2[dim3], v5[dim1], v5[dim2], v5[dim3], v6[dim1], v6[dim2], v6[dim3]);
-      v3[dim1] = domain->boxhi[dim1];
-      v3[dim2] = domain->boxhi[dim2];
-      v3[dim3] = 0.0;
-      MathExtra::scale3(2.0, v3);
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v5);
-      v3[dim1] = domain->boxhi[dim1];
-      v3[dim2] = domain->boxlo[dim2];
-      v3[dim3] = 0.0;
-      MathExtra::scale3(2.0, v3);
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v6);
-      utils::print(fp, "draw triangle {{{} {} {}}} {{{} {} {}}} {{{} {} {}}}\n", v2[dim1], v2[dim2],
-                   v2[dim3], v5[dim1], v5[dim2], v5[dim3], v6[dim1], v6[dim2], v6[dim3]);
-      v3[dim1] = domain->boxlo[dim1];
-      v3[dim2] = domain->boxhi[dim2];
-      v3[dim3] = 0.0;
-      MathExtra::scale3(2.0, v3);
-      MathExtra::cross3(v1, v3, v4);
-      MathExtra::add3(v4, v2, v6);
-      utils::print(fp, "draw triangle {{{} {} {}}} {{{} {} {}}} {{{} {} {}}}\n", v2[dim1], v2[dim2],
-                   v2[dim3], v5[dim1], v5[dim2], v5[dim3], v6[dim1], v6[dim2], v6[dim3]);
-    }
-
   } else if (regstyle == "prism") {
     const auto prism = dynamic_cast<RegPrism *>(region);
     if (!prism) {