kim-api  2.1.4-git+v2.1.3-git-1-g7847914a.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
LennardJones612Implementation.hpp
Go to the documentation of this file.
1 //
2 // CDDL HEADER START
3 //
4 // The contents of this file are subject to the terms of the Common Development
5 // and Distribution License Version 1.0 (the "License").
6 //
7 // You can obtain a copy of the license at
8 // http://www.opensource.org/licenses/CDDL-1.0. See the License for the
9 // specific language governing permissions and limitations under the License.
10 //
11 // When distributing Covered Code, include this CDDL HEADER in each file and
12 // include the License file in a prominent location with the name LICENSE.CDDL.
13 // If applicable, add the following below this CDDL HEADER, with the fields
14 // enclosed by brackets "[]" replaced with your own identifying information:
15 //
16 // Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
17 //
18 // CDDL HEADER END
19 //
20 
21 //
22 // Copyright (c) 2015, Regents of the University of Minnesota.
23 // All rights reserved.
24 //
25 // Contributors:
26 // Ryan S. Elliott
27 // Andrew Akerson
28 //
29 
30 
31 #ifndef LENNARD_JONES_612_IMPLEMENTATION_HPP_
32 #define LENNARD_JONES_612_IMPLEMENTATION_HPP_
33 
34 #include "KIM_LogMacros.hpp"
35 #include "LennardJones612.hpp"
36 #include <cmath>
37 #include <vector>
38 
39 #define DIMENSION 3
40 #define ONE 1.0
41 #define HALF 0.5
42 
43 #define MAX_PARAMETER_FILES 1
44 
45 #define PARAM_SHIFT_INDEX 0
46 #define PARAM_CUTOFFS_INDEX 1
47 #define PARAM_EPSILONS_INDEX 2
48 #define PARAM_SIGMAS_INDEX 3
49 
50 
51 //==============================================================================
52 //
53 // Type definitions, enumerations, and helper function prototypes
54 //
55 //==============================================================================
56 
57 // type declaration for get neighbor functions
58 typedef int(GetNeighborFunction)(void const * const,
59  int const,
60  int * const,
61  int const ** const);
62 // type declaration for vector of constant dimension
63 typedef double VectorOfSizeDIM[DIMENSION];
64 typedef double VectorOfSizeSix[6];
65 
66 // helper routine declarations
67 void AllocateAndInitialize2DArray(double **& arrayPtr,
68  int const extentZero,
69  int const extentOne);
70 void Deallocate2DArray(double **& arrayPtr);
71 
72 //==============================================================================
73 //
74 // Declaration of LennardJones612Implementation class
75 //
76 //==============================================================================
77 
78 //******************************************************************************
80 {
81  public:
83  KIM::ModelDriverCreate * const modelDriverCreate,
84  KIM::LengthUnit const requestedLengthUnit,
85  KIM::EnergyUnit const requestedEnergyUnit,
86  KIM::ChargeUnit const requestedChargeUnit,
87  KIM::TemperatureUnit const requestedTemperatureUnit,
88  KIM::TimeUnit const requestedTimeUnit,
89  int * const ier);
90  ~LennardJones612Implementation(); // no explicit Destroy() needed here
91 
92  int Refresh(KIM::ModelRefresh * const modelRefresh);
93  int Compute(KIM::ModelCompute const * const modelCompute,
94  KIM::ModelComputeArguments const * const modelComputeArguments);
96  modelComputeArgumentsCreate) const;
98  modelComputeArgumentsDestroy) const;
99 
100 
101  private:
102  // Constant values that never change
103  // Set in constructor (via SetConstantValues)
104  //
105  //
106  // LennardJones612Implementation: constants
107  int numberModelSpecies_;
108  std::vector<int> modelSpeciesCodeList_;
109  int numberUniqueSpeciesPairs_;
110 
111 
112  // Constant values that are read from the input files and never change
113  // Set in constructor (via functions listed below)
114  //
115  //
116  // Private Model Parameters
117  // Memory allocated in AllocatePrivateParameterMemory() (from constructor)
118  // Memory deallocated in destructor
119  // Data set in ReadParameterFile routines
120  // none
121  //
122  // KIM API: Model Parameters whose (pointer) values never change
123  // Memory allocated in AllocateParameterMemory() (from constructor)
124  // Memory deallocated in destructor
125  // Data set in ReadParameterFile routines OR by KIM Simulator
126  int shift_;
127  double * cutoffs_;
128  double * epsilons_;
129  double * sigmas_;
130 
131  // Mutable values that only change when Refresh() executes
132  // Set in Refresh (via SetRefreshMutableValues)
133  //
134  //
135  // KIM API: Model Parameters (can be changed directly by KIM Simulator)
136  // none
137  //
138  // LennardJones612Implementation: values (changed only by Refresh())
139  double influenceDistance_;
140  double ** cutoffsSq2D_;
141  int modelWillNotRequestNeighborsOfNoncontributingParticles_;
142  double ** fourEpsilonSigma6_2D_;
143  double ** fourEpsilonSigma12_2D_;
144  double ** twentyFourEpsilonSigma6_2D_;
145  double ** fortyEightEpsilonSigma12_2D_;
146  double ** oneSixtyEightEpsilonSigma6_2D_;
147  double ** sixTwentyFourEpsilonSigma12_2D_;
148  double ** shifts2D_;
149 
150 
151  // Mutable values that can change with each call to Refresh() and Compute()
152  // Memory may be reallocated on each call
153  //
154  //
155  // LennardJones612Implementation: values that change
156  int cachedNumberOfParticles_;
157 
158 
159  // Helper methods
160  //
161  //
162  // Related to constructor
163  void AllocatePrivateParameterMemory();
164  void AllocateParameterMemory();
165  static int
166  OpenParameterFiles(KIM::ModelDriverCreate * const modelDriverCreate,
167  int const numberParameterFiles,
168  FILE * parameterFilePointers[MAX_PARAMETER_FILES]);
169  int ProcessParameterFiles(
170  KIM::ModelDriverCreate * const modelDriverCreate,
171  int const numberParameterFiles,
172  FILE * const parameterFilePointers[MAX_PARAMETER_FILES]);
173  void getNextDataLine(FILE * const filePtr,
174  char * const nextLine,
175  int const maxSize,
176  int * endOfFileFlag);
177  static void
178  CloseParameterFiles(int const numberParameterFiles,
179  FILE * const parameterFilePointers[MAX_PARAMETER_FILES]);
180  int ConvertUnits(KIM::ModelDriverCreate * const modelDriverCreate,
181  KIM::LengthUnit const requestedLengthUnit,
182  KIM::EnergyUnit const requestedEnergyUnit,
183  KIM::ChargeUnit const requestedChargeUnit,
184  KIM::TemperatureUnit const requestedTemperatureUnit,
185  KIM::TimeUnit const requestedTimeUnit);
186  int RegisterKIMModelSettings(
187  KIM::ModelDriverCreate * const modelDriverCreate) const;
188  int RegisterKIMComputeArgumentsSettings(
189  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
190  const;
191  int RegisterKIMParameters(KIM::ModelDriverCreate * const modelDriverCreate);
192  int RegisterKIMFunctions(
193  KIM::ModelDriverCreate * const modelDriverCreate) const;
194  //
195  // Related to Refresh()
196  template<class ModelObj>
197  int SetRefreshMutableValues(ModelObj * const modelObj);
198 
199  //
200  // Related to Compute()
201  int SetComputeMutableValues(
202  KIM::ModelComputeArguments const * const modelComputeArguments,
203  bool & isComputeProcess_dEdr,
204  bool & isComputeProcess_d2Edr2,
205  bool & isComputeEnergy,
206  bool & isComputeForces,
207  bool & isComputeParticleEnergy,
208  bool & isComputeVirial,
209  bool & isComputeParticleVirial,
210  int const *& particleSpeciesCodes,
211  int const *& particleContributing,
212  VectorOfSizeDIM const *& coordinates,
213  double *& energy,
214  double *& particleEnergy,
215  VectorOfSizeDIM *& forces,
216  VectorOfSizeSix *& virial,
217  VectorOfSizeSix *& particleViral);
218  int CheckParticleSpeciesCodes(KIM::ModelCompute const * const modelCompute,
219  int const * const particleSpeciesCodes) const;
220  int GetComputeIndex(const bool & isComputeProcess_dEdr,
221  const bool & isComputeProcess_d2Edr2,
222  const bool & isComputeEnergy,
223  const bool & isComputeForces,
224  const bool & isComputeParticleEnergy,
225  const bool & isComputeVirial,
226  const bool & isComputeParticleVirial,
227  const bool & isShift) const;
228  void ProcessVirialTerm(const double & dEidr,
229  const double & rij,
230  const double * const r_ij,
231  const int & i,
232  const int & j,
233  VectorOfSizeSix virial) const;
234  void ProcessParticleVirialTerm(const double & dEidr,
235  const double & rij,
236  const double * const r_ij,
237  const int & i,
238  const int & j,
239  VectorOfSizeSix * const particleVirial) const;
240 
241  // compute functions
242  template<bool isComputeProcess_dEdr,
243  bool isComputeProcess_d2Edr2,
244  bool isComputeEnergy,
245  bool isComputeForces,
246  bool isComputeParticleEnergy,
247  bool isComputeVirial,
248  bool isComputeParticleVirial,
249  bool isShift>
250  int Compute(KIM::ModelCompute const * const modelCompute,
251  KIM::ModelComputeArguments const * const modelComputeArguments,
252  const int * const particleSpeciesCodes,
253  const int * const particleContributing,
254  const VectorOfSizeDIM * const coordinates,
255  double * const energy,
256  VectorOfSizeDIM * const forces,
257  double * const particleEnergy,
258  VectorOfSizeSix virial,
259  VectorOfSizeSix * const particleVirial) const;
260 };
261 
262 //==============================================================================
263 //
264 // Definition of LennardJones612Implementation::Compute functions
265 //
266 // NOTE: Here we rely on the compiler optimizations to prune dead code
267 // after the template expansions. This provides high efficiency
268 // and easy maintenance.
269 //
270 //==============================================================================
271 
272 //******************************************************************************
273 // MACRO to compute Lennard-Jones phi
274 // (used for efficiency)
275 //
276 // exshift - expression to be added to the end of the phi value
277 #define LENNARD_JONES_PHI(exshift) \
278  phi = r6iv \
279  * (constFourEpsSig12_2D[iSpecies][jSpecies] * r6iv \
280  - constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
281 
282 //******************************************************************************
283 #undef KIM_LOGGER_OBJECT_NAME
284 #define KIM_LOGGER_OBJECT_NAME modelCompute
285 //
286 template<bool isComputeProcess_dEdr,
287  bool isComputeProcess_d2Edr2,
288  bool isComputeEnergy,
289  bool isComputeForces,
290  bool isComputeParticleEnergy,
291  bool isComputeVirial,
292  bool isComputeParticleVirial,
293  bool isShift>
295  KIM::ModelCompute const * const modelCompute,
296  KIM::ModelComputeArguments const * const modelComputeArguments,
297  const int * const particleSpeciesCodes,
298  const int * const particleContributing,
299  const VectorOfSizeDIM * const coordinates,
300  double * const energy,
301  VectorOfSizeDIM * const forces,
302  double * const particleEnergy,
303  VectorOfSizeSix virial,
304  VectorOfSizeSix * const particleVirial) const
305 {
306  int ier = false;
307 
308  if ((isComputeEnergy == false) && (isComputeParticleEnergy == false)
309  && (isComputeForces == false) && (isComputeProcess_dEdr == false)
310  && (isComputeProcess_d2Edr2 == false) && (isComputeVirial == false)
311  && (isComputeParticleVirial == false))
312  return ier;
313 
314  // initialize energy and forces
315  if (isComputeEnergy == true) { *energy = 0.0; }
316  if (isComputeVirial == true)
317  {
318  for (int i = 0; i < 6; ++i) virial[i] = 0.0;
319  }
320  if (isComputeParticleEnergy == true)
321  {
322  int const cachedNumParticles = cachedNumberOfParticles_;
323  for (int i = 0; i < cachedNumParticles; ++i) { particleEnergy[i] = 0.0; }
324  }
325  if (isComputeForces == true)
326  {
327  int const cachedNumParticles = cachedNumberOfParticles_;
328  for (int i = 0; i < cachedNumParticles; ++i)
329  {
330  for (int j = 0; j < DIMENSION; ++j) forces[i][j] = 0.0;
331  }
332  }
333  if (isComputeParticleVirial == true)
334  {
335  int const cachedNumParticles = cachedNumberOfParticles_;
336  for (int i = 0; i < cachedNumParticles; ++i)
337  {
338  for (int j = 0; j < 6; ++j) particleVirial[i][j] = 0.0;
339  }
340  }
341 
342  // calculate contribution from pair function
343  //
344  // Setup loop over contributing particles
345  int ii = 0;
346  int numnei = 0;
347  int const * n1atom = NULL;
348  double const * const * const constCutoffsSq2D = cutoffsSq2D_;
349  double const * const * const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
350  double const * const * const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
351  double const * const * const constTwentyFourEpsSig6_2D
352  = twentyFourEpsilonSigma6_2D_;
353  double const * const * const constFortyEightEpsSig12_2D
354  = fortyEightEpsilonSigma12_2D_;
355  double const * const * const constOneSixtyEightEpsSig6_2D
356  = oneSixtyEightEpsilonSigma6_2D_;
357  double const * const * const constSixTwentyFourEpsSig12_2D
358  = sixTwentyFourEpsilonSigma12_2D_;
359  double const * const * const constShifts2D = shifts2D_;
360  for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
361  {
362  if (particleContributing[ii])
363  {
364  modelComputeArguments->GetNeighborList(0, ii, &numnei, &n1atom);
365  int const numNei = numnei;
366  int const * const n1Atom = n1atom;
367  int const i = ii;
368  int const iSpecies = particleSpeciesCodes[i];
369 
370  // Setup loop over neighbors of current particle
371  for (int jj = 0; jj < numNei; ++jj)
372  {
373  int const j = n1Atom[jj];
374  int const jContrib = particleContributing[j];
375 
376  if (!(jContrib && (j < i))) // effective half-list
377  {
378  int const jSpecies = particleSpeciesCodes[j];
379  double * r_ij;
380  double r_ijValue[DIMENSION];
381  // Compute r_ij
382  r_ij = r_ijValue;
383  for (int k = 0; k < DIMENSION; ++k)
384  r_ij[k] = coordinates[j][k] - coordinates[i][k];
385  double const * const r_ij_const = const_cast<double *>(r_ij);
386 
387  // compute distance squared
388  double const rij2 = r_ij_const[0] * r_ij_const[0]
389  + r_ij_const[1] * r_ij_const[1]
390  + r_ij_const[2] * r_ij_const[2];
391 
392  if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
393  { // compute contribution to energy, force, etc.
394  double phi = 0.0;
395  double dphiByR = 0.0;
396  double d2phi = 0.0;
397  double dEidrByR = 0.0;
398  double d2Eidr2 = 0.0;
399  double const r2iv = 1.0 / rij2;
400  double const r6iv = r2iv * r2iv * r2iv;
401  // Compute pair potential and its derivatives
402  if (isComputeProcess_d2Edr2 == true)
403  { // Compute d2phi
404  d2phi
405  = r6iv
406  * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies] * r6iv
407  - constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
408  * r2iv;
409  if (jContrib == 1) { d2Eidr2 = d2phi; }
410  else
411  {
412  d2Eidr2 = 0.5 * d2phi;
413  }
414  }
415 
416  if ((isComputeProcess_dEdr == true) || (isComputeForces == true)
417  || (isComputeVirial == true)
418  || (isComputeParticleVirial == true))
419  { // Compute dphi
420  dphiByR
421  = r6iv
422  * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies]
423  - constFortyEightEpsSig12_2D[iSpecies][jSpecies] * r6iv)
424  * r2iv;
425  if (jContrib == 1) { dEidrByR = dphiByR; }
426  else
427  {
428  dEidrByR = 0.5 * dphiByR;
429  }
430  }
431 
432  if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
433  { // Compute phi
434  if (isShift == true)
435  { LENNARD_JONES_PHI(-constShifts2D[iSpecies][jSpecies]); }
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  { ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial); }
491 
492  if (isComputeParticleVirial == true)
493  {
494  ProcessParticleVirialTerm(
495  dEidr, rij, r_ij_const, i, j, particleVirial);
496  }
497  }
498 
499  // Call process_d2Edr2
500  if (isComputeProcess_d2Edr2 == true)
501  {
502  double const rij = sqrt(rij2);
503  double const R_pairs[2] = {rij, rij};
504  double const * const pRs = &R_pairs[0];
505  double const Rij_pairs[6] = {r_ij_const[0],
506  r_ij_const[1],
507  r_ij_const[2],
508  r_ij_const[0],
509  r_ij_const[1],
510  r_ij_const[2]};
511  double const * const pRijConsts = &Rij_pairs[0];
512  int const i_pairs[2] = {i, i};
513  int const j_pairs[2] = {j, j};
514  int const * const pis = &i_pairs[0];
515  int const * const pjs = &j_pairs[0];
516 
517  ier = modelComputeArguments->ProcessD2EDr2Term(
518  d2Eidr2, pRs, pRijConsts, pis, pjs);
519  if (ier)
520  {
521  LOG_ERROR("process_d2Edr2");
522  return ier;
523  }
524  }
525  } // if particleContributing
526  } // if i < j or j non contributing
527  } // if particles i and j interact
528  } // end of first neighbor loop
529  } // end of loop over contributing particles
530 
531  // everything is good
532  ier = false;
533  return ier;
534 }
535 
536 #endif // LENNARD_JONES_612_IMPLEMENTATION_HPP_
int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
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&#39;s supported by the KIM API.
An Extensible Enumeration for the TimeUnit&#39;s supported by the KIM API.
#define MAX_PARAMETER_FILES
ComputeArgumentName const coordinates
The standard coordinates argument.
#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)
int ProcessDEDrTerm(double const de, double const r, double const *const dx, int const i, int const j) const
Call the Simulator&#39;s COMPUTE_CALLBACK_NAME::ProcessDEDrTerm routine.
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&#39;s COMPUTE_CALLBACK_NAME::ProcessD2EDr2Term routine.
double VectorOfSizeDIM[DIMENSION]
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
An Extensible Enumeration for the LengthUnit&#39;s supported by the KIM API.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
ComputeArgumentName const particleContributing
The standard particleContributing argument.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
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)
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
An Extensible Enumeration for the EnergyUnit&#39;s supported by the KIM API.
int Refresh(KIM::ModelRefresh *const modelRefresh)
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
ComputeArgumentName const particleSpeciesCodes
The standard particleSpeciesCodes argument.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
#define LENNARD_JONES_PHI(exshift)
An Extensible Enumeration for the ChargeUnit&#39;s supported by the KIM API.
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...
double VectorOfSizeSix[6]
int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)