kim-api 2.3.0+AppleClang.AppleClang.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
Loading...
Searching...
No Matches
LennardJones612Implementation.cpp
Go to the documentation of this file.
1//
2// KIM-API: An API for interatomic models
3// Copyright (c) 2013--2022, Regents of the University of Minnesota.
4// All rights reserved.
5//
6// Contributors:
7// Ryan S. Elliott
8// Stephen M. Whalen
9// Andrew Akerson
10//
11// SPDX-License-Identifier: LGPL-2.1-or-later
12//
13// This library is free software; you can redistribute it and/or
14// modify it under the terms of the GNU Lesser General Public
15// License as published by the Free Software Foundation; either
16// version 2.1 of the License, or (at your option) any later version.
17//
18// This library is distributed in the hope that it will be useful,
19// but WITHOUT ANY WARRANTY; without even the implied warranty of
20// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21// Lesser General Public License for more details.
22//
23// You should have received a copy of the GNU Lesser General Public License
24// along with this library; if not, write to the Free Software Foundation,
25// Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26//
27
28
29#include <cmath>
30#include <cstring>
31#include <iostream> // IWYU pragma: keep // BUG WORK-AROUND
32#include <map>
33#include <sstream>
34#include <string>
35
37#include "LennardJones612.hpp"
39
40#define MAXLINE 1024
41#define IGNORE_RESULT(fn) \
42 if (fn) {}
43
44
45//==============================================================================
46//
47// Implementation of LennardJones612Implementation public member functions
48//
49//==============================================================================
50
51//******************************************************************************
52#undef KIM_LOGGER_OBJECT_NAME
53#define KIM_LOGGER_OBJECT_NAME modelDriverCreate
54//
56 KIM::ModelDriverCreate * const modelDriverCreate,
57 KIM::LengthUnit const requestedLengthUnit,
58 KIM::EnergyUnit const requestedEnergyUnit,
59 KIM::ChargeUnit const requestedChargeUnit,
60 KIM::TemperatureUnit const requestedTemperatureUnit,
61 KIM::TimeUnit const requestedTimeUnit,
62 int * const ier) :
63 numberModelSpecies_(0),
64 numberUniqueSpeciesPairs_(0),
65 shift_(0),
66 cutoffs_(NULL),
67 epsilons_(NULL),
68 sigmas_(NULL),
69 influenceDistance_(0.0),
70 cutoffsSq2D_(NULL),
71 modelWillNotRequestNeighborsOfNoncontributingParticles_(1),
72 fourEpsilonSigma6_2D_(NULL),
73 fourEpsilonSigma12_2D_(NULL),
74 twentyFourEpsilonSigma6_2D_(NULL),
75 fortyEightEpsilonSigma12_2D_(NULL),
76 oneSixtyEightEpsilonSigma6_2D_(NULL),
77 sixTwentyFourEpsilonSigma12_2D_(NULL),
78 shifts2D_(NULL),
79 cachedNumberOfParticles_(0)
80{
81 FILE * parameterFilePointers[MAX_PARAMETER_FILES];
82 int numberParameterFiles;
83 modelDriverCreate->GetNumberOfParameterFiles(&numberParameterFiles);
84 *ier = OpenParameterFiles(
85 modelDriverCreate, numberParameterFiles, parameterFilePointers);
86 if (*ier) return;
87
88 *ier = ProcessParameterFiles(
89 modelDriverCreate, numberParameterFiles, parameterFilePointers);
90 CloseParameterFiles(numberParameterFiles, parameterFilePointers);
91 if (*ier) return;
92
93 *ier = ConvertUnits(modelDriverCreate,
94 requestedLengthUnit,
95 requestedEnergyUnit,
96 requestedChargeUnit,
97 requestedTemperatureUnit,
98 requestedTimeUnit);
99 if (*ier) return;
100
101 *ier = SetRefreshMutableValues(modelDriverCreate);
102 if (*ier) return;
103
104 *ier = RegisterKIMModelSettings(modelDriverCreate);
105 if (*ier) return;
106
107 *ier = RegisterKIMParameters(modelDriverCreate);
108 if (*ier) return;
109
110 *ier = RegisterKIMFunctions(modelDriverCreate);
111 if (*ier) return;
112
113 // everything is good
114 *ier = false;
115 return;
116}
117
118//******************************************************************************
120{ // note: it is ok to delete a null pointer and we have ensured that
121 // everything is initialized to null
122
123 delete[] cutoffs_;
124 Deallocate2DArray(cutoffsSq2D_);
125 delete[] epsilons_;
126 delete[] sigmas_;
127 Deallocate2DArray(fourEpsilonSigma6_2D_);
128 Deallocate2DArray(fourEpsilonSigma12_2D_);
129 Deallocate2DArray(twentyFourEpsilonSigma6_2D_);
130 Deallocate2DArray(fortyEightEpsilonSigma12_2D_);
131 Deallocate2DArray(oneSixtyEightEpsilonSigma6_2D_);
132 Deallocate2DArray(sixTwentyFourEpsilonSigma12_2D_);
133 Deallocate2DArray(shifts2D_);
134}
135
136//******************************************************************************
137#undef KIM_LOGGER_OBJECT_NAME
138#define KIM_LOGGER_OBJECT_NAME modelRefresh
139//
141 KIM::ModelRefresh * const modelRefresh)
142{
143 int ier;
144
145 ier = SetRefreshMutableValues(modelRefresh);
146 if (ier) return ier;
147
148 // nothing else to do for this case
149
150 // everything is good
151 ier = false;
152 return ier;
153}
154
155//******************************************************************************
157 KIM::ModelCompute const * const modelCompute,
158 KIM::ModelComputeArguments const * const modelComputeArguments)
159{
160 int ier;
161
162 // KIM API Model Input compute flags
163 bool isComputeProcess_dEdr = false;
164 bool isComputeProcess_d2Edr2 = false;
165 //
166 // KIM API Model Output compute flags
167 bool isComputeEnergy = false;
168 bool isComputeForces = false;
169 bool isComputeParticleEnergy = false;
170 bool isComputeVirial = false;
171 bool isComputeParticleVirial = false;
172 //
173 // KIM API Model Input
174 int const * particleSpeciesCodes = NULL;
175 int const * particleContributing = NULL;
176 VectorOfSizeDIM const * coordinates = NULL;
177 //
178 // KIM API Model Output
179 double * energy = NULL;
180 double * particleEnergy = NULL;
181 VectorOfSizeDIM * forces = NULL;
182 VectorOfSizeSix * virial = NULL;
183 VectorOfSizeSix * particleVirial = NULL;
184 ier = SetComputeMutableValues(modelComputeArguments,
185 isComputeProcess_dEdr,
186 isComputeProcess_d2Edr2,
187 isComputeEnergy,
188 isComputeForces,
189 isComputeParticleEnergy,
190 isComputeVirial,
191 isComputeParticleVirial,
192 particleSpeciesCodes,
193 particleContributing,
194 coordinates,
195 energy,
196 particleEnergy,
197 forces,
198 virial,
199 particleVirial);
200 if (ier) return ier;
201
202 // Skip this check for efficiency
203 //
204 // ier = CheckParticleSpecies(modelComputeArguments, particleSpeciesCodes);
205 // if (ier) return ier;
206
207 bool const isShift = (1 == shift_);
208
210 return ier;
211}
212
213//******************************************************************************
215 KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
216{
217 int ier;
218
219 ier = RegisterKIMComputeArgumentsSettings(modelComputeArgumentsCreate);
220 if (ier) return ier;
221
222 // nothing else to do for this case
223
224 // everything is good
225 ier = false;
226 return ier;
227}
228
229//******************************************************************************
232 /* modelComputeArgumentsDestroy */) const
233{
234 int ier;
235
236 // nothing else to do for this case
237
238 // everything is good
239 ier = false;
240 return ier;
241}
242
243//==============================================================================
244//
245// Implementation of LennardJones612Implementation private member functions
246//
247//==============================================================================
248
249//******************************************************************************
250void LennardJones612Implementation::AllocatePrivateParameterMemory()
251{
252 // nothing to do for this case
253}
254
255//******************************************************************************
256void LennardJones612Implementation::AllocateParameterMemory()
257{ // allocate memory for data
258 cutoffs_ = new double[numberUniqueSpeciesPairs_];
260 cutoffsSq2D_, numberModelSpecies_, numberModelSpecies_);
261
262 epsilons_ = new double[numberUniqueSpeciesPairs_];
263 sigmas_ = new double[numberUniqueSpeciesPairs_];
265 fourEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
267 fourEpsilonSigma12_2D_, numberModelSpecies_, numberModelSpecies_);
269 twentyFourEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
271 fortyEightEpsilonSigma12_2D_, numberModelSpecies_, numberModelSpecies_);
273 oneSixtyEightEpsilonSigma6_2D_, numberModelSpecies_, numberModelSpecies_);
274 AllocateAndInitialize2DArray(sixTwentyFourEpsilonSigma12_2D_,
275 numberModelSpecies_,
276 numberModelSpecies_);
277
279 shifts2D_, numberModelSpecies_, numberModelSpecies_);
280}
281
282//******************************************************************************
283#undef KIM_LOGGER_OBJECT_NAME
284#define KIM_LOGGER_OBJECT_NAME modelDriverCreate
285//
286int LennardJones612Implementation::OpenParameterFiles(
287 KIM::ModelDriverCreate * const modelDriverCreate,
288 int const numberParameterFiles,
289 FILE * parameterFilePointers[MAX_PARAMETER_FILES])
290{
291 int ier;
292
293 if (numberParameterFiles > MAX_PARAMETER_FILES)
294 {
295 ier = true;
296 LOG_ERROR("LennardJones612 given too many parameter files");
297 return ier;
298 }
299
300 std::string const * paramFileDirName;
301 modelDriverCreate->GetParameterFileDirectoryName(&paramFileDirName);
302 for (int i = 0; i < numberParameterFiles; ++i)
303 {
304 std::string const * paramFileName;
305 ier = modelDriverCreate->GetParameterFileBasename(i, &paramFileName);
306 if (ier)
307 {
308 LOG_ERROR("Unable to get parameter file name");
309 return ier;
310 }
311 std::string filename = *paramFileDirName + "/" + *paramFileName;
312 parameterFilePointers[i] = fopen(filename.c_str(), "r");
313 if (parameterFilePointers[i] == 0)
314 {
315 char message[MAXLINE];
316 sprintf(message,
317 "LennardJones612 parameter file number %d cannot be opened",
318 i);
319 ier = true;
320 LOG_ERROR(message);
321 for (int j = i - 1; j >= 0; --j) { fclose(parameterFilePointers[j]); }
322 return ier;
323 }
324 }
325
326 // everything is good
327 ier = false;
328 return ier;
329}
330
331//******************************************************************************
332#undef KIM_LOGGER_OBJECT_NAME
333#define KIM_LOGGER_OBJECT_NAME modelDriverCreate
334//
335int LennardJones612Implementation::ProcessParameterFiles(
336 KIM::ModelDriverCreate * const modelDriverCreate,
337 int const /* numberParameterFiles */,
338 FILE * const parameterFilePointers[MAX_PARAMETER_FILES])
339{
340 int N, ier;
341 int endOfFileFlag = 0;
342 char spec1[MAXLINE], spec2[MAXLINE], nextLine[MAXLINE];
343 char * nextLinePtr;
344 int iIndex, jIndex, indx, iiIndex, jjIndex;
345 double nextCutoff, nextEpsilon, nextSigma;
346
347 nextLinePtr = nextLine;
348
349 getNextDataLine(
350 parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
351 ier = sscanf(nextLine, "%d %d", &N, &shift_);
352 if (ier != 2)
353 {
354 sprintf(nextLine, "unable to read first line of the parameter file");
355 ier = true;
356 LOG_ERROR(nextLine);
357 fclose(parameterFilePointers[0]);
358 return ier;
359 }
360 numberModelSpecies_ = N;
361 numberUniqueSpeciesPairs_
362 = ((numberModelSpecies_ + 1) * numberModelSpecies_) / 2;
363 AllocateParameterMemory();
364
365 // set all values in the arrays to -1 for mixing later
366 for (int i = 0; i < ((N + 1) * N / 2); i++)
367 {
368 cutoffs_[i] = -1;
369 epsilons_[i] = -1;
370 sigmas_[i] = -1;
371 }
372
373
374 // keep track of known species
375 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
376 modelSpeciesMap;
377 std::vector<KIM::SpeciesName> speciesNameVector;
378 int index = 0;
379
380 // Read and process data lines
381 getNextDataLine(
382 parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
383 while (endOfFileFlag == 0)
384 {
385 ier = sscanf(nextLine,
386 "%s %s %lg %lg %lg",
387 spec1,
388 spec2,
389 &nextCutoff,
390 &nextEpsilon,
391 &nextSigma);
392 if (ier != 5)
393 {
394 sprintf(nextLine, "error reading lines of the parameter file");
395 LOG_ERROR(nextLine);
396 return true;
397 }
398
399 // convert species strings to proper type instances
400 KIM::SpeciesName const specName1(spec1);
401 KIM::SpeciesName const specName2(spec2);
402
403 // check for new species
404 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
405 const_iterator iIter
406 = modelSpeciesMap.find(specName1);
407 if (iIter == modelSpeciesMap.end())
408 {
409 modelSpeciesMap[specName1] = index;
410 modelSpeciesCodeList_.push_back(index);
411 speciesNameVector.push_back(specName1);
412
413 ier = modelDriverCreate->SetSpeciesCode(specName1, index);
414 if (ier) return ier;
415 iIndex = index;
416 index++;
417 }
418 else
419 {
420 iIndex = modelSpeciesMap[specName1];
421 }
422 std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>::
423 const_iterator jIter
424 = modelSpeciesMap.find(specName2);
425 if (jIter == modelSpeciesMap.end())
426 {
427 modelSpeciesMap[specName2] = index;
428 modelSpeciesCodeList_.push_back(index);
429 speciesNameVector.push_back(specName2);
430
431 ier = modelDriverCreate->SetSpeciesCode(specName2, index);
432 if (ier) return ier;
433 jIndex = index;
434 index++;
435 }
436 else
437 {
438 jIndex = modelSpeciesMap[specName2];
439 }
440
441 if (iIndex >= jIndex)
442 {
443 indx = jIndex * N + iIndex - (jIndex * jIndex + jIndex) / 2;
444 }
445 else
446 {
447 indx = iIndex * N + jIndex - (iIndex * iIndex + iIndex) / 2;
448 }
449 cutoffs_[indx] = nextCutoff;
450 epsilons_[indx] = nextEpsilon;
451 sigmas_[indx] = nextSigma;
452
453 getNextDataLine(
454 parameterFilePointers[0], nextLinePtr, MAXLINE, &endOfFileFlag);
455 }
456
457 // check that we got all like - like pairs
458 std::stringstream ss;
459 ss << "There are not values for like-like pairs of:";
460 for (int i = 0; i < N; i++)
461 {
462 if (cutoffs_[(i * N + i - (i * i + i) / 2)] == -1)
463 {
464 ss << " ";
465 ss << (speciesNameVector[i].ToString()).c_str();
466 ier = -1;
467 }
468 }
469 if (ier == -1)
470 {
471 LOG_ERROR(ss);
472 return true;
473 }
474
475 // Perform Mixing if nessisary
476 for (int jIndex = 0; jIndex < N; jIndex++)
477 {
478 jjIndex = (jIndex * N + jIndex - (jIndex * jIndex + jIndex) / 2);
479 for (int iIndex = (jIndex + 1); iIndex < N; iIndex++)
480 {
481 indx = jIndex * N + iIndex - (jIndex * jIndex + jIndex) / 2;
482 if (cutoffs_[indx] == -1)
483 {
484 iiIndex = (iIndex * N + iIndex - (iIndex * iIndex + iIndex) / 2);
485 epsilons_[indx] = sqrt(epsilons_[iiIndex] * epsilons_[jjIndex]);
486 sigmas_[indx] = (sigmas_[iiIndex] + sigmas_[jjIndex]) / 2.0;
487 cutoffs_[indx] = (cutoffs_[iiIndex] + cutoffs_[jjIndex]) / 2.0;
488 }
489 }
490 }
491
492 // everything is good
493 ier = false;
494 return ier;
495}
496
497//******************************************************************************
498void LennardJones612Implementation::getNextDataLine(FILE * const filePtr,
499 char * nextLinePtr,
500 int const maxSize,
501 int * endOfFileFlag)
502{
503 do {
504 if (fgets(nextLinePtr, maxSize, filePtr) == NULL)
505 {
506 *endOfFileFlag = 1;
507 break;
508 }
509 while ((nextLinePtr[0] == ' ' || nextLinePtr[0] == '\t')
510 || (nextLinePtr[0] == '\n' || nextLinePtr[0] == '\r'))
511 {
512 nextLinePtr = (nextLinePtr + 1);
513 }
514 } while ((strncmp("#", nextLinePtr, 1) == 0) || (strlen(nextLinePtr) == 0));
515}
516
517//******************************************************************************
518void LennardJones612Implementation::CloseParameterFiles(
519 int const numberParameterFiles,
520 FILE * const parameterFilePointers[MAX_PARAMETER_FILES])
521{
522 for (int i = 0; i < numberParameterFiles; ++i)
523 fclose(parameterFilePointers[i]);
524}
525
526//******************************************************************************
527#undef KIM_LOGGER_OBJECT_NAME
528#define KIM_LOGGER_OBJECT_NAME modelDriverCreate
529//
530int LennardJones612Implementation::ConvertUnits(
531 KIM::ModelDriverCreate * const modelDriverCreate,
532 KIM::LengthUnit const requestedLengthUnit,
533 KIM::EnergyUnit const requestedEnergyUnit,
534 KIM::ChargeUnit const requestedChargeUnit,
535 KIM::TemperatureUnit const requestedTemperatureUnit,
536 KIM::TimeUnit const requestedTimeUnit)
537{
538 int ier;
539
540 // define default base units
546
547 // changing units of cutoffs and sigmas
548 double convertLength = 1.0;
550 fromEnergy,
551 fromCharge,
552 fromTemperature,
553 fromTime,
554 requestedLengthUnit,
555 requestedEnergyUnit,
556 requestedChargeUnit,
557 requestedTemperatureUnit,
558 requestedTimeUnit,
559 1.0,
560 0.0,
561 0.0,
562 0.0,
563 0.0,
564 &convertLength);
565 if (ier)
566 {
567 LOG_ERROR("Unable to convert length unit");
568 return ier;
569 }
570 if (convertLength != ONE)
571 {
572 for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
573 {
574 cutoffs_[i] *= convertLength; // convert to active units
575 sigmas_[i] *= convertLength; // convert to active units
576 }
577 }
578 // changing units of epsilons
579 double convertEnergy = 1.0;
581 fromEnergy,
582 fromCharge,
583 fromTemperature,
584 fromTime,
585 requestedLengthUnit,
586 requestedEnergyUnit,
587 requestedChargeUnit,
588 requestedTemperatureUnit,
589 requestedTimeUnit,
590 0.0,
591 1.0,
592 0.0,
593 0.0,
594 0.0,
595 &convertEnergy);
596 if (ier)
597 {
598 LOG_ERROR("Unable to convert energy unit");
599 return ier;
600 }
601 if (convertEnergy != ONE)
602 {
603 for (int i = 0; i < numberUniqueSpeciesPairs_; ++i)
604 {
605 epsilons_[i] *= convertEnergy; // convert to active units
606 }
607 }
608
609 // register units
610 ier = modelDriverCreate->SetUnits(requestedLengthUnit,
611 requestedEnergyUnit,
615 if (ier)
616 {
617 LOG_ERROR("Unable to set units to requested values");
618 return ier;
619 }
620
621 // everything is good
622 ier = false;
623 return ier;
624}
625
626//******************************************************************************
627int LennardJones612Implementation::RegisterKIMModelSettings(
628 KIM::ModelDriverCreate * const modelDriverCreate) const
629{
630 // register numbering
631 int error = modelDriverCreate->SetModelNumbering(KIM::NUMBERING::zeroBased);
632
633 return error;
634}
635
636//******************************************************************************
637#undef KIM_LOGGER_OBJECT_NAME
638#define KIM_LOGGER_OBJECT_NAME modelComputeArgumentsCreate
639//
640int LennardJones612Implementation::RegisterKIMComputeArgumentsSettings(
641 KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate) const
642{
643 // register arguments
644 LOG_INFORMATION("Register argument supportStatus");
645 int error = modelComputeArgumentsCreate->SetArgumentSupportStatus(
648 || modelComputeArgumentsCreate->SetArgumentSupportStatus(
651 || modelComputeArgumentsCreate->SetArgumentSupportStatus(
654 || modelComputeArgumentsCreate->SetArgumentSupportStatus(
657 || modelComputeArgumentsCreate->SetArgumentSupportStatus(
660
661
662 // register callbacks
663 LOG_INFORMATION("Register callback supportStatus");
664 error = error
665 || modelComputeArgumentsCreate->SetCallbackSupportStatus(
668 || modelComputeArgumentsCreate->SetCallbackSupportStatus(
671
672 return error;
673}
674
675//******************************************************************************
676// helper macro
677#define SNUM(x) \
678 static_cast<std::ostringstream const &>(std::ostringstream() \
679 << std::dec << x) \
680 .str()
681//******************************************************************************
682#undef KIM_LOGGER_OBJECT_NAME
683#define KIM_LOGGER_OBJECT_NAME modelDriverCreate
684//
685int LennardJones612Implementation::RegisterKIMParameters(
686 KIM::ModelDriverCreate * const modelDriverCreate)
687{
688 int ier = false;
689
690 // publish parameters (order is important)
691 ier = modelDriverCreate->SetParameterPointer(
692 1,
693 &shift_,
694 "shift",
695 "If (shift == 1), all LJ potentials are shifted to zero energy "
696 "at their respective cutoff distance. Otherwise, no shifting is "
697 "performed.");
698 if (ier)
699 {
700 LOG_ERROR("set_parameter shift");
701 return ier;
702 }
703
704 ier = modelDriverCreate->SetParameterPointer(
705 numberUniqueSpeciesPairs_,
706 cutoffs_,
707 "cutoffs",
708 "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
709 + ") "
710 "in row-major storage. Ordering is according to SpeciesCode "
711 "values. "
712 "For example, to find the parameter related to SpeciesCode 'i' and "
713 "SpeciesCode 'j' (i >= j), use (zero-based) "
714 "index = (j*N + i - (j*j + j)/2).");
715 if (ier)
716 {
717 LOG_ERROR("set_parameter cutoffs");
718 return ier;
719 }
720 ier = modelDriverCreate->SetParameterPointer(
721 numberUniqueSpeciesPairs_,
722 epsilons_,
723 "epsilons",
724 "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
725 + ") "
726 "in row-major storage. Ordering is according to SpeciesCode "
727 "values. "
728 "For example, to find the parameter related to SpeciesCode 'i' and "
729 "SpeciesCode 'j' (i >= j), use (zero-based) "
730 "index = (j*N + i - (j*j + j)/2).");
731 if (ier)
732 {
733 LOG_ERROR("set_parameter epsilons");
734 return ier;
735 }
736 ier = modelDriverCreate->SetParameterPointer(
737 numberUniqueSpeciesPairs_,
738 sigmas_,
739 "sigmas",
740 "Lower-triangular matrix (of size N=" + SNUM(numberModelSpecies_)
741 + ") "
742 "in row-major storage. Ordering is according to SpeciesCode "
743 "values. "
744 "For example, to find the parameter related to SpeciesCode 'i' and "
745 "SpeciesCode 'j' (i >= j), use (zero-based) "
746 "index = (j*N + i - (j*j + j)/2).");
747 if (ier)
748 {
749 LOG_ERROR("set_parameter sigmas");
750 return ier;
751 }
752
753 // everything is good
754 ier = false;
755 return ier;
756}
757
758//******************************************************************************
759int LennardJones612Implementation::RegisterKIMFunctions(
760 KIM::ModelDriverCreate * const modelDriverCreate) const
761{
762 int error;
763
764 // Use function pointer definitions to verify correct prototypes
772
773 // register the destroy() and reinit() functions
774 error = modelDriverCreate->SetRoutinePointer(
777 true,
778 reinterpret_cast<KIM::Function *>(destroy))
779 || modelDriverCreate->SetRoutinePointer(
782 true,
783 reinterpret_cast<KIM::Function *>(refresh))
784 || modelDriverCreate->SetRoutinePointer(
787 true,
788 reinterpret_cast<KIM::Function *>(compute))
789 || modelDriverCreate->SetRoutinePointer(
792 true,
793 reinterpret_cast<KIM::Function *>(CACreate))
794 || modelDriverCreate->SetRoutinePointer(
797 true,
798 reinterpret_cast<KIM::Function *>(CADestroy));
799 return error;
800}
801
802//******************************************************************************
803template<class ModelObj>
804int LennardJones612Implementation::SetRefreshMutableValues(
805 ModelObj * const modelObj)
806{ // use (possibly) new values of parameters to compute other quantities
807 // NOTE: This function is templated because it's called with both a
808 // modelDriverCreate object during initialization and with a
809 // modelRefresh object when the Model's parameters have been altered
810 int ier;
811
812 // update cutoffsSq, epsilons, and sigmas
813 for (int i = 0; i < numberModelSpecies_; ++i)
814 {
815 for (int j = 0; j <= i; ++j)
816 {
817 int const index = j * numberModelSpecies_ + i - (j * j + j) / 2;
818 cutoffsSq2D_[i][j] = cutoffsSq2D_[j][i]
819 = (cutoffs_[index] * cutoffs_[index]);
820 fourEpsilonSigma6_2D_[i][j] = fourEpsilonSigma6_2D_[j][i]
821 = 4.0 * epsilons_[index] * pow(sigmas_[index], 6.0);
822 fourEpsilonSigma12_2D_[i][j] = fourEpsilonSigma12_2D_[j][i]
823 = 4.0 * epsilons_[index] * pow(sigmas_[index], 12.0);
824 twentyFourEpsilonSigma6_2D_[i][j] = twentyFourEpsilonSigma6_2D_[j][i]
825 = 6.0 * fourEpsilonSigma6_2D_[i][j];
826 fortyEightEpsilonSigma12_2D_[i][j] = fortyEightEpsilonSigma12_2D_[j][i]
827 = 12.0 * fourEpsilonSigma12_2D_[i][j];
828 oneSixtyEightEpsilonSigma6_2D_[i][j]
829 = oneSixtyEightEpsilonSigma6_2D_[j][i]
830 = 7.0 * twentyFourEpsilonSigma6_2D_[i][j];
831 sixTwentyFourEpsilonSigma12_2D_[i][j]
832 = sixTwentyFourEpsilonSigma12_2D_[j][i]
833 = 13.0 * fortyEightEpsilonSigma12_2D_[i][j];
834 }
835 }
836
837 // update cutoff value in KIM API object
838 influenceDistance_ = 0.0;
839
840 for (int i = 0; i < numberModelSpecies_; i++)
841 {
842 int indexI = modelSpeciesCodeList_[i];
843
844 for (int j = 0; j < numberModelSpecies_; j++)
845 {
846 int indexJ = modelSpeciesCodeList_[j];
847
848 if (influenceDistance_ < cutoffsSq2D_[indexI][indexJ])
849 {
850 influenceDistance_ = cutoffsSq2D_[indexI][indexJ];
851 }
852 }
853 }
854
855 influenceDistance_ = sqrt(influenceDistance_);
856 modelObj->SetInfluenceDistancePointer(&influenceDistance_);
857 modelObj->SetNeighborListPointers(
858 1,
859 &influenceDistance_,
860 &modelWillNotRequestNeighborsOfNoncontributingParticles_);
861
862 // update shifts
863 // compute and set shifts2D_ check if minus sign
864 double const * const * const constFourEpsSig6_2D = fourEpsilonSigma6_2D_;
865 double const * const * const constFourEpsSig12_2D = fourEpsilonSigma12_2D_;
866 if (1 == shift_)
867 {
868 double phi;
869 for (int iSpecies = 0; iSpecies < numberModelSpecies_; iSpecies++)
870 {
871 for (int jSpecies = 0; jSpecies <= iSpecies; jSpecies++)
872 {
873 int const index = jSpecies * numberModelSpecies_ + iSpecies
874 - (jSpecies * jSpecies + jSpecies) / 2;
875 double const rij2 = cutoffs_[index] * cutoffs_[index];
876 double const r2iv = 1.0 / rij2;
877 double const r6iv = r2iv * r2iv * r2iv;
879 shifts2D_[iSpecies][jSpecies] = shifts2D_[jSpecies][iSpecies] = phi;
880 }
881 }
882 }
883
884 // everything is good
885 ier = false;
886 return ier;
887}
888
889//******************************************************************************
890#undef KIM_LOGGER_OBJECT_NAME
891#define KIM_LOGGER_OBJECT_NAME modelComputeArguments
892//
893int LennardJones612Implementation::SetComputeMutableValues(
894 KIM::ModelComputeArguments const * const modelComputeArguments,
895 bool & isComputeProcess_dEdr,
896 bool & isComputeProcess_d2Edr2,
897 bool & isComputeEnergy,
898 bool & isComputeForces,
899 bool & isComputeParticleEnergy,
900 bool & isComputeVirial,
901 bool & isComputeParticleVirial,
902 int const *& particleSpeciesCodes,
903 int const *& particleContributing,
904 VectorOfSizeDIM const *& coordinates,
905 double *& energy,
906 double *& particleEnergy,
907 VectorOfSizeDIM *& forces,
908 VectorOfSizeSix *& virial,
909 VectorOfSizeSix *& particleVirial)
910{
911 int ier = true;
912
913 // get compute flags
914 int compProcess_dEdr;
915 int compProcess_d2Edr2;
916
917 modelComputeArguments->IsCallbackPresent(
919 modelComputeArguments->IsCallbackPresent(
921
922 isComputeProcess_dEdr = compProcess_dEdr;
923 isComputeProcess_d2Edr2 = compProcess_d2Edr2;
924
925 int const * numberOfParticles;
926 ier = modelComputeArguments->GetArgumentPointer(
928 || modelComputeArguments->GetArgumentPointer(
930 &particleSpeciesCodes)
931 || modelComputeArguments->GetArgumentPointer(
933 &particleContributing)
934 || modelComputeArguments->GetArgumentPointer(
936 (double const **) &coordinates)
937 || modelComputeArguments->GetArgumentPointer(
939 || modelComputeArguments->GetArgumentPointer(
941 || modelComputeArguments->GetArgumentPointer(
943 (double const **) &forces)
944 || modelComputeArguments->GetArgumentPointer(
946 (double const **) &virial)
947 || modelComputeArguments->GetArgumentPointer(
949 (double const **) &particleVirial);
950 if (ier)
951 {
952 LOG_ERROR("GetArgumentPointer");
953 return ier;
954 }
955
956 isComputeEnergy = (energy != NULL);
957 isComputeParticleEnergy = (particleEnergy != NULL);
958 isComputeForces = (forces != NULL);
959 isComputeVirial = (virial != NULL);
960 isComputeParticleVirial = (particleVirial != NULL);
961
962 // update values
963 cachedNumberOfParticles_ = *numberOfParticles;
964
965 // everything is good
966 ier = false;
967 return ier;
968}
969
970//******************************************************************************
971#undef KIM_LOGGER_OBJECT_NAME
972#define KIM_LOGGER_OBJECT_NAME modelCompute
973int LennardJones612Implementation::CheckParticleSpeciesCodes(
974 KIM::ModelCompute const * const modelCompute,
975 int const * const particleSpeciesCodes) const
976{
977 int ier;
978 for (int i = 0; i < cachedNumberOfParticles_; ++i)
979 {
980 if ((particleSpeciesCodes[i] < 0)
981 || (particleSpeciesCodes[i] >= numberModelSpecies_))
982 {
983 ier = true;
984 LOG_ERROR("unsupported particle species codes detected");
985 return ier;
986 }
987 }
988
989 // everything is good
990 ier = false;
991 return ier;
992}
993
994//******************************************************************************
995int LennardJones612Implementation::GetComputeIndex(
996 const bool & isComputeProcess_dEdr,
997 const bool & isComputeProcess_d2Edr2,
998 const bool & isComputeEnergy,
999 const bool & isComputeForces,
1000 const bool & isComputeParticleEnergy,
1001 const bool & isComputeVirial,
1002 const bool & isComputeParticleVirial,
1003 const bool & isShift) const
1004{
1005 // const int processdE = 2;
1006 const int processd2E = 2;
1007 const int energy = 2;
1008 const int force = 2;
1009 const int particleEnergy = 2;
1010 const int virial = 2;
1011 const int particleVirial = 2;
1012 const int shift = 2;
1013
1014
1015 int index = 0;
1016
1017 // processdE
1018 index += (int(isComputeProcess_dEdr)) * processd2E * energy * force
1019 * particleEnergy * virial * particleVirial * shift;
1020
1021 // processd2E
1022 index += (int(isComputeProcess_d2Edr2)) * energy * force * particleEnergy
1023 * virial * particleVirial * shift;
1024
1025 // energy
1026 index += (int(isComputeEnergy)) * force * particleEnergy * virial
1027 * particleVirial * shift;
1028
1029 // force
1030 index += (int(isComputeForces)) * particleEnergy * virial * particleVirial
1031 * shift;
1032
1033 // particleEnergy
1034 index += (int(isComputeParticleEnergy)) * virial * particleVirial * shift;
1035
1036 // virial
1037 index += (int(isComputeVirial)) * particleVirial * shift;
1038
1039 // particleVirial
1040 index += (int(isComputeParticleVirial)) * shift;
1041
1042 // shift
1043 index += (int(isShift));
1044
1045 return index;
1046}
1047
1048//******************************************************************************
1049void LennardJones612Implementation::ProcessVirialTerm(
1050 const double & dEidr,
1051 const double & rij,
1052 const double * const r_ij,
1053 const int & /* i */,
1054 const int & /* j */,
1055 VectorOfSizeSix virial) const
1056{
1057 double const v = dEidr / rij;
1058
1059 virial[0] += v * r_ij[0] * r_ij[0];
1060 virial[1] += v * r_ij[1] * r_ij[1];
1061 virial[2] += v * r_ij[2] * r_ij[2];
1062 virial[3] += v * r_ij[1] * r_ij[2];
1063 virial[4] += v * r_ij[0] * r_ij[2];
1064 virial[5] += v * r_ij[0] * r_ij[1];
1065}
1066
1067//******************************************************************************
1068void LennardJones612Implementation::ProcessParticleVirialTerm(
1069 const double & dEidr,
1070 const double & rij,
1071 const double * const r_ij,
1072 const int & i,
1073 const int & j,
1074 VectorOfSizeSix * const particleVirial) const
1075{
1076 double const v = dEidr / rij;
1077 VectorOfSizeSix vir;
1078
1079 vir[0] = 0.5 * v * r_ij[0] * r_ij[0];
1080 vir[1] = 0.5 * v * r_ij[1] * r_ij[1];
1081 vir[2] = 0.5 * v * r_ij[2] * r_ij[2];
1082 vir[3] = 0.5 * v * r_ij[1] * r_ij[2];
1083 vir[4] = 0.5 * v * r_ij[0] * r_ij[2];
1084 vir[5] = 0.5 * v * r_ij[0] * r_ij[1];
1085
1086 for (int k = 0; k < 6; ++k)
1087 {
1088 particleVirial[i][k] += vir[k];
1089 particleVirial[j][k] += vir[k];
1090 }
1091}
1092
1093//==============================================================================
1094//
1095// Implementation of helper functions
1096//
1097//==============================================================================
1098
1099//******************************************************************************
1100void AllocateAndInitialize2DArray(double **& arrayPtr,
1101 int const extentZero,
1102 int const extentOne)
1103{ // allocate memory and set pointers
1104 arrayPtr = new double *[extentZero];
1105 arrayPtr[0] = new double[extentZero * extentOne];
1106 for (int i = 1; i < extentZero; ++i)
1107 {
1108 arrayPtr[i] = arrayPtr[i - 1] + extentOne;
1109 }
1110
1111 // initialize
1112 for (int i = 0; i < extentZero; ++i)
1113 {
1114 for (int j = 0; j < extentOne; ++j) { arrayPtr[i][j] = 0.0; }
1115 }
1116}
1117
1118//******************************************************************************
1119void Deallocate2DArray(double **& arrayPtr)
1120{ // deallocate memory
1121 if (arrayPtr != NULL) delete[] arrayPtr[0];
1122 delete[] arrayPtr;
1123
1124 // nullify pointer
1125 arrayPtr = NULL;
1126}
#define LOG_ERROR(message)
Convenience macro for ERROR Log entries with compile-time optimization.
#define LOG_INFORMATION(message)
Convenience macro for INFORMATION Log entries with compile-time optimization.
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
void Deallocate2DArray(double **&arrayPtr)
void AllocateAndInitialize2DArray(double **&arrayPtr, int const extentZero, int const extentOne)
double VectorOfSizeDIM[DIMENSION]
#define MAX_PARAMETER_FILES
void Deallocate2DArray(double **&arrayPtr)
#define LENNARD_JONES_PHI(exshift)
double VectorOfSizeSix[6]
An Extensible Enumeration for the ChargeUnit's supported by the KIM API.
An Extensible Enumeration for the EnergyUnit's supported by the KIM API.
An Extensible Enumeration for the LengthUnit's supported by the KIM API.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
int SetCallbackSupportStatus(ComputeCallbackName const computeCallbackName, SupportStatus const supportStatus)
Set the SupportStatus of a ComputeCallbackName.
int SetArgumentSupportStatus(ComputeArgumentName const computeArgumentName, SupportStatus const supportStatus)
Set the SupportStatus of a ComputeArgumentName.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
int GetArgumentPointer(ComputeArgumentName const computeArgumentName, int const **const ptr) const
Get the data pointer for a ComputeArgumentName.
int IsCallbackPresent(ComputeCallbackName const computeCallbackName, int *const present) const
Determine if the Simulator has provided a non-NULL function pointer for a ComputeCallbackName of inte...
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
void GetNumberOfParameterFiles(int *const numberOfParameterFiles) const
Get the number of parameter files provided by the parameterized model.
int SetModelNumbering(Numbering const numbering)
Set the Model's particle Numbering.
void GetParameterFileDirectoryName(std::string const **const directoryName) const
Get absolute path name of the temporary directory where parameter files provided by the model are wri...
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
Set the Model's base unit values.
int SetRoutinePointer(ModelRoutineName const modelRoutineName, LanguageName const languageName, int const required, Function *const fptr)
Set the function pointer for the ModelRoutineName of interest.
int GetParameterFileBasename(int const index, std::string const **const parameterFileBasename) const
Get a particular parameter file basename. The file is located in the Model's parameter file directory...
int SetSpeciesCode(SpeciesName const speciesName, int const code)
Set integer code for supported SpeciesName.
int SetParameterPointer(int const extent, int *const ptr, std::string const &name, std::string const &description)
Set the next parameter data pointer to be provided by the model.
static int ConvertUnit(LengthUnit const fromLengthUnit, EnergyUnit const fromEnergyUnit, ChargeUnit const fromChargeUnit, TemperatureUnit const fromTemperatureUnit, TimeUnit const fromTimeUnit, LengthUnit const toLengthUnit, EnergyUnit const toEnergyUnit, ChargeUnit const toChargeUnit, TemperatureUnit const toTemperatureUnit, TimeUnit const toTimeUnit, double const lengthExponent, double const energyExponent, double const chargeExponent, double const temperatureExponent, double const timeExponent, double *const conversionFactor)
Get the multiplicative factor to convert between a derived unit represented in two different sets of ...
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::R...
An Extensible Enumeration for the SpeciesName's supported by the KIM API.
An Extensible Enumeration for the TemperatureUnit's supported by the KIM API.
An Extensible Enumeration for the TimeUnit's supported by the KIM API.
static int ComputeArgumentsCreate(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
static int Destroy(KIM::ModelDestroy *const modelDestroy)
static int Refresh(KIM::ModelRefresh *const modelRefresh)
static int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
static int ComputeArgumentsDestroy(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
int Refresh(KIM::ModelRefresh *const modelRefresh)
int ComputeArgumentsCreate(KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate) const
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)
int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
int ComputeArgumentsDestroy(KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy) const
ChargeUnit const e
The standard electron unit of charge.
ChargeUnit const unused
Indicates that a ChargeUnit is not used.
ComputeArgumentName const partialParticleEnergy
The standard partialParticleEnergy argument.
ComputeArgumentName const coordinates
The standard coordinates argument.
ComputeArgumentName const numberOfParticles
The standard numberOfParticles argument.
ComputeArgumentName const particleSpeciesCodes
The standard particleSpeciesCodes argument.
ComputeArgumentName const partialEnergy
The standard partialEnergy argument.
ComputeArgumentName const partialParticleVirial
The standard partialParticleVirial argument.
ComputeArgumentName const partialForces
The standard partialForces argument.
ComputeArgumentName const particleContributing
The standard particleContributing argument.
ComputeArgumentName const partialVirial
The standard partialVirial argument.
ComputeCallbackName const ProcessDEDrTerm
The standard ProcessDEDrTerm callback.
ComputeCallbackName const ProcessD2EDr2Term
The standard ProcessD2EDr2Term callback.
EnergyUnit const eV
The standard electronvolt unit of energy.
LanguageName const cpp
The standard cpp language.
LengthUnit const A
The standard angstrom unit of length.
ModelRoutineName const Destroy
The standard Destroy routine.
ModelRoutineName const Compute
The standard Compute routine.
ModelRoutineName const Refresh
The standard Refresh routine.
ModelRoutineName const ComputeArgumentsDestroy
The standard ComputeArgumentsDestroy routine.
ModelRoutineName const ComputeArgumentsCreate
The standard ComputeArgumentsCreate routine.
Numbering const zeroBased
The standard zeroBased numbering.
SpeciesName const N
The standard Nitrogen species.
SupportStatus const optional
The standard optional status.
TemperatureUnit const unused
Indicates that a TemperatureUnit is not used.
TemperatureUnit const K
The standard Kelvin unit of temperature.
TimeUnit const unused
Indicates that a TimeUnit is not used.
TimeUnit const ps
The standard picosecond unit of time.
int ModelComputeFunction(ModelCompute const *const modelCompute, ModelComputeArguments const *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::Compute routine.
int ModelComputeArgumentsCreateFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsCreate routine.
int ModelDestroyFunction(ModelDestroy *const modelDestroy)
Prototype for MODEL_ROUTINE_NAME::Destroy routine.
int ModelRefreshFunction(ModelRefresh *const modelRefresh)
Prototype for MODEL_ROUTINE_NAME::Refresh routine.
void() Function(void)
Generic function type.
int ModelComputeArgumentsDestroyFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsDestroy routine.
recursive subroutine, public refresh(model_refresh_handle, ierr)
recursive subroutine, public destroy(model_destroy_handle, ierr)