31 #ifndef LENNARD_JONES_612_IMPLEMENTATION_HPP_ 32 #define LENNARD_JONES_612_IMPLEMENTATION_HPP_ 42 #define MAX_PARAMETER_FILES 1 44 #define PARAM_SHIFT_INDEX 0 45 #define PARAM_CUTOFFS_INDEX 1 46 #define PARAM_EPSILONS_INDEX 2 47 #define PARAM_SIGMAS_INDEX 3 58 int *
const,
int const **
const);
96 int numberModelSpecies_;
97 std::vector<int> modelSpeciesCodeList_;
98 int numberUniqueSpeciesPairs_;
131 double influenceDistance_;
132 double** cutoffsSq2D_;
133 double** fourEpsilonSigma6_2D_;
134 double** fourEpsilonSigma12_2D_;
135 double** twentyFourEpsilonSigma6_2D_;
136 double** fortyEightEpsilonSigma12_2D_;
137 double** oneSixtyEightEpsilonSigma6_2D_;
138 double** sixTwentyFourEpsilonSigma12_2D_;
147 int cachedNumberOfParticles_;
154 void AllocateFreeParameterMemory();
155 static int OpenParameterFiles(
157 int const numberParameterFiles,
159 static void CloseParameterFiles(
160 int const numberParameterFiles,
162 int ProcessParameterFiles(
164 int const numberParameterFiles,
166 void getNextDataLine(FILE*
const filePtr,
char*
const nextLine,
167 int const maxSize,
int* endOfFileFlag);
175 int RegisterKIMModelSettings(
177 int RegisterKIMParameters(
179 int RegisterKIMFunctions(
183 template<
class ModelObj>
184 int SetReinitMutableValues(ModelObj *
const modelObj);
189 bool& isComputeProcess_dEdr,
190 bool& isComputeProcess_d2Edr2,
191 bool& isComputeEnergy,
192 bool& isComputeForces,
193 bool& isComputeParticleEnergy,
198 double*& particleEnergy,
202 int GetComputeIndex(
const bool& isComputeProcess_dEdr,
203 const bool& isComputeProcess_d2Edr2,
204 const bool& isComputeEnergy,
205 const bool& isComputeForces,
206 const bool& isComputeParticleEnergy,
207 const bool& isShift)
const;
210 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
211 bool isComputeEnergy,
bool isComputeForces,
212 bool isComputeParticleEnergy,
bool isShift >
217 double*
const energy,
219 double*
const particleEnergy);
237 #define LENNARD_JONES_PHI(exshift) \ 238 phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \ 239 constFourEpsSig6_2D[iSpecies][jSpecies]) exshift; 243 template<
bool isComputeProcess_dEdr,
bool isComputeProcess_d2Edr2,
244 bool isComputeEnergy,
bool isComputeForces,
245 bool isComputeParticleEnergy,
bool isShift >
251 double*
const energy,
253 double*
const particleEnergy)
257 if ((isComputeEnergy ==
false) &&
258 (isComputeParticleEnergy ==
false) &&
259 (isComputeForces ==
false) &&
260 (isComputeProcess_dEdr ==
false) &&
261 (isComputeProcess_d2Edr2 ==
false))
265 if (isComputeEnergy ==
true)
269 if (isComputeParticleEnergy ==
true)
271 int const cachedNumParticles = cachedNumberOfParticles_;
272 for (
int i = 0; i < cachedNumParticles; ++i)
274 particleEnergy[i] = 0.0;
277 if (isComputeForces ==
true)
279 int const cachedNumParticles = cachedNumberOfParticles_;
280 for (
int i = 0; i < cachedNumParticles; ++i)
292 int const * n1atom = 0;
293 double const*
const*
const constCutoffsSq2D = cutoffsSq2D_;
294 double const*
const*
const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
295 double const*
const*
const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
296 double const*
const*
const constTwentyFourEpsSig6_2D
297 = twentyFourEpsilonSigma6_2D_;
298 double const*
const*
const constFortyEightEpsSig12_2D
299 = fortyEightEpsilonSigma12_2D_;
300 double const*
const*
const constOneSixtyEightEpsSig6_2D
301 = oneSixtyEightEpsilonSigma6_2D_;
302 double const*
const*
const constSixTwentyFourEpsSig12_2D
303 = sixTwentyFourEpsilonSigma12_2D_;
304 double const*
const*
const constShifts2D = shifts2D_;
305 for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
310 int const numNei = numnei;
311 int const *
const n1Atom = n1atom;
316 for (
int jj = 0; jj < numNei; ++jj)
318 int const j = n1Atom[jj];
326 double const*
const r_ij_const =
const_cast<double*
>(r_ij);
330 r_ij_const[0] * r_ij_const[0] +
331 r_ij_const[1] * r_ij_const[1] +
332 r_ij_const[2] * r_ij_const[2];
334 if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
337 double dphiByR = 0.0;
339 double dEidrByR = 0.0;
340 double d2Eidr2 = 0.0;
341 double const r2iv = 1.0/rij2;
342 double const r6iv = r2iv*r2iv*r2iv;
344 if (isComputeProcess_d2Edr2 ==
true)
347 r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]*r6iv -
348 constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
353 if ((isComputeProcess_dEdr ==
true) || (isComputeForces ==
true))
356 r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
357 constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
359 dEidrByR = 0.5*dphiByR;
362 if ((isComputeEnergy ==
true) || (isComputeParticleEnergy ==
true))
375 if (isComputeEnergy ==
true)
381 if (isComputeParticleEnergy ==
true)
383 double const halfPhi = 0.5*phi;
384 particleEnergy[i] += halfPhi;
388 if (isComputeForces ==
true)
392 double const contrib = dEidrByR * r_ij_const[k];
393 forces[i][k] += contrib;
394 forces[j][k] -= contrib;
399 if (isComputeProcess_dEdr ==
true)
401 double const rij = sqrt(rij2);
402 double const dEidr = dEidrByR*rij;
412 if (isComputeProcess_d2Edr2 ==
true)
414 double const rij = sqrt(rij2);
415 double const R_pairs[2] = {rij, rij};
416 double const*
const pRs = &R_pairs[0];
417 double const Rij_pairs[6]
418 = {r_ij_const[0], r_ij_const[1], r_ij_const[2],
419 r_ij_const[0], r_ij_const[1], r_ij_const[2]};
420 double const*
const pRijConsts = &Rij_pairs[0];
421 int const i_pairs[2] = {i, i};
422 int const j_pairs[2] = {j, j};
423 int const*
const pis = &i_pairs[0];
424 int const*
const pjs = &j_pairs[0];
444 #endif // LENNARD_JONES_612_IMPLEMENTATION_HPP_
int ProcessDEDrTerm(double const de, double const r, double const *const dx, int const i, int const j) const
int GetNeighborList(int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle) const
#define MAX_PARAMETER_FILES
ArgumentName const particleContributing
int() GetNeighborFunction(void const *const, int const, int *const, int const **const)
double VectorOfSizeDIM[DIMENSION]
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
ArgumentName const coordinates
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)
ArgumentName const particleSpeciesCodes
int Refresh(KIM::ModelRefresh *const modelRefresh)
int Compute(KIM::ModelCompute const *const modelCompute)
~LennardJones612Implementation()
int ProcessD2EDr2Term(double const de, double const *const r, double const *const dx, int const *const i, int const *const j) const
#define LENNARD_JONES_PHI(exshift)
#define LOG_ERROR(message)