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 /* ----------------------------------------------------------------------
15 Contributing author: Rene Halver (JSC)
16 ------------------------------------------------------------------------- */
17
18 #include "scafacos.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "domain.h"
23 #include "error.h"
24 #include "force.h"
25 #include "memory.h"
26 #include "pair_hybrid.h"
27
28 #include <cstdlib>
29 #include <cstring>
30
31 // ScaFaCoS library
32
33 #include "fcs.h"
34
35 using namespace LAMMPS_NS;
36
37 /* ---------------------------------------------------------------------- */
38
Scafacos(LAMMPS * lmp)39 Scafacos::Scafacos(LAMMPS *lmp) : KSpace(lmp)
40 {
41 me = comm->me;
42 initialized = 0;
43
44 maxatom = 0;
45 xpbc = nullptr;
46 epot = nullptr;
47 efield = nullptr;
48 fcs = nullptr;
49 }
50
51 /* ---------------------------------------------------------------------- */
52
settings(int narg,char ** arg)53 void Scafacos::settings(int narg, char **arg)
54 {
55 if (narg != 2) error->all(FLERR, "Illegal scafacos command");
56
57 method = utils::strdup(arg[0]);
58 tolerance = utils::numeric(FLERR, arg[1], false, lmp);
59
60 // optional ScaFaCoS library setting defaults
61 // choose the correct default tolerance type for chosen method
62 // throw an error if a not yet supported solver is chosen
63 if (strcmp(method, "fmm") == 0) {
64 tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
65 fmm_tuning_flag = 0;
66 } else if (strcmp(method, "p3m") == 0 || strcmp(method, "p2nfft") == 0 ||
67 strcmp(method, "ewald") == 0) {
68 tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
69 } else if (strcmp(method, "direct") == 0) {
70 ; // direct summation has no tolerance type
71 } else {
72 error->all(FLERR, "Unsupported ScaFaCoS method");
73 }
74 }
75
76 /* ---------------------------------------------------------------------- */
77
~Scafacos()78 Scafacos::~Scafacos()
79 {
80 delete[] method;
81
82 memory->destroy(xpbc);
83 memory->destroy(epot);
84 memory->destroy(efield);
85
86 // clean up of the ScaFaCoS handle and internal arrays
87 fcs_destroy((FCS) fcs);
88 }
89
90 /* ---------------------------------------------------------------------- */
91
init()92 void Scafacos::init()
93 {
94 // error checks
95 if (me == 0) {
96 utils::logmesg(lmp, "Setting up ScaFaCoS with solver {} ...\n", method);
97
98 if (strcmp(method, "p3m") == 0)
99 error->warning(FLERR, "Virial computation for P3M not available");
100
101 if (strcmp(method, "ewald") == 0)
102 error->warning(FLERR, "Virial computation for Ewald not available");
103 }
104
105 if (!atom->q_flag) error->all(FLERR, "Kspace style requires atom attribute q");
106
107 if (domain->dimension == 2) error->all(FLERR, "Cannot use ScaFaCoS with 2d simulation");
108
109 if (domain->triclinic) error->all(FLERR, "Cannot use ScaFaCoS with triclinic domain yet");
110
111 if (atom->natoms > INT_MAX && sizeof(int) != 8)
112 error->all(FLERR, "ScaFaCoS atom count exceeds 2B");
113
114 if ((atom->molecular != Atom::ATOMIC) && (atom->nbonds + atom->nangles + atom->ndihedrals) > 0) {
115 int flag = 0;
116 if ((force->special_coul[1] == 1.0) && (force->special_coul[2] == 1.0) &&
117 (force->special_coul[3] == 1.0))
118 ++flag;
119
120 PairHybrid *ph = reinterpret_cast<PairHybrid *>(force->pair_match("^hybrid", 0));
121 if (ph) {
122 for (int isub = 0; isub < ph->nstyles; ++isub)
123 if (force->pair_match("coul/exclude", 0, isub)) ++flag;
124 } else {
125 if (force->pair_match("coul/exclude", 0)) ++flag;
126 }
127 if (!flag)
128 error->all(FLERR,
129 "Must use pair style coul/exclude or 'special_bonds coul 1.0 1.0 1.0'"
130 "for molecular charged systems with ScaFaCos");
131 }
132
133 FCSResult result;
134
135 // one-time initialization of ScaFaCoS
136
137 qqrd2e = force->qqrd2e;
138
139 if (!initialized) {
140 result = fcs_init((FCS *) &fcs, method, world);
141 check_result((void *) &result);
142
143 setup_handle();
144
145 // using other methods lead to termination of the program,
146 // since they have no tolerance tuning methods
147 if (strcmp(method, "fmm") == 0 || strcmp(method, "p3m") == 0 || strcmp(method, "p2nfft") == 0 ||
148 strcmp(method, "ewald") == 0) {
149 result = fcs_set_tolerance((FCS) fcs, tolerance_type, tolerance);
150 check_result((void *) &result);
151 }
152
153 double **x = atom->x;
154 double *q = atom->q;
155 int nlocal = atom->nlocal;
156
157 if (strcmp(method, "fmm") == 0) {
158 if (fmm_tuning_flag == 1)
159 fcs_fmm_set_internal_tuning((FCS) fcs, FCS_FMM_INHOMOGENOUS_SYSTEM);
160 else
161 fcs_fmm_set_internal_tuning((FCS) fcs, FCS_FMM_HOMOGENOUS_SYSTEM);
162 }
163
164 // for the FMM at least one particle is required per process
165 if (strcmp(method, "fmm") == 0) {
166 int empty = (nlocal == 0) ? 1 : 0;
167 MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPI_INT, MPI_SUM, world);
168 if (empty > 0)
169 fcs_set_redistribute((FCS) fcs, 1);
170 else
171 fcs_set_redistribute((FCS) fcs, 0);
172 }
173
174 result = fcs_tune((FCS) fcs, nlocal, &x[0][0], q);
175 check_result((void *) &result);
176 // more useful here, since the parameters should be tuned now
177 if (me == 0) fcs_print_parameters((FCS) fcs);
178 }
179
180 initialized = 1;
181 }
182
183 /* ---------------------------------------------------------------------- */
184
compute(int eflag,int vflag)185 void Scafacos::compute(int eflag, int vflag)
186 {
187 double **x = atom->x;
188 double *q = atom->q;
189 int nlocal = atom->nlocal;
190
191 const double qscale = qqrd2e;
192 FCSResult result;
193
194 // for the FMM at least one particle is required per process
195 if (strcmp(method, "fmm")) {
196 int empty = (nlocal == 0) ? 1 : 0;
197 MPI_Allreduce(MPI_IN_PLACE, &empty, 1, MPI_INT, MPI_SUM, world);
198 if (empty > 0)
199 fcs_set_redistribute((FCS) fcs, 1);
200 else
201 fcs_set_redistribute((FCS) fcs, 0);
202 }
203
204 ev_init(eflag, vflag);
205
206 // grow xpbc, epot, efield if necessary
207
208 if (nlocal > maxatom || maxatom == 0) {
209 memory->destroy(xpbc);
210 memory->destroy(epot);
211 memory->destroy(efield);
212 maxatom = atom->nmax;
213 memory->create(xpbc, 3 * maxatom, "scafacos:xpbc");
214 memory->create(epot, maxatom, "scafacos:epot");
215 memory->create(efield, maxatom, 3, "scafacos:efield");
216 }
217
218 if (vflag_global) {
219
220 // for P3M or Ewald we cannot compute the virial. skip it.
221
222 if ((strcmp(method, "p3m") != 0) && (strcmp(method, "ewald") != 0)) {
223 result = fcs_set_compute_virial((FCS) fcs, 1);
224 check_result((void *) &result);
225 }
226 }
227
228 // pack coords into xpbc and apply PBC
229 memcpy(xpbc, &x[0][0], 3 * nlocal * sizeof(double));
230
231 if (domain->xperiodic || domain->yperiodic || domain->zperiodic) {
232 int j = 0;
233 for (int i = 0; i < nlocal; i++) {
234 domain->remap(&xpbc[j]);
235 j += 3;
236 }
237 }
238
239 // if simulation box has changed, call fcs_tune()
240
241 if (box_has_changed()) {
242 setup_handle();
243 result = fcs_tune((FCS) fcs, nlocal, xpbc, q);
244 check_result((void *) &result);
245 }
246
247 // invoke ScaFaCoS solver
248
249 result = fcs_run((FCS) fcs, nlocal, xpbc, q, &efield[0][0], epot);
250 check_result((void *) &result);
251
252 // extract virial
253
254 if (vflag_global) {
255
256 // for P3M or Ewald we cannot compute the virial. skip it.
257
258 if ((strcmp(method, "p3m") != 0) && (strcmp(method, "ewald") != 0)) {
259 result = fcs_get_virial((FCS) fcs, virial_int);
260 check_result((void *) &result);
261
262 virial[0] = virial_int[0];
263 virial[1] = virial_int[1];
264 virial[2] = virial_int[2];
265 virial[3] = virial_int[4];
266 virial[4] = virial_int[5];
267 virial[5] = virial_int[8];
268 }
269 }
270
271 // apply Efield to each particle
272 // accumulate total energy
273
274 double **f = atom->f;
275
276 double qone;
277 double myeng = 0.0;
278
279 for (int i = 0; i < nlocal; i++) {
280 qone = q[i] * qscale;
281 f[i][0] += qone * efield[i][0];
282 f[i][1] += qone * efield[i][1];
283 f[i][2] += qone * efield[i][2];
284 myeng += 0.5 * qone * epot[i];
285 }
286
287 if (eflag_atom) {
288 for (int i = 0; i < nlocal; i++) eatom[i] = 0.5 * qscale * q[i] * epot[i];
289 }
290
291 MPI_Allreduce(&myeng, &energy, 1, MPI_DOUBLE, MPI_SUM, world);
292 }
293
294 /* ---------------------------------------------------------------------- */
295
modify_param(int narg,char ** arg)296 int Scafacos::modify_param(int narg, char **arg)
297 {
298 // add any Scafacos options here you want to expose to LAMMPS
299 // syntax: kspace_modify scafacos keyword value1 value2 ...
300 // keyword = tolerance
301 // value1 = energy, energy_rel, etc
302 // everyone of these should have a default, so user doesn't need to set
303
304 if (strcmp(arg[0], "scafacos") != 0) return 0;
305
306 if (strcmp(arg[1], "tolerance") == 0) {
307 if (narg < 3) error->all(FLERR, "Illegal kspace_modify command (tolerance)");
308 if (strcmp(arg[2], "energy") == 0)
309 tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
310 else if (strcmp(arg[2], "energy_rel") == 0)
311 tolerance_type = FCS_TOLERANCE_TYPE_ENERGY_REL;
312 else if (strcmp(arg[2], "field") == 0)
313 tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
314 else if (strcmp(arg[2], "field_rel") == 0)
315 tolerance_type = FCS_TOLERANCE_TYPE_FIELD_REL;
316 else if (strcmp(arg[2], "potential") == 0)
317 tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL;
318 else if (strcmp(arg[2], "potential_rel") == 0)
319 tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL_REL;
320 else
321 error->all(FLERR, "Illegal kspace_modify command (tolerance argument)");
322 // check if method is compatatible to chosen tolerance type
323 if ((strcmp(method, "fmm") == 0 &&
324 (tolerance_type != FCS_TOLERANCE_TYPE_ENERGY &&
325 tolerance_type != FCS_TOLERANCE_TYPE_ENERGY_REL)) ||
326 (strcmp(method, "p2nfft") == 0 &&
327 (tolerance_type != FCS_TOLERANCE_TYPE_FIELD &&
328 tolerance_type != FCS_TOLERANCE_TYPE_POTENTIAL)) ||
329 (strcmp(method, "p3m") == 0 && (tolerance_type != FCS_TOLERANCE_TYPE_FIELD)) ||
330 (strcmp(method, "ewald") == 0 && (tolerance_type != FCS_TOLERANCE_TYPE_FIELD)))
331 error->all(FLERR, "Illegal kspace_modify command \
332 (invalid tolerance / method combination)");
333 return 3;
334 }
335
336 // keyword = fmm_inhomogen_tuning
337 // value1 = 0, 1
338 // 0 -> homogenous system (default)
339 // 1 -> inhomogenous system (more internal tuning is provided (sequential!))
340 if (strcmp(arg[1], "fmm_tuning") == 0) {
341 if (me == 0) utils::logmesg(lmp, "ScaFaCoS setting fmm inhomogen tuning ...\n");
342
343 if (narg < 3) error->all(FLERR, "Illegal kspace_modify command (fmm_tuning)");
344 fmm_tuning_flag = utils::inumeric(FLERR, arg[2], false, lmp);
345 return 3;
346 }
347
348 return 0;
349 }
350
351 /* ----------------------------------------------------------------------
352 memory usage of local arrays
353 ------------------------------------------------------------------------- */
354
memory_usage()355 double Scafacos::memory_usage()
356 {
357 double bytes = 0.0;
358 bytes += (double) maxatom * sizeof(double);
359 bytes += (double) 3 * maxatom * sizeof(double);
360 return bytes;
361 }
362
363 /* ----------------------------------------------------------------------
364 setup of ScaFaCoS handle with common parameters
365 ------------------------------------------------------------------------- */
366
setup_handle()367 void Scafacos::setup_handle()
368 {
369 FCSResult result;
370
371 // store simulation box params
372
373 // setup periodicity
374 old_periodicity[0] = domain->xperiodic;
375 old_periodicity[1] = domain->yperiodic;
376 old_periodicity[2] = domain->zperiodic;
377
378 // setup box origin (lower left front corner of the system)
379 old_origin[0] = domain->boxlo[0];
380 old_origin[1] = domain->boxlo[1];
381 old_origin[2] = domain->boxlo[2];
382
383 // setup box vectors (base vectors of the system box)
384 old_box_x[0] = domain->prd[0];
385 old_box_x[1] = old_box_x[2] = 0.0;
386 old_box_y[1] = domain->prd[1];
387 old_box_y[0] = old_box_y[2] = 0.0;
388 old_box_z[2] = domain->prd[2];
389 old_box_z[1] = old_box_z[0] = 0.0;
390
391 // setup number of atoms in the system
392 old_natoms = atom->natoms;
393
394 // store parameters to ScaFaCoS handle
395 result = fcs_set_box_a((FCS) fcs, old_box_x);
396 check_result((void *) &result);
397
398 result = fcs_set_box_b((FCS) fcs, old_box_y);
399 check_result((void *) &result);
400
401 result = fcs_set_box_c((FCS) fcs, old_box_z);
402 check_result((void *) &result);
403
404 result = fcs_set_box_origin((FCS) fcs, old_origin);
405 check_result((void *) &result);
406
407 result = fcs_set_periodicity((FCS) fcs, old_periodicity);
408 check_result((void *) &result);
409
410 result = fcs_set_total_particles((FCS) fcs, old_natoms);
411 check_result((void *) &result);
412
413 // allow ScaFaCoS to calculate the near field computations for now
414 // TODO: allow the delegation of the near field computations
415 // within LAMMPS
416 // (near_field_flag = 1 -> enables the internal near field calcs
417 // 0 -> disables the internal near field calcs
418 int near_field_flag = 1;
419 result = fcs_set_near_field_flag((FCS) fcs, near_field_flag);
420 check_result((void *) &result);
421 }
422
423 /* ----------------------------------------------------------------------
424 check if box parameters changed, requiring a new call to fcs_tune
425 ------------------------------------------------------------------------- */
426
box_has_changed()427 bool Scafacos::box_has_changed()
428 {
429 int *periodicity = domain->periodicity;
430 double *prd = domain->prd;
431
432 bool changed = (periodicity[0] != old_periodicity[0]) || (periodicity[1] != old_periodicity[1]) ||
433 (periodicity[2] != old_periodicity[2]) || (domain->boundary[0][0] != old_origin[0]) ||
434 (domain->boundary[1][0] != old_origin[1]) || (domain->boundary[2][0] != old_origin[2]) ||
435 (prd[0] != old_box_x[0]) || (prd[1] != old_box_y[1]) || (prd[2] != old_box_z[2]) ||
436 (atom->natoms != old_natoms);
437
438 return changed;
439 }
440
441 /* ----------------------------------------------------------------------
442 check ScaFaCoS result for error condition
443 ------------------------------------------------------------------------- */
444
check_result(void * result_p)445 void Scafacos::check_result(void *result_p)
446 {
447 FCSResult result = *(FCSResult *) result_p;
448 if (!result) return;
449
450 error->one(FLERR, "ScaFaCoS: {}\n{}\n", fcs_result_get_function(result),
451 fcs_result_get_message(result));
452 }
453