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