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  KIM_Model * model;
103  /* model inputs */
104  int numberOfParticles_cluster = NCLUSTERPARTS;
105  int speciesIsSupported;
106  int particleSpecies_cluster_model[NCLUSTERPARTS];
107  int coordinatesShapeCluster[2];
108  double coords_cluster[NCLUSTERPARTS][DIM];
109  NeighList nl_cluster_model;
110  /* model outputs */
111  int number_of_cutoffs_cluster_model;
112  double influence_distance_cluster_model;
113  double const * cutoff_cluster_model;
114  double energy_cluster_model;
115  double forces_cluster[NCLUSTERPARTS*DIM];
116 
117  char modelname[NAMESTRLEN];
118  int requestedUnitsAccepted;
119  int modelArCode;
120  int numberOfComputeArgumentNames;
121  KIM_ComputeArguments * computeArguments;
122  KIM_ComputeArgumentName computeArgumentName;
123  KIM_SupportStatus supportStatus;
124  int numberOfComputeCallbackNames;
125  KIM_ComputeCallbackName computeCallbackName;
126 
127  /* Get KIM Model names */
128  printf("Please enter valid KIM Model name: \n");
129  error = scanf("%s", modelname);
130  if (1 != error)
131  {
132  MY_ERROR("Unable to read model name");
133  }
134 
135  /* initialize the model */
142  modelname,
143  &requestedUnitsAccepted,
144  &model);
145  if (error) MY_ERROR("KIM_create_model_interface()");
146 
147  /* Check for compatibility with the model */
148  if (!requestedUnitsAccepted) MY_ERROR("Must adapt to model units");
149 
150  /* check species */
152  model, KIM_SPECIES_NAME_Ar, &speciesIsSupported, &modelArCode);
153  if ((error) || (!speciesIsSupported))
154  {
155  MY_ERROR("Species Ar not supported");
156  }
157 
158  error = KIM_Model_ComputeArgumentsCreate(model, &computeArguments);
159  if (error)
160  {
161  MY_ERROR("KIM_Model_ComputeArgumentsCreate");
162  }
163 
164  /* check arguments */
166  &numberOfComputeArgumentNames);
167  for (i=0; i<numberOfComputeArgumentNames; ++i)
168  {
170  i, &computeArgumentName);
171  if (error) MY_ERROR("can't get argument name");
173  computeArgumentName,
174  &supportStatus);
175  if (error) MY_ERROR("can't get argument supportStatus");
176 
177  /* can only handle energy and force as a required arg */
179  {
181  computeArgumentName,
184  computeArgumentName,
186  {
187  MY_ERROR("unsupported required argument");
188  }
189  }
190 
191  /* must have energy and forces */
193  computeArgumentName,
195  ||
197  computeArgumentName,
199  {
200  if (! ((KIM_SupportStatus_Equal(supportStatus,
202  ||
203  (KIM_SupportStatus_Equal(supportStatus,
205  {
206  MY_ERROR("energy or forces not available");
207  }
208  }
209  }
210 
211  /* check call backs */
213  &numberOfComputeCallbackNames);
214  for (i=0; i<numberOfComputeCallbackNames; ++i)
215  {
217  i, &computeCallbackName);
218  if (error) MY_ERROR("can't get call back name");
220  computeCallbackName,
221  &supportStatus);
222  if (error) MY_ERROR("can't get call back supportStatus");
223 
224  /* cannot handle any "required" call backs */
226  {
227  MY_ERROR("unsupported required call back");
228  }
229  }
230 
231  /* We're compatible with the model. Let's do it. */
232 
233  /* setup particleSpecies */
234  particleSpecies_cluster_model[0] = KIM_API_get_species_code(pkim_cluster_model, "Ar", &error);
235  if (KIM_ERROR_OK > error) MY_ERROR("get_species_code");
236  for (i = 1; i < NCLUSTERPARTS; ++i)
237  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
238 
239  /* setup neighbor lists */
240  /* allocate memory for list */
241  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
242  nl_cluster_model.NNeighbors = (int*) malloc(NCLUSTERPARTS*sizeof(int));
243  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("malloc unsuccessful");
244 
245  nl_cluster_model.neighborList = (int*) malloc(NCLUSTERPARTS*NCLUSTERPARTS*sizeof(int));
246  if (NULL==nl_cluster_model.neighborList) MY_ERROR("malloc unsuccessful");
247 
248  /* ready to compute */
249  printf("This is Test : ex_test_Ar_fcc_cluster\n");
250  printf("--------------------------------------------------------------------------------\n");
251  printf("Results for KIM Model : %s\n", modelname);
252 
253  printf("%20s, %20s, %20s\n", "Energy", "Force Norm", "Lattice Spacing");
254  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
255  {
256  /* update coordinates for cluster */
257  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
258  /* compute neighbor lists */
259  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
260  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
261 
262  /* call compute functions */
263  error = KIM_API_model_compute(pkim_cluster_model);
264  if (KIM_ERROR_OK > error) MY_ERROR("compute");
265 
266  /* compute force norm */
267  force_norm = 0.0;
268  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
269  {
270  force_norm += forces_cluster[i]*forces_cluster[i];
271  }
272  force_norm = sqrt(force_norm);
273 
274  /* print the results */
275  printf("%20.10e, %20.10e, %20.10e\n",
276  energy_cluster_model,
277  force_norm,
278  CurrentSpacing);
279  }
280 
281  /* free memory of neighbor lists */
282  free(nl_cluster_model.NNeighbors);
283  free(nl_cluster_model.neighborList);
284 
285  /* call model destroy */
286  error = KIM_API_model_destroy(pkim_cluster_model);
287  if (KIM_ERROR_OK > error) MY_ERROR("destroy");
288 
289  /* free pkim objects */
290  KIM_API_free(&pkim_cluster_model, &error);
291  if (KIM_ERROR_OK > error) MY_ERROR("free");
292 
293  /* everything is great */
294  return 0;
295 }
296 
297 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
298 {
299  /* local variables */
300  double FCCshifts[4][DIM];
301  double latVec[DIM];
302  int a;
303  int i;
304  int j;
305  int k;
306  int m;
307  int n;
308 
309  /* create a cubic FCC cluster of parts */
310  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
311  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
312  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
313  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
314 
315  a = 0;
316  for (i = 0; i < nCellsPerSide; ++i)
317  {
318  latVec[0] = ((double) i)*FCCspacing;
319  for (j = 0; j < nCellsPerSide; ++j)
320  {
321  latVec[1] = ((double) j)*FCCspacing;
322  for (k = 0; k < nCellsPerSide; ++k)
323  {
324  latVec[2] = ((double) k)*FCCspacing;
325  for (m = 0; m < 4; ++m)
326  {
327  for (n = 0; n < DIM; ++n)
328  {
329  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
330  }
331  a++;
332  }
333  }
334  /* add in the remaining three faces */
335  /* pos-x face */
336  latVec[0] = NCELLSPERSIDE*FCCspacing;
337  latVec[1] = ((double) i)*FCCspacing;
338  latVec[2] = ((double) j)*FCCspacing;
339  for (n = 0; n < DIM; ++n)
340  {
341  coords[a*DIM + n] = latVec[n];
342  }
343  a++;
344  for (n = 0; n < DIM; ++n)
345  {
346  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
347  }
348  a++;
349  /* pos-y face */
350  latVec[0] = ((double) i)*FCCspacing;
351  latVec[1] = NCELLSPERSIDE*FCCspacing;
352  latVec[2] = ((double) j)*FCCspacing;
353  for (n = 0; n < DIM; ++n)
354  {
355  coords[a*DIM + n] = latVec[n];
356  }
357  a++;
358  for (n = 0; n < DIM; ++n)
359  {
360  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
361  }
362  a++;
363  /* pos-z face */
364  latVec[0] = ((double) i)*FCCspacing;
365  latVec[1] = ((double) j)*FCCspacing;
366  latVec[2] = NCELLSPERSIDE*FCCspacing;
367  for (n = 0; n < DIM; ++n)
368  {
369  coords[a*DIM + n] = latVec[n];
370  }
371  a++;
372  for (n = 0; n < DIM; ++n)
373  {
374  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
375  }
376  a++;
377  }
378  /* add in the remaining three edges */
379  latVec[0] = ((double) i)*FCCspacing;
380  latVec[1] = NCELLSPERSIDE*FCCspacing;
381  latVec[2] = NCELLSPERSIDE*FCCspacing;
382  for (n = 0; n < DIM; ++n)
383  {
384  coords[a*DIM + n] = latVec[n];
385  }
386  a++;
387  latVec[0] = NCELLSPERSIDE*FCCspacing;
388  latVec[1] = ((double) i)*FCCspacing;
389  latVec[2] = NCELLSPERSIDE*FCCspacing;
390  for (n = 0; n < DIM; ++n)
391  {
392  coords[a*DIM + n] = latVec[n];
393  }
394  a++;
395  latVec[0] = NCELLSPERSIDE*FCCspacing;
396  latVec[1] = NCELLSPERSIDE*FCCspacing;
397  latVec[2] = ((double) i)*FCCspacing;
398  for (n = 0; n < DIM; ++n)
399  {
400  coords[a*DIM + n] = latVec[n];
401  }
402  a++;
403  }
404  /* add in the remaining corner */
405  for (n = 0; n < DIM; ++n)
406  {
407  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
408  }
409  a++;
410 
411  return;
412 }
413 
414 
415 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
416  double cutoff, NeighList* nl)
417 {
418  /* local variables */
419  int i;
420  int j;
421  int k;
422  int a;
423 
424  double dx[DIM];
425  double r2;
426  double cutoff2;
427 
428  nl->cutoff = cutoff;
429 
430  cutoff2 = cutoff*cutoff;
431 
432  for (i = 0; i < numberOfParticles; ++i)
433  {
434  a = 0;
435  for (j = 0; j < numberOfParticles; ++j)
436  {
437  r2 = 0.0;
438  for (k = 0; k < DIM; ++k)
439  {
440  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
441  r2 += dx[k]*dx[k];
442  }
443 
444  if (r2 < cutoff2)
445  {
446  if ((half && i < j) || (!half && i != j))
447  {
448  /* part j is a neighbor of part i */
449  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
450  a++;
451  }
452  }
453  }
454  /* part i has `a' neighbors */
455  (*nl).NNeighbors[i] = a;
456  }
457 
458  return;
459 }
460 
461 int get_cluster_neigh(void const * const dataObject,
462  int const numberOfCutoffs, double const * const cutoffs,
463  int const neighborListIndex,
464  int const particleNumber, int * const numberOfNeighbors,
465  int const ** const neighborsOfParticle)
466 {
467  /* local variables */
468  int error = TRUE;
469  NeighList* nl = (NeighList*) dataObject;
470  int numberOfParticles = nl->numberOfParticles;
471 
472  if ((numberOfCutoffs != 1) || (cutoffs[0] > nl->cutoff)) return error;
473 
474  if (neighborListIndex != 0) return error;
475 
476  /* initialize numNeigh */
477  *numberOfNeighbors = 0;
478 
479  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
480  {
481  MY_WARNING("Invalid part ID in get_cluster_neigh");
482  return TRUE;
483  }
484 
485  /* set the returned number of neighbors for the returned part */
486  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
487 
488  /* set the location for the returned neighbor list */
489  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
490 
491  return FALSE;
492 }
struct KIM_ComputeArguments KIM_ComputeArguments
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
KIM_SupportStatus const KIM_SUPPORT_STATUS_optional
int KIM_ComputeArguments_GetCallbackSupportStatus(KIM_ComputeArguments const *const computeArguments, KIM_ComputeCallbackName const computeCallbackName, KIM_SupportStatus *const supportStatus)
void KIM_COMPUTE_ARGUMENT_NAME_GetNumberOfComputeArgumentNames(int *const numberOfComputeArgumentNames)
void KIM_COMPUTE_CALLBACK_NAME_GetNumberOfComputeCallbackNames(int *const numberOfComputeCallbackNames)
#define MY_WARNING(message)
#define NCLUSTERPARTS
KIM_ChargeUnit const KIM_CHARGE_UNIT_e
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialForces
#define FALSE
int KIM_Model_GetSpeciesSupportAndCode(KIM_Model const *const model, KIM_SpeciesName const speciesName, int *const speciesIsSupported, int *const code)
KIM_SpeciesName const KIM_SPECIES_NAME_Ar
int KIM_Model_ComputeArgumentsCreate(KIM_Model const *const model, KIM_ComputeArguments **const computeArguments)
#define MY_ERROR(message)
#define NCELLSPERSIDE
LengthUnit const m
int KIM_Model_Create(KIM_Numbering const numbering, KIM_LengthUnit const requestedLengthUnit, KIM_EnergyUnit const requestedEnergyUnit, KIM_ChargeUnit const requestedChargeUnit, KIM_TemperatureUnit const requestedTemperatureUnit, KIM_TimeUnit const requestedTimeUnit, char const *const modelName, int *const requestedUnitsAccepted, KIM_Model **const model)
#define FCCSPACING
#define TRUE
int KIM_ComputeArgumentName_NotEqual(KIM_ComputeArgumentName const left, KIM_ComputeArgumentName const right)
KIM_TimeUnit const KIM_TIME_UNIT_ps
int KIM_COMPUTE_CALLBACK_NAME_GetComputeCallbackName(int const index, KIM_ComputeCallbackName *const computeCallbackName)
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
int KIM_SupportStatus_Equal(KIM_SupportStatus const left, KIM_SupportStatus const right)
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
#define NAMESTRLEN
ComputeArgumentName const numberOfParticles
int KIM_COMPUTE_ARGUMENT_NAME_GetComputeArgumentName(int const index, KIM_ComputeArgumentName *const computeArgumentName)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialEnergy
#define DIM
KIM_SupportStatus const KIM_SUPPORT_STATUS_required
int main()
int KIM_ComputeArguments_GetArgumentSupportStatus(KIM_ComputeArguments const *const computeArguments, KIM_ComputeArgumentName const computeArgumentName, KIM_SupportStatus *const supportStatus)
struct KIM_Model KIM_Model
Definition: KIM_Model.h:104
void fcc_cluster_neighborlist(int half, int numberOfParticles, double *coords, double cutoff, NeighList *nl)
KIM_LengthUnit const KIM_LENGTH_UNIT_A
LogVerbosity const error
int KIM_ComputeArgumentName_Equal(KIM_ComputeArgumentName const left, KIM_ComputeArgumentName const right)
KIM_TemperatureUnit const KIM_TEMPERATURE_UNIT_K
KIM_Numbering const KIM_NUMBERING_zeroBased