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 numberOfNeighborLists,
87  double const * const cutoffs,
88  int const neighborListIndex, int const particleNumber,
89  int * const numberOfNeighbors,
90  int const ** const neighborsOfParticle);
91 
92 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
93 
94 
95 /* Main program */
96 int main()
97 {
98  /* Local variable declarations */
99  double const MinSpacing = 0.8*FCCSPACING;
100  double const MaxSpacing = 1.2*FCCSPACING;
101  double const SpacingIncr = 0.025*FCCSPACING;
102  double CurrentSpacing;
103  double cutpad = 0.75; /* Angstroms */
104  double force_norm;
105  int i;
106  int error;
107 
108 
109  /* model inputs */
110  int numberOfParticles_cluster = NCLUSTERPARTS;
111  int numberOfSpecies = 1;
112  int particleSpecies_cluster_model[NCLUSTERPARTS];
113  double coords_cluster[NCLUSTERPARTS][DIM];
114  NeighList nl_cluster_model;
115  /* model outputs */
116  double influence_distance_cluster_model;
117  int number_of_neighbor_lists;
118  double const * cutoff_cluster_model;
119  double energy_cluster_model;
120  double forces_cluster[NCLUSTERPARTS*DIM];
121 
122  std::string modelname;
123 
124  /* Get KIM Model names */
125  printf("Please enter valid KIM Model name: \n");
126  std::cin >> modelname;
127 
128 
129  /* initialize the model */
130  KIM::Model * kim_cluster_model;
131  int requestedUnitsAccepted;
138  modelname,
139  &requestedUnitsAccepted,
140  &kim_cluster_model);
141  if (error)
142  {
143  MY_ERROR("KIM::Model::Create()");
144  }
145 
146  // Check for compatibility with the model
147  if (!requestedUnitsAccepted)
148  {
149  MY_ERROR("Must Adapt to model units");
150  }
151 
152  // print model units
153  KIM::LengthUnit lengthUnit;
154  KIM::EnergyUnit energyUnit;
155  KIM::ChargeUnit chargeUnit;
156  KIM::TemperatureUnit temperatureUnit;
157  KIM::TimeUnit timeUnit;
158 
159  kim_cluster_model->GetUnits(&lengthUnit, &energyUnit, &chargeUnit,
160  &temperatureUnit, &timeUnit);
161 
162  std::cout << "LengthUnit is \"" << lengthUnit.String() << "\"" << std::endl
163  << "EnergyUnit is \"" << energyUnit.String() << "\"" << std::endl
164  << "ChargeUnit is \"" << chargeUnit.String() << "\"" << std::endl
165  << "TemperatureUnit is \"" << temperatureUnit.String()
166  << "\"" << std::endl
167  << "TimeUnit is \"" << timeUnit.String() << "\"" << std::endl;
168 
169  // check species
170  int speciesIsSupported;
171  int modelArCode;
172  error = kim_cluster_model->GetSpeciesSupportAndCode(
173  KIM::SPECIES_NAME::Ar, &speciesIsSupported, &modelArCode);
174  if ((error) || (!speciesIsSupported))
175  {
176  MY_ERROR("Species Ar not supported");
177  }
178 
179  KIM::ComputeArguments * computeArguments;
180  error = kim_cluster_model->ComputeArgumentsCreate(&computeArguments);
181  if (error)
182  {
183  MY_ERROR("Unable to create a ComputeArguments object.");
184  }
185 
186  // check compute arguments
187  int numberOfComputeArgumentNames;
189  &numberOfComputeArgumentNames);
190  for (int i=0; i<numberOfComputeArgumentNames; ++i)
191  {
192  KIM::ComputeArgumentName computeArgumentName;
193  KIM::SupportStatus supportStatus;
195  KIM::DataType dataType;
197  computeArgumentName, &dataType);
198  error = computeArguments->GetArgumentSupportStatus(computeArgumentName,
199  &supportStatus);
200  if (error)
201  MY_ERROR("unable to get ComputeArgument SupportStatus");
202 
203  std::cout << "ComputeArgument Name \""
204  << computeArgumentName.String() << "\""
205  << " is of type \""
206  << dataType.String() << "\""
207  << " and has supportStatus \""
208  << supportStatus.String() << "\""
209  << std::endl;
210 
211  // can only handle energy and force as a required arg
212  if (supportStatus == KIM::SUPPORT_STATUS::required)
213  {
214  if ((computeArgumentName != KIM::COMPUTE_ARGUMENT_NAME::partialEnergy)
215  &&
216  (computeArgumentName != KIM::COMPUTE_ARGUMENT_NAME::partialForces))
217  {
218  MY_ERROR("unsupported required ComputeArgument");
219  }
220  }
221 
222  // must have energy and forces
223  if ((computeArgumentName == KIM::COMPUTE_ARGUMENT_NAME::partialEnergy)
224  ||
225  (computeArgumentName == KIM::COMPUTE_ARGUMENT_NAME::partialForces))
226  {
227  if (! ((supportStatus == KIM::SUPPORT_STATUS::required)
228  ||
229  (supportStatus == KIM::SUPPORT_STATUS::optional)))
230  {
231  MY_ERROR("energy or forces not available");
232  }
233  }
234  }
235 
236  // check compute callbacks
237  int numberOfComputeCallbackNames;
239  &numberOfComputeCallbackNames);
240  for (int i=0; i<numberOfComputeCallbackNames; ++i)
241  {
242  KIM::ComputeCallbackName computeCallbackName;
244  KIM::SupportStatus supportStatus;
245  computeArguments->GetCallbackSupportStatus(computeCallbackName,
246  &supportStatus);
247 
248  std::cout << "ComputeCallback Name \""
249  << computeCallbackName.String() << "\""
250  << " has supportStatus \""
251  << supportStatus.String() << "\""
252  << std::endl;
253 
254  // cannot handle any "required" callbacks
255  if (supportStatus == KIM::SUPPORT_STATUS::required)
256  {
257  MY_ERROR("unsupported required ComputeCallback");
258  }
259  }
260 
261  int numberOfParameters;
262  kim_cluster_model->GetNumberOfParameters(&numberOfParameters);
263  for (int i=0; i<numberOfParameters; ++i)
264  {
265  KIM::DataType dataType;
266  std::string const * str;
267  int extent;
268  kim_cluster_model->GetParameterDataTypeExtentAndDescription(
269  i, &dataType, &extent, &str);
270  std::cout << "Parameter No. " << i
271  << " has data type \"" << dataType.String() << "\""
272  << " with extent " << extent
273  << " and description : " << *str << std::endl;
274  }
275 
276  // We're compatible with the model. Let's do it.
277 
278  /* setup particleSpecies */
279  particleSpecies_cluster_model[0] = kim_cluster_model.get_species_code("Ar", &error);
280  if (KIM_ERROR_OK > error) MY_ERROR("get_species_code");
281  for (i = 1; i < NCLUSTERPARTS; ++i)
282  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
283 
284  /* setup neighbor lists */
285  /* allocate memory for list */
286  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
287  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
288  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("new unsuccessful");
289 
290  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
291  if (NULL==nl_cluster_model.neighborList) MY_ERROR("new unsuccessful");
292 
293  /* ready to compute */
294  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
295  std::cout << "This is Test : ex_test_Ar_fcc_cluster_cpp\n";
296  std::cout << "--------------------------------------------------------------------------------\n";
297  std::cout << "Results for KIM Model : " << modelname << std::endl;
298 
299  std::cout << std::setw(20) << "Energy"
300  << std::setw(20) << "Force Norm"
301  << std::setw(20) << "Lattice Spacing"
302  << std::endl;
303  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
304  {
305  /* update coordinates for cluster */
306  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
307  /* compute neighbor lists */
308  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
309  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
310 
311  /* call compute functions */
312  error = kim_cluster_model.model_compute();
313  if (KIM_ERROR_OK > error) MY_ERROR("compute");
314 
315  /* compute force norm */
316  force_norm = 0.0;
317  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
318  {
319  force_norm += forces_cluster[i]*forces_cluster[i];
320  }
321  force_norm = sqrt(force_norm);
322 
323  /* print the results */
324  std::cout << std::setw(20) << energy_cluster_model
325  << std::setw(20) << force_norm
326  << std::setw(20) << CurrentSpacing
327  << std::endl;
328  }
329 
330 
331  /* call model destroy */
332  error = kim_cluster_model.model_destroy();
333  if (KIM_ERROR_OK > error) MY_ERROR("destroy");
334 
335  /* free memory of neighbor lists */
336  delete [] nl_cluster_model.NNeighbors;
337  delete [] nl_cluster_model.neighborList;
338 
339  /* everything is great */
340  return 0;
341 }
342 
343 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
344 {
345  /* local variables */
346  double FCCshifts[4][DIM];
347  double latVec[DIM];
348  int a;
349  int i;
350  int j;
351  int k;
352  int m;
353  int n;
354 
355  /* create a cubic FCC cluster of parts */
356  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
357  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
358  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
359  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
360 
361  a = 0;
362  for (i = 0; i < nCellsPerSide; ++i)
363  {
364  latVec[0] = ((double) i)*FCCspacing;
365  for (j = 0; j < nCellsPerSide; ++j)
366  {
367  latVec[1] = ((double) j)*FCCspacing;
368  for (k = 0; k < nCellsPerSide; ++k)
369  {
370  latVec[2] = ((double) k)*FCCspacing;
371  for (m = 0; m < 4; ++m)
372  {
373  for (n = 0; n < DIM; ++n)
374  {
375  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
376  }
377  a++;
378  }
379  }
380  /* add in the remaining three faces */
381  /* pos-x face */
382  latVec[0] = NCELLSPERSIDE*FCCspacing;
383  latVec[1] = ((double) i)*FCCspacing;
384  latVec[2] = ((double) j)*FCCspacing;
385  for (n = 0; n < DIM; ++n)
386  {
387  coords[a*DIM + n] = latVec[n];
388  }
389  a++;
390  for (n = 0; n < DIM; ++n)
391  {
392  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
393  }
394  a++;
395  /* pos-y face */
396  latVec[0] = ((double) i)*FCCspacing;
397  latVec[1] = NCELLSPERSIDE*FCCspacing;
398  latVec[2] = ((double) j)*FCCspacing;
399  for (n = 0; n < DIM; ++n)
400  {
401  coords[a*DIM + n] = latVec[n];
402  }
403  a++;
404  for (n = 0; n < DIM; ++n)
405  {
406  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
407  }
408  a++;
409  /* pos-z face */
410  latVec[0] = ((double) i)*FCCspacing;
411  latVec[1] = ((double) j)*FCCspacing;
412  latVec[2] = NCELLSPERSIDE*FCCspacing;
413  for (n = 0; n < DIM; ++n)
414  {
415  coords[a*DIM + n] = latVec[n];
416  }
417  a++;
418  for (n = 0; n < DIM; ++n)
419  {
420  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
421  }
422  a++;
423  }
424  /* add in the remaining three edges */
425  latVec[0] = ((double) i)*FCCspacing;
426  latVec[1] = NCELLSPERSIDE*FCCspacing;
427  latVec[2] = NCELLSPERSIDE*FCCspacing;
428  for (n = 0; n < DIM; ++n)
429  {
430  coords[a*DIM + n] = latVec[n];
431  }
432  a++;
433  latVec[0] = NCELLSPERSIDE*FCCspacing;
434  latVec[1] = ((double) i)*FCCspacing;
435  latVec[2] = NCELLSPERSIDE*FCCspacing;
436  for (n = 0; n < DIM; ++n)
437  {
438  coords[a*DIM + n] = latVec[n];
439  }
440  a++;
441  latVec[0] = NCELLSPERSIDE*FCCspacing;
442  latVec[1] = NCELLSPERSIDE*FCCspacing;
443  latVec[2] = ((double) i)*FCCspacing;
444  for (n = 0; n < DIM; ++n)
445  {
446  coords[a*DIM + n] = latVec[n];
447  }
448  a++;
449  }
450  /* add in the remaining corner */
451  for (n = 0; n < DIM; ++n)
452  {
453  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
454  }
455  a++;
456 
457  return;
458 }
459 
460 
461 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
462  double cutoff, NeighList* nl)
463 {
464  /* local variables */
465  int i;
466  int j;
467  int k;
468  int a;
469 
470  double dx[DIM];
471  double r2;
472  double cutoff2;
473 
474  nl->cutoff = cutoff;
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 numberOfNeighborLists,
509  double const * const cutoffs,
510  int const neighborListIndex, int const particleNumber,
511  int * const numberOfNeighbors,
512  int const ** const neighborsOfParticle)
513 {
514  /* local variables */
515  int error = true;
516  NeighList* nl = (NeighList*) dataObject;
518 
519  if ((numberOfNeighborLists != 1) || (cutoffs[0] > nl->cutoff)) return error;
520 
521  if (neighborListIndex != 0) return error;
522 
523  /* initialize numNeigh */
524  *numberOfNeighbors = 0;
525 
526  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
527  {
528  MY_WARNING("Invalid part ID in get_cluster_neigh");
529  return true;
530  }
531 
532  /* set the returned number of neighbors for the returned part */
533  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
534 
535  /* set the location for the returned neighbor list */
536  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
537 
538  return false;
539 }
#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