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