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 authors: Albert Bartok (Cambridge University)
16 Aidan Thompson (Sandia, athomps@sandia.gov)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_quip.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "force.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neigh_request.h"
29 #include "neighbor.h"
30 #include "potential_file_reader.h"
31 #include "update.h"
32
33 #include <cmath>
34 #include <cstring>
35
36 using namespace LAMMPS_NS;
37
38 /* ---------------------------------------------------------------------- */
39
PairQUIP(LAMMPS * lmp)40 PairQUIP::PairQUIP(LAMMPS *lmp) : Pair(lmp)
41 {
42 single_enable = 0;
43 restartinfo = 0;
44 one_coeff = 1;
45 no_virial_fdotr_compute = 1;
46 manybody_flag = 1;
47 centroidstressflag = CENTROID_NOTAVAIL;
48 unit_convert_flag = utils::NOCONVERT;
49 map = nullptr;
50 quip_potential = nullptr;
51 quip_file = nullptr;
52 quip_string = nullptr;
53 }
54
~PairQUIP()55 PairQUIP::~PairQUIP()
56 {
57 if (allocated) {
58 memory->destroy(setflag);
59 memory->destroy(cutsq);
60 delete[] map;
61 }
62 delete[] quip_potential;
63 delete[] quip_file;
64 delete[] quip_string;
65 }
66
compute(int eflag,int vflag)67 void PairQUIP::compute(int eflag, int vflag)
68 {
69 int inum, jnum, sum_num_neigh, ii, jj, i, iquip;
70 int *ilist;
71 int *jlist;
72 int *numneigh, **firstneigh;
73 int *quip_num_neigh, *quip_neigh, *atomic_numbers;
74
75 int nlocal = atom->nlocal;
76 int nghost = atom->nghost;
77 int ntotal = nlocal + nghost;
78 int *type = atom->type;
79 tagint *tag = atom->tag;
80
81 double **x = atom->x;
82 double **f = atom->f;
83
84 double *quip_local_e, *quip_force, *quip_local_virial, *quip_virial, quip_energy, *lattice;
85
86 ev_init(eflag, vflag);
87
88 inum = list->inum;
89 ilist = list->ilist;
90 numneigh = list->numneigh;
91 firstneigh = list->firstneigh;
92
93 sum_num_neigh = 0;
94 quip_num_neigh = new int[inum];
95
96 for (ii = 0; ii < inum; ii++) {
97 i = ilist[ii];
98 quip_num_neigh[ii] = numneigh[i];
99 sum_num_neigh += numneigh[i];
100 }
101
102 quip_neigh = new int[sum_num_neigh];
103 iquip = 0;
104
105 for (ii = 0; ii < inum; ii++) {
106 i = ilist[ii];
107 jlist = firstneigh[i];
108 jnum = numneigh[i];
109
110 for (jj = 0; jj < jnum; jj++) {
111 quip_neigh[iquip] = (jlist[jj] & NEIGHMASK) + 1;
112 iquip++;
113 }
114 }
115
116 atomic_numbers = new int[ntotal];
117 for (ii = 0; ii < ntotal; ii++) atomic_numbers[ii] = map[type[ii]];
118
119 quip_local_e = new double[ntotal];
120 quip_force = new double[ntotal * 3];
121 quip_local_virial = new double[ntotal * 9];
122 quip_virial = new double[9];
123
124 lattice = new double[9];
125 lattice[0] = domain->xprd;
126 lattice[1] = 0.0;
127 lattice[2] = 0.0;
128 lattice[3] = domain->xy;
129 lattice[4] = domain->yprd;
130 lattice[5] = 0.0;
131 lattice[6] = domain->xz;
132 lattice[7] = domain->yz;
133 lattice[8] = domain->zprd;
134
135 #if defined(LAMMPS_BIGBIG)
136 int *tmptag = new int[ntotal];
137 int tmplarge = 0, toolarge = 0;
138 for (ii = 0; ii < ntotal; ++ii) {
139 tmptag[ii] = tag[ii];
140 if (tag[ii] > MAXSMALLINT) tmplarge = 1;
141 }
142 MPI_Allreduce(&tmplarge, &toolarge, 1, MPI_INT, MPI_MAX, world);
143 if (toolarge > 0) error->all(FLERR, "Pair style quip does not support 64-bit atom IDs");
144
145 quip_lammps_wrapper(&nlocal, &nghost, atomic_numbers, tmptag, &inum, &sum_num_neigh, ilist,
146 quip_num_neigh, quip_neigh, lattice, quip_potential, &n_quip_potential,
147 &x[0][0], &quip_energy, quip_local_e, quip_virial, quip_local_virial,
148 quip_force);
149
150 delete[] tmptag;
151 #else
152 quip_lammps_wrapper(&nlocal, &nghost, atomic_numbers, tag, &inum, &sum_num_neigh, ilist,
153 quip_num_neigh, quip_neigh, lattice, quip_potential, &n_quip_potential,
154 &x[0][0], &quip_energy, quip_local_e, quip_virial, quip_local_virial,
155 quip_force);
156 #endif
157
158 iquip = 0;
159 for (ii = 0; ii < ntotal; ii++) {
160 for (jj = 0; jj < 3; jj++) {
161 f[ii][jj] += quip_force[iquip];
162 iquip++;
163 }
164 }
165
166 if (eflag_global) { eng_vdwl = quip_energy; }
167
168 if (eflag_atom) {
169 for (ii = 0; ii < ntotal; ii++) { eatom[ii] = quip_local_e[ii]; }
170 }
171
172 if (vflag_global) {
173 virial[0] = quip_virial[0];
174 virial[1] = quip_virial[4];
175 virial[2] = quip_virial[8];
176 virial[3] = (quip_virial[3] + quip_virial[1]) * 0.5;
177 virial[4] = (quip_virial[2] + quip_virial[6]) * 0.5;
178 virial[5] = (quip_virial[5] + quip_virial[7]) * 0.5;
179 }
180
181 if (vflag_atom) {
182 int iatom = 0;
183 for (ii = 0; ii < ntotal; ii++) {
184 vatom[ii][0] += quip_local_virial[iatom + 0];
185 vatom[ii][1] += quip_local_virial[iatom + 4];
186 vatom[ii][2] += quip_local_virial[iatom + 8];
187 vatom[ii][3] += (quip_local_virial[iatom + 3] + quip_local_virial[iatom + 1]) * 0.5;
188 vatom[ii][4] += (quip_local_virial[iatom + 2] + quip_local_virial[iatom + 6]) * 0.5;
189 vatom[ii][5] += (quip_local_virial[iatom + 5] + quip_local_virial[iatom + 7]) * 0.5;
190 iatom += 9;
191 }
192 }
193
194 delete[] atomic_numbers;
195 delete[] quip_num_neigh;
196 delete[] quip_neigh;
197 delete[] quip_local_e;
198 delete[] quip_force;
199 delete[] quip_virial;
200 delete[] quip_local_virial;
201 delete[] lattice;
202 }
203
204 /* ----------------------------------------------------------------------
205 global settings
206 ------------------------------------------------------------------------- */
207
settings(int narg,char **)208 void PairQUIP::settings(int narg, char ** /* arg */)
209 {
210 if (narg != 0) error->all(FLERR, "Illegal pair_style command");
211
212 // check if linked to the correct QUIP library API version
213 // as of 2017-07-19 this is API_VERSION 1
214 if (quip_lammps_api_version() != 1)
215 error->all(FLERR,
216 "QUIP LAMMPS wrapper API version is not compatible "
217 "with this version of LAMMPS");
218
219 // QUIP potentials are parameterized in metal units
220
221 if (strcmp("metal", update->unit_style) != 0)
222 error->all(FLERR, "QUIP potentials require 'metal' units");
223 }
224
225 /* ----------------------------------------------------------------------
226 set coeffs for one or more type pairs
227 ------------------------------------------------------------------------- */
allocate()228 void PairQUIP::allocate()
229 {
230 allocated = 1;
231 int n = atom->ntypes;
232
233 setflag = memory->create(setflag, n + 1, n + 1, "pair:setflag");
234 cutsq = memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
235 map = new int[n + 1];
236 }
237
coeff(int narg,char ** arg)238 void PairQUIP::coeff(int narg, char **arg)
239 {
240 if (!allocated) allocate();
241
242 int n = atom->ntypes;
243 if (narg != (4 + n))
244 error->all(FLERR, "Number of arguments {} is not correct, it should be {}", narg, 4 + n);
245
246 if (comm->me == 0) {
247 PotentialFileReader reader(lmp, arg[2], "QUIP", unit_convert_flag);
248 auto comment = reader.next_string();
249 }
250
251 // use expanded file name, including LAMMPS_POTENTIALS search path
252 quip_file = utils::strdup(utils::get_potential_file_path(arg[2]));
253 quip_string = utils::strdup(arg[3]);
254 n_quip_file = strlen(quip_file);
255 n_quip_string = strlen(quip_string);
256
257 for (int i = 4; i < narg; i++) {
258 if (strcmp(arg[i], "NULL") == 0)
259 map[i - 3] = -1;
260 else
261 map[i - 3] = utils::inumeric(FLERR, arg[i], false, lmp);
262 }
263
264 // clear setflag since coeff() called once with I,J = * *
265
266 n = atom->ntypes;
267 for (int i = 1; i <= n; i++)
268 for (int j = i; j <= n; j++) setflag[i][j] = 0;
269
270 // set setflag i,j for type pairs where both are mapped to elements
271
272 int count = 0;
273 for (int i = 1; i <= n; i++)
274 for (int j = i; j <= n; j++)
275 if (map[i] >= 0 && map[j] >= 0) {
276 setflag[i][j] = 1;
277 count++;
278 }
279
280 if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
281
282 // Initialise potential
283 // First call initializes potential via the fortran code in memory,
284 // and returns the necessary size of quip_potential. This behavior
285 // is invoked by setting n_potential_quip to 0.
286 n_quip_potential = 0;
287 quip_potential = new int[0];
288 quip_lammps_potential_initialise(quip_potential, &n_quip_potential, &cutoff, quip_file,
289 &n_quip_file, quip_string, &n_quip_string);
290 delete[] quip_potential;
291
292 // Allocate quip_potential integer array. This initialise call will transfer
293 // the location of the previously initialised potential to the quip_potential
294 // variable, and we will use it as a handle when calling the actual calculation
295 // routine. We return the cutoff as well.
296 quip_potential = new int[n_quip_potential];
297 quip_lammps_potential_initialise(quip_potential, &n_quip_potential, &cutoff, quip_file,
298 &n_quip_file, quip_string, &n_quip_string);
299 }
300
301 /* ----------------------------------------------------------------------
302 init specific to this pair style
303 ------------------------------------------------------------------------- */
304
init_style()305 void PairQUIP::init_style()
306 {
307 // Require newton pair on
308
309 if (force->newton_pair != 1) error->all(FLERR, "Pair style quip requires newton pair on");
310
311 // Initialise neighbor list
312 int irequest_full = neighbor->request(this);
313 neighbor->requests[irequest_full]->half = 0;
314 neighbor->requests[irequest_full]->full = 1;
315 }
316
317 /* ----------------------------------------------------------------------
318 init for one type pair i,j and corresponding j,i
319 ------------------------------------------------------------------------- */
320
init_one(int,int)321 double PairQUIP::init_one(int /*i*/, int /*j*/)
322 {
323 return cutoff;
324 }
325