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* NNeighbors;
58  int* neighborList;
59 } NeighList;
60 
61 /* Define prototypes */
62 char const * const descriptor();
63 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
64  double cutoff, NeighList* nl);
65 
66 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
67  int* numnei, int** nei1part, double** Rij);
68 
69 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
70 
71 
72 /* Main program */
73 int main()
74 {
75  /* Local variable declarations */
76  double const MinSpacing = 0.8*FCCSPACING;
77  double const MaxSpacing = 1.2*FCCSPACING;
78  double const SpacingIncr = 0.025*FCCSPACING;
79  double CurrentSpacing;
80  double cutpad = 0.75; /* Angstroms */
81  double force_norm;
82  int i;
83  int status;
84 
85 
86  /* model inputs */
87  int numberOfParticles_cluster = NCLUSTERPARTS;
88  int numberOfSpecies = 1;
89  int particleSpecies_cluster_model[NCLUSTERPARTS];
90  double coords_cluster[NCLUSTERPARTS][DIM];
91  NeighList nl_cluster_model;
92  /* model outputs */
93  double cutoff_cluster_model;
94  double energy_cluster_model;
95  double forces_cluster[NCLUSTERPARTS*DIM];
96 
97  std::string modelname;
98 
99  /* Get KIM Model names */
100  printf("Please enter valid KIM Model name: \n");
101  std::cin >> modelname;
102 
103 
104  /* initialize the model */
105  KIM_API_model kim_cluster_model;
106  status = kim_cluster_model.string_init(descriptor(), modelname.c_str());
107  if (KIM_STATUS_OK > status)
108  REPORT_ERROR(__LINE__, __FILE__,"KIM_API_string_init()",status);
109 
110  kim_cluster_model.setm_data(&status, 8*4,
111  "numberOfParticles", 1, &numberOfParticles_cluster, 1,
112  "numberOfSpecies", 1, &numberOfSpecies, 1,
113  "particleSpecies", numberOfParticles_cluster, &particleSpecies_cluster_model, 1,
114  "coordinates", DIM*numberOfParticles_cluster, coords_cluster, 1,
115  "neighObject", 1, &nl_cluster_model, 1,
116  "cutoff", 1, &cutoff_cluster_model, 1,
117  "energy", 1, &energy_cluster_model, 1,
118  "forces", DIM*numberOfParticles_cluster, forces_cluster, 1);
119  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_setm_data",status);
120  status = kim_cluster_model.set_method("get_neigh", 1, (func_ptr) &get_cluster_neigh);
121  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_set_method",status);
122 
123  /* call model init routine */
124  status = kim_cluster_model.model_init();
125  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_model_init", status);
126 
127  /* setup particleSpecies */
128  particleSpecies_cluster_model[0] = kim_cluster_model.get_species_code("Ar", &status);
129  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"get_species_code", status);
130  for (i = 1; i < NCLUSTERPARTS; ++i)
131  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
132 
133  /* setup neighbor lists */
134  /* allocate memory for list */
135  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
136  if (NULL==nl_cluster_model.NNeighbors) REPORT_ERROR(__LINE__, __FILE__,"new unsuccessful", -1);
137 
138  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
139  if (NULL==nl_cluster_model.neighborList) REPORT_ERROR(__LINE__, __FILE__,"new unsuccessful", -1);
140 
141  /* ready to compute */
142  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
143  std::cout << "This is Test : ex_test_Ar_fcc_cluster_cpp\n";
144  std::cout << "--------------------------------------------------------------------------------\n";
145  std::cout << "Results for KIM Model : " << modelname << std::endl;
146 
147  std::cout << std::setw(20) << "Energy"
148  << std::setw(20) << "Force Norm"
149  << std::setw(20) << "Lattice Spacing"
150  << std::endl;
151  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
152  {
153  /* update coordinates for cluster */
154  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
155  /* compute neighbor lists */
156  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
157  (cutoff_cluster_model + cutpad), &nl_cluster_model);
158 
159  /* call compute functions */
160  status = kim_cluster_model.model_compute();
161  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"compute", status);
162 
163  /* compute force norm */
164  force_norm = 0.0;
165  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
166  {
167  force_norm += forces_cluster[i]*forces_cluster[i];
168  }
169  force_norm = sqrt(force_norm);
170 
171  /* print the results */
172  std::cout << std::setw(20) << energy_cluster_model
173  << std::setw(20) << force_norm
174  << std::setw(20) << CurrentSpacing
175  << std::endl;
176  }
177 
178 
179  /* call model destroy */
180  status = kim_cluster_model.model_destroy();
181  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"destroy", status);
182 
183  /* free memory of neighbor lists */
184  delete [] nl_cluster_model.NNeighbors;
185  delete [] nl_cluster_model.neighborList;
186 
187  /* everything is great */
188  return 0;
189 }
190 
191 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
192 {
193  /* local variables */
194  double FCCshifts[4][DIM];
195  double latVec[DIM];
196  int a;
197  int i;
198  int j;
199  int k;
200  int m;
201  int n;
202 
203  /* create a cubic FCC cluster of parts */
204  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
205  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
206  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
207  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
208 
209  a = 0;
210  for (i = 0; i < nCellsPerSide; ++i)
211  {
212  latVec[0] = ((double) i)*FCCspacing;
213  for (j = 0; j < nCellsPerSide; ++j)
214  {
215  latVec[1] = ((double) j)*FCCspacing;
216  for (k = 0; k < nCellsPerSide; ++k)
217  {
218  latVec[2] = ((double) k)*FCCspacing;
219  for (m = 0; m < 4; ++m)
220  {
221  for (n = 0; n < DIM; ++n)
222  {
223  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
224  }
225  a++;
226  }
227  }
228  /* add in the remaining three faces */
229  /* pos-x face */
230  latVec[0] = NCELLSPERSIDE*FCCspacing;
231  latVec[1] = ((double) i)*FCCspacing;
232  latVec[2] = ((double) j)*FCCspacing;
233  for (n = 0; n < DIM; ++n)
234  {
235  coords[a*DIM + n] = latVec[n];
236  }
237  a++;
238  for (n = 0; n < DIM; ++n)
239  {
240  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
241  }
242  a++;
243  /* pos-y face */
244  latVec[0] = ((double) i)*FCCspacing;
245  latVec[1] = NCELLSPERSIDE*FCCspacing;
246  latVec[2] = ((double) j)*FCCspacing;
247  for (n = 0; n < DIM; ++n)
248  {
249  coords[a*DIM + n] = latVec[n];
250  }
251  a++;
252  for (n = 0; n < DIM; ++n)
253  {
254  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
255  }
256  a++;
257  /* pos-z face */
258  latVec[0] = ((double) i)*FCCspacing;
259  latVec[1] = ((double) j)*FCCspacing;
260  latVec[2] = NCELLSPERSIDE*FCCspacing;
261  for (n = 0; n < DIM; ++n)
262  {
263  coords[a*DIM + n] = latVec[n];
264  }
265  a++;
266  for (n = 0; n < DIM; ++n)
267  {
268  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
269  }
270  a++;
271  }
272  /* add in the remaining three edges */
273  latVec[0] = ((double) i)*FCCspacing;
274  latVec[1] = NCELLSPERSIDE*FCCspacing;
275  latVec[2] = NCELLSPERSIDE*FCCspacing;
276  for (n = 0; n < DIM; ++n)
277  {
278  coords[a*DIM + n] = latVec[n];
279  }
280  a++;
281  latVec[0] = NCELLSPERSIDE*FCCspacing;
282  latVec[1] = ((double) i)*FCCspacing;
283  latVec[2] = NCELLSPERSIDE*FCCspacing;
284  for (n = 0; n < DIM; ++n)
285  {
286  coords[a*DIM + n] = latVec[n];
287  }
288  a++;
289  latVec[0] = NCELLSPERSIDE*FCCspacing;
290  latVec[1] = NCELLSPERSIDE*FCCspacing;
291  latVec[2] = ((double) i)*FCCspacing;
292  for (n = 0; n < DIM; ++n)
293  {
294  coords[a*DIM + n] = latVec[n];
295  }
296  a++;
297  }
298  /* add in the remaining corner */
299  for (n = 0; n < DIM; ++n)
300  {
301  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
302  }
303  a++;
304 
305  return;
306 }
307 
308 
309 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
310  double cutoff, NeighList* nl)
311 {
312  /* local variables */
313  int i;
314  int j;
315  int k;
316  int a;
317 
318  double dx[DIM];
319  double r2;
320  double cutoff2;
321 
322  cutoff2 = cutoff*cutoff;
323 
324  for (i = 0; i < numberOfParticles; ++i)
325  {
326  a = 0;
327  for (j = 0; j < numberOfParticles; ++j)
328  {
329  r2 = 0.0;
330  for (k = 0; k < DIM; ++k)
331  {
332  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
333  r2 += dx[k]*dx[k];
334  }
335 
336  if (r2 < cutoff2)
337  {
338  if ((half && i < j) || (!half && i != j))
339  {
340  /* part j is a neighbor of part i */
341  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
342  a++;
343  }
344  }
345  }
346  /* part i has `a' neighbors */
347  (*nl).NNeighbors[i] = a;
348  }
349 
350  return;
351 }
352 
353 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
354  int* numnei, int** nei1part, double** Rij)
355 {
356  /* local variables */
357  KIM_API_model * pkim = *(KIM_API_model **) kimmdl;
358  int partToReturn;
359  int status;
360  int* numberOfParticles;
361  NeighList* nl;
362 
363  /* initialize numnei */
364  *numnei = 0;
365 
366  /* unpack neighbor list object */
367  numberOfParticles = (int*) pkim->get_data("numberOfParticles", &status);
368  if (KIM_STATUS_OK > status)
369  {
370  pkim->report_error(__LINE__, __FILE__,"get_data", status);
371  return status;
372  }
373 
374  nl = (NeighList*) pkim->get_data("neighObject", &status);
375  if (KIM_STATUS_OK > status)
376  {
377  pkim->report_error(__LINE__, __FILE__,"get_data", status);
378  return status;
379  }
380 
381  /* check mode and request */
382  if (1 == *mode) /* locator mode */
383  {
384  if ((*request >= *numberOfParticles) || (*request < 0)) /* invalid id */
385  {
386  pkim->report_error(__LINE__, __FILE__,"Invalid part ID in get_cluster_neigh", KIM_STATUS_PARTICLE_INVALID_ID);
387  return KIM_STATUS_PARTICLE_INVALID_ID;
388  }
389  else
390  {
391  partToReturn = *request;
392  }
393  }
394  else /* invalid mode */
395  {
396  pkim->report_error(__LINE__, __FILE__,"Invalid mode in get_cluster_neigh", KIM_STATUS_NEIGH_INVALID_MODE);
397  return KIM_STATUS_NEIGH_INVALID_MODE;
398  }
399 
400  /* set the returned part */
401  *part = partToReturn;
402 
403  /* set the returned number of neighbors for the returned part */
404  *numnei = (*nl).NNeighbors[*part];
405 
406  /* set the location for the returned neighbor list */
407  *nei1part = &((*nl).neighborList[(*part)*NCLUSTERPARTS]);
408 
409  /* set the pointer to Rij to appropriate value */
410  *Rij = 0;
411 
412  return KIM_STATUS_OK;
413 }
414 
415 char const * const descriptor()
416 {
417  return
418  "KIM_API_Version := 1.6.0\n"
419  "Unit_length := A\n"
420  "Unit_energy := eV\n"
421  "Unit_charge := e\n"
422  "Unit_temperature := K\n"
423  "Unit_time := ps\n"
424  "\n"
425  "PARTICLE_SPECIES:\n"
426  "Ar spec 1\n"
427  "\n"
428  "CONVENTIONS:\n"
429  "ZeroBasedLists flag\n"
430  "Neigh_LocaAccess flag\n"
431  "NEIGH_PURE_F flag\n"
432  "\n"
433  "MODEL_INPUT:\n"
434  "numberOfParticles integer none []\n"
435  "numberOfSpecies integer none []\n"
436  "particleSpecies integer none [numberOfParticles]\n"
437  "coordinates double length [numberOfParticles,3]\n"
438  "get_neigh method none []\n"
439  "neighObject pointer none []\n"
440  "\n"
441  "MODEL_OUTPUT:\n"
442  "destroy method none []\n"
443  "compute method none []\n"
444  "reinit method none []\n"
445  "cutoff double length []\n"
446  "energy double energy []\n"
447  "forces double force [numberOfParticles,3]\n";
448 }
#define FCCSPACING
#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)
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