1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <mpi.h>
47 #include <string.h>
48 #include <stdlib.h>
49 #include "compute_pressure.h"
50 #include "atom.h"
51 #include "update.h"
52 #include "domain.h"
53 #include "modify.h"
54 #include "fix.h"
55 #include "force.h"
56 #include "pair.h"
57 #include "bond.h"
58 #include "angle.h"
59 #include "dihedral.h"
60 #include "improper.h"
61 #include "kspace.h"
62 #include "error.h"
63
64 using namespace LAMMPS_NS;
65
66 /* ---------------------------------------------------------------------- */
67
ComputePressure(LAMMPS * lmp,int & iarg,int narg,char ** arg)68 ComputePressure::ComputePressure(LAMMPS *lmp, int &iarg, int narg, char **arg) :
69 Compute(lmp, iarg, narg, arg)
70 {
71 if (narg < iarg+1) error->all(FLERR,"Illegal compute pressure command");
72 if (igroup) error->all(FLERR,"Compute pressure must use group all");
73
74 scalar_flag = vector_flag = 1;
75 size_vector = 6;
76 extscalar = 0;
77 extvector = 0;
78 pressflag = 1;
79 timeflag = 1;
80
81 // store temperature ID used by pressure computation
82 // insure it is valid for temperature computation
83
84 int n = strlen(arg[iarg]) + 1;
85 id_temp = new char[n];
86 strcpy(id_temp,arg[iarg]);
87 iarg++;
88
89 int icompute = modify->find_compute(id_temp);
90 if (icompute < 0)
91 error->all(FLERR,"Could not find compute pressure temperature ID");
92 if (modify->compute[icompute]->tempflag == 0)
93 error->all(FLERR,
94 "Compute pressure temperature ID does not compute temperature");
95
96 // process optional args
97
98 if (narg == iarg) {
99 keflag = 1;
100 pairflag = 1;
101 bondflag = angleflag = dihedralflag = improperflag = 1;
102 kspaceflag = fixflag = 1;
103 } else {
104 keflag = 0;
105 pairflag = 0;
106 bondflag = angleflag = dihedralflag = improperflag = 0;
107 kspaceflag = fixflag = 0;
108 while (iarg < narg) {
109 if (strcmp(arg[iarg],"ke") == 0) keflag = 1;
110 else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
111 else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
112 else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
113 else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
114 else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
115 else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
116 else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
117 else if (strcmp(arg[iarg],"virial") == 0) {
118 pairflag = 1;
119 bondflag = angleflag = dihedralflag = improperflag = 1;
120 kspaceflag = fixflag = 1;
121 } else error->all(FLERR,"Illegal compute pressure command");
122 iarg++;
123 }
124 }
125
126 vector = new double[6];
127 nvirial = 0;
128 vptr = NULL;
129 }
130
131 /* ---------------------------------------------------------------------- */
132
~ComputePressure()133 ComputePressure::~ComputePressure()
134 {
135 delete [] id_temp;
136 delete [] vector;
137 delete [] vptr;
138 }
139
140 /* ---------------------------------------------------------------------- */
141
init()142 void ComputePressure::init()
143 {
144 boltz = force->boltz;
145 nktv2p = force->nktv2p;
146 dimension = domain->dimension;
147
148 // set temperature compute, must be done in init()
149 // fixes could have changed or compute_modify could have changed it
150
151 int icompute = modify->find_compute(id_temp);
152 if (icompute < 0)
153 error->all(FLERR,"Could not find compute pressure temperature ID");
154 temperature = modify->compute[icompute];
155
156 // detect contributions to virial
157 // vptr points to all virial[6] contributions
158
159 delete [] vptr;
160 nvirial = 0;
161 vptr = NULL;
162
163 if (pairflag && force->pair) nvirial++;
164 if (bondflag && atom->molecular && force->bond) nvirial++;
165 if (angleflag && atom->molecular && force->angle) nvirial++;
166 if (dihedralflag && atom->molecular && force->dihedral) nvirial++;
167 if (improperflag && atom->molecular && force->improper) nvirial++;
168 if (fixflag)
169 for (int i = 0; i < modify->nfix; i++)
170 if (modify->fix[i]->virial_flag) nvirial++;
171
172 if (nvirial) {
173 vptr = new double*[nvirial];
174 nvirial = 0;
175 if (pairflag && force->pair) vptr[nvirial++] = force->pair->virial;
176 if (bondflag && force->bond) vptr[nvirial++] = force->bond->virial;
177 if (angleflag && force->angle) vptr[nvirial++] = force->angle->virial;
178 if (dihedralflag && force->dihedral)
179 vptr[nvirial++] = force->dihedral->virial;
180 if (improperflag && force->improper)
181 vptr[nvirial++] = force->improper->virial;
182 if (fixflag)
183 for (int i = 0; i < modify->nfix; i++)
184 if (modify->fix[i]->virial_flag)
185 vptr[nvirial++] = modify->fix[i]->virial;
186 }
187
188 // flag Kspace contribution separately, since not summed across procs
189
190 if (kspaceflag && force->kspace) kspace_virial = force->kspace->virial;
191 else kspace_virial = NULL;
192 }
193
194 /* ----------------------------------------------------------------------
195 compute total pressure, averaged over Pxx, Pyy, Pzz
196 ------------------------------------------------------------------------- */
197
compute_scalar()198 double ComputePressure::compute_scalar()
199 {
200 invoked_scalar = update->ntimestep;
201 if (update->vflag_global != invoked_scalar)
202 error->all(FLERR,"Virial was not tallied on needed timestep");
203
204 // invoke temperature it it hasn't been already
205
206 double t = 0.0;
207 if (keflag) {
208 if (temperature->invoked_scalar != update->ntimestep)
209 t = temperature->compute_scalar();
210 else t = temperature->scalar;
211 }
212
213 if (dimension == 3) {
214 inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
215 virial_compute(3,3);
216 if (keflag)
217 scalar = (temperature->dof * boltz * t +
218 virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
219 else
220 scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
221 } else {
222 inv_volume = 1.0 / (domain->xprd * domain->yprd);
223 virial_compute(2,2);
224 if (keflag)
225 scalar = (temperature->dof * boltz * t +
226 virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
227 else
228 scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
229 }
230
231 return scalar;
232 }
233
234 /* ----------------------------------------------------------------------
235 compute pressure tensor
236 assume KE tensor has already been computed
237 ------------------------------------------------------------------------- */
238
compute_vector()239 void ComputePressure::compute_vector()
240 {
241 invoked_vector = update->ntimestep;
242 if (update->vflag_global != invoked_vector)
243 error->all(FLERR,"Virial was not tallied on needed timestep");
244
245 // invoke temperature if it hasn't been already
246
247 double *ke_tensor = NULL;
248 if (keflag) {
249 if (temperature->invoked_vector != update->ntimestep)
250 temperature->compute_vector();
251 ke_tensor = temperature->vector;
252 }
253
254 if (dimension == 3) {
255 inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
256 virial_compute(6,3);
257 if (keflag) {
258 for (int i = 0; i < 6; i++)
259 vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;
260 } else
261 for (int i = 0; i < 6; i++)
262 vector[i] = virial[i] * inv_volume * nktv2p;
263 } else {
264 inv_volume = 1.0 / (domain->xprd * domain->yprd);
265 virial_compute(4,2);
266 if (keflag) {
267 vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;
268 vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;
269 vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;
270 vector[2] = vector[4] = vector[5] = 0.0;
271 } else {
272 vector[0] = virial[0] * inv_volume * nktv2p;
273 vector[1] = virial[1] * inv_volume * nktv2p;
274 vector[3] = virial[3] * inv_volume * nktv2p;
275 vector[2] = vector[4] = vector[5] = 0.0;
276 }
277 }
278 }
279
280 /* ---------------------------------------------------------------------- */
281
virial_compute(int n,int ndiag)282 void ComputePressure::virial_compute(int n, int ndiag)
283 {
284 int i,j;
285 double v[6],*vcomponent;
286
287 for (i = 0; i < n; i++) v[i] = 0.0;
288
289 // sum contributions to virial from forces and fixes
290
291 for (j = 0; j < nvirial; j++) {
292 vcomponent = vptr[j];
293 for (i = 0; i < n; i++) v[i] += vcomponent[i];
294 }
295
296 // sum virial across procs
297
298 MPI_Allreduce(v,virial,n,MPI_DOUBLE,MPI_SUM,world);
299
300 // KSpace virial contribution is already summed across procs
301
302 if (kspace_virial)
303 for (i = 0; i < n; i++) virial[i] += kspace_virial[i];
304
305 // LJ long-range tail correction
306
307 if (force->pair && force->pair->tail_flag)
308 for (i = 0; i < ndiag; i++) virial[i] += force->pair->ptail * inv_volume;
309 }
310
311 /* ---------------------------------------------------------------------- */
312
reset_extra_compute_fix(const char * id_new)313 void ComputePressure::reset_extra_compute_fix(const char *id_new)
314 {
315 delete [] id_temp;
316 int n = strlen(id_new) + 1;
317 id_temp = new char[n];
318 strcpy(id_temp,id_new);
319 }
320