1 /*
2 *
3 * CDDL HEADER START
4 *
5 * The contents of this file are subject to the terms of the Common Development
6 * and Distribution License Version 1.0 (the "License").
7 *
8 * You can obtain a copy of the license at
9 * http://www.opensource.org/licenses/CDDL-1.0. See the License for the
10 * specific language governing permissions and limitations under the License.
11 *
12 * When distributing Covered Code, include this CDDL HEADER in each file and
13 * include the License file in a prominent location with the name LICENSE.CDDL.
14 * If applicable, add the following below this CDDL HEADER, with the fields
15 * enclosed by brackets "[]" replaced with your own identifying information:
16 *
17 * Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
18 *
19 * CDDL HEADER END
20 *
21
22 *
23 * Copyright (c) 2013, Regents of the University of Minnesota.
24 * All rights reserved.
25 *
26 * Contributors:
27 * Ryan S. Elliott
28 * Ellad B. Tadmor
29 * Valeriu Smirichinski
30 * Stephen M. Whalen
31 * Yaser Afshar
32 *
33 */
34
35 /*******************************************************************************
36 *
37 * Pair_Morse_Shifted
38 *
39 * Morse pair potential KIM Model Driver
40 * shifted to have zero energy at the cutoff radius
41 *
42 * Language: C
43 *
44 *******************************************************************************/
45
46
47 #include "KIM_ModelDriverHeaders.h"
48 #include <math.h>
49 #include <stdio.h>
50 #include <stdlib.h>
51 #include <string.h>
52
53 #include "KIM_LogMacros.h"
54
55 #define TRUE 1
56 #define FALSE 0
57 #define ONE 1.0
58
59
60 /******************************************************************************
61 * Below are the definitions and values of all Model parameters
62 *******************************************************************************/
63 #define DIM 3 /* dimensionality of space */
64 #define SPECCODE 1 /* internal species code */
65
66
67 /* Define prototypes for Model Driver init */
68 /**/
69 /* Define prototype for Model Driver initialization */
70 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
71 KIM_LengthUnit const requestedLengthUnit,
72 KIM_EnergyUnit const requestedEnergyUnit,
73 KIM_ChargeUnit const requestedChargeUnit,
74 KIM_TemperatureUnit const requestedTemperatureUnit,
75 KIM_TimeUnit const requestedTimeUnit);
76
77 /* Define prototype for destroy (defined as static to avoid namespace clashes
78 * with other codes) */
79 static int destroy(KIM_ModelDestroy * const modelDestroy);
80
81 /* Define prototype for compute routine */
82 static int
83 compute(KIM_ModelCompute const * const modelCompute,
84 KIM_ModelComputeArguments const * const modelComputeArguments);
85 static int compute_arguments_create(
86 KIM_ModelCompute const * const modelCompute,
87 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
88 static int compute_arguments_destroy(
89 KIM_ModelCompute const * const modelCompute,
90 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
91
92 /* Define prototype for refresh routine */
93 static int refresh(KIM_ModelRefresh * const modelRefresh);
94 /**/
95 static void calc_phi(double * epsilon,
96 double * C,
97 double * Rzero,
98 double * shift,
99 double * cutoff,
100 double r,
101 double * phi);
102
103 static void calc_phi_dphi(double * epsilon,
104 double * C,
105 double * Rzero,
106 double * shift,
107 double * cutoff,
108 double r,
109 double * phi,
110 double * dphi);
111
112 static void calc_phi_d2phi(double * epsilon,
113 double * C,
114 double * Rzero,
115 double * shift,
116 double * cutoff,
117 double r,
118 double * phi,
119 double * dphi,
120 double * d2phi);
121
122 /* Define model_buffer structure */
123 struct model_buffer
124 {
125 int modelWillNotRequestNeighborsOfNoncontributingParticles;
126
127 double influenceDistance;
128 double Pcutoff;
129 double cutsq;
130 double epsilon;
131 double C;
132 double Rzero;
133 double shift;
134 };
135
136
137 /* Calculate pair potential phi(r) */
calc_phi(double * epsilon,double * C,double * Rzero,double * shift,double * cutoff,double r,double * phi)138 static void calc_phi(double * epsilon,
139 double * C,
140 double * Rzero,
141 double * shift,
142 double * cutoff,
143 double r,
144 double * phi)
145 {
146 if (r > *cutoff)
147 {
148 /* Argument exceeds cutoff radius */
149 *phi = 0.0;
150 return;
151 }
152
153 /* local variables */
154 double ep;
155 double ep2;
156
157 ep = exp(-(*C) * (r - *Rzero));
158 ep2 = ep * ep;
159
160 *phi = (*epsilon) * (-ep2 + 2.0 * ep) + *shift;
161
162 return;
163 }
164
165 /* Calculate pair potential phi(r) and its derivative dphi(r) */
calc_phi_dphi(double * epsilon,double * C,double * Rzero,double * shift,double * cutoff,double r,double * phi,double * dphi)166 static void calc_phi_dphi(double * epsilon,
167 double * C,
168 double * Rzero,
169 double * shift,
170 double * cutoff,
171 double r,
172 double * phi,
173 double * dphi)
174 {
175 if (r > *cutoff)
176 {
177 /* Argument exceeds cutoff radius */
178 *phi = 0.0;
179 *dphi = 0.0;
180 return;
181 }
182
183 /* local variables */
184 double ep;
185 double ep2;
186
187 ep = exp(-(*C) * (r - *Rzero));
188 ep2 = ep * ep;
189
190 *phi = (*epsilon) * (-ep2 + 2.0 * ep) + *shift;
191 *dphi = 2.0 * (*epsilon) * (*C) * (-ep + ep2);
192
193 return;
194 }
195
196 /*
197 Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r),
198 d2phi(r)
199 */
calc_phi_d2phi(double * epsilon,double * C,double * Rzero,double * shift,double * cutoff,double r,double * phi,double * dphi,double * d2phi)200 static void calc_phi_d2phi(double * epsilon,
201 double * C,
202 double * Rzero,
203 double * shift,
204 double * cutoff,
205 double r,
206 double * phi,
207 double * dphi,
208 double * d2phi)
209 {
210 if (r > *cutoff)
211 {
212 /* Argument exceeds cutoff radius */
213 *phi = 0.0;
214 *dphi = 0.0;
215 *d2phi = 0.0;
216 return;
217 }
218
219 /* local variables */
220 double ep;
221 double ep2;
222
223 ep = exp(-(*C) * (r - *Rzero));
224 ep2 = ep * ep;
225
226 *phi = (*epsilon) * (-ep2 + 2.0 * ep) + *shift;
227 *dphi = 2.0 * (*epsilon) * (*C) * (-ep + ep2);
228 *d2phi = 2.0 * (*epsilon) * (*C) * (*C) * (ep - 2.0 * ep2);
229
230 return;
231 }
232
233
234 /* ######## Function Definitions ######## */
235
236 /* Primary compute function*/
237 #undef KIM_LOGGER_FUNCTION_NAME
238 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelCompute_LogEntry
239 #undef KIM_LOGGER_OBJECT_NAME
240 #define KIM_LOGGER_OBJECT_NAME modelCompute
241 static int
compute(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArguments const * const modelComputeArguments)242 compute(KIM_ModelCompute const * const modelCompute,
243 KIM_ModelComputeArguments const * const modelComputeArguments)
244 {
245 /* local variables */
246 double R;
247 double R_pairs[2];
248 double * pR_pairs = &(R_pairs[0]);
249 double Rsqij;
250 double phi;
251 double dphi;
252 double d2phi;
253 double dEidr = 0.0;
254 double d2Eidr = 0.0;
255 double Rij[DIM];
256 double * pRij = &(Rij[0]);
257 double Rij_pairs[2][3];
258 double * pRij_pairs = &(Rij_pairs[0][0]);
259 double v;
260 double vir[6];
261 int ier;
262 int i;
263 int i_pairs[2];
264 int * pi_pairs = &(i_pairs[0]);
265 int j;
266 int j_pairs[2];
267 int * pj_pairs = &(j_pairs[0]);
268 int jj;
269 int k;
270 int const * neighListOfCurrentAtom;
271 struct model_buffer * buffer;
272 int comp_energy;
273 int comp_force;
274 int comp_particleEnergy;
275 int comp_virial;
276 int comp_particleVirial;
277 int comp_process_dEdr;
278 int comp_process_d2Edr2;
279
280 int * nAtoms;
281 int * particleSpeciesCodes;
282 int * particleContributing;
283 double * cutoff;
284 double * cutsq;
285 double * epsilon;
286 double * C;
287 double * Rzero;
288 double * shift;
289 double * coords;
290 double * energy;
291 double * force;
292 double * particleEnergy;
293 double * virial;
294 double * particleVirial;
295 int numOfAtomNeigh;
296
297 /* Get model buffer */
298 KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
299
300 /* unpack info from the buffer */
301 cutoff = &(buffer->Pcutoff);
302 cutsq = &(buffer->cutsq);
303 epsilon = &(buffer->epsilon);
304 C = &(buffer->C);
305 Rzero = &(buffer->Rzero);
306 shift = &(buffer->shift);
307
308 /* check to see if we have been asked to compute the forces, */
309 /* particleEnergy, and d1Edr */
310 LOG_DEBUG("Checking if call backs are present.");
311 KIM_ModelComputeArguments_IsCallbackPresent(
312 modelComputeArguments,
313 KIM_COMPUTE_CALLBACK_NAME_ProcessDEDrTerm,
314 &comp_process_dEdr);
315 KIM_ModelComputeArguments_IsCallbackPresent(
316 modelComputeArguments,
317 KIM_COMPUTE_CALLBACK_NAME_ProcessD2EDr2Term,
318 &comp_process_d2Edr2);
319
320 LOG_DEBUG("Getting data pointers");
321 ier = KIM_ModelComputeArguments_GetArgumentPointerInteger(
322 modelComputeArguments,
323 KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles,
324 &nAtoms)
325 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
326 modelComputeArguments,
327 KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes,
328 &particleSpeciesCodes)
329 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
330 modelComputeArguments,
331 KIM_COMPUTE_ARGUMENT_NAME_particleContributing,
332 &particleContributing)
333 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
334 modelComputeArguments,
335 KIM_COMPUTE_ARGUMENT_NAME_coordinates,
336 &coords)
337 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
338 modelComputeArguments,
339 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
340 &energy)
341 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
342 modelComputeArguments,
343 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
344 &force)
345 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
346 modelComputeArguments,
347 KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
348 &particleEnergy)
349 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
350 modelComputeArguments,
351 KIM_COMPUTE_ARGUMENT_NAME_partialVirial,
352 &virial)
353 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
354 modelComputeArguments,
355 KIM_COMPUTE_ARGUMENT_NAME_partialParticleVirial,
356 &particleVirial);
357 if (ier == TRUE)
358 {
359 LOG_ERROR("Could not retrieve compute arguments");
360 return ier;
361 }
362
363 /* Check to see what we've been asked to compute */
364 comp_energy = (energy != 0);
365 comp_force = (force != 0);
366 comp_particleEnergy = (particleEnergy != 0);
367 comp_virial = (virial != 0);
368 comp_particleVirial = (particleVirial != 0);
369
370 /* Check to be sure that the species are correct */
371 ier = TRUE;
372 for (i = 0; i < *nAtoms; ++i)
373 {
374 if (SPECCODE != particleSpeciesCodes[i])
375 {
376 LOG_ERROR("Unexpected species code detected");
377 return ier;
378 }
379 }
380 ier = FALSE;
381
382 if (comp_particleEnergy)
383 {
384 for (i = 0; i < *nAtoms; ++i) { particleEnergy[i] = 0.0; }
385 }
386 if (comp_energy) { *energy = 0.0; }
387
388 if (comp_force)
389 {
390 for (i = 0; i < *nAtoms; ++i)
391 {
392 for (k = 0; k < DIM; ++k) { force[i * DIM + k] = 0.0; }
393 }
394 }
395
396 if (comp_virial)
397 {
398 for (i = 0; i < 6; ++i) virial[i] = 0.0;
399 }
400
401 if (comp_particleVirial)
402 {
403 for (i = 0; i < *nAtoms; ++i)
404 {
405 for (k = 0; k < 6; ++k) particleVirial[i * 6 + k] = 0.0;
406 }
407 }
408
409 /* Compute energy and forces */
410
411 /* loop over particles and compute enregy and forces */
412 LOG_DEBUG("Starting main compute loop");
413 for (i = 0; i < *nAtoms; ++i)
414 {
415 if (particleContributing[i])
416 {
417 ier = KIM_ModelComputeArguments_GetNeighborList(modelComputeArguments,
418 0,
419 i,
420 &numOfAtomNeigh,
421 &neighListOfCurrentAtom);
422 if (ier)
423 {
424 /* some sort of problem, exit */
425 LOG_ERROR("GetNeighborList failed");
426 ier = TRUE;
427 return ier;
428 }
429
430 /* loop over the neighbors of atom i */
431 for (jj = 0; jj < numOfAtomNeigh; ++jj)
432 {
433 /* get neighbor ID */
434 j = neighListOfCurrentAtom[jj];
435
436 if (!(particleContributing[j] && j < i)) /* short-circuit half-list */
437 {
438 /* compute relative position vector and squared distance */
439 Rsqij = 0.0;
440 for (k = 0; k < DIM; ++k)
441 {
442 Rij[k] = coords[j * DIM + k] - coords[i * DIM + k];
443
444 /* compute squared distance */
445 Rsqij += Rij[k] * Rij[k];
446 }
447
448 /* compute energy and force */
449 if (Rsqij < *cutsq) /* particles are interacting ? */
450 {
451 R = sqrt(Rsqij);
452 if (comp_process_d2Edr2)
453 {
454 /* compute pair potential and its derivatives */
455 calc_phi_d2phi(
456 epsilon, C, Rzero, shift, cutoff, R, &phi, &dphi, &d2phi);
457
458 /* compute dEidr */
459 /* Half mode -- double contribution */
460 if (particleContributing[j])
461 {
462 dEidr = dphi;
463 d2Eidr = d2phi;
464 }
465 else
466 {
467 dEidr = 0.5 * dphi;
468 d2Eidr = 0.5 * d2phi;
469 }
470 }
471 else if (comp_force || comp_process_dEdr ||
472 comp_virial || comp_particleVirial)
473 {
474 /* compute pair potential and its derivative */
475 calc_phi_dphi(epsilon, C, Rzero, shift, cutoff, R, &phi, &dphi);
476
477 /* compute dEidr */
478 /* Half mode -- double contribution */
479 if (particleContributing[j]) { dEidr = dphi; }
480 else
481 {
482 dEidr = 0.5 * dphi;
483 }
484 }
485 else
486 {
487 /* compute just pair potential */
488 calc_phi(epsilon, C, Rzero, shift, cutoff, R, &phi);
489 }
490
491 /* contribution to energy */
492 if (comp_particleEnergy)
493 {
494 particleEnergy[i] += 0.5 * phi;
495 if (particleContributing[j]) { particleEnergy[j] += 0.5 * phi; }
496 }
497 if (comp_energy)
498 {
499 if (particleContributing[j]) { *energy += phi; }
500 else
501 {
502 *energy += 0.5 * phi;
503 }
504 }
505 /* contribution to process_dEdr */
506 if (comp_process_dEdr)
507 {
508 ier = KIM_ModelComputeArguments_ProcessDEDrTerm(
509 modelComputeArguments, dEidr, R, pRij, i, j);
510 }
511 /* contribution to virial and particle virial */
512 if (comp_virial || comp_particleVirial)
513 {
514 /* Virial has 6 components and is stored as a 6-element
515 vector in the following order: xx, yy, zz, yz, xz, xy. */
516 v = dEidr / R;
517
518 vir[0] = v * Rij[0] * Rij[0];
519 vir[1] = v * Rij[1] * Rij[1];
520 vir[2] = v * Rij[2] * Rij[2];
521 vir[3] = v * Rij[1] * Rij[2];
522 vir[4] = v * Rij[0] * Rij[2];
523 vir[5] = v * Rij[0] * Rij[1];
524
525 if (comp_virial)
526 {
527 virial[0] += vir[0];
528 virial[1] += vir[1];
529 virial[2] += vir[2];
530 virial[3] += vir[3];
531 virial[4] += vir[4];
532 virial[5] += vir[5];
533 }
534
535 if (comp_particleVirial)
536 {
537 vir[0] *= 0.5;
538 vir[1] *= 0.5;
539 vir[2] *= 0.5;
540 vir[3] *= 0.5;
541 vir[4] *= 0.5;
542 vir[5] *= 0.5;
543
544 k = i * 6;
545 particleVirial[k + 0] += vir[0];
546 particleVirial[k + 1] += vir[1];
547 particleVirial[k + 2] += vir[2];
548 particleVirial[k + 3] += vir[3];
549 particleVirial[k + 4] += vir[4];
550 particleVirial[k + 5] += vir[5];
551
552 k = j * 6;
553 particleVirial[k + 0] += vir[0];
554 particleVirial[k + 1] += vir[1];
555 particleVirial[k + 2] += vir[2];
556 particleVirial[k + 3] += vir[3];
557 particleVirial[k + 4] += vir[4];
558 particleVirial[k + 5] += vir[5];
559 }
560 }
561
562 /* contribution to process_d2Edr2 */
563 if (comp_process_d2Edr2)
564 {
565 R_pairs[0] = R_pairs[1] = R;
566 Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
567 Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
568 Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
569 i_pairs[0] = i_pairs[1] = i;
570 j_pairs[0] = j_pairs[1] = j;
571
572 ier = KIM_ModelComputeArguments_ProcessD2EDr2Term(
573 modelComputeArguments,
574 d2Eidr,
575 pR_pairs,
576 pRij_pairs,
577 pi_pairs,
578 pj_pairs);
579 }
580
581 /* contribution to forces */
582 if (comp_force)
583 {
584 for (k = 0; k < DIM; ++k)
585 { /* accumulate force on atom i */
586 force[i * DIM + k] += dEidr * Rij[k] / R;
587 /* accumulate force on atom j */
588 force[j * DIM + k] -= dEidr * Rij[k] / R;
589 }
590 }
591 }
592 }
593 } /* loop on jj */
594 } /* if contributing */
595 } /* loop on i */
596 LOG_DEBUG("Finished compute loop");
597
598 /* everything is great */
599 ier = FALSE;
600
601 return ier;
602 }
603
604 /* Define prototype for Model Driver initialization */
605 #undef KIM_LOGGER_FUNCTION_NAME
606 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDriverCreate_LogEntry
607 #undef KIM_LOGGER_OBJECT_NAME
608 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
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)609 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
610 KIM_LengthUnit const requestedLengthUnit,
611 KIM_EnergyUnit const requestedEnergyUnit,
612 KIM_ChargeUnit const requestedChargeUnit,
613 KIM_TemperatureUnit const requestedTemperatureUnit,
614 KIM_TimeUnit const requestedTimeUnit)
615 {
616 /* Create Model buffer */
617 struct model_buffer * buffer;
618
619 /* KIM variables */
620 int numberOfParameterFiles;
621 char const * paramfile1name; /* This driver only actually supports a single
622 parameter file */
623 char speciesNameString[100];
624 KIM_SpeciesName speciesName;
625
626 /* Model Parameters */
627 double a;
628 double cutoff;
629 double epsilon;
630 double C;
631 double Rzero;
632 int ier;
633
634 KIM_LengthUnit fromLength;
635 KIM_EnergyUnit fromEnergy;
636 KIM_ChargeUnit fromCharge;
637 KIM_TemperatureUnit fromTemperature;
638 KIM_TimeUnit fromTime;
639
640 /* Misc */
641 FILE * fid;
642 double dummy;
643 double convertLength = 1.0;
644 double convertEnergy = 1.0;
645
646 ier = FALSE;
647
648 /* Get number of parameter files and error out if more than one is listed */
649 KIM_ModelDriverCreate_GetNumberOfParameterFiles(modelDriverCreate,
650 &numberOfParameterFiles);
651 if (numberOfParameterFiles != 1)
652 {
653 ier = TRUE;
654 LOG_ERROR(
655 "Incorrect number of parameter files. Only one parameter file is currently \
656 supported by this driver.");
657 return ier;
658 }
659
660 /* Set paramfile1name */
661 ier = KIM_ModelDriverCreate_GetParameterFileName(
662 modelDriverCreate, 0, ¶mfile1name);
663 if (ier == TRUE)
664 {
665 LOG_ERROR("Unable to get parameter file name");
666 return ier;
667 }
668
669 /* Read in Model parameters from parameter file */
670 fid = fopen(paramfile1name, "r");
671 if (fid == NULL)
672 {
673 ier = TRUE;
674 LOG_ERROR("Unable to open parameter file for Morse parameters");
675 return ier;
676 }
677
678 ier = fscanf(fid,
679 "%s \n%lf \n%lf \n%lf \n%lf",
680 speciesNameString,
681 &a, /* cutoff distance in angstroms */
682 &epsilon, /* Morse epsilon in eV */
683 &C, /* Morse C in 1/Angstroms */
684 &Rzero); /* Morse Rzero in Angstroms */
685 fclose(fid);
686
687 /* check that we read the right number of parameters */
688 if (5 != ier)
689 {
690 ier = TRUE;
691 LOG_ERROR("Unable to read all Morse parameters");
692 return ier;
693 }
694 /* Define default base units */
695 fromLength = KIM_LENGTH_UNIT_A;
696 fromEnergy = KIM_ENERGY_UNIT_eV;
697 fromCharge = KIM_CHARGE_UNIT_e;
698 fromTemperature = KIM_TEMPERATURE_UNIT_K;
699 fromTime = KIM_TIME_UNIT_ps;
700
701 /* Convert units of Model parameters */
702 ier = KIM_ModelDriverCreate_ConvertUnit(fromLength,
703 fromEnergy,
704 fromCharge,
705 fromTemperature,
706 fromTime,
707 requestedLengthUnit,
708 requestedEnergyUnit,
709 requestedChargeUnit,
710 requestedTemperatureUnit,
711 requestedTimeUnit,
712 1.0,
713 0.0,
714 0.0,
715 0.0,
716 0.0,
717 &convertLength);
718 if (ier)
719 {
720 LOG_ERROR("Unable to convert length unit");
721 return ier;
722 }
723
724 ier = KIM_ModelDriverCreate_ConvertUnit(fromLength,
725 fromEnergy,
726 fromCharge,
727 fromTemperature,
728 fromTime,
729 requestedLengthUnit,
730 requestedEnergyUnit,
731 requestedChargeUnit,
732 requestedTemperatureUnit,
733 requestedTimeUnit,
734 0.0,
735 1.0,
736 0.0,
737 0.0,
738 0.0,
739 &convertEnergy);
740 if (ier)
741 {
742 LOG_ERROR("Unable to convert energy unit");
743 return ier;
744 }
745
746 if (convertLength != ONE)
747 {
748 Rzero *= convertLength;
749 a *= convertLength;
750 C *= (1.0 / convertLength);
751 }
752
753 if (convertEnergy != ONE) { epsilon *= convertEnergy; }
754
755 /* Set cutoff equal to 'a' and work with it from now on */
756 cutoff = a;
757
758 /* Set units */
759 LOG_INFORMATION("Setting units");
760 ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
761 requestedLengthUnit,
762 requestedEnergyUnit,
763 KIM_CHARGE_UNIT_unused,
764 KIM_TEMPERATURE_UNIT_unused,
765 KIM_TIME_UNIT_unused);
766 if (ier == TRUE) { LOG_ERROR("Unable to set units to requested set"); }
767
768 /* Indicate that we use zero-based indexing for particles in this driver */
769 LOG_INFORMATION("Setting particle indexing to zero-based");
770 ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate,
771 KIM_NUMBERING_zeroBased);
772 if (ier == TRUE)
773 {
774 LOG_ERROR("Unable to set particle indexing");
775 return ier;
776 }
777
778 /* Store pointers to functions in the KIM object */
779 LOG_INFORMATION("Registering Model function pointers");
780 ier = ier
781 || KIM_ModelDriverCreate_SetRoutinePointer(
782 modelDriverCreate,
783 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsCreate,
784 KIM_LANGUAGE_NAME_c,
785 TRUE,
786 (KIM_Function *) &compute_arguments_create);
787 ier = ier
788 || KIM_ModelDriverCreate_SetRoutinePointer(
789 modelDriverCreate,
790 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsDestroy,
791 KIM_LANGUAGE_NAME_c,
792 TRUE,
793 (KIM_Function *) &compute_arguments_destroy);
794 ier = ier
795 || KIM_ModelDriverCreate_SetRoutinePointer(
796 modelDriverCreate,
797 KIM_MODEL_ROUTINE_NAME_Destroy,
798 KIM_LANGUAGE_NAME_c,
799 TRUE,
800 (KIM_Function *) &destroy);
801 ier = ier
802 || KIM_ModelDriverCreate_SetRoutinePointer(
803 modelDriverCreate,
804 KIM_MODEL_ROUTINE_NAME_Compute,
805 KIM_LANGUAGE_NAME_c,
806 TRUE,
807 (KIM_Function *) &compute);
808 ier = ier
809 || KIM_ModelDriverCreate_SetRoutinePointer(
810 modelDriverCreate,
811 KIM_MODEL_ROUTINE_NAME_Refresh,
812 KIM_LANGUAGE_NAME_c,
813 TRUE,
814 (KIM_Function *) &refresh);
815 if (ier == TRUE)
816 {
817 LOG_ERROR("Unable to register Model function pointers");
818 return ier;
819 }
820
821 /* Register species */
822 LOG_INFORMATION("Setting species code");
823 speciesName = KIM_SpeciesName_FromString(speciesNameString);
824 ier = ier
825 || KIM_ModelDriverCreate_SetSpeciesCode(
826 modelDriverCreate, speciesName, SPECCODE);
827 if (ier == TRUE)
828 {
829 LOG_ERROR("Unable to set species code");
830 return ier;
831 }
832
833 /* Allocate buffer */
834 LOG_INFORMATION("Allocating memory for Model buffer");
835 buffer = (struct model_buffer *) malloc(sizeof(struct model_buffer));
836 if (buffer == NULL)
837 {
838 ier = TRUE;
839 LOG_ERROR("Could not allocate Model buffer");
840 return ier;
841 }
842
843 /* Store model buffer in KIM object */
844 LOG_INFORMATION("Registering Model buffer");
845 KIM_ModelDriverCreate_SetModelBufferPointer(modelDriverCreate,
846 (void *) buffer);
847
848 /* Store parameter values in buffer */
849 LOG_INFORMATION("Loading Model parameters into buffer");
850 buffer->influenceDistance = cutoff;
851 buffer->cutsq = cutoff * cutoff;
852 buffer->epsilon = epsilon;
853 buffer->Rzero = Rzero;
854 buffer->C = C;
855 buffer->Pcutoff = cutoff;
856 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles = 1;
857
858 /* set value of parameter shift */
859 dummy = 0.0;
860 /* call calc_phi with r=cutoff and shift=0.0 */
861 calc_phi(&(buffer->epsilon),
862 &(buffer->C),
863 &(buffer->Rzero),
864 &dummy,
865 &cutoff,
866 cutoff,
867 &(buffer->shift));
868 /* set shift to -shift */
869 buffer->shift = -buffer->shift;
870
871 /* Register influence distance pointer */
872 LOG_INFORMATION("Registering influence distance pointer");
873 KIM_ModelDriverCreate_SetInfluenceDistancePointer(
874 modelDriverCreate, &(buffer->influenceDistance));
875
876 /* Register cutoff pointer */
877 LOG_INFORMATION("Registering cutoff pointer");
878 KIM_ModelDriverCreate_SetNeighborListPointers(
879 modelDriverCreate,
880 1,
881 &(buffer->Pcutoff),
882 (const int *) &(
883 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
884
885
886 /* Register Model parameters */
887 LOG_INFORMATION("Registering Model parameters");
888 KIM_ModelDriverCreate_SetParameterPointerDouble(
889 modelDriverCreate,
890 1,
891 &(buffer->Pcutoff),
892 "cutoff",
893 "Distance beyond which particles do not "
894 "interact with one another");
895 KIM_ModelDriverCreate_SetParameterPointerDouble(modelDriverCreate,
896 1,
897 &(buffer->Rzero),
898 "Rzero",
899 "Equilibrium bond distance.");
900 KIM_ModelDriverCreate_SetParameterPointerDouble(
901 modelDriverCreate,
902 1,
903 &(buffer->C),
904 "C",
905 "Multiplying factor in the arguments of the exponential terms.");
906 KIM_ModelDriverCreate_SetParameterPointerDouble(
907 modelDriverCreate,
908 1,
909 &(buffer->epsilon),
910 "epsilon",
911 "Multiplying factor of the traditional "
912 "Morse potential, i.e. not including the added shift. This is also the "
913 "maximum depth of the potential well if no shift were added.");
914
915
916 if (ier == TRUE)
917 {
918 free(buffer);
919 LOG_ERROR("Unable to successfully initialize Model");
920 }
921
922 return ier;
923 }
924
925 /* Refresh function */
926 #undef KIM_LOGGER_FUNCTION_NAME
927 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelRefresh_LogEntry
928 #undef KIM_LOGGER_OBJECT_NAME
929 #define KIM_LOGGER_OBJECT_NAME modelRefresh
refresh(KIM_ModelRefresh * const modelRefresh)930 static int refresh(KIM_ModelRefresh * const modelRefresh)
931 {
932 /* Local variables */
933 double cutoff;
934 double dummy;
935 struct model_buffer * buffer;
936
937 /* Get model buffer from KIM object */
938 LOG_INFORMATION("Retrieving Model buffer");
939 KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
940
941 LOG_INFORMATION("Resetting influence distance and cutoff, and shift");
942 cutoff = buffer->Pcutoff;
943
944 /* set value of parameter cutsq */
945 buffer->cutsq = cutoff*cutoff;
946
947 /* set value of influence distance */
948 buffer->influenceDistance = cutoff;
949
950 /* set value of parameter shift */
951 dummy = 0.0;
952 /* call calc_phi with r=cutoff and shift=0.0 */
953 calc_phi(&(buffer->epsilon),
954 &(buffer->C),
955 &(buffer->Rzero),
956 &dummy,
957 &cutoff,
958 cutoff,
959 &(buffer->shift));
960 /* set shift to -shift */
961 buffer->shift = -buffer->shift;
962 KIM_ModelRefresh_SetInfluenceDistancePointer(modelRefresh,
963 &(buffer->influenceDistance));
964
965 KIM_ModelRefresh_SetNeighborListPointers(
966 modelRefresh,
967 1,
968 &(buffer->Pcutoff),
969 (const int *) &(
970 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
971
972 return FALSE;
973 }
974
975 /* Destroy function */
976 #undef KIM_LOGGER_FUNCTION_NAME
977 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDestroy_LogEntry
978 #undef KIM_LOGGER_OBJECT_NAME
979 #define KIM_LOGGER_OBJECT_NAME modelDestroy
destroy(KIM_ModelDestroy * const modelDestroy)980 static int destroy(KIM_ModelDestroy * const modelDestroy)
981 {
982 /* Local variables */
983 struct model_buffer * buffer;
984
985 /* Get model buffer from KIM object */
986 KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
987
988 /* Free the buffer */
989 LOG_INFORMATION("Deallocating Model buffer");
990 free(buffer);
991
992 return FALSE;
993 }
994
995 /* compute_arguments create routine */
996 #undef KIM_LOGGER_FUNCTION_NAME
997 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelComputeArgumentsCreate_LogEntry
998 #undef KIM_LOGGER_OBJECT_NAME
999 #define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate
compute_arguments_create(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)1000 static int compute_arguments_create(
1001 KIM_ModelCompute const * const modelCompute,
1002 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
1003 {
1004 int ier;
1005
1006 (void) modelCompute; /* avoid unused argument warning */
1007
1008 ier = FALSE;
1009
1010 /* Register arguments */
1011 ier = ier
1012 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
1013 modelComputeArgumentsCreate,
1014 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
1015 KIM_SUPPORT_STATUS_optional);
1016 ier = ier
1017 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
1018 modelComputeArgumentsCreate,
1019 KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
1020 KIM_SUPPORT_STATUS_optional);
1021 ier = ier
1022 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
1023 modelComputeArgumentsCreate,
1024 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
1025 KIM_SUPPORT_STATUS_optional);
1026
1027 ier = ier
1028 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
1029 modelComputeArgumentsCreate,
1030 KIM_COMPUTE_ARGUMENT_NAME_partialVirial,
1031 KIM_SUPPORT_STATUS_optional);
1032
1033 ier = ier
1034 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
1035 modelComputeArgumentsCreate,
1036 KIM_COMPUTE_ARGUMENT_NAME_partialParticleVirial,
1037 KIM_SUPPORT_STATUS_optional);
1038
1039 /* register call backs */
1040 LOG_INFORMATION("Register call back supportStatus");
1041 ier = ier
1042 || KIM_ModelComputeArgumentsCreate_SetCallbackSupportStatus(
1043 modelComputeArgumentsCreate,
1044 KIM_COMPUTE_CALLBACK_NAME_ProcessDEDrTerm,
1045 KIM_SUPPORT_STATUS_optional);
1046 ier = ier
1047 || KIM_ModelComputeArgumentsCreate_SetCallbackSupportStatus(
1048 modelComputeArgumentsCreate,
1049 KIM_COMPUTE_CALLBACK_NAME_ProcessD2EDr2Term,
1050 KIM_SUPPORT_STATUS_optional);
1051
1052 if (ier == TRUE) { LOG_ERROR("Unable to set argument supportStatus."); }
1053 return ier;
1054 }
1055
1056 /* compute arguments destroy routine */
compute_arguments_destroy(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)1057 static int compute_arguments_destroy(
1058 KIM_ModelCompute const * const modelCompute,
1059 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
1060 {
1061 (void) modelCompute; /* avoid unused parameter warning */
1062 (void) modelComputeArgumentsDestroy;
1063
1064 /* Nothing further to do */
1065
1066 return FALSE;
1067 }
1068