48 #define IGNORE_RESULT(fn) if(fn){} 67 : numberModelSpecies_(0),
68 numberUniqueSpeciesPairs_(0),
73 influenceDistance_(0.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),
82 cachedNumberOfParticles_(0)
85 int numberParameterFiles;
87 &numberParameterFiles);
88 *ier = OpenParameterFiles(modelDriverCreate, numberParameterFiles,
89 parameterFilePointers);
92 *ier = ProcessParameterFiles(modelDriverCreate,
94 parameterFilePointers);
95 CloseParameterFiles(numberParameterFiles, parameterFilePointers);
98 *ier = ConvertUnits(modelDriverCreate,
102 requestedTemperatureUnit,
106 *ier = SetReinitMutableValues(modelDriverCreate);
109 *ier = RegisterKIMModelSettings(modelDriverCreate);
112 *ier = RegisterKIMParameters(modelDriverCreate);
115 *ier = RegisterKIMFunctions(modelDriverCreate);
147 ier = SetReinitMutableValues(modelRefresh);
164 bool isComputeProcess_dEdr =
false;
165 bool isComputeProcess_d2Edr2 =
false;
168 bool isComputeEnergy =
false;
169 bool isComputeForces =
false;
170 bool isComputeParticleEnergy =
false;
179 double* particleEnergy = 0;
181 ier = SetComputeMutableValues(modelCompute, isComputeProcess_dEdr,
182 isComputeProcess_d2Edr2, isComputeEnergy,
183 isComputeForces, isComputeParticleEnergy,
193 bool const isShift = (1 == shift_);
195 #include "LennardJones612ImplementationComputeDispatch.cpp" 207 int LennardJones612Implementation::OpenParameterFiles(
209 int const numberParameterFiles,
217 LOG_ERROR(
"LennardJones612 given too many parameter files");
221 for (
int i = 0; i < numberParameterFiles; ++i)
223 std::string paramFileName;
229 LOG_ERROR(
"Unable to get parameter file name");
232 parameterFilePointers[i] = fopen(paramFileName.c_str(),
"r");
233 if (parameterFilePointers[i] == 0)
237 "LennardJones612 parameter file number %d cannot be opened",
241 for (
int j = i - 1; i <= 0; --i)
243 fclose(parameterFilePointers[j]);
255 int LennardJones612Implementation::ProcessParameterFiles(
257 int const numberParameterFiles,
261 int endOfFileFlag = 0;
264 int iIndex, jIndex , indx, iiIndex, jjIndex;
265 double nextCutoff, nextEpsilon, nextSigma;
266 nextLinePtr = nextLine;
268 getNextDataLine(parameterFilePointers[0], nextLinePtr,
270 ier = sscanf(nextLine,
"%d %d", &
N, &shift_);
273 sprintf(nextLine,
"unable to read first line of the parameter file");
276 fclose(parameterFilePointers[0]);
279 numberModelSpecies_ =
N;
280 numberUniqueSpeciesPairs_ = ((numberModelSpecies_+1)*numberModelSpecies_)/2;
281 AllocateFreeParameterMemory();
284 for (
int i = 0; i < ((
N+1)*
N/2); i++)
293 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
295 std::vector<KIM::SpeciesName> speciesNameVector;
299 getNextDataLine(parameterFilePointers[0], nextLinePtr,
301 while (endOfFileFlag == 0)
303 ier = sscanf(nextLine,
"%s %s %lg %lg %lg",
304 spec1, spec2, &nextCutoff, &nextEpsilon, &nextSigma);
307 sprintf(nextLine,
"error reading lines of the parameter file");
317 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
318 const_iterator iIter = modelSpeciesMap.find(specName1);
319 if (iIter == modelSpeciesMap.end())
321 modelSpeciesMap[specName1] = index;
322 modelSpeciesCodeList_.push_back(index);
323 speciesNameVector.push_back(specName1);
332 iIndex = modelSpeciesMap[specName1];
334 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
335 const_iterator jIter = modelSpeciesMap.find(specName2);
336 if (jIter == modelSpeciesMap.end())
338 modelSpeciesMap[specName2] = index;
339 modelSpeciesCodeList_.push_back(index);
340 speciesNameVector.push_back(specName2);
349 jIndex = modelSpeciesMap[specName2];
352 if (iIndex >= jIndex)
354 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
358 indx = iIndex*
N + jIndex - (iIndex*iIndex + iIndex)/2;
360 cutoffs_[indx] = nextCutoff;
361 epsilons_[indx] = nextEpsilon;
362 sigmas_[indx] = nextSigma;
364 getNextDataLine(parameterFilePointers[0], nextLinePtr,
369 sprintf(nextLine,
"There are not values for like-like pairs of:");
370 for (
int i = 0; i <
N; i++)
372 if (cutoffs_[(i*
N + i - (i*i + i)/2)] == -1)
374 strcat(nextLine,
" ");
375 strcat(nextLine, (speciesNameVector[i].String()).c_str());
386 for (
int jIndex = 0; jIndex <
N; jIndex++)
388 jjIndex = (jIndex*
N + jIndex - (jIndex*jIndex + jIndex)/2);
389 for (
int iIndex = (jIndex+1) ; iIndex <
N; iIndex++)
391 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
392 if (cutoffs_[indx] == -1)
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;
408 void LennardJones612Implementation::getNextDataLine(
409 FILE*
const filePtr,
char* nextLinePtr,
int const maxSize,
414 if(fgets(nextLinePtr, maxSize, filePtr) == NULL)
419 while ((nextLinePtr[0] ==
' ' || nextLinePtr[0] ==
'\t') ||
420 (nextLinePtr[0] ==
'\n' || nextLinePtr[0] ==
'\r' ))
422 nextLinePtr = (nextLinePtr + 1);
425 while ((strncmp(
"#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
429 void LennardJones612Implementation::CloseParameterFiles(
430 int const numberParameterFiles,
433 for (
int i = 0; i < numberParameterFiles; ++i)
434 fclose(parameterFilePointers[i]);
438 void LennardJones612Implementation::AllocateFreeParameterMemory()
440 cutoffs_ =
new double[numberUniqueSpeciesPairs_];
442 numberModelSpecies_);
444 epsilons_ =
new double[numberUniqueSpeciesPairs_];
445 sigmas_ =
new double[numberUniqueSpeciesPairs_];
447 numberModelSpecies_);
449 numberModelSpecies_);
451 numberModelSpecies_);
453 numberModelSpecies_, numberModelSpecies_);
455 numberModelSpecies_, numberModelSpecies_);
457 numberModelSpecies_, numberModelSpecies_);
460 numberModelSpecies_);
465 int LennardJones612Implementation::ConvertUnits(
483 double convertLength = 1.0;
485 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
486 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
487 requestedTemperatureUnit, requestedTimeUnit,
488 1.0, 0.0, 0.0, 0.0, 0.0,
492 LOG_ERROR(
"Unable to convert length unit");
495 if (convertLength !=
ONE)
497 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
499 cutoffs_[i] *= convertLength;
500 sigmas_[i] *= convertLength;
504 double convertEnergy = 1.0;
506 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
507 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
508 requestedTemperatureUnit, requestedTimeUnit,
509 0.0, 1.0, 0.0, 0.0, 0.0,
513 LOG_ERROR(
"Unable to convert energy unit");
516 if (convertEnergy !=
ONE)
518 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
520 epsilons_[i] *= convertEnergy;
529 requestedTemperatureUnit,
533 LOG_ERROR(
"Unable to set units to requested values");
543 int LennardJones612Implementation::RegisterKIMModelSettings(
573 int LennardJones612Implementation::RegisterKIMParameters(
586 numberUniqueSpeciesPairs_, cutoffs_,
"cutoffs");
593 numberUniqueSpeciesPairs_, epsilons_,
"epsilons");
600 numberUniqueSpeciesPairs_, sigmas_,
"sigmas");
613 int LennardJones612Implementation::RegisterKIMFunctions(
631 template<
class ModelObj>
632 int LennardJones612Implementation::SetReinitMutableValues(
633 ModelObj *
const modelObj)
638 for (
int i = 0; i < numberModelSpecies_; ++i)
640 for (
int j = 0; j <= i ; ++j)
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];
663 influenceDistance_ = 0.0;
665 for (
int i = 0; i < numberModelSpecies_; i++)
667 int indexI = modelSpeciesCodeList_[i];
669 for (
int j = 0; j < numberModelSpecies_; j++)
671 int indexJ = modelSpeciesCodeList_[j];
673 if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
675 influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
680 influenceDistance_ = sqrt(influenceDistance_);
681 modelObj->SetInfluenceDistancePointer(&influenceDistance_);
682 modelObj->SetNeighborListCutoffsPointer(1, &influenceDistance_);
686 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
687 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
691 for (
int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
693 for(
int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
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;
713 int LennardJones612Implementation::SetComputeMutableValues(
715 bool& isComputeProcess_dEdr,
716 bool& isComputeProcess_d2Edr2,
717 bool& isComputeEnergy,
718 bool& isComputeForces,
719 bool& isComputeParticleEnergy,
724 double*& particleEnergy,
730 int compProcess_dEdr;
731 int compProcess_d2Edr2;
736 &compProcess_d2Edr2);
738 isComputeProcess_dEdr = compProcess_dEdr;
739 isComputeProcess_d2Edr2 = compProcess_d2Edr2;
765 (
double const **
const) &forces);
772 isComputeEnergy = (energy != 0);
773 isComputeParticleEnergy = (particleEnergy != 0);
774 isComputeForces = (forces != 0);
785 int LennardJones612Implementation::CheckParticleSpeciesCodes(
791 for (
int i = 0; i < cachedNumberOfParticles_; ++i)
797 LOG_ERROR(
"unsupported particle species codes detected");
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 817 const int processd2E = 2;
818 const int energy = 2;
820 const int particleEnergy = 2;
827 index += (int(isComputeProcess_dEdr))
828 * processd2E * energy * force * particleEnergy * shift;
831 index += (int(isComputeProcess_d2Edr2))
832 * energy * force * particleEnergy * shift;
835 index += (int(isComputeEnergy))
836 * force * particleEnergy * shift;
839 index += (int(isComputeForces))
840 * particleEnergy * shift;
843 index += (int(isComputeParticleEnergy))
847 index += (int(isShift));
862 arrayPtr =
new double*[extentZero];
863 arrayPtr[0] =
new double[extentZero * extentOne];
864 for (
int i = 1; i < extentZero; ++i)
866 arrayPtr[i] = arrayPtr[i-1] + extentOne;
870 for (
int i = 0; i < extentZero; ++i)
872 for (
int j = 0; j < extentOne; ++j)
874 arrayPtr[i][j] = 0.0;
882 if (arrayPtr != 0)
delete [] arrayPtr[0];
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
#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)
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
int Refresh(KIM::ModelRefresh *const modelRefresh)
static int Compute(KIM::ModelCompute const *const modelCompute)
int Compute(KIM::ModelCompute const *const modelCompute)
#define LOG_INFORMATION(message)
~LennardJones612Implementation()
#define LENNARD_JONES_PHI(exshift)
ArgumentName const partialParticleEnergy
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)
CallbackName const ProcessDEDrTerm
SupportStatus const optional