lammps-sys 0.6.0

Generates bindings to LAMMPS' C interface (with optional builds from source)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

/* ----------------------------------------------------------------------
   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.h"
#include <cmath>
#include <cstring>
#include "error.h"
#include "math_extra.h"
#include "math_special.h"
#include "modify.h"
#include "update.h"

using namespace LAMMPS_NS;

/* ---------------------------------------------------------------------- */

ComputeGyrationShape::ComputeGyrationShape(LAMMPS *lmp, int narg, char **arg) :
  Compute(lmp, narg, arg), id_gyration(NULL)
{
  if (narg != 4) error->all(FLERR,"Illegal compute gyration/shape command");

  vector_flag = 1;
  size_vector = 6;
  extscalar = 0;
  extvector = 0;

  // ID of compute gyration 
  int n = strlen(arg[3]) + 1;
  id_gyration = new char[n];
  strcpy(id_gyration,arg[3]);

  init();

  vector = new double[6];
}

/* ---------------------------------------------------------------------- */

ComputeGyrationShape::~ComputeGyrationShape()
{
  delete [] id_gyration;
  delete [] vector;
}

/* ---------------------------------------------------------------------- */

void ComputeGyrationShape::init()
{
  // check that the compute gyration command exist
  int icompute = modify->find_compute(id_gyration);
  if (icompute < 0)
    error->all(FLERR,"Compute gyration ID does not exist for "
               "compute gyration/shape");

  // check the id_gyration corresponds really to a compute gyration command
  c_gyration = (Compute *) modify->compute[icompute];
  if (strcmp(c_gyration->style,"gyration") != 0)
    error->all(FLERR,"Compute gyration compute ID does not point to "
               "gyration compute for compute gyration/shape");
}

/* ----------------------------------------------------------------------
   compute shape parameters based on the eigenvalues of the
   gyration tensor of group of atoms
------------------------------------------------------------------------- */

void ComputeGyrationShape::compute_vector()
{
  invoked_vector = update->ntimestep;
  c_gyration->compute_vector();
  double *gyration_tensor = c_gyration->vector;

  // call the function for the calculation of the eigenvalues
  double ione[3][3], evalues[3], evectors[3][3];

  ione[0][0] = gyration_tensor[0];
  ione[1][1] = gyration_tensor[1];
  ione[2][2] = gyration_tensor[2];
  ione[0][1] = ione[1][0] = gyration_tensor[3];
  ione[1][2] = ione[2][1] = gyration_tensor[4];
  ione[0][2] = ione[2][0] = gyration_tensor[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 = MathSpecial::square(sq_eigen_x)
    + MathSpecial::square(sq_eigen_y)
    + MathSpecial::square(sq_eigen_z);
  double denominator = MathSpecial::square(sq_eigen_x+sq_eigen_y+sq_eigen_z);

  vector[0] = evalues[0];
  vector[1] = evalues[1];
  vector[2] = evalues[2];
  vector[3] = sq_eigen_z - 0.5*(sq_eigen_x + sq_eigen_y);
  vector[4] = sq_eigen_y - sq_eigen_x;
  vector[5] = 0.5*(3*nominator/denominator - 1.0);
}