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 
41 
42 #define MAXLINE 1024
43 #define IGNORE_RESULT(fn) if(fn){}
44 
45 
46 //==============================================================================
47 //
48 // Implementation of LennardJones612Implementation public member functions
49 //
50 //==============================================================================
51 
52 //******************************************************************************
55  KIM::ModelDriverCreate * const modelDriverCreate,
56  KIM::LengthUnit const requestedLengthUnit,
57  KIM::EnergyUnit const requestedEnergyUnit,
58  KIM::ChargeUnit const requestedChargeUnit,
59  KIM::TemperatureUnit const requestedTemperatureUnit,
60  KIM::TimeUnit const requestedTimeUnit,
61  int * const ier)
62 : numberModelSpecies_(0),
63  numberUniqueSpeciesPairs_(0),
64  shift_(0),
65  cutoffs_(NULL),
66  epsilons_(NULL),
67  sigmas_(NULL),
68  influenceDistance_(0.0),
69  cutoffsSq2D_(NULL),
70  paddingNeighborHints_(1),
71  halfListHints_(1),
72  fourEpsilonSigma6_2D_(NULL),
73  fourEpsilonSigma12_2D_(NULL),
74  twentyFourEpsilonSigma6_2D_(NULL),
75  fortyEightEpsilonSigma12_2D_(NULL),
76  oneSixtyEightEpsilonSigma6_2D_(NULL),
77  sixTwentyFourEpsilonSigma12_2D_(NULL),
78  shifts2D_(NULL),
79  cachedNumberOfParticles_(0)
80 {
81  FILE* parameterFilePointers[MAX_PARAMETER_FILES];
82  int numberParameterFiles;
83  modelDriverCreate->GetNumberOfParameterFiles(
84  &numberParameterFiles);
85  *ier = OpenParameterFiles(modelDriverCreate, numberParameterFiles,
86  parameterFilePointers);
87  if (*ier) return;
88 
89  *ier = ProcessParameterFiles(modelDriverCreate,
90  numberParameterFiles,
91  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 //******************************************************************************
141  KIM::ModelRefresh * const modelRefresh)
142 {
143  int ier;
144 
145  ier = SetRefreshMutableValues(modelRefresh);
146  if (ier) return ier;
147 
148  // nothing else to do for this case
149 
150  // everything is good
151  ier = false;
152  return ier;
153 }
154 
155 //******************************************************************************
157  KIM::ModelCompute const * const modelCompute,
158  KIM::ModelComputeArguments const * const modelComputeArguments)
159 {
160  int ier;
161 
162  // KIM API Model Input compute flags
163  bool isComputeProcess_dEdr = false;
164  bool isComputeProcess_d2Edr2 = false;
165  //
166  // KIM API Model Output compute flags
167  bool isComputeEnergy = false;
168  bool isComputeForces = false;
169  bool isComputeParticleEnergy = false;
170  bool isComputeVirial = false;
171  bool isComputeParticleVirial = false;
172  //
173  // KIM API Model Input
174  int const* particleSpeciesCodes = NULL;
175  int const* particleContributing = NULL;
176  VectorOfSizeDIM const* coordinates = NULL;
177  //
178  // KIM API Model Output
179  double* energy = NULL;
180  double* particleEnergy = NULL;
181  VectorOfSizeDIM* forces = NULL;
182  VectorOfSizeSix* virial = NULL;
183  VectorOfSizeSix* particleVirial = NULL;
184  ier = SetComputeMutableValues(modelComputeArguments,
185  isComputeProcess_dEdr,
186  isComputeProcess_d2Edr2, isComputeEnergy,
187  isComputeForces, isComputeParticleEnergy,
188  isComputeVirial, isComputeParticleVirial,
190  coordinates, energy, particleEnergy, forces,
191  virial, particleVirial);
192  if (ier) return ier;
193 
194  // Skip this check for efficiency
195  //
196  //ier = CheckParticleSpecies(modelComputeArguments, particleSpeciesCodes);
197  //if (ier) return ier;
198 
199  bool const isShift = (1 == shift_);
200 
201 #include "LennardJones612ImplementationComputeDispatch.cpp"
202  return ier;
203 }
204 
205 //******************************************************************************
207  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
208 {
209  int ier;
210 
211  ier = RegisterKIMComputeArgumentsSettings(modelComputeArgumentsCreate);
212  if (ier) return ier;
213 
214  // nothing else to do for this case
215 
216  // everything is good
217  ier = false;
218  return ier;
219 }
220 
221 //******************************************************************************
223  KIM::ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
224  const
225 {
226  int ier;
227 
228  // nothing else to do for this case
229 
230  // everything is good
231  ier = false;
232  return ier;
233 }
234 
235 
236 //==============================================================================
237 //
238 // Implementation of LennardJones612Implementation private member functions
239 //
240 //==============================================================================
241 
242 //******************************************************************************
243 void LennardJones612Implementation::AllocatePrivateParameterMemory()
244 {
245  // nothing to do for this case
246 }
247 
248 //******************************************************************************
249 void LennardJones612Implementation::AllocateParameterMemory()
250 { // allocate memory for data
251  cutoffs_ = new double[numberUniqueSpeciesPairs_];
252  AllocateAndInitialize2DArray(cutoffsSq2D_, numberModelSpecies_,
253  numberModelSpecies_);
254 
255  epsilons_ = new double[numberUniqueSpeciesPairs_];
256  sigmas_ = new double[numberUniqueSpeciesPairs_];
257  AllocateAndInitialize2DArray(fourEpsilonSigma6_2D_, numberModelSpecies_,
258  numberModelSpecies_);
259  AllocateAndInitialize2DArray(fourEpsilonSigma12_2D_, numberModelSpecies_,
260  numberModelSpecies_);
261  AllocateAndInitialize2DArray(twentyFourEpsilonSigma6_2D_, numberModelSpecies_,
262  numberModelSpecies_);
263  AllocateAndInitialize2DArray(fortyEightEpsilonSigma12_2D_,
264  numberModelSpecies_, numberModelSpecies_);
265  AllocateAndInitialize2DArray(oneSixtyEightEpsilonSigma6_2D_,
266  numberModelSpecies_, numberModelSpecies_);
267  AllocateAndInitialize2DArray(sixTwentyFourEpsilonSigma12_2D_,
268  numberModelSpecies_, numberModelSpecies_);
269 
270  AllocateAndInitialize2DArray(shifts2D_, numberModelSpecies_,
271  numberModelSpecies_);
272 }
273 
274 //******************************************************************************
276 int LennardJones612Implementation::OpenParameterFiles(
277  KIM::ModelDriverCreate * const modelDriverCreate,
278  int const numberParameterFiles,
279  FILE* parameterFilePointers[MAX_PARAMETER_FILES])
280 {
281  int ier;
282 
283  if (numberParameterFiles > MAX_PARAMETER_FILES)
284  {
285  ier = true;
286  LOG_ERROR("LennardJones612 given too many parameter files");
287  return ier;
288  }
289 
290  for (int i = 0; i < numberParameterFiles; ++i)
291  {
292  std::string const * paramFileName;
293  ier = modelDriverCreate->GetParameterFileName(
294  i,
295  &paramFileName);
296  if (ier)
297  {
298  LOG_ERROR("Unable to get parameter file name");
299  return ier;
300  }
301  parameterFilePointers[i] = fopen(paramFileName->c_str(), "r");
302  if (parameterFilePointers[i] == 0)
303  {
304  char message[MAXLINE];
305  sprintf(message,
306  "LennardJones612 parameter file number %d cannot be opened",
307  i);
308  ier = true;
309  LOG_ERROR(message);
310  for (int j = i - 1; i <= 0; --i)
311  {
312  fclose(parameterFilePointers[j]);
313  }
314  return ier;
315  }
316  }
317 
318  // everything is good
319  ier = false;
320  return ier;
321 }
322 
323 //******************************************************************************
325 int LennardJones612Implementation::ProcessParameterFiles(
326  KIM::ModelDriverCreate * const modelDriverCreate,
327  int const numberParameterFiles,
328  FILE* const parameterFilePointers[MAX_PARAMETER_FILES])
329 {
330  int N, ier;
331  int endOfFileFlag = 0;
332  char spec1[MAXLINE], spec2[MAXLINE], nextLine[MAXLINE];
333  char *nextLinePtr;
334  int iIndex, jIndex , indx, iiIndex, jjIndex;
335  double nextCutoff, nextEpsilon, nextSigma;
336  nextLinePtr = nextLine;
337 
338  getNextDataLine(parameterFilePointers[0], nextLinePtr,
339  MAXLINE, &endOfFileFlag);
340  ier = sscanf(nextLine, "%d %d", &N, &shift_);
341  if (ier != 2)
342  {
343  sprintf(nextLine, "unable to read first line of the parameter file");
344  ier = true;
345  LOG_ERROR(nextLine);
346  fclose(parameterFilePointers[0]);
347  return ier;
348  }
349  numberModelSpecies_ = N;
350  numberUniqueSpeciesPairs_ = ((numberModelSpecies_+1)*numberModelSpecies_)/2;
351  AllocateParameterMemory();
352 
353  // set all values in the arrays to -1 for mixing later
354  for (int i = 0; i < ((N+1)*N/2); i++)
355  {
356  cutoffs_[i] = -1;
357  epsilons_[i] = -1;
358  sigmas_[i] = -1;
359  }
360 
361 
362  // keep track of known species
363  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
364  modelSpeciesMap;
365  std::vector<KIM::SpeciesName> speciesNameVector;
366  int index = 0;
367 
368  // Read and process data lines
369  getNextDataLine(parameterFilePointers[0], nextLinePtr,
370  MAXLINE, &endOfFileFlag);
371  while (endOfFileFlag == 0)
372  {
373  ier = sscanf(nextLine, "%s %s %lg %lg %lg",
374  spec1, spec2, &nextCutoff, &nextEpsilon, &nextSigma);
375  if (ier != 5)
376  {
377  sprintf(nextLine, "error reading lines of the parameter file");
378  LOG_ERROR(nextLine);
379  return true;
380  }
381 
382  // convert species strings to proper type instances
383  KIM::SpeciesName const specName1(spec1);
384  KIM::SpeciesName const specName2(spec2);
385 
386  // check for new species
387  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
388  const_iterator iIter = modelSpeciesMap.find(specName1);
389  if (iIter == modelSpeciesMap.end())
390  {
391  modelSpeciesMap[specName1] = index;
392  modelSpeciesCodeList_.push_back(index);
393  speciesNameVector.push_back(specName1);
394 
395  ier = modelDriverCreate->SetSpeciesCode(specName1, index);
396  if (ier) return ier;
397  iIndex = index;
398  index++;
399  }
400  else
401  {
402  iIndex = modelSpeciesMap[specName1];
403  }
404  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
405  const_iterator jIter = modelSpeciesMap.find(specName2);
406  if (jIter == modelSpeciesMap.end())
407  {
408  modelSpeciesMap[specName2] = index;
409  modelSpeciesCodeList_.push_back(index);
410  speciesNameVector.push_back(specName2);
411 
412  ier = modelDriverCreate->SetSpeciesCode(specName2, index);
413  if (ier) return ier;
414  jIndex = index;
415  index++;
416  }
417  else
418  {
419  jIndex = modelSpeciesMap[specName2];
420  }
421 
422  if (iIndex >= jIndex)
423  {
424  indx = jIndex*N + iIndex - (jIndex*jIndex + jIndex)/2;
425  }
426  else
427  {
428  indx = iIndex*N + jIndex - (iIndex*iIndex + iIndex)/2;
429  }
430  cutoffs_[indx] = nextCutoff;
431  epsilons_[indx] = nextEpsilon;
432  sigmas_[indx] = nextSigma;
433 
434  getNextDataLine(parameterFilePointers[0], nextLinePtr,
435  MAXLINE, &endOfFileFlag);
436  }
437 
438  // check that we got all like - like pairs
439  sprintf(nextLine, "There are not values for like-like pairs of:");
440  for (int i = 0; i < N; i++)
441  {
442  if (cutoffs_[(i*N + i - (i*i + i)/2)] == -1)
443  {
444  strcat(nextLine, " ");
445  strcat(nextLine, (speciesNameVector[i].String()).c_str());
446  ier = -1;
447  }
448  }
449  if (ier == -1)
450  {
451  LOG_ERROR(nextLine);
452  return true;
453  }
454 
455  // Perform Mixing if nessisary
456  for (int jIndex = 0; jIndex < N; jIndex++)
457  {
458  jjIndex = (jIndex*N + jIndex - (jIndex*jIndex + jIndex)/2);
459  for (int iIndex = (jIndex+1) ; iIndex < N; iIndex++)
460  {
461  indx = jIndex*N + iIndex - (jIndex*jIndex + jIndex)/2;
462  if (cutoffs_[indx] == -1)
463  {
464  iiIndex = (iIndex*N + iIndex - (iIndex*iIndex + iIndex)/2);
465  epsilons_[indx] = sqrt(epsilons_[iiIndex]*epsilons_[jjIndex]);
466  sigmas_[indx] = (sigmas_[iiIndex] + sigmas_[jjIndex])/2.0;
467  cutoffs_[indx] = (cutoffs_[iiIndex] + cutoffs_[jjIndex])/2.0;
468  }
469  }
470  }
471 
472  // everything is good
473  ier = false;
474  return ier;
475 }
476 
477 //******************************************************************************
478 void LennardJones612Implementation::getNextDataLine(
479  FILE* const filePtr, char* nextLinePtr, int const maxSize,
480  int *endOfFileFlag)
481 {
482  do
483  {
484  if(fgets(nextLinePtr, maxSize, filePtr) == NULL)
485  {
486  *endOfFileFlag = 1;
487  break;
488  }
489  while ((nextLinePtr[0] == ' ' || nextLinePtr[0] == '\t') ||
490  (nextLinePtr[0] == '\n' || nextLinePtr[0] == '\r' ))
491  {
492  nextLinePtr = (nextLinePtr + 1);
493  }
494  }
495  while ((strncmp("#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
496 }
497 
498 //******************************************************************************
499 void LennardJones612Implementation::CloseParameterFiles(
500  int const numberParameterFiles,
501  FILE* const parameterFilePointers[MAX_PARAMETER_FILES])
502 {
503  for (int i = 0; i < numberParameterFiles; ++i)
504  fclose(parameterFilePointers[i]);
505 }
506 
507 //******************************************************************************
509 int LennardJones612Implementation::ConvertUnits(
510  KIM::ModelDriverCreate * const modelDriverCreate,
511  KIM::LengthUnit const requestedLengthUnit,
512  KIM::EnergyUnit const requestedEnergyUnit,
513  KIM::ChargeUnit const requestedChargeUnit,
514  KIM::TemperatureUnit const requestedTemperatureUnit,
515  KIM::TimeUnit const requestedTimeUnit)
516 {
517  int ier;
518 
519  // define default base units
525 
526  // changing units of cutoffs and sigmas
527  double convertLength = 1.0;
528  ier = modelDriverCreate->ConvertUnit(
529  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
530  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
531  requestedTemperatureUnit, requestedTimeUnit,
532  1.0, 0.0, 0.0, 0.0, 0.0,
533  &convertLength);
534  if (ier)
535  {
536  LOG_ERROR("Unable to convert length unit");
537  return ier;
538  }
539  if (convertLength != ONE)
540  {
541  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
542  {
543  cutoffs_[i] *= convertLength; // convert to active units
544  sigmas_[i] *= convertLength; // convert to active units
545  }
546  }
547  // changing units of epsilons
548  double convertEnergy = 1.0;
549  ier = modelDriverCreate->ConvertUnit(
550  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
551  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
552  requestedTemperatureUnit, requestedTimeUnit,
553  0.0, 1.0, 0.0, 0.0, 0.0,
554  &convertEnergy);
555  if (ier)
556  {
557  LOG_ERROR("Unable to convert energy unit");
558  return ier;
559  }
560  if (convertEnergy != ONE)
561  {
562  for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
563  {
564  epsilons_[i] *= convertEnergy; // convert to active units
565  }
566  }
567 
568  // register units
569  ier = modelDriverCreate->SetUnits(
570  requestedLengthUnit,
571  requestedEnergyUnit,
575  if (ier)
576  {
577  LOG_ERROR("Unable to set units to requested values");
578  return ier;
579  }
580 
581  // everything is good
582  ier = false;
583  return ier;
584 }
585 
586 //******************************************************************************
587 int LennardJones612Implementation::RegisterKIMModelSettings(
588  KIM::ModelDriverCreate * const modelDriverCreate) const
589 {
590  // register numbering
591  int error = modelDriverCreate->SetModelNumbering(
593 
594  return error;
595 }
596 
597 //******************************************************************************
599 int LennardJones612Implementation::RegisterKIMComputeArgumentsSettings(
600  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
601 {
602  // register arguments
603  LOG_INFORMATION("Register argument supportStatus");
604  int error =
605  modelComputeArgumentsCreate->SetArgumentSupportStatus(
608  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
611  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
614  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
617  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
620 
621 
622  // register callbacks
623  LOG_INFORMATION("Register callback supportStatus");
624  error = error
625  || modelComputeArgumentsCreate->SetCallbackSupportStatus(
628  || modelComputeArgumentsCreate->SetCallbackSupportStatus(
631 
632  return error;
633 }
634 
635 //******************************************************************************
637 int LennardJones612Implementation::RegisterKIMParameters(
638  KIM::ModelDriverCreate * const modelDriverCreate)
639 {
640  int ier = false;
641 
642  // publish parameters (order is important)
643  ier = modelDriverCreate->SetParameterPointer(1, &shift_, "shift");
644  if (ier)
645  {
646  LOG_ERROR("set_parameter shift");
647  return ier;
648  }
649  ier = modelDriverCreate->SetParameterPointer(
650  numberUniqueSpeciesPairs_, cutoffs_, "cutoffs");
651  if (ier)
652  {
653  LOG_ERROR("set_parameter cutoffs");
654  return ier;
655  }
656  ier = modelDriverCreate->SetParameterPointer(
657  numberUniqueSpeciesPairs_, epsilons_, "epsilons");
658  if (ier)
659  {
660  LOG_ERROR("set_parameter epsilons");
661  return ier;
662  }
663  ier = modelDriverCreate->SetParameterPointer(
664  numberUniqueSpeciesPairs_, sigmas_, "sigmas");
665  if (ier)
666  {
667  LOG_ERROR("set_parameter sigmas");
668  return ier;
669  }
670 
671  // everything is good
672  ier = false;
673  return ier;
674 }
675 
676 //******************************************************************************
677 int LennardJones612Implementation::RegisterKIMFunctions(
678  KIM::ModelDriverCreate * const modelDriverCreate)
679  const
680 {
681  int error;
682 
683  // register the destroy() and reinit() functions
684  error = modelDriverCreate->SetDestroyPointer(
686  || modelDriverCreate->SetRefreshPointer(
688  || modelDriverCreate->SetComputePointer(
690  || modelDriverCreate->SetComputeArgumentsCreatePointer(
693  || modelDriverCreate->SetComputeArgumentsDestroyPointer(
696 
697  return error;
698 }
699 
700 //******************************************************************************
701 template<class ModelObj>
702 int LennardJones612Implementation::SetRefreshMutableValues(
703  ModelObj * const modelObj)
704 { // use (possibly) new values of parameters to compute other quantities
705  // NOTE: This function is templated because it's called with both a
706  // modelDriverCreate object during initialization and with a
707  // modelRefresh object when the Model's parameters have been altered
708  int ier;
709 
710  // update cutoffsSq, epsilons, and sigmas
711  for (int i = 0; i < numberModelSpecies_; ++i)
712  {
713  for (int j = 0; j <= i ; ++j)
714  {
715  int const index = j*numberModelSpecies_ + i - (j*j + j)/2;
716  cutoffsSq2D_[i][j] = cutoffsSq2D_[j][i]
717  = (cutoffs_[index]*cutoffs_[index]);
718  fourEpsilonSigma6_2D_[i][j] = fourEpsilonSigma6_2D_[j][i]
719  = 4.0*epsilons_[index]*pow(sigmas_[index],6.0);
720  fourEpsilonSigma12_2D_[i][j] = fourEpsilonSigma12_2D_[j][i]
721  = 4.0*epsilons_[index]*pow(sigmas_[index],12.0);
722  twentyFourEpsilonSigma6_2D_[i][j] = twentyFourEpsilonSigma6_2D_[j][i]
723  = 6.0*fourEpsilonSigma6_2D_[i][j];
724  fortyEightEpsilonSigma12_2D_[i][j] = fortyEightEpsilonSigma12_2D_[j][i]
725  = 12.0*fourEpsilonSigma12_2D_[i][j];
726  oneSixtyEightEpsilonSigma6_2D_[i][j]
727  = oneSixtyEightEpsilonSigma6_2D_[j][i]
728  = 7.0*twentyFourEpsilonSigma6_2D_[i][j];
729  sixTwentyFourEpsilonSigma12_2D_[i][j]
730  = sixTwentyFourEpsilonSigma12_2D_[j][i]
731  = 13.0*fortyEightEpsilonSigma12_2D_[i][j];
732  }
733  }
734 
735  // update cutoff value in KIM API object
736  influenceDistance_ = 0.0;
737 
738  for (int i = 0; i < numberModelSpecies_; i++)
739  {
740  int indexI = modelSpeciesCodeList_[i];
741 
742  for (int j = 0; j < numberModelSpecies_; j++)
743  {
744  int indexJ = modelSpeciesCodeList_[j];
745 
746  if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
747  {
748  influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
749  }
750  }
751  }
752 
753  influenceDistance_ = sqrt(influenceDistance_);
754  modelObj->SetInfluenceDistancePointer(&influenceDistance_);
755  modelObj->SetNeighborListPointers(1,
756  &influenceDistance_,
757  &paddingNeighborHints_,
758  &halfListHints_);
759 
760  // update shifts
761  // compute and set shifts2D_ check if minus sign
762  double const* const* const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
763  double const* const* const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
764  if (1 == shift_)
765  {
766  double phi;
767  for (int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
768  {
769  for(int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
770  {
771  int const index = jSpecies*numberModelSpecies_ + iSpecies
772  - (jSpecies*jSpecies + jSpecies)/2;
773  double const rij2 = cutoffs_[index]*cutoffs_[index];
774  double const r2iv = 1.0/rij2;
775  double const r6iv = r2iv*r2iv*r2iv;
777  shifts2D_[iSpecies][jSpecies] = shifts2D_[jSpecies][iSpecies] = phi;
778  }
779  }
780  }
781 
782  // everything is good
783  ier = false;
784  return ier;
785 }
786 
787 //******************************************************************************
789 int LennardJones612Implementation::SetComputeMutableValues(
790  KIM::ModelComputeArguments const * const modelComputeArguments,
791  bool& isComputeProcess_dEdr,
792  bool& isComputeProcess_d2Edr2,
793  bool& isComputeEnergy,
794  bool& isComputeForces,
795  bool& isComputeParticleEnergy,
796  bool& isComputeVirial,
797  bool& isComputeParticleVirial,
798  int const*& particleSpeciesCodes,
799  int const*& particleContributing,
801  double*& energy,
802  double*& particleEnergy,
803  VectorOfSizeDIM*& forces,
804  VectorOfSizeSix*& virial,
805  VectorOfSizeSix*& particleVirial)
806 {
807  int ier = true;
808 
809  // get compute flags
810  int compProcess_dEdr;
811  int compProcess_d2Edr2;
812 
813  modelComputeArguments->IsCallbackPresent(
815  &compProcess_dEdr);
816  modelComputeArguments->IsCallbackPresent(
818  &compProcess_d2Edr2);
819 
820  isComputeProcess_dEdr = compProcess_dEdr;
821  isComputeProcess_d2Edr2 = compProcess_d2Edr2;
822 
823  int const* numberOfParticles;
824  ier =
825  modelComputeArguments->GetArgumentPointer(
828  || modelComputeArguments->GetArgumentPointer(
831  || modelComputeArguments->GetArgumentPointer(
834  || modelComputeArguments->GetArgumentPointer(
836  (double const ** const) &coordinates)
837  || modelComputeArguments->GetArgumentPointer(
839  &energy)
840  || modelComputeArguments->GetArgumentPointer(
842  &particleEnergy)
843  || modelComputeArguments->GetArgumentPointer(
845  (double const ** const) &forces)
846  || modelComputeArguments->GetArgumentPointer(
848  (double const ** const) &virial)
849  || modelComputeArguments->GetArgumentPointer(
851  (double const ** const) &particleVirial);
852  if (ier)
853  {
854  LOG_ERROR("GetArgumentPointer");
855  return ier;
856  }
857 
858  isComputeEnergy = (energy != NULL);
859  isComputeParticleEnergy = (particleEnergy != NULL);
860  isComputeForces = (forces != NULL);
861  isComputeVirial = (virial != NULL);
862  isComputeParticleVirial = (particleVirial != NULL);
863 
864  // update values
865  cachedNumberOfParticles_ = *numberOfParticles;
866 
867  // everything is good
868  ier = false;
869  return ier;
870 }
871 
872 //******************************************************************************
874 int LennardJones612Implementation::CheckParticleSpeciesCodes(
875  KIM::ModelCompute const * const modelCompute,
876  int const* const particleSpeciesCodes)
877  const
878 {
879  int ier;
880  for (int i = 0; i < cachedNumberOfParticles_; ++i)
881  {
882  if ((particleSpeciesCodes[i] < 0) ||
883  (particleSpeciesCodes[i] >= numberModelSpecies_))
884  {
885  ier = true;
886  LOG_ERROR("unsupported particle species codes detected");
887  return ier;
888  }
889  }
890 
891  // everything is good
892  ier = false;
893  return ier;
894 }
895 
896 //******************************************************************************
897 int LennardJones612Implementation::GetComputeIndex(
898  const bool& isComputeProcess_dEdr,
899  const bool& isComputeProcess_d2Edr2,
900  const bool& isComputeEnergy,
901  const bool& isComputeForces,
902  const bool& isComputeParticleEnergy,
903  const bool& isComputeVirial,
904  const bool& isComputeParticleVirial,
905  const bool& isShift) const
906 {
907  //const int processdE = 2;
908  const int processd2E = 2;
909  const int energy = 2;
910  const int force = 2;
911  const int particleEnergy = 2;
912  const int virial = 2;
913  const int particleVirial = 2;
914  const int shift = 2;
915 
916 
917  int index = 0;
918 
919  // processdE
920  index += (int(isComputeProcess_dEdr))
921  * processd2E * energy * force * particleEnergy * virial
922  * particleVirial* shift;
923 
924  // processd2E
925  index += (int(isComputeProcess_d2Edr2))
926  * energy * force * particleEnergy * virial * particleVirial * shift;
927 
928  // energy
929  index += (int(isComputeEnergy))
930  * force * particleEnergy * virial * particleVirial * shift;
931 
932  // force
933  index += (int(isComputeForces))
934  * particleEnergy * virial * particleVirial * shift;
935 
936  // particleEnergy
937  index += (int(isComputeParticleEnergy))
938  * virial * particleVirial * shift;
939 
940  // virial
941  index += (int(isComputeVirial))
942  * particleVirial * shift;
943 
944  // particleVirial
945  index += (int(isComputeParticleVirial))
946  * shift;
947 
948  // shift
949  index += (int(isShift));
950 
951  return index;
952 }
953 
954 //******************************************************************************
955 void LennardJones612Implementation::ProcessVirialTerm(
956  const double& dEidr,
957  const double& rij,
958  const double* const r_ij,
959  const int& i,
960  const int& j,
961  VectorOfSizeSix virial) const
962 {
963  double const v = dEidr/rij;
964 
965  virial[0] += v * r_ij[0] * r_ij[0];
966  virial[1] += v * r_ij[1] * r_ij[1];
967  virial[2] += v * r_ij[2] * r_ij[2];
968  virial[3] += v * r_ij[1] * r_ij[2];
969  virial[4] += v * r_ij[0] * r_ij[2];
970  virial[5] += v * r_ij[0] * r_ij[1];
971 }
972 
973 //******************************************************************************
974 void LennardJones612Implementation::ProcessParticleVirialTerm(
975  const double& dEidr,
976  const double& rij,
977  const double* const r_ij,
978  const int& i,
979  const int& j,
980  VectorOfSizeSix* const particleVirial) const
981 {
982  double const v = dEidr/rij;
983  VectorOfSizeSix vir;
984 
985  vir[0] = 0.5 * v * r_ij[0] * r_ij[0];
986  vir[1] = 0.5 * v * r_ij[1] * r_ij[1];
987  vir[2] = 0.5 * v * r_ij[2] * r_ij[2];
988  vir[3] = 0.5 * v * r_ij[1] * r_ij[2];
989  vir[4] = 0.5 * v * r_ij[0] * r_ij[2];
990  vir[5] = 0.5 * v * r_ij[0] * r_ij[1];
991 
992  for (int k = 0; k < 6; ++k)
993  {
994  particleVirial[i][k] += vir[k];
995  particleVirial[j][k] += vir[k];
996  }
997 }
998 
999 //==============================================================================
1000 //
1001 // Implementation of helper functions
1002 //
1003 //==============================================================================
1004 
1005 //******************************************************************************
1006 void AllocateAndInitialize2DArray(double**& arrayPtr, int const extentZero,
1007  int const extentOne)
1008 { // allocate memory and set pointers
1009  arrayPtr = new double*[extentZero];
1010  arrayPtr[0] = new double[extentZero * extentOne];
1011  for (int i = 1; i < extentZero; ++i)
1012  {
1013  arrayPtr[i] = arrayPtr[i-1] + extentOne;
1014  }
1015 
1016  // initialize
1017  for (int i = 0; i < extentZero; ++i)
1018  {
1019  for (int j = 0; j < extentOne; ++j)
1020  {
1021  arrayPtr[i][j] = 0.0;
1022  }
1023  }
1024 }
1025 
1026 //******************************************************************************
1027 void Deallocate2DArray(double**& arrayPtr)
1028 { // deallocate memory
1029  if (arrayPtr != NULL) delete [] arrayPtr[0];
1030  delete [] arrayPtr;
1031 
1032  // nullify pointer
1033  arrayPtr = NULL;
1034 }
int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
TemperatureUnit const unused
ComputeArgumentName const partialParticleEnergy
int SetParameterPointer(int const extent, int *const ptr, std::string const &description)
static int Refresh(KIM::ModelRefresh *const modelRefresh)
int SetComputeArgumentsCreatePointer(LanguageName const languageName, func *const fptr)
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
ComputeArgumentName const coordinates
int SetDestroyPointer(LanguageName const languageName, func *const fptr)
int SetRefreshPointer(LanguageName const languageName, func *const fptr)
ComputeArgumentName const partialVirial
ChargeUnit const unused
ChargeUnit const e
void GetNumberOfParameterFiles(int *const numberOfParameterFiles) const
ComputeCallbackName const ProcessD2EDr2Term
double VectorOfSizeDIM[DIMENSION]
int GetArgumentPointer(ComputeArgumentName const computeArgumentName, int const **const ptr) const
int SetModelNumbering(Numbering const numbering)
ComputeArgumentName const particleContributing
static int ComputeArgumentsCreate(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
ComputeArgumentName const partialEnergy
int SetComputePointer(LanguageName const languageName, func *const fptr)
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)
int GetParameterFileName(int const index, std::string const **const parameterFileName) const
static int Destroy(KIM::ModelDestroy *const modelDestroy)
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
#define LOG_ERROR(message)
LanguageName const cpp
ComputeArgumentName const partialParticleVirial
ComputeCallbackName const ProcessDEDrTerm
int Refresh(KIM::ModelRefresh *const modelRefresh)
LengthUnit const A
ComputeArgumentName const partialForces
#define LOG_INFORMATION(message)
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
EnergyUnit const eV
static int ComputeArgumentsDestroy(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
ComputeArgumentName const particleSpeciesCodes
#define LENNARD_JONES_PHI(exshift)
int SetArgumentSupportStatus(ComputeArgumentName const clomputeArgumentName, SupportStatus const supportStatus)
ComputeArgumentName const numberOfParticles
SpeciesName const N
double VectorOfSizeSix[6]
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
static int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
int SetSpeciesCode(SpeciesName const speciesName, int const code)
TimeUnit const unused
int SetCallbackSupportStatus(ComputeCallbackName const computeCallbackName, SupportStatus const supportStatus)
int SetComputeArgumentsDestroyPointer(LanguageName const languageName, func *const fptr)
LogVerbosity const error
int IsCallbackPresent(ComputeCallbackName const computeCallbackName, int *const present) const
void() func()
Definition: KIM_func.hpp:40
SupportStatus const optional
TemperatureUnit const K