1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, 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_dipole.h"
15
16 #include "atom.h"
17 #include "domain.h"
18 #include "error.h"
19 #include "update.h"
20
21 #include <cmath>
22 #include <cstring>
23
24 using namespace LAMMPS_NS;
25
26 enum { MASSCENTER, GEOMCENTER };
27
28 /* ---------------------------------------------------------------------- */
29
ComputeDipole(LAMMPS * lmp,int narg,char ** arg)30 ComputeDipole::ComputeDipole(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
31 {
32 if ((narg < 3) || (narg > 4)) error->all(FLERR, "Illegal compute dipole command");
33
34 scalar_flag = 1;
35 vector_flag = 1;
36 size_vector = 3;
37 extscalar = 0;
38 extvector = 0;
39
40 vector = new double[size_vector];
41 vector[0] = vector[1] = vector[2] = 0.0;
42 usecenter = MASSCENTER;
43
44 if (narg == 4) {
45 if (utils::strmatch(arg[3], "^geom"))
46 usecenter = GEOMCENTER;
47 else if (strcmp(arg[3], "mass") == 0)
48 usecenter = MASSCENTER;
49 else
50 error->all(FLERR, "Illegal compute dipole command");
51 }
52 }
53
54 /* ---------------------------------------------------------------------- */
55
~ComputeDipole()56 ComputeDipole::~ComputeDipole()
57 {
58 delete[] vector;
59 }
60
61 /* ---------------------------------------------------------------------- */
62
compute_vector()63 void ComputeDipole::compute_vector()
64 {
65 invoked_vector = update->ntimestep;
66
67 const auto x = atom->x;
68 const auto mask = atom->mask;
69 const auto type = atom->type;
70 const auto image = atom->image;
71 const auto mass = atom->mass;
72 const auto rmass = atom->rmass;
73 const auto q = atom->q;
74 const auto mu = atom->mu;
75 const auto nlocal = atom->nlocal;
76
77 double dipole[3] = {0.0, 0.0, 0.0};
78 double comproc[3] = {0.0, 0.0, 0.0};
79 double com[3] = {0.0, 0.0, 0.0};
80 double masstotal = 0.0;
81 double chrgtotal = 0.0;
82 double massproc = 0.0;
83 double chrgproc = 0.0;
84
85 for (int i = 0; i < nlocal; ++i) {
86 if (mask[i] & groupbit) {
87 double unwrap[3];
88 double massone = 1.0; // for usecenter == GEOMCENTER
89 if (usecenter == MASSCENTER) {
90 if (rmass)
91 massone = rmass[i];
92 else
93 massone = mass[type[i]];
94 }
95 massproc += massone;
96 if (atom->q_flag) chrgproc += q[i];
97 domain->unmap(x[i], image[i], unwrap);
98 comproc[0] += unwrap[0] * massone;
99 comproc[1] += unwrap[1] * massone;
100 comproc[2] += unwrap[2] * massone;
101 }
102 }
103 MPI_Allreduce(&massproc, &masstotal, 1, MPI_DOUBLE, MPI_SUM, world);
104 MPI_Allreduce(&chrgproc, &chrgtotal, 1, MPI_DOUBLE, MPI_SUM, world);
105 MPI_Allreduce(comproc, com, 3, MPI_DOUBLE, MPI_SUM, world);
106
107 if (masstotal > 0.0) {
108 com[0] /= masstotal;
109 com[1] /= masstotal;
110 com[2] /= masstotal;
111 }
112
113 // compute dipole moment
114
115 for (int i = 0; i < nlocal; i++) {
116 if (mask[i] & groupbit) {
117 double unwrap[3];
118 domain->unmap(x[i], image[i], unwrap);
119 if (atom->q_flag) {
120 dipole[0] += q[i] * unwrap[0];
121 dipole[1] += q[i] * unwrap[1];
122 dipole[2] += q[i] * unwrap[2];
123 }
124 if (atom->mu_flag) {
125 dipole[0] += mu[i][0];
126 dipole[1] += mu[i][1];
127 dipole[2] += mu[i][2];
128 }
129 }
130 }
131
132 MPI_Allreduce(dipole, vector, 3, MPI_DOUBLE, MPI_SUM, world);
133
134 // correct for position dependence with a net charged group
135 vector[0] -= chrgtotal * com[0];
136 vector[1] -= chrgtotal * com[1];
137 vector[2] -= chrgtotal * com[2];
138 }
139
140 /* ---------------------------------------------------------------------- */
141
compute_scalar()142 double ComputeDipole::compute_scalar()
143 {
144 if (invoked_vector != update->ntimestep) compute_vector();
145
146 invoked_scalar = update->ntimestep;
147 scalar = sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]);
148 return scalar;
149 }
150