← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/compute_com_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_com_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_com_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
 
ComputeCOMMolecule::ComputeCOMMolecule(LAMMPS *lmp, int narg, char **arg) :
26
 
  Compute(lmp, narg, arg)
27
 
{
28
 
  if (narg != 3) error->all(FLERR,"Illegal compute com/molecule command");
29
 
 
30
 
  if (atom->molecular == 0)
31
 
    error->all(FLERR,"Compute com/molecule requires molecular atom style");
32
 
 
33
 
  array_flag = 1;
34
 
  size_array_cols = 3;
35
 
  extarray = 0;
36
 
 
37
 
  // setup molecule-based data
38
 
 
39
 
  nmolecules = molecules_in_group(idlo,idhi);
40
 
  size_array_rows = nmolecules;
41
 
 
42
 
  memory->create(massproc,nmolecules,"com/molecule:massproc");
43
 
  memory->create(masstotal,nmolecules,"com/molecule:masstotal");
44
 
  memory->create(com,nmolecules,3,"com/molecule:com");
45
 
  memory->create(comall,nmolecules,3,"com/molecule:comall");
46
 
  array = comall;
47
 
 
48
 
  // compute masstotal for each molecule
49
 
 
50
 
  int *mask = atom->mask;
51
 
  tagint *molecule = atom->molecule;
52
 
  int *type = atom->type;
53
 
  double *mass = atom->mass;
54
 
  double *rmass = atom->rmass;
55
 
  int nlocal = atom->nlocal;
56
 
 
57
 
  tagint imol;
58
 
  double massone;
59
 
 
60
 
  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
61
 
 
62
 
  for (int i = 0; i < nlocal; i++)
63
 
    if (mask[i] & groupbit) {
64
 
      if (rmass) massone = rmass[i];
65
 
      else massone = mass[type[i]];
66
 
      imol = molecule[i];
67
 
      if (molmap) imol = molmap[imol-idlo];
68
 
      else imol--;
69
 
      massproc[imol] += massone;
70
 
    }
71
 
 
72
 
  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
73
 
}
74
 
 
75
 
/* ---------------------------------------------------------------------- */
76
 
 
77
 
ComputeCOMMolecule::~ComputeCOMMolecule()
78
 
{
79
 
  memory->destroy(massproc);
80
 
  memory->destroy(masstotal);
81
 
  memory->destroy(com);
82
 
  memory->destroy(comall);
83
 
}
84
 
 
85
 
/* ---------------------------------------------------------------------- */
86
 
 
87
 
void ComputeCOMMolecule::init()
88
 
{
89
 
  int ntmp = molecules_in_group(idlo,idhi);
90
 
  if (ntmp != nmolecules)
91
 
    error->all(FLERR,"Molecule count changed in compute com/molecule");
92
 
}
93
 
 
94
 
/* ---------------------------------------------------------------------- */
95
 
 
96
 
void ComputeCOMMolecule::compute_array()
97
 
{
98
 
  tagint imol;
99
 
  double massone;
100
 
  double unwrap[3];
101
 
 
102
 
  invoked_array = update->ntimestep;
103
 
 
104
 
  for (int i = 0; i < nmolecules; i++)
105
 
    com[i][0] = com[i][1] = com[i][2] = 0.0;
106
 
 
107
 
  double **x = atom->x;
108
 
  int *mask = atom->mask;
109
 
  tagint *molecule = atom->molecule;
110
 
  int *type = atom->type;
111
 
  imageint *image = atom->image;
112
 
  double *mass = atom->mass;
113
 
  double *rmass = atom->rmass;
114
 
  int nlocal = atom->nlocal;
115
 
 
116
 
  for (int i = 0; i < nlocal; i++)
117
 
    if (mask[i] & groupbit) {
118
 
      if (rmass) massone = rmass[i];
119
 
      else massone = mass[type[i]];
120
 
      imol = molecule[i];
121
 
      if (molmap) imol = molmap[imol-idlo];
122
 
      else imol--;
123
 
      domain->unmap(x[i],image[i],unwrap);
124
 
      com[imol][0] += unwrap[0] * massone;
125
 
      com[imol][1] += unwrap[1] * massone;
126
 
      com[imol][2] += unwrap[2] * massone;
127
 
    }
128
 
 
129
 
  MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
130
 
                MPI_DOUBLE,MPI_SUM,world);
131
 
  for (int i = 0; i < nmolecules; i++) {
132
 
    comall[i][0] /= masstotal[i];
133
 
    comall[i][1] /= masstotal[i];
134
 
    comall[i][2] /= masstotal[i];
135
 
  }
136
 
}
137
 
 
138
 
/* ----------------------------------------------------------------------
139
 
   memory usage of local data
140
 
------------------------------------------------------------------------- */
141
 
 
142
 
double ComputeCOMMolecule::memory_usage()
143
 
{
144
 
  double bytes = (bigint) nmolecules * 2 * sizeof(double);
145
 
  if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
146
 
  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
147
 
  return bytes;
148
 
}