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 <stdio.h>
35 #include <iostream>
36 #include <iomanip>
37 #include <string>
38 #include <cmath>
39 #include "KIM_SimulatorHeaders.hpp"
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 #define MY_ERROR(message) \
51  { \
52  std::cout << "* Error : \"" << message << "\" : " \
53  << __LINE__ << ":" << __FILE__ << std::endl; \
54  exit(1); \
55  }
56 
57 #define MY_WARNING(message) \
58  { \
59  std::cout << "* Error : \"" << message << "\" : " \
60  << __LINE__ << ":" << __FILE__ << std::endl; \
61  }
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 
77 int get_cluster_neigh(void const * const dataObject,
78  int const numberOfNeighborLists,
79  double const * const cutoffs,
80  int const neighborListIndex, int const particleNumber,
81  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  /* model inputs */
102  int numberOfParticles_cluster = NCLUSTERPARTS;
103  int particleSpecies_cluster_model[NCLUSTERPARTS];
104  int particleContributing_cluster_model[NCLUSTERPARTS];
105  double coords_cluster[NCLUSTERPARTS][DIM];
106  NeighList nl_cluster_model;
107  /* model outputs */
108  double influence_distance_cluster_model;
109  int number_of_neighbor_lists;
110  double const * cutoff_cluster_model;
111  double energy_cluster_model;
112  double forces_cluster[NCLUSTERPARTS*DIM];
113 
114  std::string modelname;
115 
116  /* Get KIM Model names */
117  printf("Please enter valid KIM Model name: \n");
118  std::cin >> modelname;
119 
120 
121  /* initialize the model */
122  KIM::Model * kim_cluster_model;
123  int requestedUnitsAccepted;
130  modelname,
131  &requestedUnitsAccepted,
132  &kim_cluster_model);
133  if (error)
134  {
135  MY_ERROR("KIM::Model::Create()");
136  }
137 
138  // Check for compatibility with the model
139  if (!requestedUnitsAccepted)
140  {
141  MY_ERROR("Must Adapt to model units");
142  }
143 
144  // print model units
145  KIM::LengthUnit lengthUnit;
146  KIM::EnergyUnit energyUnit;
147  KIM::ChargeUnit chargeUnit;
148  KIM::TemperatureUnit temperatureUnit;
149  KIM::TimeUnit timeUnit;
150 
151  kim_cluster_model->GetUnits(&lengthUnit, &energyUnit, &chargeUnit,
152  &temperatureUnit, &timeUnit);
153 
154  std::cout << "LengthUnit is \"" << lengthUnit.String() << "\"" << std::endl
155  << "EnergyUnit is \"" << energyUnit.String() << "\"" << std::endl
156  << "ChargeUnit is \"" << chargeUnit.String() << "\"" << std::endl
157  << "TemperatureUnit is \"" << temperatureUnit.String()
158  << "\"" << std::endl
159  << "TimeUnit is \"" << timeUnit.String() << "\"" << std::endl;
160 
161  // check species
162  int speciesIsSupported;
163  int modelArCode;
164  error = kim_cluster_model->GetSpeciesSupportAndCode(
165  KIM::SPECIES_NAME::Ar, &speciesIsSupported, &modelArCode);
166  if ((error) || (!speciesIsSupported))
167  {
168  MY_ERROR("Species Ar not supported");
169  }
170 
171  KIM::ComputeArguments * computeArguments;
172  error = kim_cluster_model->ComputeArgumentsCreate(&computeArguments);
173  if (error)
174  {
175  MY_ERROR("Unable to create a ComputeArguments object.");
176  }
177 
178  // check compute arguments
179  int numberOfComputeArgumentNames;
181  &numberOfComputeArgumentNames);
182  for (int i=0; i<numberOfComputeArgumentNames; ++i)
183  {
184  KIM::ComputeArgumentName computeArgumentName;
185  KIM::SupportStatus supportStatus;
187  KIM::DataType dataType;
189  computeArgumentName, &dataType);
190  error = computeArguments->GetArgumentSupportStatus(computeArgumentName,
191  &supportStatus);
192  if (error)
193  MY_ERROR("unable to get ComputeArgument SupportStatus");
194 
195  std::cout << "ComputeArgument Name \""
196  << computeArgumentName.String() << "\""
197  << " is of type \""
198  << dataType.String() << "\""
199  << " and has supportStatus \""
200  << supportStatus.String() << "\""
201  << std::endl;
202 
203  // can only handle energy and force as a required arg
204  if (supportStatus == KIM::SUPPORT_STATUS::required)
205  {
206  if ((computeArgumentName != KIM::COMPUTE_ARGUMENT_NAME::partialEnergy)
207  &&
208  (computeArgumentName != KIM::COMPUTE_ARGUMENT_NAME::partialForces))
209  {
210  MY_ERROR("unsupported required ComputeArgument");
211  }
212  }
213 
214  // must have energy and forces
215  if ((computeArgumentName == KIM::COMPUTE_ARGUMENT_NAME::partialEnergy)
216  ||
217  (computeArgumentName == KIM::COMPUTE_ARGUMENT_NAME::partialForces))
218  {
219  if (! ((supportStatus == KIM::SUPPORT_STATUS::required)
220  ||
221  (supportStatus == KIM::SUPPORT_STATUS::optional)))
222  {
223  MY_ERROR("energy or forces not available");
224  }
225  }
226  }
227 
228  // check compute callbacks
229  int numberOfComputeCallbackNames;
231  &numberOfComputeCallbackNames);
232  for (int i=0; i<numberOfComputeCallbackNames; ++i)
233  {
234  KIM::ComputeCallbackName computeCallbackName;
236  KIM::SupportStatus supportStatus;
237  computeArguments->GetCallbackSupportStatus(computeCallbackName,
238  &supportStatus);
239 
240  std::cout << "ComputeCallback Name \""
241  << computeCallbackName.String() << "\""
242  << " has supportStatus \""
243  << supportStatus.String() << "\""
244  << std::endl;
245 
246  // cannot handle any "required" callbacks
247  if (supportStatus == KIM::SUPPORT_STATUS::required)
248  {
249  MY_ERROR("unsupported required ComputeCallback");
250  }
251  }
252 
253  int numberOfParameters;
254  kim_cluster_model->GetNumberOfParameters(&numberOfParameters);
255  for (int i=0; i<numberOfParameters; ++i)
256  {
257  KIM::DataType dataType;
258  std::string const * str;
259  int extent;
260  kim_cluster_model->GetParameterDataTypeExtentAndDescription(
261  i, &dataType, &extent, &str);
262  std::cout << "Parameter No. " << i
263  << " has data type \"" << dataType.String() << "\""
264  << " with extent " << extent
265  << " and description : " << *str << std::endl;
266  }
267 
268  // We're compatible with the model. Let's do it.
269 
270  error = computeArguments->SetArgumentPointer(
272  (int *) &numberOfParticles_cluster)
273  || computeArguments->SetArgumentPointer(
275  particleSpecies_cluster_model)
276  || computeArguments->SetArgumentPointer(
278  particleContributing_cluster_model)
279  || computeArguments->SetArgumentPointer(
280  KIM::COMPUTE_ARGUMENT_NAME::coordinates, (double*) coords_cluster)
281  || computeArguments->SetArgumentPointer(
282  KIM::COMPUTE_ARGUMENT_NAME::partialEnergy, &energy_cluster_model)
283  || computeArguments->SetArgumentPointer(
284  KIM::COMPUTE_ARGUMENT_NAME::partialForces, (double*) forces_cluster);
285  if (error) MY_ERROR("KIM_API_set_data");
286  error = computeArguments->SetCallbackPointer(
290  &nl_cluster_model);
291  if (error) MY_ERROR("set_call_back");
292 
293  kim_cluster_model->GetInfluenceDistance(&influence_distance_cluster_model);
294  int const * paddingNeighborHints;
295  int const * halfListHints;
296  kim_cluster_model->GetNeighborListPointers(&number_of_neighbor_lists,
297  &cutoff_cluster_model,
298  &paddingNeighborHints,
299  &halfListHints);
300  std::cout << "Model has influence distance of : "
301  << influence_distance_cluster_model << std::endl;
302  std::cout << "Model has numberOfNeighborLists : " << number_of_neighbor_lists
303  << std::endl;
304  for (int i=0; i<number_of_neighbor_lists; ++i)
305  {
306  std::cout << "\t" << "Neighbor list " << i << " has cutoff "
307  << cutoff_cluster_model[i] << " with paddingNeighborHint "
308  << paddingNeighborHints[i] << " and halfListHint "
309  << halfListHints[i]
310  << std::endl;
311  }
312  // ignoring hints from here on...
313  if (number_of_neighbor_lists != 1) MY_ERROR("too many neighbor lists");
314 
315  /* setup particleSpecies */
316  int isSpeciesSupported;
317  error = kim_cluster_model->GetSpeciesSupportAndCode(
319  &isSpeciesSupported,
320  &(particleSpecies_cluster_model[0]));
321  if (error) MY_ERROR("get_species_code");
322  for (i = 1; i < NCLUSTERPARTS; ++i)
323  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
324  /* setup particleContributing */
325  for (i = 0; i < NCLUSTERPARTS; ++i)
326  particleContributing_cluster_model[i] = 1; /* every particle contributes */
327 
328  /* setup neighbor lists */
329  /* allocate memory for list */
330  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
331  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
332  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("new unsuccessful");
333 
334  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
335  if (NULL==nl_cluster_model.neighborList) MY_ERROR("new unsuccessful");
336 
337  /* ready to compute */
338  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
339  std::cout << "This is Test : ex_test_Ar_fcc_cluster_cpp\n";
340  std::cout << "--------------------------------------------------------------------------------\n";
341  std::cout << "Results for KIM Model : " << modelname << std::endl;
342 
343  std::cout << std::setw(20) << "Energy"
344  << std::setw(20) << "Force Norm"
345  << std::setw(20) << "Lattice Spacing"
346  << std::endl;
347  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
348  {
349  /* update coordinates for cluster */
350  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
351  /* compute neighbor lists */
352  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
353  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
354 
355  /* call compute functions */
356  error = kim_cluster_model->Compute(computeArguments);
357  if (error) MY_ERROR("compute");
358 
359  /* compute force norm */
360  force_norm = 0.0;
361  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
362  {
363  force_norm += forces_cluster[i]*forces_cluster[i];
364  }
365  force_norm = sqrt(force_norm);
366 
367  /* print the results */
368  std::cout << std::setw(20) << energy_cluster_model
369  << std::setw(20) << force_norm
370  << std::setw(20) << CurrentSpacing
371  << std::endl;
372  }
373 
374  /* call compute arguments destroy */
375  error = kim_cluster_model->ComputeArgumentsDestroy(&computeArguments);
376  if (error)
377  {
378  MY_ERROR("Unable to destroy compute arguments");
379  }
380 
381  /* call model destroy */
382  KIM::Model::Destroy(&kim_cluster_model);
383 
384  /* free memory of neighbor lists */
385  delete [] nl_cluster_model.NNeighbors;
386  delete [] nl_cluster_model.neighborList;
387 
388  /* everything is great */
389  return 0;
390 }
391 
392 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
393 {
394  /* local variables */
395  double FCCshifts[4][DIM];
396  double latVec[DIM];
397  int a;
398  int i;
399  int j;
400  int k;
401  int m;
402  int n;
403 
404  /* create a cubic FCC cluster of parts */
405  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
406  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
407  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
408  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
409 
410  a = 0;
411  for (i = 0; i < nCellsPerSide; ++i)
412  {
413  latVec[0] = ((double) i)*FCCspacing;
414  for (j = 0; j < nCellsPerSide; ++j)
415  {
416  latVec[1] = ((double) j)*FCCspacing;
417  for (k = 0; k < nCellsPerSide; ++k)
418  {
419  latVec[2] = ((double) k)*FCCspacing;
420  for (m = 0; m < 4; ++m)
421  {
422  for (n = 0; n < DIM; ++n)
423  {
424  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
425  }
426  a++;
427  }
428  }
429  /* add in the remaining three faces */
430  /* pos-x face */
431  latVec[0] = NCELLSPERSIDE*FCCspacing;
432  latVec[1] = ((double) i)*FCCspacing;
433  latVec[2] = ((double) j)*FCCspacing;
434  for (n = 0; n < DIM; ++n)
435  {
436  coords[a*DIM + n] = latVec[n];
437  }
438  a++;
439  for (n = 0; n < DIM; ++n)
440  {
441  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
442  }
443  a++;
444  /* pos-y face */
445  latVec[0] = ((double) i)*FCCspacing;
446  latVec[1] = NCELLSPERSIDE*FCCspacing;
447  latVec[2] = ((double) j)*FCCspacing;
448  for (n = 0; n < DIM; ++n)
449  {
450  coords[a*DIM + n] = latVec[n];
451  }
452  a++;
453  for (n = 0; n < DIM; ++n)
454  {
455  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
456  }
457  a++;
458  /* pos-z face */
459  latVec[0] = ((double) i)*FCCspacing;
460  latVec[1] = ((double) j)*FCCspacing;
461  latVec[2] = NCELLSPERSIDE*FCCspacing;
462  for (n = 0; n < DIM; ++n)
463  {
464  coords[a*DIM + n] = latVec[n];
465  }
466  a++;
467  for (n = 0; n < DIM; ++n)
468  {
469  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
470  }
471  a++;
472  }
473  /* add in the remaining three edges */
474  latVec[0] = ((double) i)*FCCspacing;
475  latVec[1] = NCELLSPERSIDE*FCCspacing;
476  latVec[2] = NCELLSPERSIDE*FCCspacing;
477  for (n = 0; n < DIM; ++n)
478  {
479  coords[a*DIM + n] = latVec[n];
480  }
481  a++;
482  latVec[0] = NCELLSPERSIDE*FCCspacing;
483  latVec[1] = ((double) i)*FCCspacing;
484  latVec[2] = NCELLSPERSIDE*FCCspacing;
485  for (n = 0; n < DIM; ++n)
486  {
487  coords[a*DIM + n] = latVec[n];
488  }
489  a++;
490  latVec[0] = NCELLSPERSIDE*FCCspacing;
491  latVec[1] = NCELLSPERSIDE*FCCspacing;
492  latVec[2] = ((double) i)*FCCspacing;
493  for (n = 0; n < DIM; ++n)
494  {
495  coords[a*DIM + n] = latVec[n];
496  }
497  a++;
498  }
499  /* add in the remaining corner */
500  for (n = 0; n < DIM; ++n)
501  {
502  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
503  }
504  a++;
505 
506  return;
507 }
508 
509 
510 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
511  double cutoff, NeighList* nl)
512 {
513  /* local variables */
514  int i;
515  int j;
516  int k;
517  int a;
518 
519  double dx[DIM];
520  double r2;
521  double cutoff2;
522 
523  nl->cutoff = cutoff;
524 
525  cutoff2 = cutoff*cutoff;
526 
527  for (i = 0; i < numberOfParticles; ++i)
528  {
529  a = 0;
530  for (j = 0; j < numberOfParticles; ++j)
531  {
532  r2 = 0.0;
533  for (k = 0; k < DIM; ++k)
534  {
535  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
536  r2 += dx[k]*dx[k];
537  }
538 
539  if (r2 < cutoff2)
540  {
541  if ((half && i < j) || (!half && i != j))
542  {
543  /* part j is a neighbor of part i */
544  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
545  a++;
546  }
547  }
548  }
549  /* part i has `a' neighbors */
550  (*nl).NNeighbors[i] = a;
551  }
552 
553  return;
554 }
555 
556 int get_cluster_neigh(void const * const dataObject,
557  int const numberOfNeighborLists,
558  double const * const cutoffs,
559  int const neighborListIndex, int const particleNumber,
560  int * const numberOfNeighbors,
561  int const ** const neighborsOfParticle)
562 {
563  /* local variables */
564  int error = true;
565  NeighList* nl = (NeighList*) dataObject;
567 
568  if ((numberOfNeighborLists != 1) || (cutoffs[0] > nl->cutoff)) return error;
569 
570  if (neighborListIndex != 0) return error;
571 
572  /* initialize numNeigh */
573  *numberOfNeighbors = 0;
574 
575  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
576  {
577  MY_WARNING("Invalid part ID in get_cluster_neigh");
578  return true;
579  }
580 
581  /* set the returned number of neighbors for the returned part */
582  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
583 
584  /* set the location for the returned neighbor list */
585  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
586 
587  return false;
588 }
int SetArgumentPointer(ComputeArgumentName const computeArgumentName, int const *const ptr)
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
ComputeArgumentName const coordinates
void GetNumberOfComputeCallbackNames(int *const numberOfComputeCallbackNames)
#define MY_ERROR(message)
ChargeUnit const e
std::string const & String() const
#define MY_WARNING(message)
void GetNumberOfComputeArgumentNames(int *const numberOfComputeArgumentNames)
void GetNeighborListPointers(int *const numberOfNeighborLists, double const **const cutoffs, int const **const paddingNeighborHints, int const **const halfListHints) const
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
int ComputeArgumentsCreate(ComputeArguments **const computeArguments) const
ComputeCallbackName const GetNeighborList
SupportStatus const required
int Compute(ComputeArguments const *const computeArguments) const
std::string const & String() const
#define NCLUSTERPARTS
void GetInfluenceDistance(double *const influenceDistance) const
ComputeArgumentName const particleContributing
int GetCallbackSupportStatus(ComputeCallbackName const computeCallbackName, SupportStatus *const supportStatus) const
ComputeArgumentName const partialEnergy
LengthUnit const m
static void Destroy(Model **const model)
LanguageName const cpp
int GetArgumentSupportStatus(ComputeArgumentName const computeArgumentName, SupportStatus *const supportStatus) const
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
int SetCallbackPointer(ComputeCallbackName const computeCallbackName, LanguageName const languageName, func *const fptr, void const *const dataObject)
std::string const & String() const
ComputeArgumentName const particleSpeciesCodes
#define FCCSPACING
ComputeArgumentName const numberOfParticles
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
#define NCELLSPERSIDE
int GetSpeciesSupportAndCode(SpeciesName const speciesName, int *const speciesIsSupported, int *const code) const
int ComputeArgumentsDestroy(ComputeArguments **const computeArguments) 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
void() func()
Definition: KIM_func.hpp:40
SupportStatus const optional
TemperatureUnit const K