diff --git a/doc/src/Commands_compute.txt b/doc/src/Commands_compute.txt index e035eb8431..5c6eba4ddf 100644 --- a/doc/src/Commands_compute.txt +++ b/doc/src/Commands_compute.txt @@ -67,6 +67,7 @@ KOKKOS, o = USER-OMP, t = OPT. "gyration"_compute_gyration.html, "gyration/chunk"_compute_gyration_chunk.html, "gyration/shape"_compute_gyration_shape.html, +"gyration/shape/chunk"_compute_gyration_shape_chunk.html, "heat/flux"_compute_heat_flux.html, "heat/flux/tally"_compute_tally.html, "hexorder/atom"_compute_hexorder_atom.html, diff --git a/doc/src/compute.txt b/doc/src/compute.txt index b54d2d2e7b..2d6d26ff16 100644 --- a/doc/src/compute.txt +++ b/doc/src/compute.txt @@ -213,7 +213,8 @@ compute"_Commands_compute.html doc page are followed by one or more of "group/group"_compute_group_group.html - energy/force between two groups of atoms "gyration"_compute_gyration.html - radius of gyration of group of atoms "gyration/chunk"_compute_gyration_chunk.html - radius of gyration for each chunk -"gyration/shape"_compute_gyration_shape.html - compute shape parameters from radius of gyration tensor +"gyration/shape"_compute_gyration_shape.html - shape parameters from gyration tensor +"gyration/shape/chunk"_compute_gyration_shape_chunk.html - shape parameters from gyration tensor for each chunk "heat/flux"_compute_heat_flux.html - heat flux through a group of atoms "heat/flux/tally"_compute_tally.html - "hexorder/atom"_compute_hexorder_atom.html - bond orientational order parameter q6 diff --git a/doc/src/compute_gyration_shape_chunk.txt b/doc/src/compute_gyration_shape_chunk.txt new file mode 100644 index 0000000000..bc34594905 --- /dev/null +++ b/doc/src/compute_gyration_shape_chunk.txt @@ -0,0 +1,89 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +compute gyration/shape/chunk command :h3 + +[Syntax:] + +compute ID group-ID gyration/shape/chunk compute-ID :pre + +ID, group-ID are documented in "compute"_compute.html command +gyration/shape/chunk = style name of this compute command +compute-ID = ID of "compute gyration/chunk"_compute_gyration_chunk.html command :ul + +[Examples:] + +compute 1 molecule gyration/shape/chunk pe :pre + +[Description:] + +Define a computation that calculates the eigenvalues of the gyration tensor and +three shape parameters of multiple chunks of atoms. The computation includes +all effects due to atoms passing through periodic boundaries. + +The three computed shape parameters are the asphericity, b, the acylindricity, c, +and the relative shape anisotropy, k: + +:c,image(Eqs/compute_shape_parameters.jpg) + +where lx <= ly <= lz are the three eigenvalues of the gyration tensor. +The asphericity is always non-negative and zero only when the three principal +moments are equal. This zero condition is met when the distribution of particles +is spherically symmetric (hence the name asphericity) but also whenever the particle +distribution is symmetric with respect to the three coordinate axes, e.g., +when the particles are distributed uniformly on a cube, tetrahedron or other Platonic +solid. The acylindricity is always non-negative and zero only when the two principal +moments are equal. This zero condition is met when the distribution of particles is +cylindrically symmetric (hence the name, acylindricity), but also whenever the particle +distribution is symmetric with respect to the two coordinate axes, e.g., when the +particles are distributed uniformly on a regular prism. the relative shape anisotropy +is bounded between zero (if all points are spherically symmetric) and one +(if all points lie on a line). + +The tensor keyword must be specified in the compute gyration/chunk command. + +NOTE: The coordinates of an atom contribute to the gyration tensor in +"unwrapped" form, by using the image flags associated with each atom. +See the "dump custom"_dump.html command for a discussion of "unwrapped" +coordinates. See the Atoms section of the "read_data"_read_data.html +command for a discussion of image flags and how they are set for each +atom. You can reset the image flags (e.g. to 0) before invoking this +compute by using the "set image"_set.html command. + +[Output info:] + +This compute calculates a global array with six columns, +which can be accessed by indices 1-6. The first three columns are the +eigenvalues of the gyration tensor followed by the asphericity, the acylindricity +and the relative shape anisotropy. The computed values can be used by any command +that uses global array values from a compute as input. See the "Howto +output"_Howto_output.html doc page for an overview of LAMMPS output +options. + +The array calculated by this compute is +"intensive". The first five columns will be in +distance^2 "units"_units.html while the sixth one is dimensionless. + +[Restrictions:] + +This compute is part of the USER-MISC package. It is only enabled if +LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +[Related commands:] + +"compute gyration/chunk"_compute_gyration_chunk.html +"compute gyration/shape"_compute_gyration_shape.html + +[Default:] none + +:line + +:link(Theodorou) +[(Theodorou)] Theodorou, Suter, Macromolecules, 18, 1206 (1985). + diff --git a/src/USER-MISC/README b/src/USER-MISC/README index 3b96f3685d..18b8cfb58d 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -31,6 +31,7 @@ compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013 compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017 compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018 compute gyration/shape, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 July 2019 +compute gyration/shape/chunk, Evangelos Voyiatzis, evoyiatzis at gmail.com, 22 October 2019 compute hma, Andrew Schultz & David Kofke (UB), ajs42 at buffalo.edu & kofke at buffalo.edu, 1 Jul 2019 compute pressure/cylinder, Cody K. Addington (NCSU), , 2 Oct 2018 compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk, 28 June 2019 diff --git a/src/USER-MISC/compute_gyration_shape_chunk.cpp b/src/USER-MISC/compute_gyration_shape_chunk.cpp new file mode 100644 index 0000000000..08484d9301 --- /dev/null +++ b/src/USER-MISC/compute_gyration_shape_chunk.cpp @@ -0,0 +1,191 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + * Contributing author: Evangelos Voyiatzis (Royal DSM) + * ------------------------------------------------------------------------- */ + + +#include "compute_gyration_shape_chunk.h" +#include +#include +#include "error.h" +#include "math_extra.h" +#include "math_special.h" +#include "modify.h" +#include "memory.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGyrationShapeChunk::ComputeGyrationShapeChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), id_gyration_chunk(NULL), shape_parameters(NULL) +{ + if (narg != 4) error->all(FLERR,"Illegal compute gyration/shape/chunk command"); + + // ID of compute gyration + int n = strlen(arg[3]) + 1; + id_gyration_chunk = new char[n]; + strcpy(id_gyration_chunk,arg[3]); + + init(); + + array_flag = 1; + size_array_cols = 6; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + firstflag = 1; + former_nchunks = 0; + current_nchunks = 1; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeGyrationShapeChunk::~ComputeGyrationShapeChunk() +{ + delete [] id_gyration_chunk; + memory->destroy(shape_parameters); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGyrationShapeChunk::init() +{ + // check that the compute gyration command exist + int icompute = modify->find_compute(id_gyration_chunk); + if (icompute < 0) + error->all(FLERR,"Compute gyration/chunk ID does not exist for " + "compute gyration/shape/chunk"); + + // check the id_gyration_chunk corresponds really to a compute gyration/chunk command + c_gyration_chunk = (Compute *) modify->compute[icompute]; + if (strcmp(c_gyration_chunk->style,"gyration/chunk") != 0) + error->all(FLERR,"Compute gyration/shape/chunk does not point to " + "gyration compute/chunk"); + + // check the compute gyration/chunk command computes the whole gyration tensor + if (c_gyration_chunk->array_flag == 0) + error->all(FLERR,"Compute gyration/chunk where gyration/shape/chunk points to " + "does not calculate the gyration tensor"); + +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGyrationShapeChunk::setup() +{ + // one-time calculation of per-chunk mass + // done in setup, so that ComputeChunkAtom::setup() is already called + + if (firstflag) { + compute_array(); + firstflag = 0; + } +} + +/* ---------------------------------------------------------------------- + compute shape parameters based on the eigenvalues of the + gyration tensor of group of atoms +------------------------------------------------------------------------- */ + +void ComputeGyrationShapeChunk::compute_array() +{ + invoked_array = update->ntimestep; + c_gyration_chunk->compute_array(); + + current_nchunks = c_gyration_chunk->size_array_rows; // how to check for the number of chunks in the gyration/chunk? + if (former_nchunks != current_nchunks) allocate(); + + double **gyration_tensor = c_gyration_chunk->array; + + // call the function for the calculation of the eigenvalues + double ione[3][3], evalues[3], evectors[3][3]; + + for (int ichunk = 0; ichunk < current_nchunks; ichunk++) { + + ione[0][0] = gyration_tensor[ichunk][0]; + ione[1][1] = gyration_tensor[ichunk][1]; + ione[2][2] = gyration_tensor[ichunk][2]; + ione[0][1] = ione[1][0] = gyration_tensor[ichunk][3]; + ione[0][2] = ione[2][0] = gyration_tensor[ichunk][4]; + ione[1][2] = ione[2][1] = gyration_tensor[ichunk][5]; + + int ierror = MathExtra::jacobi(ione,evalues,evectors); + if (ierror) error->all(FLERR, "Insufficient Jacobi rotations " + "for gyration/shape"); + + // sort the eigenvalues according to their size with bubble sort + double t; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 2-i; j++) { + if (fabs(evalues[j]) < fabs(evalues[j+1])) { + t = evalues[j]; + evalues[j] = evalues[j+1]; + evalues[j+1] = t; + } + } + } + + // compute the shape parameters of the gyration tensor + double sq_eigen_x = MathSpecial::square(evalues[0]); + double sq_eigen_y = MathSpecial::square(evalues[1]); + double sq_eigen_z = MathSpecial::square(evalues[2]); + + double nominator = sq_eigen_x + sq_eigen_y + sq_eigen_z; + double denominator = MathSpecial::square(evalues[0]+evalues[1]+evalues[2]); + + shape_parameters[ichunk][0] = evalues[0]; + shape_parameters[ichunk][1] = evalues[1]; + shape_parameters[ichunk][2] = evalues[2]; + shape_parameters[ichunk][3] = evalues[0] - 0.5*(evalues[1] + evalues[2]); + shape_parameters[ichunk][4] = evalues[1] - evalues[2]; + shape_parameters[ichunk][5] = 1.5*nominator/denominator - 0.5; + } +} + +/* ---------------------------------------------------------------------- + * calculate and return # of chunks = length of vector/array + * ------------------------------------------------------------------------- */ + +int ComputeGyrationShapeChunk::lock_length() +{ + int number_of_chunks = c_gyration_chunk->size_array_rows; + return number_of_chunks; +} + +/* ---------------------------------------------------------------------- + * free and reallocate per-chunk arrays + * ---------------------------------------------------------------------- */ + +void ComputeGyrationShapeChunk::allocate() +{ + memory->destroy(shape_parameters); + former_nchunks = current_nchunks; + memory->create(shape_parameters,current_nchunks,6,"gyration/shape/chunk:shape_parameters"); + array = shape_parameters; +} + +/* ---------------------------------------------------------------------- + * memory usage of local data + * ---------------------------------------------------------------------- */ + +double ComputeGyrationShapeChunk::memory_usage() +{ + double bytes = (bigint) current_nchunks * 6 * sizeof(double); + return bytes; +} diff --git a/src/USER-MISC/compute_gyration_shape_chunk.h b/src/USER-MISC/compute_gyration_shape_chunk.h new file mode 100644 index 0000000000..1ce4b9f758 --- /dev/null +++ b/src/USER-MISC/compute_gyration_shape_chunk.h @@ -0,0 +1,75 @@ +/* -*- 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(gyration/shape/chunk,ComputeGyrationShapeChunk) + +#else + +#ifndef LMP_COMPUTE_GYRATION_SHAPE_CHUNK_H +#define LMP_COMPUTE_GYRATION_SHAPE_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGyrationShapeChunk : public Compute { + public: + char *id_gyration_chunk; // fields accessed by other classes + + ComputeGyrationShapeChunk(class LAMMPS *, int, char **); + ~ComputeGyrationShapeChunk(); + void init(); + void setup(); + void compute_array(); + + int lock_length(); + + double memory_usage(); + + private: + int current_nchunks, former_nchunks; + int firstflag; + double **shape_parameters; + class Compute *c_gyration_chunk; + + void allocate(); + +}; + +} + +#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 gyration/chunk ID does not exist for compute gyration/shape/chunk + +Self-explanatory. Provide a valid compute ID + +E: Compute gyration/shape/chunk ID does not point to a gyration/chunk compute + +Self-explanatory. Provide an ID of a compute gyration/chunk command. + +E: Compute gyration/chunk does not compute gyration tensor + +Self-explanatory. Use keyword tensor in compute gyration/chunk command +*/