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