43 #define IGNORE_RESULT(fn) if(fn){} 62 : numberModelSpecies_(0),
63 numberUniqueSpeciesPairs_(0),
68 influenceDistance_(0.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),
77 cachedNumberOfParticles_(0)
80 int numberParameterFiles;
82 &numberParameterFiles);
83 *ier = OpenParameterFiles(modelDriverCreate, numberParameterFiles,
84 parameterFilePointers);
87 *ier = ProcessParameterFiles(modelDriverCreate,
89 parameterFilePointers);
90 CloseParameterFiles(numberParameterFiles, parameterFilePointers);
93 *ier = ConvertUnits(modelDriverCreate,
97 requestedTemperatureUnit,
101 *ier = SetRefreshMutableValues(modelDriverCreate);
104 *ier = RegisterKIMModelSettings(modelDriverCreate);
107 *ier = RegisterKIMParameters(modelDriverCreate);
110 *ier = RegisterKIMFunctions(modelDriverCreate);
143 ier = SetRefreshMutableValues(modelRefresh);
161 bool isComputeProcess_dEdr =
false;
162 bool isComputeProcess_d2Edr2 =
false;
165 bool isComputeEnergy =
false;
166 bool isComputeForces =
false;
167 bool isComputeParticleEnergy =
false;
168 bool isComputeVirial =
false;
169 bool isComputeParticleVirial =
false;
178 double* particleEnergy = 0;
182 ier = SetComputeMutableValues(modelComputeArguments,
183 isComputeProcess_dEdr,
184 isComputeProcess_d2Edr2, isComputeEnergy,
185 isComputeForces, isComputeParticleEnergy,
186 isComputeVirial, isComputeParticleVirial,
189 virial, particleVirial);
197 bool const isShift = (1 == shift_);
199 #include "LennardJones612ImplementationComputeDispatch.cpp" 209 ier = RegisterKIMComputeArgumentsSettings(modelComputeArgumentsCreate);
241 void LennardJones612Implementation::AllocatePrivateParameterMemory()
247 void LennardJones612Implementation::AllocateParameterMemory()
249 cutoffs_ =
new double[numberUniqueSpeciesPairs_];
251 numberModelSpecies_);
253 epsilons_ =
new double[numberUniqueSpeciesPairs_];
254 sigmas_ =
new double[numberUniqueSpeciesPairs_];
256 numberModelSpecies_);
258 numberModelSpecies_);
260 numberModelSpecies_);
262 numberModelSpecies_, numberModelSpecies_);
264 numberModelSpecies_, numberModelSpecies_);
266 numberModelSpecies_, numberModelSpecies_);
269 numberModelSpecies_);
274 int LennardJones612Implementation::OpenParameterFiles(
276 int const numberParameterFiles,
284 LOG_ERROR(
"LennardJones612 given too many parameter files");
288 for (
int i = 0; i < numberParameterFiles; ++i)
290 std::string
const * paramFileName;
296 LOG_ERROR(
"Unable to get parameter file name");
299 parameterFilePointers[i] = fopen(paramFileName->c_str(),
"r");
300 if (parameterFilePointers[i] == 0)
304 "LennardJones612 parameter file number %d cannot be opened",
308 for (
int j = i - 1; i <= 0; --i)
310 fclose(parameterFilePointers[j]);
323 int LennardJones612Implementation::ProcessParameterFiles(
325 int const numberParameterFiles,
329 int endOfFileFlag = 0;
332 int iIndex, jIndex , indx, iiIndex, jjIndex;
333 double nextCutoff, nextEpsilon, nextSigma;
334 nextLinePtr = nextLine;
336 getNextDataLine(parameterFilePointers[0], nextLinePtr,
338 ier = sscanf(nextLine,
"%d %d", &
N, &shift_);
341 sprintf(nextLine,
"unable to read first line of the parameter file");
344 fclose(parameterFilePointers[0]);
347 numberModelSpecies_ =
N;
348 numberUniqueSpeciesPairs_ = ((numberModelSpecies_+1)*numberModelSpecies_)/2;
349 AllocateParameterMemory();
352 for (
int i = 0; i < ((
N+1)*
N/2); i++)
361 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
363 std::vector<KIM::SpeciesName> speciesNameVector;
367 getNextDataLine(parameterFilePointers[0], nextLinePtr,
369 while (endOfFileFlag == 0)
371 ier = sscanf(nextLine,
"%s %s %lg %lg %lg",
372 spec1, spec2, &nextCutoff, &nextEpsilon, &nextSigma);
375 sprintf(nextLine,
"error reading lines of the parameter file");
385 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
386 const_iterator iIter = modelSpeciesMap.find(specName1);
387 if (iIter == modelSpeciesMap.end())
389 modelSpeciesMap[specName1] = index;
390 modelSpeciesCodeList_.push_back(index);
391 speciesNameVector.push_back(specName1);
400 iIndex = modelSpeciesMap[specName1];
402 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
403 const_iterator jIter = modelSpeciesMap.find(specName2);
404 if (jIter == modelSpeciesMap.end())
406 modelSpeciesMap[specName2] = index;
407 modelSpeciesCodeList_.push_back(index);
408 speciesNameVector.push_back(specName2);
417 jIndex = modelSpeciesMap[specName2];
420 if (iIndex >= jIndex)
422 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
426 indx = iIndex*
N + jIndex - (iIndex*iIndex + iIndex)/2;
428 cutoffs_[indx] = nextCutoff;
429 epsilons_[indx] = nextEpsilon;
430 sigmas_[indx] = nextSigma;
432 getNextDataLine(parameterFilePointers[0], nextLinePtr,
437 sprintf(nextLine,
"There are not values for like-like pairs of:");
438 for (
int i = 0; i <
N; i++)
440 if (cutoffs_[(i*
N + i - (i*i + i)/2)] == -1)
442 strcat(nextLine,
" ");
443 strcat(nextLine, (speciesNameVector[i].String()).c_str());
454 for (
int jIndex = 0; jIndex <
N; jIndex++)
456 jjIndex = (jIndex*
N + jIndex - (jIndex*jIndex + jIndex)/2);
457 for (
int iIndex = (jIndex+1) ; iIndex <
N; iIndex++)
459 indx = jIndex*
N + iIndex - (jIndex*jIndex + jIndex)/2;
460 if (cutoffs_[indx] == -1)
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;
476 void LennardJones612Implementation::getNextDataLine(
477 FILE*
const filePtr,
char* nextLinePtr,
int const maxSize,
482 if(fgets(nextLinePtr, maxSize, filePtr) == NULL)
487 while ((nextLinePtr[0] ==
' ' || nextLinePtr[0] ==
'\t') ||
488 (nextLinePtr[0] ==
'\n' || nextLinePtr[0] ==
'\r' ))
490 nextLinePtr = (nextLinePtr + 1);
493 while ((strncmp(
"#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
497 void LennardJones612Implementation::CloseParameterFiles(
498 int const numberParameterFiles,
501 for (
int i = 0; i < numberParameterFiles; ++i)
502 fclose(parameterFilePointers[i]);
507 int LennardJones612Implementation::ConvertUnits(
525 double convertLength = 1.0;
527 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
528 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
529 requestedTemperatureUnit, requestedTimeUnit,
530 1.0, 0.0, 0.0, 0.0, 0.0,
534 LOG_ERROR(
"Unable to convert length unit");
537 if (convertLength !=
ONE)
539 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
541 cutoffs_[i] *= convertLength;
542 sigmas_[i] *= convertLength;
546 double convertEnergy = 1.0;
548 fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
549 requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
550 requestedTemperatureUnit, requestedTimeUnit,
551 0.0, 1.0, 0.0, 0.0, 0.0,
555 LOG_ERROR(
"Unable to convert energy unit");
558 if (convertEnergy !=
ONE)
560 for (
int i = 0; i < numberUniqueSpeciesPairs_; ++i)
562 epsilons_[i] *= convertEnergy;
575 LOG_ERROR(
"Unable to set units to requested values");
585 int LennardJones612Implementation::RegisterKIMModelSettings(
597 int LennardJones612Implementation::RegisterKIMComputeArgumentsSettings(
635 int LennardJones612Implementation::RegisterKIMParameters(
648 numberUniqueSpeciesPairs_, cutoffs_,
"cutoffs");
655 numberUniqueSpeciesPairs_, epsilons_,
"epsilons");
662 numberUniqueSpeciesPairs_, sigmas_,
"sigmas");
675 int LennardJones612Implementation::RegisterKIMFunctions(
699 template<
class ModelObj>
700 int LennardJones612Implementation::SetRefreshMutableValues(
701 ModelObj *
const modelObj)
709 for (
int i = 0; i < numberModelSpecies_; ++i)
711 for (
int j = 0; j <= i ; ++j)
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];
734 influenceDistance_ = 0.0;
736 for (
int i = 0; i < numberModelSpecies_; i++)
738 int indexI = modelSpeciesCodeList_[i];
740 for (
int j = 0; j < numberModelSpecies_; j++)
742 int indexJ = modelSpeciesCodeList_[j];
744 if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
746 influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
751 influenceDistance_ = sqrt(influenceDistance_);
752 modelObj->SetInfluenceDistancePointer(&influenceDistance_);
753 modelObj->SetNeighborListCutoffsPointer(1, &influenceDistance_);
757 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
758 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
762 for (
int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
764 for(
int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
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;
784 int LennardJones612Implementation::SetComputeMutableValues(
786 bool& isComputeProcess_dEdr,
787 bool& isComputeProcess_d2Edr2,
788 bool& isComputeEnergy,
789 bool& isComputeForces,
790 bool& isComputeParticleEnergy,
791 bool& isComputeVirial,
792 bool& isComputeParticleVirial,
797 double*& particleEnergy,
805 int compProcess_dEdr;
806 int compProcess_d2Edr2;
813 &compProcess_d2Edr2);
815 isComputeProcess_dEdr = compProcess_dEdr;
816 isComputeProcess_d2Edr2 = compProcess_d2Edr2;
840 (
double const **
const) &forces)
843 (
double const **
const) &virial)
846 (
double const **
const) &particleVirial);
853 isComputeEnergy = (energy != 0);
854 isComputeParticleEnergy = (particleEnergy != 0);
855 isComputeForces = (forces != 0);
856 isComputeVirial = (virial != 0);
857 isComputeParticleVirial = (particleVirial != 0);
869 int LennardJones612Implementation::CheckParticleSpeciesCodes(
875 for (
int i = 0; i < cachedNumberOfParticles_; ++i)
881 LOG_ERROR(
"unsupported particle species codes detected");
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 903 const int processd2E = 2;
904 const int energy = 2;
906 const int particleEnergy = 2;
907 const int virial = 2;
908 const int particleVirial = 2;
915 index += (int(isComputeProcess_dEdr))
916 * processd2E * energy * force * particleEnergy * virial
917 * particleVirial* shift;
920 index += (int(isComputeProcess_d2Edr2))
921 * energy * force * particleEnergy * virial * particleVirial * shift;
924 index += (int(isComputeEnergy))
925 * force * particleEnergy * virial * particleVirial * shift;
928 index += (int(isComputeForces))
929 * particleEnergy * virial * particleVirial * shift;
932 index += (int(isComputeParticleEnergy))
933 * virial * particleVirial * shift;
936 index += (int(isComputeVirial))
937 * particleVirial * shift;
940 index += (int(isComputeParticleVirial))
944 index += (int(isShift));
950 void LennardJones612Implementation::ProcessVirialTerm(
953 const double*
const r_ij,
958 double const v = dEidr/rij;
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];
969 void LennardJones612Implementation::ProcessParticleVirialTerm(
972 const double*
const r_ij,
977 double const v = dEidr/rij;
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];
987 for (
int k = 0; k < 6; ++k)
989 particleVirial[i][k] += vir[k];
990 particleVirial[j][k] += vir[k];
1002 int const extentOne)
1004 arrayPtr =
new double*[extentZero];
1005 arrayPtr[0] =
new double[extentZero * extentOne];
1006 for (
int i = 1; i < extentZero; ++i)
1008 arrayPtr[i] = arrayPtr[i-1] + extentOne;
1012 for (
int i = 0; i < extentZero; ++i)
1014 for (
int j = 0; j < extentOne; ++j)
1016 arrayPtr[i][j] = 0.0;
1024 if (arrayPtr != 0)
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