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) 2019, Regents of the University of Minnesota. */
25 /* All rights reserved. */
26 /* */
27 /* Contributors: */
28 /* Daniel S. Karls */
29 /* Ellad B. Tadmor */
30 /* Ryan S. Elliott */
31 /* */
32
33 #include "KIM_LogMacros.h"
34 #include "KIM_ModelDriverHeaders.h"
35 #include <math.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 #include "bondorder.inc"
41
42 #define TRUE 1
43 #define FALSE 0
44
45 /******************************************************************************/
46 /* Below are the definitions for some constants */
47 /******************************************************************************/
48 #define DIM 3 /* dimensionality of space */
49 #define SPECCODE 1 /* internal species code */
50 #define SPEC_NAME_LEN 64 /* max length of species name string */
51 #define SPEC_NAME_FMT "%63s" /* scan format of species name string */
52
53 /* Define prototype for Model Driver init */
54 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
55 KIM_LengthUnit const requestedLengthUnit,
56 KIM_EnergyUnit const requestedEnergyUnit,
57 KIM_ChargeUnit const requestedChargeUnit,
58 KIM_TemperatureUnit const requestedTemperatureUnit,
59 KIM_TimeUnit const requestedTimeUnit);
60
61 /* Define prototypes for destroy */
62 /* defined as static to avoid namespace clashes with other codes */
63 static int destroy_routine(KIM_ModelDestroy * const modelDestroy);
64
65 /* Define prototype for routines */
66 static int
67 compute_routine(KIM_ModelCompute const * const modelCompute,
68 KIM_ModelComputeArguments const * const modelComputeArguments);
69 static int compute_arguments_create(
70 KIM_ModelCompute const * const modelCompute,
71 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate);
72 static int compute_arguments_destroy(
73 KIM_ModelCompute const * const modelCompute,
74 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy);
75 static int refresh_routine(KIM_ModelRefresh * const modelRefresh);
76 static int
77 write_parameterized_model(KIM_ModelWriteParameterizedModel const * const
78 modelWriteParameterizedModel);
79
80 /* Define model_buffer structure */
81 struct model_buffer
82 {
83 double influenceDistance;
84 double cutoff;
85 double cut_sq;
86 int modelWillNotRequestNeighborsOfNoncontributingParticles;
87 char species_name[SPEC_NAME_LEN];
88 double * params;
89 };
90
91 /* compute function */
92 #undef KIM_LOGGER_FUNCTION_NAME
93 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelCompute_LogEntry
94 #undef KIM_LOGGER_OBJECT_NAME
95 #define KIM_LOGGER_OBJECT_NAME modelCompute
96 static int
compute_routine(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArguments const * const modelComputeArguments)97 compute_routine(KIM_ModelCompute const * const modelCompute,
98 KIM_ModelComputeArguments const * const modelComputeArguments)
99 {
100 /* local variables */
101 double tmp;
102 double Rij;
103 double Rik;
104 double Rjk;
105 double Ril;
106 double Rsqij;
107 double Rsqik;
108 double Rsqjk;
109 double Rsqil;
110 double Z;
111 double Z_i;
112 double *dZi_dR;
113 double bij;
114 double zeta_ij;
115 double phi2;
116 double dphi2_dRij;
117 double dphi2_dbij;
118 double dbij_dzeta_ij;
119 double dbij_dZi;
120 double phi2_force_contrib;
121 double coordination_force_contrib;
122 double phi3;
123 double dphi3_dRij;
124 double dphi3_dRik;
125 double dphi3_dRjk;
126 double phi3_force_contrib_ij;
127 double phi3_force_contrib_ik;
128 double phi3_force_contrib_jk;
129 double r_ij[DIM];
130 double r_ik[DIM];
131 double r_jk[DIM];
132 double r_il[DIM];
133 int ier;
134 int atom_i;
135 int neigh_count;
136 int neigh_count2;
137 int atom_j;
138 int atom_k;
139 int atom_l;
140 int dim;
141 int dZi_dR_size;
142 int const * neigh_list_of_current_part;
143 struct model_buffer const * buffer;
144 int comp_energy;
145 int comp_force;
146
147 int const * n_parts;
148 int const * particle_species_codes;
149 int const * particle_contributing;
150 double const * cut_sq;
151 double const * params;
152 double const * coords;
153 double * energy;
154 double * force;
155 int numb_of_part_neigh;
156
157 /* get buffer from KIM object */
158 KIM_ModelCompute_GetModelBufferPointer(modelCompute, (void **) &buffer);
159
160 /* unpack info from the buffer */
161 cut_sq = &(buffer->cut_sq);
162 params = buffer->params;
163
164 ier = KIM_ModelComputeArguments_GetArgumentPointerInteger(
165 modelComputeArguments,
166 KIM_COMPUTE_ARGUMENT_NAME_numberOfParticles,
167 (int **) &n_parts)
168 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
169 modelComputeArguments,
170 KIM_COMPUTE_ARGUMENT_NAME_particleSpeciesCodes,
171 (int **) &particle_species_codes)
172 || KIM_ModelComputeArguments_GetArgumentPointerInteger(
173 modelComputeArguments,
174 KIM_COMPUTE_ARGUMENT_NAME_particleContributing,
175 (int **) &particle_contributing)
176 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
177 modelComputeArguments,
178 KIM_COMPUTE_ARGUMENT_NAME_coordinates,
179 (double **) &coords)
180 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
181 modelComputeArguments,
182 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
183 (double **) &energy)
184 || KIM_ModelComputeArguments_GetArgumentPointerDouble(
185 modelComputeArguments,
186 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
187 (double **) &force);
188 if (ier)
189 {
190 LOG_ERROR("Unable to get argument pointer.");
191 return ier;
192 }
193
194 comp_energy = (energy != NULL);
195 comp_force = (force != NULL);
196
197 /* Check to be sure that the species are correct */
198 /**/
199 ier = TRUE; /* assume an error */
200 for (atom_i = 0; atom_i < *n_parts; ++atom_i)
201 {
202 if (SPECCODE != particle_species_codes[atom_i])
203 {
204 LOG_ERROR("Unexpected species code detected.");
205 return ier;
206 }
207 }
208 ier = FALSE; /* everything is ok */
209
210 /* initialize potential energies, forces, and virial term */
211 if (comp_energy) { *energy = 0.0; }
212
213 if (comp_force)
214 {
215 for (atom_i = 0; atom_i < *n_parts; ++atom_i)
216 {
217 for (dim = 0; dim < DIM; ++dim) { force[atom_i * DIM + dim] = 0.0; }
218 }
219 }
220
221 /* perform initial allocation of array containing derivatives of Z_i
222 * w.r.t. Rij for each neighbor j of atom i */
223 dZi_dR_size = 100;
224 dZi_dR = (double*) malloc(sizeof(double) * dZi_dR_size);
225
226 /* Compute energy and forces */
227
228 /* loop over particles and compute energy and forces */
229 for (atom_i = 0; atom_i < *n_parts; ++atom_i)
230 {
231 if (particle_contributing[atom_i])
232 {
233 ier = KIM_ModelComputeArguments_GetNeighborList(
234 modelComputeArguments,
235 0,
236 atom_i,
237 &numb_of_part_neigh,
238 &neigh_list_of_current_part);
239 if (ier)
240 {
241 /* some sort of problem, exit */
242 LOG_ERROR("Unable to get neighbor list.");
243
244 /* free dynamic memory */
245 free(dZi_dR);
246 ier = TRUE;
247 return ier;
248 }
249
250 /* reset coordination */
251 Z_i = 0.0;
252
253 /* check if number of neighbors exceeds our current array size for
254 * storing derivatives of coordination w.r.t. each neighbor */
255 if (numb_of_part_neigh > dZi_dR_size)
256 {
257 free(dZi_dR);
258 while (dZi_dR_size < numb_of_part_neigh)
259 {
260 dZi_dR_size *= 2;
261 }
262 dZi_dR = (double*) malloc(sizeof(double) * dZi_dR_size);
263 }
264
265 /* first loop over the neighbors of particle atom_i to calculate
266 * coordination */
267 for (neigh_count = 0; neigh_count < numb_of_part_neigh; ++neigh_count)
268 {
269 atom_j = neigh_list_of_current_part[neigh_count]; /* get ID */
270
271 /* compute relative position vector between atoms i and j, and squared
272 * distance */
273 Rsqij = 0.0;
274 for (dim = 0; dim < DIM; ++dim)
275 {
276 r_ij[dim] = coords[atom_j * DIM + dim] - coords[atom_i * DIM + dim];
277 /* compute squared distance */
278 Rsqij += r_ij[dim] * r_ij[dim];
279 }
280
281 /* compute energy and force */
282 if (Rsqij < *cut_sq)
283 {
284 /* particles i and j are interacting */
285 Rij = sqrt(Rsqij);
286
287 if (comp_force)
288 {
289 /* compute coordination and its derivative w.r.t. Rij */
290 calc_Zi_dZi(params, Rij, &Z, &dZi_dR[neigh_count]);
291 }
292 else
293 {
294 /* compute just coordination */
295 calc_Zi_dZi(params, Rij, &Z, NULL);
296 }
297 Z_i += Z;
298 }
299 }
300
301 /* Second loop over the neighbors of particle atom_i to calculate energy
302 * & forces */
303 for (neigh_count = 0; neigh_count < numb_of_part_neigh; ++neigh_count)
304 {
305 atom_j = neigh_list_of_current_part[neigh_count]; /* get ID */
306
307 /* compute relative position vector between atoms i and j, and squared
308 * distance */
309 Rsqij = 0.0;
310 for (dim = 0; dim < DIM; ++dim)
311 {
312 r_ij[dim] = coords[atom_j * DIM + dim] - coords[atom_i * DIM + dim];
313 /* compute squared distance */
314 Rsqij += r_ij[dim] * r_ij[dim];
315 }
316
317 /* reset bond order */
318 zeta_ij = 0.0;
319
320 /* compute energy and force */
321 if (Rsqij < *cut_sq)
322 {
323 /* particles i and j are interacting */
324 Rij = sqrt(Rsqij);
325
326 /* First secondary loop over neighbors to calculate bond order */
327 for (neigh_count2 = 0; neigh_count2 < numb_of_part_neigh;
328 ++neigh_count2)
329 {
330 atom_k = neigh_list_of_current_part[neigh_count2]; /* get ID */
331
332 if (atom_k == atom_j) { continue; }
333
334 /* compute relative position vector between atoms i and k, and
335 * squared distance */
336 Rsqik = 0.0;
337 for (dim = 0; dim < DIM; ++dim)
338 {
339 r_ik[dim]
340 = coords[atom_k * DIM + dim] - coords[atom_i * DIM + dim];
341 /* compute squared distance */
342 Rsqik += r_ik[dim] * r_ik[dim];
343 }
344
345 /* compute energy and force */
346 if (Rsqik < *cut_sq)
347 {
348 /* particles i and k are interacting */
349 Rik = sqrt(Rsqik);
350
351 /* compute relative position vector between atoms j and k, and
352 * squared distance */
353 Rsqjk = 0.0;
354 for (dim = 0; dim < DIM; ++dim)
355 {
356 r_jk[dim]
357 = coords[atom_k * DIM + dim] - coords[atom_j * DIM + dim];
358 /* compute squared distance */
359 Rsqjk += r_jk[dim] * r_jk[dim];
360 }
361
362 if (SKIP_CHECK_ON_RJK || Rsqjk < *cut_sq)
363 {
364 Rjk = sqrt(Rsqjk);
365
366 /* compute just three-body interaction */
367 calc_phi3_dphi3(params, Rij, Rik, Rjk, &phi3, NULL, NULL, NULL);
368
369 /* contribution to bond order */
370 zeta_ij += phi3;
371 }
372 }
373 } /* loop on neigh_count2 */
374
375 if (comp_force)
376 {
377 /* compute bond order and its derivative w.r.t. zeta_ij */
378 calc_bij_dbij(params, Z_i, zeta_ij, &bij, &dbij_dZi, &dbij_dzeta_ij);
379
380 /* compute two-body energy and its derivatives w.r.t. Rij and bij */
381 calc_phi2_dphi2(params, Rij, bij, &phi2, &dphi2_dRij, &dphi2_dbij);
382 }
383 else
384 {
385 /* compute just bond order */
386 calc_bij_dbij(params, Z_i, zeta_ij, &bij, NULL, NULL);
387
388 /* compute just two-body energy */
389 calc_phi2_dphi2(params, Rij, bij, &phi2, NULL, NULL);
390 }
391
392 /* contribution to energy */
393 if (comp_energy) { *energy += 0.5 * phi2; }
394
395 /* contribution to forces */
396 if (comp_force)
397 {
398 for (dim = 0; dim < DIM; ++dim)
399 {
400 phi2_force_contrib = 0.5 * dphi2_dRij * r_ij[dim] / Rij;
401
402 force[atom_i * DIM + dim]
403 += phi2_force_contrib; /* accumulate force on atom_i */
404 force[atom_j * DIM + dim]
405 -= phi2_force_contrib; /* accumulate force on atom_j */
406 }
407
408 /* Another secondary loop to apply coordination forces */
409 for (neigh_count2 = 0; neigh_count2 < numb_of_part_neigh;
410 ++neigh_count2)
411 {
412 atom_l = neigh_list_of_current_part[neigh_count2]; /* get ID */
413
414 /* compute relative position vector between atoms i and l, and
415 * squared distance */
416 Rsqil = 0.0;
417 for (dim = 0; dim < DIM; ++dim)
418 {
419 r_il[dim]
420 = coords[atom_l * DIM + dim] - coords[atom_i * DIM + dim];
421 /* compute squared distance */
422 Rsqil += r_il[dim] * r_il[dim];
423 }
424
425 /* compute energy and force */
426 if (Rsqil < *cut_sq)
427 {
428 /* particles i and l are interacting */
429 Ril = sqrt(Rsqil);
430
431 /* contribution to forces */
432 for (dim = 0; dim < DIM; ++dim)
433 {
434 coordination_force_contrib = 0.5 * dphi2_dbij * dbij_dZi *
435 dZi_dR[neigh_count2] * r_il[dim]/Ril;
436
437 force[atom_i * DIM + dim]
438 += coordination_force_contrib;
439
440 force[atom_l * DIM + dim]
441 -= coordination_force_contrib;
442 }
443 }
444 }
445
446 /* Another, final, secondary loop over neighbors to calculate three-body
447 * energy and forces */
448 for (neigh_count2 = 0; neigh_count2 < numb_of_part_neigh;
449 ++neigh_count2)
450 {
451 atom_k = neigh_list_of_current_part[neigh_count2]; /* get ID */
452
453 if (atom_k == atom_j) { continue; }
454
455 /* compute relative position vector between atoms i and k, and
456 * squared distance */
457 Rsqik = 0.0;
458 for (dim = 0; dim < DIM; ++dim)
459 {
460 r_ik[dim]
461 = coords[atom_k * DIM + dim] - coords[atom_i * DIM + dim];
462 /* compute squared distance */
463 Rsqik += r_ik[dim] * r_ik[dim];
464 }
465
466 /* compute energy and force */
467 if (Rsqik < *cut_sq)
468 {
469 /* particles i and k are interacting */
470 Rik = sqrt(Rsqik);
471
472 /* compute relative position vector between atoms j and k, and
473 * squared distance */
474 Rsqjk = 0.0;
475 for (dim = 0; dim < DIM; ++dim)
476 {
477 r_jk[dim]
478 = coords[atom_k * DIM + dim] - coords[atom_j * DIM + dim];
479 /* compute squared distance */
480 Rsqjk += r_jk[dim] * r_jk[dim];
481 }
482
483 if (SKIP_CHECK_ON_RJK || Rsqjk < *cut_sq)
484 {
485 Rjk = sqrt(Rsqjk);
486
487 /* compute three-body interaction and its derivatives with
488 * respect to Rij and Rik */
489 calc_phi3_dphi3(params,
490 Rij,
491 Rik,
492 Rjk,
493 &phi3,
494 &dphi3_dRij,
495 &dphi3_dRik,
496 &dphi3_dRjk);
497
498 /* contribution to forces */
499 for (dim = 0; dim < DIM; ++dim)
500 {
501 tmp = 0.5 * dphi2_dbij * dbij_dzeta_ij;
502 phi3_force_contrib_ij = dphi3_dRij * r_ij[dim] / Rij;
503 phi3_force_contrib_ik = dphi3_dRik * r_ik[dim] / Rik;
504 phi3_force_contrib_jk = dphi3_dRjk * r_jk[dim] / Rjk;
505
506 force[atom_i * DIM + dim]
507 += tmp * (phi3_force_contrib_ij
508 + phi3_force_contrib_ik); /* accumulate phi3
509 derivatives for
510 atom_i */
511 force[atom_j * DIM + dim]
512 += tmp * (-phi3_force_contrib_ij
513 + phi3_force_contrib_jk); /* accumulate three-body
514 force on atom_j */
515 force[atom_k * DIM + dim]
516 += tmp * (-phi3_force_contrib_ik
517 - phi3_force_contrib_jk); /* accumulate partial
518 three-body force on
519 atom_k */
520 }
521 } /* if Rsqjk < *cut_sq */
522 } /* if Rsqik < *cut_sq */
523 } /* loop on neigh_count2 */
524 } /* check on whether forces were asked for */
525 } /* if Rsqij < *cut_sq */
526 } /* loop on neigh_count */
527 } /* if particle_contributing */
528 } /* loop on atom_i */
529
530 /* free the dynamic array used to store derivatives of coordination w.r.t.
531 * Rij for each neighbor */
532 free(dZi_dR);
533
534 /* everything is great */
535 ier = FALSE;
536
537 return ier;
538 }
539
540 /* Create function */
541 #undef KIM_LOGGER_FUNCTION_NAME
542 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelDriverCreate_LogEntry
543 #undef KIM_LOGGER_OBJECT_NAME
544 #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)545 int model_driver_create(KIM_ModelDriverCreate * const modelDriverCreate,
546 KIM_LengthUnit const requestedLengthUnit,
547 KIM_EnergyUnit const requestedEnergyUnit,
548 KIM_ChargeUnit const requestedChargeUnit,
549 KIM_TemperatureUnit const requestedTemperatureUnit,
550 KIM_TimeUnit const requestedTimeUnit)
551 {
552 /* KIM variables */
553 int number_of_parameter_files;
554 char const * param_file_1_name;
555
556 /* Local variables */
557 FILE * fid;
558 char species_name_string[SPEC_NAME_LEN];
559 char strbuf[256];
560 KIM_SpeciesName species_name;
561
562 int ier;
563 int param_count;
564 double * params;
565 double conversion_factor;
566 struct model_buffer * buffer;
567
568 /* Use function pointer definitions to verify prototypes */
569 KIM_ModelDriverCreateFunction * create = model_driver_create;
570 KIM_ModelComputeArgumentsCreateFunction * ca_create
571 = compute_arguments_create;
572 KIM_ModelComputeFunction * compute = compute_routine;
573 KIM_ModelRefreshFunction * refresh = refresh_routine;
574 KIM_ModelWriteParameterizedModelFunction * write_model
575 = write_parameterized_model;
576 KIM_ModelComputeArgumentsDestroyFunction * ca_destroy
577 = compute_arguments_destroy;
578 KIM_ModelDestroyFunction * destroy = destroy_routine;
579
580
581 (void) create; /* avoid unused parameter warnings */
582 (void) requestedChargeUnit;
583 (void) requestedTemperatureUnit;
584 (void) requestedTimeUnit;
585
586 /* using requested units */
587 ier = KIM_ModelDriverCreate_SetUnits(modelDriverCreate,
588 requestedLengthUnit,
589 requestedEnergyUnit,
590 KIM_CHARGE_UNIT_unused,
591 KIM_TEMPERATURE_UNIT_unused,
592 KIM_TIME_UNIT_unused);
593 if (ier == TRUE)
594 {
595 LOG_ERROR("Unable to set units.");
596 return ier;
597 }
598
599 ier = KIM_ModelDriverCreate_SetModelNumbering(modelDriverCreate,
600 KIM_NUMBERING_zeroBased);
601 if (ier == TRUE)
602 {
603 LOG_ERROR("Unable to set numbering.");
604 return ier;
605 }
606
607 /* store pointer to functions in KIM object */
608 ier = KIM_ModelDriverCreate_SetRoutinePointer(
609 modelDriverCreate,
610 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsCreate,
611 KIM_LANGUAGE_NAME_c,
612 TRUE,
613 (KIM_Function *) ca_create)
614 || KIM_ModelDriverCreate_SetRoutinePointer(
615 modelDriverCreate,
616 KIM_MODEL_ROUTINE_NAME_Compute,
617 KIM_LANGUAGE_NAME_c,
618 TRUE,
619 (KIM_Function *) compute)
620 || KIM_ModelDriverCreate_SetRoutinePointer(
621 modelDriverCreate,
622 KIM_MODEL_ROUTINE_NAME_Refresh,
623 KIM_LANGUAGE_NAME_c,
624 TRUE,
625 (KIM_Function *) refresh)
626 || KIM_ModelDriverCreate_SetRoutinePointer(
627 modelDriverCreate,
628 KIM_MODEL_ROUTINE_NAME_WriteParameterizedModel,
629 KIM_LANGUAGE_NAME_c,
630 FALSE,
631 (KIM_Function *) write_model)
632 || KIM_ModelDriverCreate_SetRoutinePointer(
633 modelDriverCreate,
634 KIM_MODEL_ROUTINE_NAME_ComputeArgumentsDestroy,
635 KIM_LANGUAGE_NAME_c,
636 TRUE,
637 (KIM_Function *) ca_destroy)
638 || KIM_ModelDriverCreate_SetRoutinePointer(
639 modelDriverCreate,
640 KIM_MODEL_ROUTINE_NAME_Destroy,
641 KIM_LANGUAGE_NAME_c,
642 TRUE,
643 (KIM_Function *) destroy);
644
645 if (ier == TRUE)
646 {
647 LOG_ERROR("Unable to register model function pointers.");
648 return ier;
649 }
650
651 /* get number of parameter files */
652 KIM_ModelDriverCreate_GetNumberOfParameterFiles(modelDriverCreate,
653 &number_of_parameter_files);
654 /* set param_file_1_name */
655 if (number_of_parameter_files != 1)
656 {
657 ier = TRUE;
658 LOG_ERROR("Incorrect number of parameter files.");
659 return ier;
660 }
661 ier = KIM_ModelDriverCreate_GetParameterFileName(
662 modelDriverCreate, 0, ¶m_file_1_name);
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(param_file_1_name, "r");
671 if (fid == NULL)
672 {
673 ier = TRUE;
674 LOG_ERROR("Unable to open parameter file.");
675 return ier;
676 }
677
678 ier = fscanf(fid, SPEC_NAME_FMT, species_name_string); /* species symbol */
679
680 /* Check that we successfully read in this line */
681 if (1 != ier)
682 {
683 ier = TRUE;
684 fclose(fid);
685 LOG_ERROR("Unable to read species from parameter file.");
686 }
687
688 /* Consume the rest of the species line */
689 if (fgets(strbuf, 256, fid) == NULL)
690 {
691 if (ferror(fid))
692 {
693 /* A read error was encountered */
694 ier = TRUE;
695 fclose(fid);
696 LOG_ERROR("Error occurred while trying to read parameter file (ferror).");
697 return ier;
698 }
699 else if (feof(fid))
700 {
701 /* Reached the EOF mark before reading a character */
702 ier = TRUE;
703 fclose(fid);
704 LOG_ERROR("Error occurred while trying to read parameter file (reached EOF).");
705 return ier;
706 }
707 }
708
709 params = (double *) malloc(sizeof(double) * NUM_PARAMS);
710 if (params == NULL)
711 {
712 ier = TRUE;
713 fclose(fid);
714 LOG_ERROR("Unable to malloc memory for parameters.");
715 return ier;
716 }
717 for (param_count = 0; param_count < NUM_PARAMS; ++param_count)
718 {
719 if (fgets(strbuf, 256, fid) == NULL)
720 {
721 if (ferror(fid))
722 {
723 /* A read error was encountered */
724 ier = TRUE;
725 fclose(fid);
726 LOG_ERROR("Error occurred while trying to read parameter file (ferror).");
727 return ier;
728 }
729 else if (feof(fid))
730 {
731 /* Reached the EOF mark before reading a character */
732 ier = TRUE;
733 fclose(fid);
734 LOG_ERROR("Error occurred while trying to read parameter file (reached EOF).");
735 return ier;
736 }
737 }
738
739 /* Cast string read in to double (this will throw out comments) */
740 params[param_count] = strtod(strbuf, NULL);
741
742 /* convert units */
743 if ((param_units[param_count][0] != 0.0)
744 || (param_units[param_count][1] != 0.0))
745 {
746 ier = KIM_ModelDriverCreate_ConvertUnit(*length_unit,
747 *energy_unit,
748 KIM_CHARGE_UNIT_unused,
749 KIM_TEMPERATURE_UNIT_unused,
750 KIM_TIME_UNIT_unused,
751 requestedLengthUnit,
752 requestedEnergyUnit,
753 KIM_CHARGE_UNIT_unused,
754 KIM_TEMPERATURE_UNIT_unused,
755 KIM_TIME_UNIT_unused,
756 param_units[param_count][0],
757 param_units[param_count][1],
758 0.0,
759 0.0,
760 0.0,
761 &conversion_factor);
762 if (ier == TRUE)
763 {
764 fclose(fid);
765 LOG_ERROR("Unable to convert units of parameter.");
766 return ier;
767 }
768 params[param_count] *= conversion_factor;
769 }
770 }
771 fclose(fid);
772
773 /* register species */
774 species_name = KIM_SpeciesName_FromString(species_name_string);
775 ier = KIM_ModelDriverCreate_SetSpeciesCode(
776 modelDriverCreate, species_name, SPECCODE);
777 if (ier == TRUE)
778 {
779 LOG_ERROR("Unable to set species code.");
780 return ier;
781 }
782
783 /* allocate buffer */
784 buffer = (struct model_buffer *) malloc(sizeof(struct model_buffer));
785 if (NULL == buffer)
786 {
787 ier = TRUE;
788 LOG_ERROR("Unable to malloc memory for buffer.");
789 return ier;
790 }
791
792 /* setup buffer */
793 buffer->params = params;
794 buffer->influenceDistance = get_influence_distance(buffer->params);
795 buffer->cutoff = buffer->influenceDistance;
796 buffer->cut_sq = buffer->cutoff * buffer->cutoff;
797 buffer->modelWillNotRequestNeighborsOfNoncontributingParticles = 1;
798 sprintf(buffer->species_name,
799 "%s", /* no buffer overrun possible, due to above code */
800 species_name_string);
801 /* end setup buffer */
802
803 /* store in model buffer */
804 KIM_ModelDriverCreate_SetModelBufferPointer(modelDriverCreate,
805 (void *) buffer);
806
807 /* publish model parameters */
808 for (param_count = 0; param_count < NUM_PARAMS; ++param_count)
809 {
810 ier = KIM_ModelDriverCreate_SetParameterPointerDouble(
811 modelDriverCreate,
812 1,
813 &(buffer->params[param_count]),
814 param_strings[param_count][0],
815 param_strings[param_count][1]);
816
817 if (ier == TRUE)
818 {
819 LOG_ERROR("Unable to set parameter pointer(s).");
820 return TRUE;
821 }
822 }
823
824 /* store model cutoff in KIM object */
825 KIM_ModelDriverCreate_SetInfluenceDistancePointer(
826 modelDriverCreate, &(buffer->influenceDistance));
827 KIM_ModelDriverCreate_SetNeighborListPointers(
828 modelDriverCreate,
829 1,
830 &(buffer->cutoff),
831 &(buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
832
833 return FALSE;
834 }
835
836 /* Refresh function */
837 #undef KIM_LOGGER_FUNCTION_NAME
838 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelRefresh_LogEntry
839 #undef KIM_LOGGER_OBJECT_NAME
840 #define KIM_LOGGER_OBJECT_NAME modelRefresh
refresh_routine(KIM_ModelRefresh * const modelRefresh)841 int refresh_routine(KIM_ModelRefresh * const modelRefresh)
842 {
843 struct model_buffer * buffer;
844
845 /* get model buffer from KIM object */
846 KIM_ModelRefresh_GetModelBufferPointer(modelRefresh, (void **) &buffer);
847
848 /* update cutoff and cutoff square */
849 buffer->influenceDistance = get_influence_distance(buffer->params);
850 buffer->cutoff = buffer->influenceDistance;
851 buffer->cut_sq = buffer->cutoff * buffer->cutoff;
852
853 /* store model cutoff in KIM object */
854 KIM_ModelRefresh_SetInfluenceDistancePointer(modelRefresh,
855 &(buffer->influenceDistance));
856 KIM_ModelRefresh_SetNeighborListPointers(
857 modelRefresh,
858 1,
859 &(buffer->cutoff),
860 &(buffer->modelWillNotRequestNeighborsOfNoncontributingParticles));
861
862 return FALSE;
863 }
864
865
866 /* destroy function */
destroy_routine(KIM_ModelDestroy * const modelDestroy)867 static int destroy_routine(KIM_ModelDestroy * const modelDestroy)
868 {
869 /* Local variables */
870 struct model_buffer * buffer;
871 int ier;
872
873 /* get model buffer from KIM object */
874 KIM_ModelDestroy_GetModelBufferPointer(modelDestroy, (void **) &buffer);
875
876 /* free the buffer */
877 free(buffer->params);
878 free(buffer);
879
880 ier = FALSE;
881 return ier;
882 }
883
884 /* compute arguments create routine */
885 #undef KIM_LOGGER_FUNCTION_NAME
886 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelComputeArgumentsCreate_LogEntry
887 #undef KIM_LOGGER_OBJECT_NAME
888 #define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate
compute_arguments_create(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)889 static int compute_arguments_create(
890 KIM_ModelCompute const * const modelCompute,
891 KIM_ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
892 {
893 int ier;
894
895 (void) modelCompute; /* avoid unused parameter warning */
896
897 /* register arguments */
898 ier = KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
899 modelComputeArgumentsCreate,
900 KIM_COMPUTE_ARGUMENT_NAME_partialEnergy,
901 KIM_SUPPORT_STATUS_optional)
902 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
903 modelComputeArgumentsCreate,
904 KIM_COMPUTE_ARGUMENT_NAME_partialParticleEnergy,
905 KIM_SUPPORT_STATUS_optional)
906 || KIM_ModelComputeArgumentsCreate_SetArgumentSupportStatus(
907 modelComputeArgumentsCreate,
908 KIM_COMPUTE_ARGUMENT_NAME_partialForces,
909 KIM_SUPPORT_STATUS_optional);
910 if (ier == TRUE)
911 {
912 LOG_ERROR("Unable to set argument supportStatus.");
913 return TRUE;
914 }
915 else
916 {
917 return FALSE;
918 }
919 }
920
921 /* compute arguments destroy routine */
compute_arguments_destroy(KIM_ModelCompute const * const modelCompute,KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)922 static int compute_arguments_destroy(
923 KIM_ModelCompute const * const modelCompute,
924 KIM_ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
925 {
926 (void) modelCompute; /* avoid unused parameter warning */
927 (void) modelComputeArgumentsDestroy;
928
929 /* Nothing further to do */
930
931 return FALSE;
932 }
933
934 /* write parameterized model routine */
935 #undef KIM_LOGGER_FUNCTION_NAME
936 #define KIM_LOGGER_FUNCTION_NAME KIM_ModelWriteParameterizedModel_LogEntry
937 #undef KIM_LOGGER_OBJECT_NAME
938 #define KIM_LOGGER_OBJECT_NAME modelWriteParameterizedModel
939 /* */
940 #define STR_LEN 2048
write_parameterized_model(KIM_ModelWriteParameterizedModel const * const modelWriteParameterizedModel)941 static int write_parameterized_model(
942 KIM_ModelWriteParameterizedModel const * const modelWriteParameterizedModel)
943 {
944 FILE * fp;
945 char string_buffer[STR_LEN];
946 struct model_buffer const * buffer;
947 char const * path;
948 char const * model_name;
949 int param_count;
950 int max_str_len;
951
952 /* get buffer from KIM object */
953 KIM_ModelWriteParameterizedModel_GetModelBufferPointer(
954 modelWriteParameterizedModel, (void **) &buffer);
955
956 KIM_ModelWriteParameterizedModel_GetPath(modelWriteParameterizedModel, &path);
957 KIM_ModelWriteParameterizedModel_GetModelName(modelWriteParameterizedModel,
958 &model_name);
959 max_str_len = strlen(path) + strlen(model_name) + 9;
960 if (max_str_len >= STR_LEN)
961 {
962 LOG_ERROR("Path and model name too long for internal buffers.");
963 return TRUE;
964 }
965
966 sprintf(string_buffer, "%s.params", model_name);
967 KIM_ModelWriteParameterizedModel_SetParameterFileName(
968 modelWriteParameterizedModel, string_buffer);
969 sprintf(string_buffer, "%s/%s.params", path, model_name);
970 fp = fopen(string_buffer, "w");
971 if (NULL == fp)
972 {
973 LOG_ERROR("Unable to open parameter file for write.");
974 return TRUE;
975 }
976
977 fprintf(fp, "%s\n", buffer->species_name);
978
979 for (param_count = 0; param_count < NUM_PARAMS; ++param_count)
980 { fprintf(fp, "%20.15f\n", buffer->params[param_count]); }
981
982 fclose(fp);
983
984 return FALSE;
985 }
986