KIM API V2
ex_model_driver_P_Morse.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_driver_P_Morse pair potential KIM Model Driver */
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_ModelDriverHeaders.h"
48 
49 #define TRUE 1
50 #define FALSE 0
51 
52 /******************************************************************************/
53 /* Below are the definitions for some constants */
54 /******************************************************************************/
55 #define DIM 3 /* dimensionality of space */
56 #define SPECCODE 1 /* internal species code */
57 
58 
59 /* Define prototype for Model Driver init */
61  KIM_ModelDriverCreate * const modelDriverCreate,
62  KIM_LengthUnit const requestedLengthUnit,
63  KIM_EnergyUnit const requestedEnergyUnit,
64  KIM_ChargeUnit const requestedChargeUnit,
65  KIM_TemperatureUnit const requestedTemperatureUnit,
66  KIM_TimeUnit const requestedTimeUnit);
67 
68 /* Define prototypes for destroy */
69 /* defined as static to avoid namespace clashes with other codes */
70 static int destroy(KIM_ModelDestroy * const modelDestroy);
71 
72 /* Define prototype for compute routine */
73 static int compute(
74  KIM_ModelCompute const * const modelCompute,
75  KIM_ModelComputeArguments const * const modelComputeArguments);
76 static int compute_arguments_create(
77  KIM_ModelCompute const * const modelCompute,
78  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
79 static int compute_arguments_destroy(
80  KIM_ModelCompute const * const modelCompute,
81  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
82 
83 /* Define prototype for refresh routine */
84 static int refresh(KIM_ModelRefresh * const modelRefresh);
85 
86 /* Define prototypes for pair potential calculations */
87 static void calc_phi(double const * epsilon,
88  double const * C,
89  double const * Rzero,
90  double const * shift,
91  double const cutoff, double const r, double* phi);
92 
93 static void calc_phi_dphi(double const* epsilon,
94  double const* C,
95  double const* Rzero,
96  double const* shift,
97  double const cutoff, double const r,
98  double* phi, double* dphi);
99 
100 /* Define model_buffer structure */
102 {
104  double cutoff;
105  double cutsq;
108  double epsilon;
109  double C;
110  double Rzero;
111  double shift;
112 };
113 
114 
115 /* Calculate pair potential phi(r) */
116 static void calc_phi(double const* epsilon,
117  double const* C,
118  double const* Rzero,
119  double const* shift,
120  double const cutoff, double const r, double* phi)
121 {
122  /* local variables */
123  double ep;
124  double ep2;
125 
126  ep = exp(-(*C)*(r-*Rzero));
127  ep2 = ep*ep;
128 
129  if (r > cutoff)
130  {
131  /* Argument exceeds cutoff radius */
132  *phi = 0.0;
133  }
134  else
135  {
136  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
137  }
138 
139  return;
140 }
141 
142 /* Calculate pair potential phi(r) and its derivative dphi(r) */
143 static void calc_phi_dphi(double const* epsilon,
144  double const* C,
145  double const* Rzero,
146  double const* shift,
147  double const cutoff, double const r,
148  double* phi, double* dphi)
149 {
150  /* local variables */
151  double ep;
152  double ep2;
153 
154  ep = exp(-(*C)*(r-*Rzero));
155  ep2 = ep*ep;
156 
157  if (r > cutoff)
158  {
159  /* Argument exceeds cutoff radius */
160  *phi = 0.0;
161  *dphi = 0.0;
162  }
163  else
164  {
165  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
166  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
167  }
168 
169  return;
170 }
171 
172 /* compute function */
174 static int compute(
175  KIM_ModelCompute const * const modelCompute,
176  KIM_ModelComputeArguments const * const modelComputeArguments)
177 {
178  /* local variables */
179  double R;
180  double Rsqij;
181  double phi;
182  double dphi;
183  double dEidr;
184  double Rij[DIM];
185  int ier;
186  int i;
187  int j;
188  int jj;
189  int k;
190  int const * neighListOfCurrentPart;
191  struct model_buffer* buffer;
192  int comp_energy;
193  int comp_force;
194  int comp_particleEnergy;
195 
196  int* nParts;
199  double cutoff;
200  double* cutsq;
201  double* epsilon;
202  double* C;
203  double* Rzero;
204  double* shift;
205  double* coords;
206  double* energy;
207  double* force;
208  double* particleEnergy;
209  int numOfPartNeigh;
210 
211  /* get buffer from KIM object */
212  KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
213 
214  /* unpack info from the buffer */
216  cutsq = &(buffer->cutsq);
217  epsilon = &(buffer->epsilon);
218  C = &(buffer->C);
219  Rzero = &(buffer->Rzero);
220  shift = &(buffer->shift);
221 
222  ier =
224  modelComputeArguments,
226  &nParts)
227  ||
229  modelComputeArguments,
232  ||
234  modelComputeArguments,
237  ||
239  modelComputeArguments,
241  &coords)
242  ||
244  modelComputeArguments,
246  &energy)
247  ||
249  modelComputeArguments,
251  &force)
252  ||
254  modelComputeArguments,
256  &particleEnergy);
257  if (ier)
258  {
259  LOG_ERROR("GetArgumentPointer");
260  return ier;
261  }
262 
263  comp_energy = (energy != NULL);
264  comp_force = (force != NULL);
265  comp_particleEnergy = (particleEnergy != NULL);
266 
267  /* Check to be sure that the species are correct */
268 
269  ier = TRUE; /* assume an error */
270  for (i = 0; i < *nParts; ++i)
271  {
272  if ( SPECCODE != particleSpeciesCodes[i])
273  {
274  LOG_ERROR("Unexpected species code detected");
275  return ier;
276  }
277  }
278  ier = FALSE; /* everything is ok */
279 
280  /* initialize potential energies, forces, and virial term */
281  if (comp_particleEnergy)
282  {
283  for (i = 0; i < *nParts; ++i)
284  {
285  particleEnergy[i] = 0.0;
286  }
287  }
288  if (comp_energy)
289  {
290  *energy = 0.0;
291  }
292 
293  if (comp_force)
294  {
295  for (i = 0; i < *nParts; ++i)
296  {
297  for (k = 0; k < DIM; ++k)
298  {
299  force[i*DIM + k] = 0.0;
300  }
301  }
302  }
303 
304  /* Compute energy and forces */
305 
306  /* loop over particles and compute enregy and forces */
307  for (i = 0; i < *nParts; ++i)
308  {
309  if (particleContributing[i])
310  {
312  modelComputeArguments,
313  0, i, &numOfPartNeigh, &neighListOfCurrentPart);
314  if (ier)
315  {
316  /* some sort of problem, exit */
317  LOG_ERROR("KIM_get_neigh");
318  ier = TRUE;
319  return ier;
320  }
321 
322  /* loop over the neighbors of particle i */
323  for (jj = 0; jj < numOfPartNeigh; ++ jj)
324  {
325  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
326 
327  /* compute relative position vector and squared distance */
328  Rsqij = 0.0;
329  for (k = 0; k < DIM; ++k)
330  {
331  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
332  /* compute squared distance */
333  Rsqij += Rij[k]*Rij[k];
334  }
335 
336  /* compute energy and force */
337  if (Rsqij < *cutsq)
338  {
339  /* particles are interacting ? */
340  R = sqrt(Rsqij);
341  if (comp_force)
342  {
343  /* compute pair potential and its derivative */
345  C,
346  Rzero,
347  shift,
348  cutoff, R, &phi, &dphi);
349 
350  /* compute dEidr */
351  dEidr = 0.5*dphi;
352  }
353  else
354  {
355  /* compute just pair potential */
357  C,
358  Rzero,
359  shift,
360  cutoff, R, &phi);
361  }
362 
363  /* contribution to energy */
364  if (comp_particleEnergy)
365  {
366  particleEnergy[i] += 0.5*phi;
367  }
368  if (comp_energy)
369  {
370  *energy += 0.5*phi;
371  }
372 
373  /* contribution to forces */
374  if (comp_force)
375  {
376  for (k = 0; k < DIM; ++k)
377  {
378  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
379  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */
380  }
381  }
382  } /* if Rsqij */
383  } /* loop on jj */
384  } /* if particleContributing */
385  } /* infinite while loop (terminated by break statements above) */
386 
387  /* everything is great */
388  ier = FALSE;
389 
390  return ier;
391 }
392 
393 /* Create function */
396  KIM_ModelDriverCreate * const modelDriverCreate,
397  KIM_LengthUnit const requestedLengthUnit,
398  KIM_EnergyUnit const requestedEnergyUnit,
399  KIM_ChargeUnit const requestedChargeUnit,
400  KIM_TemperatureUnit const requestedTemperatureUnit,
401  KIM_TimeUnit const requestedTimeUnit)
402 {
403  /* KIM variables */
404  int numberOfParameterFiles;
405  char const * paramfile1name;
406 
407  /* Local variables */
408  FILE* fid;
409  char speciesNameString[100];
410  KIM_SpeciesName speciesName;
411  double cutoff;
412  double epsilon;
413  double C;
414  double Rzero;
415  int ier;
416  double dummy;
417  struct model_buffer* buffer;
418 
419 
420  /* using fixed units */
421  ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
427  if (ier == TRUE)
428  {
429  LOG_ERROR("Problem setting units");
430  return ier;
431  }
432 
433  ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate,
435  if (ier == TRUE)
436  {
437  LOG_ERROR("Unable to set numbering");
438  return ier;
439  }
440 
441  /* store pointer to functions in KIM object */
442  KIM_ModelDriverCreate_SetDestroyPointer(modelDriverCreate,
444  (func *) destroy);
451  KIM_ModelDriverCreate_SetComputePointer(modelDriverCreate,
453  (func *) compute);
455  modelDriverCreate, KIM_LANGUAGE_NAME_c, (func *) refresh);
456 
457  /* get number of parameter files */
459  modelDriverCreate, &numberOfParameterFiles);
460  /* set paramfile1name */
461  if (numberOfParameterFiles != 1)
462  {
463  ier = TRUE;
464  LOG_ERROR("Incorrect number of parameter files.");
465  return ier;
466  }
468  modelDriverCreate,
469  0,
470  &paramfile1name);
471  if (ier == TRUE)
472  {
473  LOG_ERROR("Unable to get parameter file name.");
474  return ier;
475  }
476 
477  /* Read in model parameters from parameter file */
478  fid = fopen(paramfile1name, "r");
479  if (fid == NULL)
480  {
481  ier = TRUE;
482  LOG_ERROR("Unable to open parameter file for Morse parameters");
483  return ier;
484  }
485 
486  ier = fscanf(fid, "%s \n%lf \n%lf \n%lf \n%lf",
487  speciesNameString, /* element symbol */
488  &cutoff, /* cutoff distance in angstroms */
489  &epsilon, /* Morse epsilon in eV */
490  &C, /* Morse C in 1/Angstroms */
491  &Rzero /* Morse Rzero in Angstroms */
492  );
493  fclose(fid);
494 
495  /* check that we read the right number of parameters */
496  if (5 != ier)
497  {
498  ier = TRUE;
499  LOG_ERROR("Unable to read all parameters");
500  return ier;
501  }
502 
503  /* register species */
504  speciesName = KIM_SpeciesName_FromString(speciesNameString);
506  modelDriverCreate,
507  speciesName,
508  SPECCODE);
509  if (ier == TRUE)
510  {
511  LOG_ERROR("Unable to set species code for Ar.");
512  return ier;
513  }
514 
515 
516  /* allocate buffer */
517  buffer = (struct model_buffer*) malloc(sizeof(struct model_buffer));
518  if (NULL == buffer)
519  {
520  ier = TRUE;
521  LOG_ERROR("malloc");
522  return ier;
523  }
524 
525  /* setup buffer */
526  /* set value of parameters */
528  buffer->cutoff = cutoff;
529  buffer->cutsq = (cutoff)*(cutoff);
531  buffer->halfListHint = 0;
532  buffer->epsilon = epsilon;
533  buffer->C = C;
534  buffer->Rzero = Rzero;
535  /* set value of parameter shift */
536  dummy = 0.0;
537  /* call calc_phi with r=cutoff and shift=0.0 */
538  calc_phi(&(buffer->epsilon),
539  &(buffer->C),
540  &(buffer->Rzero),
541  &dummy,
542  cutoff, cutoff, &(buffer->shift));
543  /* set shift to -shift */
544  buffer->shift = -buffer->shift;
545 
546  /* end setup buffer */
547 
548  /* store in model buffer */
550  (void*) buffer);
551 
552  /* store model cutoff in KIM object */
554  modelDriverCreate,
557  modelDriverCreate,
558  1,
559  &(buffer->cutoff),
561  &(buffer->halfListHint));
562 
563  return FALSE;
564 }
565 
566 /* refresh function */
567 static int refresh(KIM_ModelRefresh * const modelRefresh)
568 {
569  /* Local variables */
570  struct model_buffer* buffer;
571 
572  /* get model buffer from KIM object */
573  KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
574 
576  modelRefresh, &(buffer->influenceDistance));
578  modelRefresh,
579  1,
580  &(buffer->cutoff),
582  &(buffer->halfListHint));
583 
584  return FALSE;
585 }
586 
587 /* destroy function */
588 static int destroy(KIM_ModelDestroy * const modelDestroy)
589 {
590  /* Local variables */
591  struct model_buffer* buffer;
592  int ier;
593 
594  /* get model buffer from KIM object */
595  KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
596 
597  /* free the buffer */
598  free(buffer);
599 
600  ier = FALSE;
601  return ier;
602 }
603 
604 /* compute arguments create routine */
607  KIM_ModelCompute const * const modelCompute,
608  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
609 {
610  int ier;
611  /* register arguments */
613  modelComputeArgumentsCreate,
616  ||
618  modelComputeArgumentsCreate,
621  ||
623  modelComputeArgumentsCreate,
626  if (ier == TRUE)
627  {
628  LOG_ERROR("Unable to set argument supportStatus.");
629  return TRUE;
630  }
631  else
632  {
633  return FALSE;
634  }
635 }
636 
637 /* compue arguments destroy routine */
640  KIM_ModelCompute const * const modelCompute,
641  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
642 {
643  /* nothing to be done */
644 
645  return FALSE;
646 }
void KIM_ModelDriverCreate_GetNumberOfParameterFiles(KIM_ModelDriverCreate *const modelDriverCreate, int *const numberOfParameterFiles)
#define SPECCODE
KIM_SupportStatus const KIM_SUPPORT_STATUS_optional
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_coordinates
int KIM_ModelDriverCreate_SetUnits(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LengthUnit const lengthUnit, KIM_EnergyUnit const energyUnit, KIM_ChargeUnit const chargeUnit, KIM_TemperatureUnit const temperatureUnit, KIM_TimeUnit const timeUnit)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy
void KIM_ModelDestroy_GetModelBufferPointer(KIM_ModelDestroy const *const modelDestroy, void **const ptr)
static int compute_arguments_destroy(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
KIM_LanguageName const KIM_LANGUAGE_NAME_c
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialForces
struct KIM_ModelRefresh KIM_ModelRefresh
int KIM_ModelComputeArguments_GetArgumentPointerDouble(KIM_ModelComputeArguments const *const modelComputeArguments, KIM_ComputeArgumentName const computeArgumentName, double **const ptr)
KIM_SpeciesName KIM_SpeciesName_FromString(char const *const str)
void KIM_ModelDriverCreate_SetInfluenceDistancePointer(KIM_ModelDriverCreate *const modelDriverCreate, double *const influenceDistance)
void KIM_ModelRefresh_GetModelBufferPointer(KIM_ModelRefresh const *const modelRefresh, void **const ptr)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes
static void calc_phi_dphi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi, double *dphi)
KIM_TemperatureUnit const KIM_TEMPERATURE_UNIT_unused
int KIM_ModelDriverCreate_SetComputeArgumentsCreatePointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
ChargeUnit const C
KIM_TimeUnit const KIM_TIME_UNIT_unused
#define TRUE
struct KIM_ModelDestroy KIM_ModelDestroy
void KIM_ModelDriverCreate_SetNeighborListPointers(KIM_ModelDriverCreate *const modelDriverCreate, int const numberOfNeighborLists, double const *const cutoffs, int const *const paddingNeighborHints, int const *const halfListHints)
ComputeArgumentName const particleContributing
int KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate, KIM_ComputeArgumentName const computeArgumentName, KIM_SupportStatus const supportStatus)
struct KIM_ModelDriverCreate KIM_ModelDriverCreate
void KIM_ModelCompute_GetModelBufferPointer(KIM_ModelCompute const *const modelCompute, void **const ptr)
void() func()
Definition: KIM_func.h:39
static int compute(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArguments const *const modelComputeArguments)
#define LOG_ERROR(message)
int model_driver_create(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LengthUnit const requestedLengthUnit, KIM_EnergyUnit const requestedEnergyUnit, KIM_ChargeUnit const requestedChargeUnit, KIM_TemperatureUnit const requestedTemperatureUnit, KIM_TimeUnit const requestedTimeUnit)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles
struct KIM_ModelComputeArguments KIM_ModelComputeArguments
struct KIM_ModelComputeArgumentsCreate KIM_ModelComputeArgumentsCreate
int KIM_ModelDriverCreate_GetParameterFileName(KIM_ModelDriverCreate *const modelDriverCreate, int const index, char const **const parameterFileName)
int KIM_ModelDriverCreate_SetSpeciesCode(KIM_ModelDriverCreate *const modelDriverCreate, KIM_SpeciesName const speciesName, int const code)
#define FALSE
static int compute_arguments_create(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
int KIM_ModelDriverCreate_SetComputePointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
int KIM_ModelDriverCreate_SetComputeArgumentsDestroyPointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
int KIM_ModelDriverCreate_SetDestroyPointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
static void calc_phi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi)
int KIM_ModelComputeArguments_GetNeighborList(KIM_ModelComputeArguments const *const modelComputeArguments, int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle)
#define DIM
double influenceDistance
struct KIM_ModelComputeArgumentsDestroy KIM_ModelComputeArgumentsDestroy
struct buffer buffer
KIM_EnergyUnit const KIM_ENERGY_UNIT_eV
static int refresh(KIM_ModelRefresh *const modelRefresh)
void KIM_ModelRefresh_SetInfluenceDistancePointer(KIM_ModelRefresh *const modelRefresh, double *const influenceDistance)
KIM_ChargeUnit const KIM_CHARGE_UNIT_unused
ComputeArgumentName const particleSpeciesCodes
static int destroy(KIM_ModelDestroy *const modelDestroy)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_partialEnergy
void KIM_ModelDriverCreate_SetModelBufferPointer(KIM_ModelDriverCreate *const modelDriverCreate, void *const ptr)
void KIM_ModelRefresh_SetNeighborListPointers(KIM_ModelRefresh *const modelRefresh, int const numberOfNeighborLists, double const *const cutoffs, int const *const paddingNeighborHints, int const *const halfListHints)
int KIM_ModelDriverCreate_SetModelNumbering(KIM_ModelDriverCreate *const modelDriverCreate, KIM_Numbering const numbering)
KIM_LengthUnit const KIM_LENGTH_UNIT_A
int KIM_ModelDriverCreate_SetRefreshPointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
KIM_Numbering const KIM_NUMBERING_zeroBased
int KIM_ModelComputeArguments_GetArgumentPointerInteger(KIM_ModelComputeArguments const *const modelComputeArguments, KIM_ComputeArgumentName const computeArgumentName, int **const ptr)
KIM_ComputeArgumentName const KIM_COMPUTE_ARGUMENT_NAME_particleContributing
struct KIM_ModelCompute KIM_ModelCompute