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