1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_KSPACE_H
#define LMP_KSPACE_H
#include "pointers.h" // IWYU pragma: export
#ifdef FFT_SINGLE
typedef float FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_FLOAT
#else
typedef double FFT_SCALAR;
#define MPI_FFT_SCALAR MPI_DOUBLE
#endif
namespace LAMMPS_NS {
class KSpace : protected Pointers {
friend class ThrOMP;
friend class FixOMP;
public:
double energy; // accumulated energies
double energy_1,energy_6;
double virial[6]; // accumulated virial
double *eatom,**vatom; // accumulated per-atom energy/virial
double e2group; // accumulated group-group energy
double f2group[3]; // accumulated group-group force
int triclinic_support; // 1 if supports triclinic geometries
int ewaldflag; // 1 if a Ewald solver
int pppmflag; // 1 if a PPPM solver
int msmflag; // 1 if a MSM solver
int dispersionflag; // 1 if a LJ/dispersion solver
int tip4pflag; // 1 if a TIP4P solver
int dipoleflag; // 1 if a dipole solver
int spinflag; // 1 if a spin solver
int differentiation_flag;
int neighrequest_flag; // used to avoid obsolete construction
// of neighbor lists
int mixflag; // 1 if geometric mixing rules are enforced
// for LJ coefficients
int slabflag;
int scalar_pressure_flag; // 1 if using MSM fast scalar pressure
double slab_volfactor;
int warn_nonneutral; // 0 = error if non-neutral system
// 1 = warn once if non-neutral system
// 2 = warn, but already warned
int warn_nocharge; // 0 = already warned
// 1 = warn if zero charge
int order,order_6,order_allocated;
double accuracy; // accuracy of KSpace solver (force units)
double accuracy_absolute; // user-specified accuracy in force units
double accuracy_relative; // user-specified dimensionless accuracy
// accurary = acc_rel * two_charge_force
double accuracy_real_6; // real space accuracy for
// dispersion solver (force units)
double accuracy_kspace_6; // reciprocal space accuracy for
// dispersion solver (force units)
int auto_disp_flag; // use automatic parameter generation for pppm/disp
double two_charge_force; // force in user units of two point
// charges separated by 1 Angstrom
double g_ewald,g_ewald_6;
int nx_pppm,ny_pppm,nz_pppm; // global FFT grid for Coulombics
int nx_pppm_6,ny_pppm_6,nz_pppm_6; // global FFT grid for dispersion
int nx_msm_max,ny_msm_max,nz_msm_max;
int group_group_enable; // 1 if style supports group/group calculation
// KOKKOS host/device flag and data masks
ExecutionSpace execution_space;
unsigned int datamask_read,datamask_modify;
int copymode;
int compute_flag; // 0 if skip compute()
int fftbench; // 0 if skip FFT timing
int collective_flag; // 1 if use MPI collectives for FFT/remap
int stagger_flag; // 1 if using staggered PPPM grids
double splittol; // tolerance for when to truncate splitting
KSpace(class LAMMPS *);
virtual ~KSpace();
void triclinic_check();
void modify_params(int, char **);
void *extract(const char *);
void compute_dummy(int, int);
// triclinic
void x2lamdaT(double *, double *);
void lamda2xT(double *, double *);
void lamda2xvector(double *, double *);
void kspacebbox(double, double *);
// public so can be called by commands that change charge
void qsum_qsq(int warning_flag = 1);
// general child-class methods
virtual void settings(int, char **) {};
virtual void init() = 0;
virtual void setup() = 0;
virtual void setup_grid() {};
virtual void compute(int, int) = 0;
virtual void compute_group_group(int, int, int) {};
virtual void pack_forward(int, FFT_SCALAR *, int, int *) {};
virtual void unpack_forward(int, FFT_SCALAR *, int, int *) {};
virtual void pack_reverse(int, FFT_SCALAR *, int, int *) {};
virtual void unpack_reverse(int, FFT_SCALAR *, int, int *) {};
virtual int timing(int, double &, double &) {return 0;}
virtual int timing_1d(int, double &) {return 0;}
virtual int timing_3d(int, double &) {return 0;}
virtual int modify_param(int, char **) {return 0;}
virtual double memory_usage() {return 0.0;}
/* ----------------------------------------------------------------------
compute gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164-177
------------------------------------------------------------------------- */
double gamma(const double &rho) const
{
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double g = gcons[split_order][0];
double rho_n = rho2;
for (int n = 1; n <= split_order; n++) {
g += gcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return g;
} else return (1.0/rho);
}
/* ----------------------------------------------------------------------
compute the derivative of gamma for MSM and pair styles
see Eq 4 from Parallel Computing 35 (2009) 164-177
------------------------------------------------------------------------- */
double dgamma(const double &rho) const
{
if (rho <= 1.0) {
const int split_order = order/2;
const double rho2 = rho*rho;
double dg = dgcons[split_order][0]*rho;
double rho_n = rho*rho2;
for (int n = 1; n < split_order; n++) {
dg += dgcons[split_order][n]*rho_n;
rho_n *= rho2;
}
return dg;
} else return (-1.0/rho/rho);
}
double **get_gcons() { return gcons; }
double **get_dgcons() { return dgcons; }
protected:
int gridflag,gridflag_6;
int gewaldflag,gewaldflag_6;
int minorder,overlap_allowed;
int adjust_cutoff_flag;
int suffix_flag; // suffix compatibility flag
bigint natoms_original;
double scale,qqrd2e;
double qsum,qsqsum,q2;
double **gcons,**dgcons; // accumulated per-atom energy/virial
int evflag,evflag_atom;
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int maxeatom,maxvatom;
int kewaldflag; // 1 if kspace range set for Ewald sum
int kx_ewald,ky_ewald,kz_ewald; // kspace settings for Ewald sum
void pair_check();
void ev_init(int eflag, int vflag, int alloc = 1) {
if (eflag||vflag) ev_setup(eflag, vflag, alloc);
else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = 0;
}
void ev_setup(int, int, int alloc = 1);
double estimate_table_accuracy(double, double);
};
}
#endif
/* ERROR/WARNING messages:
E: KSpace style does not yet support triclinic geometries
The specified kspace style does not allow for non-orthogonal
simulation boxes.
E: KSpace solver requires a pair style
No pair style is defined.
E: KSpace style is incompatible with Pair style
Setting a kspace style requires that a pair style with matching
long-range Coulombic or dispersion components be used.
W: Using kspace solver on system with no charge
Self-explanatory.
E: System is not charge neutral, net charge = %g
The total charge on all atoms on the system is not 0.0.
For some KSpace solvers this is an error.
W: System is not charge neutral, net charge = %g
The total charge on all atoms on the system is not 0.0.
For some KSpace solvers this is only a warning.
W: For better accuracy use 'pair_modify table 0'
The user-specified force accuracy cannot be achieved unless the table
feature is disabled by using 'pair_modify table 0'.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Bad kspace_modify slab parameter
Kspace_modify value for the slab/volume keyword must be >= 2.0.
E: Kspace_modify mesh parameter must be all zero or all positive
Valid kspace mesh parameters are >0. The code will try to auto-detect
suitable values when all three mesh sizes are set to zero (the default).
E: Kspace_modify mesh/disp parameter must be all zero or all positive
Valid kspace mesh/disp parameters are >0. The code will try to auto-detect
suitable values when all three mesh sizes are set to zero [and]
the required accuracy via {force/disp/real} as well as
{force/disp/kspace} is set.
W: Kspace_modify slab param < 2.0 may cause unphysical behavior
The kspace_modify slab parameter should be larger to insure periodic
grids padded with empty space do not overlap.
E: Bad kspace_modify kmax/ewald parameter
Kspace_modify values for the kmax/ewald keyword must be integers > 0
E: Kspace_modify eigtol must be smaller than one
Self-explanatory.
*/