KIM API V2
ex_test_Ar_fcc_cluster.c
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 <stdio.h>
35 #include <math.h>
36 #include "KIM_API_C.h"
37 #include "KIM_API_status.h"
38 
39 #define NAMESTRLEN 128
40 
41 #define FCCSPACING 5.260
42 #define DIM 3
43 #define NCELLSPERSIDE 2
44 #define NCLUSTERPARTS (4*(NCELLSPERSIDE*NCELLSPERSIDE*NCELLSPERSIDE) + \
45  6*(NCELLSPERSIDE*NCELLSPERSIDE) \
46  + 3*(NCELLSPERSIDE) + 1)
47 
48 #define REPORT_ERROR(LN, FL, MSG, STAT) { \
49  KIM_API_report_error(LN, FL, MSG, STAT); \
50  exit(STAT); \
51  }
52 
53 /* Define neighborlist structure */
54 typedef struct
55 {
56  int iteratorId;
57  int* NNeighbors;
58  int* neighborList;
59 } NeighList;
60 
61 /* Define prototypes */
62 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
63  double cutoff, NeighList* nl);
64 
65 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
66  int* numnei, int** nei1part, double** Rij);
67 
68 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
69 
70 
71 /* Main program */
72 int main()
73 {
74  /* Local variable declarations */
75  double const MinSpacing = 0.8*FCCSPACING;
76  double const MaxSpacing = 1.2*FCCSPACING;
77  double const SpacingIncr = 0.025*FCCSPACING;
78  double CurrentSpacing;
79  double cutpad = 0.75; /* Angstroms */
80  double force_norm;
81  int i;
82  int status;
83 
84 
85  /* KIM variable declarations */
86  void* pkim_cluster_model;
87  /* model inputs */
88  int numberOfParticles_cluster = NCLUSTERPARTS;
89  int numberOfSpecies = 1;
90  int particleSpeciesShapeCluster[1];
91  int particleSpecies_cluster_model[NCLUSTERPARTS];
92  int coordinatesShapeCluster[2];
93  double coords_cluster[NCLUSTERPARTS][DIM];
94  NeighList nl_cluster_model;
95  /* model outputs */
96  double cutoff_cluster_model;
97  double energy_cluster_model;
98  double forces_cluster[NCLUSTERPARTS*DIM];
99 
100  char testkimfile[] = "descriptor.kim";
101  char modelname[NAMESTRLEN];
102 
103  /* Get KIM Model names */
104  printf("Please enter valid KIM Model name: \n");
105  status = scanf("%s", modelname);
106  if (1 != status)
107  {
108  REPORT_ERROR(__LINE__, __FILE__, "Unable to read model name",
109  status);
110  }
111 
112  /* initialize the model */
113  status = KIM_API_file_init(&pkim_cluster_model, testkimfile, modelname);
114  if (KIM_STATUS_OK > status)
115  REPORT_ERROR(__LINE__, __FILE__,"KIM_API_file_init()",status);
116 
117  /* Assign shapes & register memory */
118  coordinatesShapeCluster[0] = numberOfParticles_cluster;
119  coordinatesShapeCluster[1] = 3;
120 
121  KIM_API_set_shape(pkim_cluster_model, "particleSpecies", particleSpeciesShapeCluster, 1, &status);
122  if (KIM_STATUS_OK > status)
123  REPORT_ERROR(__LINE__, __FILE__,"KIM_API_set_shape",status);
124  KIM_API_set_shape(pkim_cluster_model, "coordinates", coordinatesShapeCluster, 2, &status);
125  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_set_shape",status);
126  KIM_API_setm_data(pkim_cluster_model, &status, 8*4,
127  "numberOfParticles", 1, &numberOfParticles_cluster, 1,
128  "numberOfSpecies", 1, &numberOfSpecies, 1,
129  "particleSpecies", numberOfParticles_cluster, &particleSpecies_cluster_model, 1,
130  "coordinates", DIM*numberOfParticles_cluster, coords_cluster, 1,
131  "neighObject", 1, &nl_cluster_model, 1,
132  "cutoff", 1, &cutoff_cluster_model, 1,
133  "energy", 1, &energy_cluster_model, 1,
134  "forces", DIM*numberOfParticles_cluster, forces_cluster, 1);
135  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_setm_data",status);
136  status = KIM_API_set_method(pkim_cluster_model, "get_neigh", 1, (func_ptr) &get_cluster_neigh);
137  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_set_method",status);
138 
139  /* call model init routine */
140  status = KIM_API_model_init(pkim_cluster_model);
141  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"KIM_API_model_init", status);
142 
143  /* setup particleSpecies */
144  particleSpecies_cluster_model[0] = KIM_API_get_species_code(pkim_cluster_model, "Ar", &status);
145  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"get_species_code", status);
146  for (i = 1; i < NCLUSTERPARTS; ++i)
147  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
148 
149  /* setup neighbor lists */
150  /* allocate memory for list */
151  nl_cluster_model.NNeighbors = (int*) malloc(NCLUSTERPARTS*sizeof(int));
152  if (NULL==nl_cluster_model.NNeighbors) REPORT_ERROR(__LINE__, __FILE__,"malloc unsuccessful", -1);
153 
154  nl_cluster_model.neighborList = (int*) malloc(NCLUSTERPARTS*NCLUSTERPARTS*sizeof(int));
155  if (NULL==nl_cluster_model.neighborList) REPORT_ERROR(__LINE__, __FILE__,"malloc unsuccessful", -1);
156 
157  /* ready to compute */
158  printf("This is Test : ex_test_Ar_fcc_cluster\n");
159  printf("--------------------------------------------------------------------------------\n");
160  printf("Results for KIM Model : %s\n", modelname);
161 
162  printf("%20s, %20s, %20s\n", "Energy", "Force Norm", "Lattice Spacing");
163  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
164  {
165  /* update coordinates for cluster */
166  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
167  /* compute neighbor lists */
168  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
169  (cutoff_cluster_model + cutpad), &nl_cluster_model);
170 
171  /* call compute functions */
172  status = KIM_API_model_compute(pkim_cluster_model);
173  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"compute", status);
174 
175  /* compute force norm */
176  force_norm = 0.0;
177  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
178  {
179  force_norm += forces_cluster[i]*forces_cluster[i];
180  }
181  force_norm = sqrt(force_norm);
182 
183  /* print the results */
184  printf("%20.10e, %20.10e, %20.10e\n",
185  energy_cluster_model,
186  force_norm,
187  CurrentSpacing);
188  }
189 
190  /* free memory of neighbor lists */
191  free(nl_cluster_model.NNeighbors);
192  free(nl_cluster_model.neighborList);
193 
194  /* call model destroy */
195  status = KIM_API_model_destroy(pkim_cluster_model);
196  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"destroy", status);
197 
198  /* free pkim objects */
199  KIM_API_free(&pkim_cluster_model, &status);
200  if (KIM_STATUS_OK > status) REPORT_ERROR(__LINE__, __FILE__,"free", status);
201 
202  /* everything is great */
203  return 0;
204 }
205 
206 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
207 {
208  /* local variables */
209  double FCCshifts[4][DIM];
210  double latVec[DIM];
211  int a;
212  int i;
213  int j;
214  int k;
215  int m;
216  int n;
217 
218  /* create a cubic FCC cluster of parts */
219  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
220  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
221  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
222  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
223 
224  a = 0;
225  for (i = 0; i < nCellsPerSide; ++i)
226  {
227  latVec[0] = ((double) i)*FCCspacing;
228  for (j = 0; j < nCellsPerSide; ++j)
229  {
230  latVec[1] = ((double) j)*FCCspacing;
231  for (k = 0; k < nCellsPerSide; ++k)
232  {
233  latVec[2] = ((double) k)*FCCspacing;
234  for (m = 0; m < 4; ++m)
235  {
236  for (n = 0; n < DIM; ++n)
237  {
238  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
239  }
240  a++;
241  }
242  }
243  /* add in the remaining three faces */
244  /* pos-x face */
245  latVec[0] = NCELLSPERSIDE*FCCspacing;
246  latVec[1] = ((double) i)*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[3][n];
256  }
257  a++;
258  /* pos-y face */
259  latVec[0] = ((double) i)*FCCspacing;
260  latVec[1] = NCELLSPERSIDE*FCCspacing;
261  latVec[2] = ((double) j)*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[2][n];
270  }
271  a++;
272  /* pos-z face */
273  latVec[0] = ((double) i)*FCCspacing;
274  latVec[1] = ((double) j)*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  for (n = 0; n < DIM; ++n)
282  {
283  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
284  }
285  a++;
286  }
287  /* add in the remaining three edges */
288  latVec[0] = ((double) i)*FCCspacing;
289  latVec[1] = NCELLSPERSIDE*FCCspacing;
290  latVec[2] = NCELLSPERSIDE*FCCspacing;
291  for (n = 0; n < DIM; ++n)
292  {
293  coords[a*DIM + n] = latVec[n];
294  }
295  a++;
296  latVec[0] = NCELLSPERSIDE*FCCspacing;
297  latVec[1] = ((double) i)*FCCspacing;
298  latVec[2] = NCELLSPERSIDE*FCCspacing;
299  for (n = 0; n < DIM; ++n)
300  {
301  coords[a*DIM + n] = latVec[n];
302  }
303  a++;
304  latVec[0] = NCELLSPERSIDE*FCCspacing;
305  latVec[1] = NCELLSPERSIDE*FCCspacing;
306  latVec[2] = ((double) i)*FCCspacing;
307  for (n = 0; n < DIM; ++n)
308  {
309  coords[a*DIM + n] = latVec[n];
310  }
311  a++;
312  }
313  /* add in the remaining corner */
314  for (n = 0; n < DIM; ++n)
315  {
316  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
317  }
318  a++;
319 
320  return;
321 }
322 
323 
324 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
325  double cutoff, NeighList* nl)
326 {
327  /* local variables */
328  int i;
329  int j;
330  int k;
331  int a;
332 
333  double dx[DIM];
334  double r2;
335  double cutoff2;
336 
337  cutoff2 = cutoff*cutoff;
338 
339  for (i = 0; i < numberOfParticles; ++i)
340  {
341  a = 0;
342  for (j = 0; j < numberOfParticles; ++j)
343  {
344  r2 = 0.0;
345  for (k = 0; k < DIM; ++k)
346  {
347  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
348  r2 += dx[k]*dx[k];
349  }
350 
351  if (r2 < cutoff2)
352  {
353  if ((half && i < j) || (!half && i != j))
354  {
355  /* part j is a neighbor of part i */
356  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
357  a++;
358  }
359  }
360  }
361  /* part i has `a' neighbors */
362  (*nl).NNeighbors[i] = a;
363  }
364 
365  return;
366 }
367 
368 int get_cluster_neigh(void* kimmdl, int *mode, int *request, int* part,
369  int* numnei, int** nei1part, double** Rij)
370 {
371  /* local variables */
372  intptr_t* pkim = *((intptr_t**) kimmdl);
373  int partToReturn;
374  int status;
375  int* numberOfParticles;
376  NeighList* nl;
377 
378  /* initialize numnei */
379  *numnei = 0;
380 
381  /* unpack neighbor list object */
382  numberOfParticles = (int*) KIM_API_get_data(pkim, "numberOfParticles", &status);
383  if (KIM_STATUS_OK > status)
384  {
385  KIM_API_report_error(__LINE__, __FILE__,"get_data", status);
386  return status;
387  }
388 
389  nl = (NeighList*) KIM_API_get_data(pkim, "neighObject", &status);
390  if (KIM_STATUS_OK > status)
391  {
392  KIM_API_report_error(__LINE__, __FILE__,"get_data", status);
393  return status;
394  }
395 
396  /* check mode and request */
397  if (0 == *mode) /* iterator mode */
398  {
399  if (0 == *request) /* reset iterator */
400  {
401  (*nl).iteratorId = -1;
402  return KIM_STATUS_NEIGH_ITER_INIT_OK;
403  }
404  else if (1 == *request) /* increment iterator */
405  {
406  (*nl).iteratorId++;
407  if ((*nl).iteratorId >= *numberOfParticles)
408  {
409  return KIM_STATUS_NEIGH_ITER_PAST_END;
410  }
411  else
412  {
413  partToReturn = (*nl).iteratorId;
414  }
415  }
416  else /* invalid request value */
417  {
418  KIM_API_report_error(__LINE__, __FILE__,"Invalid request in get_cluster_neigh", KIM_STATUS_NEIGH_INVALID_REQUEST);
419  return KIM_STATUS_NEIGH_INVALID_REQUEST;
420  }
421  }
422  else if (1 == *mode) /* locator mode */
423  {
424  if ((*request >= *numberOfParticles) || (*request < 0)) /* invalid id */
425  {
426  KIM_API_report_error(__LINE__, __FILE__,"Invalid part ID in get_cluster_neigh", KIM_STATUS_PARTICLE_INVALID_ID);
427  return KIM_STATUS_PARTICLE_INVALID_ID;
428  }
429  else
430  {
431  partToReturn = *request;
432  }
433  }
434  else /* invalid mode */
435  {
436  KIM_API_report_error(__LINE__, __FILE__,"Invalid mode in get_cluster_neigh", KIM_STATUS_NEIGH_INVALID_MODE);
437  return KIM_STATUS_NEIGH_INVALID_MODE;
438  }
439 
440  /* set the returned part */
441  *part = partToReturn;
442 
443  /* set the returned number of neighbors for the returned part */
444  *numnei = (*nl).NNeighbors[*part];
445 
446  /* set the location for the returned neighbor list */
447  *nei1part = &((*nl).neighborList[(*part)*NCLUSTERPARTS]);
448 
449  /* set the pointer to Rij to appropriate value */
450  *Rij = 0;
451 
452  return KIM_STATUS_OK;
453 }
#define REPORT_ERROR(LN, FL, MSG, STAT)
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
#define FCCSPACING
#define DIM
#define NAMESTRLEN
LengthUnit const m
#define NCELLSPERSIDE
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
#define NCLUSTERPARTS
ComputeArgumentName const numberOfParticles
int main()
void fcc_cluster_neighborlist(int half, int numberOfParticles, double *coords, double cutoff, NeighList *nl)