31 #ifndef LENNARD_JONES_612_IMPLEMENTATION_HPP_ 32 #define LENNARD_JONES_612_IMPLEMENTATION_HPP_ 43 #define MAX_PARAMETER_FILES 1 45 #define PARAM_SHIFT_INDEX 0 46 #define PARAM_CUTOFFS_INDEX 1 47 #define PARAM_EPSILONS_INDEX 2 48 #define PARAM_SIGMAS_INDEX 3 59 int *
const,
int const **
const);
107 int numberModelSpecies_;
108 std::vector<int> modelSpeciesCodeList_;
109 int numberUniqueSpeciesPairs_;
139 double influenceDistance_;
140 double** cutoffsSq2D_;
141 double** fourEpsilonSigma6_2D_;
142 double** fourEpsilonSigma12_2D_;
143 double** twentyFourEpsilonSigma6_2D_;
144 double** fortyEightEpsilonSigma12_2D_;
145 double** oneSixtyEightEpsilonSigma6_2D_;
146 double** sixTwentyFourEpsilonSigma12_2D_;
155 int cachedNumberOfParticles_;
162 void AllocatePrivateParameterMemory();
163 void AllocateParameterMemory();
164 static int OpenParameterFiles(
166 int const numberParameterFiles,
168 int ProcessParameterFiles(
170 int const numberParameterFiles,
172 void getNextDataLine(FILE*
const filePtr,
char*
const nextLine,
173 int const maxSize,
int* endOfFileFlag);
174 static void CloseParameterFiles(
175 int const numberParameterFiles,
184 int RegisterKIMModelSettings(
186 int RegisterKIMComputeArgumentsSettings(
189 int RegisterKIMParameters(
191 int RegisterKIMFunctions(
195 template<
class ModelObj>
196 int SetRefreshMutableValues(ModelObj *
const modelObj);
200 int SetComputeMutableValues(
202 bool& isComputeProcess_dEdr,
203 bool& isComputeProcess_d2Edr2,
204 bool& isComputeEnergy,
205 bool& isComputeForces,
206 bool& isComputeParticleEnergy,
207 bool& isComputeVirial,
208 bool& isComputeParticleVirial,
213 double*& particleEnergy,
219 int GetComputeIndex(
const bool& isComputeProcess_dEdr,
220 const bool& isComputeProcess_d2Edr2,
221 const bool& isComputeEnergy,
222 const bool& isComputeForces,
223 const bool& isComputeParticleEnergy,
224 const bool& isComputeVirial,
225 const bool& isComputeParticleVirial,
226 const bool& isShift)
const;
227 void ProcessVirialTerm(
const double& dEidr,
229 const double*
const r_ij,
233 void ProcessParticleVirialTerm(
const double& dEidr,
235 const double*
const r_ij,
241 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
242 bool isComputeEnergy,
bool isComputeForces,
243 bool isComputeParticleEnergy,
bool isComputeVirial,
244 bool isComputeParticleVirial,
bool isShift >
250 double*
const energy,
252 double*
const particleEnergy,
272 #define LENNARD_JONES_PHI(exshift) \ 273 phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \ 274 constFourEpsSig6_2D[iSpecies][jSpecies]) exshift; 278 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
279 bool isComputeEnergy,
bool isComputeForces,
280 bool isComputeParticleEnergy,
bool isComputeVirial,
281 bool isComputeParticleVirial,
bool isShift >
288 double*
const energy,
290 double*
const particleEnergy,
296 if ((isComputeEnergy ==
false) &&
297 (isComputeParticleEnergy ==
false) &&
298 (isComputeForces ==
false) &&
299 (isComputeProcess_dEdr ==
false) &&
300 (isComputeProcess_d2Edr2 ==
false) &&
301 (isComputeVirial ==
false) &&
302 (isComputeParticleVirial ==
false))
306 if (isComputeEnergy ==
true)
310 if (isComputeVirial ==
true)
312 for (
int i = 0; i < 6; ++i) virial[i] = 0.0;
314 if (isComputeParticleEnergy ==
true)
316 int const cachedNumParticles = cachedNumberOfParticles_;
317 for (
int i = 0; i < cachedNumParticles; ++i)
319 particleEnergy[i] = 0.0;
322 if (isComputeForces ==
true)
324 int const cachedNumParticles = cachedNumberOfParticles_;
325 for (
int i = 0; i < cachedNumParticles; ++i)
331 if (isComputeParticleVirial ==
true)
333 int const cachedNumParticles = cachedNumberOfParticles_;
334 for (
int i = 0; i < cachedNumParticles; ++i)
336 for (
int j = 0; j < 6; ++j)
337 particleVirial[i][j] = 0.0;
346 int const * n1atom = 0;
347 double const*
const*
const constCutoffsSq2D = cutoffsSq2D_;
348 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
349 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
350 double const*
const*
const constTwentyFourEpsSig6_2D
351 = twentyFourEpsilonSigma6_2D_;
352 double const*
const*
const constFortyEightEpsSig12_2D
353 = fortyEightEpsilonSigma12_2D_;
354 double const*
const*
const constOneSixtyEightEpsSig6_2D
355 = oneSixtyEightEpsilonSigma6_2D_;
356 double const*
const*
const constSixTwentyFourEpsSig12_2D
357 = sixTwentyFourEpsilonSigma12_2D_;
358 double const*
const*
const constShifts2D = shifts2D_;
359 for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
364 int const numNei = numnei;
365 int const *
const n1Atom = n1atom;
370 for (
int jj = 0; jj < numNei; ++jj)
372 int const j = n1Atom[jj];
380 double const*
const r_ij_const =
const_cast<double*
>(r_ij);
384 r_ij_const[0] * r_ij_const[0] +
385 r_ij_const[1] * r_ij_const[1] +
386 r_ij_const[2] * r_ij_const[2];
388 if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
391 double dphiByR = 0.0;
393 double dEidrByR = 0.0;
394 double d2Eidr2 = 0.0;
395 double const r2iv = 1.0/rij2;
396 double const r6iv = r2iv*r2iv*r2iv;
398 if (isComputeProcess_d2Edr2 ==
true)
401 r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]*r6iv -
402 constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
407 if ((isComputeProcess_dEdr ==
true) || (isComputeForces ==
true) ||
408 (isComputeVirial ==
true) || (isComputeParticleVirial ==
true))
411 r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
412 constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
414 dEidrByR = 0.5*dphiByR;
417 if ((isComputeEnergy ==
true) || (isComputeParticleEnergy ==
true))
430 if (isComputeEnergy ==
true)
436 if (isComputeParticleEnergy ==
true)
438 double const halfPhi = 0.5*phi;
439 particleEnergy[i] += halfPhi;
443 if (isComputeForces ==
true)
447 double const contrib = dEidrByR * r_ij_const[k];
448 forces[i][k] += contrib;
449 forces[j][k] -= contrib;
454 if ((isComputeProcess_dEdr ==
true) ||
455 (isComputeVirial ==
true) ||
456 (isComputeParticleVirial ==
true))
458 double const rij = sqrt(rij2);
459 double const dEidr = dEidrByR*rij;
461 if (isComputeProcess_dEdr ==
true)
463 ier = modelComputeArguments
472 if (isComputeVirial ==
true)
474 ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
477 if (isComputeParticleVirial ==
true)
479 ProcessParticleVirialTerm(dEidr, rij, r_ij_const, i, j,
485 if (isComputeProcess_d2Edr2 ==
true)
487 double const rij = sqrt(rij2);
488 double const R_pairs[2] = {rij, rij};
489 double const*
const pRs = &R_pairs[0];
490 double const Rij_pairs[6]
491 = {r_ij_const[0], r_ij_const[1], r_ij_const[2],
492 r_ij_const[0], r_ij_const[1], r_ij_const[2]};
493 double const*
const pRijConsts = &Rij_pairs[0];
494 int const i_pairs[2] = {i, i};
495 int const j_pairs[2] = {j, j};
496 int const*
const pis = &i_pairs[0];
497 int const*
const pjs = &j_pairs[0];
499 ier = modelComputeArguments
517 #endif // LENNARD_JONES_612_IMPLEMENTATION_HPP_ int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
#define MAX_PARAMETER_FILES
ComputeArgumentName const coordinates
int() GetNeighborFunction(void const *const, int const, int *const, int const **const)
int ProcessDEDrTerm(double const de, double const r, double const *const dx, int const i, int const j) const
int ProcessD2EDr2Term(double const de, double const *const r, double const *const dx, int const *const i, int const *const j) const
double VectorOfSizeDIM[DIMENSION]
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
ComputeArgumentName const particleContributing
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)
void Deallocate2DArray(double **&arrayPtr)
#define LOG_ERROR(message)
int Refresh(KIM::ModelRefresh *const modelRefresh)
~LennardJones612Implementation()
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
ComputeArgumentName const particleSpeciesCodes
#define LENNARD_JONES_PHI(exshift)
int GetNeighborList(int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle) const
double VectorOfSizeSix[6]
int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)