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;
67 };
68 typedef struct buffer buffer;
69 
70 /* Define prototype for Model create */
71 int model_create(KIM_ModelCreate * const modelCreate,
72  KIM_LengthUnit const requestedLengthUnit,
73  KIM_EnergyUnit const requestedEnergyUnit,
74  KIM_ChargeUnit const requestedChargeUnit,
75  KIM_TemperatureUnit const requestedTemperatureUnit,
76  KIM_TimeUnit const requestedTimeUnit);
77 
78 /* Define prototype for other routines */
79 static int compute(void* km);
80 static int compute_arguments_create(
81  KIM_ModelCompute const * const modelCompute,
82  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
83 static int compute_arguments_destroy(
84  KIM_ModelCompute const * const modelCompute,
85  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
86 
87 /* Define prototypes for pair potential calculations */
88 static void calc_phi(double* epsilon,
89  double* C,
90  double* Rzero,
91  double* shift,
92  double* cutoff, double r, double* phi);
93 
94 static void calc_phi_dphi(double* epsilon,
95  double* C,
96  double* Rzero,
97  double* shift,
98  double* cutoff, double r, double* phi, double* dphi);
99 
100 static void calc_phi_d2phi(double* epsilon,
101  double* C,
102  double* Rzero,
103  double* shift,
104  double* cutoff, double r, double* phi, double* dphi,
105  double* d2phi);
106 
107 /* Calculate pair potential phi(r) */
108 static void calc_phi(double* epsilon,
109  double* C,
110  double* Rzero,
111  double* shift,
112  double* cutoff, double r, double* phi) {
113  /* local variables */
114  double ep;
115  double ep2;
116 
117  ep = exp(-(*C)*(r-*Rzero));
118  ep2 = ep*ep;
119 
120  if (r > *cutoff) {
121  /* Argument exceeds cutoff radius */
122  *phi = 0.0; }
123  else {
124  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift; }
125 
126  return; }
127 
128 /* Calculate pair potential phi(r) and its derivative dphi(r) */
129 static void calc_phi_dphi(double* epsilon,
130  double* C,
131  double* Rzero,
132  double* shift,
133  double* cutoff, double r, double* phi, double* dphi) {
134  /* local variables */
135  double ep;
136  double ep2;
137 
138  ep = exp(-(*C)*(r-*Rzero));
139  ep2 = ep*ep;
140 
141  if (r > *cutoff) {
142  /* Argument exceeds cutoff radius */
143  *phi = 0.0;
144  *dphi = 0.0; }
145  else {
146  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
147  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 ); }
148 
149  return; }
150 
151 /* Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r), */
152 /* d2phi(r) */
153 static void calc_phi_d2phi(double* epsilon,
154  double* C,
155  double* Rzero,
156  double* shift,
157  double* cutoff, double r, double* phi, double* dphi,
158  double* d2phi) {
159  /* local variables */
160  double ep;
161  double ep2;
162 
163  ep = exp(-(*C)*(r-*Rzero));
164  ep2 = ep*ep;
165 
166  if (r > *cutoff) {
167  /* Argument exceeds cutoff radius */
168  *phi = 0.0;
169  *dphi = 0.0;
170  *d2phi = 0.0; }
171  else {
172  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
173  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
174  *d2phi = 2.0*(*epsilon)*(*C)*(*C)*(ep - 2.0*ep2); }
175 
176  return; }
177 
179 //static int compute(void* km) {
180 // /* local variables */
181 // intptr_t* pkim = *((intptr_t**) km);
182 // double R;
183 // double R_pairs[2];
184 // double *pR_pairs = &(R_pairs[0]);
185 // double Rsqij;
186 // double phi;
187 // double dphi;
188 // double d2phi;
189 // double dEidr;
190 // double d2Eidr;
191 // double Rij[DIM];
192 // double *pRij = &(Rij[0]);
193 // double Rij_pairs[2][3];
194 // double *pRij_pairs = &(Rij_pairs[0][0]);
195 // int ier;
196 // int i;
197 // int i_pairs[2];
198 // int *pi_pairs = &(i_pairs[0]);
199 // int j;
200 // int j_pairs[2];
201 // int *pj_pairs = &(j_pairs[0]);
202 // int jj;
203 // int k;
204 // int currentPart;
205 // int* neighListOfCurrentPart;
206 // int comp_energy;
207 // int comp_force;
208 // int comp_particleEnergy;
209 // int comp_process_dEdr;
210 // int comp_process_d2Edr2;
211 // int one = 1;
212 // int request;
213 //
214 // int* nParts;
215 // int* particleSpecies;
216 // double* cutoff;
217 // double cutsq;
218 // double epsilon;
219 // double C;
220 // double Rzero;
221 // double shift;
222 // double* Rij_list;
223 // double* coords;
224 // double* energy;
225 // double* force;
226 // double* particleEnergy;
227 // int numOfPartNeigh;
228 // double dummy;
229 //
230 // /* check to see if we have been asked to compute the forces, */
231 // /* particleEnergy, and d1Edr */
232 // KIM_API_getm_compute(pkim, &ier, 5*3,
233 // "energy", &comp_energy, 1,
234 // "forces", &comp_force, 1,
235 // "particleEnergy", &comp_particleEnergy, 1,
236 // "process_dEdr", &comp_process_dEdr, 1,
237 // "process_d2Edr2", &comp_process_d2Edr2, 1
238 // );
239 // if (KIM_STATUS_OK > ier) {
240 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_compute", ier);
241 // return ier; }
242 //
243 // KIM_API_getm_data(
244 // pkim, &ier, 7*3,
245 // "cutoff", &cutoff, 1,
246 // "numberOfParticles", &nParts, 1,
247 // "particleSpecies", &particleSpecies,1,
248 // "coordinates", &coords, 1,
249 // "energy", &energy, comp_energy,
250 // "forces", &force, comp_force,
251 // "particleEnergy", &particleEnergy, comp_particleEnergy
252 // );
253 // if (KIM_STATUS_OK > ier) {
254 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_data", ier);
255 // return ier; }
256 //
257 // /* set value of parameters */
258 // cutsq = (*cutoff)*(*cutoff);
259 // epsilon = EPSILON;
260 // C = PARAM_C;
261 // Rzero = RZERO;
262 // /* set value of parameter shift */
263 // dummy = 0.0;
264 // /* call calc_phi with r=cutoff and shift=0.0 */
265 // calc_phi(&epsilon, &C, &Rzero, &dummy, cutoff, *cutoff, &shift);
266 // /* set shift to -shift */
267 // shift = -(shift);
268 //
269 // /* Check to be sure that the species are correct */
270 // /**/
271 // ier = KIM_STATUS_FAIL; /* assume an error */
272 // for (i = 0; i < *nParts; ++i) {
273 // if ( SPECCODE != particleSpecies[i]) {
274 // KIM_API_report_error(__LINE__, __FILE__,
275 // "Unexpected species detected", ier);
276 // return ier; } }
277 // ier = KIM_STATUS_OK; /* everything is ok */
278 //
279 // /* initialize potential energies, forces, and virial term */
280 // if (comp_particleEnergy) {
281 // for (i = 0; i < *nParts; ++i) {
282 // particleEnergy[i] = 0.0; } }
283 // if (comp_energy) {
284 // *energy = 0.0; }
285 //
286 // if (comp_force) {
287 // for (i = 0; i < *nParts; ++i) {
288 // for (k = 0; k < DIM; ++k) {
289 // force[i*DIM + k] = 0.0; } } }
290 //
291 // /* Compute energy and forces */
292 //
293 // /* loop over particles and compute enregy and forces */
294 // for (i = 0; i < *nParts; ++i) {
295 // request = i;
296 // ier = KIM_API_get_neigh(pkim, one, request, &currentPart,
297 // &numOfPartNeigh, &neighListOfCurrentPart,
298 // &Rij_list);
299 // if (KIM_STATUS_OK != ier) {
300 // /* some sort of problem, exit */
301 // KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
302 // ier = KIM_STATUS_FAIL;
303 // return ier; }
304 //
305 // /* loop over the neighbors of particle i */
306 // for (jj = 0; jj < numOfPartNeigh; ++ jj) {
307 // j = neighListOfCurrentPart[jj]; /* get neighbor ID */
308 //
309 // /* compute relative position vector and squared distance */
310 // Rsqij = 0.0;
311 // for (k = 0; k < DIM; ++k) {
312 // Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
313 //
314 // /* compute squared distance */
315 // Rsqij += Rij[k]*Rij[k]; }
316 //
317 // /* compute energy and force */
318 // if (Rsqij < cutsq) {
319 // /* particles are interacting ? */
320 // R = sqrt(Rsqij);
321 // if (comp_process_d2Edr2) {
322 // /* compute pair potential and its derivatives */
323 // calc_phi_d2phi(&epsilon,
324 // &C,
325 // &Rzero,
326 // &shift,
327 // cutoff, R, &phi, &dphi, &d2phi);
328 //
329 // /* compute dEidr */
330 // dEidr = 0.5*dphi;
331 // d2Eidr = 0.5*d2phi; }
332 // else if (comp_force || comp_process_dEdr) {
333 // /* compute pair potential and its derivative */
334 // calc_phi_dphi(&epsilon,
335 // &C,
336 // &Rzero,
337 // &shift,
338 // cutoff, R, &phi, &dphi);
339 //
340 // /* compute dEidr */
341 // dEidr = 0.5*dphi; }
342 // else {
343 // /* compute just pair potential */
344 // calc_phi(&epsilon,
345 // &C,
346 // &Rzero,
347 // &shift,
348 // cutoff, R, &phi); }
349 //
350 // /* contribution to energy */
351 // if (comp_particleEnergy) {
352 // particleEnergy[i] += 0.5*phi; }
353 // if (comp_energy) {
354 // *energy += 0.5*phi; }
355 //
356 // /* contribution to process_dEdr */
357 // if (comp_process_dEdr) {
358 // ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
359 //
360 // /* contribution to process_d2Edr2 */
361 // if (comp_process_d2Edr2) {
362 // R_pairs[0] = R_pairs[1] = R;
363 // Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
364 // Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
365 // Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
366 // i_pairs[0] = i_pairs[1] = i;
367 // j_pairs[0] = j_pairs[1] = j;
368 //
369 // ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
370 // &pi_pairs, &pj_pairs); }
371 //
372 // /* contribution to forces */
373 // if (comp_force) {
374 // for (k = 0; k < DIM; ++k) {
375 // force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
376 // force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */ } } }
377 // } /* loop on jj */
378 // } /* loop on i */
379 //
380 // /* everything is great */
381 // ier = KIM_STATUS_OK;
382 //
383 // return ier; }
384 
385 /* Create function */
387 int model_create(KIM_ModelCreate * const modelCreate,
388  KIM_LengthUnit const requestedLengthUnit,
389  KIM_EnergyUnit const requestedEnergyUnit,
390  KIM_ChargeUnit const requestedChargeUnit,
391  KIM_TemperatureUnit const requestedTemperatureUnit,
392  KIM_TimeUnit const requestedTimeUnit)
393 {
394  buffer * bufferPointer;
395  int error;
396 
397  /* set units */
398  LOG_INFORMATION("Set model units");
400  modelCreate, /* ignoring requested units */
406 
407  /* register species */
408  LOG_INFORMATION("Setting species code");
409  error = error ||
410  KIM_ModelCreate_SetSpeciesCode(modelCreate,
412 
413  /* register numbering */
414  LOG_INFORMATION("Setting model numbering");
417 
418  /* register function pointers */
419  LOG_INFORMATION("Register model function pointers");
420  error = error ||
422  modelCreate,
424  (func *) 2);
425  error = error ||
427  modelCreate,
430  error = error ||
432  modelCreate,
435  error = error ||
437  modelCreate,
439  (func *) 2);
440  error = error ||
442  modelCreate,
444  (func *) 2);
445 
446  /* allocate buffer */
447  bufferPointer = (buffer *) malloc(sizeof(buffer));
448 
449  /* store model buffer in KIM object */
450  LOG_INFORMATION("Set influence distance and cutoffs");
452  bufferPointer);
453 
454  /* set buffer values */
455  bufferPointer->influenceDistance = CUTOFF;
456  bufferPointer->cutoff = CUTOFF;
457 
458  /* register influence distance */
460  modelCreate,
461  &(bufferPointer->influenceDistance));
462 
463  /* register cutoff */
465  &(bufferPointer->cutoff));
466 
467  if (error)
468  {
469  free(bufferPointer);
470  LOG_ERROR("Unable to successfully initialize model");
471  return TRUE;
472  }
473  else
474  {
475  printf("%s\n",KIM_ModelCreate_String(modelCreate));
476  return FALSE;
477  }
478 }
479 
480 /* compute arguments create routine */
483  KIM_ModelCompute const * const modelCompute,
484  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
485 {
486  int error;
487  /* register arguments */
488  LOG_INFORMATION("Register argument supportStatus");
489  error =
491  modelComputeArgumentsCreate,
493  error = error ||
495  modelComputeArgumentsCreate,
497  error = error ||
499  modelComputeArgumentsCreate,
502 
503  /* register call backs */
504  LOG_INFORMATION("Register call back supportStatus");
505  error = error ||
507  modelComputeArgumentsCreate,
510  error = error ||
512  modelComputeArgumentsCreate,
515 
516  if (error)
517  {
518  LOG_ERROR("Unable to successfully initialize compute arguments");
519  return TRUE;
520  }
521  else
522  {
524  modelComputeArgumentsCreate));
525  return FALSE;
526  }
527 }
528 
529 /* compue arguments destroy routine */
532  KIM_ModelCompute const * const modelCompute,
533  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
534 {
535  /* nothing to be done */
536 
537  return FALSE;
538 }
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
void KIM_ModelCreate_SetNeighborListCutoffsPointer(KIM_ModelCreate *const modelCreate, int const numberOfCutoffs, double const *const cutoffs)
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)
struct KIM_ModelComputeArgumentsDestroy KIM_ModelComputeArgumentsDestroy
struct buffer buffer
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
KIM_ChargeUnit const KIM_CHARGE_UNIT_unused
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialEnergy
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