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 */
101 struct model_buffer
102 {
103  double influenceDistance;
104  double cutsq;
105  double epsilon;
106  double C;
107  double Rzero;
108  double shift;
109 };
110 
111 
112 /* Calculate pair potential phi(r) */
113 static void calc_phi(double const* epsilon,
114  double const* C,
115  double const* Rzero,
116  double const* shift,
117  double const cutoff, double const r, double* phi)
118 {
119  /* local variables */
120  double ep;
121  double ep2;
122 
123  ep = exp(-(*C)*(r-*Rzero));
124  ep2 = ep*ep;
125 
126  if (r > cutoff)
127  {
128  /* Argument exceeds cutoff radius */
129  *phi = 0.0;
130  }
131  else
132  {
133  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
134  }
135 
136  return;
137 }
138 
139 /* Calculate pair potential phi(r) and its derivative dphi(r) */
140 static void calc_phi_dphi(double const* epsilon,
141  double const* C,
142  double const* Rzero,
143  double const* shift,
144  double const cutoff, double const r,
145  double* phi, double* dphi)
146 {
147  /* local variables */
148  double ep;
149  double ep2;
150 
151  ep = exp(-(*C)*(r-*Rzero));
152  ep2 = ep*ep;
153 
154  if (r > cutoff)
155  {
156  /* Argument exceeds cutoff radius */
157  *phi = 0.0;
158  *dphi = 0.0;
159  }
160  else
161  {
162  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
163  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
164  }
165 
166  return;
167 }
168 
169 /* compute function */
171 static int compute(
172  KIM_ModelCompute const * const modelCompute,
173  KIM_ModelComputeArguments const * const modelComputeArguments)
174 {
175  /* local variables */
176  double R;
177  double Rsqij;
178  double phi;
179  double dphi;
180  double dEidr;
181  double Rij[DIM];
182  int ier;
183  int i;
184  int j;
185  int jj;
186  int k;
187  int const * neighListOfCurrentPart;
188  struct model_buffer* buffer;
189  int comp_energy;
190  int comp_force;
191  int comp_particleEnergy;
192 
193  int* nParts;
196  double cutoff;
197  double* cutsq;
198  double* epsilon;
199  double* C;
200  double* Rzero;
201  double* shift;
202  double* coords;
203  double* energy;
204  double* force;
205  double* particleEnergy;
206  int numOfPartNeigh;
207 
208  /* get buffer from KIM object */
209  KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
210 
211  /* unpack info from the buffer */
212  cutoff = buffer->influenceDistance;
213  cutsq = &(buffer->cutsq);
214  epsilon = &(buffer->epsilon);
215  C = &(buffer->C);
216  Rzero = &(buffer->Rzero);
217  shift = &(buffer->shift);
218 
219  ier =
221  modelComputeArguments,
223  &nParts)
224  ||
226  modelComputeArguments,
229  ||
231  modelComputeArguments,
234  ||
236  modelComputeArguments,
238  &coords)
239  ||
241  modelComputeArguments,
243  &energy)
244  ||
246  modelComputeArguments,
248  &force)
249  ||
251  modelComputeArguments,
253  &particleEnergy);
254  if (ier)
255  {
256  LOG_ERROR("GetArgumentPointer");
257  return ier;
258  }
259 
260  comp_energy = (energy != 0);
261  comp_force = (force != 0);
262  comp_particleEnergy = (particleEnergy != 0);
263 
264  /* Check to be sure that the species are correct */
265 
266  ier = TRUE; /* assume an error */
267  for (i = 0; i < *nParts; ++i)
268  {
269  if ( SPECCODE != particleSpeciesCodes[i])
270  {
271  LOG_ERROR("Unexpected species code detected");
272  return ier;
273  }
274  }
275  ier = FALSE; /* everything is ok */
276 
277  /* initialize potential energies, forces, and virial term */
278  if (comp_particleEnergy)
279  {
280  for (i = 0; i < *nParts; ++i)
281  {
282  particleEnergy[i] = 0.0;
283  }
284  }
285  if (comp_energy)
286  {
287  *energy = 0.0;
288  }
289 
290  if (comp_force)
291  {
292  for (i = 0; i < *nParts; ++i)
293  {
294  for (k = 0; k < DIM; ++k)
295  {
296  force[i*DIM + k] = 0.0;
297  }
298  }
299  }
300 
301  /* Compute energy and forces */
302 
303  /* loop over particles and compute enregy and forces */
304  for (i = 0; i < *nParts; ++i)
305  {
306  if (particleContributing[i])
307  {
309  modelComputeArguments,
310  0, i, &numOfPartNeigh, &neighListOfCurrentPart);
311  if (ier)
312  {
313  /* some sort of problem, exit */
314  LOG_ERROR("KIM_get_neigh");
315  ier = TRUE;
316  return ier;
317  }
318 
319  /* loop over the neighbors of particle i */
320  for (jj = 0; jj < numOfPartNeigh; ++ jj)
321  {
322  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
323 
324  /* compute relative position vector and squared distance */
325  Rsqij = 0.0;
326  for (k = 0; k < DIM; ++k)
327  {
328  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
329  /* compute squared distance */
330  Rsqij += Rij[k]*Rij[k];
331  }
332 
333  /* compute energy and force */
334  if (Rsqij < *cutsq)
335  {
336  /* particles are interacting ? */
337  R = sqrt(Rsqij);
338  if (comp_force)
339  {
340  /* compute pair potential and its derivative */
341  calc_phi_dphi(epsilon,
342  C,
343  Rzero,
344  shift,
345  cutoff, R, &phi, &dphi);
346 
347  /* compute dEidr */
348  dEidr = 0.5*dphi;
349  }
350  else
351  {
352  /* compute just pair potential */
353  calc_phi(epsilon,
354  C,
355  Rzero,
356  shift,
357  cutoff, R, &phi);
358  }
359 
360  /* contribution to energy */
361  if (comp_particleEnergy)
362  {
363  particleEnergy[i] += 0.5*phi;
364  }
365  if (comp_energy)
366  {
367  *energy += 0.5*phi;
368  }
369 
370  /* contribution to forces */
371  if (comp_force)
372  {
373  for (k = 0; k < DIM; ++k)
374  {
375  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
376  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */
377  }
378  }
379  } /* if Rsqij */
380  } /* loop on jj */
381  } /* if particleContributing */
382  } /* infinite while loop (terminated by break statements above) */
383 
384  /* everything is great */
385  ier = FALSE;
386 
387  return ier;
388 }
389 
390 /* Create function */
393  KIM_ModelDriverCreate * const modelDriverCreate,
394  KIM_LengthUnit const requestedLengthUnit,
395  KIM_EnergyUnit const requestedEnergyUnit,
396  KIM_ChargeUnit const requestedChargeUnit,
397  KIM_TemperatureUnit const requestedTemperatureUnit,
398  KIM_TimeUnit const requestedTimeUnit)
399 {
400  /* KIM variables */
401  int numberOfParameterFiles;
402  char const * paramfile1name;
403 
404  /* Local variables */
405  FILE* fid;
406  char speciesNameString[100];
407  KIM_SpeciesName speciesName;
408  double cutoff;
409  double epsilon;
410  double C;
411  double Rzero;
412  int ier;
413  double dummy;
414  struct model_buffer* buffer;
415 
416 
417  /* using fixed units */
418  ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
424  if (ier == TRUE)
425  {
426  LOG_ERROR("Problem setting units");
427  return ier;
428  }
429 
430  ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate,
432  if (ier == TRUE)
433  {
434  LOG_ERROR("Unable to set numbering");
435  return ier;
436  }
437 
438  /* store pointer to functions in KIM object */
439  KIM_ModelDriverCreate_SetDestroyPointer(modelDriverCreate,
441  (func *) destroy);
448  KIM_ModelDriverCreate_SetComputePointer(modelDriverCreate,
450  (func *) compute);
452  modelDriverCreate, KIM_LANGUAGE_NAME_c, (func *) refresh);
453 
454  /* get number of parameter files */
456  modelDriverCreate, &numberOfParameterFiles);
457  /* set paramfile1name */
458  if (numberOfParameterFiles != 1)
459  {
460  ier = TRUE;
461  LOG_ERROR("Incorrect number of parameter files.");
462  return ier;
463  }
465  modelDriverCreate,
466  0,
467  &paramfile1name);
468  if (ier == TRUE)
469  {
470  LOG_ERROR("Unable to get parameter file name.");
471  return ier;
472  }
473 
474  /* Read in model parameters from parameter file */
475  fid = fopen(paramfile1name, "r");
476  if (fid == NULL)
477  {
478  ier = TRUE;
479  LOG_ERROR("Unable to open parameter file for Morse parameters");
480  return ier;
481  }
482 
483  ier = fscanf(fid, "%s \n%lf \n%lf \n%lf \n%lf",
484  speciesNameString, /* element symbol */
485  &cutoff, /* cutoff distance in angstroms */
486  &epsilon, /* Morse epsilon in eV */
487  &C, /* Morse C in 1/Angstroms */
488  &Rzero /* Morse Rzero in Angstroms */
489  );
490  fclose(fid);
491 
492  /* check that we read the right number of parameters */
493  if (5 != ier)
494  {
495  ier = TRUE;
496  LOG_ERROR("Unable to read all parameters");
497  return ier;
498  }
499 
500  /* register species */
501  speciesName = KIM_SpeciesName_FromString(speciesNameString);
503  modelDriverCreate,
504  speciesName,
505  SPECCODE);
506  if (ier == TRUE)
507  {
508  LOG_ERROR("Unable to set species code for Ar.");
509  return ier;
510  }
511 
512 
513  /* allocate buffer */
514  buffer = (struct model_buffer*) malloc(sizeof(struct model_buffer));
515  if (NULL == buffer)
516  {
517  ier = TRUE;
518  LOG_ERROR("malloc");
519  return ier;
520  }
521 
522  /* setup buffer */
523  /* set value of parameters */
524  buffer->influenceDistance = cutoff;
525  buffer->cutsq = (cutoff)*(cutoff);
526  buffer->epsilon = epsilon;
527  buffer->C = C;
528  buffer->Rzero = Rzero;
529  /* set value of parameter shift */
530  dummy = 0.0;
531  /* call calc_phi with r=cutoff and shift=0.0 */
532  calc_phi(&(buffer->epsilon),
533  &(buffer->C),
534  &(buffer->Rzero),
535  &dummy,
536  cutoff, cutoff, &(buffer->shift));
537  /* set shift to -shift */
538  buffer->shift = -buffer->shift;
539 
540  /* end setup buffer */
541 
542  /* store in model buffer */
544  (void*) buffer);
545 
546  /* store model cutoff in KIM object */
548  modelDriverCreate,
549  &(buffer->influenceDistance));
551  modelDriverCreate, 1,
552  &(buffer->influenceDistance));
553 
554  return FALSE;
555 }
556 
557 /* refresh function */
558 static int refresh(KIM_ModelRefresh * const modelRefresh)
559 {
560  /* Local variables */
561  struct model_buffer* buffer;
562 
563  /* get model buffer from KIM object */
564  KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
565 
567  modelRefresh, &(buffer->influenceDistance));
569  modelRefresh, 1, &(buffer->influenceDistance));
570 
571  return FALSE;
572 }
573 
574 /* destroy function */
575 static int destroy(KIM_ModelDestroy * const modelDestroy)
576 {
577  /* Local variables */
578  struct model_buffer* buffer;
579  int ier;
580 
581  /* get model buffer from KIM object */
582  KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
583 
584  /* free the buffer */
585  free(buffer);
586 
587  ier = FALSE;
588  return ier;
589 }
590 
591 /* compute arguments create routine */
594  KIM_ModelCompute const * const modelCompute,
595  KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
596 {
597  int ier;
598  /* register arguments */
600  modelComputeArgumentsCreate,
603  ||
605  modelComputeArgumentsCreate,
608  ||
610  modelComputeArgumentsCreate,
613  if (ier == TRUE)
614  {
615  LOG_ERROR("Unable to set argument supportStatus.");
616  return TRUE;
617  }
618  else
619  {
620  return FALSE;
621  }
622 }
623 
624 /* compue arguments destroy routine */
627  KIM_ModelCompute const * const modelCompute,
628  KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
629 {
630  /* nothing to be done */
631 
632  return FALSE;
633 }
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
void KIM_ModelRefresh_SetNeighborListCutoffsPointer(KIM_ModelRefresh *const modelRefresh, int const numberOfCutoffs, double const *const cutoffs)
struct KIM_ModelDestroy KIM_ModelDestroy
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 KIM_ModelDriverCreate_SetNeighborListCutoffsPointer(KIM_ModelDriverCreate *const modelDriverCreate, int const numberOfCutoffs, double const *const cutoffs)
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
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)
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