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 // C or Fortran style library interface to LAMMPS
47 // customize by adding new LAMMPS-specific functions
48
49 #include "lmptype.h"
50 #include <mpi.h>
51 #include <string.h>
52 #include <stdlib.h>
53 #include "library.h"
54 #include "lammps.h"
55 #include "input.h"
56 #include "atom.h"
57 #include "domain.h"
58 #include "update.h"
59 #include "group.h"
60 #include "input.h"
61 #include "variable.h"
62 #include "modify.h"
63 #include "compute.h"
64 #include "fix.h"
65 #include "comm.h"
66 #include "memory.h"
67 #include "error.h"
68
69 using namespace LAMMPS_NS;
70
71 /* ----------------------------------------------------------------------
72 create an instance of LAMMPS and return pointer to it
73 pass in command-line args and MPI communicator to run on
74 ------------------------------------------------------------------------- */
75
lammps_open(int argc,char ** argv,MPI_Comm communicator,void ** ptr)76 void lammps_open(int argc, char **argv, MPI_Comm communicator, void **ptr)
77 {
78 LAMMPS *lmp = new LAMMPS(argc,argv,communicator);
79 *ptr = (void *) lmp;
80 }
81
82 /* ----------------------------------------------------------------------
83 create an instance of LAMMPS and return pointer to it
84 caller doesn't know MPI communicator, so use MPI_COMM_WORLD
85 intialize MPI if needed
86 ------------------------------------------------------------------------- */
87
lammps_open_no_mpi(int argc,char ** argv,void ** ptr)88 void lammps_open_no_mpi(int argc, char **argv, void **ptr)
89 {
90 int flag;
91 MPI_Initialized(&flag);
92
93 if (!flag) {
94 int argc = 0;
95 char **argv = NULL;
96 MPI_Init(&argc,&argv);
97 }
98
99 MPI_Comm communicator = MPI_COMM_WORLD;
100
101 LAMMPS *lmp = new LAMMPS(argc,argv,communicator);
102 *ptr = (void *) lmp;
103 }
104
105 /* ----------------------------------------------------------------------
106 destruct an instance of LAMMPS
107 ------------------------------------------------------------------------- */
108
lammps_close(void * ptr)109 void lammps_close(void *ptr)
110 {
111 LAMMPS *lmp = (LAMMPS *) ptr;
112 delete lmp;
113 }
114
115 /* ----------------------------------------------------------------------
116 process an input script in filename str
117 ------------------------------------------------------------------------- */
118
lammps_file(void * ptr,const char * str)119 void lammps_file(void *ptr, const char *str)
120 {
121 LAMMPS *lmp = (LAMMPS *) ptr;
122 lmp->input->file(str);
123 }
124
125 /* ----------------------------------------------------------------------
126 process a single input command in str
127 ------------------------------------------------------------------------- */
128
lammps_command(void * ptr,const char * str)129 char *lammps_command(void *ptr, const char *str)
130 {
131 LAMMPS *lmp = (LAMMPS *) ptr;
132 return lmp->input->one(str);
133 }
134
135 /* ----------------------------------------------------------------------
136 clean-up function to free memory allocated by lib and returned to caller
137 ------------------------------------------------------------------------- */
138
lammps_free(void * ptr)139 void lammps_free(void *ptr)
140 {
141 free(ptr);
142 }
143
144 /* ----------------------------------------------------------------------
145 add LAMMPS-specific library functions
146 all must receive LAMMPS pointer as argument
147 customize by adding a function here and in library.h header file
148 ------------------------------------------------------------------------- */
149
150 /* ----------------------------------------------------------------------
151 extract a pointer to an internal LAMMPS global entity
152 name = desired quantity, e.g. dt or boxyhi or natoms
153 returns a void pointer to the entity
154 which the caller can cast to the proper data type
155 returns a NULL if name not listed below
156 customize by adding names
157 ------------------------------------------------------------------------- */
158
lammps_extract_global(void * ptr,const char * name)159 void *lammps_extract_global(void *ptr, const char *name)
160 {
161 LAMMPS *lmp = (LAMMPS *) ptr;
162
163 if (strcmp(name,"dt") == 0) return (void *) &lmp->update->dt;
164 if (strcmp(name,"atime") == 0) return (void *)&lmp->update->atime;
165 if (strcmp(name,"boxxlo") == 0) return (void *) &lmp->domain->boxlo[0];
166 if (strcmp(name,"boxxhi") == 0) return (void *) &lmp->domain->boxhi[0];
167 if (strcmp(name,"boxylo") == 0) return (void *) &lmp->domain->boxlo[1];
168 if (strcmp(name,"boxyhi") == 0) return (void *) &lmp->domain->boxhi[1];
169 if (strcmp(name,"boxzlo") == 0) return (void *) &lmp->domain->boxlo[2];
170 if (strcmp(name,"boxzhi") == 0) return (void *) &lmp->domain->boxhi[2];
171 if (strcmp(name,"subxlo") == 0) return (void *) &lmp->domain->sublo[0];
172 if (strcmp(name,"subxhi") == 0) return (void *) &lmp->domain->subhi[0];
173 if (strcmp(name,"subylo") == 0) return (void *) &lmp->domain->sublo[1];
174 if (strcmp(name,"subyhi") == 0) return (void *) &lmp->domain->subhi[1];
175 if (strcmp(name,"subzlo") == 0) return (void *) &lmp->domain->sublo[2];
176 if (strcmp(name,"subzhi") == 0) return (void *) &lmp->domain->subhi[2];
177 if (strcmp(name,"procx") == 0) return (void *) &lmp->comm->procgrid[0];
178 if (strcmp(name,"procy") == 0) return (void *) &lmp->comm->procgrid[1];
179 if (strcmp(name,"procz") == 0) return (void *) &lmp->comm->procgrid[2];
180 if (strcmp(name,"procneighxleft") == 0) return (void *) &lmp->comm->procneigh[0][0];
181 if (strcmp(name,"procneighxright") == 0) return (void *) &lmp->comm->procneigh[0][1];
182 if (strcmp(name,"procneighyleft") == 0) return (void *) &lmp->comm->procneigh[1][0];
183 if (strcmp(name,"procneighyright") == 0) return (void *) &lmp->comm->procneigh[1][1];
184 if (strcmp(name,"procneighzleft") == 0) return (void *) &lmp->comm->procneigh[2][0];
185 if (strcmp(name,"procneighzright") == 0) return (void *) &lmp->comm->procneigh[2][1];
186 if (strcmp(name,"mylocx") == 0) return (void *) &lmp->comm->myloc[0];
187 if (strcmp(name,"mylocy") == 0) return (void *) &lmp->comm->myloc[1];
188 if (strcmp(name,"mylocz") == 0) return (void *) &lmp->comm->myloc[2];
189 if (strcmp(name,"natoms") == 0) return (void *) &lmp->atom->natoms;
190 if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal;
191 if (strcmp(name,"nghost") == 0) return (void *) &lmp->atom->nghost;
192 if (strcmp(name,"ago") == 0) return (void *) &lmp->neighbor->ago; //INT, time steps since last neighbor->decide,
193 if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep;
194 if (strcmp(name,"firststep") == 0) return (void *) &lmp->update->firststep;
195 return NULL;
196 }
197
198 /* ----------------------------------------------------------------------
199 extract a pointer to an internal LAMMPS atom-based entity
200 name = desired quantity, e.g. x or mass
201 returns a void pointer to the entity
202 which the caller can cast to the proper data type
203 returns a NULL if Atom::extract() does not recognize the name
204 customize by adding names to Atom::extract()
205 ------------------------------------------------------------------------- */
206
lammps_extract_atom(void * ptr,const char * name)207 void *lammps_extract_atom(void *ptr, const char *name)
208 {
209 LAMMPS *lmp = (LAMMPS *) ptr;
210 return lmp->atom->extract(name);
211 }
212
213 /* ----------------------------------------------------------------------
214 extract a pointer to an internal LAMMPS compute-based entity
215 id = compute ID
216 style = 0 for global data, 1 for per-atom data, 2 for local data
217 type = 0 for scalar, 1 for vector, 2 for array
218 for global data, returns a pointer to the
219 compute's internal data structure for the entity
220 caller should cast it to (double *) for a scalar or vector
221 caller should cast it to (double **) for an array
222 for per-atom or local data, returns a pointer to the
223 compute's internal data structure for the entity
224 caller should cast it to (double *) for a vector
225 caller should cast it to (double **) for an array
226 returns a void pointer to the compute's internal data structure
227 for the entity which the caller can cast to the proper data type
228 returns a NULL if id is not recognized or style/type not supported
229 IMPORTANT: if the compute is not current it will be invoked
230 LAMMPS cannot easily check here if it is valid to invoke the compute,
231 so caller must insure that it is OK
232 ------------------------------------------------------------------------- */
233
lammps_extract_compute(void * ptr,const char * id,int style,int type)234 void *lammps_extract_compute(void *ptr, const char *id, int style, int type)
235 {
236 LAMMPS *lmp = (LAMMPS *) ptr;
237
238 int icompute = lmp->modify->find_compute(id);
239 if (icompute < 0) return NULL;
240 Compute *compute = lmp->modify->compute[icompute];
241
242 if (style == 0) {
243 if (type == 0) {
244 if (!compute->scalar_flag) return NULL;
245 if (compute->invoked_scalar != lmp->update->ntimestep)
246 compute->compute_scalar();
247 return (void *) &compute->scalar;
248 }
249 if (type == 1) {
250 if (!compute->vector_flag) return NULL;
251 if (compute->invoked_vector != lmp->update->ntimestep)
252 compute->compute_vector();
253 return (void *) compute->vector;
254 }
255 if (type == 2) {
256 if (!compute->array_flag) return NULL;
257 if (compute->invoked_array != lmp->update->ntimestep)
258 compute->compute_array();
259 return (void *) compute->array;
260 }
261 }
262
263 if (style == 1) {
264 if (!compute->peratom_flag) return NULL;
265 if (type == 1) {
266 if (compute->invoked_peratom != lmp->update->ntimestep)
267 compute->compute_peratom();
268 return (void *) compute->vector_atom;
269 }
270 if (type == 2) {
271 if (compute->invoked_peratom != lmp->update->ntimestep)
272 compute->compute_peratom();
273 return (void *) compute->array_atom;
274 }
275 }
276
277 if (style == 2) {
278 if (!compute->local_flag) return NULL;
279 if (type == 1) {
280 if (compute->invoked_local != lmp->update->ntimestep)
281 compute->compute_local();
282 return (void *) compute->vector_local;
283 }
284 if (type == 2) {
285 if (compute->invoked_local != lmp->update->ntimestep)
286 compute->compute_local();
287 return (void *) compute->array_local;
288 }
289 }
290
291 return NULL;
292 }
293
294 /* ----------------------------------------------------------------------
295 extract a pointer to an internal LAMMPS fix-based entity
296 id = fix ID
297 style = 0 for global data, 1 for per-atom data, 2 for local data
298 type = 0 for scalar, 1 for vector, 2 for array
299 i,j = indices needed only to specify which global vector or array value
300 for global data, returns a pointer to a memory location
301 which is allocated by this function
302 which the caller can cast to a (double *) which points to the value
303 for per-atom or local data, returns a pointer to the
304 fix's internal data structure for the entity
305 caller should cast it to (double *) for a vector
306 caller should cast it to (double **) for an array
307 returns a NULL if id is not recognized or style/type not supported
308 IMPORTANT: for global data,
309 this function allocates a double to store the value in,
310 so the caller must free this memory to avoid a leak, e.g.
311 double *dptr = (double *) lammps_extract_fix();
312 double value = *dptr;
313 free(dptr);
314 IMPORTANT: LAMMPS cannot easily check here when info extracted from
315 the fix is valid, so caller must insure that it is OK
316 ------------------------------------------------------------------------- */
317
lammps_extract_fix(void * ptr,const char * id,int style,int type,int i,int j)318 void *lammps_extract_fix(void *ptr, const char *id, int style, int type,
319 int i, int j)
320 {
321 LAMMPS *lmp = (LAMMPS *) ptr;
322
323 int ifix = lmp->modify->find_fix(id);
324 if (ifix < 0) return NULL;
325 Fix *fix = lmp->modify->fix[ifix];
326
327 if (style == 0) {
328 double *dptr = (double *) malloc(sizeof(double));
329 if (type == 0) {
330 if (!fix->scalar_flag) return NULL;
331 *dptr = fix->compute_scalar();
332 return (void *) dptr;
333 }
334 if (type == 1) {
335 if (!fix->vector_flag) return NULL;
336 *dptr = fix->compute_vector(i);
337 return (void *) dptr;
338 }
339 if (type == 2) {
340 if (!fix->array_flag) return NULL;
341 *dptr = fix->compute_array(i,j);
342 return (void *) dptr;
343 }
344 }
345
346 if (style == 1) {
347 if (!fix->peratom_flag) return NULL;
348 if (type == 1) return (void *) fix->vector_atom;
349 if (type == 2) return (void *) fix->array_atom;
350 }
351
352 if (style == 2) {
353 if (!fix->local_flag) return NULL;
354 if (type == 1) return (void *) fix->vector_local;
355 if (type == 2) return (void *) fix->array_local;
356 }
357
358 return NULL;
359 }
360
361 /* ----------------------------------------------------------------------
362 extract a pointer to an internal LAMMPS evaluated variable
363 name = variable name, must be equal-style or atom-style variable
364 group = group ID for evaluating an atom-style variable, else NULL
365 for equal-style variable, returns a pointer to a memory location
366 which is allocated by this function
367 which the caller can cast to a (double *) which points to the value
368 for atom-style variable, returns a pointer to the
369 vector of per-atom values on each processor,
370 which the caller can cast to a (double *) which points to the values
371 returns a NULL if name is not recognized or not equal-style or atom-style
372 IMPORTANT: for both equal-style and atom-style variables,
373 this function allocates memory to store the variable data in
374 so the caller must free this memory to avoid a leak
375 e.g. for equal-style variables
376 double *dptr = (double *) lammps_extract_variable();
377 double value = *dptr;
378 free(dptr);
379 e.g. for atom-style variables
380 double *vector = (double *) lammps_extract_variable();
381 use the vector values
382 free(vector);
383 IMPORTANT: LAMMPS cannot easily check here when it is valid to evaluate
384 the variable or any fixes or computes or thermodynamic info it references,
385 so caller must insure that it is OK
386 ------------------------------------------------------------------------- */
387
lammps_extract_variable(void * ptr,char * name,char * group)388 void *lammps_extract_variable(void *ptr, char *name, char *group)
389 {
390 LAMMPS *lmp = (LAMMPS *) ptr;
391
392 int ivar = lmp->input->variable->find(name);
393 if (ivar < 0) return NULL;
394
395 if (lmp->input->variable->equalstyle(ivar)) {
396 double *dptr = (double *) malloc(sizeof(double));
397 *dptr = lmp->input->variable->compute_equal(ivar);
398 return (void *) dptr;
399 }
400
401 if (lmp->input->variable->atomstyle(ivar)) {
402 int igroup = lmp->group->find(group);
403 if (igroup < 0) return NULL;
404 int nlocal = lmp->atom->nlocal;
405 double *vector = (double *) malloc(nlocal*sizeof(double));
406 lmp->input->variable->compute_atom(ivar,igroup,vector,1,0);
407 return (void *) vector;
408 }
409
410 return NULL;
411 }
412
413 /* ----------------------------------------------------------------------
414 return the total number of atoms in the system
415 useful before call to lammps_get_atoms() so can pre-allocate vector
416 ------------------------------------------------------------------------- */
417
lammps_get_natoms(void * ptr)418 int lammps_get_natoms(void *ptr)
419 {
420 LAMMPS *lmp = (LAMMPS *) ptr;
421 if (lmp->atom->natoms > MAXSMALLINT) return 0;
422 int natoms = static_cast<int> (lmp->atom->natoms);
423 return natoms;
424 }
425
426 /* ----------------------------------------------------------------------
427 gather the named atom-based entity across all processors
428 name = desired quantity, e.g. x or charge
429 type = 0 for integer values, 1 for double values
430 count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
431 return atom-based values in data, ordered by count, then by atom ID
432 e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
433 data must be pre-allocated by caller to correct length
434 ------------------------------------------------------------------------- */
435
lammps_gather_atoms(void * ptr,const char * name,int type,int count,void * data)436 void lammps_gather_atoms(void *ptr, const char *name,
437 int type, int count, void *data)
438 {
439 LAMMPS *lmp = (LAMMPS *) ptr;
440
441 // error if tags are not defined or not consecutive
442
443 int flag = 0;
444 if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
445 if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
446 if (flag && lmp->comm->me == 0) {
447 lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
448 return;
449 }
450
451 int natoms = static_cast<int> (lmp->atom->natoms);
452
453 int i,j,offset;
454 void *vptr = lmp->atom->extract(name);
455
456 // copy = Natom length vector of per-atom values
457 // use atom ID to insert each atom's values into copy
458 // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
459
460 if (type == 0) {
461 int *vector = NULL;
462 int **array = NULL;
463 if (count == 1) vector = (int *) vptr;
464 else array = (int **) vptr;
465
466 int *copy;
467 lmp->memory->create(copy,count*natoms,"lib/gather:copy");
468 for (i = 0; i < count*natoms; i++) copy[i] = 0;
469
470 int *tag = lmp->atom->tag;
471 int nlocal = lmp->atom->nlocal;
472
473 if (count == 1)
474 for (i = 0; i < nlocal; i++)
475 copy[tag[i]-1] = vector[i];
476 else
477 for (i = 0; i < nlocal; i++) {
478 offset = count*(tag[i]-1);
479 for (j = 0; j < count; j++)
480 copy[offset++] = array[i][0];
481 }
482
483 MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
484 lmp->memory->destroy(copy);
485
486 } else {
487 double *vector = NULL;
488 double **array = NULL;
489 if (count == 1) vector = (double *) vptr;
490 else array = (double **) vptr;
491
492 double *copy;
493 lmp->memory->create(copy,count*natoms,"lib/gather:copy");
494 for (i = 0; i < count*natoms; i++) copy[i] = 0.0;
495
496 int *tag = lmp->atom->tag;
497 int nlocal = lmp->atom->nlocal;
498
499 if (count == 1) {
500 for (i = 0; i < nlocal; i++)
501 copy[tag[i]-1] = vector[i];
502 } else {
503 for (i = 0; i < nlocal; i++) {
504 offset = count*(tag[i]-1);
505 for (j = 0; j < count; j++)
506 copy[offset++] = array[i][j];
507 }
508 }
509
510 MPI_Allreduce(copy,data,count*natoms,MPI_DOUBLE,MPI_SUM,lmp->world);
511 lmp->memory->destroy(copy);
512 }
513 }
514
515 /* ----------------------------------------------------------------------
516 scatter the named atom-based entity across all processors
517 name = desired quantity, e.g. x or charge
518 type = 0 for integer values, 1 for double values
519 count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
520 data = atom-based values in data, ordered by count, then by atom ID
521 e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
522 ------------------------------------------------------------------------- */
523
lammps_scatter_atoms(void * ptr,const char * name,int type,int count,void * data)524 void lammps_scatter_atoms(void *ptr, const char *name,
525 int type, int count, void *data)
526 {
527 LAMMPS *lmp = (LAMMPS *) ptr;
528
529 // error if tags are not defined or not consecutive or no atom map
530
531 int flag = 0;
532 if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
533 if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
534 if (lmp->atom->map_style == 0) flag = 1;
535 if (flag && lmp->comm->me == 0) {
536 lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms");
537 return;
538 }
539
540 int natoms = static_cast<int> (lmp->atom->natoms);
541
542 int i,j,m,offset;
543 void *vptr = lmp->atom->extract(name);
544
545 // copy = Natom length vector of per-atom values
546 // use atom ID to insert each atom's values into copy
547 // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
548
549 if (type == 0) {
550 int *vector = NULL;
551 int **array = NULL;
552 if (count == 1) vector = (int *) vptr;
553 else array = (int **) vptr;
554 int *dptr = (int *) data;
555
556 if (count == 1) {
557 for (i = 0; i < natoms; i++)
558 if ((m = lmp->atom->map(i+1)) >= 0)
559 vector[m] = dptr[i];
560 } else {
561 for (i = 0; i < natoms; i++)
562 if ((m = lmp->atom->map(i+1)) >= 0) {
563 offset = count*i;
564 for (j = 0; j < count; j++)
565 array[m][j] = dptr[offset++];
566 }
567 }
568 } else {
569 double *vector = NULL;
570 double **array = NULL;
571 if (count == 1) vector = (double *) vptr;
572 else array = (double **) vptr;
573 double *dptr = (double *) data;
574
575 if (count == 1) {
576 for (i = 0; i < natoms; i++)
577 if ((m = lmp->atom->map(i+1)) >= 0)
578 vector[m] = dptr[i];
579 } else {
580 for (i = 0; i < natoms; i++) {
581 if ((m = lmp->atom->map(i+1)) >= 0) {
582 offset = count*i;
583 for (j = 0; j < count; j++)
584 array[m][j] = dptr[offset++];
585 }
586 }
587 }
588 }
589 }
590