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 "KIM_ModelDriverHeaders.hpp"
36 #include <cmath>
37 #include <cstdio>
38 #include <vector>
39
40 #define DIMENSION 3
41 #define ONE 1.0
42 #define HALF 0.5
43
44 #define MAX_PARAMETER_FILES 1
45
46 #define PARAM_SHIFT_INDEX 0
47 #define PARAM_CUTOFFS_INDEX 1
48 #define PARAM_EPSILONS_INDEX 2
49 #define PARAM_SIGMAS_INDEX 3
50
51
52 //==============================================================================
53 //
54 // Type definitions, enumerations, and helper function prototypes
55 //
56 //==============================================================================
57
58 // type declaration for get neighbor functions
59 typedef int(GetNeighborFunction)(void const * const,
60 int const,
61 int * const,
62 int const ** const);
63 // type declaration for vector of constant dimension
64 typedef double VectorOfSizeDIM[DIMENSION];
65 typedef double VectorOfSizeSix[6];
66
67 // helper routine declarations
68 void AllocateAndInitialize2DArray(double **& arrayPtr,
69 int const extentZero,
70 int const extentOne);
71 void Deallocate2DArray(double **& arrayPtr);
72
73 //==============================================================================
74 //
75 // Declaration of LennardJones612Implementation class
76 //
77 //==============================================================================
78
79 //******************************************************************************
80 class LennardJones612Implementation
81 {
82 public:
83 LennardJones612Implementation(
84 KIM::ModelDriverCreate * const modelDriverCreate,
85 KIM::LengthUnit const requestedLengthUnit,
86 KIM::EnergyUnit const requestedEnergyUnit,
87 KIM::ChargeUnit const requestedChargeUnit,
88 KIM::TemperatureUnit const requestedTemperatureUnit,
89 KIM::TimeUnit const requestedTimeUnit,
90 int * const ier);
91 ~LennardJones612Implementation(); // no explicit Destroy() needed here
92
93 int Refresh(KIM::ModelRefresh * const modelRefresh);
94 int Compute(KIM::ModelCompute const * const modelCompute,
95 KIM::ModelComputeArguments const * const modelComputeArguments);
96 int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate * const
97 modelComputeArgumentsCreate) const;
98 int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy * const
99 modelComputeArgumentsDestroy) const;
100
101
102 private:
103 // Constant values that never change
104 // Set in constructor (via SetConstantValues)
105 //
106 //
107 // LennardJones612Implementation: constants
108 int numberModelSpecies_;
109 std::vector<int> modelSpeciesCodeList_;
110 int numberUniqueSpeciesPairs_;
111
112
113 // Constant values that are read from the input files and never change
114 // Set in constructor (via functions listed below)
115 //
116 //
117 // Private Model Parameters
118 // Memory allocated in AllocatePrivateParameterMemory() (from constructor)
119 // Memory deallocated in destructor
120 // Data set in ReadParameterFile routines
121 // none
122 //
123 // KIM API: Model Parameters whose (pointer) values never change
124 // Memory allocated in AllocateParameterMemory() (from constructor)
125 // Memory deallocated in destructor
126 // Data set in ReadParameterFile routines OR by KIM Simulator
127 int shift_;
128 double * cutoffs_;
129 double * epsilons_;
130 double * sigmas_;
131
132 // Mutable values that only change when Refresh() executes
133 // Set in Refresh (via SetRefreshMutableValues)
134 //
135 //
136 // KIM API: Model Parameters (can be changed directly by KIM Simulator)
137 // none
138 //
139 // LennardJones612Implementation: values (changed only by Refresh())
140 double influenceDistance_;
141 double ** cutoffsSq2D_;
142 int modelWillNotRequestNeighborsOfNoncontributingParticles_;
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
167 OpenParameterFiles(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,
175 char * const nextLine,
176 int const maxSize,
177 int * endOfFileFlag);
178 static void
179 CloseParameterFiles(int const numberParameterFiles,
180 FILE * const parameterFilePointers[MAX_PARAMETER_FILES]);
181 int ConvertUnits(KIM::ModelDriverCreate * const modelDriverCreate,
182 KIM::LengthUnit const requestedLengthUnit,
183 KIM::EnergyUnit const requestedEnergyUnit,
184 KIM::ChargeUnit const requestedChargeUnit,
185 KIM::TemperatureUnit const requestedTemperatureUnit,
186 KIM::TimeUnit const requestedTimeUnit);
187 int RegisterKIMModelSettings(
188 KIM::ModelDriverCreate * const modelDriverCreate) const;
189 int RegisterKIMComputeArgumentsSettings(
190 KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
191 const;
192 int RegisterKIMParameters(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,
213 VectorOfSizeDIM const *& coordinates,
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,
244 bool isComputeProcess_d2Edr2,
245 bool isComputeEnergy,
246 bool isComputeForces,
247 bool isComputeParticleEnergy,
248 bool isComputeVirial,
249 bool isComputeParticleVirial,
250 bool isShift>
251 int Compute(KIM::ModelCompute const * const modelCompute,
252 KIM::ModelComputeArguments const * const modelComputeArguments,
253 const int * const particleSpeciesCodes,
254 const int * const particleContributing,
255 const VectorOfSizeDIM * const coordinates,
256 double * const energy,
257 VectorOfSizeDIM * const forces,
258 double * const particleEnergy,
259 VectorOfSizeSix virial,
260 VectorOfSizeSix * const particleVirial) const;
261 };
262
263 //==============================================================================
264 //
265 // Definition of LennardJones612Implementation::Compute functions
266 //
267 // NOTE: Here we rely on the compiler optimizations to prune dead code
268 // after the template expansions. This provides high efficiency
269 // and easy maintenance.
270 //
271 //==============================================================================
272
273 //******************************************************************************
274 // MACRO to compute Lennard-Jones phi
275 // (used for efficiency)
276 //
277 // exshift - expression to be added to the end of the phi value
278 #define LENNARD_JONES_PHI(exshift) \
279 phi = r6iv \
280 * (constFourEpsSig12_2D[iSpecies][jSpecies] * r6iv \
281 - constFourEpsSig6_2D[iSpecies][jSpecies]) exshift;
282
283 //******************************************************************************
284 #undef KIM_LOGGER_OBJECT_NAME
285 #define KIM_LOGGER_OBJECT_NAME modelCompute
286 //
287 template<bool isComputeProcess_dEdr,
288 bool isComputeProcess_d2Edr2,
289 bool isComputeEnergy,
290 bool isComputeForces,
291 bool isComputeParticleEnergy,
292 bool isComputeVirial,
293 bool isComputeParticleVirial,
294 bool isShift>
Compute(KIM::ModelCompute const * const modelCompute,KIM::ModelComputeArguments const * const modelComputeArguments,const int * const particleSpeciesCodes,const int * const particleContributing,const VectorOfSizeDIM * const coordinates,double * const energy,VectorOfSizeDIM * const forces,double * const particleEnergy,VectorOfSizeSix virial,VectorOfSizeSix * const particleVirial) const295 int LennardJones612Implementation::Compute(
296 KIM::ModelCompute const * const modelCompute,
297 KIM::ModelComputeArguments const * const modelComputeArguments,
298 const int * const particleSpeciesCodes,
299 const int * const particleContributing,
300 const VectorOfSizeDIM * const coordinates,
301 double * const energy,
302 VectorOfSizeDIM * const forces,
303 double * const particleEnergy,
304 VectorOfSizeSix virial,
305 VectorOfSizeSix * const particleVirial) const
306 {
307 int ier = false;
308
309 if ((isComputeEnergy == false) && (isComputeParticleEnergy == false)
310 && (isComputeForces == false) && (isComputeProcess_dEdr == false)
311 && (isComputeProcess_d2Edr2 == false) && (isComputeVirial == false)
312 && (isComputeParticleVirial == false))
313 return ier;
314
315 // initialize energy and forces
316 if (isComputeEnergy == true) { *energy = 0.0; }
317 if (isComputeVirial == true)
318 {
319 for (int i = 0; i < 6; ++i) virial[i] = 0.0;
320 }
321 if (isComputeParticleEnergy == true)
322 {
323 int const cachedNumParticles = cachedNumberOfParticles_;
324 for (int i = 0; i < cachedNumParticles; ++i) { particleEnergy[i] = 0.0; }
325 }
326 if (isComputeForces == true)
327 {
328 int const cachedNumParticles = cachedNumberOfParticles_;
329 for (int i = 0; i < cachedNumParticles; ++i)
330 {
331 for (int j = 0; j < DIMENSION; ++j) forces[i][j] = 0.0;
332 }
333 }
334 if (isComputeParticleVirial == true)
335 {
336 int const cachedNumParticles = cachedNumberOfParticles_;
337 for (int i = 0; i < cachedNumParticles; ++i)
338 {
339 for (int j = 0; j < 6; ++j) 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 int const jContrib = particleContributing[j];
376
377 if (!(jContrib && (j < i))) // effective half-list
378 {
379 int const jSpecies = particleSpeciesCodes[j];
380 double * r_ij;
381 double r_ijValue[DIMENSION];
382 // Compute r_ij
383 r_ij = r_ijValue;
384 for (int k = 0; k < DIMENSION; ++k)
385 r_ij[k] = coordinates[j][k] - coordinates[i][k];
386 double const * const r_ij_const = const_cast<double *>(r_ij);
387
388 // compute distance squared
389 double const rij2 = 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
407 * (constSixTwentyFourEpsSig12_2D[iSpecies][jSpecies] * r6iv
408 - constOneSixtyEightEpsSig6_2D[iSpecies][jSpecies])
409 * r2iv;
410 if (jContrib == 1) { d2Eidr2 = d2phi; }
411 else
412 {
413 d2Eidr2 = 0.5 * d2phi;
414 }
415 }
416
417 if ((isComputeProcess_dEdr == true) || (isComputeForces == true)
418 || (isComputeVirial == true)
419 || (isComputeParticleVirial == true))
420 { // Compute dphi
421 dphiByR
422 = r6iv
423 * (constTwentyFourEpsSig6_2D[iSpecies][jSpecies]
424 - constFortyEightEpsSig12_2D[iSpecies][jSpecies] * r6iv)
425 * r2iv;
426 if (jContrib == 1) { dEidrByR = dphiByR; }
427 else
428 {
429 dEidrByR = 0.5 * dphiByR;
430 }
431 }
432
433 if ((isComputeEnergy == true) || (isComputeParticleEnergy == true))
434 { // Compute phi
435 if (isShift == true)
436 {
437 LENNARD_JONES_PHI(-constShifts2D[iSpecies][jSpecies]);
438 }
439 else
440 {
441 LENNARD_JONES_PHI(;);
442 }
443 }
444
445 // Contribution to energy
446 if (isComputeEnergy == true)
447 {
448 if (jContrib == 1) { *energy += phi; }
449 else
450 {
451 *energy += 0.5 * phi;
452 }
453 }
454
455 // Contribution to particleEnergy
456 if (isComputeParticleEnergy == true)
457 {
458 double const halfPhi = 0.5 * phi;
459 particleEnergy[i] += halfPhi;
460 if (jContrib == 1) { particleEnergy[j] += halfPhi; }
461 }
462
463 // Contribution to forces
464 if (isComputeForces == true)
465 {
466 for (int k = 0; k < DIMENSION; ++k)
467 {
468 double const contrib = dEidrByR * r_ij_const[k];
469 forces[i][k] += contrib;
470 forces[j][k] -= contrib;
471 }
472 }
473
474 // Call process_dEdr
475 if ((isComputeProcess_dEdr == true) || (isComputeVirial == true)
476 || (isComputeParticleVirial == true))
477 {
478 double const rij = sqrt(rij2);
479 double const dEidr = dEidrByR * rij;
480
481 if (isComputeProcess_dEdr == true)
482 {
483 ier = modelComputeArguments->ProcessDEDrTerm(
484 dEidr, rij, r_ij_const, i, j);
485 if (ier)
486 {
487 LOG_ERROR("process_dEdr");
488 return ier;
489 }
490 }
491
492 if (isComputeVirial == true)
493 {
494 ProcessVirialTerm(dEidr, rij, r_ij_const, i, j, virial);
495 }
496
497 if (isComputeParticleVirial == true)
498 {
499 ProcessParticleVirialTerm(
500 dEidr, rij, r_ij_const, i, j, particleVirial);
501 }
502 }
503
504 // Call process_d2Edr2
505 if (isComputeProcess_d2Edr2 == true)
506 {
507 double const rij = sqrt(rij2);
508 double const R_pairs[2] = {rij, rij};
509 double const * const pRs = &R_pairs[0];
510 double const Rij_pairs[6] = {r_ij_const[0],
511 r_ij_const[1],
512 r_ij_const[2],
513 r_ij_const[0],
514 r_ij_const[1],
515 r_ij_const[2]};
516 double const * const pRijConsts = &Rij_pairs[0];
517 int const i_pairs[2] = {i, i};
518 int const j_pairs[2] = {j, j};
519 int const * const pis = &i_pairs[0];
520 int const * const pjs = &j_pairs[0];
521
522 ier = modelComputeArguments->ProcessD2EDr2Term(
523 d2Eidr2, pRs, pRijConsts, pis, pjs);
524 if (ier)
525 {
526 LOG_ERROR("process_d2Edr2");
527 return ier;
528 }
529 }
530 } // if particleContributing
531 } // if i < j or j non contributing
532 } // if particles i and j interact
533 } // end of first neighbor loop
534 } // end of loop over contributing particles
535
536 // everything is good
537 ier = false;
538 return ier;
539 }
540
541 #endif // LENNARD_JONES_612_IMPLEMENTATION_HPP_
542