kim-api  2.1.4-git+v2.1.3-git-3-g4c859c7f.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
LennardJones_Ar.cpp
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) 2018--2019, Regents of the University of Minnesota.
23 // All rights reserved.
24 //
25 // Contributors:
26 // Ryan S. Elliott
27 //
28 
29 #include "KIM_LogMacros.hpp"
30 #include "KIM_ModelHeaders.hpp"
31 #include <math.h>
32 
33 #define DIMENSION 3
34 
35 
36 namespace
37 {
39 {
40  public:
41  //****************************************************************************
42  LennardJones_Ar(KIM::ModelCreate * const modelCreate,
43  KIM::LengthUnit const requestedLengthUnit,
44  KIM::EnergyUnit const requestedEnergyUnit,
45  KIM::ChargeUnit const requestedChargeUnit,
46  KIM::TemperatureUnit const requestedTemperatureUnit,
47  KIM::TimeUnit const requestedTimeUnit,
48  int * const error) :
49  epsilon_(0.0104),
50  sigma_(3.4000),
51  influenceDistance_(8.1500),
52  cutoff_(influenceDistance_),
53  cutoffSq_(cutoff_ * cutoff_),
54  modelWillNotRequestNeighborsOfNoncontributingParticles_(1)
55  {
56  *error = ConvertUnits(modelCreate,
57  requestedLengthUnit,
58  requestedEnergyUnit,
59  requestedChargeUnit,
60  requestedTemperatureUnit,
61  requestedTimeUnit);
62  if (*error) return;
63 
65 
66  modelCreate->SetInfluenceDistancePointer(&influenceDistance_);
67  modelCreate->SetNeighborListPointers(
68  1, &cutoff_, &modelWillNotRequestNeighborsOfNoncontributingParticles_);
69 
70  modelCreate->SetSpeciesCode(KIM::SPECIES_NAME::Ar, 0);
71 
72  // use function pointer declarations to verify prototypes
79 
80  *error = modelCreate->SetRoutinePointer(
83  true,
84  reinterpret_cast<KIM::Function *>(CACreate))
85  || modelCreate->SetRoutinePointer(
88  true,
89  reinterpret_cast<KIM::Function *>(compute))
90  || modelCreate->SetRoutinePointer(
93  true,
94  reinterpret_cast<KIM::Function *>(CADestroy))
95  || modelCreate->SetRoutinePointer(
98  true,
99  reinterpret_cast<KIM::Function *>(destroy));
100  if (*error) return;
101 
102  // everything is good
103  *error = false;
104  return;
105  };
106 
107  //****************************************************************************
109 
110  //****************************************************************************
111  // no need to make these "extern" since KIM will only access them
112  // via function pointers. "static" is required so that there is not
113  // an implicit this pointer added to the prototype by the C++ compiler
114  static int Destroy(KIM::ModelDestroy * const modelDestroy)
115  {
116  LennardJones_Ar * model;
117  modelDestroy->GetModelBufferPointer(reinterpret_cast<void **>(&model));
118 
119  if (model != NULL)
120  {
121  // delete object itself
122  delete model;
123  }
124 
125  // everything is good
126  return false;
127  }
128 
129  //****************************************************************************
130 #undef KIM_LOGGER_OBJECT_NAME
131 #define KIM_LOGGER_OBJECT_NAME modelCompute
132  //
133  static int
134  Compute(KIM::ModelCompute const * const modelCompute,
135  KIM::ModelComputeArguments const * const modelComputeArguments)
136  {
137  int const * numberOfParticlesPointer;
138  int const * particleSpeciesCodes;
139  int const * particleContributing;
140  double const * coordinates;
141  double * partialEnergy;
142  double * partialForces;
143 
144  LennardJones_Ar * lj;
145  modelCompute->GetModelBufferPointer(reinterpret_cast<void **>(&lj));
146  double const epsilon = lj->epsilon_;
147  double const sigma = lj->sigma_;
148  double const cutoffSq = lj->cutoffSq_;
149 
150  int error = modelComputeArguments->GetArgumentPointer(
152  &numberOfParticlesPointer)
153  || modelComputeArguments->GetArgumentPointer(
155  &particleSpeciesCodes)
156  || modelComputeArguments->GetArgumentPointer(
158  &particleContributing)
159  || modelComputeArguments->GetArgumentPointer(
161  (double const **) &coordinates)
162  || modelComputeArguments->GetArgumentPointer(
164  || modelComputeArguments->GetArgumentPointer(
166  (double const **) &partialForces);
167  if (error)
168  {
169  LOG_ERROR("Unable to get argument pointers");
170  return error;
171  }
172 
173  int const numberOfParticles = *numberOfParticlesPointer;
174 
175  // initialize energy and forces
176  *partialEnergy = 0.0;
177  int const extent = numberOfParticles * DIMENSION;
178  for (int i = 0; i < extent; ++i) { partialForces[i] = 0.0; }
179 
180  int jContributing;
181  int i, j, jj, numberOfNeighbors;
182  int const * neighbors;
183  double phi;
184  double xcoord, ycoord, zcoord;
185  double xrij, yrij, zrij;
186  double rij2;
187  double r2inv, r6inv, dphiByR, dEidrByR;
188  double xdf, ydf, zdf;
189  double const fortyEight = 48.0 * epsilon * pow(sigma, 12.0);
190  double const twentyFour = 24.0 * epsilon * pow(sigma, 6.0);
191  double const four12 = 4.0 * epsilon * pow(sigma, 12.0);
192  double const four6 = 4.0 * epsilon * pow(sigma, 6.0);
193  for (i = 0; i < numberOfParticles; ++i)
194  {
195  if (particleContributing[i])
196  {
197  xcoord = coordinates[i * DIMENSION + 0];
198  ycoord = coordinates[i * DIMENSION + 1];
199  zcoord = coordinates[i * DIMENSION + 2];
200 
201  modelComputeArguments->GetNeighborList(
202  0, i, &numberOfNeighbors, &neighbors);
203 
204  for (jj = 0; jj < numberOfNeighbors; ++jj)
205  {
206  j = neighbors[jj];
207  jContributing = particleContributing[j];
208 
209  if (!(jContributing && (j < i)))
210  {
211  xrij = coordinates[j * DIMENSION + 0] - xcoord;
212  yrij = coordinates[j * DIMENSION + 1] - ycoord;
213  zrij = coordinates[j * DIMENSION + 2] - zcoord;
214 
215  rij2 = xrij * xrij + yrij * yrij + zrij * zrij;
216 
217  if (rij2 < cutoffSq)
218  {
219  r2inv = 1.0 / rij2;
220  r6inv = r2inv * r2inv * r2inv;
221  phi = 0.5 * r6inv * (four12 * r6inv - four6);
222  dphiByR = r6inv * (twentyFour - fortyEight * r6inv) * r2inv;
223 
224  *partialEnergy += phi;
225  if (jContributing)
226  {
227  *partialEnergy += phi;
228  dEidrByR = dphiByR;
229  }
230  else
231  {
232  dEidrByR = 0.5 * dphiByR;
233  }
234 
235  xdf = dEidrByR * xrij;
236  ydf = dEidrByR * yrij;
237  zdf = dEidrByR * zrij;
238  partialForces[i * DIMENSION + 0] += xdf;
239  partialForces[i * DIMENSION + 1] += ydf;
240  partialForces[i * DIMENSION + 2] += zdf;
241  partialForces[j * DIMENSION + 0] -= xdf;
242  partialForces[j * DIMENSION + 1] -= ydf;
243  partialForces[j * DIMENSION + 2] -= zdf;
244  } // if (rij2 < cutoffSq_)
245  } // if (i < j)
246  } // for jj
247  } // if (particleContributing[i])
248  } // for i
249 
250  // everything is good
251  return false;
252  };
253 
254  //****************************************************************************
256  KIM::ModelCompute const * const /* modelCompute */,
257  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
258  {
259  // register arguments
260  int error = modelComputeArgumentsCreate->SetArgumentSupportStatus(
263  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
266 
267  // register callbacks
268  //
269  // none
270 
271  return error;
272  }
273 
274  //****************************************************************************
275  static int
276  ComputeArgumentsDestroy(KIM::ModelCompute const * const /* modelCompute */,
278  /* modelComputeArgumentsDestroy */)
279  {
280  // nothing further to do
281 
282  return false;
283  }
284 
285  private:
286  //****************************************************************************
287  // Member variables
288  double epsilon_;
289  double sigma_;
290  double influenceDistance_;
291  double cutoff_;
292  double cutoffSq_;
293  int const modelWillNotRequestNeighborsOfNoncontributingParticles_;
294 
295  //****************************************************************************
296 #undef KIM_LOGGER_OBJECT_NAME
297 #define KIM_LOGGER_OBJECT_NAME modelCreate
298  //
299  int ConvertUnits(KIM::ModelCreate * const modelCreate,
300  KIM::LengthUnit const requestedLengthUnit,
301  KIM::EnergyUnit const requestedEnergyUnit,
302  KIM::ChargeUnit const requestedChargeUnit,
303  KIM::TemperatureUnit const requestedTemperatureUnit,
304  KIM::TimeUnit const requestedTimeUnit)
305  {
306  int ier;
307 
308  // define default base units
314 
315  // changing units of cutoffs and sigmas
316  double convertLength = 1.0;
317  ier = KIM::ModelCreate::ConvertUnit(fromLength,
318  fromEnergy,
319  fromCharge,
320  fromTemperature,
321  fromTime,
322  requestedLengthUnit,
323  requestedEnergyUnit,
324  requestedChargeUnit,
325  requestedTemperatureUnit,
326  requestedTimeUnit,
327  1.0,
328  0.0,
329  0.0,
330  0.0,
331  0.0,
332  &convertLength);
333  if (ier)
334  {
335  LOG_ERROR("Unable to convert length unit");
336  return ier;
337  }
338  influenceDistance_ *= convertLength; // convert to active units
339  cutoff_ = influenceDistance_;
340  cutoffSq_ = cutoff_ * cutoff_;
341  sigma_ *= convertLength; // convert to active units
342 
343  // changing units of epsilons
344  double convertEnergy = 1.0;
345  ier = KIM::ModelCreate::ConvertUnit(fromLength,
346  fromEnergy,
347  fromCharge,
348  fromTemperature,
349  fromTime,
350  requestedLengthUnit,
351  requestedEnergyUnit,
352  requestedChargeUnit,
353  requestedTemperatureUnit,
354  requestedTimeUnit,
355  0.0,
356  1.0,
357  0.0,
358  0.0,
359  0.0,
360  &convertEnergy);
361  if (ier)
362  {
363  LOG_ERROR("Unable to convert energy unit");
364  return ier;
365  }
366  epsilon_ *= convertEnergy; // convert to active units
367 
368  // register units
369  ier = modelCreate->SetUnits(requestedLengthUnit,
370  requestedEnergyUnit,
374  if (ier)
375  {
376  LOG_ERROR("Unable to set units to requested values");
377  return ier;
378  }
379 
380  // everything is good
381  ier = false;
382  return ier;
383  }
384 };
385 
386 } // namespace
387 
388 extern "C" {
389 //******************************************************************************
390 int model_create(KIM::ModelCreate * const modelCreate,
391  KIM::LengthUnit const requestedLengthUnit,
392  KIM::EnergyUnit const requestedEnergyUnit,
393  KIM::ChargeUnit const requestedChargeUnit,
394  KIM::TemperatureUnit const requestedTemperatureUnit,
395  KIM::TimeUnit const requestedTimeUnit)
396 {
397  int error;
398 
399  LennardJones_Ar * const model = new LennardJones_Ar(modelCreate,
400  requestedLengthUnit,
401  requestedEnergyUnit,
402  requestedChargeUnit,
403  requestedTemperatureUnit,
404  requestedTimeUnit,
405  &error);
406  if (error)
407  {
408  // constructor already reported the error
409  delete model;
410  return error;
411  }
412 
413  // register pointer to LennardJones_Ar object in moedelCreate object
414  modelCreate->SetModelBufferPointer(reinterpret_cast<void *>(model));
415 
416  // everything is good
417  return false;
418 }
419 } // extern "C"
int ModelComputeArgumentsDestroyFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsDestroy routine.
TemperatureUnit const unused
Indicates that a TemperatureUnit is not used.
An Extensible Enumeration for the TemperatureUnit&#39;s supported by the KIM API.
SpeciesName const Ar
The standard Argon species.
An Extensible Enumeration for the TimeUnit&#39;s supported by the KIM API.
int ModelDestroyFunction(ModelDestroy *const modelDestroy)
Prototype for MODEL_ROUTINE_NAME::Destroy routine.
int SetSpeciesCode(SpeciesName const speciesName, int const code)
Set integer code for supported SpeciesName.
int SetModelNumbering(Numbering const numbering)
Set the Model&#39;s particle Numbering.
recursive subroutine, public destroy(model_destroy_handle, ierr)
ComputeArgumentName const coordinates
The standard coordinates argument.
#define LOG_ERROR(message)
Convenience macro for ERROR Log entries with compile-time optimization.
ModelRoutineName const Destroy
The standard Destroy routine.
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
Set the Model&#39;s base unit values.
ChargeUnit const unused
Indicates that a ChargeUnit is not used.
void GetModelBufferPointer(void **const ptr) const
Get the Model&#39;s buffer pointer within the Model object.
void SetInfluenceDistancePointer(double const *const influenceDistance)
Set the Model&#39;s influence distance data pointer.
int SetArgumentSupportStatus(ComputeArgumentName const computeArgumentName, SupportStatus const supportStatus)
Set the SupportStatus of a ComputeArgumentName.
int ModelComputeFunction(ModelCompute const *const modelCompute, ModelComputeArguments const *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::Compute routine.
void SetModelBufferPointer(void *const ptr)
Set the Model&#39;s buffer pointer within the Model object.
SupportStatus const required
The standard required status.
An Extensible Enumeration for the LengthUnit&#39;s supported by the KIM API.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::D...
int GetArgumentPointer(ComputeArgumentName const computeArgumentName, int const **const ptr) const
Get the data pointer for a ComputeArgumentName.
int ModelComputeArgumentsCreateFunction(ModelCompute const *const modelCompute, ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
Prototype for MODEL_ROUTINE_NAME::ComputeArgumentsCreate routine.
ComputeArgumentName const particleContributing
The standard particleContributing argument.
ComputeArgumentName const partialEnergy
The standard partialEnergy argument.
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
ModelRoutineName const ComputeArgumentsCreate
The standard ComputeArgumentsCreate routine.
ModelRoutineName const ComputeArgumentsDestroy
The standard ComputeArgumentsDestroy routine.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
void SetNeighborListPointers(int const numberOfNeighborLists, double const *const cutoffs, int const *const modelWillNotRequestNeighborsOfNoncontributingParticles)
Set the Model&#39;s neighbor list data pointers.
#define DIMENSION
An Extensible Enumeration for the EnergyUnit&#39;s supported by the KIM API.
LanguageName const cpp
The standard cpp language.
ModelRoutineName const Compute
The standard Compute routine.
LengthUnit const A
The standard angstrom unit of length.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
ComputeArgumentName const partialForces
The standard partialForces argument.
void GetModelBufferPointer(void **const ptr) const
Get the Model&#39;s buffer pointer within the Model object.
EnergyUnit const eV
The standard electronvolt unit of energy.
ComputeArgumentName const particleSpeciesCodes
The standard particleSpeciesCodes argument.
Provides the interface to a KIM API ComputeArguments object for use by models within their MODEL_ROUT...
An Extensible Enumeration for the ChargeUnit&#39;s supported by the KIM API.
ComputeArgumentName const numberOfParticles
The standard numberOfParticles argument.
int GetNeighborList(int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle) const
Get the neighbor list for a particle of interest corresponding to a particular neighbor list cutoff d...
int model_create(KIM::ModelCreate *const modelCreate, KIM::LengthUnit const requestedLengthUnit, KIM::EnergyUnit const requestedEnergyUnit, KIM::ChargeUnit const requestedChargeUnit, KIM::TemperatureUnit const requestedTemperatureUnit, KIM::TimeUnit const requestedTimeUnit)
LennardJones_Ar(KIM::ModelCreate *const modelCreate, KIM::LengthUnit const requestedLengthUnit, KIM::EnergyUnit const requestedEnergyUnit, KIM::ChargeUnit const requestedChargeUnit, KIM::TemperatureUnit const requestedTemperatureUnit, KIM::TimeUnit const requestedTimeUnit, int *const error)
Provides the interface to a KIM API Model object for use by models within their MODEL_ROUTINE_NAME::C...
static int Destroy(KIM::ModelDestroy *const modelDestroy)
Numbering const zeroBased
The standard zeroBased numbering.
TimeUnit const unused
Indicates that a TimeUnit is not used.
static int ComputeArgumentsDestroy(KIM::ModelCompute const *const, KIM::ModelComputeArgumentsDestroy *const)
int SetRoutinePointer(ModelRoutineName const modelRoutineName, LanguageName const languageName, int const required, Function *const fptr)
Set the function pointer for the ModelRoutineName of interest.
LogVerbosity const error
The standard error verbosity.
static int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
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 ...
static int ComputeArgumentsCreate(KIM::ModelCompute const *const, KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)