KIM API V2
ex_test_Ar_fcc_cluster_cpp.cpp
Go to the documentation of this file.
1 //
2 //
3 // CDDL HEADER START
4 //
5 // The contents of this file are subject to the terms of the Common Development
6 // and Distribution License Version 1.0 (the "License").
7 //
8 // You can obtain a copy of the license at
9 // http://www.opensource.org/licenses/CDDL-1.0. See the License for the
10 // specific language governing permissions and limitations under the License.
11 //
12 // When distributing Covered Code, include this CDDL HEADER in each file and
13 // include the License file in a prominent location with the name LICENSE.CDDL.
14 // If applicable, add the following below this CDDL HEADER, with the fields
15 // enclosed by brackets "[]" replaced with your own identifying information:
16 //
17 // Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
18 //
19 // CDDL HEADER END
20 //
21 //
22 //
23 // Copyright (c) 2013--2018, Regents of the University of Minnesota.
24 // All rights reserved.
25 //
26 // Contributors:
27 // Ryan S. Elliott
28 // Stephen M. Whalen
29 //
30 //
31 
32 
33 #include <stdlib.h>
34 #include <iomanip>
35 #include <string>
36 #include <cmath>
37 #include "KIM_API.h"
38 #include "KIM_API_status.h"
39 
40 #define NAMESTRLEN 128
41 
42 #define FCCSPACING 5.260
43 #define DIM 3
44 #define NCELLSPERSIDE 2
45 #define NCLUSTERPARTS (4*(NCELLSPERSIDE*NCELLSPERSIDE*NCELLSPERSIDE) + \
46  6*(NCELLSPERSIDE*NCELLSPERSIDE) \
47  + 3*(NCELLSPERSIDE) + 1)
48 
49 #define REPORT_ERROR(LN, FL, MSG, STAT) { \
50  kim_cluster_model.report_error(LN, FL, MSG, STAT); \
51  exit(STAT); \
52  }
53 
54 /* Define neighborlist structure */
55 typedef struct
56 {
57  int iteratorId;
58  int* NNeighbors;
59  int* neighborList;
60 } NeighList;
61 
62 /* Define prototypes */
63 char const * const descriptor();
64 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
65  double cutoff, NeighList* nl);
66 
67 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
68  int* numnei, int** nei1part, double** Rij);
69 
70 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
71 
72 
73 /* Main program */
74 int main()
75 {
76  /* Local variable declarations */
77  double const MinSpacing = 0.8*FCCSPACING;
78  double const MaxSpacing = 1.2*FCCSPACING;
79  double const SpacingIncr = 0.025*FCCSPACING;
80  double CurrentSpacing;
81  double cutpad = 0.75; /* Angstroms */
82  double force_norm;
83  int i;
84  int status;
85 
86 
87  /* model inputs */
88  int numberOfParticles_cluster = NCLUSTERPARTS;
89  int numberOfSpecies = 1;
90  int particleSpecies_cluster_model[NCLUSTERPARTS];
91  double coords_cluster[NCLUSTERPARTS][DIM];
92  NeighList nl_cluster_model;
93  /* model outputs */
94  double cutoff_cluster_model;
95  double energy_cluster_model;
96  double forces_cluster[NCLUSTERPARTS*DIM];
97 
98  std::string modelname;
99 
100  /* Get KIM Model names */
101  printf("Please enter valid KIM Model name: \n");
102  std::cin >> modelname;
103 
104 
105  /* initialize the model */
106  KIM_API_model kim_cluster_model;
107  status = kim_cluster_model.string_init(descriptor(), modelname.c_str());
108  if (KIM_STATUS_OK > status)
109  REPORT_ERROR(__LINE__, __FILE__,"KIM_API_string_init()",status);
110 
111  kim_cluster_model.setm_data(&status, 8*4,
112  "numberOfParticles", 1, &numberOfParticles_cluster, 1,
113  "numberOfSpecies", 1, &numberOfSpecies, 1,
114  "particleSpecies", numberOfParticles_cluster, &particleSpecies_cluster_model, 1,
115  "coordinates", DIM*numberOfParticles_cluster, coords_cluster, 1,
116  "neighObject", 1, &nl_cluster_model, 1,
117  "cutoff", 1, &cutoff_cluster_model, 1,
118  "energy", 1, &energy_cluster_model, 1,
119  "forces", DIM*numberOfParticles_cluster, forces_cluster, 1);
120  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_setm_data",status);
121  status = kim_cluster_model.set_method("get_neigh", 1, (func_ptr) &get_cluster_neigh);
122  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_set_method",status);
123 
124  /* call model init routine */
125  status = kim_cluster_model.model_init();
126  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_model_init", status);
127 
128  /* setup particleSpecies */
129  particleSpecies_cluster_model[0] = kim_cluster_model.get_species_code("Ar", &status);
130  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"get_species_code", status);
131  for (i = 1; i < NCLUSTERPARTS; ++i)
132  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
133 
134  /* setup neighbor lists */
135  /* allocate memory for list */
136  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
137  if (NULL==nl_cluster_model.NNeighbors) REPORT_ERROR(__LINE__, __FILE__,"new unsuccessful", -1);
138 
139  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
140  if (NULL==nl_cluster_model.neighborList) REPORT_ERROR(__LINE__, __FILE__,"new unsuccessful", -1);
141 
142  /* ready to compute */
143  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
144  std::cout << "This is Test : ex_test_Ar_fcc_cluster_cpp\n";
145  std::cout << "--------------------------------------------------------------------------------\n";
146  std::cout << "Results for KIM Model : " << modelname << std::endl;
147 
148  std::cout << std::setw(20) << "Energy"
149  << std::setw(20) << "Force Norm"
150  << std::setw(20) << "Lattice Spacing"
151  << std::endl;
152  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
153  {
154  /* update coordinates for cluster */
155  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
156  /* compute neighbor lists */
157  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
158  (cutoff_cluster_model + cutpad), &nl_cluster_model);
159 
160  /* call compute functions */
161  status = kim_cluster_model.model_compute();
162  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"compute", status);
163 
164  /* compute force norm */
165  force_norm = 0.0;
166  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
167  {
168  force_norm += forces_cluster[i]*forces_cluster[i];
169  }
170  force_norm = sqrt(force_norm);
171 
172  /* print the results */
173  std::cout << std::setw(20) << energy_cluster_model
174  << std::setw(20) << force_norm
175  << std::setw(20) << CurrentSpacing
176  << std::endl;
177  }
178 
179 
180  /* call model destroy */
181  status = kim_cluster_model.model_destroy();
182  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"destroy", status);
183 
184  /* free memory of neighbor lists */
185  delete [] nl_cluster_model.NNeighbors;
186  delete [] nl_cluster_model.neighborList;
187 
188  /* everything is great */
189  return 0;
190 }
191 
192 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
193 {
194  /* local variables */
195  double FCCshifts[4][DIM];
196  double latVec[DIM];
197  int a;
198  int i;
199  int j;
200  int k;
201  int m;
202  int n;
203 
204  /* create a cubic FCC cluster of parts */
205  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
206  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
207  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
208  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
209 
210  a = 0;
211  for (i = 0; i < nCellsPerSide; ++i)
212  {
213  latVec[0] = ((double) i)*FCCspacing;
214  for (j = 0; j < nCellsPerSide; ++j)
215  {
216  latVec[1] = ((double) j)*FCCspacing;
217  for (k = 0; k < nCellsPerSide; ++k)
218  {
219  latVec[2] = ((double) k)*FCCspacing;
220  for (m = 0; m < 4; ++m)
221  {
222  for (n = 0; n < DIM; ++n)
223  {
224  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
225  }
226  a++;
227  }
228  }
229  /* add in the remaining three faces */
230  /* pos-x face */
231  latVec[0] = NCELLSPERSIDE*FCCspacing;
232  latVec[1] = ((double) i)*FCCspacing;
233  latVec[2] = ((double) j)*FCCspacing;
234  for (n = 0; n < DIM; ++n)
235  {
236  coords[a*DIM + n] = latVec[n];
237  }
238  a++;
239  for (n = 0; n < DIM; ++n)
240  {
241  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
242  }
243  a++;
244  /* pos-y face */
245  latVec[0] = ((double) i)*FCCspacing;
246  latVec[1] = NCELLSPERSIDE*FCCspacing;
247  latVec[2] = ((double) j)*FCCspacing;
248  for (n = 0; n < DIM; ++n)
249  {
250  coords[a*DIM + n] = latVec[n];
251  }
252  a++;
253  for (n = 0; n < DIM; ++n)
254  {
255  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
256  }
257  a++;
258  /* pos-z face */
259  latVec[0] = ((double) i)*FCCspacing;
260  latVec[1] = ((double) j)*FCCspacing;
261  latVec[2] = NCELLSPERSIDE*FCCspacing;
262  for (n = 0; n < DIM; ++n)
263  {
264  coords[a*DIM + n] = latVec[n];
265  }
266  a++;
267  for (n = 0; n < DIM; ++n)
268  {
269  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
270  }
271  a++;
272  }
273  /* add in the remaining three edges */
274  latVec[0] = ((double) i)*FCCspacing;
275  latVec[1] = NCELLSPERSIDE*FCCspacing;
276  latVec[2] = NCELLSPERSIDE*FCCspacing;
277  for (n = 0; n < DIM; ++n)
278  {
279  coords[a*DIM + n] = latVec[n];
280  }
281  a++;
282  latVec[0] = NCELLSPERSIDE*FCCspacing;
283  latVec[1] = ((double) i)*FCCspacing;
284  latVec[2] = NCELLSPERSIDE*FCCspacing;
285  for (n = 0; n < DIM; ++n)
286  {
287  coords[a*DIM + n] = latVec[n];
288  }
289  a++;
290  latVec[0] = NCELLSPERSIDE*FCCspacing;
291  latVec[1] = NCELLSPERSIDE*FCCspacing;
292  latVec[2] = ((double) i)*FCCspacing;
293  for (n = 0; n < DIM; ++n)
294  {
295  coords[a*DIM + n] = latVec[n];
296  }
297  a++;
298  }
299  /* add in the remaining corner */
300  for (n = 0; n < DIM; ++n)
301  {
302  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
303  }
304  a++;
305 
306  return;
307 }
308 
309 
310 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
311  double cutoff, NeighList* nl)
312 {
313  /* local variables */
314  int i;
315  int j;
316  int k;
317  int a;
318 
319  double dx[DIM];
320  double r2;
321  double cutoff2;
322 
323  cutoff2 = cutoff*cutoff;
324 
325  for (i = 0; i < numberOfParticles; ++i)
326  {
327  a = 0;
328  for (j = 0; j < numberOfParticles; ++j)
329  {
330  r2 = 0.0;
331  for (k = 0; k < DIM; ++k)
332  {
333  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
334  r2 += dx[k]*dx[k];
335  }
336 
337  if (r2 < cutoff2)
338  {
339  if ((half && i < j) || (!half && i != j))
340  {
341  /* part j is a neighbor of part i */
342  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
343  a++;
344  }
345  }
346  }
347  /* part i has `a' neighbors */
348  (*nl).NNeighbors[i] = a;
349  }
350 
351  return;
352 }
353 
354 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
355  int* numnei, int** nei1part, double** Rij)
356 {
357  /* local variables */
358  KIM_API_model * pkim = *(KIM_API_model **) kimmdl;
359  int partToReturn;
360  int status;
361  int* numberOfParticles;
362  NeighList* nl;
363 
364  /* initialize numnei */
365  *numnei = 0;
366 
367  /* unpack neighbor list object */
368  numberOfParticles = (int*) pkim->get_data("numberOfParticles", &status);
369  if (KIM_STATUS_OK > status)
370  {
371  pkim->report_error(__LINE__, __FILE__,"get_data", status);
372  return status;
373  }
374 
375  nl = (NeighList*) pkim->get_data("neighObject", &status);
376  if (KIM_STATUS_OK > status)
377  {
378  pkim->report_error(__LINE__, __FILE__,"get_data", status);
379  return status;
380  }
381 
382  /* check mode and request */
383  if (0 == *mode) /* iterator mode */
384  {
385  if (0 == *request) /* reset iterator */
386  {
387  (*nl).iteratorId = -1;
388  return KIM_STATUS_NEIGH_ITER_INIT_OK;
389  }
390  else if (1 == *request) /* increment iterator */
391  {
392  (*nl).iteratorId++;
393  if ((*nl).iteratorId >= *numberOfParticles)
394  {
395  return KIM_STATUS_NEIGH_ITER_PAST_END;
396  }
397  else
398  {
399  partToReturn = (*nl).iteratorId;
400  }
401  }
402  else /* invalid request value */
403  {
404  pkim->report_error(__LINE__, __FILE__,"Invalid request in get_cluster_neigh", KIM_STATUS_NEIGH_INVALID_REQUEST);
405  return KIM_STATUS_NEIGH_INVALID_REQUEST;
406  }
407  }
408  else if (1 == *mode) /* locator mode */
409  {
410  if ((*request >= *numberOfParticles) || (*request < 0)) /* invalid id */
411  {
412  pkim->report_error(__LINE__, __FILE__,"Invalid part ID in get_cluster_neigh", KIM_STATUS_PARTICLE_INVALID_ID);
413  return KIM_STATUS_PARTICLE_INVALID_ID;
414  }
415  else
416  {
417  partToReturn = *request;
418  }
419  }
420  else /* invalid mode */
421  {
422  pkim->report_error(__LINE__, __FILE__,"Invalid mode in get_cluster_neigh", KIM_STATUS_NEIGH_INVALID_MODE);
423  return KIM_STATUS_NEIGH_INVALID_MODE;
424  }
425 
426  /* set the returned part */
427  *part = partToReturn;
428 
429  /* set the returned number of neighbors for the returned part */
430  *numnei = (*nl).NNeighbors[*part];
431 
432  /* set the location for the returned neighbor list */
433  *nei1part = &((*nl).neighborList[(*part)*NCLUSTERPARTS]);
434 
435  /* set the pointer to Rij to appropriate value */
436  *Rij = 0;
437 
438  return KIM_STATUS_OK;
439 }
440 
441 char const * const descriptor()
442 {
443  return
444  "KIM_API_Version := 1.6.0\n"
445  "Unit_length := A\n"
446  "Unit_energy := eV\n"
447  "Unit_charge := e\n"
448  "Unit_temperature := K\n"
449  "Unit_time := ps\n"
450  "\n"
451  "PARTICLE_SPECIES:\n"
452  "Ar spec 1\n"
453  "\n"
454  "CONVENTIONS:\n"
455  "ZeroBasedLists flag\n"
456  "Neigh_BothAccess flag\n"
457  "NEIGH_PURE_F flag\n"
458  "\n"
459  "MODEL_INPUT:\n"
460  "numberOfParticles integer none []\n"
461  "numberOfSpecies integer none []\n"
462  "particleSpecies integer none [numberOfParticles]\n"
463  "coordinates double length [numberOfParticles,3]\n"
464  "get_neigh method none []\n"
465  "neighObject pointer none []\n"
466  "\n"
467  "MODEL_OUTPUT:\n"
468  "destroy method none []\n"
469  "compute method none []\n"
470  "reinit method none []\n"
471  "cutoff double length []\n"
472  "energy double energy []\n"
473  "forces double force [numberOfParticles,3]\n";
474 }
#define NCLUSTERPARTS
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
LengthUnit const m
void fcc_cluster_neighborlist(int half, int numberOfParticles, double *coords, double cutoff, NeighList *nl)
#define FCCSPACING
char const *const descriptor()
ComputeArgumentName const numberOfParticles
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
#define REPORT_ERROR(LN, FL, MSG, STAT)
#define NCELLSPERSIDE