1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 #include "kspace.h"
16
17 #include "atom.h"
18 #include "atom_masks.h"
19 #include "comm.h"
20 #include "domain.h"
21 #include "error.h"
22 #include "force.h"
23 #include "memory.h"
24 #include "pair.h"
25 #include "suffix.h"
26
27 #include <cmath>
28 #include <cstring>
29
30 using namespace LAMMPS_NS;
31
32 #define SMALL 0.00001
33
34 /* ---------------------------------------------------------------------- */
35
KSpace(LAMMPS * lmp)36 KSpace::KSpace(LAMMPS *lmp) : Pointers(lmp)
37 {
38 order_allocated = 0;
39 energy = 0.0;
40 virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
41
42 triclinic_support = 1;
43 ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag =
44 dipoleflag = spinflag = 0;
45 compute_flag = 1;
46 group_group_enable = 0;
47 stagger_flag = 0;
48
49 order = 5;
50 gridflag = 0;
51 gewaldflag = 0;
52 minorder = 2;
53 overlap_allowed = 1;
54 fftbench = 0;
55
56 // default to using MPI collectives for FFT/remap only on IBM BlueGene
57
58 #ifdef __bg__
59 collective_flag = 1;
60 #else
61 collective_flag = 0;
62 #endif
63
64 kewaldflag = 0;
65
66 order_6 = 5;
67 gridflag_6 = 0;
68 gewaldflag_6 = 0;
69 auto_disp_flag = 0;
70
71 slabflag = 0;
72 differentiation_flag = 0;
73 slab_volfactor = 1;
74 suffix_flag = Suffix::NONE;
75 adjust_cutoff_flag = 1;
76 scalar_pressure_flag = 0;
77 warn_nonneutral = 1;
78 warn_nocharge = 1;
79
80 accuracy_absolute = -1.0;
81 accuracy_real_6 = -1.0;
82 accuracy_kspace_6 = -1.0;
83
84 neighrequest_flag = 1;
85 mixflag = 0;
86
87 splittol = 1.0e-6;
88
89 maxeatom = maxvatom = 0;
90 eatom = nullptr;
91 vatom = nullptr;
92 centroidstressflag = CENTROID_NOTAVAIL;
93
94 execution_space = Host;
95 datamask_read = ALL_MASK;
96 datamask_modify = ALL_MASK;
97 copymode = 0;
98
99 memory->create(gcons,7,7,"kspace:gcons");
100 gcons[2][0] = 15.0 / 8.0;
101 gcons[2][1] = -5.0 / 4.0;
102 gcons[2][2] = 3.0 / 8.0;
103 gcons[3][0] = 35.0 / 16.0;
104 gcons[3][1] = -35.0 / 16.0;
105 gcons[3][2] = 21.0 / 16.0;
106 gcons[3][3] = -5.0 / 16.0;
107 gcons[4][0] = 315.0 / 128.0;
108 gcons[4][1] = -105.0 / 32.0;
109 gcons[4][2] = 189.0 / 64.0;
110 gcons[4][3] = -45.0 / 32.0;
111 gcons[4][4] = 35.0 / 128.0;
112 gcons[5][0] = 693.0 / 256.0;
113 gcons[5][1] = -1155.0 / 256.0;
114 gcons[5][2] = 693.0 / 128.0;
115 gcons[5][3] = -495.0 / 128.0;
116 gcons[5][4] = 385.0 / 256.0;
117 gcons[5][5] = -63.0 / 256.0;
118 gcons[6][0] = 3003.0 / 1024.0;
119 gcons[6][1] = -3003.0 / 512.0;
120 gcons[6][2] = 9009.0 / 1024.0;
121 gcons[6][3] = -2145.0 / 256.0;
122 gcons[6][4] = 5005.0 / 1024.0;
123 gcons[6][5] = -819.0 / 512.0;
124 gcons[6][6] = 231.0 / 1024.0;
125
126 memory->create(dgcons,7,6,"kspace:dgcons");
127 dgcons[2][0] = -5.0 / 2.0;
128 dgcons[2][1] = 3.0 / 2.0;
129 dgcons[3][0] = -35.0 / 8.0;
130 dgcons[3][1] = 21.0 / 4.0;
131 dgcons[3][2] = -15.0 / 8.0;
132 dgcons[4][0] = -105.0 / 16.0;
133 dgcons[4][1] = 189.0 / 16.0;
134 dgcons[4][2] = -135.0 / 16.0;
135 dgcons[4][3] = 35.0 / 16.0;
136 dgcons[5][0] = -1155.0 / 128.0;
137 dgcons[5][1] = 693.0 / 32.0;
138 dgcons[5][2] = -1485.0 / 64.0;
139 dgcons[5][3] = 385.0 / 32.0;
140 dgcons[5][4] = -315.0 / 128.0;
141 dgcons[6][0] = -3003.0 / 256.0;
142 dgcons[6][1] = 9009.0 / 256.0;
143 dgcons[6][2] = -6435.0 / 128.0;
144 dgcons[6][3] = 5005.0 / 128.0;
145 dgcons[6][4] = -4095.0 / 256.0;
146 dgcons[6][5] = 693.0 / 256.0;
147 }
148
149 /* ---------------------------------------------------------------------- */
150
~KSpace()151 KSpace::~KSpace()
152 {
153 if (copymode) return;
154
155 memory->destroy(eatom);
156 memory->destroy(vatom);
157 memory->destroy(gcons);
158 memory->destroy(dgcons);
159 }
160
161 /* ----------------------------------------------------------------------
162 calculate this in init() so that units are finalized
163 ------------------------------------------------------------------------- */
164
two_charge()165 void KSpace::two_charge()
166 {
167 two_charge_force = force->qqr2e *
168 (force->qelectron * force->qelectron) /
169 (force->angstrom * force->angstrom);
170 }
171
172 /* ---------------------------------------------------------------------- */
173
triclinic_check()174 void KSpace::triclinic_check()
175 {
176 if (domain->triclinic && triclinic_support != 1)
177 error->all(FLERR,"KSpace style does not yet support triclinic geometries");
178 }
179
180 /* ---------------------------------------------------------------------- */
181
compute_dummy(int eflag,int vflag)182 void KSpace::compute_dummy(int eflag, int vflag)
183 {
184 ev_init(eflag,vflag);
185 }
186
187 /* ----------------------------------------------------------------------
188 check that pair style is compatible with long-range solver
189 ------------------------------------------------------------------------- */
190
pair_check()191 void KSpace::pair_check()
192 {
193 if (force->pair == nullptr)
194 error->all(FLERR,"KSpace solver requires a pair style");
195
196 if (ewaldflag && !force->pair->ewaldflag)
197 error->all(FLERR,"KSpace style is incompatible with Pair style");
198 if (pppmflag && !force->pair->pppmflag)
199 error->all(FLERR,"KSpace style is incompatible with Pair style");
200 if (msmflag && !force->pair->msmflag)
201 error->all(FLERR,"KSpace style is incompatible with Pair style");
202 if (dispersionflag && !force->pair->dispersionflag)
203 error->all(FLERR,"KSpace style is incompatible with Pair style");
204 if (dipoleflag && !force->pair->dipoleflag)
205 error->all(FLERR,"KSpace style is incompatible with Pair style");
206 if (spinflag && !force->pair->spinflag)
207 error->all(FLERR,"KSpace style is incompatible with Pair style");
208 if (tip4pflag && !force->pair->tip4pflag)
209 error->all(FLERR,"KSpace style is incompatible with Pair style");
210
211 if (force->pair->dispersionflag && !dispersionflag)
212 error->all(FLERR,"KSpace style is incompatible with Pair style");
213 if (force->pair->tip4pflag && !tip4pflag)
214 error->all(FLERR,"KSpace style is incompatible with Pair style");
215 }
216
217 /* ----------------------------------------------------------------------
218 setup for energy, virial computation
219 see integrate::ev_set() for bitwise settings of eflag/vflag
220 set the following flags, values are otherwise set to 0:
221 evflag != 0 if any bits of eflag or vflag are set
222 eflag_global != 0 if ENERGY_GLOBAL bit of eflag set
223 eflag_atom != 0 if ENERGY_ATOM bit of eflag set
224 eflag_either != 0 if eflag_global or eflag_atom is set
225 vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set
226 vflag_atom != 0 if VIRIAL_ATOM bit of vflag set
227 no current support for centroid stress
228 vflag_either != 0 if vflag_global or vflag_atom is set
229 evflag_atom != 0 if eflag_atom or vflag_atom is set
230 ------------------------------------------------------------------------- */
231
ev_setup(int eflag,int vflag,int alloc)232 void KSpace::ev_setup(int eflag, int vflag, int alloc)
233 {
234 int i,n;
235
236 evflag = 1;
237
238 eflag_either = eflag;
239 eflag_global = eflag & ENERGY_GLOBAL;
240 eflag_atom = eflag & ENERGY_ATOM;
241
242 vflag_either = vflag;
243 vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR);
244 vflag_atom = vflag & VIRIAL_ATOM;
245
246 if (eflag_atom || vflag_atom) evflag_atom = 1;
247 else evflag_atom = 0;
248
249 // reallocate per-atom arrays if necessary
250
251 if (eflag_atom && atom->nmax > maxeatom) {
252 maxeatom = atom->nmax;
253 if (alloc) {
254 memory->destroy(eatom);
255 memory->create(eatom,maxeatom,"kspace:eatom");
256 }
257 }
258 if (vflag_atom && atom->nmax > maxvatom) {
259 maxvatom = atom->nmax;
260 if (alloc) {
261 memory->destroy(vatom);
262 memory->create(vatom,maxvatom,6,"kspace:vatom");
263 }
264 }
265
266 // zero accumulators
267
268 if (eflag_global) energy = 0.0;
269 if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
270 if (eflag_atom && alloc) {
271 n = atom->nlocal;
272 if (tip4pflag) n += atom->nghost;
273 for (i = 0; i < n; i++) eatom[i] = 0.0;
274 }
275 if (vflag_atom && alloc) {
276 n = atom->nlocal;
277 if (tip4pflag) n += atom->nghost;
278 for (i = 0; i < n; i++) {
279 vatom[i][0] = 0.0;
280 vatom[i][1] = 0.0;
281 vatom[i][2] = 0.0;
282 vatom[i][3] = 0.0;
283 vatom[i][4] = 0.0;
284 vatom[i][5] = 0.0;
285 }
286 }
287 }
288
289 /* ----------------------------------------------------------------------
290 compute qsum,qsqsum,q2 and give error/warning if not charge neutral
291 called initially, when particle count changes, when charges are changed
292 ------------------------------------------------------------------------- */
293
qsum_qsq(int warning_flag)294 void KSpace::qsum_qsq(int warning_flag)
295 {
296 const double * const q = atom->q;
297 const int nlocal = atom->nlocal;
298 double qsum_local(0.0), qsqsum_local(0.0);
299
300 #if defined(_OPENMP)
301 #pragma omp parallel for default(shared) reduction(+:qsum_local,qsqsum_local)
302 #endif
303 for (int i = 0; i < nlocal; i++) {
304 qsum_local += q[i];
305 qsqsum_local += q[i]*q[i];
306 }
307
308 MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
309 MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
310
311 if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge && warning_flag) {
312 error->warning(FLERR,"Using kspace solver on system with no charge");
313 warn_nocharge = 0;
314 }
315
316 q2 = qsqsum * force->qqrd2e;
317
318 // not yet sure of the correction needed for non-neutral systems
319 // so issue warning or error
320
321 if (fabs(qsum) > SMALL) {
322 std::string message = fmt::format("System is not charge neutral, net "
323 "charge = {:.8}",qsum);
324 if (!warn_nonneutral) error->all(FLERR,message);
325 if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,message);
326 warn_nonneutral = 2;
327 }
328 }
329
330 /* ----------------------------------------------------------------------
331 estimate the accuracy of the short-range coulomb tables
332 ------------------------------------------------------------------------- */
333
estimate_table_accuracy(double q2_over_sqrt,double spr)334 double KSpace::estimate_table_accuracy(double q2_over_sqrt, double spr)
335 {
336 double table_accuracy = 0.0;
337 int nctb = force->pair->ncoultablebits;
338 if (comm->me == 0) {
339 if (nctb)
340 error->message(FLERR," using {}-bit tables for long-range coulomb",nctb);
341 else
342 error->message(FLERR," using polynomial approximation for long-range coulomb");
343 }
344
345 if (nctb) {
346 double empirical_precision[17];
347 empirical_precision[6] = 6.99E-03;
348 empirical_precision[7] = 1.78E-03;
349 empirical_precision[8] = 4.72E-04;
350 empirical_precision[9] = 1.17E-04;
351 empirical_precision[10] = 2.95E-05;
352 empirical_precision[11] = 7.41E-06;
353 empirical_precision[12] = 1.76E-06;
354 empirical_precision[13] = 9.28E-07;
355 empirical_precision[14] = 7.46E-07;
356 empirical_precision[15] = 7.32E-07;
357 empirical_precision[16] = 7.30E-07;
358 if (nctb <= 6) table_accuracy = empirical_precision[6];
359 else if (nctb <= 16) table_accuracy = empirical_precision[nctb];
360 else table_accuracy = empirical_precision[16];
361 table_accuracy *= q2_over_sqrt;
362 if ((table_accuracy > spr) && (comm->me == 0))
363 error->warning(FLERR,"For better accuracy use 'pair_modify table 0'");
364 }
365
366 return table_accuracy;
367 }
368
369 /* ----------------------------------------------------------------------
370 convert box coords vector to transposed triclinic lamda (0-1) coords
371 vector, lamda = [(H^-1)^T] * v, does not preserve vector magnitude
372 v and lamda can point to same 3-vector
373 ------------------------------------------------------------------------- */
374
x2lamdaT(double * v,double * lamda)375 void KSpace::x2lamdaT(double *v, double *lamda)
376 {
377 double *h_inv = domain->h_inv;
378 double lamda_tmp[3];
379
380 lamda_tmp[0] = h_inv[0]*v[0];
381 lamda_tmp[1] = h_inv[5]*v[0] + h_inv[1]*v[1];
382 lamda_tmp[2] = h_inv[4]*v[0] + h_inv[3]*v[1] + h_inv[2]*v[2];
383
384 lamda[0] = lamda_tmp[0];
385 lamda[1] = lamda_tmp[1];
386 lamda[2] = lamda_tmp[2];
387 }
388
389 /* ----------------------------------------------------------------------
390 convert lamda (0-1) coords vector to transposed box coords vector
391 lamda = (H^T) * v, does not preserve vector magnitude
392 v and lamda can point to same 3-vector
393 ------------------------------------------------------------------------- */
394
lamda2xT(double * lamda,double * v)395 void KSpace::lamda2xT(double *lamda, double *v)
396 {
397 double h[6];
398 h[0] = domain->h[0];
399 h[1] = domain->h[1];
400 h[2] = domain->h[2];
401 h[3] = fabs(domain->h[3]);
402 h[4] = fabs(domain->h[4]);
403 h[5] = fabs(domain->h[5]);
404 double v_tmp[3];
405
406 v_tmp[0] = h[0]*lamda[0];
407 v_tmp[1] = h[5]*lamda[0] + h[1]*lamda[1];
408 v_tmp[2] = h[4]*lamda[0] + h[3]*lamda[1] + h[2]*lamda[2];
409
410 v[0] = v_tmp[0];
411 v[1] = v_tmp[1];
412 v[2] = v_tmp[2];
413 }
414
415 /* ----------------------------------------------------------------------
416 convert triclinic lamda (0-1) coords vector to box coords vector
417 v = H * lamda, does not preserve vector magnitude
418 lamda and v can point to same 3-vector
419 ------------------------------------------------------------------------- */
420
lamda2xvector(double * lamda,double * v)421 void KSpace::lamda2xvector(double *lamda, double *v)
422 {
423 double *h = domain->h;
424
425 v[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2];
426 v[1] = h[1]*lamda[1] + h[3]*lamda[2];
427 v[2] = h[2]*lamda[2];
428 }
429
430 /* ----------------------------------------------------------------------
431 convert a sphere in box coords to an ellipsoid in lamda (0-1)
432 coords and return the tight (axis-aligned) bounding box, does not
433 preserve vector magnitude see:
434 http://www.loria.fr/~shornus/ellipsoid-bbox.html (no longer online) and
435 https://yiningkarlli.blogspot.com/2013/02/bounding-boxes-for-ellipsoidsfigure.html
436 ------------------------------------------------------------------------- */
437
kspacebbox(double r,double * b)438 void KSpace::kspacebbox(double r, double *b)
439 {
440 double *h = domain->h;
441 double lx,ly,lz,xy,xz,yz;
442 lx = h[0];
443 ly = h[1];
444 lz = h[2];
445 yz = h[3];
446 xz = h[4];
447 xy = h[5];
448
449 b[0] = r*sqrt(ly*ly*lz*lz + ly*ly*xz*xz - 2.0*ly*xy*xz*yz + lz*lz*xy*xy +
450 xy*xy*yz*yz)/(lx*ly*lz);
451 b[1] = r*sqrt(lz*lz + yz*yz)/(ly*lz);
452 b[2] = r/lz;
453 }
454
455 /* ----------------------------------------------------------------------
456 modify parameters of the KSpace style
457 ------------------------------------------------------------------------- */
458
modify_params(int narg,char ** arg)459 void KSpace::modify_params(int narg, char **arg)
460 {
461 int iarg = 0;
462 while (iarg < narg) {
463 if (strcmp(arg[iarg],"mesh") == 0) {
464 if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
465 nx_pppm = nx_msm_max = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
466 ny_pppm = ny_msm_max = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
467 nz_pppm = nz_msm_max = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
468 if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0)
469 gridflag = 0;
470 else if (nx_pppm <= 0 || ny_pppm <= 0 || nz_pppm <= 0)
471 error->all(FLERR,"Kspace_modify mesh parameters must be all "
472 "zero or all positive");
473 else gridflag = 1;
474 iarg += 4;
475 } else if (strcmp(arg[iarg],"mesh/disp") == 0) {
476 if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
477 nx_pppm_6 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
478 ny_pppm_6 = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
479 nz_pppm_6 = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
480 if (nx_pppm_6 == 0 && ny_pppm_6 == 0 && nz_pppm_6 == 0)
481 gridflag_6 = 0;
482 else if (nx_pppm_6 <= 0 || ny_pppm_6 <= 0 || nz_pppm_6 == 0)
483 error->all(FLERR,"Kspace_modify mesh/disp parameters must be all "
484 "zero or all positive");
485 else gridflag_6 = 1;
486 iarg += 4;
487 } else if (strcmp(arg[iarg],"order") == 0) {
488 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
489 order = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
490 iarg += 2;
491 } else if (strcmp(arg[iarg],"order/disp") == 0) {
492 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
493 order_6 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
494 iarg += 2;
495 } else if (strcmp(arg[iarg],"minorder") == 0) {
496 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
497 minorder = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
498 if (minorder < 2) error->all(FLERR,"Illegal kspace_modify command");
499 iarg += 2;
500 } else if (strcmp(arg[iarg],"overlap") == 0) {
501 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
502 if (strcmp(arg[iarg+1],"yes") == 0) overlap_allowed = 1;
503 else if (strcmp(arg[iarg+1],"no") == 0) overlap_allowed = 0;
504 else error->all(FLERR,"Illegal kspace_modify command");
505 iarg += 2;
506 } else if (strcmp(arg[iarg],"force") == 0) {
507 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
508 accuracy_absolute = utils::numeric(FLERR,arg[iarg+1],false,lmp);
509 iarg += 2;
510 } else if (strcmp(arg[iarg],"gewald") == 0) {
511 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
512 g_ewald = utils::numeric(FLERR,arg[iarg+1],false,lmp);
513 if (g_ewald == 0.0) gewaldflag = 0;
514 else gewaldflag = 1;
515 iarg += 2;
516 } else if (strcmp(arg[iarg],"gewald/disp") == 0) {
517 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
518 g_ewald_6 = utils::numeric(FLERR,arg[iarg+1],false,lmp);
519 if (g_ewald_6 == 0.0) gewaldflag_6 = 0;
520 else gewaldflag_6 = 1;
521 iarg += 2;
522 } else if (strcmp(arg[iarg],"slab") == 0) {
523 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
524 if (strcmp(arg[iarg+1],"nozforce") == 0) {
525 slabflag = 2;
526 } else {
527 slabflag = 1;
528 slab_volfactor = utils::numeric(FLERR,arg[iarg+1],false,lmp);
529 if (slab_volfactor <= 1.0)
530 error->all(FLERR,"Bad kspace_modify slab parameter");
531 if (slab_volfactor < 2.0 && comm->me == 0)
532 error->warning(FLERR,"Kspace_modify slab param < 2.0 may "
533 "cause unphysical behavior");
534 }
535 iarg += 2;
536 } else if (strcmp(arg[iarg],"compute") == 0) {
537 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
538 if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1;
539 else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0;
540 else error->all(FLERR,"Illegal kspace_modify command");
541 iarg += 2;
542 } else if (strcmp(arg[iarg],"fftbench") == 0) {
543 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
544 if (strcmp(arg[iarg+1],"yes") == 0) fftbench = 1;
545 else if (strcmp(arg[iarg+1],"no") == 0) fftbench = 0;
546 else error->all(FLERR,"Illegal kspace_modify command");
547 iarg += 2;
548 } else if (strcmp(arg[iarg],"collective") == 0) {
549 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
550 if (strcmp(arg[iarg+1],"yes") == 0) collective_flag = 1;
551 else if (strcmp(arg[iarg+1],"no") == 0) collective_flag = 0;
552 else error->all(FLERR,"Illegal kspace_modify command");
553 iarg += 2;
554 } else if (strcmp(arg[iarg],"diff") == 0) {
555 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
556 if (strcmp(arg[iarg+1],"ad") == 0) differentiation_flag = 1;
557 else if (strcmp(arg[iarg+1],"ik") == 0) differentiation_flag = 0;
558 else error->all(FLERR, "Illegal kspace_modify command");
559 iarg += 2;
560 } else if (strcmp(arg[iarg],"cutoff/adjust") == 0) {
561 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
562 if (strcmp(arg[iarg+1],"yes") == 0) adjust_cutoff_flag = 1;
563 else if (strcmp(arg[iarg+1],"no") == 0) adjust_cutoff_flag = 0;
564 else error->all(FLERR,"Illegal kspace_modify command");
565 iarg += 2;
566 } else if (strcmp(arg[iarg],"kmax/ewald") == 0) {
567 if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command");
568 kx_ewald = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
569 ky_ewald = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
570 kz_ewald = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
571 if (kx_ewald < 0 || ky_ewald < 0 || kz_ewald < 0)
572 error->all(FLERR,"Bad kspace_modify kmax/ewald parameter");
573 if (kx_ewald > 0 && ky_ewald > 0 && kz_ewald > 0)
574 kewaldflag = 1;
575 else
576 kewaldflag = 0;
577 iarg += 4;
578 } else if (strcmp(arg[iarg],"mix/disp") == 0) {
579 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
580 if (strcmp(arg[iarg+1],"pair") == 0) mixflag = 0;
581 else if (strcmp(arg[iarg+1],"geom") == 0) mixflag = 1;
582 else if (strcmp(arg[iarg+1],"none") == 0) mixflag = 2;
583 else error->all(FLERR,"Illegal kspace_modify command");
584 iarg += 2;
585 } else if (strcmp(arg[iarg],"force/disp/real") == 0) {
586 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
587 accuracy_real_6 = utils::numeric(FLERR,arg[iarg+1],false,lmp);
588 iarg += 2;
589 } else if (strcmp(arg[iarg],"force/disp/kspace") == 0) {
590 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
591 accuracy_kspace_6 = utils::numeric(FLERR,arg[iarg+1],false,lmp);
592 iarg += 2;
593 } else if (strcmp(arg[iarg],"eigtol") == 0) {
594 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
595 splittol = utils::numeric(FLERR,arg[iarg+1],false,lmp);
596 if (splittol >= 1.0)
597 error->all(FLERR,"Kspace_modify eigtol must be smaller than one");
598 iarg += 2;
599 } else if (strcmp(arg[iarg],"pressure/scalar") == 0) {
600 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
601 if (strcmp(arg[iarg+1],"yes") == 0) scalar_pressure_flag = 1;
602 else if (strcmp(arg[iarg+1],"no") == 0) scalar_pressure_flag = 0;
603 else error->all(FLERR,"Illegal kspace_modify command");
604 iarg += 2;
605 } else if (strcmp(arg[iarg],"disp/auto") == 0) {
606 if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
607 if (strcmp(arg[iarg+1],"yes") == 0) auto_disp_flag = 1;
608 else if (strcmp(arg[iarg+1],"no") == 0) auto_disp_flag = 0;
609 else error->all(FLERR,"Illegal kspace_modify command");
610 iarg += 2;
611 } else {
612 int n = modify_param(narg-iarg,&arg[iarg]);
613 if (n == 0) error->all(FLERR,"Illegal kspace_modify command");
614 iarg += n;
615 }
616 }
617 }
618
619 /* ---------------------------------------------------------------------- */
620
extract(const char * str)621 void *KSpace::extract(const char *str)
622 {
623 if (strcmp(str,"scale") == 0) return (void *) &scale;
624 return nullptr;
625 }
626