← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to src/compute_inertia_molecule.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/compute_inertia_molecule.cpp
1
 
/* ----------------------------------------------------------------------
2
 
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3
 
   http://lammps.sandia.gov, Sandia National Laboratories
4
 
   Steve Plimpton, sjplimp@sandia.gov
5
 
 
6
 
   Copyright (2003) Sandia Corporation.  Under the terms of Contract
7
 
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8
 
   certain rights in this software.  This software is distributed under
9
 
   the GNU General Public License.
10
 
 
11
 
   See the README file in the top-level LAMMPS directory.
12
 
------------------------------------------------------------------------- */
13
 
 
14
 
#include "compute_inertia_molecule.h"
15
 
#include "atom.h"
16
 
#include "update.h"
17
 
#include "domain.h"
18
 
#include "memory.h"
19
 
#include "error.h"
20
 
 
21
 
using namespace LAMMPS_NS;
22
 
 
23
 
/* ---------------------------------------------------------------------- */
24
 
 
25
 
ComputeInertiaMolecule::
26
 
ComputeInertiaMolecule(LAMMPS *lmp, int narg, char **arg) : 
27
 
  Compute(lmp, narg, arg)
28
 
{
29
 
  if (narg != 3) error->all(FLERR,"Illegal compute inertia/molecule command");
30
 
 
31
 
  if (atom->molecular == 0)
32
 
    error->all(FLERR,"Compute inertia/molecule requires molecular atom style");
33
 
 
34
 
  array_flag = 1;
35
 
  size_array_cols = 6;
36
 
  extarray = 0;
37
 
 
38
 
  // setup molecule-based data
39
 
 
40
 
  nmolecules = molecules_in_group(idlo,idhi);
41
 
  size_array_rows = nmolecules;
42
 
 
43
 
  memory->create(massproc,nmolecules,"inertia/molecule:massproc");
44
 
  memory->create(masstotal,nmolecules,"inertia/molecule:masstotal");
45
 
  memory->create(com,nmolecules,3,"inertia/molecule:com");
46
 
  memory->create(comall,nmolecules,3,"inertia/molecule:comall");
47
 
  memory->create(inertia,nmolecules,6,"inertia/molecule:inertia");
48
 
  memory->create(inertiaall,nmolecules,6,"inertia/molecule:inertiaall");
49
 
  array = inertiaall;
50
 
 
51
 
  // compute masstotal for each molecule
52
 
 
53
 
  int *mask = atom->mask;
54
 
  tagint *molecule = atom->molecule;
55
 
  int *type = atom->type;
56
 
  double *mass = atom->mass;
57
 
  double *rmass = atom->rmass;
58
 
  int nlocal = atom->nlocal;
59
 
 
60
 
  tagint imol;
61
 
  double massone;
62
 
 
63
 
  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
64
 
 
65
 
  for (int i = 0; i < nlocal; i++)
66
 
    if (mask[i] & groupbit) {
67
 
      if (rmass) massone = rmass[i];
68
 
      else massone = mass[type[i]];
69
 
      imol = molecule[i];
70
 
      if (molmap) imol = molmap[imol-idlo];
71
 
      else imol--;
72
 
      massproc[imol] += massone;
73
 
    }
74
 
 
75
 
  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
76
 
}
77
 
 
78
 
/* ---------------------------------------------------------------------- */
79
 
 
80
 
ComputeInertiaMolecule::~ComputeInertiaMolecule()
81
 
{
82
 
  memory->destroy(massproc);
83
 
  memory->destroy(masstotal);
84
 
  memory->destroy(com);
85
 
  memory->destroy(comall);
86
 
  memory->destroy(inertia);
87
 
  memory->destroy(inertiaall);
88
 
}
89
 
 
90
 
/* ---------------------------------------------------------------------- */
91
 
 
92
 
void ComputeInertiaMolecule::init()
93
 
{
94
 
  int ntmp = molecules_in_group(idlo,idhi);
95
 
  if (ntmp != nmolecules)
96
 
    error->all(FLERR,"Molecule count changed in compute inertia/molecule");
97
 
}
98
 
 
99
 
/* ---------------------------------------------------------------------- */
100
 
 
101
 
void ComputeInertiaMolecule::compute_array()
102
 
{
103
 
  int i,j;
104
 
  tagint imol;
105
 
  double dx,dy,dz,massone;
106
 
  double unwrap[3];
107
 
 
108
 
  invoked_array = update->ntimestep;
109
 
 
110
 
  double **x = atom->x;
111
 
  int *mask = atom->mask;
112
 
  tagint *molecule = atom->molecule;
113
 
  int *type = atom->type;
114
 
  imageint *image = atom->image;
115
 
  double *mass = atom->mass;
116
 
  double *rmass = atom->rmass;
117
 
  int nlocal = atom->nlocal;
118
 
 
119
 
  // center-of-mass for each molecule
120
 
 
121
 
  for (i = 0; i < nmolecules; i++)
122
 
    com[i][0] = com[i][1] = com[i][2] = 0.0;
123
 
 
124
 
  for (int i = 0; i < nlocal; i++)
125
 
    if (mask[i] & groupbit) {
126
 
      if (rmass) massone = rmass[i];
127
 
      else massone = mass[type[i]];
128
 
      imol = molecule[i];
129
 
      if (molmap) imol = molmap[imol-idlo];
130
 
      else imol--;
131
 
      domain->unmap(x[i],image[i],unwrap);
132
 
      com[imol][0] += unwrap[0] * massone;
133
 
      com[imol][1] += unwrap[1] * massone;
134
 
      com[imol][2] += unwrap[2] * massone;
135
 
    }
136
 
 
137
 
  MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
138
 
                MPI_DOUBLE,MPI_SUM,world);
139
 
  for (i = 0; i < nmolecules; i++) {
140
 
    comall[i][0] /= masstotal[i];
141
 
    comall[i][1] /= masstotal[i];
142
 
    comall[i][2] /= masstotal[i];
143
 
  }
144
 
 
145
 
  // inertia tensor for each molecule
146
 
 
147
 
  for (i = 0; i < nmolecules; i++)
148
 
    for (j = 0; j < 6; j++)
149
 
      inertia[i][j] = 0.0;
150
 
 
151
 
  for (i = 0; i < nlocal; i++)
152
 
    if (mask[i] & groupbit) {
153
 
      if (rmass) massone = rmass[i];
154
 
      else massone = mass[type[i]];
155
 
      imol = molecule[i];
156
 
      if (molmap) imol = molmap[imol-idlo];
157
 
      else imol--;
158
 
      domain->unmap(x[i],image[i],unwrap);
159
 
      dx = unwrap[0] - comall[imol][0];
160
 
      dy = unwrap[1] - comall[imol][1];
161
 
      dz = unwrap[2] - comall[imol][2];
162
 
      inertia[imol][0] += massone * (dy*dy + dz*dz);
163
 
      inertia[imol][1] += massone * (dx*dx + dz*dz);
164
 
      inertia[imol][2] += massone * (dx*dx + dy*dy);
165
 
      inertia[imol][3] -= massone * dx*dy;
166
 
      inertia[imol][4] -= massone * dy*dz;
167
 
      inertia[imol][5] -= massone * dx*dz;
168
 
    }
169
 
 
170
 
  MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nmolecules,
171
 
                MPI_DOUBLE,MPI_SUM,world);
172
 
}
173
 
 
174
 
/* ----------------------------------------------------------------------
175
 
   memory usage of local data
176
 
------------------------------------------------------------------------- */
177
 
 
178
 
double ComputeInertiaMolecule::memory_usage()
179
 
{
180
 
  double bytes = (bigint) nmolecules * 2 * sizeof(double);
181
 
  if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
182
 
  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
183
 
  bytes += (bigint) nmolecules * 2*6 * sizeof(double);
184
 
  return bytes;
185
 
}