43 #define IGNORE_RESULT(fn) if(fn){} 62 : numberModelSpecies_(0),
63 numberUniqueSpeciesPairs_(0),
68 influenceDistance_(0.0),
70 paddingNeighborHints_(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),
79 cachedNumberOfParticles_(0)
82 int numberParameterFiles;
84 &numberParameterFiles);
85 *ier = OpenParameterFiles(modelDriverCreate, numberParameterFiles,
86 parameterFilePointers);
89 *ier = ProcessParameterFiles(modelDriverCreate,
91 parameterFilePointers);
92 CloseParameterFiles(numberParameterFiles, parameterFilePointers);
95 *ier = ConvertUnits(modelDriverCreate,
99 requestedTemperatureUnit,
103 *ier = SetRefreshMutableValues(modelDriverCreate);
106 *ier = RegisterKIMModelSettings(modelDriverCreate);
109 *ier = RegisterKIMParameters(modelDriverCreate);
112 *ier = RegisterKIMFunctions(modelDriverCreate);
145 ier = SetRefreshMutableValues(modelRefresh);
163 bool isComputeProcess_dEdr =
false;
164 bool isComputeProcess_d2Edr2 =
false;
167 bool isComputeEnergy =
false;
168 bool isComputeForces =
false;
169 bool isComputeParticleEnergy =
false;
170 bool isComputeVirial =
false;
171 bool isComputeParticleVirial =
false;
179 double* energy = NULL;
180 double* particleEnergy = NULL;
184 ier = SetComputeMutableValues(modelComputeArguments,
185 isComputeProcess_dEdr,
186 isComputeProcess_d2Edr2, isComputeEnergy,
187 isComputeForces, isComputeParticleEnergy,
188 isComputeVirial, isComputeParticleVirial,
191 virial, particleVirial);
199 bool const isShift = (1 == shift_);
201 #include "LennardJones612ImplementationComputeDispatch.cpp" 211 ier = RegisterKIMComputeArgumentsSettings(modelComputeArgumentsCreate);
243 void LennardJones612Implementation::AllocatePrivateParameterMemory()
249 void LennardJones612Implementation::AllocateParameterMemory()
251 cutoffs_ =
new double[numberUniqueSpeciesPairs_];
253 numberModelSpecies_);
255 epsilons_ =
new double[numberUniqueSpeciesPairs_];
256 sigmas_ =
new double[numberUniqueSpeciesPairs_];
258 numberModelSpecies_);
260 numberModelSpecies_);
262 numberModelSpecies_);
264 numberModelSpecies_, numberModelSpecies_);
266 numberModelSpecies_, numberModelSpecies_);
268 numberModelSpecies_, numberModelSpecies_);
271 numberModelSpecies_);
276 int LennardJones612Implementation::OpenParameterFiles(
278 int const numberParameterFiles,
286 LOG_ERROR(
"LennardJones612 given too many parameter files");
290 for (
int i = 0; i < numberParameterFiles; ++i)
292 std::string
const * paramFileName;
298 LOG_ERROR(
"Unable to get parameter file name");
301 parameterFilePointers[i] = fopen(paramFileName->c_str(),
"r");
302 if (parameterFilePointers[i] == 0)
306 "LennardJones612 parameter file number %d cannot be opened",
310 for (
int j = i - 1; i <= 0; --i)
312 fclose(parameterFilePointers[j]);
325 int LennardJones612Implementation::ProcessParameterFiles(
327 int const numberParameterFiles,
331 int endOfFileFlag = 0;
334 int iIndex, jIndex , indx, iiIndex, jjIndex;
335 double nextCutoff, nextEpsilon, nextSigma;
336 nextLinePtr = nextLine;
338 getNextDataLine(parameterFilePointers[0], nextLinePtr,
340 ier = sscanf(nextLine,
"%d %d", &
N, &shift_);
343 sprintf(nextLine,
"unable to read first line of the parameter file");
346 fclose(parameterFilePointers[0]);
349 numberModelSpecies_ =
N;
350 numberUniqueSpeciesPairs_ = ((numberModelSpecies_+1)*numberModelSpecies_)/2;
351 AllocateParameterMemory();
354 for (
int i = 0; i < ((
N+1)*
N/2); i++)
363 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
365 std::vector<KIM::SpeciesName> speciesNameVector;
369 getNextDataLine(parameterFilePointers[0], nextLinePtr,
371 while (endOfFileFlag == 0)
373 ier = sscanf(nextLine,
"%s %s %lg %lg %lg",
374 spec1, spec2, &nextCutoff, &nextEpsilon, &nextSigma);
377 sprintf(nextLine,
"error reading lines of the parameter file");
387 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
388 const_iterator iIter = modelSpeciesMap.find(specName1);
389 if (iIter == modelSpeciesMap.end())
391 modelSpeciesMap[specName1] = index;
392 modelSpeciesCodeList_.push_back(index);
393 speciesNameVector.push_back(specName1);
402 iIndex = modelSpeciesMap[specName1];
404 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
405 const_iterator jIter = modelSpeciesMap.find(specName2);
406 if (jIter == modelSpeciesMap.end())
408 modelSpeciesMap[specName2] = index;
409 modelSpeciesCodeList_.push_back(index);
410 speciesNameVector.push_back(specName2);
419 jIndex = modelSpeciesMap[specName2];
422 if (iIndex >= jIndex)
424 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
428 indx = iIndex*
N + jIndex - (iIndex*iIndex + iIndex)/2;
430 cutoffs_[indx] = nextCutoff;
431 epsilons_[indx] = nextEpsilon;
432 sigmas_[indx] = nextSigma;
434 getNextDataLine(parameterFilePointers[0], nextLinePtr,
439 sprintf(nextLine,
"There are not values for like-like pairs of:");
440 for (
int i = 0; i <
N; i++)
442 if (cutoffs_[(i*
N + i - (i*i + i)/2)] == -1)
444 strcat(nextLine,
" ");
445 strcat(nextLine, (speciesNameVector[i].String()).c_str());
456 for (
int jIndex = 0; jIndex <
N; jIndex++)
458 jjIndex = (jIndex*
N + jIndex - (jIndex*jIndex + jIndex)/2);
459 for (
int iIndex = (jIndex+1) ; iIndex <
N; iIndex++)
461 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
462 if (cutoffs_[indx] == -1)
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;
478 void LennardJones612Implementation::getNextDataLine(
479 FILE*
const filePtr,
char* nextLinePtr,
int const maxSize,
484 if(fgets(nextLinePtr, maxSize, filePtr) == NULL)
489 while ((nextLinePtr[0] ==
' ' || nextLinePtr[0] ==
'\t') ||
490 (nextLinePtr[0] ==
'\n' || nextLinePtr[0] ==
'\r' ))
492 nextLinePtr = (nextLinePtr + 1);
495 while ((strncmp(
"#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
499 void LennardJones612Implementation::CloseParameterFiles(
500 int const numberParameterFiles,
503 for (
int i = 0; i < numberParameterFiles; ++i)
504 fclose(parameterFilePointers[i]);
509 int LennardJones612Implementation::ConvertUnits(
527 double convertLength = 1.0;
529 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
530 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
531 requestedTemperatureUnit, requestedTimeUnit,
532 1.0, 0.0, 0.0, 0.0, 0.0,
536 LOG_ERROR(
"Unable to convert length unit");
539 if (convertLength !=
ONE)
541 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
543 cutoffs_[i] *= convertLength;
544 sigmas_[i] *= convertLength;
548 double convertEnergy = 1.0;
550 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
551 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
552 requestedTemperatureUnit, requestedTimeUnit,
553 0.0, 1.0, 0.0, 0.0, 0.0,
557 LOG_ERROR(
"Unable to convert energy unit");
560 if (convertEnergy !=
ONE)
562 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
564 epsilons_[i] *= convertEnergy;
577 LOG_ERROR(
"Unable to set units to requested values");
587 int LennardJones612Implementation::RegisterKIMModelSettings(
599 int LennardJones612Implementation::RegisterKIMComputeArgumentsSettings(
637 int LennardJones612Implementation::RegisterKIMParameters(
650 numberUniqueSpeciesPairs_, cutoffs_,
"cutoffs");
657 numberUniqueSpeciesPairs_, epsilons_,
"epsilons");
664 numberUniqueSpeciesPairs_, sigmas_,
"sigmas");
677 int LennardJones612Implementation::RegisterKIMFunctions(
701 template<
class ModelObj>
702 int LennardJones612Implementation::SetRefreshMutableValues(
703 ModelObj *
const modelObj)
711 for (
int i = 0; i < numberModelSpecies_; ++i)
713 for (
int j = 0; j <= i ; ++j)
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];
736 influenceDistance_ = 0.0;
738 for (
int i = 0; i < numberModelSpecies_; i++)
740 int indexI = modelSpeciesCodeList_[i];
742 for (
int j = 0; j < numberModelSpecies_; j++)
744 int indexJ = modelSpeciesCodeList_[j];
746 if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
748 influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
753 influenceDistance_ = sqrt(influenceDistance_);
754 modelObj->SetInfluenceDistancePointer(&influenceDistance_);
755 modelObj->SetNeighborListPointers(1,
757 &paddingNeighborHints_,
762 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
763 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
767 for (
int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
769 for(
int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
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;
789 int LennardJones612Implementation::SetComputeMutableValues(
791 bool& isComputeProcess_dEdr,
792 bool& isComputeProcess_d2Edr2,
793 bool& isComputeEnergy,
794 bool& isComputeForces,
795 bool& isComputeParticleEnergy,
796 bool& isComputeVirial,
797 bool& isComputeParticleVirial,
802 double*& particleEnergy,
810 int compProcess_dEdr;
811 int compProcess_d2Edr2;
818 &compProcess_d2Edr2);
820 isComputeProcess_dEdr = compProcess_dEdr;
821 isComputeProcess_d2Edr2 = compProcess_d2Edr2;
845 (
double const **
const) &forces)
848 (
double const **
const) &virial)
851 (
double const **
const) &particleVirial);
858 isComputeEnergy = (energy != NULL);
859 isComputeParticleEnergy = (particleEnergy != NULL);
860 isComputeForces = (forces != NULL);
861 isComputeVirial = (virial != NULL);
862 isComputeParticleVirial = (particleVirial != NULL);
874 int LennardJones612Implementation::CheckParticleSpeciesCodes(
880 for (
int i = 0; i < cachedNumberOfParticles_; ++i)
886 LOG_ERROR(
"unsupported particle species codes detected");
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 908 const int processd2E = 2;
909 const int energy = 2;
911 const int particleEnergy = 2;
912 const int virial = 2;
913 const int particleVirial = 2;
920 index += (int(isComputeProcess_dEdr))
921 * processd2E * energy * force * particleEnergy * virial
922 * particleVirial* shift;
925 index += (int(isComputeProcess_d2Edr2))
926 * energy * force * particleEnergy * virial * particleVirial * shift;
929 index += (int(isComputeEnergy))
930 * force * particleEnergy * virial * particleVirial * shift;
933 index += (int(isComputeForces))
934 * particleEnergy * virial * particleVirial * shift;
937 index += (int(isComputeParticleEnergy))
938 * virial * particleVirial * shift;
941 index += (int(isComputeVirial))
942 * particleVirial * shift;
945 index += (int(isComputeParticleVirial))
949 index += (int(isShift));
955 void LennardJones612Implementation::ProcessVirialTerm(
958 const double*
const r_ij,
963 double const v = dEidr/rij;
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];
974 void LennardJones612Implementation::ProcessParticleVirialTerm(
977 const double*
const r_ij,
982 double const v = dEidr/rij;
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];
992 for (
int k = 0; k < 6; ++k)
994 particleVirial[i][k] += vir[k];
995 particleVirial[j][k] += vir[k];
1007 int const extentOne)
1009 arrayPtr =
new double*[extentZero];
1010 arrayPtr[0] =
new double[extentZero * extentOne];
1011 for (
int i = 1; i < extentZero; ++i)
1013 arrayPtr[i] = arrayPtr[i-1] + extentOne;
1017 for (
int i = 0; i < extentZero; ++i)
1019 for (
int j = 0; j < extentOne; ++j)
1021 arrayPtr[i][j] = 0.0;
1029 if (arrayPtr != NULL)
delete [] arrayPtr[0];
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
#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
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)
ComputeArgumentName const partialParticleVirial
ComputeCallbackName const ProcessDEDrTerm
int Refresh(KIM::ModelRefresh *const modelRefresh)
ComputeArgumentName const partialForces
#define LOG_INFORMATION(message)
~LennardJones612Implementation()
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
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
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)
int SetCallbackSupportStatus(ComputeCallbackName const computeCallbackName, SupportStatus const supportStatus)
int SetComputeArgumentsDestroyPointer(LanguageName const languageName, func *const fptr)
int IsCallbackPresent(ComputeCallbackName const computeCallbackName, int *const present) const
SupportStatus const optional