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