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 char const * const descriptor();
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 numberOfNeighborLists,
88  double const * const cutoffs,
89  int const neighborListIndex, int const particleNumber,
90  int * const numberOfNeighbors,
91  int const ** const neighborsOfParticle);
92 
93 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords);
94 
95 
96 /* Main program */
97 int main()
98 {
99  /* Local variable declarations */
100  double const MinSpacing = 0.8*FCCSPACING;
101  double const MaxSpacing = 1.2*FCCSPACING;
102  double const SpacingIncr = 0.025*FCCSPACING;
103  double CurrentSpacing;
104  double cutpad = 0.75; /* Angstroms */
105  double force_norm;
106  int i;
107  int error;
108 
109 
110  /* model inputs */
111  int numberOfParticles_cluster = NCLUSTERPARTS;
112  int numberOfSpecies = 1;
113  int particleSpecies_cluster_model[NCLUSTERPARTS];
114  double coords_cluster[NCLUSTERPARTS][DIM];
115  NeighList nl_cluster_model;
116  /* model outputs */
117  double 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_API_model kim_cluster_model;
130  error = kim_cluster_model.string_init(descriptor(), modelname.c_str());
131  if (KIM_ERROR_OK > error)
132  MY_ERROR("KIM_API_string_init()");
133 
134  kim_cluster_model.setm_data(&error, 8*4,
135  "numberOfParticles", 1, &numberOfParticles_cluster, 1,
136  "numberOfSpecies", 1, &numberOfSpecies, 1,
137  "particleSpecies", numberOfParticles_cluster, &particleSpecies_cluster_model, 1,
138  "coordinates", DIM*numberOfParticles_cluster, coords_cluster, 1,
139  "neighObject", 1, &nl_cluster_model, 1,
140  "cutoff", 1, &cutoff_cluster_model, 1,
141  "energy", 1, &energy_cluster_model, 1,
142  "forces", DIM*numberOfParticles_cluster, forces_cluster, 1);
143  if (KIM_ERROR_OK > error) MY_ERROR("KIM_API_setm_data");
144  error = kim_cluster_model.set_method("get_neigh", 1, (func_ptr) &get_cluster_neigh);
145  if (KIM_ERROR_OK > error) MY_ERROR("KIM_API_set_method");
146 
147  /* call model init routine */
148  error = kim_cluster_model.model_init();
149  if (KIM_ERROR_OK > error) MY_ERROR("KIM_API_model_init");
150 
151  /* setup particleSpecies */
152  particleSpecies_cluster_model[0] = kim_cluster_model.get_species_code("Ar", &error);
153  if (KIM_ERROR_OK > error) MY_ERROR("get_species_code");
154  for (i = 1; i < NCLUSTERPARTS; ++i)
155  particleSpecies_cluster_model[i] = particleSpecies_cluster_model[0];
156 
157  /* setup neighbor lists */
158  /* allocate memory for list */
159  nl_cluster_model.numberOfParticles = NCLUSTERPARTS;
160  nl_cluster_model.NNeighbors = new int[NCLUSTERPARTS];
161  if (NULL==nl_cluster_model.NNeighbors) MY_ERROR("new unsuccessful");
162 
163  nl_cluster_model.neighborList = new int[NCLUSTERPARTS*NCLUSTERPARTS];
164  if (NULL==nl_cluster_model.neighborList) MY_ERROR("new unsuccessful");
165 
166  /* ready to compute */
167  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(10);
168  std::cout << "This is Test : ex_test_Ar_fcc_cluster_cpp\n";
169  std::cout << "--------------------------------------------------------------------------------\n";
170  std::cout << "Results for KIM Model : " << modelname << std::endl;
171 
172  std::cout << std::setw(20) << "Energy"
173  << std::setw(20) << "Force Norm"
174  << std::setw(20) << "Lattice Spacing"
175  << std::endl;
176  for (CurrentSpacing = MinSpacing; CurrentSpacing < MaxSpacing; CurrentSpacing += SpacingIncr)
177  {
178  /* update coordinates for cluster */
179  create_FCC_cluster(CurrentSpacing, NCELLSPERSIDE, &(coords_cluster[0][0]));
180  /* compute neighbor lists */
181  fcc_cluster_neighborlist(0, NCLUSTERPARTS, &(coords_cluster[0][0]),
182  (*cutoff_cluster_model + cutpad), &nl_cluster_model);
183 
184  /* call compute functions */
185  error = kim_cluster_model.model_compute();
186  if (KIM_ERROR_OK > error) MY_ERROR("compute");
187 
188  /* compute force norm */
189  force_norm = 0.0;
190  for (i=0; i < DIM*numberOfParticles_cluster; ++i)
191  {
192  force_norm += forces_cluster[i]*forces_cluster[i];
193  }
194  force_norm = sqrt(force_norm);
195 
196  /* print the results */
197  std::cout << std::setw(20) << energy_cluster_model
198  << std::setw(20) << force_norm
199  << std::setw(20) << CurrentSpacing
200  << std::endl;
201  }
202 
203 
204  /* call model destroy */
205  error = kim_cluster_model.model_destroy();
206  if (KIM_ERROR_OK > error) MY_ERROR("destroy");
207 
208  /* free memory of neighbor lists */
209  delete [] nl_cluster_model.NNeighbors;
210  delete [] nl_cluster_model.neighborList;
211 
212  /* everything is great */
213  return 0;
214 }
215 
216 void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
217 {
218  /* local variables */
219  double FCCshifts[4][DIM];
220  double latVec[DIM];
221  int a;
222  int i;
223  int j;
224  int k;
225  int m;
226  int n;
227 
228  /* create a cubic FCC cluster of parts */
229  FCCshifts[0][0] = 0.0; FCCshifts[0][1] = 0.0; FCCshifts[0][2] = 0.0;
230  FCCshifts[1][0] = 0.5*FCCspacing; FCCshifts[1][1] = 0.5*FCCspacing; FCCshifts[1][2] = 0.0;
231  FCCshifts[2][0] = 0.5*FCCspacing; FCCshifts[2][1] = 0.0; FCCshifts[2][2] = 0.5*FCCspacing;
232  FCCshifts[3][0] = 0.0; FCCshifts[3][1] = 0.5*FCCspacing; FCCshifts[3][2] = 0.5*FCCspacing;
233 
234  a = 0;
235  for (i = 0; i < nCellsPerSide; ++i)
236  {
237  latVec[0] = ((double) i)*FCCspacing;
238  for (j = 0; j < nCellsPerSide; ++j)
239  {
240  latVec[1] = ((double) j)*FCCspacing;
241  for (k = 0; k < nCellsPerSide; ++k)
242  {
243  latVec[2] = ((double) k)*FCCspacing;
244  for (m = 0; m < 4; ++m)
245  {
246  for (n = 0; n < DIM; ++n)
247  {
248  coords[a*DIM + n] = latVec[n] + FCCshifts[m][n];
249  }
250  a++;
251  }
252  }
253  /* add in the remaining three faces */
254  /* pos-x face */
255  latVec[0] = NCELLSPERSIDE*FCCspacing;
256  latVec[1] = ((double) i)*FCCspacing;
257  latVec[2] = ((double) j)*FCCspacing;
258  for (n = 0; n < DIM; ++n)
259  {
260  coords[a*DIM + n] = latVec[n];
261  }
262  a++;
263  for (n = 0; n < DIM; ++n)
264  {
265  coords[a*DIM + n] = latVec[n] + FCCshifts[3][n];
266  }
267  a++;
268  /* pos-y face */
269  latVec[0] = ((double) i)*FCCspacing;
270  latVec[1] = NCELLSPERSIDE*FCCspacing;
271  latVec[2] = ((double) j)*FCCspacing;
272  for (n = 0; n < DIM; ++n)
273  {
274  coords[a*DIM + n] = latVec[n];
275  }
276  a++;
277  for (n = 0; n < DIM; ++n)
278  {
279  coords[a*DIM + n] = latVec[n] + FCCshifts[2][n];
280  }
281  a++;
282  /* pos-z face */
283  latVec[0] = ((double) i)*FCCspacing;
284  latVec[1] = ((double) j)*FCCspacing;
285  latVec[2] = NCELLSPERSIDE*FCCspacing;
286  for (n = 0; n < DIM; ++n)
287  {
288  coords[a*DIM + n] = latVec[n];
289  }
290  a++;
291  for (n = 0; n < DIM; ++n)
292  {
293  coords[a*DIM + n] = latVec[n] + FCCshifts[1][n];
294  }
295  a++;
296  }
297  /* add in the remaining three edges */
298  latVec[0] = ((double) i)*FCCspacing;
299  latVec[1] = NCELLSPERSIDE*FCCspacing;
300  latVec[2] = NCELLSPERSIDE*FCCspacing;
301  for (n = 0; n < DIM; ++n)
302  {
303  coords[a*DIM + n] = latVec[n];
304  }
305  a++;
306  latVec[0] = NCELLSPERSIDE*FCCspacing;
307  latVec[1] = ((double) i)*FCCspacing;
308  latVec[2] = NCELLSPERSIDE*FCCspacing;
309  for (n = 0; n < DIM; ++n)
310  {
311  coords[a*DIM + n] = latVec[n];
312  }
313  a++;
314  latVec[0] = NCELLSPERSIDE*FCCspacing;
315  latVec[1] = NCELLSPERSIDE*FCCspacing;
316  latVec[2] = ((double) i)*FCCspacing;
317  for (n = 0; n < DIM; ++n)
318  {
319  coords[a*DIM + n] = latVec[n];
320  }
321  a++;
322  }
323  /* add in the remaining corner */
324  for (n = 0; n < DIM; ++n)
325  {
326  coords[a*DIM + n] = NCELLSPERSIDE*FCCspacing;
327  }
328  a++;
329 
330  return;
331 }
332 
333 
334 void fcc_cluster_neighborlist(int half, int numberOfParticles, double* coords,
335  double cutoff, NeighList* nl)
336 {
337  /* local variables */
338  int i;
339  int j;
340  int k;
341  int a;
342 
343  double dx[DIM];
344  double r2;
345  double cutoff2;
346 
347  nl->cutoff = cutoff;
348 
349  cutoff2 = cutoff*cutoff;
350 
351  for (i = 0; i < numberOfParticles; ++i)
352  {
353  a = 0;
354  for (j = 0; j < numberOfParticles; ++j)
355  {
356  r2 = 0.0;
357  for (k = 0; k < DIM; ++k)
358  {
359  dx[k] = coords[j*DIM + k] - coords[i*DIM + k];
360  r2 += dx[k]*dx[k];
361  }
362 
363  if (r2 < cutoff2)
364  {
365  if ((half && i < j) || (!half && i != j))
366  {
367  /* part j is a neighbor of part i */
368  (*nl).neighborList[i*NCLUSTERPARTS + a] = j;
369  a++;
370  }
371  }
372  }
373  /* part i has `a' neighbors */
374  (*nl).NNeighbors[i] = a;
375  }
376 
377  return;
378 }
379 
380 int get_cluster_neigh(void const * const dataObject,
381  int const numberOfNeighborLists,
382  double const * const cutoffs,
383  int const neighborListIndex, int const particleNumber,
384  int * const numberOfNeighbors,
385  int const ** const neighborsOfParticle)
386 {
387  /* local variables */
388  int error = true;
389  NeighList* nl = (NeighList*) dataObject;
391 
392  if ((numberOfNeighborLists != 1) || (cutoffs[0] > nl->cutoff)) return error;
393 
394  if (neighborListIndex != 0) return error;
395 
396  /* initialize numNeigh */
397  *numberOfNeighbors = 0;
398 
399  if ((particleNumber >= numberOfParticles) || (particleNumber < 0)) /* invalid id */
400  {
401  MY_WARNING("Invalid part ID in get_cluster_neigh");
402  return true;
403  }
404 
405  /* set the returned number of neighbors for the returned part */
406  *numberOfNeighbors = (*nl).NNeighbors[particleNumber];
407 
408  /* set the location for the returned neighbor list */
409  *neighborsOfParticle = &((*nl).neighborList[(particleNumber)*numberOfParticles]);
410 
411  return false;
412 }
413 
414 char const * const descriptor()
415 {
416  return
417  "KIM_API_Version := 1.6.0\n"
418  "Unit_length := A\n"
419  "Unit_energy := eV\n"
420  "Unit_charge := e\n"
421  "Unit_temperature := K\n"
422  "Unit_time := ps\n"
423  "\n"
424  "PARTICLE_SPECIES:\n"
425  "Ar spec 1\n"
426  "\n"
427  "CONVENTIONS:\n"
428  "ZeroBasedLists flag\n"
429  "Neigh_LocaAccess flag\n"
430  "NEIGH_PURE_F flag\n"
431  "\n"
432  "MODEL_INPUT:\n"
433  "numberOfParticles integer none []\n"
434  "numberOfSpecies integer none []\n"
435  "particleSpecies integer none [numberOfParticles]\n"
436  "coordinates double length [numberOfParticles,3]\n"
437  "get_neigh method none []\n"
438  "neighObject pointer none []\n"
439  "\n"
440  "MODEL_OUTPUT:\n"
441  "destroy method none []\n"
442  "compute method none []\n"
443  "reinit method none []\n"
444  "cutoff double length []\n"
445  "energy double energy []\n"
446  "forces double force [numberOfParticles,3]\n";
447 }
#define NCELLSPERSIDE
#define MY_ERROR(message)
int get_cluster_neigh(void *kimmdl, int *mode, int *request, int *part, int *numnei, int **nei1part, double **Rij)
LengthUnit const m
void fcc_cluster_neighborlist(int half, int numberOfParticles, double *coords, double cutoff, NeighList *nl)
#define MY_WARNING(message)
char const *const descriptor()
ComputeArgumentName const numberOfParticles
void create_FCC_cluster(double FCCspacing, int nCellsPerSide, double *coords)
#define FCCSPACING
#define NCLUSTERPARTS
LogVerbosity const error