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_API_C.h"
48 #include "KIM_API_status.h"
49 
50 /******************************************************************************/
51 /* Below are the definitions and values of all Model parameters */
52 /******************************************************************************/
53 #define DIM 3 /* dimensionality of space */
54 #define SPECCODE 1 /* internal species code */
55 #define CUTOFF 8.15 /* Angstroms */
56 #define EPSILON -0.0134783698072604 /* eV */
57 #define PARAM_C 1.545 /* 1/Angstroms */
58 #define RZERO 3.786 /* Angstroms */
59 
60 
61 /* Define prototype for Model init */
62 int model_init(void* km);
63 
64 /* Define prototype for compute routine */
65 static int compute(void* km);
66 
67 /* Define prototypes for pair potential calculations */
68 static void calc_phi(double* epsilon,
69  double* C,
70  double* Rzero,
71  double* shift,
72  double* cutoff, double r, double* phi);
73 
74 static void calc_phi_dphi(double* epsilon,
75  double* C,
76  double* Rzero,
77  double* shift,
78  double* cutoff, double r, double* phi, double* dphi);
79 
80 static void calc_phi_d2phi(double* epsilon,
81  double* C,
82  double* Rzero,
83  double* shift,
84  double* cutoff, double r, double* phi, double* dphi,
85  double* d2phi);
86 
87 /* Calculate pair potential phi(r) */
88 static void calc_phi(double* epsilon,
89  double* C,
90  double* Rzero,
91  double* shift,
92  double* cutoff, double r, double* phi) {
93  /* local variables */
94  double ep;
95  double ep2;
96 
97  ep = exp(-(*C)*(r-*Rzero));
98  ep2 = ep*ep;
99 
100  if (r > *cutoff) {
101  /* Argument exceeds cutoff radius */
102  *phi = 0.0; }
103  else {
104  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift; }
105 
106  return; }
107 
108 /* Calculate pair potential phi(r) and its derivative dphi(r) */
109 static void calc_phi_dphi(double* epsilon,
110  double* C,
111  double* Rzero,
112  double* shift,
113  double* cutoff, double r, double* phi, double* dphi) {
114  /* local variables */
115  double ep;
116  double ep2;
117 
118  ep = exp(-(*C)*(r-*Rzero));
119  ep2 = ep*ep;
120 
121  if (r > *cutoff) {
122  /* Argument exceeds cutoff radius */
123  *phi = 0.0;
124  *dphi = 0.0; }
125  else {
126  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
127  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 ); }
128 
129  return; }
130 
131 /* Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r), */
132 /* d2phi(r) */
133 static void calc_phi_d2phi(double* epsilon,
134  double* C,
135  double* Rzero,
136  double* shift,
137  double* cutoff, double r, double* phi, double* dphi,
138  double* d2phi) {
139  /* local variables */
140  double ep;
141  double ep2;
142 
143  ep = exp(-(*C)*(r-*Rzero));
144  ep2 = ep*ep;
145 
146  if (r > *cutoff) {
147  /* Argument exceeds cutoff radius */
148  *phi = 0.0;
149  *dphi = 0.0;
150  *d2phi = 0.0; }
151  else {
152  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
153  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
154  *d2phi = 2.0*(*epsilon)*(*C)*(*C)*(ep - 2.0*ep2); }
155 
156  return; }
157 
158 /* compute function */
159 static int compute(void* km) {
160  /* local variables */
161  intptr_t* pkim = *((intptr_t**) km);
162  double R;
163  double R_pairs[2];
164  double *pR_pairs = &(R_pairs[0]);
165  double Rsqij;
166  double phi;
167  double dphi;
168  double d2phi;
169  double dEidr;
170  double d2Eidr;
171  double Rij[DIM];
172  double *pRij = &(Rij[0]);
173  double Rij_pairs[2][3];
174  double *pRij_pairs = &(Rij_pairs[0][0]);
175  int ier;
176  int i;
177  int i_pairs[2];
178  int *pi_pairs = &(i_pairs[0]);
179  int j;
180  int j_pairs[2];
181  int *pj_pairs = &(j_pairs[0]);
182  int jj;
183  int k;
184  int currentPart;
185  int* neighListOfCurrentPart;
186  int comp_energy;
187  int comp_force;
188  int comp_particleEnergy;
189  int comp_process_dEdr;
190  int comp_process_d2Edr2;
191  int NBC;
192  int HalfOrFull;
193  int IterOrLoca;
194  int zero = 0;
195  int one = 1;
196  int request;
197 
198  int* nParts;
199  int* particleSpecies;
200  double* cutoff;
201  double cutsq;
202  double epsilon;
203  double C;
204  double Rzero;
205  double shift;
206  double* Rij_list;
207  double* coords;
208  double* energy;
209  double* force;
210  double* particleEnergy;
211  int* numContrib;
212  int numberContrib;
213  int numOfPartNeigh;
214  const char* NBCstr;
215  double dummy;
216 
217  /* Determine neighbor list boundary condition (NBC) */
218  ier = KIM_API_get_NBC_method(pkim, &NBCstr);
219  if (KIM_STATUS_OK > ier) {
220  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_NBC_method", ier);
221  return ier; }
222  if ((!strcmp("NEIGH_PURE_H",NBCstr)) || (!strcmp("NEIGH_PURE_F",NBCstr))) {
223  NBC = 1; }
224  else if (!strcmp("CLUSTER",NBCstr)) {
225  NBC = 3; }
226  else {
227  ier = KIM_STATUS_FAIL;
228  KIM_API_report_error(__LINE__, __FILE__, "Unknown NBC method", ier);
229  return ier; }
230 
231  /* Determine if Half or Full neighbor lists are being used */
232  /*****************************/
233  /* HalfOrFull = 1 -- Half */
234  /* = 2 -- Full */
235  /*****************************/
236  if (KIM_API_is_half_neighbors(pkim, &ier)) {
237  HalfOrFull = 1; }
238  else {
239  HalfOrFull = 2; }
240 
241  /* determine neighbor list handling mode */
242  if (NBC != 3) {
243  /******************************/
244  /* IterOrLoca = 1 -- Iterator */
245  /* = 2 -- Locator */
246  /******************************/
247  IterOrLoca = KIM_API_get_neigh_mode(pkim, &ier);
248  if (KIM_STATUS_OK > ier) {
249  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh_mode", ier);
250  return ier; }
251  if ((IterOrLoca != 1) && (IterOrLoca != 2)) {
252  ier = KIM_STATUS_FAIL;
253  KIM_API_report_error(__LINE__, __FILE__,
254  "Unsupported IterOrLoca mode", ier);
255  return ier; } }
256  else {
257  IterOrLoca = 2; /* for CLUSTER NBC */ }
258 
259  /* check to see if we have been asked to compute the forces, */
260  /* particleEnergy, and d1Edr */
261  KIM_API_getm_compute(pkim, &ier, 5*3,
262  "energy", &comp_energy, 1,
263  "forces", &comp_force, 1,
264  "particleEnergy", &comp_particleEnergy, 1,
265  "process_dEdr", &comp_process_dEdr, 1,
266  "process_d2Edr2", &comp_process_d2Edr2, 1
267  );
268  if (KIM_STATUS_OK > ier) {
269  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_compute", ier);
270  return ier; }
271 
272  KIM_API_getm_data(
273  pkim, &ier, 8*3,
274  "cutoff", &cutoff, 1,
275  "numberOfParticles", &nParts, 1,
276  "particleSpecies", &particleSpecies,1,
277  "coordinates", &coords, 1,
278  "numberContributingParticles", &numContrib, (HalfOrFull==1),
279  "energy", &energy, comp_energy,
280  "forces", &force, comp_force,
281  "particleEnergy", &particleEnergy, comp_particleEnergy
282  );
283  if (KIM_STATUS_OK > ier) {
284  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_data", ier);
285  return ier; }
286 
287  /* set value of parameters */
288  cutsq = (*cutoff)*(*cutoff);
289  epsilon = EPSILON;
290  C = PARAM_C;
291  Rzero = RZERO;
292  /* set value of parameter shift */
293  dummy = 0.0;
294  /* call calc_phi with r=cutoff and shift=0.0 */
295  calc_phi(&epsilon, &C, &Rzero, &dummy, cutoff, *cutoff, &shift);
296  /* set shift to -shift */
297  shift = -(shift);
298 
299  if (HalfOrFull == 1) {
300  if (3 != NBC) {
301  /* non-CLUSTER cases */
302  numberContrib = *numContrib; }
303  else {
304  /* CLUSTER cases */
305  numberContrib = *nParts; } }
306  else {
307  /* provide initialization even if not used */
308  numberContrib = *nParts; }
309 
310  /* Check to be sure that the species are correct */
311 
312  ier = KIM_STATUS_FAIL; /* assume an error */
313  for (i = 0; i < *nParts; ++i) {
314  if ( SPECCODE != particleSpecies[i]) {
315  KIM_API_report_error(__LINE__, __FILE__,
316  "Unexpected species detected", ier);
317  return ier; } }
318  ier = KIM_STATUS_OK; /* everything is ok */
319 
320  /* initialize potential energies, forces, and virial term */
321  if (comp_particleEnergy) {
322  for (i = 0; i < *nParts; ++i) {
323  particleEnergy[i] = 0.0; } }
324  if (comp_energy) {
325  *energy = 0.0; }
326 
327  if (comp_force) {
328  for (i = 0; i < *nParts; ++i) {
329  for (k = 0; k < DIM; ++k) {
330  force[i*DIM + k] = 0.0; } } }
331 
332  /* Initialize neighbor handling for CLUSTER NBC */
333  if (3 == NBC) {
334  /* CLUSTER */
335  neighListOfCurrentPart = (int *) malloc((*nParts)*sizeof(int)); }
336 
337  /* Initialize neighbor handling for Iterator mode */
338 
339  if (1 == IterOrLoca) {
340  ier = KIM_API_get_neigh(pkim, zero, zero, &currentPart, &numOfPartNeigh,
341  &neighListOfCurrentPart, &Rij_list);
342  /* check for successful initialization */
343  if (KIM_STATUS_NEIGH_ITER_INIT_OK != ier) {
344  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
345  ier = KIM_STATUS_FAIL;
346  return ier; } }
347 
348  /* Compute energy and forces */
349 
350  /* loop over particles and compute enregy and forces */
351  i = -1;
352  while( 1 ) {
353  /* Set up neighbor list for next part for all NBC methods */
354  if (1 == IterOrLoca) {
355  /* ITERATOR mode */
356  ier = KIM_API_get_neigh(pkim, zero, one, &currentPart, &numOfPartNeigh,
357  &neighListOfCurrentPart, &Rij_list);
358  if (KIM_STATUS_NEIGH_ITER_PAST_END == ier) {
359  /* the end of the list, terminate loop */
360  break; }
361  if (KIM_STATUS_OK > ier) {
362  /* some sort of problem, exit */
363  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
364  return ier; }
365 
366  i = currentPart; }
367  else {
368  i++;
369  if (*nParts <= i) {
370  /* incremented past end of list, terminate loop */
371  break; }
372 
373  if (3 == NBC) {
374  /* CLUSTER NBC method */
375  numOfPartNeigh = *nParts - (i + 1);
376  for (k = 0; k < numOfPartNeigh; ++k) {
377  neighListOfCurrentPart[k] = i + k + 1; }
378  ier = KIM_STATUS_OK; }
379  else {
380  request = i;
381  ier = KIM_API_get_neigh(pkim, one, request, &currentPart,
382  &numOfPartNeigh, &neighListOfCurrentPart,
383  &Rij_list);
384  if (KIM_STATUS_OK != ier) {
385  /* some sort of problem, exit */
386  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
387  ier = KIM_STATUS_FAIL;
388  return ier; } } }
389 
390  /* loop over the neighbors of particle i */
391  for (jj = 0; jj < numOfPartNeigh; ++ jj) {
392  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
393 
394  /* compute relative position vector and squared distance */
395  Rsqij = 0.0;
396  for (k = 0; k < DIM; ++k) {
397  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
398 
399  /* compute squared distance */
400  Rsqij += Rij[k]*Rij[k]; }
401 
402  /* compute energy and force */
403  if (Rsqij < cutsq) {
404  /* particles are interacting ? */
405  R = sqrt(Rsqij);
406  if (comp_process_d2Edr2) {
407  /* compute pair potential and its derivatives */
408  calc_phi_d2phi(&epsilon,
409  &C,
410  &Rzero,
411  &shift,
412  cutoff, R, &phi, &dphi, &d2phi);
413 
414  /* compute dEidr */
415  if ((1 == HalfOrFull) && (j < numberContrib)) {
416  /* Half mode -- double contribution */
417  dEidr = dphi;
418  d2Eidr = d2phi; }
419  else {
420  /* Full mode -- regular contribution */
421  dEidr = 0.5*dphi;
422  d2Eidr = 0.5*d2phi; } }
423  else if (comp_force || comp_process_dEdr) {
424  /* compute pair potential and its derivative */
425  calc_phi_dphi(&epsilon,
426  &C,
427  &Rzero,
428  &shift,
429  cutoff, R, &phi, &dphi);
430 
431  /* compute dEidr */
432  if ((1 == HalfOrFull) && (j < numberContrib)) {
433  /* Half mode -- double contribution */
434  dEidr = dphi; }
435  else {
436  /* Full mode -- regular contribution */
437  dEidr = 0.5*dphi; } }
438  else {
439  /* compute just pair potential */
440  calc_phi(&epsilon,
441  &C,
442  &Rzero,
443  &shift,
444  cutoff, R, &phi); }
445 
446  /* contribution to energy */
447  if (comp_particleEnergy) {
448  particleEnergy[i] += 0.5*phi;
449  /* if half list add energy for the other particle in the pair */
450  if ((1 == HalfOrFull) && (j < numberContrib)) {
451  particleEnergy[j] += 0.5*phi; } }
452  if (comp_energy) {
453  if ((1 == HalfOrFull) && (j < numberContrib)) {
454  /* Half mode -- add v to total energy */
455  *energy += phi; }
456  else {
457  /* Full mode -- add half v to total energy */
458  *energy += 0.5*phi; } }
459 
460  /* contribution to process_dEdr */
461  if (comp_process_dEdr) {
462  ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
463 
464  /* contribution to process_d2Edr2 */
465  if (comp_process_d2Edr2) {
466  R_pairs[0] = R_pairs[1] = R;
467  Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
468  Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
469  Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
470  i_pairs[0] = i_pairs[1] = i;
471  j_pairs[0] = j_pairs[1] = j;
472 
473  ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
474  &pi_pairs, &pj_pairs); }
475 
476  /* contribution to forces */
477  if (comp_force) {
478  for (k = 0; k < DIM; ++k) {
479  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
480  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */ } } }
481  } /* loop on jj */
482  } /* infinite while loop (terminated by break statements above) */
483 
484  /* Free temporary storage */
485  if (3 == NBC) {
486  free(neighListOfCurrentPart); }
487 
488  /* everything is great */
489  ier = KIM_STATUS_OK;
490 
491  return ier; }
492 
493 /* Initialization function */
494 int model_init(void* km) {
495  /* KIM variables */
496  intptr_t* pkim = *((intptr_t**) km);
497 
498  /* Local variables */
499  double* model_cutoff;
500  int ier;
501 
502  /* store pointer to function in KIM object */
503  ier = KIM_API_set_method(pkim, "compute", 1, (func_ptr) &compute);
504  if (KIM_STATUS_OK > ier) {
505  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_set_data", ier);
506  return ier; }
507 
508  /* store model cutoff in KIM object */
509  model_cutoff = (double*) KIM_API_get_data(pkim, "cutoff", &ier);
510  if (KIM_STATUS_OK > ier) {
511  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_data", ier);
512  return ier; }
513  *model_cutoff = CUTOFF;
514 
515  return KIM_STATUS_OK; }
#define EPSILON
#define PARAM_C
static int compute(void *km)
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
static void calc_phi_d2phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi, double *d2phi)
ChargeUnit const C
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
int model_init(void *km)
#define DIM
#define SPECCODE
#define RZERO
#define CUTOFF
real(c_double), parameter, public model_cutoff