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