kim-api 2.3.0+AppleClang.AppleClang.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
Loading...
Searching...
No Matches
LennardJones612Implementation.hpp
Go to the documentation of this file.
1//
2// KIM-API: An API for interatomic models
3// Copyright (c) 2013--2022, Regents of the University of Minnesota.
4// All rights reserved.
5//
6// Contributors:
7// Ryan S. Elliott
8// Andrew Akerson
9//
10// SPDX-License-Identifier: LGPL-2.1-or-later
11//
12// This library is free software; you can redistribute it and/or
13// modify it under the terms of the GNU Lesser General Public
14// License as published by the Free Software Foundation; either
15// version 2.1 of the License, or (at your option) any later version.
16//
17// This library is distributed in the hope that it will be useful,
18// but WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20// Lesser General Public License for more details.
21//
22// You should have received a copy of the GNU Lesser General Public License
23// along with this library; if not, write to the Free Software Foundation,
24// Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25//
26
27
28#ifndef LENNARD_JONES_612_IMPLEMENTATION_HPP_
29#define LENNARD_JONES_612_IMPLEMENTATION_HPP_
30
31#include "KIM_LogMacros.hpp"
33#include <cmath>
34#include <cstdio>
35#include <vector>
36
37#define DIMENSION 3
38#define ONE 1.0
39#define HALF 0.5
40
41#define MAX_PARAMETER_FILES 1
42
43#define PARAM_SHIFT_INDEX 0
44#define PARAM_CUTOFFS_INDEX 1
45#define PARAM_EPSILONS_INDEX 2
46#define PARAM_SIGMAS_INDEX 3
47
48
49//==============================================================================
50//
51// Type definitions, enumerations, and helper function prototypes
52//
53//==============================================================================
54
55// type declaration for get neighbor functions
56typedef int(GetNeighborFunction)(void const * const,
57 int const,
58 int * const,
59 int const ** const);
60// type declaration for vector of constant dimension
61typedef double VectorOfSizeDIM[DIMENSION];
62typedef double VectorOfSizeSix[6];
63
64// helper routine declarations
65void AllocateAndInitialize2DArray(double **& arrayPtr,
66 int const extentZero,
67 int const extentOne);
68void Deallocate2DArray(double **& arrayPtr);
69
70//==============================================================================
71//
72// Declaration of LennardJones612Implementation class
73//
74//==============================================================================
75
76//******************************************************************************
78{
79 public:
81 KIM::ModelDriverCreate * const modelDriverCreate,
82 KIM::LengthUnit const requestedLengthUnit,
83 KIM::EnergyUnit const requestedEnergyUnit,
84 KIM::ChargeUnit const requestedChargeUnit,
85 KIM::TemperatureUnit const requestedTemperatureUnit,
86 KIM::TimeUnit const requestedTimeUnit,
87 int * const ier);
88 ~LennardJones612Implementation(); // no explicit Destroy() needed here
89
90 int Refresh(KIM::ModelRefresh * const modelRefresh);
91 int Compute(KIM::ModelCompute const * const modelCompute,
92 KIM::ModelComputeArguments const * const modelComputeArguments);
94 modelComputeArgumentsCreate) const;
96 modelComputeArgumentsDestroy) const;
97
98
99 private:
100 // Constant values that never change
101 // Set in constructor (via SetConstantValues)
102 //
103 //
104 // LennardJones612Implementation: constants
105 int numberModelSpecies_;
106 std::vector<int> modelSpeciesCodeList_;
107 int numberUniqueSpeciesPairs_;
108
109
110 // Constant values that are read from the input files and never change
111 // Set in constructor (via functions listed below)
112 //
113 //
114 // Private Model Parameters
115 // Memory allocated in AllocatePrivateParameterMemory() (from constructor)
116 // Memory deallocated in destructor
117 // Data set in ReadParameterFile routines
118 // none
119 //
120 // KIM API: Model Parameters whose (pointer) values never change
121 // Memory allocated in AllocateParameterMemory() (from constructor)
122 // Memory deallocated in destructor
123 // Data set in ReadParameterFile routines OR by KIM Simulator
124 int shift_;
125 double * cutoffs_;
126 double * epsilons_;
127 double * sigmas_;
128
129 // Mutable values that only change when Refresh() executes
130 // Set in Refresh (via SetRefreshMutableValues)
131 //
132 //
133 // KIM API: Model Parameters (can be changed directly by KIM Simulator)
134 // none
135 //
136 // LennardJones612Implementation: values (changed only by Refresh())
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_;
146 double ** shifts2D_;
147
148
149 // Mutable values that can change with each call to Refresh() and Compute()
150 // Memory may be reallocated on each call
151 //
152 //
153 // LennardJones612Implementation: values that change
154 int cachedNumberOfParticles_;
155
156
157 // Helper methods
158 //
159 //
160 // Related to constructor
161 void AllocatePrivateParameterMemory();
162 void AllocateParameterMemory();
163 static int
164 OpenParameterFiles(KIM::ModelDriverCreate * const modelDriverCreate,
165 int const numberParameterFiles,
166 FILE * parameterFilePointers[MAX_PARAMETER_FILES]);
167 int ProcessParameterFiles(
168 KIM::ModelDriverCreate * const modelDriverCreate,
169 int const numberParameterFiles,
170 FILE * const parameterFilePointers[MAX_PARAMETER_FILES]);
171 void getNextDataLine(FILE * const filePtr,
172 char * const nextLine,
173 int const maxSize,
174 int * endOfFileFlag);
175 static void
176 CloseParameterFiles(int const numberParameterFiles,
177 FILE * const parameterFilePointers[MAX_PARAMETER_FILES]);
178 int ConvertUnits(KIM::ModelDriverCreate * const modelDriverCreate,
179 KIM::LengthUnit const requestedLengthUnit,
180 KIM::EnergyUnit const requestedEnergyUnit,
181 KIM::ChargeUnit const requestedChargeUnit,
182 KIM::TemperatureUnit const requestedTemperatureUnit,
183 KIM::TimeUnit const requestedTimeUnit);
184 int RegisterKIMModelSettings(
185 KIM::ModelDriverCreate * const modelDriverCreate) const;
186 int RegisterKIMComputeArgumentsSettings(
187 KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
188 const;
189 int RegisterKIMParameters(KIM::ModelDriverCreate * const modelDriverCreate);
190 int RegisterKIMFunctions(
191 KIM::ModelDriverCreate * const modelDriverCreate) const;
192 //
193 // Related to Refresh()
194 template<class ModelObj>
195 int SetRefreshMutableValues(ModelObj * const modelObj);
196
197 //
198 // Related to Compute()
199 int SetComputeMutableValues(
200 KIM::ModelComputeArguments const * const modelComputeArguments,
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,
210 VectorOfSizeDIM const *& coordinates,
211 double *& energy,
212 double *& particleEnergy,
213 VectorOfSizeDIM *& forces,
214 VectorOfSizeSix *& virial,
215 VectorOfSizeSix *& particleViral);
216 int CheckParticleSpeciesCodes(KIM::ModelCompute const * const modelCompute,
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,
227 const double & rij,
228 const double * const r_ij,
229 const int & i,
230 const int & j,
231 VectorOfSizeSix virial) const;
232 void ProcessParticleVirialTerm(const double & dEidr,
233 const double & rij,
234 const double * const r_ij,
235 const int & i,
236 const int & j,
237 VectorOfSizeSix * const particleVirial) const;
238
239 // compute functions
240 template<bool isComputeProcess_dEdr,
241 bool isComputeProcess_d2Edr2,
242 bool isComputeEnergy,
243 bool isComputeForces,
244 bool isComputeParticleEnergy,
245 bool isComputeVirial,
246 bool isComputeParticleVirial,
247 bool isShift>
248 int Compute(KIM::ModelCompute const * const modelCompute,
249 KIM::ModelComputeArguments const * const modelComputeArguments,
250 const int * const particleSpeciesCodes,
251 const int * const particleContributing,
252 const VectorOfSizeDIM * const coordinates,
253 double * const energy,
254 VectorOfSizeDIM * const forces,
255 double * const particleEnergy,
256 VectorOfSizeSix virial,
257 VectorOfSizeSix * const particleVirial) const;
258};
259
260//==============================================================================
261//
262// Definition of LennardJones612Implementation::Compute functions
263//
264// NOTE: Here we rely on the compiler optimizations to prune dead code
265// after the template expansions. This provides high efficiency
266// and easy maintenance.
267//
268//==============================================================================
269
270//******************************************************************************
271// MACRO to compute Lennard-Jones phi
272// (used for efficiency)
273//
274// exshift - expression to be added to the end of the phi value
275#define LENNARD_JONES_PHI(exshift) \
276 phi = r6iv \
277 * (constFourEpsSig12_2D[iSpecies][jSpecies] * r6iv \
278 - constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
279
280//******************************************************************************
281#undef KIM_LOGGER_OBJECT_NAME
282#define KIM_LOGGER_OBJECT_NAME modelCompute
283//
284template<bool isComputeProcess_dEdr,
285 bool isComputeProcess_d2Edr2,
286 bool isComputeEnergy,
287 bool isComputeForces,
288 bool isComputeParticleEnergy,
289 bool isComputeVirial,
290 bool isComputeParticleVirial,
291 bool isShift>
293 KIM::ModelCompute const * const modelCompute,
294 KIM::ModelComputeArguments const * const modelComputeArguments,
295 const int * const particleSpeciesCodes,
296 const int * const particleContributing,
297 const VectorOfSizeDIM * const coordinates,
298 double * const energy,
299 VectorOfSizeDIM * const forces,
300 double * const particleEnergy,
301 VectorOfSizeSix virial,
302 VectorOfSizeSix * const particleVirial) const
303{
304 int ier = false;
305
306 if ((isComputeEnergy == false) && (isComputeParticleEnergy == false)
307 && (isComputeForces == false) && (isComputeProcess_dEdr == false)
308 && (isComputeProcess_d2Edr2 == false) && (isComputeVirial == false)
309 && (isComputeParticleVirial == false))
310 return ier;
311
312 // initialize energy and forces
313 if (isComputeEnergy == true) { *energy = 0.0; }
314 if (isComputeVirial == true)
315 {
316 for (int i = 0; i < 6; ++i) virial[i] = 0.0;
317 }
318 if (isComputeParticleEnergy == true)
319 {
320 int const cachedNumParticles = cachedNumberOfParticles_;
321 for (int i = 0; i < cachedNumParticles; ++i) { particleEnergy[i] = 0.0; }
322 }
323 if (isComputeForces == true)
324 {
325 int const cachedNumParticles = cachedNumberOfParticles_;
326 for (int i = 0; i < cachedNumParticles; ++i)
327 {
328 for (int j = 0; j < DIMENSION; ++j) forces[i][j] = 0.0;
329 }
330 }
331 if (isComputeParticleVirial == true)
332 {
333 int const cachedNumParticles = cachedNumberOfParticles_;
334 for (int i = 0; i < cachedNumParticles; ++i)
335 {
336 for (int j = 0; j < 6; ++j) particleVirial[i][j] = 0.0;
337 }
338 }
339
340 // calculate contribution from pair function
341 //
342 // Setup loop over contributing particles
343 int ii = 0;
344 int numnei = 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)
359 {
360 if (particleContributing[ii])
361 {
362 modelComputeArguments->GetNeighborList(0, ii, &numnei, &n1atom);
363 int const numNei = numnei;
364 int const * const n1Atom = n1atom;
365 int const i = ii;
366 int const iSpecies = particleSpeciesCodes[i];
367
368 // Setup loop over neighbors of current particle
369 for (int jj = 0; jj < numNei; ++jj)
370 {
371 int const j = n1Atom[jj];
372 int const jContrib = particleContributing[j];
373
374 if (!(jContrib && (j < i))) // effective half-list
375 {
376 int const jSpecies = particleSpeciesCodes[j];
377 double * r_ij;
378 double r_ijValue[DIMENSION];
379 // Compute r_ij
380 r_ij = r_ijValue;
381 for (int k = 0; k < DIMENSION; ++k)
382 r_ij[k] = coordinates[j][k] - coordinates[i][k];
383 double const * const r_ij_const = const_cast<double *>(r_ij);
384
385 // compute distance squared
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];
389
390 if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
391 { // compute contribution to energy, force, etc.
392 double phi = 0.0;
393 double dphiByR = 0.0;
394 double d2phi = 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;
399 // Compute pair potential and its derivatives
400 if (isComputeProcess_d2Edr2 == true)
401 { // Compute d2phi
402 d2phi
403 = r6iv
404 * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies] * r6iv
405 - constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
406 * r2iv;
407 if (jContrib == 1) { d2Eidr2 = d2phi; }
408 else
409 {
410 d2Eidr2 = 0.5 * d2phi;
411 }
412 }
413
414 if ((isComputeProcess_dEdr == true) || (isComputeForces == true)
415 || (isComputeVirial == true)
416 || (isComputeParticleVirial == true))
417 { // Compute dphi
418 dphiByR
419 = r6iv
420 * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies]
421 - constFortyEightEpsSig12_2D[iSpecies][jSpecies] * r6iv)
422 * r2iv;
423 if (jContrib == 1) { dEidrByR = dphiByR; }
424 else
425 {
426 dEidrByR = 0.5 * dphiByR;
427 }
428 }
429
430 if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
431 { // Compute phi
432 if (isShift == true)
433 {
434 LENNARD_JONES_PHI(-constShifts2D[iSpecies][jSpecies]);
435 }
436 else
437 {
439 }
440 }
441
442 // Contribution to energy
443 if (isComputeEnergy == true)
444 {
445 if (jContrib == 1) { *energy += phi; }
446 else
447 {
448 *energy += 0.5 * phi;
449 }
450 }
451
452 // Contribution to particleEnergy
453 if (isComputeParticleEnergy == true)
454 {
455 double const halfPhi = 0.5 * phi;
456 particleEnergy[i] += halfPhi;
457 if (jContrib == 1) { particleEnergy[j] += halfPhi; }
458 }
459
460 // Contribution to forces
461 if (isComputeForces == true)
462 {
463 for (int k = 0; k < DIMENSION; ++k)
464 {
465 double const contrib = dEidrByR * r_ij_const[k];
466 forces[i][k] += contrib;
467 forces[j][k] -= contrib;
468 }
469 }
470
471 // Call process_dEdr
472 if ((isComputeProcess_dEdr == true) || (isComputeVirial == true)
473 || (isComputeParticleVirial == true))
474 {
475 double const rij = sqrt(rij2);
476 double const dEidr = dEidrByR * rij;
477
478 if (isComputeProcess_dEdr == true)
479 {
480 ier = modelComputeArguments->ProcessDEDrTerm(
481 dEidr, rij, r_ij_const, i, j);
482 if (ier)
483 {
484 LOG_ERROR("process_dEdr");
485 return ier;
486 }
487 }
488
489 if (isComputeVirial == true)
490 {
491 ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
492 }
493
494 if (isComputeParticleVirial == true)
495 {
496 ProcessParticleVirialTerm(
497 dEidr, rij, r_ij_const, i, j, particleVirial);
498 }
499 }
500
501 // Call process_d2Edr2
502 if (isComputeProcess_d2Edr2 == true)
503 {
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],
508 r_ij_const[1],
509 r_ij_const[2],
510 r_ij_const[0],
511 r_ij_const[1],
512 r_ij_const[2]};
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];
518
519 ier = modelComputeArguments->ProcessD2EDr2Term(
520 d2Eidr2, pRs, pRijConsts, pis, pjs);
521 if (ier)
522 {
523 LOG_ERROR("process_d2Edr2");
524 return ier;
525 }
526 }
527 } // if particleContributing
528 } // if i < j or j non contributing
529 } // if particles i and j interact
530 } // end of first neighbor loop
531 } // end of loop over contributing particles
532
533 // everything is good
534 ier = false;
535 return ier;
536}
537
538#endif // LENNARD_JONES_612_IMPLEMENTATION_HPP_
#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.
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