KIM API V2
ex_test_Ar_fcc_cluster_cpp.cpp
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 // Release: This file is part of the kim-api-v2.0.0-alpha.0 package.
34 //
35 
36 #include <stdlib.h>
37 #include <iostream>
38 #include <iomanip>
39 #include "KIM_LogVerbosity.hpp"
40 #include "KIM_LanguageName.hpp"
41 #include "KIM_DataType.hpp"
42 #include "KIM_SpeciesName.hpp"
43 #include "KIM_Numbering.hpp"
44 #include "KIM_Model.hpp"
45 #include "KIM_ArgumentName.hpp"
46 #include "KIM_CallbackName.hpp"
47 #include "KIM_SupportStatus.hpp"
48 #include "KIM_UnitSystem.hpp"
49 
50 #define NAMESTRLEN 128
51 
52 #define FCCSPACING 5.260
53 #define DIM 3
54 #define NCELLSPERSIDE 2
55 #define NCLUSTERPARTS (4*(NCELLSPERSIDE*NCELLSPERSIDE*NCELLSPERSIDE) + \
56  6*(NCELLSPERSIDE*NCELLSPERSIDE) \
57  + 3*(NCELLSPERSIDE) + 1)
58 
59 #define MY_ERROR(message) \
60  { \
61  std::cout << "* Error : \"" << message << "\" : " \
62  << __LINE__ << ":" << __FILE__ << std::endl; \
63  exit(1); \
64  }
65 
66 #define MY_WARNING(message) \
67  { \
68  std::cout << "* Error : \"" << message << "\" : " \
69  << __LINE__ << ":" << __FILE__ << std::endl; \
70  }
71 
72 
73 /* Define neighborlist structure */
74 typedef struct
75 {
77  int iteratorId;
78  int* NNeighbors;
79  int* neighborList;
80 } NeighList;
81 
82 /* Define prototypes */
83 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
84  double cutoff, NeighList* nl);
85 
86 int get_cluster_neigh(void const * const dataObject,
87  int const neighborListIndex, int const particleNumber,
88  int * const numberOfNeighbors,
89  int const ** const neighborsOfParticle);
90 
91 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
92 
93 
94 /* Main program */
95 int main()
96 {
97  /* Local variable declarations */
98  double const MinSpacing = 0.8*FCCSPACING;
99  double const MaxSpacing = 1.2*FCCSPACING;
100  double const SpacingIncr = 0.025*FCCSPACING;
101  double CurrentSpacing;
102  double cutpad = 0.75; /* Angstroms */
103  int i;
104  int error;
105 
106 
107  /* model inputs */
108  int numberOfParticles_cluster = NCLUSTERPARTS;
109  int particleSpecies_cluster_model[NCLUSTERPARTS];
110  int particleContributing_cluster_model[NCLUSTERPARTS];
111  double coords_cluster[NCLUSTERPARTS][DIM];
112  NeighList nl_cluster_model;
113  /* model outputs */
114  double influence_distance_cluster_model;
115  int number_of_cutoffs;
116  double const * cutoff_cluster_model;
117  double energy_cluster_model;
118 
119  std::string modelname;
120 
121  /* Get KIM Model names */
122  std::cout << "Please enter valid KIM Model name: \n";
123  std::cin >> modelname;
124 
125 
126  /* initialize the model */
127  KIM::Model * kim_cluster_model;
128  int requestedUnitsAccepted;
136  modelname,
137  &requestedUnitsAccepted,
138  &kim_cluster_model);
139  if (error)
140  {
141  MY_ERROR("KIM_create_model_interface()");
142  }
143 
144  // Check for compatibility with the model
145  if (!requestedUnitsAccepted)
146  {
147  MY_ERROR("Must Adapt to model units");
148  }
149 
150  // print model units
151  KIM::LengthUnit lengthUnit;
152  KIM::EnergyUnit energyUnit;
153  KIM::ChargeUnit chargeUnit;
154  KIM::TemperatureUnit temperatureUnit;
155  KIM::TimeUnit timeUnit;
156 
157  kim_cluster_model->GetUnits(&lengthUnit, &energyUnit, &chargeUnit,
158  &temperatureUnit, &timeUnit);
159 
160  std::cout << "LengthUnit is \"" << lengthUnit.String() << "\"" << std::endl
161  << "EnergyUnit is \"" << energyUnit.String() << "\"" << std::endl
162  << "ChargeUnit is \"" << chargeUnit.String() << "\"" << std::endl
163  << "TemperatureUnit is \"" << temperatureUnit.String()
164  << "\"" << std::endl
165  << "TimeUnit is \"" << timeUnit.String() << "\"" << std::endl;
166 
167  // check species
168  int speciesIsSupported;
169  int modelArCode;
170  error = kim_cluster_model->GetSpeciesSupportAndCode(
171  KIM::SPECIES_NAME::Ar, &speciesIsSupported, &modelArCode);
172  if ((error) || (!speciesIsSupported))
173  {
174  MY_ERROR("Species Ar not supported");
175  }
176 
177  // check arguments
178  int numberOfArguments;
179  KIM::ARGUMENT_NAME::GetNumberOfArguments(&numberOfArguments);
180  for (int i=0; i<numberOfArguments; ++i)
181  {
182  KIM::ArgumentName argumentName;
183  KIM::SupportStatus supportStatus;
184  KIM::ARGUMENT_NAME::GetArgumentName(i, &argumentName);
185  KIM::DataType dataType;
186  KIM::ARGUMENT_NAME::GetArgumentDataType(argumentName, &dataType);
187  error = kim_cluster_model->GetArgumentSupportStatus(argumentName,
188  &supportStatus);
189  if (error)
190  MY_ERROR("unable to get argument supportStatus");
191 
192  std::cout << "Argument Name \""
193  << argumentName.String() << "\""
194  << " is of type \""
195  << dataType.String() << "\""
196  << " and has supportStatus \""
197  << supportStatus.String() << "\""
198  << std::endl;
199 
200  // can only handle energy as a required arg
201  if (supportStatus == KIM::SUPPORT_STATUS::required)
202  {
203  if (argumentName != KIM::ARGUMENT_NAME::partialEnergy)
204  {
205  MY_ERROR("unsupported required argument");
206  }
207  }
208 
209  // must have energy
210  if (argumentName == KIM::ARGUMENT_NAME::partialEnergy)
211  {
212  if (! ((supportStatus == KIM::SUPPORT_STATUS::required)
213  ||
214  (supportStatus == KIM::SUPPORT_STATUS::optional)))
215  {
216  MY_ERROR("energy not available");
217  }
218  }
219  }
220 
221  // check call backs
222  int numberOfCallbacks;
223  KIM::CALLBACK_NAME::GetNumberOfCallbacks(&numberOfCallbacks);
224  for (int i=0; i<numberOfCallbacks; ++i)
225  {
226  KIM::CallbackName callbackName;
227  KIM::CALLBACK_NAME::GetCallbackName(i, &callbackName);
228  KIM::SupportStatus supportStatus;
229  kim_cluster_model->GetCallbackSupportStatus(callbackName, &supportStatus);
230 
231  std::cout << "Callback Name \""
232  << callbackName.String() << "\""
233  << " has supportStatus \""
234  << supportStatus.String() << "\""
235  << std::endl;
236 
237  // cannot handle any "required" call backs
238  if (supportStatus == KIM::SUPPORT_STATUS::required)
239  {
240  MY_ERROR("unsupported required call back");
241  }
242  }
243 
244  // We're compatible with the model. Let's do it.
245 
246  int numberOfParameters;
247  kim_cluster_model->GetNumberOfParameters(&numberOfParameters);
248  for (int i=0; i<numberOfParameters; ++i)
249  {
250  KIM::DataType dataType;
251  std::string str;
252  int extent;
253  kim_cluster_model->GetParameterDataTypeExtentAndDescription(
254  i, &dataType, &extent, &str);
255  std::cout << "Parameter No. " << i
256  << " has data type \"" << dataType.String() << "\""
257  << " with extent " << extent
258  << " and description : " << str << std::endl;
259  }
260 
261  error = kim_cluster_model->SetArgumentPointer(
262  KIM::ARGUMENT_NAME::numberOfParticles, (int *) &numberOfParticles_cluster)
263  || kim_cluster_model->SetArgumentPointer(
265  particleSpecies_cluster_model)
266  || kim_cluster_model->SetArgumentPointer(
268  particleContributing_cluster_model)
269  || kim_cluster_model->SetArgumentPointer(
270  KIM::ARGUMENT_NAME::coordinates, (double*) coords_cluster)
271  || kim_cluster_model->SetArgumentPointer(
272  KIM::ARGUMENT_NAME::partialEnergy, &energy_cluster_model);
273  if (error) MY_ERROR("KIM_API_set_data");
274  error = kim_cluster_model->SetCallbackPointer(
278  &nl_cluster_model);
279  if (error) MY_ERROR("set_call_back");
280 
281  kim_cluster_model->GetInfluenceDistance(&influence_distance_cluster_model);
282  kim_cluster_model->GetNeighborListCutoffsPointer(&number_of_cutoffs,
283  &cutoff_cluster_model);
284  if (number_of_cutoffs != 1) MY_ERROR("too many cutoffs");
285 
286  /* setup particleSpecies */
287  int isSpeciesSupported;
288  error = kim_cluster_model->GetSpeciesSupportAndCode(
290  &isSpeciesSupported,
291  &(particleSpecies_cluster_model[0]));
292  if (error) MY_ERROR("get_species_code");
293  for (i = 1; i < NCLUSTERPARTS; ++i)
294  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
295  /* setup particleContributing */
296  for (i = 0; i < NCLUSTERPARTS; ++i)
297  particleContributing_cluster_model[i] = 1; /* every particle contributes */
298 
299  /* setup neighbor lists */
300  /* allocate memory for list */
301  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
302  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
303  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("new unsuccessful");
304 
305  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
306  if (NULL==nl_cluster_model.neighborList) MY_ERROR("new unsuccessful");
307 
308  /* ready to compute */
309  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
310  std::cout << "--------------------------------------------------------------------------------\n";
311  std::cout << "This is Test : ex_test_Ar_fcc_cluster\n";
312  std::cout << "MODEL is : " << modelname << std::endl;
313 
314  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
315  {
316  /* update coordinates for cluster */
317  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
318  /* compute neighbor lists */
319  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
320  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
321 
322  /* call compute functions */
323  error = kim_cluster_model->Compute();
324  if (error) MY_ERROR("compute");
325 
326  /* print the results */
327  std::cout << "Energy for " << NCLUSTERPARTS << " parts = "
328  << std::setw(20) << energy_cluster_model
329  << std::setw(20) << CurrentSpacing
330  << std::endl;
331  }
332 
333 
334  /* call model destroy */
335  KIM::Model::Destroy(&kim_cluster_model);
336 
337  /* free memory of neighbor lists */
338  delete [] nl_cluster_model.NNeighbors;
339  delete [] nl_cluster_model.neighborList;
340 
341  /* everything is great */
342  return 0;
343 }
344 
345 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
346 {
347  /* local variables */
348  double FCCshifts[4][DIM];
349  double latVec[DIM];
350  int a;
351  int i;
352  int j;
353  int k;
354  int m;
355  int n;
356 
357  /* create a cubic FCC cluster of parts */
358  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
359  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
360  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
361  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
362 
363  a = 0;
364  for (i = 0; i < nCellsPerSide; ++i)
365  {
366  latVec[0] = ((double) i)*FCCspacing;
367  for (j = 0; j < nCellsPerSide; ++j)
368  {
369  latVec[1] = ((double) j)*FCCspacing;
370  for (k = 0; k < nCellsPerSide; ++k)
371  {
372  latVec[2] = ((double) k)*FCCspacing;
373  for (m = 0; m < 4; ++m)
374  {
375  for (n = 0; n < DIM; ++n)
376  {
377  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
378  }
379  a++;
380  }
381  }
382  /* add in the remaining three faces */
383  /* pos-x face */
384  latVec[0] = NCELLSPERSIDE*FCCspacing;
385  latVec[1] = ((double) i)*FCCspacing;
386  latVec[2] = ((double) j)*FCCspacing;
387  for (n = 0; n < DIM; ++n)
388  {
389  coords[a*DIM + n] = latVec[n];
390  }
391  a++;
392  for (n = 0; n < DIM; ++n)
393  {
394  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
395  }
396  a++;
397  /* pos-y face */
398  latVec[0] = ((double) i)*FCCspacing;
399  latVec[1] = NCELLSPERSIDE*FCCspacing;
400  latVec[2] = ((double) j)*FCCspacing;
401  for (n = 0; n < DIM; ++n)
402  {
403  coords[a*DIM + n] = latVec[n];
404  }
405  a++;
406  for (n = 0; n < DIM; ++n)
407  {
408  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
409  }
410  a++;
411  /* pos-z face */
412  latVec[0] = ((double) i)*FCCspacing;
413  latVec[1] = ((double) j)*FCCspacing;
414  latVec[2] = NCELLSPERSIDE*FCCspacing;
415  for (n = 0; n < DIM; ++n)
416  {
417  coords[a*DIM + n] = latVec[n];
418  }
419  a++;
420  for (n = 0; n < DIM; ++n)
421  {
422  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
423  }
424  a++;
425  }
426  /* add in the remaining three edges */
427  latVec[0] = ((double) i)*FCCspacing;
428  latVec[1] = NCELLSPERSIDE*FCCspacing;
429  latVec[2] = NCELLSPERSIDE*FCCspacing;
430  for (n = 0; n < DIM; ++n)
431  {
432  coords[a*DIM + n] = latVec[n];
433  }
434  a++;
435  latVec[0] = NCELLSPERSIDE*FCCspacing;
436  latVec[1] = ((double) i)*FCCspacing;
437  latVec[2] = NCELLSPERSIDE*FCCspacing;
438  for (n = 0; n < DIM; ++n)
439  {
440  coords[a*DIM + n] = latVec[n];
441  }
442  a++;
443  latVec[0] = NCELLSPERSIDE*FCCspacing;
444  latVec[1] = NCELLSPERSIDE*FCCspacing;
445  latVec[2] = ((double) i)*FCCspacing;
446  for (n = 0; n < DIM; ++n)
447  {
448  coords[a*DIM + n] = latVec[n];
449  }
450  a++;
451  }
452  /* add in the remaining corner */
453  for (n = 0; n < DIM; ++n)
454  {
455  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
456  }
457  a++;
458 
459  return;
460 }
461 
462 
463 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
464  double cutoff, NeighList* nl)
465 {
466  /* local variables */
467  int i;
468  int j;
469  int k;
470  int a;
471 
472  double dx[DIM];
473  double r2;
474  double cutoff2;
475 
476  cutoff2 = cutoff*cutoff;
477 
478  for (i = 0; i < numberOfParticles; ++i)
479  {
480  a = 0;
481  for (j = 0; j < numberOfParticles; ++j)
482  {
483  r2 = 0.0;
484  for (k = 0; k < DIM; ++k)
485  {
486  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
487  r2 += dx[k]*dx[k];
488  }
489 
490  if (r2 < cutoff2)
491  {
492  if ((half && i < j) || (!half && i != j))
493  {
494  /* part j is a neighbor of part i */
495  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
496  a++;
497  }
498  }
499  }
500  /* part i has `a' neighbors */
501  (*nl).NNeighbors[i] = a;
502  }
503 
504  return;
505 }
506 
507 int get_cluster_neigh(void const * const dataObject,
508  int const neighborListIndex, int const particleNumber,
509  int * const numberOfNeighbors,
510  int const ** const neighborsOfParticle)
511 {
512  /* local variables */
513  int error = true;
514  NeighList* nl = (NeighList*) dataObject;
515  int numberOfParticles = nl->numberOfParticles;
516 
517  if (neighborListIndex != 0) return error;
518 
519  /* initialize numNeigh */
520  *numberOfNeighbors = 0;
521 
522  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
523  {
524  MY_WARNING("Invalid part ID in get_cluster_neigh");
525  return true;
526  }
527 
528  /* set the returned number of neighbors for the returned part */
529  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
530 
531  /* set the location for the returned neighbor list */
532  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
533 
534  return false;
535 }
static int Create(Numbering const numbering, LengthUnit const requestedLengthUnit, EnergyUnit const requestedEnergyUnit, ChargeUnit const requestedChargeUnit, TemperatureUnit const requestedTemperatureUnit, TimeUnit const requestedTimeUnit, std::string const &modelName, int *const requestedUnitsAccepted, Model **const model)
void GetUnits(LengthUnit *const lengthUnit, EnergyUnit *const energyUnit, ChargeUnit *const chargeUnit, TemperatureUnit *const temperatureUnit, TimeUnit *const timeUnit) const
SpeciesName const Ar
int GetArgumentName(int const index, ArgumentName *const argumentName)
TimeUnit const ps
#define FCCSPACING
int Compute() const
ArgumentName const particleContributing
std::string String() const
std::string String() const
ArgumentName const numberOfParticles
void GetNeighborListCutoffsPointer(int *const numberOfCutoffs, double const **const cutoffs) const
int GetParameterDataTypeExtentAndDescription(int const index, DataType *const dataType, int *extent, std::string *const description) const
ChargeUnit const e
int GetArgumentDataType(ArgumentName const argumentName, DataType *const dataType)
std::string String() const
ArgumentName const partialEnergy
std::string String() const
SupportStatus const required
std::string String() const
std::string String() const
int SetCallbackPointer(CallbackName const callbackName, LanguageName const languageName, func *const fptr, void const *const dataObject)
ArgumentName const coordinates
void GetInfluenceDistance(double *const influenceDistance) const
void GetNumberOfCallbacks(int *const numberOfCallbacks)
int GetCallbackName(int const index, CallbackName *const callbackName)
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
LengthUnit const m
static void Destroy(Model **const model)
ArgumentName const particleSpeciesCodes
int SetArgumentPointer(ArgumentName const argumentName, int const *const ptr)
LanguageName const cpp
std::string String() const
std::string String() const
void GetNumberOfArguments(int *const numberOfArguments)
LengthUnit const A
#define NCELLSPERSIDE
int get_cluster_neigh(void const *const dataObject, int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle)
int GetArgumentSupportStatus(ArgumentName const argumentName, SupportStatus *const supportStatus) const
void GetNumberOfParameters(int *const numberOfParameters) const
#define NCLUSTERPARTS
EnergyUnit const eV
#define MY_WARNING(message)
int GetCallbackSupportStatus(CallbackName const callbackName, SupportStatus *const supportStatus) const
int GetSpeciesSupportAndCode(KIM::SpeciesName const speciesName, int *const speciesIsSupported, int *const code) const
#define MY_ERROR(message)
Numbering const zeroBased
CallbackName const GetNeighborList
std::string String() const
LogVerbosity const error
void() func()
Definition: KIM_func.hpp:40
void fcc_cluster_neighborlist(int half, int numberOfParticles, double *coords, double cutoff, NeighList *nl)
SupportStatus const optional
TemperatureUnit const K