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--2016, Regents of the University of Minnesota.
24 * All rights reserved.
25 *
26 * Contributors:
27 * Ryan S. Elliott
28 * Andrew Akerson
29 * Ellad B. Tadmor
30 * Valeriu Smirichinski
31 * Stephen M. Whalen
32 *
33 */
34
35
36 /*******************************************************************************
37 *
38 * MorseEIP_GuthikondaElliott_2009
39 *
40 * Temperature-dependent Morse pair potential KIM Model Driver
41 * Shifted to have zero energy at the cutoff radius
42 *
43 * Language: C
44 *
45 ******************************************************************************/
46
47 #include "KIM_ModelDriverHeaders.h"
48 #include <math.h>
49 #include <stdarg.h>
50 #include <stdio.h>
51 #include <stdlib.h>
52 #include <string.h>
53
54 #include "KIM_LogMacros.h"
55
56 #define DIM 3 /* dimensionality of space */
57 #define MAXLINE 1024 /* max characters in line */
58 #define TRUE 1
59 #define FALSE 0
60 #define ONE 1.0
61
62 /* Define prototypes for Model Driver create and unit conversion function */
63 /**/
64 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
65 KIM_LengthUnit const requestedLengthUnit,
66 KIM_EnergyUnit const requestedEnergyUnit,
67 KIM_ChargeUnit const requestedChargeUnit,
68 KIM_TemperatureUnit const requestedTemperatureUnit,
69 KIM_TimeUnit const requestedTimeUnit);
70
71 int ConvertUnits(KIM_ModelDriverCreate * const modelDriverCreate,
72 KIM_LengthUnit const requestedLengthUnit,
73 KIM_EnergyUnit const requestedEnergyUnit,
74 KIM_ChargeUnit const requestedChargeUnit,
75 KIM_TemperatureUnit const requestedTemperatureUnit,
76 KIM_TimeUnit const requestedTimeUnit,
77 int numberUniqueSpeciesPairs,
78 double * const cutoffs,
79 double * const A1s,
80 double * const A2s,
81 double * const r1s,
82 double * const r2s);
83
84 /* Define prototypes for Model (Driver) refresh and destroy */
85 /* defined as static to avoid namespace clashes with other Models */
86 /**/
87 static int refresh(KIM_ModelRefresh * const modelRefresh);
88 static int destroy(KIM_ModelDestroy * const modelDestroy);
89
90 /* Define prototypes for compute, compute_arguments_create, and
91 * compute_arguments_destroy */
92 /**/
93 static int
94 compute(KIM_ModelCompute const * const modelCompute,
95 KIM_ModelComputeArguments const * const modelComputeArguments);
96 static int compute_arguments_create(
97 KIM_ModelCompute const * const modelCompute,
98 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
99 static int compute_arguments_destroy(
100 KIM_ModelCompute const * const modelCompute,
101 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
102
103 /* Define prototypes for internal functions used in this driver */
104 static void calc_phi(double const * const A,
105 double const * const B,
106 double const * const rHat,
107 double const * const shift,
108 double const * const cutoff,
109 double const r,
110 double * const phi);
111
112 static void calc_phi_dphi(double const * const A,
113 double const * const B,
114 double const * const rHat,
115 double const * const shift,
116 double const * const cutoff,
117 double const r,
118 double * const phi,
119 double * const dphi);
120
121 static void calc_phi_d2phi(double const * const A,
122 double const * const B,
123 double const * const rHat,
124 double const * const shift,
125 double const * const cutoff,
126 double const r,
127 double * const phi,
128 double * const dphi,
129 double * const d2phi);
130
131 double calc_A(double const * const A1,
132 double const * const A2,
133 double const * const A3,
134 double const theta);
135
136 double calc_B(double const * const B1,
137 double const * const B2,
138 double const * const B3,
139 double const theta);
140
141 double calc_rHat(double const * const r1,
142 double const * const r2,
143 double const * const r3,
144 double const theta);
145
146 double ** AllocateAndInitialize2DArray(int const extentZero,
147 int const extentOne);
148
149 double * AllocateAndInitialize1DArray(int const numberModelSpecies);
150 void Deallocate2DArrays(int const n, ...);
151 void Deallocate1DArrays(int const n, ...);
152 void getNextDataLine(FILE * const filePtr,
153 char * nextLinePtr,
154 int const maxSize,
155 int * const endOfFileFlag);
156
157 /* Define model_buffer structure */
158 struct model_buffer
159 {
160 int modelWillNotRequestNeighborsOfNoncontributingParticles;
161
162 double influenceDistance;
163 double * cutoffs;
164
165 int numberModelSpecies;
166
167 double temperature;
168
169 double * A1s;
170 double * A2s;
171 double * A3s;
172 double * B1s;
173 double * B2s;
174 double * B3s;
175 double * r1s;
176 double * r2s;
177 double * r3s;
178
179 double ** cutsq2D;
180 double ** As2D;
181 double ** Bs2D;
182 double ** rHats2D;
183 double ** shifts2D;
184 };
185
186 /* Calculate Guthikonda Elliott Paramaters */
calc_A(double const * const A1,double const * const A2,double const * const A3,double const theta)187 double calc_A(double const * const A1,
188 double const * const A2,
189 double const * const A3,
190 double const theta)
191 {
192 double A;
193
194 A = *A1 + (*A2) * (pow(theta, *A3) - 1.0);
195 return A;
196 }
197
calc_B(double const * const B1,double const * const B2,double const * const B3,double const theta)198 double calc_B(double const * const B1,
199 double const * const B2,
200 double const * const B3,
201 double const theta)
202 {
203 double B;
204
205 B = *B1 + (*B2) * (pow(theta, *B3) - 1.0);
206 return B;
207 }
208
calc_rHat(double const * const r1,double const * const r2,double const * const r3,double const theta)209 double calc_rHat(double const * const r1,
210 double const * const r2,
211 double const * const r3,
212 double const theta)
213 {
214 double rHat;
215
216 rHat = *r1 + (*r2) * (exp((*r3) * (theta - 1.0)) - 1.0);
217 return rHat;
218 }
219
220 /* Calculate pair potential phi(r) */
calc_phi(double const * const A,double const * const B,double const * const rHat,double const * const shift,double const * const cutoff,double const r,double * const phi)221 static void calc_phi(double const * const A,
222 double const * const B,
223 double const * const rHat,
224 double const * const shift,
225 double const * const cutoff,
226 double const r,
227 double * const phi)
228 {
229 /* local variables */
230 double ep;
231 double ep2;
232 double epsilon;
233 double C;
234
235 epsilon = -(*A);
236 C = (*B) / (*rHat);
237
238 ep = exp(-(C) * (r - *rHat));
239 ep2 = ep * ep;
240
241 if (r > *cutoff)
242 {
243 /* Argument exceeds cutoff radius */
244 *phi = 0.0;
245 }
246 else
247 {
248 *phi = (epsilon) * (-ep2 + 2.0 * ep) + *shift;
249 }
250
251 return;
252 }
253
254 /* Calculate pair potential phi(r) and its derivative dphi(r) */
calc_phi_dphi(double const * const A,double const * const B,double const * const rHat,double const * const shift,double const * const cutoff,double const r,double * const phi,double * const dphi)255 static void calc_phi_dphi(double const * const A,
256 double const * const B,
257 double const * const rHat,
258 double const * const shift,
259 double const * const cutoff,
260 double const r,
261 double * const phi,
262 double * const dphi)
263 {
264 /* local variables */
265 double ep;
266 double ep2;
267 double epsilon;
268 double C;
269
270 epsilon = -(*A);
271 C = (*B) / (*rHat);
272
273 ep = exp(-(C) * (r - *rHat));
274 ep2 = ep * ep;
275
276 if (r > *cutoff)
277 {
278 /* Argument exceeds cutoff radius */
279 *phi = 0.0;
280 *dphi = 0.0;
281 }
282 else
283 {
284 *phi = (epsilon) * (-ep2 + 2.0 * ep) + *shift;
285 *dphi = 2.0 * (epsilon) * (C) * (-ep + ep2);
286 }
287
288 return;
289 }
290
291 /*
292 Calculate pair potential phi(r) and its 1st & 2nd derivatives dphi(r),
293 d2phi(r)
294 */
calc_phi_d2phi(double const * const A,double const * const B,double const * const rHat,double const * const shift,double const * const cutoff,double const r,double * const phi,double * const dphi,double * const d2phi)295 static void calc_phi_d2phi(double const * const A,
296 double const * const B,
297 double const * const rHat,
298 double const * const shift,
299 double const * const cutoff,
300 double const r,
301 double * const phi,
302 double * const dphi,
303 double * const d2phi)
304 {
305 /* local variables */
306 double ep;
307 double ep2;
308 double epsilon;
309 double C;
310
311 epsilon = -(*A);
312 C = (*B) / (*rHat);
313 ep = exp(-(C) * (r - *rHat));
314 ep2 = ep * ep;
315
316 if (r > *cutoff)
317 {
318 /* Argument exceeds cutoff radius */
319 *phi = 0.0;
320 *dphi = 0.0;
321 *d2phi = 0.0;
322 }
323 else
324 {
325 *phi = (epsilon) * (-ep2 + 2.0 * ep) + *shift;
326 *dphi = 2.0 * (epsilon) * (C) * (-ep + ep2);
327 *d2phi = 2.0 * (epsilon) * (C) * (C) * (ep - 2.0 * ep2);
328 }
329
330 return;
331 }
332
333 /*****************************************************************************/
334 /* Compute function */
335 /*****************************************************************************/
336 #undef KIM_LOGGER_FUNCTION_NAME
337 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelCompute_LogEntry
338 #undef KIM_LOGGER_OBJECT_NAME
339 #define KIM_LOGGER_OBJECT_NAME modelCompute
340 static int
compute(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArguments const * const modelComputeArguments)341 compute(KIM_ModelCompute const * const modelCompute,
342 KIM_ModelComputeArguments const * const modelComputeArguments)
343 {
344 /* local variables */
345 double R;
346 double R_pairs[2];
347 double * pR_pairs = &(R_pairs[0]);
348 double Rsqij;
349 double phi;
350 double dphi;
351 double d2phi;
352 double dEidr = 0.0;
353 double d2Eidr = 0.0;
354 double Rij[DIM];
355 double * pRij = &(Rij[0]);
356 double Rij_pairs[2][3];
357 double * pRij_pairs = &(Rij_pairs[0][0]);
358 int ier;
359 int i;
360 int i_pairs[2];
361 int * pi_pairs = &(i_pairs[0]);
362 int j;
363 int j_pairs[2];
364 int * pj_pairs = &(j_pairs[0]);
365 int jj;
366 int k;
367 int iSpecies, jSpecies;
368 int const * neighListOfCurrentAtom;
369 struct model_buffer * buffer;
370 int comp_energy;
371 int comp_force;
372 int comp_particleEnergy;
373 int comp_process_dEdr;
374 int comp_process_d2Edr2;
375
376 int * nAtoms;
377 int * particleSpeciesCodes;
378 int * particleContributing;
379 double ijcutoff;
380 double ** cutsq2D;
381 double ** As2D;
382 double ** Bs2D;
383 double ** rHats2D;
384 double ** shifts2D;
385 double * coords;
386 double * energy;
387 double * force;
388 double * particleEnergy;
389 int numOfAtomNeigh;
390
391 /* get buffer from KIM object */
392 KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
393
394 /* unpack info from the buffer */
395 cutsq2D = (buffer->cutsq2D);
396 As2D = (buffer->As2D);
397 Bs2D = (buffer->Bs2D);
398 rHats2D = (buffer->rHats2D);
399 shifts2D = (buffer->shifts2D);
400
401 /*
402 check to see if we have been asked to compute the forces, particleEnergy,
403 process_dEdr, and process_d2Edr2
404 */
405 ier = KIM_ModelComputeArguments_GetArgumentPointerInteger(
406 modelComputeArguments,
407 KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles,
408 &nAtoms)
409 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
410 modelComputeArguments,
411 KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes,
412 &particleSpeciesCodes)
413 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
414 modelComputeArguments,
415 KIM_COMPUTE_ARGUMENT_NAME_particleContributing,
416 &particleContributing)
417 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
418 modelComputeArguments,
419 KIM_COMPUTE_ARGUMENT_NAME_coordinates,
420 &coords)
421 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
422 modelComputeArguments,
423 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
424 &energy)
425 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
426 modelComputeArguments,
427 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
428 &force)
429 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
430 modelComputeArguments,
431 KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
432 &particleEnergy);
433 if (ier == TRUE)
434 {
435 LOG_ERROR("GetArgumentPointer");
436 return ier;
437 }
438
439 comp_energy = (energy != NULL);
440 comp_force = (force != NULL);
441 comp_particleEnergy = (particleEnergy != NULL);
442
443 ier = KIM_ModelComputeArguments_IsCallbackPresent(
444 modelComputeArguments,
445 KIM_COMPUTE_CALLBACK_NAME_ProcessDEDrTerm,
446 &comp_process_dEdr)
447 || KIM_ModelComputeArguments_IsCallbackPresent(
448 modelComputeArguments,
449 KIM_COMPUTE_CALLBACK_NAME_ProcessD2EDr2Term,
450 &comp_process_d2Edr2);
451 if (ier == TRUE)
452 {
453 LOG_ERROR("IsCallbackPresent");
454 return ier;
455 }
456
457 /* initialize potential energies and forces */
458 if (comp_particleEnergy)
459 {
460 for (i = 0; i < *nAtoms; ++i) { particleEnergy[i] = 0.0; }
461 }
462 if (comp_energy) { *energy = 0.0; }
463
464 if (comp_force)
465 {
466 for (i = 0; i < *nAtoms; ++i)
467 {
468 for (k = 0; k < DIM; ++k) { force[i * DIM + k] = 0.0; }
469 }
470 }
471
472 /* Compute energy and forces */
473
474 /* loop over particles and compute energy and forces */
475 for (i = 0; i < *nAtoms; i++)
476 {
477 if (particleContributing[i])
478 {
479 ier = KIM_ModelComputeArguments_GetNeighborList(modelComputeArguments,
480 0,
481 i,
482 &numOfAtomNeigh,
483 &neighListOfCurrentAtom);
484 if (ier)
485 {
486 /* some sort of problem, exit */
487 LOG_ERROR("KIM_get_neigh");
488 ier = TRUE;
489 return ier;
490 }
491
492 iSpecies = particleSpeciesCodes[i];
493
494 /* loop over the neighbors of atom i */
495 for (jj = 0; jj < numOfAtomNeigh; ++jj)
496 {
497 /* get neighbor ID */
498 j = neighListOfCurrentAtom[jj];
499
500 if (!(particleContributing[j] && j < i))
501 {
502 jSpecies = particleSpeciesCodes[j];
503 /* compute relative position vector and squared distance */
504 Rsqij = 0.0;
505 for (k = 0; k < DIM; ++k)
506 {
507 Rij[k] = coords[j * DIM + k] - coords[i * DIM + k];
508
509 /* compute squared distance */
510 Rsqij += Rij[k] * Rij[k];
511 }
512
513 /* compute energy and force */
514 if (Rsqij
515 < cutsq2D[iSpecies][jSpecies]) /* particles are interacting? */
516 {
517 ijcutoff = sqrt(cutsq2D[iSpecies][jSpecies]);
518 R = sqrt(Rsqij);
519 if (comp_process_d2Edr2)
520 {
521 /* compute pair potential and its derivatives */
522 calc_phi_d2phi(&As2D[iSpecies][jSpecies],
523 &Bs2D[iSpecies][jSpecies],
524 &rHats2D[iSpecies][jSpecies],
525 &shifts2D[iSpecies][jSpecies],
526 &ijcutoff,
527 R,
528 &phi,
529 &dphi,
530 &d2phi);
531
532 /* compute dEidr and d2Eidr */
533 dEidr = 0.5 * dphi;
534 d2Eidr = 0.5 * d2phi;
535 }
536 else if (comp_force || comp_process_dEdr)
537 {
538 /* compute pair potential and its derivative */
539 calc_phi_dphi(&(As2D[iSpecies][jSpecies]),
540 &(Bs2D[iSpecies][jSpecies]),
541 &(rHats2D[iSpecies][jSpecies]),
542 &(shifts2D[iSpecies][jSpecies]),
543 &ijcutoff,
544 R,
545 &phi,
546 &dphi);
547
548 /* compute dEidr */
549 dEidr = 0.5 * dphi;
550 }
551 else
552 {
553 /* compute just pair potential */
554 calc_phi(&(As2D[iSpecies][jSpecies]),
555 &(Bs2D[iSpecies][jSpecies]),
556 &(rHats2D[iSpecies][jSpecies]),
557 &(shifts2D[iSpecies][jSpecies]),
558 &ijcutoff,
559 R,
560 &phi);
561 }
562
563 /* contribution to energy */
564 if (comp_particleEnergy) { particleEnergy[i] += 0.5 * phi; }
565 if (comp_energy)
566 {
567 /* Full mode -- add half v to total energy */
568 *energy += 0.5 * phi;
569 }
570
571 /* contribution to process_dEdr */
572 if (comp_process_dEdr)
573 {
574 ier = KIM_ModelComputeArguments_ProcessDEDrTerm(
575 modelComputeArguments, dEidr, R, pRij, i, j);
576 }
577
578 /* contribution to process_d2Edr2 */
579 if (comp_process_d2Edr2)
580 {
581 R_pairs[0] = R_pairs[1] = R;
582 Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
583 Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
584 Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
585 i_pairs[0] = i_pairs[1] = i;
586 j_pairs[0] = j_pairs[1] = j;
587
588 ier = KIM_ModelComputeArguments_ProcessD2EDr2Term(
589 modelComputeArguments,
590 d2Eidr,
591 pR_pairs,
592 pRij_pairs,
593 pi_pairs,
594 pj_pairs);
595 }
596
597 /* contribution to forces */
598 if (comp_force)
599 {
600 for (k = 0; k < DIM; ++k)
601 { /* accumulate force on atom i */
602 force[i * DIM + k] += dEidr * Rij[k] / R;
603 /* accumulate force on atom j */
604 force[j * DIM + k] -= dEidr * Rij[k] / R;
605 }
606 }
607 }
608 } /* End effective half-list check (i < j) */
609 } /* loop on jj */
610 } /* Check on whether particle i is contributing */
611 } /* Outer loop over all atoms */
612
613 /* everything is great */
614 ier = FALSE;
615 return ier;
616 }
617
618
619 /*****************************************************************************/
620 /* Miscellaneous helper functions */
621 /*****************************************************************************/
getNextDataLine(FILE * const filePtr,char * nextLinePtr,int const maxSize,int * const endOfFileFlag)622 void getNextDataLine(FILE * const filePtr,
623 char * nextLinePtr,
624 int const maxSize,
625 int * const endOfFileFlag)
626 {
627 *endOfFileFlag = 0;
628 do
629 {
630 if (fgets(nextLinePtr, maxSize, filePtr) == NULL)
631 {
632 *endOfFileFlag = 1;
633 break;
634 }
635 while ((nextLinePtr[0] == ' ' || nextLinePtr[0] == '\t')
636 || (nextLinePtr[0] == '\n' || nextLinePtr[0] == '\r'))
637 { nextLinePtr = (nextLinePtr + 1); }
638 } while ((strncmp("#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
639 }
640
641 /*****************************************************************************/
AllocateAndInitialize2DArray(int const extentZero,int const extentOne)642 double ** AllocateAndInitialize2DArray(int const extentZero,
643 int const extentOne)
644 {
645 double ** arrayPtr;
646 int i, j;
647 /* allocate memory and set pointers */
648 arrayPtr = malloc(extentZero * sizeof(double *));
649 arrayPtr[0] = malloc((extentZero * extentOne) * sizeof(double));
650 for (i = 1; i < extentZero; ++i)
651 { arrayPtr[i] = arrayPtr[i - 1] + extentOne; }
652
653 /* initialize */
654 for (i = 0; i < extentZero; ++i)
655 {
656 for (j = 0; j < extentOne; ++j) { arrayPtr[i][j] = 0.0; }
657 }
658 return arrayPtr;
659 }
660
661 /*****************************************************************************/
AllocateAndInitialize1DArray(int const numberModelSpecies)662 double * AllocateAndInitialize1DArray(int const numberModelSpecies)
663 {
664 double * arrayPtr;
665 int numberUniqueSpeciesPairs;
666
667 /* allocate memory and set pointers */
668 numberUniqueSpeciesPairs
669 = ((numberModelSpecies + 1) * numberModelSpecies) / 2;
670 arrayPtr = calloc(numberUniqueSpeciesPairs, sizeof(double));
671 return arrayPtr;
672 }
673
674 /*****************************************************************************/
Deallocate1DArrays(int const n,...)675 void Deallocate1DArrays(int const n, ...)
676 {
677 int i;
678 va_list pointerArgs;
679
680 va_start(pointerArgs, n);
681 for (i = 0; i < n; i++) { free(va_arg(pointerArgs, double *)); }
682 va_end(pointerArgs);
683 }
684
685 /*****************************************************************************/
Deallocate2DArrays(int const n,...)686 void Deallocate2DArrays(int const n, ...)
687 {
688 double ** nextArray;
689 int i;
690 va_list doublePointerArgs;
691
692 va_start(doublePointerArgs, n);
693 for (i = 0; i < n; i++)
694 {
695 nextArray = va_arg(doublePointerArgs, double **);
696 free(nextArray[0]);
697 free(nextArray);
698 }
699 va_end(doublePointerArgs);
700 }
701
702 /*****************************************************************************/
703 /* Unit conversion function */
704 /*****************************************************************************/
705 #undef KIM_LOGGER_FUNCTION_NAME
706 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDriverCreate_LogEntry
707 #undef KIM_LOGGER_OBJECT_NAME
708 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
ConvertUnits(KIM_ModelDriverCreate * const modelDriverCreate,KIM_LengthUnit const requestedLengthUnit,KIM_EnergyUnit const requestedEnergyUnit,KIM_ChargeUnit const requestedChargeUnit,KIM_TemperatureUnit const requestedTemperatureUnit,KIM_TimeUnit const requestedTimeUnit,int numberUniqueSpeciesPairs,double * const cutoffs,double * const A1s,double * const A2s,double * const r1s,double * const r2s)709 int ConvertUnits(KIM_ModelDriverCreate * const modelDriverCreate,
710 KIM_LengthUnit const requestedLengthUnit,
711 KIM_EnergyUnit const requestedEnergyUnit,
712 KIM_ChargeUnit const requestedChargeUnit,
713 KIM_TemperatureUnit const requestedTemperatureUnit,
714 KIM_TimeUnit const requestedTimeUnit,
715 int numberUniqueSpeciesPairs,
716 double * const cutoffs,
717 double * const A1s,
718 double * const A2s,
719 double * const r1s,
720 double * const r2s)
721 {
722 int ier;
723 int i;
724
725 double convertLength = 1.0;
726 double convertEnergy = 1.0;
727
728 /* define default base units */
729 KIM_LengthUnit const fromLength = KIM_LENGTH_UNIT_A;
730 KIM_EnergyUnit const fromEnergy = KIM_ENERGY_UNIT_eV;
731 KIM_ChargeUnit const fromCharge = KIM_CHARGE_UNIT_e;
732 KIM_TemperatureUnit const fromTemperature = KIM_TEMPERATURE_UNIT_K;
733 KIM_TimeUnit const fromTime = KIM_TIME_UNIT_ps;
734
735 ier = KIM_ModelDriverCreate_ConvertUnit(fromLength,
736 fromEnergy,
737 fromCharge,
738 fromTemperature,
739 fromTime,
740 requestedLengthUnit,
741 requestedEnergyUnit,
742 requestedChargeUnit,
743 requestedTemperatureUnit,
744 requestedTimeUnit,
745 1.0,
746 0.0,
747 0.0,
748 0.0,
749 0.0,
750 &convertLength);
751 if (ier)
752 {
753 LOG_ERROR("Unable to convert length unit");
754 return ier;
755 }
756
757 if (convertLength != ONE)
758 {
759 for (i = 0; i < numberUniqueSpeciesPairs; ++i)
760 {
761 cutoffs[i] *= convertLength;
762 r1s[i] *= convertLength;
763 r2s[i] *= convertLength;
764 }
765 }
766
767 /* Energy */
768 ier = KIM_ModelDriverCreate_ConvertUnit(fromLength,
769 fromEnergy,
770 fromCharge,
771 fromTemperature,
772 fromTime,
773 requestedLengthUnit,
774 requestedEnergyUnit,
775 requestedChargeUnit,
776 requestedTemperatureUnit,
777 requestedTimeUnit,
778 0.0,
779 1.0,
780 0.0,
781 0.0,
782 0.0,
783 &convertEnergy);
784 if (ier)
785 {
786 LOG_ERROR("Unable to convert energy unit");
787 return ier;
788 }
789
790 if (convertEnergy != ONE)
791 {
792 for (i = 0; i < numberUniqueSpeciesPairs; ++i)
793 {
794 A1s[i] *= convertEnergy;
795 A2s[i] *= convertEnergy;
796 }
797 }
798
799 /* register units */
800 ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
801 requestedLengthUnit,
802 requestedEnergyUnit,
803 KIM_CHARGE_UNIT_unused,
804 KIM_TEMPERATURE_UNIT_unused,
805 requestedTimeUnit);
806 if (ier)
807 {
808 LOG_ERROR("Unable to set units to requested values");
809 return ier;
810 }
811
812 ier = FALSE;
813 return ier;
814 }
815
816
817 /*****************************************************************************/
818 /* Refresh function */
819 /*****************************************************************************/
refresh(KIM_ModelRefresh * const modelRefresh)820 static int refresh(KIM_ModelRefresh * const modelRefresh)
821 {
822 /* Local variables */
823 int ier;
824 double influenceDistance;
825 double numberModelSpecies;
826 double * cutoffs;
827 double *A1s, *A2s, *A3s;
828 double *B1s, *B2s, *B3s;
829 double *r1s, *r2s, *r3s;
830 double temperature;
831 double theta;
832
833 double nextShift;
834 double dummy;
835 int indx;
836 int i, j;
837
838 struct model_buffer * buffer;
839
840 /* get buffer from KIM object */
841 KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
842
843 /* set value for 2D parameters */
844 numberModelSpecies = buffer->numberModelSpecies;
845 cutoffs = buffer->cutoffs;
846 temperature = buffer->temperature;
847 A1s = buffer->A1s;
848 A2s = buffer->A2s;
849 A3s = buffer->A3s;
850 B1s = buffer->B1s;
851 B2s = buffer->B2s;
852 B3s = buffer->B3s;
853 r1s = buffer->r1s;
854 r2s = buffer->r2s;
855 r3s = buffer->r3s;
856
857 theta = temperature / 333.15;
858 influenceDistance = 0.0;
859 for (i = 0; i < numberModelSpecies; ++i)
860 {
861 for (j = i; j < numberModelSpecies; ++j)
862 {
863 indx = i * numberModelSpecies + j - (i * i + i) / 2;
864
865 if (influenceDistance < cutoffs[indx])
866 { influenceDistance = cutoffs[indx]; }
867
868 buffer->cutsq2D[i][j] = buffer->cutsq2D[j][i]
869 = (cutoffs[indx] * cutoffs[indx]);
870
871 buffer->As2D[i][j] = buffer->As2D[j][i]
872 = calc_A(&A1s[indx], &A2s[indx], &A3s[indx], theta);
873
874 buffer->Bs2D[i][j] = buffer->Bs2D[j][i]
875 = calc_B(&B1s[indx], &B2s[indx], &B3s[indx], theta);
876
877 buffer->rHats2D[i][j] = buffer->rHats2D[j][i]
878 = calc_rHat(&r1s[indx], &r2s[indx], &r3s[indx], theta);
879 }
880 }
881 buffer->influenceDistance = influenceDistance;
882
883 /* Set influence distance pointer and cutoffs pointer */
884 KIM_ModelRefresh_SetInfluenceDistancePointer(modelRefresh,
885 &(buffer->influenceDistance));
886 KIM_ModelRefresh_SetNeighborListPointers(
887 modelRefresh,
888 1, /* Use only a single neighbor list */
889 &(buffer->influenceDistance),
890 (const int *) &(
891 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
892
893 /* Set Values for Shifts */
894 dummy = 0.0;
895 influenceDistance += 1.0; /* add a bit to avoid rounding problem */
896 for (i = 0; i < numberModelSpecies; i++)
897 {
898 for (j = i; j < numberModelSpecies; j++)
899 {
900 indx = i * numberModelSpecies + j - (i * i + i) / 2;
901 /* call calc_phi with r=cutoff and shift=0.0 */
902 calc_phi(&(buffer->As2D[i][j]),
903 &(buffer->Bs2D[i][j]),
904 &(buffer->rHats2D[i][j]),
905 &dummy,
906 &influenceDistance,
907 cutoffs[indx],
908 &nextShift);
909 /* set shift to -shift */
910 buffer->shifts2D[i][j] = buffer->shifts2D[j][i] = -nextShift;
911 }
912 }
913
914 ier = FALSE;
915 return ier;
916 }
917
918 /*****************************************************************************/
919 /* Destroy function */
920 /*****************************************************************************/
destroy(KIM_ModelDestroy * const modelDestroy)921 static int destroy(KIM_ModelDestroy * const modelDestroy)
922 {
923 /* Local variables */
924 int ier;
925 struct model_buffer * buffer;
926
927 KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
928
929 /* destroy the buffer */
930 Deallocate2DArrays(5,
931 buffer->cutsq2D,
932 buffer->As2D,
933 buffer->Bs2D,
934 buffer->rHats2D,
935 buffer->shifts2D);
936 Deallocate1DArrays(10,
937 buffer->cutoffs,
938 buffer->A1s,
939 buffer->A2s,
940 buffer->A3s,
941 buffer->B1s,
942 buffer->B2s,
943 buffer->B3s,
944 buffer->r1s,
945 buffer->r2s,
946 buffer->r3s);
947 free(buffer);
948
949 ier = FALSE;
950 return ier;
951 }
952
953 /*****************************************************************************/
954 /* Compute arguments create function */
955 /*****************************************************************************/
956 #undef KIM_LOGGER_FUNCTION_NAME
957 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelComputeArgumentsCreate_LogEntry
958 #undef KIM_LOGGER_OBJECT_NAME
959 #define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate
compute_arguments_create(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)960 static int compute_arguments_create(
961 KIM_ModelCompute const * const modelCompute,
962 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
963 {
964 int ier;
965
966 (void) modelCompute; /* avoid unused parameter warning */
967
968 /* register arguments */
969 ier = KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
970 modelComputeArgumentsCreate,
971 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
972 KIM_SUPPORT_STATUS_optional)
973 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
974 modelComputeArgumentsCreate,
975 KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
976 KIM_SUPPORT_STATUS_optional)
977 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
978 modelComputeArgumentsCreate,
979 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
980 KIM_SUPPORT_STATUS_optional);
981
982 /* register callbacks */
983 ier = ier
984 || KIM_ModelComputeArgumentsCreate_SetCallbackSupportStatus(
985 modelComputeArgumentsCreate,
986 KIM_COMPUTE_CALLBACK_NAME_ProcessDEDrTerm,
987 KIM_SUPPORT_STATUS_optional)
988 || KIM_ModelComputeArgumentsCreate_SetCallbackSupportStatus(
989 modelComputeArgumentsCreate,
990 KIM_COMPUTE_CALLBACK_NAME_ProcessD2EDr2Term,
991 KIM_SUPPORT_STATUS_optional);
992
993 if (ier == TRUE)
994 {
995 LOG_ERROR("Unable to set argument supportStatus.");
996 return TRUE;
997 }
998 else
999 {
1000 return FALSE;
1001 }
1002 }
1003
1004 /*****************************************************************************/
1005 /* Compute arguments destroy function */
1006 /*****************************************************************************/
compute_arguments_destroy(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)1007 static int compute_arguments_destroy(
1008 KIM_ModelCompute const * const modelCompute,
1009 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
1010 {
1011 (void) modelCompute; /* avoid unused parameter warning */
1012 (void) modelComputeArgumentsDestroy;
1013
1014 /* Nothing further to do */
1015
1016 return FALSE;
1017 }
1018
1019 /*****************************************************************************/
1020 /* Create function */
1021 /*****************************************************************************/
1022 #undef KIM_LOGGER_FUNCTION_NAME
1023 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDriverCreate_LogEntry
1024 #undef KIM_LOGGER_OBJECT_NAME
1025 #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)1026 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
1027 KIM_LengthUnit const requestedLengthUnit,
1028 KIM_EnergyUnit const requestedEnergyUnit,
1029 KIM_ChargeUnit const requestedChargeUnit,
1030 KIM_TemperatureUnit const requestedTemperatureUnit,
1031 KIM_TimeUnit const requestedTimeUnit)
1032 {
1033 /* KIM variables */
1034 int numberOfParameterFiles;
1035 char const * paramfile1name;
1036
1037 /* Local variables */
1038 FILE * fid;
1039 KIM_SpeciesName speciesName;
1040
1041 double influenceDistance;
1042 double * cutoffs;
1043 double *A1s, *A2s, *A3s;
1044 double *B1s, *B2s, *B3s;
1045 double *r1s, *r2s, *r3s;
1046
1047 double theta;
1048 double dbldummy;
1049 double nextShift;
1050
1051 int ier;
1052 int i, j;
1053 struct model_buffer * buffer;
1054
1055 int numberModelSpecies;
1056 int numberUniqueSpeciesPairs;
1057
1058 int endOfFileFlag;
1059 int indx;
1060 char spec1[MAXLINE], spec2[MAXLINE], nextLine[MAXLINE];
1061 char dummy[12];
1062 char * nextLinePtr;
1063 double initialTemp;
1064 double nextCutoff;
1065 double nextA1, nextA2, nextA3;
1066 double nextB1, nextB2, nextB3;
1067 double nextr1, nextr2, nextr3;
1068
1069 nextLinePtr = nextLine;
1070
1071 /* using fixed units */
1072 ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
1073 KIM_LENGTH_UNIT_A,
1074 KIM_ENERGY_UNIT_eV,
1075 KIM_CHARGE_UNIT_unused,
1076 KIM_TEMPERATURE_UNIT_unused,
1077 KIM_TIME_UNIT_unused);
1078 if (ier == TRUE)
1079 {
1080 LOG_ERROR("Problem setting units");
1081 return ier;
1082 }
1083
1084 ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate,
1085 KIM_NUMBERING_zeroBased);
1086 if (ier == TRUE)
1087 {
1088 LOG_ERROR("Unable to set numbering");
1089 return ier;
1090 }
1091
1092 /* store pointer to functions in KIM object */
1093 KIM_ModelDriverCreate_SetRoutinePointer(modelDriverCreate,
1094 KIM_MODEL_ROUTINE_NAME_Destroy,
1095 KIM_LANGUAGE_NAME_c,
1096 TRUE,
1097 (KIM_Function *) destroy);
1098 KIM_ModelDriverCreate_SetRoutinePointer(
1099 modelDriverCreate,
1100 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsCreate,
1101 KIM_LANGUAGE_NAME_c,
1102 TRUE,
1103 (KIM_Function *) compute_arguments_create);
1104 KIM_ModelDriverCreate_SetRoutinePointer(
1105 modelDriverCreate,
1106 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsDestroy,
1107 KIM_LANGUAGE_NAME_c,
1108 TRUE,
1109 (KIM_Function *) compute_arguments_destroy);
1110 KIM_ModelDriverCreate_SetRoutinePointer(modelDriverCreate,
1111 KIM_MODEL_ROUTINE_NAME_Compute,
1112 KIM_LANGUAGE_NAME_c,
1113 TRUE,
1114 (KIM_Function *) compute);
1115 KIM_ModelDriverCreate_SetRoutinePointer(modelDriverCreate,
1116 KIM_MODEL_ROUTINE_NAME_Refresh,
1117 KIM_LANGUAGE_NAME_c,
1118 TRUE,
1119 (KIM_Function *) refresh);
1120
1121 /* get number of parameter files */
1122 KIM_ModelDriverCreate_GetNumberOfParameterFiles(modelDriverCreate,
1123 &numberOfParameterFiles);
1124 if (numberOfParameterFiles != 1)
1125 {
1126 ier = TRUE;
1127 LOG_ERROR("Incorrect number of parameter files. This driver expects only "
1128 "one Model parameter file.");
1129 return ier;
1130 }
1131
1132 /* get parameter file name */
1133 ier = KIM_ModelDriverCreate_GetParameterFileName(
1134 modelDriverCreate, 0, ¶mfile1name);
1135 if (ier == TRUE)
1136 {
1137 LOG_ERROR("Unable to get Model parameter file name.");
1138 return ier;
1139 }
1140
1141 fid = fopen(paramfile1name, "r");
1142 if (fid == NULL)
1143 {
1144 ier = TRUE;
1145 LOG_ERROR("Unable to open parameter file for Morse parameters");
1146 return ier;
1147 }
1148
1149 /* Read line 0 of parameter file */
1150 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1151 ier = sscanf(nextLine, "%d %lg", &numberModelSpecies, &initialTemp);
1152 if (ier != 2)
1153 {
1154 sprintf(nextLine, "unable to read first line of the parameter file");
1155 ier = TRUE;
1156 LOG_ERROR(nextLine);
1157 fclose(fid);
1158 return ier;
1159 }
1160
1161 /* Allocate memory used to store parameters read in */
1162 cutoffs = AllocateAndInitialize1DArray(numberModelSpecies);
1163 A1s = AllocateAndInitialize1DArray(numberModelSpecies);
1164 A2s = AllocateAndInitialize1DArray(numberModelSpecies);
1165 A3s = AllocateAndInitialize1DArray(numberModelSpecies);
1166 B1s = AllocateAndInitialize1DArray(numberModelSpecies);
1167 B2s = AllocateAndInitialize1DArray(numberModelSpecies);
1168 B3s = AllocateAndInitialize1DArray(numberModelSpecies);
1169 r1s = AllocateAndInitialize1DArray(numberModelSpecies);
1170 r2s = AllocateAndInitialize1DArray(numberModelSpecies);
1171 r3s = AllocateAndInitialize1DArray(numberModelSpecies);
1172
1173 numberUniqueSpeciesPairs
1174 = ((numberModelSpecies + 1) * numberModelSpecies) / 2;
1175 /* set all values of cutoffs1D to -1 for check later */
1176 for (i = 0; i < numberUniqueSpeciesPairs; i++) { cutoffs[i] = -1.0; }
1177
1178 /* Read all pure species lines. The species codes used corresponding to */
1179 /* the order in which they are read (starting with 0) */
1180 for (i = 0; i < numberModelSpecies; i++)
1181 {
1182 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1183 ier = sscanf(nextLine,
1184 "%s %s %lf %lf %lf %lf",
1185 spec1,
1186 spec2,
1187 &nextCutoff,
1188 &nextA1,
1189 &nextA2,
1190 &nextA3);
1191 if (ier != 6)
1192 {
1193 ier = TRUE;
1194 sprintf(nextLine, "error reading lines of the parameter file");
1195 LOG_ERROR(nextLine);
1196 fclose(fid);
1197 Deallocate1DArrays(
1198 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1199 return ier;
1200 }
1201
1202 if (strcmp(spec1, spec2) != 0)
1203 {
1204 ier = TRUE;
1205 sprintf(nextLine,
1206 "not all pure species interactions were given at "
1207 "the top of the parameter file");
1208 LOG_ERROR(nextLine);
1209 fclose(fid);
1210 Deallocate1DArrays(
1211 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1212 return ier;
1213 }
1214 else
1215 {
1216 /* Register species code in API */
1217 speciesName = KIM_SpeciesName_FromString(spec1);
1218 ier = KIM_ModelDriverCreate_SetSpeciesCode(
1219 modelDriverCreate, speciesName, i);
1220 if (ier == TRUE)
1221 {
1222 LOG_ERROR("Unable to set species code");
1223 fclose(fid);
1224 Deallocate1DArrays(
1225 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1226 return ier;
1227 }
1228
1229 /* Calculate the appropriate array indices to store these params at */
1230 indx = i * numberModelSpecies + i - (i * i + i) / 2;
1231
1232 /* Store this line of parameters in the appropriate places */
1233 cutoffs[indx] = nextCutoff;
1234 A1s[indx] = nextA1;
1235 A2s[indx] = nextA2;
1236 A3s[indx] = nextA3;
1237
1238 /* Now read the second line of parameters in this block and store */
1239 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1240 ier = sscanf(nextLine, "%lf %lf %lf", &nextB1, &nextB2, &nextB3);
1241 if (ier != 3)
1242 {
1243 ier = TRUE;
1244 sprintf(nextLine, "error reading lines of the parameter file");
1245 LOG_ERROR(nextLine);
1246 fclose(fid);
1247 Deallocate1DArrays(
1248 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1249 return ier;
1250 }
1251 B1s[indx] = nextB1;
1252 B2s[indx] = nextB2;
1253 B3s[indx] = nextB3;
1254
1255 /* Now read the third line of parameters in this block and store */
1256 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1257 ier = sscanf(nextLine, "%lf %lf %lf", &nextr1, &nextr2, &nextr3);
1258 if (ier != 3)
1259 {
1260 ier = TRUE;
1261 sprintf(nextLine, "error reading lines of the parameter file");
1262 LOG_ERROR(nextLine);
1263 fclose(fid);
1264 Deallocate1DArrays(
1265 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1266 return ier;
1267 }
1268 r1s[indx] = nextr1;
1269 r2s[indx] = nextr2;
1270 r3s[indx] = nextr3;
1271 } /* End check on whether species listed are the same */
1272 } /* End loop over number of species */
1273
1274
1275 /* Now read all of the interspecies parameters, assuming they are */
1276 /* ordered correctly: */
1277 /* (00, 01, ... ,0N, 10, 11, ... 1N, ..., N0, N1, ... NN) */
1278 for (i = 0; i < numberModelSpecies; i++)
1279 {
1280 for (j = i + 1; j < numberModelSpecies; j++)
1281 {
1282 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1283 ier = sscanf(nextLine,
1284 "%s %s %lf %lf %lf %lf",
1285 spec1,
1286 spec2,
1287 &nextCutoff,
1288 &nextA1,
1289 &nextA2,
1290 &nextA3);
1291 if (ier != 6)
1292 {
1293 ier = TRUE;
1294 sprintf(nextLine, "error reading lines of the parameter file");
1295 LOG_ERROR(nextLine);
1296 fclose(fid);
1297 Deallocate1DArrays(
1298 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1299 return ier;
1300 }
1301
1302 if (strcmp(spec1, spec2) == 0)
1303 {
1304 ier = TRUE;
1305 sprintf(nextLine,
1306 "more than one listing found for pure species "
1307 "interactions for species code %s",
1308 spec1);
1309 LOG_ERROR(nextLine);
1310 fclose(fid);
1311 Deallocate1DArrays(
1312 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1313 return ier;
1314 }
1315 else
1316 {
1317 /* Calculate the appropriate array indices to store these params at */
1318 indx = i * numberModelSpecies + j - (i * i + i) / 2;
1319
1320 /* Store this line of parameters in the appropriate places */
1321 cutoffs[indx] = nextCutoff;
1322 A1s[indx] = nextA1;
1323 A2s[indx] = nextA2;
1324 A3s[indx] = nextA3;
1325
1326 /* Now read the second line of parameters in this block and store */
1327 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1328 ier = sscanf(nextLine, "%lf %lf %lf", &nextB1, &nextB2, &nextB3);
1329 if (ier != 3)
1330 {
1331 ier = TRUE;
1332 sprintf(nextLine, "error reading lines of the parameter file");
1333 LOG_ERROR(nextLine);
1334 fclose(fid);
1335 Deallocate1DArrays(
1336 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1337 return ier;
1338 }
1339 B1s[indx] = nextB1;
1340 B2s[indx] = nextB2;
1341 B3s[indx] = nextB3;
1342
1343 /* Now read the third line of parameters in this block and store */
1344 getNextDataLine(fid, nextLinePtr, MAXLINE, &endOfFileFlag);
1345 ier = sscanf(nextLine, "%lf %lf %lf", &nextr1, &nextr2, &nextr3);
1346 if (ier != 3)
1347 {
1348 ier = TRUE;
1349 sprintf(nextLine, "error reading lines of the parameter file");
1350 LOG_ERROR(nextLine);
1351 fclose(fid);
1352 Deallocate1DArrays(
1353 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1354 return ier;
1355 }
1356 r1s[indx] = nextr1;
1357 r2s[indx] = nextr2;
1358 r3s[indx] = nextr3;
1359
1360 } /* End check on whether species listed are the same */
1361 } /* End inner loop over unique species combinations */
1362 } /* End outer loop over number of species */
1363
1364 /* close parameter file */
1365 fclose(fid);
1366
1367 /* Check that we got parameters for all pairs */
1368 ier = FALSE;
1369 sprintf(nextLine, "There are not values for the following pairs: \n");
1370 for (i = 0; i < numberModelSpecies; i++)
1371 {
1372 for (j = i; j < numberModelSpecies; j++)
1373 {
1374 indx = i * numberModelSpecies + j - (i * i + i) / 2;
1375 if (cutoffs[indx] == -1.0)
1376 {
1377 sprintf(dummy, "%d and %d\n", i, j);
1378 strcat(nextLine, dummy);
1379 ier = TRUE;
1380 }
1381 }
1382 }
1383 if (ier == TRUE)
1384 {
1385 LOG_ERROR(nextLine);
1386 Deallocate1DArrays(
1387 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1388 return ier;
1389 }
1390
1391 /* convert parameters to appropriate units (in-place) */
1392 ier = ConvertUnits(modelDriverCreate,
1393 requestedLengthUnit,
1394 requestedEnergyUnit,
1395 requestedChargeUnit,
1396 requestedTemperatureUnit,
1397 requestedTimeUnit,
1398 numberUniqueSpeciesPairs,
1399 cutoffs,
1400 A1s,
1401 A2s,
1402 r1s,
1403 r2s);
1404 if (ier == TRUE)
1405 {
1406 sprintf(nextLine, "failed to convert units");
1407 LOG_ERROR(nextLine);
1408 Deallocate1DArrays(
1409 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1410 return ier;
1411 }
1412
1413 /* allocate buffer */
1414 buffer = (struct model_buffer *) malloc(sizeof(struct model_buffer));
1415 if (NULL == buffer)
1416 {
1417 ier = TRUE;
1418 LOG_ERROR("Coul not allocate memory for Model buffer");
1419 Deallocate1DArrays(
1420 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1421 return ier;
1422 }
1423
1424 /* register model buffer */
1425 KIM_ModelDriverCreate_SetModelBufferPointer(modelDriverCreate,
1426 (void *) buffer);
1427
1428 /* Set number of species in buffer */
1429 buffer->numberModelSpecies = numberModelSpecies;
1430
1431 /* Allocate arrays within buffer */
1432 buffer->cutsq2D
1433 = AllocateAndInitialize2DArray(numberModelSpecies, numberModelSpecies);
1434 buffer->As2D
1435 = AllocateAndInitialize2DArray(numberModelSpecies, numberModelSpecies);
1436 buffer->Bs2D
1437 = AllocateAndInitialize2DArray(numberModelSpecies, numberModelSpecies);
1438 buffer->rHats2D
1439 = AllocateAndInitialize2DArray(numberModelSpecies, numberModelSpecies);
1440 buffer->shifts2D
1441 = AllocateAndInitialize2DArray(numberModelSpecies, numberModelSpecies);
1442
1443 /* set buffer parameters to those read in */
1444 buffer->temperature = initialTemp;
1445 buffer->cutoffs = cutoffs;
1446 buffer->A1s = A1s;
1447 buffer->A2s = A2s;
1448 buffer->A3s = A3s;
1449 buffer->B1s = B1s;
1450 buffer->B2s = B2s;
1451 buffer->B3s = B3s;
1452 buffer->r1s = r1s;
1453 buffer->r2s = r2s;
1454 buffer->r3s = r3s;
1455
1456 /* Request that simulator not provide neighbors of padding atoms */
1457 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles = 1;
1458
1459 /* Register params for mutability by simulator */
1460 ier = KIM_ModelDriverCreate_SetParameterPointerDouble(
1461 modelDriverCreate,
1462 1,
1463 &(buffer->temperature),
1464 "temperature",
1465 "Pseudo-temperature at which to evaluate the potential parameters.")
1466 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1467 modelDriverCreate,
1468 numberUniqueSpeciesPairs,
1469 buffer->cutoffs,
1470 "cutoffs",
1471 "Distances used to determine whether two-body interactions for "
1472 "a pair of "
1473 "atoms occur. "
1474 "This array corresponds to an upper-triangular matrix in "
1475 "row-major storage. "
1476 "Ordering is according to SpeciesCode values. For example, to "
1477 "find "
1478 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1479 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1480 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1481 modelDriverCreate,
1482 numberUniqueSpeciesPairs,
1483 buffer->A1s,
1484 "A1s",
1485 "Additive constant in the calculation of A(theta). "
1486 "This array corresponds to an upper-triangular matrix in "
1487 "row-major storage. "
1488 "Ordering is according to SpeciesCode values. For example, to "
1489 "find "
1490 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1491 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1492 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1493 modelDriverCreate,
1494 numberUniqueSpeciesPairs,
1495 buffer->A2s,
1496 "A2s",
1497 "Multiplicative constant in the calculation of A(theta). "
1498 "This array corresponds to an upper-triangular matrix in "
1499 "row-major storage. "
1500 "Ordering is according to SpeciesCode values. For example, to "
1501 "find "
1502 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1503 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1504 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1505 modelDriverCreate,
1506 numberUniqueSpeciesPairs,
1507 buffer->A3s,
1508 "A3s",
1509 "Argument of the exponential in the calculation of A(theta). "
1510 "This array corresponds to an upper-triangular matrix in "
1511 "row-major storage. "
1512 "Ordering is according to SpeciesCode values. For example, to "
1513 "find "
1514 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1515 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1516 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1517 modelDriverCreate,
1518 numberUniqueSpeciesPairs,
1519 buffer->B1s,
1520 "B1s",
1521 "Additive constant in the calculation of B(theta). "
1522 "This array corresponds to an upper-triangular matrix in "
1523 "row-major storage. "
1524 "Ordering is according to SpeciesCode values. For example, to "
1525 "find "
1526 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1527 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1528 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1529 modelDriverCreate,
1530 numberUniqueSpeciesPairs,
1531 buffer->B2s,
1532 "B2s",
1533 "Multiplicative constant in the calculation of B(theta). "
1534 "This array corresponds to an upper-triangular matrix in "
1535 "row-major storage. "
1536 "Ordering is according to SpeciesCode values. For example, to "
1537 "find "
1538 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1539 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1540 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1541 modelDriverCreate,
1542 numberUniqueSpeciesPairs,
1543 buffer->B3s,
1544 "B3s",
1545 "Argument of the exponential in the calculation of B(theta). "
1546 "This array corresponds to an upper-triangular matrix in "
1547 "row-major storage. "
1548 "Ordering is according to SpeciesCode values. For example, to "
1549 "find "
1550 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1551 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1552 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1553 modelDriverCreate,
1554 numberUniqueSpeciesPairs,
1555 buffer->r1s,
1556 "r1s",
1557 "Additive constant in the calculation of rHat(theta). "
1558 "This array corresponds to an upper-triangular matrix in "
1559 "row-major storage. "
1560 "Ordering is according to SpeciesCode values. For example, to "
1561 "find "
1562 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1563 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1564 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1565 modelDriverCreate,
1566 numberUniqueSpeciesPairs,
1567 buffer->r2s,
1568 "r2s",
1569 "Multiplicative constant in the calculation of rHat(theta). "
1570 "This array corresponds to an upper-triangular matrix in "
1571 "row-major storage. "
1572 "Ordering is according to SpeciesCode values. For example, to "
1573 "find "
1574 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1575 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).")
1576 || KIM_ModelDriverCreate_SetParameterPointerDouble(
1577 modelDriverCreate,
1578 numberUniqueSpeciesPairs,
1579 buffer->r3s,
1580 "r3s",
1581 "Multiplicative constant contained inside of the exponential "
1582 "term "
1583 "in the calculation of rHat(theta). "
1584 "This array corresponds to an upper-triangular matrix in "
1585 "row-major storage. "
1586 "Ordering is according to SpeciesCode values. For example, to "
1587 "find "
1588 "the parameter related to SpeciesCode 'i' and SpeciesCode 'j' "
1589 "(i <= j), use (zero-based) index = (i*N + j - (i*i + i)/2).");
1590 if (ier == TRUE)
1591 {
1592 LOG_ERROR("Could not register parameter pointers");
1593 Deallocate1DArrays(
1594 10, cutoffs, A1s, A2s, A3s, B1s, B2s, B3s, r1s, r2s, r3s);
1595 Deallocate2DArrays(5,
1596 buffer->cutsq2D,
1597 buffer->As2D,
1598 buffer->Bs2D,
1599 buffer->rHats2D,
1600 buffer->shifts2D);
1601 free(buffer);
1602 return ier;
1603 }
1604
1605 /* Do some processing to set up 2D arrays (same stuff as is done in */
1606 /* refresh() */
1607 theta = buffer->temperature / 333.15;
1608 influenceDistance = 0.0;
1609 for (i = 0; i < numberModelSpecies; ++i)
1610 {
1611 for (j = i; j < numberModelSpecies; ++j)
1612 {
1613 indx = i * numberModelSpecies + j - (i * i + i) / 2;
1614
1615 if (influenceDistance < cutoffs[indx])
1616 { influenceDistance = cutoffs[indx]; }
1617
1618 buffer->cutsq2D[i][j] = buffer->cutsq2D[j][i]
1619 = (cutoffs[indx] * cutoffs[indx]);
1620
1621 buffer->As2D[i][j] = buffer->As2D[j][i]
1622 = calc_A(&A1s[indx], &A2s[indx], &A3s[indx], theta);
1623
1624 buffer->Bs2D[i][j] = buffer->Bs2D[j][i]
1625 = calc_B(&B1s[indx], &B2s[indx], &B3s[indx], theta);
1626
1627 buffer->rHats2D[i][j] = buffer->rHats2D[j][i]
1628 = calc_rHat(&r1s[indx], &r2s[indx], &r3s[indx], theta);
1629 }
1630 }
1631 /* Set influence distance */
1632 buffer->influenceDistance = influenceDistance;
1633
1634 /* Set influence distance pointer and cutoff list pointer*/
1635 KIM_ModelDriverCreate_SetInfluenceDistancePointer(
1636 modelDriverCreate, &(buffer->influenceDistance));
1637 KIM_ModelDriverCreate_SetNeighborListPointers(
1638 modelDriverCreate,
1639 1, /* Use only a single neighbor list */
1640 &(buffer->influenceDistance),
1641 (const int *) &(
1642 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
1643
1644 /* Set Values for Shifts */
1645 dbldummy = 0.0;
1646 influenceDistance += 1.0; /* add a bit to avoid rounding problem */
1647 for (i = 0; i < numberModelSpecies; i++)
1648 {
1649 for (j = i; j < numberModelSpecies; j++)
1650 {
1651 indx = i * numberModelSpecies + j - (i * i + i) / 2;
1652 /* call calc_phi with r=cutoff and shift=0.0 */
1653 calc_phi(&(buffer->As2D[i][j]),
1654 &(buffer->Bs2D[i][j]),
1655 &(buffer->rHats2D[i][j]),
1656 &dbldummy,
1657 &influenceDistance,
1658 cutoffs[indx],
1659 &nextShift);
1660 /* set shift to -shift */
1661 buffer->shifts2D[i][j] = buffer->shifts2D[j][i] = -nextShift;
1662 }
1663 }
1664
1665 ier = FALSE;
1666 return ier;
1667 }
1668