KIM API V2
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, Regents of the University of Minnesota.
23 // All rights reserved.
24 //
25 // Contributors:
26 // Ryan S. Elliott
27 //
28 
29 #include <math.h>
30 #include "KIM_ModelHeaders.hpp"
31 
32 #define DIMENSION 3
33 
34 
35 namespace
36 {
37 
39 {
40  public:
41  //****************************************************************************
43  KIM::ModelCreate * const modelCreate,
44  KIM::LengthUnit const requestedLengthUnit,
45  KIM::EnergyUnit const requestedEnergyUnit,
46  KIM::ChargeUnit const requestedChargeUnit,
47  KIM::TemperatureUnit const requestedTemperatureUnit,
48  KIM::TimeUnit const requestedTimeUnit,
49  int * const error) :
50  epsilon_(0.0104),
51  sigma_(3.4000),
52  influenceDistance_(8.1500),
53  cutoff_(influenceDistance_),
54  cutoffSq_(cutoff_*cutoff_),
55  paddingNeighborHints_(1),
56  halfListHints_(1)
57  {
58  *error = ConvertUnits(modelCreate,
59  requestedLengthUnit,
60  requestedEnergyUnit,
61  requestedChargeUnit,
62  requestedTemperatureUnit,
63  requestedTimeUnit);
64  if (*error) return;
65 
67 
68  modelCreate->SetInfluenceDistancePointer(&influenceDistance_);
69  modelCreate->SetNeighborListPointers(1,
70  &cutoff_,
71  &paddingNeighborHints_,
72  &halfListHints_);
73 
74  modelCreate->SetSpeciesCode(KIM::SPECIES_NAME::Ar, 0);
75 
76  *error = *error || modelCreate->SetDestroyPointer(
78  reinterpret_cast<KIM::func *>(&(LennardJones_Ar::Destroy)));
79  *error = *error || modelCreate->SetRefreshPointer(
81  reinterpret_cast<KIM::func *>(&(LennardJones_Ar::Refresh)));
82  *error = *error || modelCreate->SetComputePointer(
84  reinterpret_cast<KIM::func *>(&(LennardJones_Ar::Compute)));
85  *error = *error || modelCreate->SetComputeArgumentsCreatePointer(
87  reinterpret_cast<KIM::func *>
88  (&(LennardJones_Ar::ComputeArgumentsCreate)));
89  *error = *error || modelCreate->SetComputeArgumentsDestroyPointer(
91  reinterpret_cast<KIM::func *>
92  (&(LennardJones_Ar::ComputeArgumentsDestroy)));
93  if (*error) return;
94 
95  // everything is good
96  *error = false;
97  return;
98  };
99 
100  //****************************************************************************
102 
103  //****************************************************************************
104  // no need to make these "extern" since KIM will only access them
105  // via function pointers. "static" is required so that there is not
106  // an implicit this pointer added to the prototype by the C++ compiler
107  static int Destroy(KIM::ModelDestroy * const modelDestroy)
108  {
109  LennardJones_Ar * model;
110  modelDestroy->GetModelBufferPointer(reinterpret_cast<void **>(&model));
111 
112  if (model != NULL)
113  {
114  // delete object itself
115  delete model;
116  }
117 
118  // everything is good
119  return false;
120  }
121 
122  //****************************************************************************
123  static int Refresh(KIM::ModelRefresh * const modelRefresh)
124  {
125  LennardJones_Ar * model;
126  modelRefresh->GetModelBufferPointer(reinterpret_cast<void **>(&model));
127 
128  // nothing to do
129 
130  modelRefresh->SetInfluenceDistancePointer(&(model->influenceDistance_));
131  modelRefresh->SetNeighborListPointers(1,
132  &(model->cutoff_),
133  &(model->paddingNeighborHints_),
134  &(model->halfListHints_));
135 
136  // everything is good
137  return false;
138  };
139 
140  //****************************************************************************
142  static int Compute(
143  KIM::ModelCompute const * const modelCompute,
144  KIM::ModelComputeArguments const * const modelComputeArguments)
145  {
146  int const * numberOfParticlesPointer;
147  int const * particleSpeciesCodes;
148  int const * particleContributing;
149  double const * coordinates;
150  double * partialEnergy;
151  double * partialForces;
152 
153  LennardJones_Ar * lj;
154  modelCompute->GetModelBufferPointer(reinterpret_cast<void **>(&lj));
155  double const epsilon = lj->epsilon_;
156  double const sigma = lj->sigma_;
157  double const cutoffSq = lj->cutoffSq_;
158 
159  int error =
160  modelComputeArguments->GetArgumentPointer(
162  &numberOfParticlesPointer)
163  || modelComputeArguments->GetArgumentPointer(
166  || modelComputeArguments->GetArgumentPointer(
169  || modelComputeArguments->GetArgumentPointer(
171  (double const ** const) &coordinates)
172  || modelComputeArguments->GetArgumentPointer(
174  &partialEnergy)
175  || modelComputeArguments->GetArgumentPointer(
177  (double const ** const) &partialForces);
178  if (error)
179  {
180  LOG_ERROR("Unable to get argument pointers");
181  return error;
182  }
183 
184  int const numberOfParticles = *numberOfParticlesPointer;
185 
186  // initialize energy and forces
187  *partialEnergy = 0.0;
188  int const extent = numberOfParticles*DIMENSION;
189  for (int i = 0; i < extent; ++i)
190  {
191  partialForces[i] = 0.0;
192  }
193 
194  int jContributing;
195  int i, j, jj, numberOfNeighbors;
196  int const * neighbors;
197  double phi;
198  double xcoord, ycoord, zcoord;
199  double xrij, yrij, zrij;
200  double rij2;
201  double r2inv, r6inv, dphiByR, dEidrByR;
202  double xdf, ydf, zdf;
203  double const fortyEight = 48.0 * epsilon * pow(sigma, 12.0);
204  double const twentyFour = 24.0 * epsilon * pow(sigma, 6.0);
205  double const four12 = 4.0 * epsilon * pow(sigma, 12.0);
206  double const four6 = 4.0 * epsilon * pow(sigma, 6.0);
207  for (i = 0; i < numberOfParticles; ++i)
208  {
209  if (particleContributing[i])
210  {
211  xcoord = coordinates[i*DIMENSION + 0];
212  ycoord = coordinates[i*DIMENSION + 1];
213  zcoord = coordinates[i*DIMENSION + 2];
214 
215  modelComputeArguments->GetNeighborList(
216  0, i, &numberOfNeighbors, &neighbors);
217 
218  for (jj = 0; jj < numberOfNeighbors; ++jj)
219  {
220  j = neighbors[jj];
221  jContributing = particleContributing[j];
222 
223  if (i < j)
224  {
225  xrij = coordinates[j*DIMENSION + 0] - xcoord;
226  yrij = coordinates[j*DIMENSION + 1] - ycoord;
227  zrij = coordinates[j*DIMENSION + 2] - zcoord;
228 
229  rij2 = xrij*xrij + yrij*yrij + zrij*zrij;
230 
231  if (rij2 < cutoffSq)
232  {
233  r2inv = 1.0/rij2;
234  r6inv = r2inv*r2inv*r2inv;
235  phi = 0.5 * r6inv * (four12*r6inv - four6);
236  dphiByR = r6inv * (twentyFour - fortyEight*r6inv) * r2inv;
237 
238  *partialEnergy += phi;
239  if (jContributing)
240  {
241  *partialEnergy += phi;
242  dEidrByR = dphiByR;
243  }
244  else
245  {
246  dEidrByR = 0.5*dphiByR;
247  }
248 
249  xdf = dEidrByR * xrij;
250  ydf = dEidrByR * yrij;
251  zdf = dEidrByR * zrij;
252  partialForces[i*DIMENSION + 0] += xdf;
253  partialForces[i*DIMENSION + 1] += ydf;
254  partialForces[i*DIMENSION + 2] += zdf;
255  partialForces[j*DIMENSION + 0] -= xdf;
256  partialForces[j*DIMENSION + 1] -= ydf;
257  partialForces[j*DIMENSION + 2] -= zdf;
258  } // if (rij2 < cutoffSq_)
259  } // if (i < j)
260  } // for jj
261  } // if (particleContributing[i])
262  } // for i
263 
264  // everything is good
265  return false;
266  };
267 
268  //****************************************************************************
270  KIM::ModelCompute const * const modelCompute,
271  KIM::ModelComputeArgumentsCreate * const modelComputeArgumentsCreate)
272  {
273  // register arguments
274  int error =
275  modelComputeArgumentsCreate->SetArgumentSupportStatus(
278  || modelComputeArgumentsCreate->SetArgumentSupportStatus(
281 
282  // register callbacks
283  //
284  // none
285 
286  return error;
287  }
288  //****************************************************************************
290  KIM::ModelCompute const * const modelCompute,
291  KIM::ModelComputeArgumentsDestroy * const modelComputeArgumentsDestroy)
292  {
293  // noting to do
294 
295  // everything is good
296  return false;
297  };
298 
299  private:
300  //****************************************************************************
301  // Member variables
302  double epsilon_;
303  double sigma_;
304  double influenceDistance_;
305  double cutoff_;
306  double cutoffSq_;
307  int const paddingNeighborHints_;
308  int const halfListHints_;
309 
310  //****************************************************************************
312  int ConvertUnits(
313  KIM::ModelCreate * const modelCreate,
314  KIM::LengthUnit const requestedLengthUnit,
315  KIM::EnergyUnit const requestedEnergyUnit,
316  KIM::ChargeUnit const requestedChargeUnit,
317  KIM::TemperatureUnit const requestedTemperatureUnit,
318  KIM::TimeUnit const requestedTimeUnit)
319  {
320  int ier;
321 
322  // define default base units
328 
329  // changing units of cutoffs and sigmas
330  double convertLength = 1.0;
331  ier = modelCreate->ConvertUnit(
332  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
333  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
334  requestedTemperatureUnit, requestedTimeUnit,
335  1.0, 0.0, 0.0, 0.0, 0.0,
336  &convertLength);
337  if (ier)
338  {
339  LOG_ERROR("Unable to convert length unit");
340  return ier;
341  }
342  influenceDistance_ *= convertLength; // convert to active units
343  cutoff_ = influenceDistance_;
344  cutoffSq_ = cutoff_*cutoff_;
345  sigma_ *= convertLength; // convert to active units
346 
347  // changing units of epsilons
348  double convertEnergy = 1.0;
349  ier = modelCreate->ConvertUnit(
350  fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
351  requestedLengthUnit, requestedEnergyUnit, requestedChargeUnit,
352  requestedTemperatureUnit, requestedTimeUnit,
353  0.0, 1.0, 0.0, 0.0, 0.0,
354  &convertEnergy);
355  if (ier)
356  {
357  LOG_ERROR("Unable to convert energy unit");
358  return ier;
359  }
360  epsilon_ *= convertEnergy; // convert to active units
361 
362  // register units
363  ier = modelCreate->SetUnits(
364  requestedLengthUnit,
365  requestedEnergyUnit,
369  if (ier)
370  {
371  LOG_ERROR("Unable to set units to requested values");
372  return ier;
373  }
374 
375  // everything is good
376  ier = false;
377  return ier;
378  }
379 };
380 
381 } // namespace
382 
383 extern "C"
384 {
385 //******************************************************************************
387  KIM::ModelCreate * const modelCreate,
388  KIM::LengthUnit const requestedLengthUnit,
389  KIM::EnergyUnit const requestedEnergyUnit,
390  KIM::ChargeUnit const requestedChargeUnit,
391  KIM::TemperatureUnit const requestedTemperatureUnit,
392  KIM::TimeUnit const requestedTimeUnit)
393 {
394  int error;
395 
396  LennardJones_Ar * const model = new LennardJones_Ar(modelCreate,
397  requestedLengthUnit,
398  requestedEnergyUnit,
399  requestedChargeUnit,
400  requestedTemperatureUnit,
401  requestedTimeUnit,
402  &error);
403  if (error)
404  {
405  // constructor already reported the error
406  delete model;
407  return error;
408  }
409 
410  // register pointer to LennardJones_Ar object in moedelCreate object
411  modelCreate->SetModelBufferPointer(reinterpret_cast<void *>(model));
412 
413  // everything is good
414  return false;
415 }
416 } // extern "C"
TemperatureUnit const unused
SpeciesName const Ar
int SetSpeciesCode(SpeciesName const speciesName, int const code)
int SetModelNumbering(Numbering const numbering)
ComputeArgumentName const coordinates
int SetComputePointer(LanguageName const languageName, func *const fptr)
static int ComputeArgumentsDestroy(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)
int SetUnits(LengthUnit const lengthUnit, EnergyUnit const energyUnit, ChargeUnit const chargeUnit, TemperatureUnit const temperatureUnit, TimeUnit const timeUnit)
ChargeUnit const unused
void GetModelBufferPointer(void **const ptr) const
void SetInfluenceDistancePointer(double const *const influenceDistance)
static int Refresh(KIM::ModelRefresh *const modelRefresh)
void SetModelBufferPointer(void *const ptr)
SupportStatus const required
int GetArgumentPointer(ComputeArgumentName const computeArgumentName, int const **const ptr) const
ComputeArgumentName const particleContributing
ComputeArgumentName const partialEnergy
int SetComputeArgumentsDestroyPointer(LanguageName const languageName, func *const fptr)
void SetNeighborListPointers(int const numberOfNeighborLists, double const *const cutoffs, int const *const paddingNeighborHints, int const *const halfListHints)
static int ComputeArgumentsCreate(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
#define DIMENSION
#define LOG_ERROR(message)
int SetComputeArgumentsCreatePointer(LanguageName const languageName, func *const fptr)
void SetNeighborListPointers(int const numberOfNeighborLists, double const *const cutoffs, int const *const paddingNeighborHints, int const *const halfListHints)
LanguageName const cpp
int SetRefreshPointer(LanguageName const languageName, func *const fptr)
LengthUnit const A
ComputeArgumentName const partialForces
void GetModelBufferPointer(void **const ptr) const
EnergyUnit const eV
ComputeArgumentName const particleSpeciesCodes
int SetArgumentSupportStatus(ComputeArgumentName const clomputeArgumentName, SupportStatus const supportStatus)
ComputeArgumentName const numberOfParticles
void SetInfluenceDistancePointer(double const *const influenceDistance)
int GetNeighborList(int const neighborListIndex, int const particleNumber, int *const numberOfNeighbors, int const **const neighborsOfParticle) const
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)
void GetModelBufferPointer(void **const ptr) const
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)
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) const
static int Destroy(KIM::ModelDestroy *const modelDestroy)
Numbering const zeroBased
TimeUnit const unused
LogVerbosity const error
static int Compute(KIM::ModelCompute const *const modelCompute, KIM::ModelComputeArguments const *const modelComputeArguments)
int SetDestroyPointer(LanguageName const languageName, func *const fptr)