KIM API V2
ex_model_Ar_P_Morse_07C.c
Go to the documentation of this file.
1 /* */
2 /* CDDL HEADER START */
3 /* */
4 /* The contents of this file are subject to the terms of the Common */
5 /* Development and Distribution License Version 1.0 (the "License"). */
6 /* */
7 /* You can obtain a copy of the license at */
8 /* http://www.opensource.org/licenses/CDDL-1.0. See the License for the */
9 /* specific language governing permissions and limitations under the License. */
10 /* */
11 /* When distributing Covered Code, include this CDDL HEADER in each file and */
12 /* include the License file in a prominent location with the name */
13 /* LICENSE.CDDL. If applicable, add the following below this CDDL HEADER, */
14 /* with the fields enclosed by brackets "[]" replaced with your own */
15 /* identifying information: */
16 /* */
17 /* Portions Copyright (c) [yyyy] [name of copyright owner]. */
18 /* All rights reserved. */
19 /* */
20 /* CDDL HEADER END */
21 /* */
22 
23 /* */
24 /* Copyright (c) 2013--2018, Regents of the University of Minnesota. */
25 /* All rights reserved. */
26 /* */
27 /* Contributors: */
28 /* Ryan S. Elliott */
29 /* Ellad B. Tadmor */
30 /* Stephen M. Whalen */
31 /* */
32 
33 /******************************************************************************/
34 /* */
35 /* ex_model_Ar_P_Morse_07C pair potential KIM Model */
36 /* shifted to have zero energy at the cutoff radius */
37 /* */
38 /* Language: C */
39 /* */
40 /******************************************************************************/
41 
42 
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <string.h>
46 #include <math.h>
47 #include "KIM_ModelHeaders.h"
48 
49 #define TRUE 1
50 #define FALSE 0
51 
52 /******************************************************************************/
53 /* Below are the definitions and values of all Model parameters */
54 /******************************************************************************/
55 #define DIM 3 /* dimensionality of space */
56 #define SPECCODE 1 /* internal species code */
57 #define CUTOFF 8.15 /* Angstroms */
58 #define EPSILON -0.0134783698072604 /* eV */
59 #define PARAM_C 1.545 /* 1/Angstroms */
60 #define RZERO 3.786 /* Angstroms */
61 
62 /* Model buffer definition */
63 struct buffer
64 {
65  double influenceDistance;
66  double cutoff;
68  int halfListHint;
69 };
70 typedef struct buffer buffer;
71 
72 /* Define prototype for Model create */
73 int model_create(KIM_ModelCreate * const modelCreate,
74  KIM_LengthUnit const requestedLengthUnit,
75  KIM_EnergyUnit const requestedEnergyUnit,
76  KIM_ChargeUnit const requestedChargeUnit,
77  KIM_TemperatureUnit const requestedTemperatureUnit,
78  KIM_TimeUnit const requestedTimeUnit);
79 
80 /* Define prototype for other routines */
81 static int compute(void* km);
82 static int compute_arguments_create(
83  KIM_ModelCompute const * const modelCompute,
84  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
85 static int compute_arguments_destroy(
86  KIM_ModelCompute const * const modelCompute,
87  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
88 
89 /* Define prototypes for pair potential calculations */
90 static void calc_phi(double* epsilon,
91  double* C,
92  double* Rzero,
93  double* shift,
94  double* cutoff, double r, double* phi);
95 
96 static void calc_phi_dphi(double* epsilon,
97  double* C,
98  double* Rzero,
99  double* shift,
100  double* cutoff, double r, double* phi, double* dphi);
101 
102 static void calc_phi_d2phi(double* epsilon,
103  double* C,
104  double* Rzero,
105  double* shift,
106  double* cutoff, double r, double* phi, double* dphi,
107  double* d2phi);
108 
109 /* Calculate pair potential phi(r) */
110 static void calc_phi(double* epsilon,
111  double* C,
112  double* Rzero,
113  double* shift,
114  double* cutoff, double r, double* phi) {
115  /* local variables */
116  double ep;
117  double ep2;
118 
119  ep = exp(-(*C)*(r-*Rzero));
120  ep2 = ep*ep;
121 
122  if (r > *cutoff) {
123  /* Argument exceeds cutoff radius */
124  *phi = 0.0; }
125  else {
126  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift; }
127 
128  return; }
129 
130 /* Calculate pair potential phi(r) and its derivative dphi(r) */
131 static void calc_phi_dphi(double* epsilon,
132  double* C,
133  double* Rzero,
134  double* shift,
135  double* cutoff, double r, double* phi, double* dphi) {
136  /* local variables */
137  double ep;
138  double ep2;
139 
140  ep = exp(-(*C)*(r-*Rzero));
141  ep2 = ep*ep;
142 
143  if (r > *cutoff) {
144  /* Argument exceeds cutoff radius */
145  *phi = 0.0;
146  *dphi = 0.0; }
147  else {
148  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
149  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 ); }
150 
151  return; }
152 
153 /* Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r), */
154 /* d2phi(r) */
155 static void calc_phi_d2phi(double* epsilon,
156  double* C,
157  double* Rzero,
158  double* shift,
159  double* cutoff, double r, double* phi, double* dphi,
160  double* d2phi) {
161  /* local variables */
162  double ep;
163  double ep2;
164 
165  ep = exp(-(*C)*(r-*Rzero));
166  ep2 = ep*ep;
167 
168  if (r > *cutoff) {
169  /* Argument exceeds cutoff radius */
170  *phi = 0.0;
171  *dphi = 0.0;
172  *d2phi = 0.0; }
173  else {
174  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
175  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
176  *d2phi = 2.0*(*epsilon)*(*C)*(*C)*(ep - 2.0*ep2); }
177 
178  return; }
179 
181 //static int compute(void* km) {
182 // /* local variables */
183 // intptr_t* pkim = *((intptr_t**) km);
184 // double R;
185 // double R_pairs[2];
186 // double *pR_pairs = &(R_pairs[0]);
187 // double Rsqij;
188 // double phi;
189 // double dphi;
190 // double d2phi;
191 // double dEidr;
192 // double d2Eidr;
193 // double Rij[DIM];
194 // double *pRij = &(Rij[0]);
195 // double Rij_pairs[2][3];
196 // double *pRij_pairs = &(Rij_pairs[0][0]);
197 // int ier;
198 // int i;
199 // int i_pairs[2];
200 // int *pi_pairs = &(i_pairs[0]);
201 // int j;
202 // int j_pairs[2];
203 // int *pj_pairs = &(j_pairs[0]);
204 // int jj;
205 // int k;
206 // int currentPart;
207 // int* neighListOfCurrentPart;
208 // int comp_energy;
209 // int comp_force;
210 // int comp_particleEnergy;
211 // int comp_process_dEdr;
212 // int comp_process_d2Edr2;
213 // int one = 1;
214 // int request;
215 //
216 // int* nParts;
217 // int* particleSpecies;
218 // double* cutoff;
219 // double cutsq;
220 // double epsilon;
221 // double C;
222 // double Rzero;
223 // double shift;
224 // double* Rij_list;
225 // double* coords;
226 // double* energy;
227 // double* force;
228 // double* particleEnergy;
229 // int numOfPartNeigh;
230 // double dummy;
231 //
232 // /* check to see if we have been asked to compute the forces, */
233 // /* particleEnergy, and d1Edr */
234 // KIM_API_getm_compute(pkim, &ier, 5*3,
235 // "energy", &comp_energy, 1,
236 // "forces", &comp_force, 1,
237 // "particleEnergy", &comp_particleEnergy, 1,
238 // "process_dEdr", &comp_process_dEdr, 1,
239 // "process_d2Edr2", &comp_process_d2Edr2, 1
240 // );
241 // if (KIM_STATUS_OK > ier) {
242 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_compute", ier);
243 // return ier; }
244 //
245 // KIM_API_getm_data(
246 // pkim, &ier, 7*3,
247 // "cutoff", &cutoff, 1,
248 // "numberOfParticles", &nParts, 1,
249 // "particleSpecies", &particleSpecies,1,
250 // "coordinates", &coords, 1,
251 // "energy", &energy, comp_energy,
252 // "forces", &force, comp_force,
253 // "particleEnergy", &particleEnergy, comp_particleEnergy
254 // );
255 // if (KIM_STATUS_OK > ier) {
256 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_data", ier);
257 // return ier; }
258 //
259 // /* set value of parameters */
260 // cutsq = (*cutoff)*(*cutoff);
261 // epsilon = EPSILON;
262 // C = PARAM_C;
263 // Rzero = RZERO;
264 // /* set value of parameter shift */
265 // dummy = 0.0;
266 // /* call calc_phi with r=cutoff and shift=0.0 */
267 // calc_phi(&epsilon, &C, &Rzero, &dummy, cutoff, *cutoff, &shift);
268 // /* set shift to -shift */
269 // shift = -(shift);
270 //
271 // /* Check to be sure that the species are correct */
272 // /**/
273 // ier = KIM_STATUS_FAIL; /* assume an error */
274 // for (i = 0; i < *nParts; ++i) {
275 // if ( SPECCODE != particleSpecies[i]) {
276 // KIM_API_report_error(__LINE__, __FILE__,
277 // "Unexpected species detected", ier);
278 // return ier; } }
279 // ier = KIM_STATUS_OK; /* everything is ok */
280 //
281 // /* initialize potential energies, forces, and virial term */
282 // if (comp_particleEnergy) {
283 // for (i = 0; i < *nParts; ++i) {
284 // particleEnergy[i] = 0.0; } }
285 // if (comp_energy) {
286 // *energy = 0.0; }
287 //
288 // if (comp_force) {
289 // for (i = 0; i < *nParts; ++i) {
290 // for (k = 0; k < DIM; ++k) {
291 // force[i*DIM + k] = 0.0; } } }
292 //
293 // /* Compute energy and forces */
294 //
295 // /* loop over particles and compute enregy and forces */
296 // for (i = 0; i < *nParts; ++i) {
297 // request = i;
298 // ier = KIM_API_get_neigh(pkim, one, request, &currentPart,
299 // &numOfPartNeigh, &neighListOfCurrentPart,
300 // &Rij_list);
301 // if (KIM_STATUS_OK != ier) {
302 // /* some sort of problem, exit */
303 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
304 // ier = KIM_STATUS_FAIL;
305 // return ier; }
306 //
307 // /* loop over the neighbors of particle i */
308 // for (jj = 0; jj < numOfPartNeigh; ++ jj) {
309 // j = neighListOfCurrentPart[jj]; /* get neighbor ID */
310 //
311 // /* compute relative position vector and squared distance */
312 // Rsqij = 0.0;
313 // for (k = 0; k < DIM; ++k) {
314 // Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
315 //
316 // /* compute squared distance */
317 // Rsqij += Rij[k]*Rij[k]; }
318 //
319 // /* compute energy and force */
320 // if (Rsqij < cutsq) {
321 // /* particles are interacting ? */
322 // R = sqrt(Rsqij);
323 // if (comp_process_d2Edr2) {
324 // /* compute pair potential and its derivatives */
325 // calc_phi_d2phi(&epsilon,
326 // &C,
327 // &Rzero,
328 // &shift,
329 // cutoff, R, &phi, &dphi, &d2phi);
330 //
331 // /* compute dEidr */
332 // dEidr = 0.5*dphi;
333 // d2Eidr = 0.5*d2phi; }
334 // else if (comp_force || comp_process_dEdr) {
335 // /* compute pair potential and its derivative */
336 // calc_phi_dphi(&epsilon,
337 // &C,
338 // &Rzero,
339 // &shift,
340 // cutoff, R, &phi, &dphi);
341 //
342 // /* compute dEidr */
343 // dEidr = 0.5*dphi; }
344 // else {
345 // /* compute just pair potential */
346 // calc_phi(&epsilon,
347 // &C,
348 // &Rzero,
349 // &shift,
350 // cutoff, R, &phi); }
351 //
352 // /* contribution to energy */
353 // if (comp_particleEnergy) {
354 // particleEnergy[i] += 0.5*phi; }
355 // if (comp_energy) {
356 // *energy += 0.5*phi; }
357 //
358 // /* contribution to process_dEdr */
359 // if (comp_process_dEdr) {
360 // ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
361 //
362 // /* contribution to process_d2Edr2 */
363 // if (comp_process_d2Edr2) {
364 // R_pairs[0] = R_pairs[1] = R;
365 // Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
366 // Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
367 // Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
368 // i_pairs[0] = i_pairs[1] = i;
369 // j_pairs[0] = j_pairs[1] = j;
370 //
371 // ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
372 // &pi_pairs, &pj_pairs); }
373 //
374 // /* contribution to forces */
375 // if (comp_force) {
376 // for (k = 0; k < DIM; ++k) {
377 // force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
378 // force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */ } } }
379 // } /* loop on jj */
380 // } /* loop on i */
381 //
382 // /* everything is great */
383 // ier = KIM_STATUS_OK;
384 //
385 // return ier; }
386 
387 /* Create function */
389 int model_create(KIM_ModelCreate * const modelCreate,
390  KIM_LengthUnit const requestedLengthUnit,
391  KIM_EnergyUnit const requestedEnergyUnit,
392  KIM_ChargeUnit const requestedChargeUnit,
393  KIM_TemperatureUnit const requestedTemperatureUnit,
394  KIM_TimeUnit const requestedTimeUnit)
395 {
396  buffer * bufferPointer;
397  int error;
398 
399  /* set units */
400  LOG_INFORMATION("Set model units");
402  modelCreate, /* ignoring requested units */
408 
409  /* register species */
410  LOG_INFORMATION("Setting species code");
411  error = error ||
412  KIM_ModelCreate_SetSpeciesCode(modelCreate,
414 
415  /* register numbering */
416  LOG_INFORMATION("Setting model numbering");
419 
420  /* register function pointers */
421  LOG_INFORMATION("Register model function pointers");
422  error = error ||
424  modelCreate,
426  (func *) 2);
427  error = error ||
429  modelCreate,
432  error = error ||
434  modelCreate,
437  error = error ||
439  modelCreate,
441  (func *) 2);
442  error = error ||
444  modelCreate,
446  (func *) 2);
447 
448  /* allocate buffer */
449  bufferPointer = (buffer *) malloc(sizeof(buffer));
450 
451  /* store model buffer in KIM object */
452  LOG_INFORMATION("Set influence distance and cutoffs");
454  bufferPointer);
455 
456  /* set buffer values */
457  bufferPointer->influenceDistance = CUTOFF;
458  bufferPointer->cutoff = CUTOFF;
459  bufferPointer->paddingNeighborHint = 1;
460  bufferPointer->halfListHint = 0;
461 
462  /* register influence distance */
464  modelCreate,
465  &(bufferPointer->influenceDistance));
466 
467  /* register cutoff */
469  modelCreate,
470  1,
471  &(bufferPointer->cutoff),
472  &(bufferPointer->paddingNeighborHint),
473  &(bufferPointer->halfListHint));
474 
475  if (error)
476  {
477  free(bufferPointer);
478  LOG_ERROR("Unable to successfully initialize model");
479  return TRUE;
480  }
481  else
482  {
483  printf("%s\n",KIM_ModelCreate_String(modelCreate));
484  return FALSE;
485  }
486 }
487 
488 /* compute arguments create routine */
491  KIM_ModelCompute const * const modelCompute,
492  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
493 {
494  int error;
495  /* register arguments */
496  LOG_INFORMATION("Register argument supportStatus");
497  error =
499  modelComputeArgumentsCreate,
501  error = error ||
503  modelComputeArgumentsCreate,
505  error = error ||
507  modelComputeArgumentsCreate,
510 
511  /* register call backs */
512  LOG_INFORMATION("Register call back supportStatus");
513  error = error ||
515  modelComputeArgumentsCreate,
518  error = error ||
520  modelComputeArgumentsCreate,
523 
524  if (error)
525  {
526  LOG_ERROR("Unable to successfully initialize compute arguments");
527  return TRUE;
528  }
529  else
530  {
532  modelComputeArgumentsCreate));
533  exit(0);
534  return FALSE;
535  }
536 }
537 
538 /* compue arguments destroy routine */
541  KIM_ModelCompute const * const modelCompute,
542  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
543 {
544  /* nothing to be done */
545 
546  return FALSE;
547 }
KIM_SupportStatus const KIM_SUPPORT_STATUS_optional
KIM_ComputeCallbackName const KIM_COMPUTE_CALLBACK_NAME_ProcessD2EDr2Term
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy
#define TRUE
KIM_LanguageName const KIM_LANGUAGE_NAME_c
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialForces
int KIM_ModelCreate_SetComputePointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
#define CUTOFF
int KIM_ModelComputeArgumentsCreate_SetCallbackSupportStatus(KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate, KIM_ComputeCallbackName const computeCallbackName, KIM_SupportStatus const supportStatus)
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
KIM_TemperatureUnit const KIM_TEMPERATURE_UNIT_unused
#define SPECCODE
static void calc_phi_d2phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi, double *d2phi)
KIM_SpeciesName const KIM_SPECIES_NAME_Ar
ChargeUnit const C
static int compute_arguments_create(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
KIM_TimeUnit const KIM_TIME_UNIT_unused
int KIM_ModelCreate_SetComputeArgumentsCreatePointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
char const *const KIM_ModelCreate_String(KIM_ModelCreate const *const modelCreate)
#define FALSE
void KIM_ModelCreate_SetModelBufferPointer(KIM_ModelCreate *const modelCreate, void *const ptr)
int KIM_ModelCreate_SetRefreshPointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
int KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate, KIM_ComputeArgumentName const computeArgumentName, KIM_SupportStatus const supportStatus)
void() func()
Definition: KIM_func.h:39
int KIM_ModelCreate_SetComputeArgumentsDestroyPointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
#define LOG_ERROR(message)
int KIM_ModelCreate_SetDestroyPointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
static int compute(void *km)
struct KIM_ModelComputeArgumentsCreate KIM_ModelComputeArgumentsCreate
void KIM_ModelCreate_SetInfluenceDistancePointer(KIM_ModelCreate *const modelCreate, double *const influenceDistance)
int KIM_ModelCreate_SetUnits(KIM_ModelCreate *const modelCreate, KIM_LengthUnit const lengthUnit, KIM_EnergyUnit const energyUnit, KIM_ChargeUnit const chargeUnit, KIM_TemperatureUnit const temperatureUnit, KIM_TimeUnit const timeUnit)
int KIM_ModelCreate_SetModelNumbering(KIM_ModelCreate *const modelCreate, KIM_Numbering const numbering)
#define LOG_INFORMATION(message)
int KIM_ModelCreate_SetSpeciesCode(KIM_ModelCreate *const modelCreate, KIM_SpeciesName const speciesName, int const code)
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
double influenceDistance
struct KIM_ModelComputeArgumentsDestroy KIM_ModelComputeArgumentsDestroy
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
KIM_ChargeUnit const KIM_CHARGE_UNIT_unused
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialEnergy
void KIM_ModelCreate_SetNeighborListPointers(KIM_ModelCreate *const modelCreate, int const numberOfNeighborLists, double const *const cutoffs, int const *const paddingNeighborHints, int const *const halfListHints)
int model_create(KIM_ModelCreate *const modelCreate, KIM_LengthUnit const requestedLengthUnit, KIM_EnergyUnit const requestedEnergyUnit, KIM_ChargeUnit const requestedChargeUnit, KIM_TemperatureUnit const requestedTemperatureUnit, KIM_TimeUnit const requestedTimeUnit)
KIM_ComputeCallbackName const KIM_COMPUTE_CALLBACK_NAME_ProcessDEDrTerm
struct KIM_ModelCreate KIM_ModelCreate
char const *const KIM_ModelComputeArgumentsCreate_String(KIM_ModelComputeArgumentsCreate const *const modelComputeArgumentsCreate)
static int compute_arguments_destroy(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
KIM_LengthUnit const KIM_LENGTH_UNIT_A
LogVerbosity const error
KIM_Numbering const KIM_NUMBERING_zeroBased
struct KIM_ModelCompute KIM_ModelCompute