KIM API V2
LennardJones612Implementation.cpp
Go to the documentation of this file.
1 //
2 // CDDL HEADER START
3 //
4 // The contents of this file are subject to the terms of the Common Development
5 // 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 LICENSE.CDDL.
13 // If applicable, add the following below this CDDL HEADER, with the fields
14 // enclosed by brackets "[]" replaced with your own identifying information:
15 //
16 // Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
17 //
18 // CDDL HEADER END
19 //
20 
21 //
22 // Copyright (c) 2013--2015, Regents of the University of Minnesota.
23 // All rights reserved.
24 //
25 // Contributors:
26 // Ryan S. Elliott
27 // Stephen M. Whalen
28 // Andrew Akerson
29 //
30 
31 
32 #include <cmath>
33 #include <cstdlib>
34 #include <cstring>
35 #include <fstream>
36 #include <iostream>
37 #include <map>
38 
40 #include "KIM_Numbering.hpp"
41 #include "KIM_LanguageName.hpp"
42 #include "KIM_SpeciesName.hpp"
43 #include "KIM_SupportStatus.hpp"
44 #include "KIM_ArgumentName.hpp"
45 #include "KIM_CallbackName.hpp"
46 
47 #define MAXLINE 1024
48 #define IGNORE_RESULT(fn) if(fn){}
49 
50 
51 //==============================================================================
52 //
53 // Implementation of LennardJones612Implementation public member functions
54 //
55 //==============================================================================
56 
57 //******************************************************************************
60  KIM::ModelDriverCreate * const modelDriverCreate,
61  KIM::LengthUnit const requestedLengthUnit,
62  KIM::EnergyUnit const requestedEnergyUnit,
63  KIM::ChargeUnit const requestedChargeUnit,
64  KIM::TemperatureUnit const requestedTemperatureUnit,
65  KIM::TimeUnit const requestedTimeUnit,
66  int * const ier)
67 : numberModelSpecies_(0),
68  numberUniqueSpeciesPairs_(0),
69  shift_(0),
70  cutoffs_(0),
71  epsilons_(0),
72  sigmas_(0),
73  influenceDistance_(0.0),
74  cutoffsSq2D_(0),
75  fourEpsilonSigma6_2D_(0),
76  fourEpsilonSigma12_2D_(0),
77  twentyFourEpsilonSigma6_2D_(0),
78  fortyEightEpsilonSigma12_2D_(0),
79  oneSixtyEightEpsilonSigma6_2D_(0),
80  sixTwentyFourEpsilonSigma12_2D_(0),
81  shifts2D_(0),
82  cachedNumberOfParticles_(0)
83 {
84  FILE* parameterFilePointers[MAX_PARAMETER_FILES];
85  int numberParameterFiles;
86  modelDriverCreate->GetNumberOfParameterFiles(
87  &numberParameterFiles);
88  *ier = OpenParameterFiles(modelDriverCreate, numberParameterFiles,
89  parameterFilePointers);
90  if (*ier) return;
91 
92  *ier = ProcessParameterFiles(modelDriverCreate,
93  numberParameterFiles,
94  parameterFilePointers);
95  CloseParameterFiles(numberParameterFiles, parameterFilePointers);
96  if (*ier) return;
97 
98  *ier = ConvertUnits(modelDriverCreate,
99  requestedLengthUnit,
100  requestedEnergyUnit,
101  requestedChargeUnit,
102  requestedTemperatureUnit,
103  requestedTimeUnit);
104  if (*ier) return;
105 
106  *ier = SetReinitMutableValues(modelDriverCreate);
107  if (*ier) return;
108 
109  *ier = RegisterKIMModelSettings(modelDriverCreate);
110  if (*ier) return;
111 
112  *ier = RegisterKIMParameters(modelDriverCreate);
113  if (*ier) return;
114 
115  *ier = RegisterKIMFunctions(modelDriverCreate);
116  if (*ier) return;
117 
118  // everything is good
119  *ier = false;
120  return;
121 }
122 
123 //******************************************************************************
125 { // note: it is ok to delete a null pointer and we have ensured that
126  // everything is initialized to null
127 
128  delete [] cutoffs_;
129  Deallocate2DArray(cutoffsSq2D_);
130  delete [] epsilons_;
131  delete [] sigmas_;
132  Deallocate2DArray(fourEpsilonSigma6_2D_);
133  Deallocate2DArray(fourEpsilonSigma12_2D_);
134  Deallocate2DArray(twentyFourEpsilonSigma6_2D_);
135  Deallocate2DArray(fortyEightEpsilonSigma12_2D_);
136  Deallocate2DArray(oneSixtyEightEpsilonSigma6_2D_);
137  Deallocate2DArray(sixTwentyFourEpsilonSigma12_2D_);
138  Deallocate2DArray(shifts2D_);
139 }
140 
141 //******************************************************************************
143  KIM::ModelRefresh * const modelRefresh)
144 {
145  int ier;
146 
147  ier = SetReinitMutableValues(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 {
161  int ier;
162 
163  // KIM API Model Input compute flags
164  bool isComputeProcess_dEdr = false;
165  bool isComputeProcess_d2Edr2 = false;
166  //
167  // KIM API Model Output compute flags
168  bool isComputeEnergy = false;
169  bool isComputeForces = false;
170  bool isComputeParticleEnergy = false;
171  //
172  // KIM API Model Input
173  int const* particleSpeciesCodes = 0;
174  int const* particleContributing = 0;
175  VectorOfSizeDIM const* coordinates = 0;
176  //
177  // KIM API Model Output
178  double* energy = 0;
179  double* particleEnergy = 0;
180  VectorOfSizeDIM* forces = 0;
181  ier = SetComputeMutableValues(modelCompute, isComputeProcess_dEdr,
182  isComputeProcess_d2Edr2, isComputeEnergy,
183  isComputeForces, isComputeParticleEnergy,
185  coordinates, energy, particleEnergy, forces);
186  if (ier) return ier;
187 
188  // Skip this check for efficiency
189  //
190  // ier = CheckParticleSpecies(pkim, particleSpeciesCodes);
191  // if (ier) return ier;
192 
193  bool const isShift = (1 == shift_);
194 
195 #include "LennardJones612ImplementationComputeDispatch.cpp"
196  return ier;
197 }
198 
199 //==============================================================================
200 //
201 // Implementation of LennardJones612Implementation private member functions
202 //
203 //==============================================================================
204 
205 //******************************************************************************
207 int LennardJones612Implementation::OpenParameterFiles(
208  KIM::ModelDriverCreate * const modelDriverCreate,
209  int const numberParameterFiles,
210  FILE* parameterFilePointers[MAX_PARAMETER_FILES])
211 {
212  int ier;
213 
214  if (numberParameterFiles > MAX_PARAMETER_FILES)
215  {
216  ier = true;
217  LOG_ERROR("LennardJones612 given too many parameter files");
218  return ier;
219  }
220 
221  for (int i = 0; i < numberParameterFiles; ++i)
222  {
223  std::string paramFileName;
224  ier = modelDriverCreate->GetParameterFileName(
225  i,
226  &paramFileName);
227  if (ier)
228  {
229  LOG_ERROR("Unable to get parameter file name");
230  return ier;
231  }
232  parameterFilePointers[i] = fopen(paramFileName.c_str(), "r");
233  if (parameterFilePointers[i] == 0)
234  {
235  char message[MAXLINE];
236  sprintf(message,
237  "LennardJones612 parameter file number %d cannot be opened",
238  i);
239  ier = true;
240  LOG_ERROR(message);
241  for (int j = i - 1; i <= 0; --i)
242  {
243  fclose(parameterFilePointers[j]);
244  }
245  return ier;
246  }
247  }
248 
249  // everything is good
250  ier = false;
251  return ier;
252 }
253 
254 //******************************************************************************
255 int LennardJones612Implementation::ProcessParameterFiles(
256  KIM::ModelDriverCreate * const modelDriverCreate,
257  int const numberParameterFiles,
258  FILE* const parameterFilePointers[MAX_PARAMETER_FILES])
259 {
260  int N, ier;
261  int endOfFileFlag = 0;
262  char spec1[MAXLINE], spec2[MAXLINE], nextLine[MAXLINE];
263  char *nextLinePtr;
264  int iIndex, jIndex , indx, iiIndex, jjIndex;
265  double nextCutoff, nextEpsilon, nextSigma;
266  nextLinePtr = nextLine;
267 
268  getNextDataLine(parameterFilePointers[0], nextLinePtr,
269  MAXLINE, &endOfFileFlag);
270  ier = sscanf(nextLine, "%d %d", &N, &shift_);
271  if (ier != 2)
272  {
273  sprintf(nextLine, "unable to read first line of the parameter file");
274  ier = true;
275  LOG_ERROR(nextLine);
276  fclose(parameterFilePointers[0]);
277  return ier;
278  }
279  numberModelSpecies_ = N;
280  numberUniqueSpeciesPairs_ = ((numberModelSpecies_+1)*numberModelSpecies_)/2;
281  AllocateFreeParameterMemory();
282 
283  // set all values in the arrays to -1 for mixing later
284  for (int i = 0; i < ((N+1)*N/2); i++)
285  {
286  cutoffs_[i] = -1;
287  epsilons_[i] = -1;
288  sigmas_[i] = -1;
289  }
290 
291 
292  // keep track of known species
293  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
294  modelSpeciesMap;
295  std::vector<KIM::SpeciesName> speciesNameVector;
296  int index = 0;
297 
298  // Read and process data lines
299  getNextDataLine(parameterFilePointers[0], nextLinePtr,
300  MAXLINE, &endOfFileFlag);
301  while (endOfFileFlag == 0)
302  {
303  ier = sscanf(nextLine, "%s %s %lg %lg %lg",
304  spec1, spec2, &nextCutoff, &nextEpsilon, &nextSigma);
305  if (ier != 5)
306  {
307  sprintf(nextLine, "error reading lines of the parameter file");
308  LOG_ERROR(nextLine);
309  return true;
310  }
311 
312  // convert species strings to proper type instances
313  KIM::SpeciesName const specName1(spec1);
314  KIM::SpeciesName const specName2(spec2);
315 
316  // check for new species
317  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
318  const_iterator iIter = modelSpeciesMap.find(specName1);
319  if (iIter == modelSpeciesMap.end())
320  {
321  modelSpeciesMap[specName1] = index;
322  modelSpeciesCodeList_.push_back(index);
323  speciesNameVector.push_back(specName1);
324 
325  ier = modelDriverCreate->SetSpeciesCode(specName1, index);
326  if (ier) return ier;
327  iIndex = index;
328  index++;
329  }
330  else
331  {
332  iIndex = modelSpeciesMap[specName1];
333  }
334  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
335  const_iterator jIter = modelSpeciesMap.find(specName2);
336  if (jIter == modelSpeciesMap.end())
337  {
338  modelSpeciesMap[specName2] = index;
339  modelSpeciesCodeList_.push_back(index);
340  speciesNameVector.push_back(specName2);
341 
342  ier = modelDriverCreate->SetSpeciesCode(specName2, index);
343  if (ier) return ier;
344  jIndex = index;
345  index++;
346  }
347  else
348  {
349  jIndex = modelSpeciesMap[specName2];
350  }
351 
352  if (iIndex >= jIndex)
353  {
354  indx = jIndex*N + iIndex - (jIndex*jIndex + jIndex)/2;
355  }
356  else
357  {
358  indx = iIndex*N + jIndex - (iIndex*iIndex + iIndex)/2;
359  }
360  cutoffs_[indx] = nextCutoff;
361  epsilons_[indx] = nextEpsilon;
362  sigmas_[indx] = nextSigma;
363 
364  getNextDataLine(parameterFilePointers[0], nextLinePtr,
365  MAXLINE, &endOfFileFlag);
366  }
367 
368  // check that we got all like - like pairs
369  sprintf(nextLine, "There are not values for like-like pairs of:");
370  for (int i = 0; i < N; i++)
371  {
372  if (cutoffs_[(i*N + i - (i*i + i)/2)] == -1)
373  {
374  strcat(nextLine, " ");
375  strcat(nextLine, (speciesNameVector[i].String()).c_str());
376  ier = -1;
377  }
378  }
379  if (ier == -1)
380  {
381  LOG_ERROR(nextLine);
382  return true;
383  }
384 
385  // Perform Mixing if nessisary
386  for (int jIndex = 0; jIndex < N; jIndex++)
387  {
388  jjIndex = (jIndex*N + jIndex - (jIndex*jIndex + jIndex)/2);
389  for (int iIndex = (jIndex+1) ; iIndex < N; iIndex++)
390  {
391  indx = jIndex*N + iIndex - (jIndex*jIndex + jIndex)/2;
392  if (cutoffs_[indx] == -1)
393  {
394  iiIndex = (iIndex*N + iIndex - (iIndex*iIndex + iIndex)/2);
395  epsilons_[indx] = sqrt(epsilons_[iiIndex]*epsilons_[jjIndex]);
396  sigmas_[indx] = (sigmas_[iiIndex] + sigmas_[jjIndex])/2.0;
397  cutoffs_[indx] = (cutoffs_[iiIndex] + cutoffs_[jjIndex])/2.0;
398  }
399  }
400  }
401 
402  // everything is good
403  ier = false;
404  return ier;
405 }
406 
407 //******************************************************************************
408 void LennardJones612Implementation::getNextDataLine(
409  FILE* const filePtr, char* nextLinePtr, int const maxSize,
410  int *endOfFileFlag)
411 {
412  do
413  {
414  if(fgets(nextLinePtr, maxSize, filePtr) == NULL)
415  {
416  *endOfFileFlag = 1;
417  break;
418  }
419  while ((nextLinePtr[0] == ' ' || nextLinePtr[0] == '\t') ||
420  (nextLinePtr[0] == '\n' || nextLinePtr[0] == '\r' ))
421  {
422  nextLinePtr = (nextLinePtr + 1);
423  }
424  }
425  while ((strncmp("#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
426 }
427 
428 //******************************************************************************
429 void LennardJones612Implementation::CloseParameterFiles(
430  int const numberParameterFiles,
431  FILE* const parameterFilePointers[MAX_PARAMETER_FILES])
432 {
433  for (int i = 0; i < numberParameterFiles; ++i)
434  fclose(parameterFilePointers[i]);
435 }
436 
437 //******************************************************************************
438 void LennardJones612Implementation::AllocateFreeParameterMemory()
439 { // allocate memory for data
440  cutoffs_ = new double[numberUniqueSpeciesPairs_];
441  AllocateAndInitialize2DArray(cutoffsSq2D_, numberModelSpecies_,
442  numberModelSpecies_);
443 
444  epsilons_ = new double[numberUniqueSpeciesPairs_];
445  sigmas_ = new double[numberUniqueSpeciesPairs_];
446  AllocateAndInitialize2DArray(fourEpsilonSigma6_2D_, numberModelSpecies_,
447  numberModelSpecies_);
448  AllocateAndInitialize2DArray(fourEpsilonSigma12_2D_, numberModelSpecies_,
449  numberModelSpecies_);
450  AllocateAndInitialize2DArray(twentyFourEpsilonSigma6_2D_, numberModelSpecies_,
451  numberModelSpecies_);
452  AllocateAndInitialize2DArray(fortyEightEpsilonSigma12_2D_,
453  numberModelSpecies_, numberModelSpecies_);
454  AllocateAndInitialize2DArray(oneSixtyEightEpsilonSigma6_2D_,
455  numberModelSpecies_, numberModelSpecies_);
456  AllocateAndInitialize2DArray(sixTwentyFourEpsilonSigma12_2D_,
457  numberModelSpecies_, numberModelSpecies_);
458 
459  AllocateAndInitialize2DArray(shifts2D_, numberModelSpecies_,
460  numberModelSpecies_);
461 }
462 
463 //******************************************************************************
465 int LennardJones612Implementation::ConvertUnits(
466  KIM::ModelDriverCreate * const modelDriverCreate,
467  KIM::LengthUnit const requestedLengthUnit,
468  KIM::EnergyUnit const requestedEnergyUnit,
469  KIM::ChargeUnit const requestedChargeUnit,
470  KIM::TemperatureUnit const requestedTemperatureUnit,
471  KIM::TimeUnit const requestedTimeUnit)
472 {
473  int ier;
474 
475  // define default base units
481 
482  // changing units of cutoffs and sigmas
483  double convertLength = 1.0;
484  ier = modelDriverCreate->ConvertUnit(
485  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
486  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
487  requestedTemperatureUnit, requestedTimeUnit,
488  1.0, 0.0, 0.0, 0.0, 0.0,
489  &convertLength);
490  if (ier)
491  {
492  LOG_ERROR("Unable to convert length unit");
493  return ier;
494  }
495  if (convertLength != ONE)
496  {
497  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
498  {
499  cutoffs_[i] *= convertLength; // convert to active units
500  sigmas_[i] *= convertLength; // convert to active units
501  }
502  }
503  // changing units of epsilons
504  double convertEnergy = 1.0;
505  ier = modelDriverCreate->ConvertUnit(
506  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
507  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
508  requestedTemperatureUnit, requestedTimeUnit,
509  0.0, 1.0, 0.0, 0.0, 0.0,
510  &convertEnergy);
511  if (ier)
512  {
513  LOG_ERROR("Unable to convert energy unit");
514  return ier;
515  }
516  if (convertEnergy != ONE)
517  {
518  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
519  {
520  epsilons_[i] *= convertEnergy; // convert to active units
521  }
522  }
523 
524  // register units
525  ier = modelDriverCreate->SetUnits(
526  requestedLengthUnit,
527  requestedEnergyUnit,
528  requestedChargeUnit,
529  requestedTemperatureUnit,
530  requestedTimeUnit);
531  if (ier)
532  {
533  LOG_ERROR("Unable to set units to requested values");
534  return ier;
535  }
536 
537  // everything is good
538  ier = false;
539  return ier;
540 }
541 
542 //******************************************************************************
543 int LennardJones612Implementation::RegisterKIMModelSettings(
544  KIM::ModelDriverCreate * const modelDriverCreate)
545 {
546  // register numbering
547  int error = modelDriverCreate->SetModelNumbering(
549 
550  // register arguments
551  LOG_INFORMATION("Register argument supportStatus");
552  error = error
553  || modelDriverCreate->SetArgumentSupportStatus(
555  || modelDriverCreate->SetArgumentSupportStatus(
557  || modelDriverCreate->SetArgumentSupportStatus(
560 
561  // register callbacks
562  LOG_INFORMATION("Register callback supportStatus");
563  error = error
564  || modelDriverCreate->SetCallbackSupportStatus(
566  || modelDriverCreate->SetCallbackSupportStatus(
568 
569  return error;
570 }
571 
572 //******************************************************************************
573 int LennardJones612Implementation::RegisterKIMParameters(
574  KIM::ModelDriverCreate * const modelDriverCreate)
575 {
576  int ier = false;
577 
578  // publish parameters (order is important)
579  ier = modelDriverCreate->SetParameterPointer(1, &shift_, "shift");
580  if (ier)
581  {
582  LOG_ERROR("set_parameter shift");
583  return ier;
584  }
585  ier = modelDriverCreate->SetParameterPointer(
586  numberUniqueSpeciesPairs_, cutoffs_, "cutoffs");
587  if (ier)
588  {
589  LOG_ERROR("set_parameter cutoffs");
590  return ier;
591  }
592  ier = modelDriverCreate->SetParameterPointer(
593  numberUniqueSpeciesPairs_, epsilons_, "epsilons");
594  if (ier)
595  {
596  LOG_ERROR("set_parameter epsilons");
597  return ier;
598  }
599  ier = modelDriverCreate->SetParameterPointer(
600  numberUniqueSpeciesPairs_, sigmas_, "sigmas");
601  if (ier)
602  {
603  LOG_ERROR("set_parameter sigmas");
604  return ier;
605  }
606 
607  // everything is good
608  ier = false;
609  return ier;
610 }
611 
612 //******************************************************************************
613 int LennardJones612Implementation::RegisterKIMFunctions(
614  KIM::ModelDriverCreate * const modelDriverCreate)
615  const
616 {
617  int error;
618 
619  // register the destroy() and reinit() functions
620  error = modelDriverCreate->SetDestroyPointer(
622  || modelDriverCreate->SetRefreshPointer(
624  || modelDriverCreate->SetComputePointer(
626 
627  return error;
628 }
629 
630 //******************************************************************************
631 template<class ModelObj>
632 int LennardJones612Implementation::SetReinitMutableValues(
633  ModelObj * const modelObj)
634 { // use (possibly) new values of free parameters to compute other quantities
635  int ier;
636 
637  // update cutoffsSq, epsilons, and sigmas
638  for (int i = 0; i < numberModelSpecies_; ++i)
639  {
640  for (int j = 0; j <= i ; ++j)
641  {
642  int const index = j*numberModelSpecies_ + i - (j*j + j)/2;
643  cutoffsSq2D_[i][j] = cutoffsSq2D_[j][i]
644  = (cutoffs_[index]*cutoffs_[index]);
645  fourEpsilonSigma6_2D_[i][j] = fourEpsilonSigma6_2D_[j][i]
646  = 4.0*epsilons_[index]*pow(sigmas_[index],6.0);
647  fourEpsilonSigma12_2D_[i][j] = fourEpsilonSigma12_2D_[j][i]
648  = 4.0*epsilons_[index]*pow(sigmas_[index],12.0);
649  twentyFourEpsilonSigma6_2D_[i][j] = twentyFourEpsilonSigma6_2D_[j][i]
650  = 6.0*fourEpsilonSigma6_2D_[i][j];
651  fortyEightEpsilonSigma12_2D_[i][j] = fortyEightEpsilonSigma12_2D_[j][i]
652  = 12.0*fourEpsilonSigma12_2D_[i][j];
653  oneSixtyEightEpsilonSigma6_2D_[i][j]
654  = oneSixtyEightEpsilonSigma6_2D_[j][i]
655  = 7.0*twentyFourEpsilonSigma6_2D_[i][j];
656  sixTwentyFourEpsilonSigma12_2D_[i][j]
657  = sixTwentyFourEpsilonSigma12_2D_[j][i]
658  = 13.0*fortyEightEpsilonSigma12_2D_[i][j];
659  }
660  }
661 
662  // update cutoff value in KIM API object
663  influenceDistance_ = 0.0;
664 
665  for (int i = 0; i < numberModelSpecies_; i++)
666  {
667  int indexI = modelSpeciesCodeList_[i];
668 
669  for (int j = 0; j < numberModelSpecies_; j++)
670  {
671  int indexJ = modelSpeciesCodeList_[j];
672 
673  if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
674  {
675  influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
676  }
677  }
678  }
679 
680  influenceDistance_ = sqrt(influenceDistance_);
681  modelObj->SetInfluenceDistancePointer(&influenceDistance_);
682  modelObj->SetNeighborListCutoffsPointer(1, &influenceDistance_);
683 
684  // update shifts
685  // compute and set shifts2D_ check if minus sign
686  double const* const* const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
687  double const* const* const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
688  if (1 == shift_)
689  {
690  double phi;
691  for (int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
692  {
693  for(int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
694  {
695  int const index = jSpecies*numberModelSpecies_ + iSpecies
696  - (jSpecies*jSpecies + jSpecies)/2;
697  double const rij2 = cutoffs_[index]*cutoffs_[index];
698  double const r2iv = 1.0/rij2;
699  double const r6iv = r2iv*r2iv*r2iv;
701  shifts2D_[iSpecies][jSpecies] = shifts2D_[jSpecies][iSpecies] = phi;
702  }
703  }
704  }
705 
706  // everything is good
707  ier = false;
708  return ier;
709 }
710 
711 //******************************************************************************
713 int LennardJones612Implementation::SetComputeMutableValues(
714  KIM::ModelCompute const * const modelCompute,
715  bool& isComputeProcess_dEdr,
716  bool& isComputeProcess_d2Edr2,
717  bool& isComputeEnergy,
718  bool& isComputeForces,
719  bool& isComputeParticleEnergy,
720  int const*& particleSpeciesCodes,
721  int const*& particleContributing,
723  double*& energy,
724  double*& particleEnergy,
725  VectorOfSizeDIM*& forces)
726 {
727  int ier = true;
728 
729  // get compute flags
730  int compProcess_dEdr;
731  int compProcess_d2Edr2;
732 
734  &compProcess_dEdr);
736  &compProcess_d2Edr2);
737 
738  isComputeProcess_dEdr = compProcess_dEdr;
739  isComputeProcess_d2Edr2 = compProcess_d2Edr2;
740 
741  // double const* cutoff; // currently unused
742  // int const* numberOfSpeciesCodes; // currently unused
743  int const* numberOfParticles;
744  ier =
745  modelCompute->GetArgumentPointer(
748  || modelCompute->GetArgumentPointer(
751  || modelCompute->GetArgumentPointer(
754  || modelCompute->GetArgumentPointer(
756  (double const ** const) &coordinates)
757  || modelCompute->GetArgumentPointer(
759  &energy)
760  || modelCompute->GetArgumentPointer(
762  &particleEnergy)
763  || modelCompute->GetArgumentPointer(
765  (double const ** const) &forces);
766  if (ier)
767  {
768  LOG_ERROR("GetArgumentPointer");
769  return ier;
770  }
771 
772  isComputeEnergy = (energy != 0);
773  isComputeParticleEnergy = (particleEnergy != 0);
774  isComputeForces = (forces != 0);
775 
776  // update values
777  cachedNumberOfParticles_ = *numberOfParticles;
778 
779  // everything is good
780  ier = false;
781  return ier;
782 }
783 
784 //******************************************************************************
785 int LennardJones612Implementation::CheckParticleSpeciesCodes(
786  KIM::ModelCompute const * const modelCompute,
787  int const* const particleSpeciesCodes)
788  const
789 {
790  int ier;
791  for (int i = 0; i < cachedNumberOfParticles_; ++i)
792  {
793  if ((particleSpeciesCodes[i] < 0) ||
794  (particleSpeciesCodes[i] >= numberModelSpecies_))
795  {
796  ier = true;
797  LOG_ERROR("unsupported particle species codes detected");
798  return ier;
799  }
800  }
801 
802  // everything is good
803  ier = false;
804  return ier;
805 }
806 
807 //******************************************************************************
808 int LennardJones612Implementation::GetComputeIndex(
809  const bool& isComputeProcess_dEdr,
810  const bool& isComputeProcess_d2Edr2,
811  const bool& isComputeEnergy,
812  const bool& isComputeForces,
813  const bool& isComputeParticleEnergy,
814  const bool& isShift) const
815 {
816  //const int processdE = 2;
817  const int processd2E = 2;
818  const int energy = 2;
819  const int force = 2;
820  const int particleEnergy = 2;
821  const int shift = 2;
822 
823 
824  int index = 0;
825 
826  // processdE
827  index += (int(isComputeProcess_dEdr))
828  * processd2E * energy * force * particleEnergy * shift;
829 
830  // processd2E
831  index += (int(isComputeProcess_d2Edr2))
832  * energy * force * particleEnergy * shift;
833 
834  // energy
835  index += (int(isComputeEnergy))
836  * force * particleEnergy * shift;
837 
838  // force
839  index += (int(isComputeForces))
840  * particleEnergy * shift;
841 
842  // particleEnergy
843  index += (int(isComputeParticleEnergy))
844  * shift;
845 
846  // shift
847  index += (int(isShift));
848 
849  return index;
850 }
851 
852 //==============================================================================
853 //
854 // Implementation of helper functions
855 //
856 //==============================================================================
857 
858 //******************************************************************************
859 void AllocateAndInitialize2DArray(double**& arrayPtr, int const extentZero,
860  int const extentOne)
861 { // allocate memory and set pointers
862  arrayPtr = new double*[extentZero];
863  arrayPtr[0] = new double[extentZero * extentOne];
864  for (int i = 1; i < extentZero; ++i)
865  {
866  arrayPtr[i] = arrayPtr[i-1] + extentOne;
867  }
868 
869  // initialize
870  for (int i = 0; i < extentZero; ++i)
871  {
872  for (int j = 0; j < extentOne; ++j)
873  {
874  arrayPtr[i][j] = 0.0;
875  }
876  }
877 }
878 
879 //******************************************************************************
880 void Deallocate2DArray(double**& arrayPtr)
881 { // deallocate memory
882  if (arrayPtr != 0) delete [] arrayPtr[0];
883  delete [] arrayPtr;
884 
885  // nullify pointer
886  arrayPtr = 0;
887 }
CallbackName const ProcessD2EDr2Term
int SetParameterPointer(int const extent, int *const ptr, std::string const &description)
static int Refresh(KIM::ModelRefresh *const modelRefresh)
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) const
TimeUnit const ps
#define MAX_PARAMETER_FILES
ArgumentName const particleContributing
ArgumentName const numberOfParticles
int SetDestroyPointer(LanguageName const languageName, func *const fptr)
ArgumentName const partialForces
int SetRefreshPointer(LanguageName const languageName, func *const fptr)
ChargeUnit const e
void GetNumberOfParameterFiles(int *const numberOfParameterFiles) const
ArgumentName const partialEnergy
int GetArgumentPointer(ArgumentName const argumentName, int const **const ptr) const
double VectorOfSizeDIM[DIMENSION]
int IsCallbackPresent(CallbackName const callbackName, int *const present) const
int SetArgumentSupportStatus(ArgumentName const argumentName, SupportStatus const supportStatus)
ArgumentName const coordinates
int SetModelNumbering(Numbering const numbering)
int SetComputePointer(LanguageName const languageName, func *const fptr)
int GetParameterFileName(int const index, std::string *const parameterFileName) const
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)
static int Destroy(KIM::ModelDestroy *const modelDestroy)
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
ArgumentName const particleSpeciesCodes
LanguageName const cpp
int Refresh(KIM::ModelRefresh *const modelRefresh)
LengthUnit const A
static int Compute(KIM::ModelCompute const *const modelCompute)
int Compute(KIM::ModelCompute const *const modelCompute)
#define LOG_INFORMATION(message)
EnergyUnit const eV
#define LENNARD_JONES_PHI(exshift)
ArgumentName const partialParticleEnergy
SpeciesName const N
int SetCallbackSupportStatus(CallbackName const callbackName, SupportStatus const supportStatus)
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
void Deallocate2DArray(double **&arrayPtr)
Numbering const zeroBased
int SetSpeciesCode(SpeciesName const speciesName, int const code)
#define LOG_ERROR(message)
LogVerbosity const error
void() func()
Definition: KIM_func.hpp:40
CallbackName const ProcessDEDrTerm
SupportStatus const optional
TemperatureUnit const K