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