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 one = 1;
192  int request;
193 
194  int* nParts;
195  int* particleSpecies;
196  double* cutoff;
197  double cutsq;
198  double epsilon;
199  double C;
200  double Rzero;
201  double shift;
202  double* Rij_list;
203  double* coords;
204  double* energy;
205  double* force;
206  double* particleEnergy;
207  int numOfPartNeigh;
208  double dummy;
209 
210  /* check to see if we have been asked to compute the forces, */
211  /* particleEnergy, and d1Edr */
212  KIM_API_getm_compute(pkim, &ier, 5*3,
213  "energy", &comp_energy, 1,
214  "forces", &comp_force, 1,
215  "particleEnergy", &comp_particleEnergy, 1,
216  "process_dEdr", &comp_process_dEdr, 1,
217  "process_d2Edr2", &comp_process_d2Edr2, 1
218  );
219  if (KIM_STATUS_OK > ier) {
220  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_compute", ier);
221  return ier; }
222 
223  KIM_API_getm_data(
224  pkim, &ier, 7*3,
225  "cutoff", &cutoff, 1,
226  "numberOfParticles", &nParts, 1,
227  "particleSpecies", &particleSpecies,1,
228  "coordinates", &coords, 1,
229  "energy", &energy, comp_energy,
230  "forces", &force, comp_force,
231  "particleEnergy", &particleEnergy, comp_particleEnergy
232  );
233  if (KIM_STATUS_OK > ier) {
234  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_getm_data", ier);
235  return ier; }
236 
237  /* set value of parameters */
238  cutsq = (*cutoff)*(*cutoff);
239  epsilon = EPSILON;
240  C = PARAM_C;
241  Rzero = RZERO;
242  /* set value of parameter shift */
243  dummy = 0.0;
244  /* call calc_phi with r=cutoff and shift=0.0 */
245  calc_phi(&epsilon, &C, &Rzero, &dummy, cutoff, *cutoff, &shift);
246  /* set shift to -shift */
247  shift = -(shift);
248 
249  /* Check to be sure that the species are correct */
250 
251  ier = KIM_STATUS_FAIL; /* assume an error */
252  for (i = 0; i < *nParts; ++i) {
253  if ( SPECCODE != particleSpecies[i]) {
254  KIM_API_report_error(__LINE__, __FILE__,
255  "Unexpected species detected", ier);
256  return ier; } }
257  ier = KIM_STATUS_OK; /* everything is ok */
258 
259  /* initialize potential energies, forces, and virial term */
260  if (comp_particleEnergy) {
261  for (i = 0; i < *nParts; ++i) {
262  particleEnergy[i] = 0.0; } }
263  if (comp_energy) {
264  *energy = 0.0; }
265 
266  if (comp_force) {
267  for (i = 0; i < *nParts; ++i) {
268  for (k = 0; k < DIM; ++k) {
269  force[i*DIM + k] = 0.0; } } }
270 
271  /* Compute energy and forces */
272 
273  /* loop over particles and compute enregy and forces */
274  for (i = 0; i < *nParts; ++i) {
275  request = i;
276  ier = KIM_API_get_neigh(pkim, one, request, &currentPart,
277  &numOfPartNeigh, &neighListOfCurrentPart,
278  &Rij_list);
279  if (KIM_STATUS_OK != ier) {
280  /* some sort of problem, exit */
281  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_neigh", ier);
282  ier = KIM_STATUS_FAIL;
283  return ier; }
284 
285  /* loop over the neighbors of particle i */
286  for (jj = 0; jj < numOfPartNeigh; ++ jj) {
287  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
288 
289  /* compute relative position vector and squared distance */
290  Rsqij = 0.0;
291  for (k = 0; k < DIM; ++k) {
292  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
293 
294  /* compute squared distance */
295  Rsqij += Rij[k]*Rij[k]; }
296 
297  /* compute energy and force */
298  if (Rsqij < cutsq) {
299  /* particles are interacting ? */
300  R = sqrt(Rsqij);
301  if (comp_process_d2Edr2) {
302  /* compute pair potential and its derivatives */
303  calc_phi_d2phi(&epsilon,
304  &C,
305  &Rzero,
306  &shift,
307  cutoff, R, &phi, &dphi, &d2phi);
308 
309  /* compute dEidr */
310  dEidr = 0.5*dphi;
311  d2Eidr = 0.5*d2phi; }
312  else if (comp_force || comp_process_dEdr) {
313  /* compute pair potential and its derivative */
314  calc_phi_dphi(&epsilon,
315  &C,
316  &Rzero,
317  &shift,
318  cutoff, R, &phi, &dphi);
319 
320  /* compute dEidr */
321  dEidr = 0.5*dphi; }
322  else {
323  /* compute just pair potential */
324  calc_phi(&epsilon,
325  &C,
326  &Rzero,
327  &shift,
328  cutoff, R, &phi); }
329 
330  /* contribution to energy */
331  if (comp_particleEnergy) {
332  particleEnergy[i] += 0.5*phi; }
333  if (comp_energy) {
334  *energy += 0.5*phi; }
335 
336  /* contribution to process_dEdr */
337  if (comp_process_dEdr) {
338  ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
339 
340  /* contribution to process_d2Edr2 */
341  if (comp_process_d2Edr2) {
342  R_pairs[0] = R_pairs[1] = R;
343  Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
344  Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
345  Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
346  i_pairs[0] = i_pairs[1] = i;
347  j_pairs[0] = j_pairs[1] = j;
348 
349  ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
350  &pi_pairs, &pj_pairs); }
351 
352  /* contribution to forces */
353  if (comp_force) {
354  for (k = 0; k < DIM; ++k) {
355  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
356  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */ } } }
357  } /* loop on jj */
358  } /* loop on i */
359 
360  /* everything is great */
361  ier = KIM_STATUS_OK;
362 
363  return ier; }
364 
365 /* Initialization function */
366 int model_init(void* km) {
367  /* KIM variables */
368  intptr_t* pkim = *((intptr_t**) km);
369 
370  /* Local variables */
371  double* model_cutoff;
372  int ier;
373 
374  /* store pointer to function in KIM object */
375  ier = KIM_API_set_method(pkim, "compute", 1, (func_ptr) &compute);
376  if (KIM_STATUS_OK > ier) {
377  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_set_data", ier);
378  return ier; }
379 
380  /* store model cutoff in KIM object */
381  model_cutoff = (double*) KIM_API_get_data(pkim, "cutoff", &ier);
382  if (KIM_STATUS_OK > ier) {
383  KIM_API_report_error(__LINE__, __FILE__, "KIM_API_get_data", ier);
384  return ier; }
385  *model_cutoff = CUTOFF;
386 
387  return KIM_STATUS_OK; }
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
#define RZERO
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
int model_init(void *km)
#define SPECCODE
#define CUTOFF
real(c_double), parameter, public model_cutoff
#define EPSILON
static int compute(void *km)
#define DIM
#define PARAM_C