28#ifndef LENNARD_JONES_612_IMPLEMENTATION_HPP_
29#define LENNARD_JONES_612_IMPLEMENTATION_HPP_
41#define MAX_PARAMETER_FILES 1
43#define PARAM_SHIFT_INDEX 0
44#define PARAM_CUTOFFS_INDEX 1
45#define PARAM_EPSILONS_INDEX 2
46#define PARAM_SIGMAS_INDEX 3
94 modelComputeArgumentsCreate)
const;
96 modelComputeArgumentsDestroy)
const;
105 int numberModelSpecies_;
106 std::vector<int> modelSpeciesCodeList_;
107 int numberUniqueSpeciesPairs_;
137 double influenceDistance_;
138 double ** cutoffsSq2D_;
139 int modelWillNotRequestNeighborsOfNoncontributingParticles_;
140 double ** fourEpsilonSigma6_2D_;
141 double ** fourEpsilonSigma12_2D_;
142 double ** twentyFourEpsilonSigma6_2D_;
143 double ** fortyEightEpsilonSigma12_2D_;
144 double ** oneSixtyEightEpsilonSigma6_2D_;
145 double ** sixTwentyFourEpsilonSigma12_2D_;
154 int cachedNumberOfParticles_;
161 void AllocatePrivateParameterMemory();
162 void AllocateParameterMemory();
165 int const numberParameterFiles,
167 int ProcessParameterFiles(
169 int const numberParameterFiles,
171 void getNextDataLine(FILE *
const filePtr,
172 char *
const nextLine,
174 int * endOfFileFlag);
176 CloseParameterFiles(
int const numberParameterFiles,
184 int RegisterKIMModelSettings(
186 int RegisterKIMComputeArgumentsSettings(
190 int RegisterKIMFunctions(
194 template<
class ModelObj>
195 int SetRefreshMutableValues(ModelObj *
const modelObj);
199 int SetComputeMutableValues(
201 bool & isComputeProcess_dEdr,
202 bool & isComputeProcess_d2Edr2,
203 bool & isComputeEnergy,
204 bool & isComputeForces,
205 bool & isComputeParticleEnergy,
206 bool & isComputeVirial,
207 bool & isComputeParticleVirial,
208 int const *& particleSpeciesCodes,
209 int const *& particleContributing,
212 double *& particleEnergy,
217 int const *
const particleSpeciesCodes)
const;
218 int GetComputeIndex(
const bool & isComputeProcess_dEdr,
219 const bool & isComputeProcess_d2Edr2,
220 const bool & isComputeEnergy,
221 const bool & isComputeForces,
222 const bool & isComputeParticleEnergy,
223 const bool & isComputeVirial,
224 const bool & isComputeParticleVirial,
225 const bool & isShift)
const;
226 void ProcessVirialTerm(
const double & dEidr,
228 const double *
const r_ij,
232 void ProcessParticleVirialTerm(
const double & dEidr,
234 const double *
const r_ij,
240 template<
bool isComputeProcess_dEdr,
241 bool isComputeProcess_d2Edr2,
242 bool isComputeEnergy,
243 bool isComputeForces,
244 bool isComputeParticleEnergy,
245 bool isComputeVirial,
246 bool isComputeParticleVirial,
250 const int *
const particleSpeciesCodes,
251 const int *
const particleContributing,
253 double *
const energy,
255 double *
const particleEnergy,
275#define LENNARD_JONES_PHI(exshift) \
277 * (constFourEpsSig12_2D[iSpecies][jSpecies] * r6iv \
278 - constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
281#undef KIM_LOGGER_OBJECT_NAME
282#define KIM_LOGGER_OBJECT_NAME modelCompute
284template<
bool isComputeProcess_dEdr,
285 bool isComputeProcess_d2Edr2,
286 bool isComputeEnergy,
287 bool isComputeForces,
288 bool isComputeParticleEnergy,
289 bool isComputeVirial,
290 bool isComputeParticleVirial,
295 const int *
const particleSpeciesCodes,
296 const int *
const particleContributing,
298 double *
const energy,
300 double *
const particleEnergy,
306 if ((isComputeEnergy ==
false) && (isComputeParticleEnergy ==
false)
307 && (isComputeForces ==
false) && (isComputeProcess_dEdr ==
false)
308 && (isComputeProcess_d2Edr2 ==
false) && (isComputeVirial ==
false)
309 && (isComputeParticleVirial ==
false))
313 if (isComputeEnergy ==
true) { *energy = 0.0; }
314 if (isComputeVirial ==
true)
316 for (
int i = 0; i < 6; ++i) virial[i] = 0.0;
318 if (isComputeParticleEnergy ==
true)
320 int const cachedNumParticles = cachedNumberOfParticles_;
321 for (
int i = 0; i < cachedNumParticles; ++i) { particleEnergy[i] = 0.0; }
323 if (isComputeForces ==
true)
325 int const cachedNumParticles = cachedNumberOfParticles_;
326 for (
int i = 0; i < cachedNumParticles; ++i)
328 for (
int j = 0; j <
DIMENSION; ++j) forces[i][j] = 0.0;
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) particleVirial[i][j] = 0.0;
345 int const * n1atom = NULL;
346 double const *
const *
const constCutoffsSq2D = cutoffsSq2D_;
347 double const *
const *
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
348 double const *
const *
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
349 double const *
const *
const constTwentyFourEpsSig6_2D
350 = twentyFourEpsilonSigma6_2D_;
351 double const *
const *
const constFortyEightEpsSig12_2D
352 = fortyEightEpsilonSigma12_2D_;
353 double const *
const *
const constOneSixtyEightEpsSig6_2D
354 = oneSixtyEightEpsilonSigma6_2D_;
355 double const *
const *
const constSixTwentyFourEpsSig12_2D
356 = sixTwentyFourEpsilonSigma12_2D_;
357 double const *
const *
const constShifts2D = shifts2D_;
358 for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
360 if (particleContributing[ii])
363 int const numNei = numnei;
364 int const *
const n1Atom = n1atom;
366 int const iSpecies = particleSpeciesCodes[i];
369 for (
int jj = 0; jj < numNei; ++jj)
371 int const j = n1Atom[jj];
372 int const jContrib = particleContributing[j];
374 if (!(jContrib && (j < i)))
376 int const jSpecies = particleSpeciesCodes[j];
382 r_ij[k] = coordinates[j][k] - coordinates[i][k];
383 double const *
const r_ij_const =
const_cast<double *
>(r_ij);
386 double const rij2 = r_ij_const[0] * r_ij_const[0]
387 + r_ij_const[1] * r_ij_const[1]
388 + r_ij_const[2] * r_ij_const[2];
390 if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
393 double dphiByR = 0.0;
395 double dEidrByR = 0.0;
396 double d2Eidr2 = 0.0;
397 double const r2iv = 1.0 / rij2;
398 double const r6iv = r2iv * r2iv * r2iv;
400 if (isComputeProcess_d2Edr2 ==
true)
404 * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies] * r6iv
405 - constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
407 if (jContrib == 1) { d2Eidr2 = d2phi; }
410 d2Eidr2 = 0.5 * d2phi;
414 if ((isComputeProcess_dEdr ==
true) || (isComputeForces ==
true)
415 || (isComputeVirial ==
true)
416 || (isComputeParticleVirial ==
true))
420 * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies]
421 - constFortyEightEpsSig12_2D[iSpecies][jSpecies] * r6iv)
423 if (jContrib == 1) { dEidrByR = dphiByR; }
426 dEidrByR = 0.5 * dphiByR;
430 if ((isComputeEnergy ==
true) || (isComputeParticleEnergy ==
true))
443 if (isComputeEnergy ==
true)
445 if (jContrib == 1) { *energy += phi; }
448 *energy += 0.5 * phi;
453 if (isComputeParticleEnergy ==
true)
455 double const halfPhi = 0.5 * phi;
456 particleEnergy[i] += halfPhi;
457 if (jContrib == 1) { particleEnergy[j] += halfPhi; }
461 if (isComputeForces ==
true)
465 double const contrib = dEidrByR * r_ij_const[k];
466 forces[i][k] += contrib;
467 forces[j][k] -= contrib;
472 if ((isComputeProcess_dEdr ==
true) || (isComputeVirial ==
true)
473 || (isComputeParticleVirial ==
true))
475 double const rij = sqrt(rij2);
476 double const dEidr = dEidrByR * rij;
478 if (isComputeProcess_dEdr ==
true)
481 dEidr, rij, r_ij_const, i, j);
489 if (isComputeVirial ==
true)
491 ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
494 if (isComputeParticleVirial ==
true)
496 ProcessParticleVirialTerm(
497 dEidr, rij, r_ij_const, i, j, particleVirial);
502 if (isComputeProcess_d2Edr2 ==
true)
504 double const rij = sqrt(rij2);
505 double const R_pairs[2] = {rij, rij};
506 double const *
const pRs = &R_pairs[0];
507 double const Rij_pairs[6] = {r_ij_const[0],
513 double const *
const pRijConsts = &Rij_pairs[0];
514 int const i_pairs[2] = {i, i};
515 int const j_pairs[2] = {j, j};
516 int const *
const pis = &i_pairs[0];
517 int const *
const pjs = &j_pairs[0];
520 d2Eidr2, pRs, pRijConsts, pis, pjs);
#define LOG_ERROR(message)
Convenience macro for ERROR Log entries with compile-time optimization.
int() GetNeighborFunction(void const *const, int const, int *const, int const **const)
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
double VectorOfSizeDIM[DIMENSION]
#define MAX_PARAMETER_FILES
void Deallocate2DArray(double **&arrayPtr)
#define LENNARD_JONES_PHI(exshift)
double VectorOfSizeSix[6]
An Extensible Enumeration for the ChargeUnit's supported by the KIM API.
An Extensible Enumeration for the EnergyUnit's supported by the KIM API.
An Extensible Enumeration for the LengthUnit's supported by the KIM API.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
int ProcessD2EDr2Term(double const de, double const *const r, double const *const dx, int const *const i, int const *const j) const
Call the Simulator's COMPUTE_CALLBACK_NAME::ProcessD2EDr2Term routine.
int GetNeighborList(int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle) const
Get the neighbor list for a particle of interest corresponding to a particular neighbor list cutoff d...
int ProcessDEDrTerm(double const de, double const r, double const *const dx, int const i, int const j) const
Call the Simulator's COMPUTE_CALLBACK_NAME::ProcessDEDrTerm routine.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::R...
An Extensible Enumeration for the TemperatureUnit's supported by the KIM API.
An Extensible Enumeration for the TimeUnit's supported by the KIM API.
~LennardJones612Implementation()
int Refresh(KIM::ModelRefresh *const modelRefresh)
int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const