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