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  double** fourEpsilonSigma6_2D_;
142  double** fourEpsilonSigma12_2D_;
143  double** twentyFourEpsilonSigma6_2D_;
144  double** fortyEightEpsilonSigma12_2D_;
145  double** oneSixtyEightEpsilonSigma6_2D_;
146  double** sixTwentyFourEpsilonSigma12_2D_;
147  double** shifts2D_;
148 
149 
150  // Mutable values that can change with each call to Refresh() and Compute()
151  // Memory may be reallocated on each call
152  //
153  //
154  // LennardJones612Implementation: values that change
155  int cachedNumberOfParticles_;
156 
157 
158  // Helper methods
159  //
160  //
161  // Related to constructor
162  void AllocatePrivateParameterMemory();
163  void AllocateParameterMemory();
164  static int OpenParameterFiles(
165  KIM::ModelDriverCreate * const modelDriverCreate,
166  int const numberParameterFiles,
167  FILE* parameterFilePointers[MAX_PARAMETER_FILES]);
168  int ProcessParameterFiles(
169  KIM::ModelDriverCreate * const modelDriverCreate,
170  int const numberParameterFiles,
171  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
172  void getNextDataLine(FILE* const filePtr, char* const nextLine,
173  int const maxSize, int* endOfFileFlag);
174  static void CloseParameterFiles(
175  int const numberParameterFiles,
176  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
177  int ConvertUnits(
178  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(
190  KIM::ModelDriverCreate * const modelDriverCreate);
191  int RegisterKIMFunctions(
192  KIM::ModelDriverCreate * const modelDriverCreate) const;
193  //
194  // Related to Refresh()
195  template<class ModelObj>
196  int SetRefreshMutableValues(ModelObj * const modelObj);
197 
198  //
199  // Related to Compute()
200  int SetComputeMutableValues(
201  KIM::ModelComputeArguments const * const modelComputeArguments,
202  bool& isComputeProcess_dEdr,
203  bool& isComputeProcess_d2Edr2,
204  bool& isComputeEnergy,
205  bool& isComputeForces,
206  bool& isComputeParticleEnergy,
207  bool& isComputeVirial,
208  bool& isComputeParticleVirial,
209  int const*& particleSpeciesCodes,
210  int const*& particleContributing,
212  double*& energy,
213  double*& particleEnergy,
214  VectorOfSizeDIM*& forces,
215  VectorOfSizeSix*& virial,
216  VectorOfSizeSix*& particleViral);
217  int CheckParticleSpeciesCodes(KIM::ModelCompute const * const modelCompute,
218  int const* const particleSpeciesCodes) const;
219  int GetComputeIndex(const bool& isComputeProcess_dEdr,
220  const bool& isComputeProcess_d2Edr2,
221  const bool& isComputeEnergy,
222  const bool& isComputeForces,
223  const bool& isComputeParticleEnergy,
224  const bool& isComputeVirial,
225  const bool& isComputeParticleVirial,
226  const bool& isShift) const;
227  void ProcessVirialTerm(const double& dEidr,
228  const double& rij,
229  const double* const r_ij,
230  const int& i,
231  const int& j,
232  VectorOfSizeSix virial) const;
233  void ProcessParticleVirialTerm(const double& dEidr,
234  const double& rij,
235  const double* const r_ij,
236  const int& i,
237  const int& j,
238  VectorOfSizeSix* const particleVirial) const;
239 
240  // compute functions
241  template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
242  bool isComputeEnergy, bool isComputeForces,
243  bool isComputeParticleEnergy, bool isComputeVirial,
244  bool isComputeParticleVirial, bool isShift >
245  int Compute(KIM::ModelCompute const * const modelCompute,
246  KIM::ModelComputeArguments const * const modelComputeArguments,
247  const int* const particleSpeciesCodes,
248  const int* const particleContributing,
249  const VectorOfSizeDIM* const coordinates,
250  double* const energy,
251  VectorOfSizeDIM* const forces,
252  double* const particleEnergy,
253  VectorOfSizeSix virial,
254  VectorOfSizeSix* const particleVirial) const;
255 };
256 
257 //==============================================================================
258 //
259 // Definition of LennardJones612Implementation::Compute functions
260 //
261 // NOTE: Here we rely on the compiler optimizations to prune dead code
262 // after the template expansions. This provides high efficiency
263 // and easy maintenance.
264 //
265 //==============================================================================
266 
267 //******************************************************************************
268 // MACRO to compute Lennard-Jones phi
269 // (used for efficiency)
270 //
271 // exshift - expression to be added to the end of the phi value
272 #define LENNARD_JONES_PHI(exshift) \
273  phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \
274  constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
275 
276 //******************************************************************************
278  template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
279  bool isComputeEnergy, bool isComputeForces,
280  bool isComputeParticleEnergy, bool isComputeVirial,
281  bool isComputeParticleVirial, bool isShift >
283  KIM::ModelCompute const * const modelCompute,
284  KIM::ModelComputeArguments const * const modelComputeArguments,
285  const int* const particleSpeciesCodes,
286  const int* const particleContributing,
287  const VectorOfSizeDIM* const coordinates,
288  double* const energy,
289  VectorOfSizeDIM* const forces,
290  double* const particleEnergy,
291  VectorOfSizeSix virial,
292  VectorOfSizeSix* const particleVirial) const
293 {
294  int ier = false;
295 
296  if ((isComputeEnergy == false) &&
297  (isComputeParticleEnergy == false) &&
298  (isComputeForces == false) &&
299  (isComputeProcess_dEdr == false) &&
300  (isComputeProcess_d2Edr2 == false) &&
301  (isComputeVirial == false) &&
302  (isComputeParticleVirial == false))
303  return ier;
304 
305  // initialize energy and forces
306  if (isComputeEnergy == true)
307  {
308  *energy = 0.0;
309  }
310  if (isComputeVirial == true)
311  {
312  for (int i = 0; i < 6; ++i) virial[i] = 0.0;
313  }
314  if (isComputeParticleEnergy == true)
315  {
316  int const cachedNumParticles = cachedNumberOfParticles_;
317  for (int i = 0; i < cachedNumParticles; ++i)
318  {
319  particleEnergy[i] = 0.0;
320  }
321  }
322  if (isComputeForces == true)
323  {
324  int const cachedNumParticles = cachedNumberOfParticles_;
325  for (int i = 0; i < cachedNumParticles; ++i)
326  {
327  for (int j = 0; j < DIMENSION; ++j)
328  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)
337  particleVirial[i][j] = 0.0;
338  }
339  }
340 
341  // calculate contribution from pair function
342  //
343  // Setup loop over contributing particles
344  int ii = 0;
345  int numnei = 0;
346  int const * n1atom = 0;
347  double const* const* const constCutoffsSq2D = cutoffsSq2D_;
348  double const* const* const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
349  double const* const* const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
350  double const* const* const constTwentyFourEpsSig6_2D
351  = twentyFourEpsilonSigma6_2D_;
352  double const* const* const constFortyEightEpsSig12_2D
353  = fortyEightEpsilonSigma12_2D_;
354  double const* const* const constOneSixtyEightEpsSig6_2D
355  = oneSixtyEightEpsilonSigma6_2D_;
356  double const* const* const constSixTwentyFourEpsSig12_2D
357  = sixTwentyFourEpsilonSigma12_2D_;
358  double const* const* const constShifts2D = shifts2D_;
359  for (ii = 0; ii < cachedNumberOfParticles_; ++ii)
360  {
361  if (particleContributing[ii])
362  {
363  modelComputeArguments->GetNeighborList(0, ii, &numnei, &n1atom);
364  int const numNei = numnei;
365  int const * const n1Atom = n1atom;
366  int const i = ii;
367  int const iSpecies = particleSpeciesCodes[i];
368 
369  // Setup loop over neighbors of current particle
370  for (int jj = 0; jj < numNei; ++jj)
371  {
372  int const j = n1Atom[jj];
373  int const jSpecies = particleSpeciesCodes[j];
374  double* r_ij;
375  double r_ijValue[DIMENSION];
376  // Compute r_ij
377  r_ij = r_ijValue;
378  for (int k = 0; k < DIMENSION; ++k)
379  r_ij[k] = coordinates[j][k] - coordinates[i][k];
380  double const* const r_ij_const = const_cast<double*>(r_ij);
381 
382  // compute distance squared
383  double const rij2 =
384  r_ij_const[0] * r_ij_const[0] +
385  r_ij_const[1] * r_ij_const[1] +
386  r_ij_const[2] * r_ij_const[2];
387 
388  if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
389  { // compute contribution to energy, force, etc.
390  double phi = 0.0;
391  double dphiByR = 0.0;
392  double d2phi = 0.0;
393  double dEidrByR = 0.0;
394  double d2Eidr2 = 0.0;
395  double const r2iv = 1.0/rij2;
396  double const r6iv = r2iv*r2iv*r2iv;
397  // Compute pair potential and its derivatives
398  if (isComputeProcess_d2Edr2 == true)
399  { // Compute d2phi
400  d2phi =
401  r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]*r6iv -
402  constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
403  * r2iv;
404  d2Eidr2 = 0.5*d2phi;
405  }
406 
407  if ((isComputeProcess_dEdr == true) || (isComputeForces == true) ||
408  (isComputeVirial == true) || (isComputeParticleVirial == true))
409  { // Compute dphi
410  dphiByR =
411  r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
412  constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
413  * r2iv;
414  dEidrByR = 0.5*dphiByR;
415  }
416 
417  if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
418  { // Compute phi
419  if (isShift == true)
420  {
421  LENNARD_JONES_PHI(- constShifts2D[iSpecies][jSpecies]);
422  }
423  else
424  {
426  }
427  }
428 
429  // Contribution to energy
430  if (isComputeEnergy == true)
431  {
432  *energy += 0.5*phi;
433  }
434 
435  // Contribution to particleEnergy
436  if (isComputeParticleEnergy == true)
437  {
438  double const halfPhi = 0.5*phi;
439  particleEnergy[i] += halfPhi;
440  }
441 
442  // Contribution to forces
443  if (isComputeForces == true)
444  {
445  for (int k = 0; k < DIMENSION; ++k)
446  {
447  double const contrib = dEidrByR * r_ij_const[k];
448  forces[i][k] += contrib;
449  forces[j][k] -= contrib;
450  }
451  }
452 
453  // Call process_dEdr
454  if ((isComputeProcess_dEdr == true) ||
455  (isComputeVirial == true) ||
456  (isComputeParticleVirial == true))
457  {
458  double const rij = sqrt(rij2);
459  double const dEidr = dEidrByR*rij;
460 
461  if (isComputeProcess_dEdr == true)
462  {
463  ier = modelComputeArguments
464  ->ProcessDEDrTerm(dEidr, rij, r_ij_const, i, j);
465  if (ier)
466  {
467  LOG_ERROR("process_dEdr");
468  return ier;
469  }
470  }
471 
472  if (isComputeVirial == true)
473  {
474  ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
475  }
476 
477  if (isComputeParticleVirial == true)
478  {
479  ProcessParticleVirialTerm(dEidr, rij, r_ij_const, i, j,
480  particleVirial);
481  }
482  }
483 
484  // Call process_d2Edr2
485  if (isComputeProcess_d2Edr2 == true)
486  {
487  double const rij = sqrt(rij2);
488  double const R_pairs[2] = {rij, rij};
489  double const* const pRs = &R_pairs[0];
490  double const Rij_pairs[6]
491  = {r_ij_const[0], r_ij_const[1], r_ij_const[2],
492  r_ij_const[0], r_ij_const[1], r_ij_const[2]};
493  double const* const pRijConsts = &Rij_pairs[0];
494  int const i_pairs[2] = {i, i};
495  int const j_pairs[2] = {j, j};
496  int const* const pis = &i_pairs[0];
497  int const* const pjs = &j_pairs[0];
498 
499  ier = modelComputeArguments
500  ->ProcessD2EDr2Term(d2Eidr2, pRs, pRijConsts, pis, pjs);
501  if (ier)
502  {
503  LOG_ERROR("process_d2Edr2");
504  return ier;
505  }
506  }
507  } // if particleContributing
508  } // if particles i and j interact
509  } // end of first neighbor loop
510  } // end of loop over contributing particles
511 
512  // everything is good
513  ier = false;
514  return ier;
515 }
516 
517 #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)