KIM API V2
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 <vector>
35 #include <cmath>
36 #include "KIM_LogVerbosity.hpp"
37 #include "LennardJones612.hpp"
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, int const,
59  int * const, int const ** const);
60 // type declaration for vector of constant dimension
61 typedef double VectorOfSizeDIM[DIMENSION];
62 typedef double VectorOfSizeSix[6];
63 
64 // helper routine declarations
65 void AllocateAndInitialize2DArray(double**& arrayPtr, int const extentZero,
66  int const extentOne);
67 void Deallocate2DArray(double**& arrayPtr);
68 
69 //==============================================================================
70 //
71 // Declaration of LennardJones612Implementation class
72 //
73 //==============================================================================
74 
75 //******************************************************************************
77 {
78  public:
80  KIM::ModelDriverCreate * const modelDriverCreate,
81  KIM::LengthUnit const requestedLengthUnit,
82  KIM::EnergyUnit const requestedEnergyUnit,
83  KIM::ChargeUnit const requestedChargeUnit,
84  KIM::TemperatureUnit const requestedTemperatureUnit,
85  KIM::TimeUnit const requestedTimeUnit,
86  int * const ier);
87  ~LennardJones612Implementation(); // no explicit Destroy() needed here
88 
89  int Refresh(KIM::ModelRefresh * const modelRefresh);
90  int Compute(
91  KIM::ModelCompute const * const modelCompute,
92  KIM::ModelComputeArguments const * const modelComputeArguments);
94  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
95  const;
97  KIM::ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
98  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 paddingNeighborHints_;
142  int halfListHints_;
143  double** fourEpsilonSigma6_2D_;
144  double** fourEpsilonSigma12_2D_;
145  double** twentyFourEpsilonSigma6_2D_;
146  double** fortyEightEpsilonSigma12_2D_;
147  double** oneSixtyEightEpsilonSigma6_2D_;
148  double** sixTwentyFourEpsilonSigma12_2D_;
149  double** shifts2D_;
150 
151 
152  // Mutable values that can change with each call to Refresh() and Compute()
153  // Memory may be reallocated on each call
154  //
155  //
156  // LennardJones612Implementation: values that change
157  int cachedNumberOfParticles_;
158 
159 
160  // Helper methods
161  //
162  //
163  // Related to constructor
164  void AllocatePrivateParameterMemory();
165  void AllocateParameterMemory();
166  static int OpenParameterFiles(
167  KIM::ModelDriverCreate * const modelDriverCreate,
168  int const numberParameterFiles,
169  FILE* parameterFilePointers[MAX_PARAMETER_FILES]);
170  int ProcessParameterFiles(
171  KIM::ModelDriverCreate * const modelDriverCreate,
172  int const numberParameterFiles,
173  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
174  void getNextDataLine(FILE* const filePtr, char* const nextLine,
175  int const maxSize, int* endOfFileFlag);
176  static void CloseParameterFiles(
177  int const numberParameterFiles,
178  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
179  int ConvertUnits(
180  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(
192  KIM::ModelDriverCreate * const modelDriverCreate);
193  int RegisterKIMFunctions(
194  KIM::ModelDriverCreate * const modelDriverCreate) const;
195  //
196  // Related to Refresh()
197  template<class ModelObj>
198  int SetRefreshMutableValues(ModelObj * const modelObj);
199 
200  //
201  // Related to Compute()
202  int SetComputeMutableValues(
203  KIM::ModelComputeArguments const * const modelComputeArguments,
204  bool& isComputeProcess_dEdr,
205  bool& isComputeProcess_d2Edr2,
206  bool& isComputeEnergy,
207  bool& isComputeForces,
208  bool& isComputeParticleEnergy,
209  bool& isComputeVirial,
210  bool& isComputeParticleVirial,
211  int const*& particleSpeciesCodes,
212  int const*& particleContributing,
214  double*& energy,
215  double*& particleEnergy,
216  VectorOfSizeDIM*& forces,
217  VectorOfSizeSix*& virial,
218  VectorOfSizeSix*& particleViral);
219  int CheckParticleSpeciesCodes(KIM::ModelCompute const * const modelCompute,
220  int const* const particleSpeciesCodes) const;
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,
230  const double& rij,
231  const double* const r_ij,
232  const int& i,
233  const int& j,
234  VectorOfSizeSix virial) const;
235  void ProcessParticleVirialTerm(const double& dEidr,
236  const double& rij,
237  const double* const r_ij,
238  const int& i,
239  const int& j,
240  VectorOfSizeSix* const particleVirial) const;
241 
242  // compute functions
243  template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
244  bool isComputeEnergy, bool isComputeForces,
245  bool isComputeParticleEnergy, bool isComputeVirial,
246  bool isComputeParticleVirial, bool isShift >
247  int Compute(KIM::ModelCompute const * const modelCompute,
248  KIM::ModelComputeArguments const * const modelComputeArguments,
249  const int* const particleSpeciesCodes,
250  const int* const particleContributing,
251  const VectorOfSizeDIM* const coordinates,
252  double* const energy,
253  VectorOfSizeDIM* const forces,
254  double* const particleEnergy,
255  VectorOfSizeSix virial,
256  VectorOfSizeSix* const particleVirial) const;
257 };
258 
259 //==============================================================================
260 //
261 // Definition of LennardJones612Implementation::Compute functions
262 //
263 // NOTE: Here we rely on the compiler optimizations to prune dead code
264 // after the template expansions. This provides high efficiency
265 // and easy maintenance.
266 //
267 //==============================================================================
268 
269 //******************************************************************************
270 // MACRO to compute Lennard-Jones phi
271 // (used for efficiency)
272 //
273 // exshift - expression to be added to the end of the phi value
274 #define LENNARD_JONES_PHI(exshift) \
275  phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \
276  constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
277 
278 //******************************************************************************
280 template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
281  bool isComputeEnergy, bool isComputeForces,
282  bool isComputeParticleEnergy, bool isComputeVirial,
283  bool isComputeParticleVirial, bool isShift >
285  KIM::ModelCompute const * const modelCompute,
286  KIM::ModelComputeArguments const * const modelComputeArguments,
287  const int* const particleSpeciesCodes,
288  const int* const particleContributing,
289  const VectorOfSizeDIM* const coordinates,
290  double* const energy,
291  VectorOfSizeDIM* const forces,
292  double* const particleEnergy,
293  VectorOfSizeSix virial,
294  VectorOfSizeSix* const particleVirial) const
295 {
296  int ier = false;
297 
298  if ((isComputeEnergy == false) &&
299  (isComputeParticleEnergy == false) &&
300  (isComputeForces == false) &&
301  (isComputeProcess_dEdr == false) &&
302  (isComputeProcess_d2Edr2 == false) &&
303  (isComputeVirial == false) &&
304  (isComputeParticleVirial == false))
305  return ier;
306 
307  // initialize energy and forces
308  if (isComputeEnergy == true)
309  {
310  *energy = 0.0;
311  }
312  if (isComputeVirial == true)
313  {
314  for (int i = 0; i < 6; ++i) virial[i] = 0.0;
315  }
316  if (isComputeParticleEnergy == true)
317  {
318  int const cachedNumParticles = cachedNumberOfParticles_;
319  for (int i = 0; i < cachedNumParticles; ++i)
320  {
321  particleEnergy[i] = 0.0;
322  }
323  }
324  if (isComputeForces == true)
325  {
326  int const cachedNumParticles = cachedNumberOfParticles_;
327  for (int i = 0; i < cachedNumParticles; ++i)
328  {
329  for (int j = 0; j < DIMENSION; ++j)
330  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)
339  particleVirial[i][j] = 0.0;
340  }
341  }
342 
343  // calculate contribution from pair function
344  //
345  // Setup loop over contributing particles
346  int ii = 0;
347  int numnei = 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)
362  {
363  if (particleContributing[ii])
364  {
365  modelComputeArguments->GetNeighborList(0, ii, &numnei, &n1atom);
366  int const numNei = numnei;
367  int const * const n1Atom = n1atom;
368  int const i = ii;
369  int const iSpecies = particleSpeciesCodes[i];
370 
371  // Setup loop over neighbors of current particle
372  for (int jj = 0; jj < numNei; ++jj)
373  {
374  int const j = n1Atom[jj];
375 
376  if (i < j) // 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 =
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];
392 
393  if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
394  { // compute contribution to energy, force, etc.
395  double phi = 0.0;
396  double dphiByR = 0.0;
397  double d2phi = 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;
402  // Compute pair potential and its derivatives
403  if (isComputeProcess_d2Edr2 == true)
404  { // Compute d2phi
405  d2phi =
406  r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]
407  *r6iv -
408  constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
409  * r2iv;
410  if (particleContributing[j] == 1)
411  {
412  d2Eidr2 = d2phi;
413  }
414  else
415  {
416  d2Eidr2 = 0.5*d2phi;
417  }
418  }
419 
420  if ((isComputeProcess_dEdr == true) || (isComputeForces == true) ||
421  (isComputeVirial == true) || (isComputeParticleVirial == true))
422  { // Compute dphi
423  dphiByR =
424  r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
425  constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
426  * r2iv;
427  if (particleContributing[j] == 1)
428  {
429  dEidrByR = dphiByR;
430  }
431  else
432  {
433  dEidrByR = 0.5*dphiByR;
434  }
435  }
436 
437  if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
438  { // Compute phi
439  if (isShift == true)
440  {
441  LENNARD_JONES_PHI(- constShifts2D[iSpecies][jSpecies]);
442  }
443  else
444  {
446  }
447  }
448 
449  // Contribution to energy
450  if (isComputeEnergy == true)
451  {
452  if (particleContributing[j] == 1)
453  {
454  *energy += phi;
455  }
456  else
457  {
458  *energy += 0.5*phi;
459  }
460  }
461 
462  // Contribution to particleEnergy
463  if (isComputeParticleEnergy == true)
464  {
465  double const halfPhi = 0.5*phi;
466  particleEnergy[i] += halfPhi;
467  if (particleContributing[j] == 1)
468  {
469  particleEnergy[j] += halfPhi;
470  }
471  }
472 
473  // Contribution to forces
474  if (isComputeForces == true)
475  {
476  for (int k = 0; k < DIMENSION; ++k)
477  {
478  double const contrib = dEidrByR * r_ij_const[k];
479  forces[i][k] += contrib;
480  forces[j][k] -= contrib;
481  }
482  }
483 
484  // Call process_dEdr
485  if ((isComputeProcess_dEdr == true) ||
486  (isComputeVirial == true) ||
487  (isComputeParticleVirial == true))
488  {
489  double const rij = sqrt(rij2);
490  double const dEidr = dEidrByR*rij;
491 
492  if (isComputeProcess_dEdr == true)
493  {
494  ier = modelComputeArguments
495  ->ProcessDEDrTerm(dEidr, rij, r_ij_const, i, j);
496  if (ier)
497  {
498  LOG_ERROR("process_dEdr");
499  return ier;
500  }
501  }
502 
503  if (isComputeVirial == true)
504  {
505  ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
506  }
507 
508  if (isComputeParticleVirial == true)
509  {
510  ProcessParticleVirialTerm(dEidr, rij, r_ij_const, i, j,
511  particleVirial);
512  }
513  }
514 
515  // Call process_d2Edr2
516  if (isComputeProcess_d2Edr2 == true)
517  {
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];
529 
530  ier = modelComputeArguments
531  ->ProcessD2EDr2Term(d2Eidr2, pRs, pRijConsts, pis, pjs);
532  if (ier)
533  {
534  LOG_ERROR("process_d2Edr2");
535  return ier;
536  }
537  }
538  } // if particleContributing
539  } // if i < j
540  } // if particles i and j interact
541  } // end of first neighbor loop
542  } // end of loop over contributing particles
543 
544  // everything is good
545  ier = false;
546  return ier;
547 }
548 
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)
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)