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 int paddingNeighborHints_;
143 double** fourEpsilonSigma6_2D_;
144 double** fourEpsilonSigma12_2D_;
145 double** twentyFourEpsilonSigma6_2D_;
146 double** fortyEightEpsilonSigma12_2D_;
147 double** oneSixtyEightEpsilonSigma6_2D_;
148 double** sixTwentyFourEpsilonSigma12_2D_;
157 int cachedNumberOfParticles_;
164 void AllocatePrivateParameterMemory();
165 void AllocateParameterMemory();
166 static int OpenParameterFiles(
168 int const numberParameterFiles,
170 int ProcessParameterFiles(
172 int const numberParameterFiles,
174 void getNextDataLine(FILE*
const filePtr,
char*
const nextLine,
175 int const maxSize,
int* endOfFileFlag);
176 static void CloseParameterFiles(
177 int const numberParameterFiles,
186 int RegisterKIMModelSettings(
188 int RegisterKIMComputeArgumentsSettings(
191 int RegisterKIMParameters(
193 int RegisterKIMFunctions(
197 template<
class ModelObj>
198 int SetRefreshMutableValues(ModelObj *
const modelObj);
202 int SetComputeMutableValues(
204 bool& isComputeProcess_dEdr,
205 bool& isComputeProcess_d2Edr2,
206 bool& isComputeEnergy,
207 bool& isComputeForces,
208 bool& isComputeParticleEnergy,
209 bool& isComputeVirial,
210 bool& isComputeParticleVirial,
215 double*& particleEnergy,
221 int GetComputeIndex(
const bool& isComputeProcess_dEdr,
222 const bool& isComputeProcess_d2Edr2,
223 const bool& isComputeEnergy,
224 const bool& isComputeForces,
225 const bool& isComputeParticleEnergy,
226 const bool& isComputeVirial,
227 const bool& isComputeParticleVirial,
228 const bool& isShift)
const;
229 void ProcessVirialTerm(
const double& dEidr,
231 const double*
const r_ij,
235 void ProcessParticleVirialTerm(
const double& dEidr,
237 const double*
const r_ij,
243 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
244 bool isComputeEnergy,
bool isComputeForces,
245 bool isComputeParticleEnergy,
bool isComputeVirial,
246 bool isComputeParticleVirial,
bool isShift >
252 double*
const energy,
254 double*
const particleEnergy,
274 #define LENNARD_JONES_PHI(exshift) \ 275 phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \ 276 constFourEpsSig6_2D[iSpecies][jSpecies]) exshift; 280 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
281 bool isComputeEnergy,
bool isComputeForces,
282 bool isComputeParticleEnergy,
bool isComputeVirial,
283 bool isComputeParticleVirial,
bool isShift >
290 double*
const energy,
292 double*
const particleEnergy,
298 if ((isComputeEnergy ==
false) &&
299 (isComputeParticleEnergy ==
false) &&
300 (isComputeForces ==
false) &&
301 (isComputeProcess_dEdr ==
false) &&
302 (isComputeProcess_d2Edr2 ==
false) &&
303 (isComputeVirial ==
false) &&
304 (isComputeParticleVirial ==
false))
308 if (isComputeEnergy ==
true)
312 if (isComputeVirial ==
true)
314 for (
int i = 0; i < 6; ++i) virial[i] = 0.0;
316 if (isComputeParticleEnergy ==
true)
318 int const cachedNumParticles = cachedNumberOfParticles_;
319 for (
int i = 0; i < cachedNumParticles; ++i)
321 particleEnergy[i] = 0.0;
324 if (isComputeForces ==
true)
326 int const cachedNumParticles = cachedNumberOfParticles_;
327 for (
int i = 0; i < cachedNumParticles; ++i)
333 if (isComputeParticleVirial ==
true)
335 int const cachedNumParticles = cachedNumberOfParticles_;
336 for (
int i = 0; i < cachedNumParticles; ++i)
338 for (
int j = 0; j < 6; ++j)
339 particleVirial[i][j] = 0.0;
348 int const * n1atom = NULL;
349 double const*
const*
const constCutoffsSq2D = cutoffsSq2D_;
350 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
351 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
352 double const*
const*
const constTwentyFourEpsSig6_2D
353 = twentyFourEpsilonSigma6_2D_;
354 double const*
const*
const constFortyEightEpsSig12_2D
355 = fortyEightEpsilonSigma12_2D_;
356 double const*
const*
const constOneSixtyEightEpsSig6_2D
357 = oneSixtyEightEpsilonSigma6_2D_;
358 double const*
const*
const constSixTwentyFourEpsSig12_2D
359 = sixTwentyFourEpsilonSigma12_2D_;
360 double const*
const*
const constShifts2D = shifts2D_;
361 for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
366 int const numNei = numnei;
367 int const *
const n1Atom = n1atom;
372 for (
int jj = 0; jj < numNei; ++jj)
374 int const j = n1Atom[jj];
385 double const*
const r_ij_const =
const_cast<double*
>(r_ij);
389 r_ij_const[0] * r_ij_const[0] +
390 r_ij_const[1] * r_ij_const[1] +
391 r_ij_const[2] * r_ij_const[2];
393 if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
396 double dphiByR = 0.0;
398 double dEidrByR = 0.0;
399 double d2Eidr2 = 0.0;
400 double const r2iv = 1.0/rij2;
401 double const r6iv = r2iv*r2iv*r2iv;
403 if (isComputeProcess_d2Edr2 ==
true)
406 r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]
408 constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
420 if ((isComputeProcess_dEdr ==
true) || (isComputeForces ==
true) ||
421 (isComputeVirial ==
true) || (isComputeParticleVirial ==
true))
424 r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
425 constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
433 dEidrByR = 0.5*dphiByR;
437 if ((isComputeEnergy ==
true) || (isComputeParticleEnergy ==
true))
450 if (isComputeEnergy ==
true)
463 if (isComputeParticleEnergy ==
true)
465 double const halfPhi = 0.5*phi;
466 particleEnergy[i] += halfPhi;
469 particleEnergy[j] += halfPhi;
474 if (isComputeForces ==
true)
478 double const contrib = dEidrByR * r_ij_const[k];
479 forces[i][k] += contrib;
480 forces[j][k] -= contrib;
485 if ((isComputeProcess_dEdr ==
true) ||
486 (isComputeVirial ==
true) ||
487 (isComputeParticleVirial ==
true))
489 double const rij = sqrt(rij2);
490 double const dEidr = dEidrByR*rij;
492 if (isComputeProcess_dEdr ==
true)
494 ier = modelComputeArguments
503 if (isComputeVirial ==
true)
505 ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
508 if (isComputeParticleVirial ==
true)
510 ProcessParticleVirialTerm(dEidr, rij, r_ij_const, i, j,
516 if (isComputeProcess_d2Edr2 ==
true)
518 double const rij = sqrt(rij2);
519 double const R_pairs[2] = {rij, rij};
520 double const*
const pRs = &R_pairs[0];
521 double const Rij_pairs[6]
522 = {r_ij_const[0], r_ij_const[1], r_ij_const[2],
523 r_ij_const[0], r_ij_const[1], r_ij_const[2]};
524 double const*
const pRijConsts = &Rij_pairs[0];
525 int const i_pairs[2] = {i, i};
526 int const j_pairs[2] = {j, j};
527 int const*
const pis = &i_pairs[0];
528 int const*
const pjs = &j_pairs[0];
530 ier = modelComputeArguments
549 #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)