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 particleContributing_cluster_model[NCLUSTERPARTS];
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  error =
235  computeArguments,
237  &numberOfParticles_cluster)
238  ||
240  computeArguments,
242  particleSpecies_cluster_model)
243  ||
245  computeArguments,
247  particleContributing_cluster_model)
248  ||
250  computeArguments, KIM_COMPUTE_ARGUMENT_NAME_coordinates,
251  (double*) &coords_cluster)
252  ||
254  computeArguments, KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
255  &energy_cluster_model)
256  ||
258  computeArguments, KIM_COMPUTE_ARGUMENT_NAME_partialForces,
259  (double*) &forces_cluster);
260  if (error) MY_ERROR("KIM_setm_data");
262  computeArguments,
265  (func *) &get_cluster_neigh,
266  &nl_cluster_model);
267 
268  KIM_Model_GetInfluenceDistance(model, &influence_distance_cluster_model);
270  &number_of_cutoffs_cluster_model,
271  &cutoff_cluster_model);
272  if (number_of_cutoffs_cluster_model != 1) MY_ERROR("too many cutoffs");
273 
274  /* setup particleSpecies */
276  model,
278  &speciesIsSupported,
279  &(particleSpecies_cluster_model[0]));
280  if (error) MY_ERROR("KIM_get_species_code");
281  for (i = 1; i < NCLUSTERPARTS; ++i)
282  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
283 
284  /* setup particleContributing */
285  for (i = 0; i < NCLUSTERPARTS; ++i)
286  {
287  particleContributing_cluster_model[i] = 1; /* all particles contribute */
288  }
289 
290  /* setup neighbor lists */
291  /* allocate memory for list */
292  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
293  nl_cluster_model.NNeighbors = (int*) malloc(NCLUSTERPARTS*sizeof(int));
294  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("malloc unsuccessful");
295 
296  nl_cluster_model.neighborList = (int*) malloc(NCLUSTERPARTS*NCLUSTERPARTS*sizeof(int));
297  if (NULL==nl_cluster_model.neighborList) MY_ERROR("malloc unsuccessful");
298 
299  /* ready to compute */
300  printf("This is Test : ex_test_Ar_fcc_cluster\n");
301  printf("--------------------------------------------------------------------------------\n");
302  printf("Results for KIM Model : %s\n", modelname);
303 
304  printf("%20s, %20s, %20s\n", "Energy", "Force Norm", "Lattice Spacing");
305  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
306  {
307  /* update coordinates for cluster */
308  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
309  /* compute neighbor lists */
310  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
311  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
312 
313  /* call compute functions */
314  error = KIM_Model_Compute(model, computeArguments);
315  if (error) MY_ERROR("KIM_model_compute");
316 
317  /* compute force norm */
318  force_norm = 0.0;
319  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
320  {
321  force_norm += forces_cluster[i]*forces_cluster[i];
322  }
323  force_norm = sqrt(force_norm);
324 
325  /* print the results */
326  printf("%20.10e, %20.10e, %20.10e\n",
327  energy_cluster_model,
328  force_norm,
329  CurrentSpacing);
330  }
331 
332  error = KIM_Model_ComputeArgumentsDestroy(model, &computeArguments);
333  if (error)
334  {
335  MY_ERROR("Unable to destroy compute arguments");
336  }
337 
338  /* free memory of neighbor lists */
339  free(nl_cluster_model.NNeighbors);
340  free(nl_cluster_model.neighborList);
341 
342  /* free pkim objects */
343  KIM_Model_Destroy(&model);
344 
345  /* everything is great */
346  return 0;
347 }
348 
349 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
350 {
351  /* local variables */
352  double FCCshifts[4][DIM];
353  double latVec[DIM];
354  int a;
355  int i;
356  int j;
357  int k;
358  int m;
359  int n;
360 
361  /* create a cubic FCC cluster of parts */
362  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
363  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
364  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
365  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
366 
367  a = 0;
368  for (i = 0; i < nCellsPerSide; ++i)
369  {
370  latVec[0] = ((double) i)*FCCspacing;
371  for (j = 0; j < nCellsPerSide; ++j)
372  {
373  latVec[1] = ((double) j)*FCCspacing;
374  for (k = 0; k < nCellsPerSide; ++k)
375  {
376  latVec[2] = ((double) k)*FCCspacing;
377  for (m = 0; m < 4; ++m)
378  {
379  for (n = 0; n < DIM; ++n)
380  {
381  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
382  }
383  a++;
384  }
385  }
386  /* add in the remaining three faces */
387  /* pos-x face */
388  latVec[0] = NCELLSPERSIDE*FCCspacing;
389  latVec[1] = ((double) i)*FCCspacing;
390  latVec[2] = ((double) j)*FCCspacing;
391  for (n = 0; n < DIM; ++n)
392  {
393  coords[a*DIM + n] = latVec[n];
394  }
395  a++;
396  for (n = 0; n < DIM; ++n)
397  {
398  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
399  }
400  a++;
401  /* pos-y face */
402  latVec[0] = ((double) i)*FCCspacing;
403  latVec[1] = NCELLSPERSIDE*FCCspacing;
404  latVec[2] = ((double) j)*FCCspacing;
405  for (n = 0; n < DIM; ++n)
406  {
407  coords[a*DIM + n] = latVec[n];
408  }
409  a++;
410  for (n = 0; n < DIM; ++n)
411  {
412  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
413  }
414  a++;
415  /* pos-z face */
416  latVec[0] = ((double) i)*FCCspacing;
417  latVec[1] = ((double) j)*FCCspacing;
418  latVec[2] = NCELLSPERSIDE*FCCspacing;
419  for (n = 0; n < DIM; ++n)
420  {
421  coords[a*DIM + n] = latVec[n];
422  }
423  a++;
424  for (n = 0; n < DIM; ++n)
425  {
426  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
427  }
428  a++;
429  }
430  /* add in the remaining three edges */
431  latVec[0] = ((double) i)*FCCspacing;
432  latVec[1] = NCELLSPERSIDE*FCCspacing;
433  latVec[2] = NCELLSPERSIDE*FCCspacing;
434  for (n = 0; n < DIM; ++n)
435  {
436  coords[a*DIM + n] = latVec[n];
437  }
438  a++;
439  latVec[0] = NCELLSPERSIDE*FCCspacing;
440  latVec[1] = ((double) i)*FCCspacing;
441  latVec[2] = NCELLSPERSIDE*FCCspacing;
442  for (n = 0; n < DIM; ++n)
443  {
444  coords[a*DIM + n] = latVec[n];
445  }
446  a++;
447  latVec[0] = NCELLSPERSIDE*FCCspacing;
448  latVec[1] = NCELLSPERSIDE*FCCspacing;
449  latVec[2] = ((double) i)*FCCspacing;
450  for (n = 0; n < DIM; ++n)
451  {
452  coords[a*DIM + n] = latVec[n];
453  }
454  a++;
455  }
456  /* add in the remaining corner */
457  for (n = 0; n < DIM; ++n)
458  {
459  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
460  }
461  a++;
462 
463  return;
464 }
465 
466 
467 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
468  double cutoff, NeighList* nl)
469 {
470  /* local variables */
471  int i;
472  int j;
473  int k;
474  int a;
475 
476  double dx[DIM];
477  double r2;
478  double cutoff2;
479 
480  nl->cutoff = cutoff;
481 
482  cutoff2 = cutoff*cutoff;
483 
484  for (i = 0; i < numberOfParticles; ++i)
485  {
486  a = 0;
487  for (j = 0; j < numberOfParticles; ++j)
488  {
489  r2 = 0.0;
490  for (k = 0; k < DIM; ++k)
491  {
492  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
493  r2 += dx[k]*dx[k];
494  }
495 
496  if (r2 < cutoff2)
497  {
498  if ((half && i < j) || (!half && i != j))
499  {
500  /* part j is a neighbor of part i */
501  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
502  a++;
503  }
504  }
505  }
506  /* part i has `a' neighbors */
507  (*nl).NNeighbors[i] = a;
508  }
509 
510  return;
511 }
512 
513 int get_cluster_neigh(void const * const dataObject,
514  int const numberOfCutoffs, double const * const cutoffs,
515  int const neighborListIndex,
516  int const particleNumber, int * const numberOfNeighbors,
517  int const ** const neighborsOfParticle)
518 {
519  /* local variables */
520  int error = TRUE;
521  NeighList* nl = (NeighList*) dataObject;
522  int numberOfParticles = nl->numberOfParticles;
523 
524  if ((numberOfCutoffs != 1) || (cutoffs[0] > nl->cutoff)) return error;
525 
526  if (neighborListIndex != 0) return error;
527 
528  /* initialize numNeigh */
529  *numberOfNeighbors = 0;
530 
531  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
532  {
533  MY_WARNING("Invalid part ID in get_cluster_neigh");
534  return TRUE;
535  }
536 
537  /* set the returned number of neighbors for the returned part */
538  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
539 
540  /* set the location for the returned neighbor list */
541  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
542 
543  return FALSE;
544 }
struct KIM_ComputeArguments KIM_ComputeArguments
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
#define FCCSPACING
void KIM_Model_GetInfluenceDistance(KIM_Model const *const model, double *const influenceDistance)
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)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_coordinates
void KIM_COMPUTE_CALLBACK_NAME_GetNumberOfComputeCallbackNames(int *const numberOfComputeCallbackNames)
#define TRUE
void KIM_Model_GetNeighborListCutoffsPointer(KIM_Model const *const model, int *const numberOfCutoffs, double const **const cutoffs)
void KIM_Model_Destroy(KIM_Model **const model)
KIM_ChargeUnit const KIM_CHARGE_UNIT_e
KIM_LanguageName const KIM_LANGUAGE_NAME_c
int KIM_ComputeArguments_SetCallbackPointer(KIM_ComputeArguments *const computeArguments, KIM_ComputeCallbackName const computeCallbackName, KIM_LanguageName const languageName, func *const fptr, void const *const dataObject)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialForces
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes
int KIM_Model_GetSpeciesSupportAndCode(KIM_Model const *const model, KIM_SpeciesName const speciesName, int *const speciesIsSupported, int *const code)
#define NAMESTRLEN
#define DIM
int KIM_ComputeArguments_SetArgumentPointerInteger(KIM_ComputeArguments *const computeArguments, KIM_ComputeArgumentName const computeArgumentName, int const *const ptr)
KIM_SpeciesName const KIM_SPECIES_NAME_Ar
int KIM_Model_ComputeArgumentsCreate(KIM_Model const *const model, KIM_ComputeArguments **const computeArguments)
void() func()
Definition: KIM_func.h:39
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)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles
int KIM_ComputeArgumentName_NotEqual(KIM_ComputeArgumentName const left, KIM_ComputeArgumentName const right)
KIM_ComputeCallbackName const KIM_COMPUTE_CALLBACK_NAME_GetNeighborList
KIM_TimeUnit const KIM_TIME_UNIT_ps
#define MY_WARNING(message)
#define NCLUSTERPARTS
int KIM_ComputeArguments_SetArgumentPointerDouble(KIM_ComputeArguments *const computeArguments, KIM_ComputeArgumentName const computeArgumentName, double const *const ptr)
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_Model_ComputeArgumentsDestroy(KIM_Model const *const model, KIM_ComputeArguments **const computeArguments)
int KIM_SupportStatus_Equal(KIM_SupportStatus const left, KIM_SupportStatus const right)
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
#define FALSE
ComputeArgumentName const numberOfParticles
int KIM_Model_Compute(KIM_Model const *const model, KIM_ComputeArguments const *const computeArguments)
int KIM_COMPUTE_ARGUMENT_NAME_GetComputeArgumentName(int const index, KIM_ComputeArgumentName *const computeArgumentName)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialEnergy
#define MY_ERROR(message)
KIM_SupportStatus const KIM_SUPPORT_STATUS_required
#define NCELLSPERSIDE
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
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_particleContributing