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 /* Release: This file is part of the kim-api-v2.0.0-alpha.0 package. */
41 /* */
42 /******************************************************************************/
43 
44 
45 #include <stdlib.h>
46 #include <stdio.h>
47 #include <string.h>
48 #include <math.h>
49 #include "KIM_Numbering.h"
50 #include "KIM_LanguageName.h"
51 #include "KIM_SpeciesName.h"
52 #include "KIM_SupportStatus.h"
53 #include "KIM_ArgumentName.h"
54 #include "KIM_CallbackName.h"
55 #include "KIM_LogVerbosity.h"
56 #include "KIM_UnitSystem.h"
57 #include "KIM_ModelCreate.h"
58 #include "KIM_ModelRefresh.h"
59 #include "KIM_ModelCompute.h"
60 #include "KIM_ModelDestroy.h"
61 
62 #define TRUE 1
63 #define FALSE 0
64 
65 /******************************************************************************/
66 /* Below are the definitions and values of all Model parameters */
67 /******************************************************************************/
68 #define DIM 3 /* dimensionality of space */
69 #define SPECCODE 1 /* internal species code */
70 #define CUTOFF 8.15 /* Angstroms */
71 #define EPSILON -0.0134783698072604 /* eV */
72 #define PARAM_C 1.545 /* 1/Angstroms */
73 #define RZERO 3.786 /* Angstroms */
74 
75 
76 /* Define prototype for Model init */
77 int model_create(KIM_ModelCreate * const modelCreate,
78  KIM_LengthUnit const requestedLengthUnit,
79  KIM_EnergyUnit const requestedEnergyUnit,
80  KIM_ChargeUnit const requestedChargeUnit,
81  KIM_TemperatureUnit const requestedTemperatureUnit,
82  KIM_TimeUnit const requestedTimeUnit);
83 
84 /* Define prototype for compute routine */
85 static int compute(KIM_ModelCompute const * const modelCompute);
86 static int model_refresh(KIM_ModelRefresh * const modelRefresh);
87 static int model_destroy(KIM_ModelDestroy * const modelDestroy);
88 
89 /* Define prototypes for pair potential calculations */
90 static void calc_phi(double* epsilon,
91  double* C,
92  double* Rzero,
93  double* shift,
94  double* cutoff, double r, double* phi);
95 
96 static void calc_phi_dphi(double* epsilon,
97  double* C,
98  double* Rzero,
99  double* shift,
100  double* cutoff, double r, double* phi, double* dphi);
101 
102 static void calc_phi_d2phi(double* epsilon,
103  double* C,
104  double* Rzero,
105  double* shift,
106  double* cutoff, double r, double* phi, double* dphi,
107  double* d2phi);
108 
109 /* Calculate pair potential phi(r) */
110 static void calc_phi(double* epsilon,
111  double* C,
112  double* Rzero,
113  double* shift,
114  double* cutoff, double r, double* phi) {
115  /* local variables */
116  double ep;
117  double ep2;
118 
119  ep = exp(-(*C)*(r-*Rzero));
120  ep2 = ep*ep;
121 
122  if (r > *cutoff) {
123  /* Argument exceeds cutoff radius */
124  *phi = 0.0; }
125  else {
126  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift; }
127 
128  return; }
129 
130 /* Calculate pair potential phi(r) and its derivative dphi(r) */
131 static void calc_phi_dphi(double* epsilon,
132  double* C,
133  double* Rzero,
134  double* shift,
135  double* cutoff, double r, double* phi, double* dphi) {
136  /* local variables */
137  double ep;
138  double ep2;
139 
140  ep = exp(-(*C)*(r-*Rzero));
141  ep2 = ep*ep;
142 
143  if (r > *cutoff) {
144  /* Argument exceeds cutoff radius */
145  *phi = 0.0;
146  *dphi = 0.0; }
147  else {
148  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
149  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 ); }
150 
151  return; }
152 
153 /* Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r), */
154 /* d2phi(r) */
155 static void calc_phi_d2phi(double* epsilon,
156  double* C,
157  double* Rzero,
158  double* shift,
159  double* cutoff, double r, double* phi, double* dphi,
160  double* d2phi) {
161  /* local variables */
162  double ep;
163  double ep2;
164 
165  ep = exp(-(*C)*(r-*Rzero));
166  ep2 = ep*ep;
167 
168  if (r > *cutoff) {
169  /* Argument exceeds cutoff radius */
170  *phi = 0.0;
171  *dphi = 0.0;
172  *d2phi = 0.0; }
173  else {
174  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
175  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
176  *d2phi = 2.0*(*epsilon)*(*C)*(*C)*(ep - 2.0*ep2); }
177 
178  return; }
179 
180 /* compute function */
182 static int compute(KIM_ModelCompute const * const modelCompute)
183 {
184  /* local variables */
185  double R;
186  double R_pairs[2];
187  double *pR_pairs = &(R_pairs[0]);
188  double Rsqij;
189  double phi;
190  double dphi;
191  double d2phi;
192  double dEidr;
193  double d2Eidr;
194  double Rij[DIM];
195  double *pRij = &(Rij[0]);
196  double Rij_pairs[2][3];
197  double const * pRij_pairs = &(Rij_pairs[0][0]);
198  int ier;
199  int i;
200  int i_pairs[2];
201  int *pi_pairs = &(i_pairs[0]);
202  int j;
203  int j_pairs[2];
204  int *pj_pairs = &(j_pairs[0]);
205  int jj;
206  int k;
207  int const * neighListOfCurrentPart;
208  int comp_energy;
209  int comp_force;
210  int comp_particleEnergy;
211  int comp_process_dEdr;
212  int comp_process_d2Edr2;
213 
214  int* nParts;
217  double* cutoff;
218  double cutsq;
219  double epsilon;
220  double C;
221  double Rzero;
222  double shift;
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  LOG_INFORMATION("Checking if call backs are present.");
235  &comp_process_dEdr);
238  &comp_process_d2Edr2);
239 
240  LOG_INFORMATION("Getting data pointers");
241  ier =
243  modelCompute,
245  &nParts)
246  ||
248  modelCompute,
251  ||
253  modelCompute,
256  ||
258  modelCompute,
260  &coords)
261  ||
263  modelCompute,
265  &energy)
266  ||
268  modelCompute,
270  &force)
271  ||
273  modelCompute,
275  &particleEnergy);
276  if (ier) {
277  LOG_ERROR("get data pointers failed");
278  return ier; }
279 
280  comp_energy = (energy != 0);
281  comp_force = (force != 0);
282  comp_particleEnergy = (particleEnergy != 0);
283 
284  /* set value of parameters */
285  KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void**) &cutoff);
286  cutsq = (*cutoff)*(*cutoff);
287  epsilon = EPSILON;
288  C = PARAM_C;
289  Rzero = RZERO;
290  /* set value of parameter shift */
291  dummy = 0.0;
292  /* call calc_phi with r=cutoff and shift=0.0 */
293  calc_phi(&epsilon, &C, &Rzero, &dummy, cutoff, *cutoff, &shift);
294  /* set shift to -shift */
295  shift = -(shift);
296 
297  /* Check to be sure that the species are correct */
298 
299  ier = TRUE; /* assume an error */
300  for (i = 0; i < *nParts; ++i) {
301  if ( SPECCODE != particleSpeciesCodes[i]) {
302  LOG_ERROR("Unexpected species code detected");
303  return ier; } }
304  ier = FALSE; /* everything is ok */
305 
306  /* initialize potential energies, forces, and virial term */
307  LOG_INFORMATION("Initializing data");
308  if (comp_particleEnergy) {
309  for (i = 0; i < *nParts; ++i) {
310  particleEnergy[i] = 0.0; } }
311  if (comp_energy) {
312  *energy = 0.0; }
313 
314  if (comp_force) {
315  for (i = 0; i < *nParts; ++i) {
316  for (k = 0; k < DIM; ++k) {
317  force[i*DIM + k] = 0.0; } } }
318 
319  /* Compute energy and forces */
320 
321  /* loop over particles and compute enregy and forces */
322  LOG_INFORMATION("Starting main compute loop");
323  for (i = 0; i< *nParts; ++i) {
324  if (particleContributing[i])
325  {
327  modelCompute, 0, i, &numOfPartNeigh, &neighListOfCurrentPart);
328  if (ier) {
329  /* some sort of problem, exit */
330  LOG_ERROR("GetNeighborList failed");
331  ier = TRUE;
332  return ier; }
333 
334  /* loop over the neighbors of particle i */
335  for (jj = 0; jj < numOfPartNeigh; ++ jj) {
336  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
337 
338  /* compute relative position vector and squared distance */
339  Rsqij = 0.0;
340  for (k = 0; k < DIM; ++k) {
341  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
342 
343  /* compute squared distance */
344  Rsqij += Rij[k]*Rij[k]; }
345 
346  /* compute energy and force */
347  if (Rsqij < cutsq) {
348  /* particles are interacting ? */
349  R = sqrt(Rsqij);
350  if (comp_process_d2Edr2) {
351  /* compute pair potential and its derivatives */
352  calc_phi_d2phi(&epsilon,
353  &C,
354  &Rzero,
355  &shift,
356  cutoff, R, &phi, &dphi, &d2phi);
357 
358  /* compute dEidr */
359  /* Full mode -- regular contribution */
360  dEidr = 0.5*dphi;
361  d2Eidr = 0.5*d2phi; }
362  else if (comp_force || comp_process_dEdr) {
363  /* compute pair potential and its derivative */
364  calc_phi_dphi(&epsilon,
365  &C,
366  &Rzero,
367  &shift,
368  cutoff, R, &phi, &dphi);
369 
370  /* compute dEidr */
371  /* Full mode -- regular contribution */
372  dEidr = 0.5*dphi; }
373  else {
374  /* compute just pair potential */
375  calc_phi(&epsilon,
376  &C,
377  &Rzero,
378  &shift,
379  cutoff, R, &phi); }
380 
381  /* contribution to energy */
382  if (comp_particleEnergy) {
383  particleEnergy[i] += 0.5*phi; }
384  if (comp_energy) {
385  /* Full mode -- add half v to total energy */
386  *energy += 0.5*phi; }
387 
388  /* contribution to process_dEdr */
389  if (comp_process_dEdr) {
391  modelCompute, dEidr, R, pRij, i, j); }
392 
393  /* contribution to process_d2Edr2 */
394  if (comp_process_d2Edr2) {
395  R_pairs[0] = R_pairs[1] = R;
396  Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
397  Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
398  Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
399  i_pairs[0] = i_pairs[1] = i;
400  j_pairs[0] = j_pairs[1] = j;
401 
403  modelCompute, d2Eidr, pR_pairs, pRij_pairs, pi_pairs, pj_pairs);
404  }
405 
406  /* contribution to forces */
407  if (comp_force) {
408  for (k = 0; k < DIM; ++k) {
409  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
410  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */ } } }
411  } /* loop on jj */
412  } /* if contributing */
413  } /* loop on i */
414  LOG_INFORMATION("Finished compute loop");
415 
416  /* everything is great */
417  ier = FALSE;
418 
419  return ier; }
420 
421 /* Create function */
423 int model_create(KIM_ModelCreate * const modelCreate,
424  KIM_LengthUnit const requestedLengthUnit,
425  KIM_EnergyUnit const requestedEnergyUnit,
426  KIM_ChargeUnit const requestedChargeUnit,
427  KIM_TemperatureUnit const requestedTemperatureUnit,
428  KIM_TimeUnit const requestedTimeUnit)
429 {
430  double* model_cutoff;
431  int error;
432 
433  /* register numbering */
434  LOG_INFORMATION("Setting model numbering");
437 
438  /* register species */
439  LOG_INFORMATION("Setting species code");
440  error = error ||
441  KIM_ModelCreate_SetSpeciesCode(modelCreate,
443 
444  /* register arguments */
445  LOG_INFORMATION("Register argument supportStatus");
446  error = error ||
448  modelCreate,
450  error = error ||
452  modelCreate,
454  error = error ||
456  modelCreate,
458 
459  /* register call backs */
460  LOG_INFORMATION("Register call back supportStatus");
461  error = error ||
463  modelCreate,
465  error = error ||
467  modelCreate,
469 
470  /* register function pointers */
471  LOG_INFORMATION("Register model function pointers");
472  error = error ||
475  (func *) &compute);
476  error = error ||
479  (func *) &model_destroy);
480  error = error ||
483  (func *) &model_refresh);
484 
485  /* set units */
486  LOG_INFORMATION("Set model units");
488  modelCreate, /* ignoring requested units */
494 
495  /* store model cutoff in KIM object */
496  LOG_INFORMATION("Set influence distance and cutoffs");
497  model_cutoff = (double*) malloc(sizeof(double));
498  *model_cutoff = CUTOFF;
500  model_cutoff);
501 
502  /* register influence distance and cutoffs */
504  model_cutoff);
506  model_cutoff);
507 
508  if (error)
509  {
510  LOG_ERROR("Unable to successfully initialize model");
511  return TRUE;
512  }
513  else
514  return FALSE;
515 }
516 
517 /* refresh function */
519 static int model_refresh(KIM_ModelRefresh * const modelRefresh)
520 {
521  /* Local variables */
522  double * model_cutoff;
523 
524  /* get model buffer from KIM object */
525  LOG_INFORMATION("Getting model buffer");
527  (void **) &model_cutoff);
528 
529  LOG_INFORMATION("Resetting influence distance and cutoffs");
531  modelRefresh, model_cutoff);
533  model_cutoff);
534 
535  return FALSE;
536 }
537 
538 /* Initialization function */
540 int model_destroy(KIM_ModelDestroy * const modelDestroy) {
541  double* model_cutoff;
542 
543  LOG_INFORMATION("Getting buffer");
544  KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &model_cutoff);
545  LOG_INFORMATION("Freeing model memory");
546  free(model_cutoff);
547 
548  return FALSE; }
KIM_SupportStatus const KIM_SUPPORT_STATUS_optional
KIM_ArgumentName const KIM_ARGUMENT_NAME_particleContributing
KIM_ArgumentName const KIM_ARGUMENT_NAME_numberOfParticles
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)
void KIM_ModelDestroy_GetModelBufferPointer(KIM_ModelDestroy const *const modelDestroy, void **const ptr)
KIM_LanguageName const KIM_LANGUAGE_NAME_c
ArgumentName const particleContributing
int KIM_ModelCreate_SetComputePointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
struct KIM_ModelRefresh KIM_ModelRefresh
static int model_refresh(KIM_ModelRefresh *const modelRefresh)
#define PARAM_C
void KIM_ModelRefresh_GetModelBufferPointer(KIM_ModelRefresh const *const modelRefresh, void **const ptr)
#define EPSILON
KIM_TemperatureUnit const KIM_TEMPERATURE_UNIT_unused
KIM_SpeciesName const KIM_SPECIES_NAME_Ar
ChargeUnit const C
KIM_TimeUnit const KIM_TIME_UNIT_unused
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
void KIM_ModelCreate_SetNeighborListCutoffsPointer(KIM_ModelCreate *const modelCreate, int const numberOfCutoffs, double const *const cutoffs)
void KIM_ModelRefresh_SetNeighborListCutoffsPointer(KIM_ModelRefresh *const modelRefresh, int const numberOfCutoffs, double const *const cutoffs)
#define TRUE
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialParticleEnergy
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialForces
struct KIM_ModelDestroy KIM_ModelDestroy
#define FALSE
void KIM_ModelCreate_SetModelBufferPointer(KIM_ModelCreate *const modelCreate, void *const ptr)
static void calc_phi_d2phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi, double *d2phi)
int KIM_ModelCreate_SetRefreshPointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
int KIM_ModelCompute_GetNeighborList(KIM_ModelCompute const *const modelCompute, int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle)
void KIM_ModelCompute_GetModelBufferPointer(KIM_ModelCompute const *const modelCompute, void **const ptr)
void() func()
Definition: KIM_func.h:39
ArgumentName const particleSpeciesCodes
int KIM_ModelCreate_SetDestroyPointer(KIM_ModelCreate *const modelCreate, KIM_LanguageName const languageName, func *const fptr)
int KIM_ModelCompute_IsCallbackPresent(KIM_ModelCompute const *const modelCompute, KIM_CallbackName const callbackName, int *const present)
#define RZERO
#define CUTOFF
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_ModelCompute_GetArgumentPointerInteger(KIM_ModelCompute const *const modelCompute, KIM_ArgumentName const argumentName, int **const ptr)
#define LOG_INFORMATION(message)
int KIM_ModelCreate_SetModelNumbering(KIM_ModelCreate *const modelCreate, KIM_Numbering const numbering)
static int compute(KIM_ModelCompute const *const modelCompute)
int KIM_ModelCreate_SetArgumentSupportStatus(KIM_ModelCreate *const modelCreate, KIM_ArgumentName const argumentName, KIM_SupportStatus const supportStatus)
int KIM_ModelCreate_SetSpeciesCode(KIM_ModelCreate *const modelCreate, KIM_SpeciesName const speciesName, int const code)
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
void KIM_ModelRefresh_SetInfluenceDistancePointer(KIM_ModelRefresh *const modelRefresh, double *const influenceDistance)
KIM_CallbackName const KIM_CALLBACK_NAME_ProcessDEDrTerm
KIM_CallbackName const KIM_CALLBACK_NAME_ProcessD2EDr2Term
KIM_ChargeUnit const KIM_CHARGE_UNIT_unused
int KIM_ModelCompute_ProcessDEDrTerm(KIM_ModelCompute const *const modelCompute, double const de, double const r, double const *const dx, int const i, int const j)
KIM_ArgumentName const KIM_ARGUMENT_NAME_coordinates
int KIM_ModelCompute_ProcessD2EDr2Term(KIM_ModelCompute const *const modelCompute, double const de, double const *const r, double const *const dx, int const *const i, int const *const j)
KIM_ArgumentName const KIM_ARGUMENT_NAME_particleSpeciesCodes
int KIM_ModelCompute_GetArgumentPointerDouble(KIM_ModelCompute const *const modelCompute, KIM_ArgumentName const argumentName, double **const ptr)
static int model_destroy(KIM_ModelDestroy *const modelDestroy)
real(c_double), parameter, public model_cutoff
struct KIM_ModelCreate KIM_ModelCreate
int KIM_ModelCreate_SetCallbackSupportStatus(KIM_ModelCreate *const modelCreate, KIM_CallbackName const callbackName, KIM_SupportStatus const supportStatus)
#define DIM
#define LOG_ERROR(message)
KIM_LengthUnit const KIM_LENGTH_UNIT_A
#define SPECCODE
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialEnergy
LogVerbosity const error
KIM_Numbering const KIM_NUMBERING_zeroBased
struct KIM_ModelCompute KIM_ModelCompute