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 "KIM_LogVerbosity.hpp"
36 #include "LennardJones612.hpp"
37 
38 #define DIMENSION 3
39 #define ONE 1.0
40 #define HALF 0.5
41 
42 #define MAX_PARAMETER_FILES 1
43 
44 #define PARAM_SHIFT_INDEX 0
45 #define PARAM_CUTOFFS_INDEX 1
46 #define PARAM_EPSILONS_INDEX 2
47 #define PARAM_SIGMAS_INDEX 3
48 
49 
50 //==============================================================================
51 //
52 // Type definitions, enumerations, and helper function prototypes
53 //
54 //==============================================================================
55 
56 // type declaration for get neighbor functions
57 typedef int (GetNeighborFunction)(void const * const, int const,
58  int * const, int const ** const);
59 // type declaration for vector of constant dimension
60 typedef double VectorOfSizeDIM[DIMENSION];
61 
62 // helper routine declarations
63 void AllocateAndInitialize2DArray(double**& arrayPtr, int const extentZero,
64  int const extentOne);
65 void Deallocate2DArray(double**& arrayPtr);
66 
67 //==============================================================================
68 //
69 // Declaration of LennardJones612Implementation class
70 //
71 //==============================================================================
72 
73 //******************************************************************************
75 {
76  public:
78  KIM::ModelDriverCreate * const modelDriverCreate,
79  KIM::LengthUnit const requestedLengthUnit,
80  KIM::EnergyUnit const requestedEnergyUnit,
81  KIM::ChargeUnit const requestedChargeUnit,
82  KIM::TemperatureUnit const requestedTemperatureUnit,
83  KIM::TimeUnit const requestedTimeUnit,
84  int * const ier);
85  ~LennardJones612Implementation(); // no explicit Destroy() needed here
86 
87  int Refresh(KIM::ModelRefresh * const modelRefresh);
88  int Compute(KIM::ModelCompute const * const modelCompute);
89 
90  private:
91  // Constant values that never change
92  // Set in constructor (via SetConstantValues)
93  //
94  //
95  // LennardJones612Implementation: constants
96  int numberModelSpecies_;
97  std::vector<int> modelSpeciesCodeList_;
98  int numberUniqueSpeciesPairs_;
99 
100 
101  // Constant values that are read from the input files and never change
102  // Set in constructor (via functions listed below)
103  //
104  //
105  // KIM API: Model Fixed Parameters
106  // Memory allocated in AllocateFixedParameterMemory()
107  // Memory deallocated in destructor
108  // Data set in ReadParameterFile routines
109  // none
110  //
111  // KIM API: Model Free Parameters whose (pointer) values never change
112  // Memory allocated in AllocateFreeParameterMemory() (from constructor)
113  // Memory deallocated in destructor
114  // Data set in ReadParameterFile routines OR by KIM Simulator
115  int shift_;
116  double* cutoffs_;
117  double* epsilons_;
118  double* sigmas_;
119 
120  // Mutable values that only change when reinit() executes
121  // Set in Reinit (via SetReinitMutableValues)
122  //
123  //
124  // KIM API: Model Fixed Parameters
125  // none
126  //
127  // KIM API: Model Free Parameters
128  // none
129  //
130  // LennardJones612Implementation: values
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_;
139  double** shifts2D_;
140 
141 
142  // Mutable values that can change with each call to Reinit() and Compute()
143  // Memory may be reallocated on each call
144  //
145  //
146  // LennardJones612Implementation: values that change
147  int cachedNumberOfParticles_;
148 
149 
150  // Helper methods
151  //
152  //
153  // Related to constructor
154  void AllocateFreeParameterMemory();
155  static int OpenParameterFiles(
156  KIM::ModelDriverCreate * const modelDriverCreate,
157  int const numberParameterFiles,
158  FILE* parameterFilePointers[MAX_PARAMETER_FILES]);
159  static void CloseParameterFiles(
160  int const numberParameterFiles,
161  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
162  int ProcessParameterFiles(
163  KIM::ModelDriverCreate * const modelDriverCreate,
164  int const numberParameterFiles,
165  FILE* const parameterFilePointers[MAX_PARAMETER_FILES]);
166  void getNextDataLine(FILE* const filePtr, char* const nextLine,
167  int const maxSize, int* endOfFileFlag);
168  int ConvertUnits(
169  KIM::ModelDriverCreate * const modelDriverCreate,
170  KIM::LengthUnit const requestedLengthUnit,
171  KIM::EnergyUnit const requestedEnergyUnit,
172  KIM::ChargeUnit const requestedChargeUnit,
173  KIM::TemperatureUnit const requestedTemperatureUnit,
174  KIM::TimeUnit const requestedTimeUnit);
175  int RegisterKIMModelSettings(
176  KIM::ModelDriverCreate * const modelDriverCreate);
177  int RegisterKIMParameters(
178  KIM::ModelDriverCreate * const modelDriverCreate);
179  int RegisterKIMFunctions(
180  KIM::ModelDriverCreate * const modelDriverCreate) const;
181  //
182  // Related to Reinit()
183  template<class ModelObj>
184  int SetReinitMutableValues(ModelObj * const modelObj);
185 
186  //
187  // Related to Compute()
188  int SetComputeMutableValues(KIM::ModelCompute const * const modelCompute,
189  bool& isComputeProcess_dEdr,
190  bool& isComputeProcess_d2Edr2,
191  bool& isComputeEnergy,
192  bool& isComputeForces,
193  bool& isComputeParticleEnergy,
194  int const*& particleSpeciesCodes,
195  int const*& particleContributing,
197  double*& energy,
198  double*& particleEnergy,
199  VectorOfSizeDIM*& forces);
200  int CheckParticleSpeciesCodes(KIM::ModelCompute const * const modelCompute,
201  int const* const particleSpeciesCodes) const;
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;
208 
209  // compute functions
210  template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
211  bool isComputeEnergy, bool isComputeForces,
212  bool isComputeParticleEnergy, bool isShift >
213  int Compute(KIM::ModelCompute const * const modelCompute,
214  const int* const particleSpeciesCodes,
215  const int* const particleContributing,
216  const VectorOfSizeDIM* const coordinates,
217  double* const energy,
218  VectorOfSizeDIM* const forces,
219  double* const particleEnergy);
220 };
221 
222 //==============================================================================
223 //
224 // Definition of LennardJones612Implementation::Compute functions
225 //
226 // NOTE: Here we rely on the compiler optimizations to prune dead code
227 // after the template expansions. This provides high efficiency
228 // and easy maintenance.
229 //
230 //==============================================================================
231 
232 //******************************************************************************
233 // MACRO to compute Lennard-Jones phi
234 // (used for efficiency)
235 //
236 // exshift - expression to be added to the end of the phi value
237 #define LENNARD_JONES_PHI(exshift) \
238  phi = r6iv * (constFourEpsSig12_2D[iSpecies][jSpecies]*r6iv - \
239  constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
240 
241 //******************************************************************************
243 template< bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
244  bool isComputeEnergy, bool isComputeForces,
245  bool isComputeParticleEnergy, bool isShift >
247  KIM::ModelCompute const * const modelCompute,
248  const int* const particleSpeciesCodes,
249  const int* const particleContributing,
250  const VectorOfSizeDIM* const coordinates,
251  double* const energy,
252  VectorOfSizeDIM* const forces,
253  double* const particleEnergy)
254 {
255  int ier = false;
256 
257  if ((isComputeEnergy == false) &&
258  (isComputeParticleEnergy == false) &&
259  (isComputeForces == false) &&
260  (isComputeProcess_dEdr == false) &&
261  (isComputeProcess_d2Edr2 == false))
262  return ier;
263 
264  // initialize energy and forces
265  if (isComputeEnergy == true)
266  {
267  *energy = 0.0;
268  }
269  if (isComputeParticleEnergy == true)
270  {
271  int const cachedNumParticles = cachedNumberOfParticles_;
272  for (int i = 0; i < cachedNumParticles; ++i)
273  {
274  particleEnergy[i] = 0.0;
275  }
276  }
277  if (isComputeForces == true)
278  {
279  int const cachedNumParticles = cachedNumberOfParticles_;
280  for (int i = 0; i < cachedNumParticles; ++i)
281  {
282  for (int j = 0; j < DIMENSION; ++j)
283  forces[i][j] = 0.0;
284  }
285  }
286 
287  // calculate contribution from pair function
288  //
289  // Setup loop over contributing particles
290  int ii = 0;
291  int numnei = 0;
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)
306  {
307  if (particleContributing[ii])
308  {
309  modelCompute->GetNeighborList(0, ii, &numnei, &n1atom);
310  int const numNei = numnei;
311  int const * const n1Atom = n1atom;
312  int const i = ii;
313  int const iSpecies = particleSpeciesCodes[i];
314 
315  // Setup loop over neighbors of current particle
316  for (int jj = 0; jj < numNei; ++jj)
317  {
318  int const j = n1Atom[jj];
319  int const jSpecies = particleSpeciesCodes[j];
320  double* r_ij;
321  double r_ijValue[DIMENSION];
322  // Compute r_ij
323  r_ij = r_ijValue;
324  for (int k = 0; k < DIMENSION; ++k)
325  r_ij[k] = coordinates[j][k] - coordinates[i][k];
326  double const* const r_ij_const = const_cast<double*>(r_ij);
327 
328  // compute distance squared
329  double const rij2 =
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];
333 
334  if (rij2 <= constCutoffsSq2D[iSpecies][jSpecies])
335  { // compute contribution to energy, force, etc.
336  double phi = 0.0;
337  double dphiByR = 0.0;
338  double d2phi = 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;
343  // Compute pair potential and its derivatives
344  if (isComputeProcess_d2Edr2 == true)
345  { // Compute d2phi
346  d2phi =
347  r6iv * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies]*r6iv -
348  constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
349  * r2iv;
350  d2Eidr2 = 0.5*d2phi;
351  }
352 
353  if ((isComputeProcess_dEdr == true) || (isComputeForces == true))
354  { // Compute dphi
355  dphiByR =
356  r6iv * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies] -
357  constFortyEightEpsSig12_2D[iSpecies][jSpecies]*r6iv)
358  * r2iv;
359  dEidrByR = 0.5*dphiByR;
360  }
361 
362  if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
363  { // Compute phi
364  if (isShift == true)
365  {
366  LENNARD_JONES_PHI(- constShifts2D[iSpecies][jSpecies]);
367  }
368  else
369  {
371  }
372  }
373 
374  // Contribution to energy
375  if (isComputeEnergy == true)
376  {
377  *energy += 0.5*phi;
378  }
379 
380  // Contribution to particleEnergy
381  if (isComputeParticleEnergy == true)
382  {
383  double const halfPhi = 0.5*phi;
384  particleEnergy[i] += halfPhi;
385  }
386 
387  // Contribution to forces
388  if (isComputeForces == true)
389  {
390  for (int k = 0; k < DIMENSION; ++k)
391  {
392  double const contrib = dEidrByR * r_ij_const[k];
393  forces[i][k] += contrib;
394  forces[j][k] -= contrib;
395  }
396  }
397 
398  // Call process_dEdr
399  if (isComputeProcess_dEdr == true)
400  {
401  double const rij = sqrt(rij2);
402  double const dEidr = dEidrByR*rij;
403  ier = modelCompute->ProcessDEDrTerm(dEidr, rij, r_ij_const, i, j);
404  if (ier)
405  {
406  LOG_ERROR("process_dEdr");
407  return ier;
408  }
409  }
410 
411  // Call process_d2Edr2
412  if (isComputeProcess_d2Edr2 == true)
413  {
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];
425 
426  ier = modelCompute->ProcessD2EDr2Term(d2Eidr2, pRs, pRijConsts, pis,
427  pjs);
428  if (ier)
429  {
430  LOG_ERROR("process_d2Edr2");
431  return ier;
432  }
433  }
434  } // if particleContributing
435  } // if particles i and j interact
436  } // end of first neighbor loop
437  } // end of loop over contributing particles
438 
439  // everything is good
440  ier = false;
441  return ier;
442 }
443 
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)
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)