47 #include "KIM_API_C.h" 48 #include "KIM_API_status.h" 56 #define EPSILON -0.0134783698072604 68 static void calc_phi(
double* epsilon,
72 double* cutoff,
double r,
double* phi);
78 double* cutoff,
double r,
double* phi,
double* dphi);
84 double* cutoff,
double r,
double* phi,
double* dphi,
92 double* cutoff,
double r,
double* phi) {
97 ep = exp(-(*
C)*(r-*Rzero));
104 *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift; }
113 double* cutoff,
double r,
double* phi,
double* dphi) {
118 ep = exp(-(*
C)*(r-*Rzero));
126 *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
127 *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 ); }
137 double* cutoff,
double r,
double* phi,
double* dphi,
143 ep = exp(-(*
C)*(r-*Rzero));
152 *phi = (*epsilon)*( -ep2 + 2.0*ep ) + *shift;
153 *dphi = 2.0*(*epsilon)*(*C)*( -ep + ep2 );
154 *d2phi = 2.0*(*epsilon)*(*C)*(*C)*(ep - 2.0*ep2); }
161 intptr_t* pkim = *((intptr_t**) km);
164 double *pR_pairs = &(R_pairs[0]);
172 double *pRij = &(Rij[0]);
173 double Rij_pairs[2][3];
174 double *pRij_pairs = &(Rij_pairs[0][0]);
178 int *pi_pairs = &(i_pairs[0]);
181 int *pj_pairs = &(j_pairs[0]);
185 int* neighListOfCurrentPart;
188 int comp_particleEnergy;
189 int comp_process_dEdr;
190 int comp_process_d2Edr2;
199 int* particleSpecies;
210 double* particleEnergy;
218 ier = KIM_API_get_NBC_method(pkim, &NBCstr);
219 if (KIM_STATUS_OK > ier) {
220 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_NBC_method", ier);
222 if ((!strcmp(
"NEIGH_PURE_H",NBCstr)) || (!strcmp(
"NEIGH_PURE_F",NBCstr))) {
224 else if (!strcmp(
"CLUSTER",NBCstr)) {
227 ier = KIM_STATUS_FAIL;
228 KIM_API_report_error(__LINE__, __FILE__,
"Unknown NBC method", ier);
236 if (KIM_API_is_half_neighbors(pkim, &ier)) {
247 IterOrLoca = KIM_API_get_neigh_mode(pkim, &ier);
248 if (KIM_STATUS_OK > ier) {
249 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_neigh_mode", ier);
251 if ((IterOrLoca != 1) && (IterOrLoca != 2)) {
252 ier = KIM_STATUS_FAIL;
253 KIM_API_report_error(__LINE__, __FILE__,
254 "Unsupported IterOrLoca mode", ier);
261 KIM_API_getm_compute(pkim, &ier, 5*3,
262 "energy", &comp_energy, 1,
263 "forces", &comp_force, 1,
264 "particleEnergy", &comp_particleEnergy, 1,
265 "process_dEdr", &comp_process_dEdr, 1,
266 "process_d2Edr2", &comp_process_d2Edr2, 1
268 if (KIM_STATUS_OK > ier) {
269 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_getm_compute", ier);
274 "cutoff", &cutoff, 1,
275 "numberOfParticles", &nParts, 1,
276 "particleSpecies", &particleSpecies,1,
277 "coordinates", &coords, 1,
278 "numberContributingParticles", &numContrib, (HalfOrFull==1),
279 "energy", &energy, comp_energy,
280 "forces", &force, comp_force,
281 "particleEnergy", &particleEnergy, comp_particleEnergy
283 if (KIM_STATUS_OK > ier) {
284 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_getm_data", ier);
288 cutsq = (*cutoff)*(*cutoff);
295 calc_phi(&epsilon, &
C, &Rzero, &dummy, cutoff, *cutoff, &shift);
299 if (HalfOrFull == 1) {
302 numberContrib = *numContrib; }
305 numberContrib = *nParts; } }
308 numberContrib = *nParts; }
312 ier = KIM_STATUS_FAIL;
313 for (i = 0; i < *nParts; ++i) {
314 if (
SPECCODE != particleSpecies[i]) {
315 KIM_API_report_error(__LINE__, __FILE__,
316 "Unexpected species detected", ier);
321 if (comp_particleEnergy) {
322 for (i = 0; i < *nParts; ++i) {
323 particleEnergy[i] = 0.0; } }
328 for (i = 0; i < *nParts; ++i) {
329 for (k = 0; k <
DIM; ++k) {
330 force[i*
DIM + k] = 0.0; } } }
335 neighListOfCurrentPart = (
int *) malloc((*nParts)*
sizeof(int)); }
339 if (1 == IterOrLoca) {
340 ier = KIM_API_get_neigh(pkim, zero, zero, ¤tPart, &numOfPartNeigh,
341 &neighListOfCurrentPart, &Rij_list);
343 if (KIM_STATUS_NEIGH_ITER_INIT_OK != ier) {
344 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_neigh", ier);
345 ier = KIM_STATUS_FAIL;
354 if (1 == IterOrLoca) {
356 ier = KIM_API_get_neigh(pkim, zero, one, ¤tPart, &numOfPartNeigh,
357 &neighListOfCurrentPart, &Rij_list);
358 if (KIM_STATUS_NEIGH_ITER_PAST_END == ier) {
361 if (KIM_STATUS_OK > ier) {
363 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_neigh", ier);
375 numOfPartNeigh = *nParts - (i + 1);
376 for (k = 0; k < numOfPartNeigh; ++k) {
377 neighListOfCurrentPart[k] = i + k + 1; }
378 ier = KIM_STATUS_OK; }
381 ier = KIM_API_get_neigh(pkim, one, request, ¤tPart,
382 &numOfPartNeigh, &neighListOfCurrentPart,
384 if (KIM_STATUS_OK != ier) {
386 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_neigh", ier);
387 ier = KIM_STATUS_FAIL;
391 for (jj = 0; jj < numOfPartNeigh; ++ jj) {
392 j = neighListOfCurrentPart[jj];
396 for (k = 0; k <
DIM; ++k) {
397 Rij[k] = coords[j*
DIM + k] - coords[i*
DIM + k];
400 Rsqij += Rij[k]*Rij[k]; }
406 if (comp_process_d2Edr2) {
412 cutoff, R, &phi, &dphi, &d2phi);
415 if ((1 == HalfOrFull) && (j < numberContrib)) {
422 d2Eidr = 0.5*d2phi; } }
423 else if (comp_force || comp_process_dEdr) {
429 cutoff, R, &phi, &dphi);
432 if ((1 == HalfOrFull) && (j < numberContrib)) {
437 dEidr = 0.5*dphi; } }
447 if (comp_particleEnergy) {
448 particleEnergy[i] += 0.5*phi;
450 if ((1 == HalfOrFull) && (j < numberContrib)) {
451 particleEnergy[j] += 0.5*phi; } }
453 if ((1 == HalfOrFull) && (j < numberContrib)) {
458 *energy += 0.5*phi; } }
461 if (comp_process_dEdr) {
462 ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
465 if (comp_process_d2Edr2) {
466 R_pairs[0] = R_pairs[1] = R;
467 Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
468 Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
469 Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
470 i_pairs[0] = i_pairs[1] = i;
471 j_pairs[0] = j_pairs[1] = j;
473 ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
474 &pi_pairs, &pj_pairs); }
478 for (k = 0; k <
DIM; ++k) {
479 force[i*
DIM + k] += dEidr*Rij[k]/R;
480 force[j*
DIM + k] -= dEidr*Rij[k]/R; } } }
486 free(neighListOfCurrentPart); }
496 intptr_t* pkim = *((intptr_t**) km);
503 ier = KIM_API_set_method(pkim,
"compute", 1, (func_ptr) &
compute);
504 if (KIM_STATUS_OK > ier) {
505 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_set_data", ier);
509 model_cutoff = (
double*) KIM_API_get_data(pkim,
"cutoff", &ier);
510 if (KIM_STATUS_OK > ier) {
511 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_data", ier);
515 return KIM_STATUS_OK; }
static int compute(void *km)
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
static void calc_phi_d2phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi, double *d2phi)
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
real(c_double), parameter, public model_cutoff