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 /* 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_LogVerbosity.h"
50 #include "KIM_LengthUnit.h"
51 #include "KIM_EnergyUnit.h"
52 #include "KIM_ChargeUnit.h"
53 #include "KIM_TemperatureUnit.h"
54 #include "KIM_TimeUnit.h"
55 #include "KIM_LanguageName.h"
56 #include "KIM_SpeciesName.h"
57 #include "KIM_SupportStatus.h"
58 #include "KIM_ArgumentName.h"
59 #include "KIM_CallbackName.h"
60 #include "KIM_ModelDriverCreate.h"
61 #include "KIM_ModelRefresh.h"
62 #include "KIM_ModelCompute.h"
63 #include "KIM_ModelDestroy.h"
64 
65 #define TRUE 1
66 #define FALSE 0
67 
68 /******************************************************************************/
69 /* Below are the definitions for some constants */
70 /******************************************************************************/
71 #define DIM 3 /* dimensionality of space */
72 #define SPECCODE 1 /* internal species code */
73 
74 
75 /* Define prototype for Model Driver init */
77  KIM_ModelDriverCreate * const modelDriverCreate,
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 prototypes for destroy */
85 /* defined as static to avoid namespace clashes with other codes */
86 static int destroy(KIM_ModelDestroy * const modelDestroy);
87 
88 /* Define prototype for compute routine */
89 static int compute(KIM_ModelCompute const * const modelCompute);
90 
91 /* Define prototype for refresh routine */
92 static int refresh(KIM_ModelRefresh * const modelRefresh);
93 
94 /* Define prototypes for pair potential calculations */
95 static void calc_phi(double const * epsilon,
96  double const * C,
97  double const * Rzero,
98  double const * shift,
99  double const cutoff, double const r, double* phi);
100 
101 static void calc_phi_dphi(double const* epsilon,
102  double const* C,
103  double const* Rzero,
104  double const* shift,
105  double const cutoff, double const r,
106  double* phi, double* dphi);
107 
108 /* Define model_buffer structure */
109 struct model_buffer
110 {
111  double influenceDistance;
112  double cutsq;
113  double epsilon;
114  double C;
115  double Rzero;
116  double shift;
117 };
118 
119 
120 /* Calculate pair potential phi(r) */
121 static void calc_phi(double const* epsilon,
122  double const* C,
123  double const* Rzero,
124  double const* shift,
125  double const cutoff, double const r, double* phi)
126 {
127  /* local variables */
128  double ep;
129  double ep2;
130 
131  ep = exp(-(*C)*(r-*Rzero));
132  ep2 = ep*ep;
133 
134  if (r > cutoff)
135  {
136  /* Argument exceeds cutoff radius */
137  *phi = 0.0;
138  }
139  else
140  {
141  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
142  }
143 
144  return;
145 }
146 
147 /* Calculate pair potential phi(r) and its derivative dphi(r) */
148 static void calc_phi_dphi(double const* epsilon,
149  double const* C,
150  double const* Rzero,
151  double const* shift,
152  double const cutoff, double const r,
153  double* phi, double* dphi)
154 {
155  /* local variables */
156  double ep;
157  double ep2;
158 
159  ep = exp(-(*C)*(r-*Rzero));
160  ep2 = ep*ep;
161 
162  if (r > cutoff)
163  {
164  /* Argument exceeds cutoff radius */
165  *phi = 0.0;
166  *dphi = 0.0;
167  }
168  else
169  {
170  *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
171  *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
172  }
173 
174  return;
175 }
176 
177 /* compute function */
179 static int compute(KIM_ModelCompute const * const modelCompute)
180 {
181  /* local variables */
182  double R;
183  double Rsqij;
184  double phi;
185  double dphi;
186  double dEidr;
187  double Rij[DIM];
188  int ier;
189  int i;
190  int j;
191  int jj;
192  int k;
193  int const * neighListOfCurrentPart;
194  struct model_buffer* buffer;
195  int comp_energy;
196  int comp_force;
197  int comp_particleEnergy;
198 
199  int* nParts;
202  double cutoff;
203  double* cutsq;
204  double* epsilon;
205  double* C;
206  double* Rzero;
207  double* shift;
208  double* coords;
209  double* energy;
210  double* force;
211  double* particleEnergy;
212  int numOfPartNeigh;
213 
214  /* get buffer from KIM object */
215  KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
216 
217  /* unpack info from the buffer */
218  cutoff = buffer->influenceDistance;
219  cutsq = &(buffer->cutsq);
220  epsilon = &(buffer->epsilon);
221  C = &(buffer->C);
222  Rzero = &(buffer->Rzero);
223  shift = &(buffer->shift);
224 
225  ier =
227  modelCompute,
229  &nParts)
230  ||
232  modelCompute,
235  ||
237  modelCompute,
240  ||
242  modelCompute,
244  &coords)
245  ||
247  modelCompute,
249  &energy)
250  ||
252  modelCompute,
254  &force)
255  ||
257  modelCompute,
259  &particleEnergy);
260  if (ier)
261  {
262  LOG_ERROR("GetArgumentPointer");
263  return ier;
264  }
265 
266  comp_energy = (energy != 0);
267  comp_force = (force != 0);
268  comp_particleEnergy = (particleEnergy != 0);
269 
270  /* Check to be sure that the species are correct */
271 
272  ier = TRUE; /* assume an error */
273  for (i = 0; i < *nParts; ++i)
274  {
275  if ( SPECCODE != particleSpeciesCodes[i])
276  {
277  LOG_ERROR("Unexpected species code detected");
278  return ier;
279  }
280  }
281  ier = FALSE; /* everything is ok */
282 
283  /* initialize potential energies, forces, and virial term */
284  if (comp_particleEnergy)
285  {
286  for (i = 0; i < *nParts; ++i)
287  {
288  particleEnergy[i] = 0.0;
289  }
290  }
291  if (comp_energy)
292  {
293  *energy = 0.0;
294  }
295 
296  if (comp_force)
297  {
298  for (i = 0; i < *nParts; ++i)
299  {
300  for (k = 0; k < DIM; ++k)
301  {
302  force[i*DIM + k] = 0.0;
303  }
304  }
305  }
306 
307  /* Compute energy and forces */
308 
309  /* loop over particles and compute enregy and forces */
310  for (i = 0; i < *nParts; ++i)
311  {
312  if (particleContributing[i])
313  {
315  modelCompute, 0, i, &numOfPartNeigh, &neighListOfCurrentPart);
316  if (ier)
317  {
318  /* some sort of problem, exit */
319  LOG_ERROR("KIM_get_neigh");
320  ier = TRUE;
321  return ier;
322  }
323 
324  /* loop over the neighbors of particle i */
325  for (jj = 0; jj < numOfPartNeigh; ++ jj)
326  {
327  j = neighListOfCurrentPart[jj]; /* get neighbor ID */
328 
329  /* compute relative position vector and squared distance */
330  Rsqij = 0.0;
331  for (k = 0; k < DIM; ++k)
332  {
333  Rij[k] = coords[j*DIM + k] - coords[i*DIM + k];
334  /* compute squared distance */
335  Rsqij += Rij[k]*Rij[k];
336  }
337 
338  /* compute energy and force */
339  if (Rsqij < *cutsq)
340  {
341  /* particles are interacting ? */
342  R = sqrt(Rsqij);
343  if (comp_force)
344  {
345  /* compute pair potential and its derivative */
346  calc_phi_dphi(epsilon,
347  C,
348  Rzero,
349  shift,
350  cutoff, R, &phi, &dphi);
351 
352  /* compute dEidr */
353  dEidr = 0.5*dphi;
354  }
355  else
356  {
357  /* compute just pair potential */
358  calc_phi(epsilon,
359  C,
360  Rzero,
361  shift,
362  cutoff, R, &phi);
363  }
364 
365  /* contribution to energy */
366  if (comp_particleEnergy)
367  {
368  particleEnergy[i] += 0.5*phi;
369  }
370  if (comp_energy)
371  {
372  *energy += 0.5*phi;
373  }
374 
375  /* contribution to forces */
376  if (comp_force)
377  {
378  for (k = 0; k < DIM; ++k)
379  {
380  force[i*DIM + k] += dEidr*Rij[k]/R; /* accumulate force on i */
381  force[j*DIM + k] -= dEidr*Rij[k]/R; /* accumulate force on j */
382  }
383  }
384  } /* if Rsqij */
385  } /* loop on jj */
386  } /* if particleContributing */
387  } /* infinite while loop (terminated by break statements above) */
388 
389  /* everything is great */
390  ier = FALSE;
391 
392  return ier;
393 }
394 
395 /* Create function */
398  KIM_ModelDriverCreate * const modelDriverCreate,
399  KIM_LengthUnit const requestedLengthUnit,
400  KIM_EnergyUnit const requestedEnergyUnit,
401  KIM_ChargeUnit const requestedChargeUnit,
402  KIM_TemperatureUnit const requestedTemperatureUnit,
403  KIM_TimeUnit const requestedTimeUnit)
404 {
405  /* KIM variables */
406  int numberOfParameterFiles;
407  char const * paramfile1name;
408 
409  /* Local variables */
410  FILE* fid;
411  char speciesNameString[100];
412  KIM_SpeciesName speciesName;
413  double cutoff;
414  double epsilon;
415  double C;
416  double Rzero;
417  int ier;
418  double dummy;
419  struct model_buffer* buffer;
420 
421 
422  /* using fixed units */
423  ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
429  if (ier == TRUE)
430  {
431  LOG_ERROR("Problem setting units");
432  return ier;
433  }
434 
435  /* register arguments */
437  modelDriverCreate,
440  ||
442  modelDriverCreate,
445  ||
447  modelDriverCreate,
450  if (ier == TRUE)
451  {
452  LOG_ERROR("Unable to set argument supportStatus.");
453  return ier;
454  }
455 
456  /* store pointer to functions in KIM object */
457  KIM_ModelDriverCreate_SetDestroyPointer(modelDriverCreate,
459  (func *) destroy);
460  KIM_ModelDriverCreate_SetComputePointer(modelDriverCreate,
462  (func *) compute);
464  modelDriverCreate, KIM_LANGUAGE_NAME_c, (func *) refresh);
465 
466  /* get number of parameter files */
468  modelDriverCreate, &numberOfParameterFiles);
469  /* set paramfile1name */
470  if (numberOfParameterFiles != 1)
471  {
472  ier = TRUE;
473  LOG_ERROR("Incorrect number of parameter files.");
474  return ier;
475  }
477  modelDriverCreate,
478  0,
479  &paramfile1name);
480  if (ier == TRUE)
481  {
482  LOG_ERROR("Unable to get parameter file name.");
483  return ier;
484  }
485 
486  /* Read in model parameters from parameter file */
487  fid = fopen(paramfile1name, "r");
488  if (fid == NULL)
489  {
490  ier = TRUE;
491  LOG_ERROR("Unable to open parameter file for Morse parameters");
492  return ier;
493  }
494 
495  ier = fscanf(fid, "%s \n%lf \n%lf \n%lf \n%lf",
496  speciesNameString, /* element symbol */
497  &cutoff, /* cutoff distance in angstroms */
498  &epsilon, /* Morse epsilon in eV */
499  &C, /* Morse C in 1/Angstroms */
500  &Rzero /* Morse Rzero in Angstroms */
501  );
502  fclose(fid);
503 
504  /* check that we read the right number of parameters */
505  if (5 != ier)
506  {
507  ier = TRUE;
508  LOG_ERROR("Unable to read all parameters");
509  return ier;
510  }
511 
512  /* register species */
513  speciesName = KIM_SpeciesNameFromString(speciesNameString);
515  modelDriverCreate,
516  speciesName,
517  SPECCODE);
518  if (ier == TRUE)
519  {
520  LOG_ERROR("Unable to set species code for Ar.");
521  return ier;
522  }
523 
524 
525  /* allocate buffer */
526  buffer = (struct model_buffer*) malloc(sizeof(struct model_buffer));
527  if (NULL == buffer)
528  {
529  ier = TRUE;
530  LOG_ERROR("malloc");
531  return ier;
532  }
533 
534  /* setup buffer */
535  /* set value of parameters */
536  buffer->influenceDistance = cutoff;
537  buffer->cutsq = (cutoff)*(cutoff);
538  buffer->epsilon = epsilon;
539  buffer->C = C;
540  buffer->Rzero = Rzero;
541  /* set value of parameter shift */
542  dummy = 0.0;
543  /* call calc_phi with r=cutoff and shift=0.0 */
544  calc_phi(&(buffer->epsilon),
545  &(buffer->C),
546  &(buffer->Rzero),
547  &dummy,
548  cutoff, cutoff, &(buffer->shift));
549  /* set shift to -shift */
550  buffer->shift = -buffer->shift;
551 
552  /* end setup buffer */
553 
554  /* store in model buffer */
556  (void*) buffer);
557 
558  /* store model cutoff in KIM object */
560  modelDriverCreate,
561  &(buffer->influenceDistance));
563  modelDriverCreate, 1,
564  &(buffer->influenceDistance));
565 
566  return FALSE;
567 }
568 
569 /* refresh function */
570 static int refresh(KIM_ModelRefresh * const modelRefresh)
571 {
572  /* Local variables */
573  struct model_buffer* buffer;
574 
575  /* get model buffer from KIM object */
576  KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
577 
579  modelRefresh, &(buffer->influenceDistance));
581  modelRefresh, 1, &(buffer->influenceDistance));
582 
583  return FALSE;
584 }
585 
586 /* destroy function */
587 static int destroy(KIM_ModelDestroy * const modelDestroy)
588 {
589  /* Local variables */
590  struct model_buffer* buffer;
591  int ier;
592 
593  /* get model buffer from KIM object */
594  KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
595 
596  /* free the buffer */
597  free(buffer);
598 
599  ier = FALSE;
600  return ier;
601 }
KIM_SpeciesName KIM_SpeciesNameFromString(char const *const str)
void KIM_ModelDriverCreate_GetNumberOfParameterFiles(KIM_ModelDriverCreate *const modelDriverCreate, int *const numberOfParameterFiles)
#define SPECCODE
KIM_SupportStatus const KIM_SUPPORT_STATUS_optional
KIM_ArgumentName const KIM_ARGUMENT_NAME_particleContributing
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)
static int compute(KIM_ModelCompute const *const modelCompute)
KIM_ArgumentName const KIM_ARGUMENT_NAME_numberOfParticles
KIM_ChargeUnit const KIM_CHARGE_UNIT_e
void KIM_ModelDestroy_GetModelBufferPointer(KIM_ModelDestroy const *const modelDestroy, void **const ptr)
KIM_LanguageName const KIM_LANGUAGE_NAME_c
ArgumentName const particleContributing
struct KIM_ModelRefresh KIM_ModelRefresh
void KIM_ModelDriverCreate_SetInfluenceDistancePointer(KIM_ModelDriverCreate *const modelDriverCreate, double *const influenceDistance)
void KIM_ModelRefresh_GetModelBufferPointer(KIM_ModelRefresh const *const modelRefresh, void **const ptr)
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)
ChargeUnit const C
#define TRUE
int KIM_ModelDriverCreate_SetArgumentSupportStatus(KIM_ModelDriverCreate *const modelDriverCreate, KIM_ArgumentName const argumentName, KIM_SupportStatus const supportStatus)
void KIM_ModelRefresh_SetNeighborListCutoffsPointer(KIM_ModelRefresh *const modelRefresh, int const numberOfCutoffs, double const *const cutoffs)
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialParticleEnergy
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialForces
struct KIM_ModelDestroy KIM_ModelDestroy
struct KIM_ModelDriverCreate KIM_ModelDriverCreate
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 KIM_ModelDriverCreate_SetNeighborListCutoffsPointer(KIM_ModelDriverCreate *const modelDriverCreate, int const numberOfCutoffs, double const *const cutoffs)
void() func()
Definition: KIM_func.h:39
ArgumentName const particleSpeciesCodes
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_TimeUnit const KIM_TIME_UNIT_ps
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
int KIM_ModelCompute_GetArgumentPointerInteger(KIM_ModelCompute const *const modelCompute, KIM_ArgumentName const argumentName, int **const ptr)
int KIM_ModelDriverCreate_SetComputePointer(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)
#define DIM
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_ArgumentName const KIM_ARGUMENT_NAME_coordinates
KIM_ArgumentName const KIM_ARGUMENT_NAME_particleSpeciesCodes
static int destroy(KIM_ModelDestroy *const modelDestroy)
int KIM_ModelCompute_GetArgumentPointerDouble(KIM_ModelCompute const *const modelCompute, KIM_ArgumentName const argumentName, double **const ptr)
void KIM_ModelDriverCreate_SetModelBufferPointer(KIM_ModelDriverCreate *const modelDriverCreate, void *const ptr)
#define LOG_ERROR(message)
KIM_LengthUnit const KIM_LENGTH_UNIT_A
KIM_ArgumentName const KIM_ARGUMENT_NAME_partialEnergy
KIM_TemperatureUnit const KIM_TEMPERATURE_UNIT_K
int KIM_ModelDriverCreate_SetRefreshPointer(KIM_ModelDriverCreate *const modelDriverCreate, KIM_LanguageName const languageName, func *const fptr)
struct KIM_ModelCompute KIM_ModelCompute