kim-api  2.1.4-git+v2.1.3-git-3-g4c859c7f.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
LennardJones612Implementation.cpp
Go to the documentation of this file.
1 // CDDL HEADER START
2 //
3 // The contents of this file are subject to the terms of the Common Development
4 // and Distribution License Version 1.0 (the "License").
5 //
6 // You can obtain a copy of the license at
7 // http://www.opensource.org/licenses/CDDL-1.0. See the License for the
8 // specific language governing permissions and limitations under the License.
9 //
10 // When distributing Covered Code, include this CDDL HEADER in each file and
11 // include the License file in a prominent location with the name LICENSE.CDDL.
12 // If applicable, add the following below this CDDL HEADER, with the fields
13 // enclosed by brackets "[]" replaced with your own identifying information:
14 //
15 // Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
16 //
17 // CDDL HEADER END
18 //
19 
20 //
21 // Copyright (c) 2013--2015, Regents of the University of Minnesota.
22 // All rights reserved.
23 //
24 // Contributors:
25 // Ryan S. Elliott
26 // Stephen M. Whalen
27 // Andrew Akerson
28 //
29 
30 
31 #include <cmath>
32 #include <cstdlib>
33 #include <cstring>
34 #include <fstream>
35 #include <iostream>
36 #include <map>
37 #include <sstream>
38 
41 
42 #define MAXLINE 1024
43 #define IGNORE_RESULT(fn) \
44  if (fn) {}
45 
46 
47 //==============================================================================
48 //
49 // Implementation of LennardJones612Implementation public member functions
50 //
51 //==============================================================================
52 
53 //******************************************************************************
54 #undef KIM_LOGGER_OBJECT_NAME
55 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
56 //
58  KIM::ModelDriverCreate * const modelDriverCreate,
59  KIM::LengthUnit const requestedLengthUnit,
60  KIM::EnergyUnit const requestedEnergyUnit,
61  KIM::ChargeUnit const requestedChargeUnit,
62  KIM::TemperatureUnit const requestedTemperatureUnit,
63  KIM::TimeUnit const requestedTimeUnit,
64  int * const ier) :
65  numberModelSpecies_(0),
66  numberUniqueSpeciesPairs_(0),
67  shift_(0),
68  cutoffs_(NULL),
69  epsilons_(NULL),
70  sigmas_(NULL),
71  influenceDistance_(0.0),
72  cutoffsSq2D_(NULL),
73  modelWillNotRequestNeighborsOfNoncontributingParticles_(1),
74  fourEpsilonSigma6_2D_(NULL),
75  fourEpsilonSigma12_2D_(NULL),
76  twentyFourEpsilonSigma6_2D_(NULL),
77  fortyEightEpsilonSigma12_2D_(NULL),
78  oneSixtyEightEpsilonSigma6_2D_(NULL),
79  sixTwentyFourEpsilonSigma12_2D_(NULL),
80  shifts2D_(NULL),
81  cachedNumberOfParticles_(0)
82 {
83  FILE * parameterFilePointers[MAX_PARAMETER_FILES];
84  int numberParameterFiles;
85  modelDriverCreate->GetNumberOfParameterFiles(&numberParameterFiles);
86  *ier = OpenParameterFiles(
87  modelDriverCreate, numberParameterFiles, parameterFilePointers);
88  if (*ier) return;
89 
90  *ier = ProcessParameterFiles(
91  modelDriverCreate, numberParameterFiles, parameterFilePointers);
92  CloseParameterFiles(numberParameterFiles, parameterFilePointers);
93  if (*ier) return;
94 
95  *ier = ConvertUnits(modelDriverCreate,
96  requestedLengthUnit,
97  requestedEnergyUnit,
98  requestedChargeUnit,
99  requestedTemperatureUnit,
100  requestedTimeUnit);
101  if (*ier) return;
102 
103  *ier = SetRefreshMutableValues(modelDriverCreate);
104  if (*ier) return;
105 
106  *ier = RegisterKIMModelSettings(modelDriverCreate);
107  if (*ier) return;
108 
109  *ier = RegisterKIMParameters(modelDriverCreate);
110  if (*ier) return;
111 
112  *ier = RegisterKIMFunctions(modelDriverCreate);
113  if (*ier) return;
114 
115  // everything is good
116  *ier = false;
117  return;
118 }
119 
120 //******************************************************************************
122 { // note: it is ok to delete a null pointer and we have ensured that
123  // everything is initialized to null
124 
125  delete[] cutoffs_;
126  Deallocate2DArray(cutoffsSq2D_);
127  delete[] epsilons_;
128  delete[] sigmas_;
129  Deallocate2DArray(fourEpsilonSigma6_2D_);
130  Deallocate2DArray(fourEpsilonSigma12_2D_);
131  Deallocate2DArray(twentyFourEpsilonSigma6_2D_);
132  Deallocate2DArray(fortyEightEpsilonSigma12_2D_);
133  Deallocate2DArray(oneSixtyEightEpsilonSigma6_2D_);
134  Deallocate2DArray(sixTwentyFourEpsilonSigma12_2D_);
135  Deallocate2DArray(shifts2D_);
136 }
137 
138 //******************************************************************************
139 #undef KIM_LOGGER_OBJECT_NAME
140 #define KIM_LOGGER_OBJECT_NAME modelRefresh
141 //
143  KIM::ModelRefresh * const modelRefresh)
144 {
145  int ier;
146 
147  ier = SetRefreshMutableValues(modelRefresh);
148  if (ier) return ier;
149 
150  // nothing else to do for this case
151 
152  // everything is good
153  ier = false;
154  return ier;
155 }
156 
157 //******************************************************************************
159  KIM::ModelCompute const * const modelCompute,
160  KIM::ModelComputeArguments const * const modelComputeArguments)
161 {
162  int ier;
163 
164  // KIM API Model Input compute flags
165  bool isComputeProcess_dEdr = false;
166  bool isComputeProcess_d2Edr2 = false;
167  //
168  // KIM API Model Output compute flags
169  bool isComputeEnergy = false;
170  bool isComputeForces = false;
171  bool isComputeParticleEnergy = false;
172  bool isComputeVirial = false;
173  bool isComputeParticleVirial = false;
174  //
175  // KIM API Model Input
176  int const * particleSpeciesCodes = NULL;
177  int const * particleContributing = NULL;
178  VectorOfSizeDIM const * coordinates = NULL;
179  //
180  // KIM API Model Output
181  double * energy = NULL;
182  double * particleEnergy = NULL;
183  VectorOfSizeDIM * forces = NULL;
184  VectorOfSizeSix * virial = NULL;
185  VectorOfSizeSix * particleVirial = NULL;
186  ier = SetComputeMutableValues(modelComputeArguments,
187  isComputeProcess_dEdr,
188  isComputeProcess_d2Edr2,
189  isComputeEnergy,
190  isComputeForces,
191  isComputeParticleEnergy,
192  isComputeVirial,
193  isComputeParticleVirial,
194  particleSpeciesCodes,
195  particleContributing,
196  coordinates,
197  energy,
198  particleEnergy,
199  forces,
200  virial,
201  particleVirial);
202  if (ier) return ier;
203 
204  // Skip this check for efficiency
205  //
206  // ier = CheckParticleSpecies(modelComputeArguments, particleSpeciesCodes);
207  // if (ier) return ier;
208 
209  bool const isShift = (1 == shift_);
210 
212  return ier;
213 }
214 
215 //******************************************************************************
217  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
218 {
219  int ier;
220 
221  ier = RegisterKIMComputeArgumentsSettings(modelComputeArgumentsCreate);
222  if (ier) return ier;
223 
224  // nothing else to do for this case
225 
226  // everything is good
227  ier = false;
228  return ier;
229 }
230 
231 //******************************************************************************
234  /* modelComputeArgumentsDestroy */) const
235 {
236  int ier;
237 
238  // nothing else to do for this case
239 
240  // everything is good
241  ier = false;
242  return ier;
243 }
244 
245 //==============================================================================
246 //
247 // Implementation of LennardJones612Implementation private member functions
248 //
249 //==============================================================================
250 
251 //******************************************************************************
252 void LennardJones612Implementation::AllocatePrivateParameterMemory()
253 {
254  // nothing to do for this case
255 }
256 
257 //******************************************************************************
258 void LennardJones612Implementation::AllocateParameterMemory()
259 { // allocate memory for data
260  cutoffs_ = new double[numberUniqueSpeciesPairs_];
262  cutoffsSq2D_, numberModelSpecies_, numberModelSpecies_);
263 
264  epsilons_ = new double[numberUniqueSpeciesPairs_];
265  sigmas_ = new double[numberUniqueSpeciesPairs_];
267  fourEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
269  fourEpsilonSigma12_2D_, numberModelSpecies_, numberModelSpecies_);
271  twentyFourEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
273  fortyEightEpsilonSigma12_2D_, numberModelSpecies_, numberModelSpecies_);
275  oneSixtyEightEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
276  AllocateAndInitialize2DArray(sixTwentyFourEpsilonSigma12_2D_,
277  numberModelSpecies_,
278  numberModelSpecies_);
279 
281  shifts2D_, numberModelSpecies_, numberModelSpecies_);
282 }
283 
284 //******************************************************************************
285 #undef KIM_LOGGER_OBJECT_NAME
286 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
287 //
288 int LennardJones612Implementation::OpenParameterFiles(
289  KIM::ModelDriverCreate * const modelDriverCreate,
290  int const numberParameterFiles,
291  FILE * parameterFilePointers[MAX_PARAMETER_FILES])
292 {
293  int ier;
294 
295  if (numberParameterFiles > MAX_PARAMETER_FILES)
296  {
297  ier = true;
298  LOG_ERROR("LennardJones612 given too many parameter files");
299  return ier;
300  }
301 
302  for (int i = 0; i < numberParameterFiles; ++i)
303  {
304  std::string const * paramFileName;
305  ier = modelDriverCreate->GetParameterFileName(i, &paramFileName);
306  if (ier)
307  {
308  LOG_ERROR("Unable to get parameter file name");
309  return ier;
310  }
311  parameterFilePointers[i] = fopen(paramFileName->c_str(), "r");
312  if (parameterFilePointers[i] == 0)
313  {
314  char message[MAXLINE];
315  sprintf(message,
316  "LennardJones612 parameter file number %d cannot be opened",
317  i);
318  ier = true;
319  LOG_ERROR(message);
320  for (int j = i - 1; j >= 0; --j) { fclose(parameterFilePointers[j]); }
321  return ier;
322  }
323  }
324 
325  // everything is good
326  ier = false;
327  return ier;
328 }
329 
330 //******************************************************************************
331 #undef KIM_LOGGER_OBJECT_NAME
332 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
333 //
334 int LennardJones612Implementation::ProcessParameterFiles(
335  KIM::ModelDriverCreate * const modelDriverCreate,
336  int const /* numberParameterFiles */,
337  FILE * const parameterFilePointers[MAX_PARAMETER_FILES])
338 {
339  int N, ier;
340  int endOfFileFlag = 0;
341  char spec1[MAXLINE], spec2[MAXLINE], nextLine[MAXLINE];
342  char * nextLinePtr;
343  int iIndex, jIndex, indx, iiIndex, jjIndex;
344  double nextCutoff, nextEpsilon, nextSigma;
345 
346  nextLinePtr = nextLine;
347 
348  getNextDataLine(
349  parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
350  ier = sscanf(nextLine, "%d %d", &N, &shift_);
351  if (ier != 2)
352  {
353  sprintf(nextLine, "unable to read first line of the parameter file");
354  ier = true;
355  LOG_ERROR(nextLine);
356  fclose(parameterFilePointers[0]);
357  return ier;
358  }
359  numberModelSpecies_ = N;
360  numberUniqueSpeciesPairs_
361  = ((numberModelSpecies_ + 1) * numberModelSpecies_) / 2;
362  AllocateParameterMemory();
363 
364  // set all values in the arrays to -1 for mixing later
365  for (int i = 0; i < ((N + 1) * N / 2); i++)
366  {
367  cutoffs_[i] = -1;
368  epsilons_[i] = -1;
369  sigmas_[i] = -1;
370  }
371 
372 
373  // keep track of known species
374  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
375  modelSpeciesMap;
376  std::vector<KIM::SpeciesName> speciesNameVector;
377  int index = 0;
378 
379  // Read and process data lines
380  getNextDataLine(
381  parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
382  while (endOfFileFlag == 0)
383  {
384  ier = sscanf(nextLine,
385  "%s %s %lg %lg %lg",
386  spec1,
387  spec2,
388  &nextCutoff,
389  &nextEpsilon,
390  &nextSigma);
391  if (ier != 5)
392  {
393  sprintf(nextLine, "error reading lines of the parameter file");
394  LOG_ERROR(nextLine);
395  return true;
396  }
397 
398  // convert species strings to proper type instances
399  KIM::SpeciesName const specName1(spec1);
400  KIM::SpeciesName const specName2(spec2);
401 
402  // check for new species
403  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
404  const_iterator iIter
405  = modelSpeciesMap.find(specName1);
406  if (iIter == modelSpeciesMap.end())
407  {
408  modelSpeciesMap[specName1] = index;
409  modelSpeciesCodeList_.push_back(index);
410  speciesNameVector.push_back(specName1);
411 
412  ier = modelDriverCreate->SetSpeciesCode(specName1, index);
413  if (ier) return ier;
414  iIndex = index;
415  index++;
416  }
417  else
418  {
419  iIndex = modelSpeciesMap[specName1];
420  }
421  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
422  const_iterator jIter
423  = modelSpeciesMap.find(specName2);
424  if (jIter == modelSpeciesMap.end())
425  {
426  modelSpeciesMap[specName2] = index;
427  modelSpeciesCodeList_.push_back(index);
428  speciesNameVector.push_back(specName2);
429 
430  ier = modelDriverCreate->SetSpeciesCode(specName2, index);
431  if (ier) return ier;
432  jIndex = index;
433  index++;
434  }
435  else
436  {
437  jIndex = modelSpeciesMap[specName2];
438  }
439 
440  if (iIndex >= jIndex)
441  { indx = jIndex * N + iIndex - (jIndex * jIndex + jIndex) / 2; }
442  else
443  {
444  indx = iIndex * N + jIndex - (iIndex * iIndex + iIndex) / 2;
445  }
446  cutoffs_[indx] = nextCutoff;
447  epsilons_[indx] = nextEpsilon;
448  sigmas_[indx] = nextSigma;
449 
450  getNextDataLine(
451  parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
452  }
453 
454  // check that we got all like - like pairs
455  std::stringstream ss;
456  ss << "There are not values for like-like pairs of:";
457  for (int i = 0; i < N; i++)
458  {
459  if (cutoffs_[(i * N + i - (i * i + i) / 2)] == -1)
460  {
461  ss << " ";
462  ss << (speciesNameVector[i].ToString()).c_str();
463  ier = -1;
464  }
465  }
466  if (ier == -1)
467  {
468  LOG_ERROR(ss);
469  return true;
470  }
471 
472  // Perform Mixing if nessisary
473  for (int jIndex = 0; jIndex < N; jIndex++)
474  {
475  jjIndex = (jIndex * N + jIndex - (jIndex * jIndex + jIndex) / 2);
476  for (int iIndex = (jIndex + 1); iIndex < N; iIndex++)
477  {
478  indx = jIndex * N + iIndex - (jIndex * jIndex + jIndex) / 2;
479  if (cutoffs_[indx] == -1)
480  {
481  iiIndex = (iIndex * N + iIndex - (iIndex * iIndex + iIndex) / 2);
482  epsilons_[indx] = sqrt(epsilons_[iiIndex] * epsilons_[jjIndex]);
483  sigmas_[indx] = (sigmas_[iiIndex] + sigmas_[jjIndex]) / 2.0;
484  cutoffs_[indx] = (cutoffs_[iiIndex] + cutoffs_[jjIndex]) / 2.0;
485  }
486  }
487  }
488 
489  // everything is good
490  ier = false;
491  return ier;
492 }
493 
494 //******************************************************************************
495 void LennardJones612Implementation::getNextDataLine(FILE * const filePtr,
496  char * nextLinePtr,
497  int const maxSize,
498  int * endOfFileFlag)
499 {
500  do
501  {
502  if (fgets(nextLinePtr, maxSize, filePtr) == NULL)
503  {
504  *endOfFileFlag = 1;
505  break;
506  }
507  while ((nextLinePtr[0] == ' ' || nextLinePtr[0] == '\t')
508  || (nextLinePtr[0] == '\n' || nextLinePtr[0] == '\r'))
509  { nextLinePtr = (nextLinePtr + 1); }
510  } while ((strncmp("#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
511 }
512 
513 //******************************************************************************
514 void LennardJones612Implementation::CloseParameterFiles(
515  int const numberParameterFiles,
516  FILE * const parameterFilePointers[MAX_PARAMETER_FILES])
517 {
518  for (int i = 0; i < numberParameterFiles; ++i)
519  fclose(parameterFilePointers[i]);
520 }
521 
522 //******************************************************************************
523 #undef KIM_LOGGER_OBJECT_NAME
524 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
525 //
526 int LennardJones612Implementation::ConvertUnits(
527  KIM::ModelDriverCreate * const modelDriverCreate,
528  KIM::LengthUnit const requestedLengthUnit,
529  KIM::EnergyUnit const requestedEnergyUnit,
530  KIM::ChargeUnit const requestedChargeUnit,
531  KIM::TemperatureUnit const requestedTemperatureUnit,
532  KIM::TimeUnit const requestedTimeUnit)
533 {
534  int ier;
535 
536  // define default base units
542 
543  // changing units of cutoffs and sigmas
544  double convertLength = 1.0;
545  ier = KIM::ModelDriverCreate::ConvertUnit(fromLength,
546  fromEnergy,
547  fromCharge,
548  fromTemperature,
549  fromTime,
550  requestedLengthUnit,
551  requestedEnergyUnit,
552  requestedChargeUnit,
553  requestedTemperatureUnit,
554  requestedTimeUnit,
555  1.0,
556  0.0,
557  0.0,
558  0.0,
559  0.0,
560  &convertLength);
561  if (ier)
562  {
563  LOG_ERROR("Unable to convert length unit");
564  return ier;
565  }
566  if (convertLength != ONE)
567  {
568  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
569  {
570  cutoffs_[i] *= convertLength; // convert to active units
571  sigmas_[i] *= convertLength; // convert to active units
572  }
573  }
574  // changing units of epsilons
575  double convertEnergy = 1.0;
576  ier = KIM::ModelDriverCreate::ConvertUnit(fromLength,
577  fromEnergy,
578  fromCharge,
579  fromTemperature,
580  fromTime,
581  requestedLengthUnit,
582  requestedEnergyUnit,
583  requestedChargeUnit,
584  requestedTemperatureUnit,
585  requestedTimeUnit,
586  0.0,
587  1.0,
588  0.0,
589  0.0,
590  0.0,
591  &convertEnergy);
592  if (ier)
593  {
594  LOG_ERROR("Unable to convert energy unit");
595  return ier;
596  }
597  if (convertEnergy != ONE)
598  {
599  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
600  {
601  epsilons_[i] *= convertEnergy; // convert to active units
602  }
603  }
604 
605  // register units
606  ier = modelDriverCreate->SetUnits(requestedLengthUnit,
607  requestedEnergyUnit,
611  if (ier)
612  {
613  LOG_ERROR("Unable to set units to requested values");
614  return ier;
615  }
616 
617  // everything is good
618  ier = false;
619  return ier;
620 }
621 
622 //******************************************************************************
623 int LennardJones612Implementation::RegisterKIMModelSettings(
624  KIM::ModelDriverCreate * const modelDriverCreate) const
625 {
626  // register numbering
627  int error = modelDriverCreate->SetModelNumbering(KIM::NUMBERING::zeroBased);
628 
629  return error;
630 }
631 
632 //******************************************************************************
633 #undef KIM_LOGGER_OBJECT_NAME
634 #define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate
635 //
636 int LennardJones612Implementation::RegisterKIMComputeArgumentsSettings(
637  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
638 {
639  // register arguments
640  LOG_INFORMATION("Register argument supportStatus");
641  int error = modelComputeArgumentsCreate->SetArgumentSupportStatus(
644  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
647  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
650  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
653  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
656 
657 
658  // register callbacks
659  LOG_INFORMATION("Register callback supportStatus");
660  error = error
661  || modelComputeArgumentsCreate->SetCallbackSupportStatus(
664  || modelComputeArgumentsCreate->SetCallbackSupportStatus(
667 
668  return error;
669 }
670 
671 //******************************************************************************
672 // helper macro
673 #define SNUM(x) \
674  static_cast<std::ostringstream const &>(std::ostringstream() \
675  << std::dec << x) \
676  .str()
677 //******************************************************************************
678 #undef KIM_LOGGER_OBJECT_NAME
679 #define KIM_LOGGER_OBJECT_NAME modelDriverCreate
680 //
681 int LennardJones612Implementation::RegisterKIMParameters(
682  KIM::ModelDriverCreate * const modelDriverCreate)
683 {
684  int ier = false;
685 
686  // publish parameters (order is important)
687  ier = modelDriverCreate->SetParameterPointer(
688  1,
689  &shift_,
690  "shift",
691  "If (shift == 1), all LJ potentials are shifted to zero energy "
692  "at their respective cutoff distance. Otherwise, no shifting is "
693  "performed.");
694  if (ier)
695  {
696  LOG_ERROR("set_parameter shift");
697  return ier;
698  }
699 
700  ier = modelDriverCreate->SetParameterPointer(
701  numberUniqueSpeciesPairs_,
702  cutoffs_,
703  "cutoffs",
704  "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
705  + ") "
706  "in row-major storage. Ordering is according to SpeciesCode "
707  "values. "
708  "For example, to find the parameter related to SpeciesCode 'i' and "
709  "SpeciesCode 'j' (i >= j), use (zero-based) "
710  "index = (j*N + i - (j*j + j)/2).");
711  if (ier)
712  {
713  LOG_ERROR("set_parameter cutoffs");
714  return ier;
715  }
716  ier = modelDriverCreate->SetParameterPointer(
717  numberUniqueSpeciesPairs_,
718  epsilons_,
719  "epsilons",
720  "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
721  + ") "
722  "in row-major storage. Ordering is according to SpeciesCode "
723  "values. "
724  "For example, to find the parameter related to SpeciesCode 'i' and "
725  "SpeciesCode 'j' (i >= j), use (zero-based) "
726  "index = (j*N + i - (j*j + j)/2).");
727  if (ier)
728  {
729  LOG_ERROR("set_parameter epsilons");
730  return ier;
731  }
732  ier = modelDriverCreate->SetParameterPointer(
733  numberUniqueSpeciesPairs_,
734  sigmas_,
735  "sigmas",
736  "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
737  + ") "
738  "in row-major storage. Ordering is according to SpeciesCode "
739  "values. "
740  "For example, to find the parameter related to SpeciesCode 'i' and "
741  "SpeciesCode 'j' (i >= j), use (zero-based) "
742  "index = (j*N + i - (j*j + j)/2).");
743  if (ier)
744  {
745  LOG_ERROR("set_parameter sigmas");
746  return ier;
747  }
748 
749  // everything is good
750  ier = false;
751  return ier;
752 }
753 
754 //******************************************************************************
755 int LennardJones612Implementation::RegisterKIMFunctions(
756  KIM::ModelDriverCreate * const modelDriverCreate) const
757 {
758  int error;
759 
760  // Use function pointer definitions to verify correct prototypes
768 
769  // register the destroy() and reinit() functions
770  error = modelDriverCreate->SetRoutinePointer(
773  true,
774  reinterpret_cast<KIM::Function *>(destroy))
775  || modelDriverCreate->SetRoutinePointer(
778  true,
779  reinterpret_cast<KIM::Function *>(refresh))
780  || modelDriverCreate->SetRoutinePointer(
783  true,
784  reinterpret_cast<KIM::Function *>(compute))
785  || modelDriverCreate->SetRoutinePointer(
788  true,
789  reinterpret_cast<KIM::Function *>(CACreate))
790  || modelDriverCreate->SetRoutinePointer(
793  true,
794  reinterpret_cast<KIM::Function *>(CADestroy));
795  return error;
796 }
797 
798 //******************************************************************************
799 template<class ModelObj>
800 int LennardJones612Implementation::SetRefreshMutableValues(
801  ModelObj * const modelObj)
802 { // use (possibly) new values of parameters to compute other quantities
803  // NOTE: This function is templated because it's called with both a
804  // modelDriverCreate object during initialization and with a
805  // modelRefresh object when the Model's parameters have been altered
806  int ier;
807 
808  // update cutoffsSq, epsilons, and sigmas
809  for (int i = 0; i < numberModelSpecies_; ++i)
810  {
811  for (int j = 0; j <= i; ++j)
812  {
813  int const index = j * numberModelSpecies_ + i - (j * j + j) / 2;
814  cutoffsSq2D_[i][j] = cutoffsSq2D_[j][i]
815  = (cutoffs_[index] * cutoffs_[index]);
816  fourEpsilonSigma6_2D_[i][j] = fourEpsilonSigma6_2D_[j][i]
817  = 4.0 * epsilons_[index] * pow(sigmas_[index], 6.0);
818  fourEpsilonSigma12_2D_[i][j] = fourEpsilonSigma12_2D_[j][i]
819  = 4.0 * epsilons_[index] * pow(sigmas_[index], 12.0);
820  twentyFourEpsilonSigma6_2D_[i][j] = twentyFourEpsilonSigma6_2D_[j][i]
821  = 6.0 * fourEpsilonSigma6_2D_[i][j];
822  fortyEightEpsilonSigma12_2D_[i][j] = fortyEightEpsilonSigma12_2D_[j][i]
823  = 12.0 * fourEpsilonSigma12_2D_[i][j];
824  oneSixtyEightEpsilonSigma6_2D_[i][j]
825  = oneSixtyEightEpsilonSigma6_2D_[j][i]
826  = 7.0 * twentyFourEpsilonSigma6_2D_[i][j];
827  sixTwentyFourEpsilonSigma12_2D_[i][j]
828  = sixTwentyFourEpsilonSigma12_2D_[j][i]
829  = 13.0 * fortyEightEpsilonSigma12_2D_[i][j];
830  }
831  }
832 
833  // update cutoff value in KIM API object
834  influenceDistance_ = 0.0;
835 
836  for (int i = 0; i < numberModelSpecies_; i++)
837  {
838  int indexI = modelSpeciesCodeList_[i];
839 
840  for (int j = 0; j < numberModelSpecies_; j++)
841  {
842  int indexJ = modelSpeciesCodeList_[j];
843 
844  if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
845  { influenceDistance_ = cutoffsSq2D_[indexI][indexJ]; }
846  }
847  }
848 
849  influenceDistance_ = sqrt(influenceDistance_);
850  modelObj->SetInfluenceDistancePointer(&influenceDistance_);
851  modelObj->SetNeighborListPointers(
852  1,
853  &influenceDistance_,
854  &modelWillNotRequestNeighborsOfNoncontributingParticles_);
855 
856  // update shifts
857  // compute and set shifts2D_ check if minus sign
858  double const * const * const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
859  double const * const * const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
860  if (1 == shift_)
861  {
862  double phi;
863  for (int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
864  {
865  for (int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
866  {
867  int const index = jSpecies * numberModelSpecies_ + iSpecies
868  - (jSpecies * jSpecies + jSpecies) / 2;
869  double const rij2 = cutoffs_[index] * cutoffs_[index];
870  double const r2iv = 1.0 / rij2;
871  double const r6iv = r2iv * r2iv * r2iv;
873  shifts2D_[iSpecies][jSpecies] = shifts2D_[jSpecies][iSpecies] = phi;
874  }
875  }
876  }
877 
878  // everything is good
879  ier = false;
880  return ier;
881 }
882 
883 //******************************************************************************
884 #undef KIM_LOGGER_OBJECT_NAME
885 #define KIM_LOGGER_OBJECT_NAME modelComputeArguments
886 //
887 int LennardJones612Implementation::SetComputeMutableValues(
888  KIM::ModelComputeArguments const * const modelComputeArguments,
889  bool & isComputeProcess_dEdr,
890  bool & isComputeProcess_d2Edr2,
891  bool & isComputeEnergy,
892  bool & isComputeForces,
893  bool & isComputeParticleEnergy,
894  bool & isComputeVirial,
895  bool & isComputeParticleVirial,
896  int const *& particleSpeciesCodes,
897  int const *& particleContributing,
898  VectorOfSizeDIM const *& coordinates,
899  double *& energy,
900  double *& particleEnergy,
901  VectorOfSizeDIM *& forces,
902  VectorOfSizeSix *& virial,
903  VectorOfSizeSix *& particleVirial)
904 {
905  int ier = true;
906 
907  // get compute flags
908  int compProcess_dEdr;
909  int compProcess_d2Edr2;
910 
911  modelComputeArguments->IsCallbackPresent(
913  modelComputeArguments->IsCallbackPresent(
914  KIM::COMPUTE_CALLBACK_NAME::ProcessD2EDr2Term, &compProcess_d2Edr2);
915 
916  isComputeProcess_dEdr = compProcess_dEdr;
917  isComputeProcess_d2Edr2 = compProcess_d2Edr2;
918 
919  int const * numberOfParticles;
920  ier = modelComputeArguments->GetArgumentPointer(
922  || modelComputeArguments->GetArgumentPointer(
925  || modelComputeArguments->GetArgumentPointer(
928  || modelComputeArguments->GetArgumentPointer(
930  (double const **) &coordinates)
931  || modelComputeArguments->GetArgumentPointer(
933  || modelComputeArguments->GetArgumentPointer(
935  || modelComputeArguments->GetArgumentPointer(
937  (double const **) &forces)
938  || modelComputeArguments->GetArgumentPointer(
940  (double const **) &virial)
941  || modelComputeArguments->GetArgumentPointer(
943  (double const **) &particleVirial);
944  if (ier)
945  {
946  LOG_ERROR("GetArgumentPointer");
947  return ier;
948  }
949 
950  isComputeEnergy = (energy != NULL);
951  isComputeParticleEnergy = (particleEnergy != NULL);
952  isComputeForces = (forces != NULL);
953  isComputeVirial = (virial != NULL);
954  isComputeParticleVirial = (particleVirial != NULL);
955 
956  // update values
957  cachedNumberOfParticles_ = *numberOfParticles;
958 
959  // everything is good
960  ier = false;
961  return ier;
962 }
963 
964 //******************************************************************************
965 #undef KIM_LOGGER_OBJECT_NAME
966 #define KIM_LOGGER_OBJECT_NAME modelCompute
967 int LennardJones612Implementation::CheckParticleSpeciesCodes(
968  KIM::ModelCompute const * const modelCompute,
969  int const * const particleSpeciesCodes) const
970 {
971  int ier;
972  for (int i = 0; i < cachedNumberOfParticles_; ++i)
973  {
974  if ((particleSpeciesCodes[i] < 0)
975  || (particleSpeciesCodes[i] >= numberModelSpecies_))
976  {
977  ier = true;
978  LOG_ERROR("unsupported particle species codes detected");
979  return ier;
980  }
981  }
982 
983  // everything is good
984  ier = false;
985  return ier;
986 }
987 
988 //******************************************************************************
989 int LennardJones612Implementation::GetComputeIndex(
990  const bool & isComputeProcess_dEdr,
991  const bool & isComputeProcess_d2Edr2,
992  const bool & isComputeEnergy,
993  const bool & isComputeForces,
994  const bool & isComputeParticleEnergy,
995  const bool & isComputeVirial,
996  const bool & isComputeParticleVirial,
997  const bool & isShift) const
998 {
999  // const int processdE = 2;
1000  const int processd2E = 2;
1001  const int energy = 2;
1002  const int force = 2;
1003  const int particleEnergy = 2;
1004  const int virial = 2;
1005  const int particleVirial = 2;
1006  const int shift = 2;
1007 
1008 
1009  int index = 0;
1010 
1011  // processdE
1012  index += (int(isComputeProcess_dEdr)) * processd2E * energy * force
1013  * particleEnergy * virial * particleVirial * shift;
1014 
1015  // processd2E
1016  index += (int(isComputeProcess_d2Edr2)) * energy * force * particleEnergy
1017  * virial * particleVirial * shift;
1018 
1019  // energy
1020  index += (int(isComputeEnergy)) * force * particleEnergy * virial
1021  * particleVirial * shift;
1022 
1023  // force
1024  index += (int(isComputeForces)) * particleEnergy * virial * particleVirial
1025  * shift;
1026 
1027  // particleEnergy
1028  index += (int(isComputeParticleEnergy)) * virial * particleVirial * shift;
1029 
1030  // virial
1031  index += (int(isComputeVirial)) * particleVirial * shift;
1032 
1033  // particleVirial
1034  index += (int(isComputeParticleVirial)) * shift;
1035 
1036  // shift
1037  index += (int(isShift));
1038 
1039  return index;
1040 }
1041 
1042 //******************************************************************************
1043 void LennardJones612Implementation::ProcessVirialTerm(
1044  const double & dEidr,
1045  const double & rij,
1046  const double * const r_ij,
1047  const int & /* i */,
1048  const int & /* j */,
1049  VectorOfSizeSix virial) const
1050 {
1051  double const v = dEidr / rij;
1052 
1053  virial[0] += v * r_ij[0] * r_ij[0];
1054  virial[1] += v * r_ij[1] * r_ij[1];
1055  virial[2] += v * r_ij[2] * r_ij[2];
1056  virial[3] += v * r_ij[1] * r_ij[2];
1057  virial[4] += v * r_ij[0] * r_ij[2];
1058  virial[5] += v * r_ij[0] * r_ij[1];
1059 }
1060 
1061 //******************************************************************************
1062 void LennardJones612Implementation::ProcessParticleVirialTerm(
1063  const double & dEidr,
1064  const double & rij,
1065  const double * const r_ij,
1066  const int & i,
1067  const int & j,
1068  VectorOfSizeSix * const particleVirial) const
1069 {
1070  double const v = dEidr / rij;
1071  VectorOfSizeSix vir;
1072 
1073  vir[0] = 0.5 * v * r_ij[0] * r_ij[0];
1074  vir[1] = 0.5 * v * r_ij[1] * r_ij[1];
1075  vir[2] = 0.5 * v * r_ij[2] * r_ij[2];
1076  vir[3] = 0.5 * v * r_ij[1] * r_ij[2];
1077  vir[4] = 0.5 * v * r_ij[0] * r_ij[2];
1078  vir[5] = 0.5 * v * r_ij[0] * r_ij[1];
1079 
1080  for (int k = 0; k < 6; ++k)
1081  {
1082  particleVirial[i][k] += vir[k];
1083  particleVirial[j][k] += vir[k];
1084  }
1085 }
1086 
1087 //==============================================================================
1088 //
1089 // Implementation of helper functions
1090 //
1091 //==============================================================================
1092 
1093 //******************************************************************************
1094 void AllocateAndInitialize2DArray(double **& arrayPtr,
1095  int const extentZero,
1096  int const extentOne)
1097 { // allocate memory and set pointers
1098  arrayPtr = new double *[extentZero];
1099  arrayPtr[0] = new double[extentZero * extentOne];
1100  for (int i = 1; i < extentZero; ++i)
1101  { arrayPtr[i] = arrayPtr[i - 1] + extentOne; }
1102 
1103  // initialize
1104  for (int i = 0; i < extentZero; ++i)
1105  {
1106  for (int j = 0; j < extentOne; ++j) { arrayPtr[i][j] = 0.0; }
1107  }
1108 }
1109 
1110 //******************************************************************************
1111 void Deallocate2DArray(double **& arrayPtr)
1112 { // deallocate memory
1113  if (arrayPtr != NULL) delete[] arrayPtr[0];
1114  delete[] arrayPtr;
1115 
1116  // nullify pointer
1117  arrayPtr = NULL;
1118 }
int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
int ModelComputeArgumentsDestroyFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsDestroy routine.
TemperatureUnit const unused
Indicates that a TemperatureUnit is not used.
ComputeArgumentName const partialParticleEnergy
The standard partialParticleEnergy argument.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::R...
An Extensible Enumeration for the TemperatureUnit&#39;s supported by the KIM API.
static int Refresh(KIM::ModelRefresh *const modelRefresh)
An Extensible Enumeration for the TimeUnit&#39;s supported by the KIM API.
int ModelDestroyFunction(ModelDestroy *const modelDestroy)
Prototype for MODEL_ROUTINE_NAME::Destroy routine.
TimeUnit const ps
The standard picosecond unit of time.
recursive subroutine, public destroy(model_destroy_handle, ierr)
#define MAX_PARAMETER_FILES
ComputeArgumentName const coordinates
The standard coordinates argument.
#define LOG_ERROR(message)
Convenience macro for ERROR Log entries with compile-time optimization.
ModelRoutineName const Destroy
The standard Destroy routine.
ComputeArgumentName const partialVirial
The standard partialVirial argument.
ChargeUnit const unused
Indicates that a ChargeUnit is not used.
ChargeUnit const e
The standard electron unit of charge.
recursive subroutine, public refresh(model_refresh_handle, ierr)
void GetNumberOfParameterFiles(int *const numberOfParameterFiles) const
Get the number of parameter files provided by the parameterized model.
int SetArgumentSupportStatus(ComputeArgumentName const computeArgumentName, SupportStatus const supportStatus)
Set the SupportStatus of a ComputeArgumentName.
int ModelComputeFunction(ModelCompute const *const modelCompute, ModelComputeArguments const *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::Compute routine.
ComputeCallbackName const ProcessD2EDr2Term
The standard ProcessD2EDr2Term callback.
double VectorOfSizeDIM[DIMENSION]
An Extensible Enumeration for the LengthUnit&#39;s supported by the KIM API.
int GetArgumentPointer(ComputeArgumentName const computeArgumentName, int const **const ptr) const
Get the data pointer for a ComputeArgumentName.
int ModelComputeArgumentsCreateFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsCreate routine.
int SetModelNumbering(Numbering const numbering)
Set the Model&#39;s particle Numbering.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
ComputeArgumentName const particleContributing
The standard particleContributing argument.
static int ComputeArgumentsCreate(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
ComputeArgumentName const partialEnergy
The standard partialEnergy argument.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
ModelRoutineName const Refresh
The standard Refresh routine.
LennardJones612Implementation(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 *const ier)
ModelRoutineName const ComputeArgumentsCreate
The standard ComputeArgumentsCreate routine.
ModelRoutineName const ComputeArgumentsDestroy
The standard ComputeArgumentsDestroy routine.
int GetParameterFileName(int const index, std::string const **const parameterFileName) const
Get a particular parameter file name.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
static int Destroy(KIM::ModelDestroy *const modelDestroy)
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
Set the Model&#39;s base unit values.
An Extensible Enumeration for the EnergyUnit&#39;s supported by the KIM API.
LanguageName const cpp
The standard cpp language.
ModelRoutineName const Compute
The standard Compute routine.
ComputeArgumentName const partialParticleVirial
The standard partialParticleVirial argument.
ComputeCallbackName const ProcessDEDrTerm
The standard ProcessDEDrTerm callback.
int SetRoutinePointer(ModelRoutineName const modelRoutineName, LanguageName const languageName, int const required, Function *const fptr)
Set the function pointer for the ModelRoutineName of interest.
int Refresh(KIM::ModelRefresh *const modelRefresh)
LengthUnit const A
The standard angstrom unit of length.
int ModelRefreshFunction(ModelRefresh *const modelRefresh)
Prototype for MODEL_ROUTINE_NAME::Refresh routine.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
ComputeArgumentName const partialForces
The standard partialForces argument.
int SetParameterPointer(int const extent, int *const ptr, std::string const &name, std::string const &description)
Set the next parameter data pointer to be provided by the model.
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
EnergyUnit const eV
The standard electronvolt unit of energy.
static int ComputeArgumentsDestroy(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
ComputeArgumentName const particleSpeciesCodes
The standard particleSpeciesCodes argument.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
#define LENNARD_JONES_PHI(exshift)
An Extensible Enumeration for the ChargeUnit&#39;s supported by the KIM API.
ComputeArgumentName const numberOfParticles
The standard numberOfParticles argument.
SpeciesName const N
The standard Nitrogen species.
double VectorOfSizeSix[6]
#define LOG_INFORMATION(message)
Convenience macro for INFORMATION Log entries with compile-time optimization.
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
void Deallocate2DArray(double **&arrayPtr)
int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
Numbering const zeroBased
The standard zeroBased numbering.
static int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
int SetSpeciesCode(SpeciesName const speciesName, int const code)
Set integer code for supported SpeciesName.
TimeUnit const unused
Indicates that a TimeUnit is not used.
An Extensible Enumeration for the SpeciesName&#39;s supported by the KIM API.
int SetCallbackSupportStatus(ComputeCallbackName const computeCallbackName, SupportStatus const supportStatus)
Set the SupportStatus of a ComputeCallbackName.
LogVerbosity const error
The standard error verbosity.
int IsCallbackPresent(ComputeCallbackName const computeCallbackName, int *const present) const
Determine if the Simulator has provided a non-NULL function pointer for a ComputeCallbackName of inte...
static int ConvertUnit(LengthUnit const fromLengthUnit, EnergyUnit const fromEnergyUnit, ChargeUnit const fromChargeUnit, TemperatureUnit const fromTemperatureUnit, TimeUnit const fromTimeUnit, LengthUnit const toLengthUnit, EnergyUnit const toEnergyUnit, ChargeUnit const toChargeUnit, TemperatureUnit const toTemperatureUnit, TimeUnit const toTimeUnit, double const lengthExponent, double const energyExponent, double const chargeExponent, double const temperatureExponent, double const timeExponent, double *const conversionFactor)
Get the multiplicative factor to convert between a derived unit represented in two different sets of ...
SupportStatus const optional
The standard optional status.
TemperatureUnit const K
The standard Kelvin unit of temperature.