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;
195 int* particleSpecies;
206 double* particleEnergy;
212 KIM_API_getm_compute(pkim, &ier, 5*3,
213 "energy", &comp_energy, 1,
214 "forces", &comp_force, 1,
215 "particleEnergy", &comp_particleEnergy, 1,
216 "process_dEdr", &comp_process_dEdr, 1,
217 "process_d2Edr2", &comp_process_d2Edr2, 1
219 if (KIM_STATUS_OK > ier) {
220 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_getm_compute", ier);
225 "cutoff", &cutoff, 1,
226 "numberOfParticles", &nParts, 1,
227 "particleSpecies", &particleSpecies,1,
228 "coordinates", &coords, 1,
229 "energy", &energy, comp_energy,
230 "forces", &force, comp_force,
231 "particleEnergy", &particleEnergy, comp_particleEnergy
233 if (KIM_STATUS_OK > ier) {
234 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_getm_data", ier);
238 cutsq = (*cutoff)*(*cutoff);
245 calc_phi(&epsilon, &
C, &Rzero, &dummy, cutoff, *cutoff, &shift);
251 ier = KIM_STATUS_FAIL;
252 for (i = 0; i < *nParts; ++i) {
253 if (
SPECCODE != particleSpecies[i]) {
254 KIM_API_report_error(__LINE__, __FILE__,
255 "Unexpected species detected", ier);
260 if (comp_particleEnergy) {
261 for (i = 0; i < *nParts; ++i) {
262 particleEnergy[i] = 0.0; } }
267 for (i = 0; i < *nParts; ++i) {
268 for (k = 0; k <
DIM; ++k) {
269 force[i*
DIM + k] = 0.0; } } }
274 for (i = 0; i < *nParts; ++i) {
276 ier = KIM_API_get_neigh(pkim, one, request, ¤tPart,
277 &numOfPartNeigh, &neighListOfCurrentPart,
279 if (KIM_STATUS_OK != ier) {
281 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_neigh", ier);
282 ier = KIM_STATUS_FAIL;
286 for (jj = 0; jj < numOfPartNeigh; ++ jj) {
287 j = neighListOfCurrentPart[jj];
291 for (k = 0; k <
DIM; ++k) {
292 Rij[k] = coords[j*
DIM + k] - coords[i*
DIM + k];
295 Rsqij += Rij[k]*Rij[k]; }
301 if (comp_process_d2Edr2) {
307 cutoff, R, &phi, &dphi, &d2phi);
311 d2Eidr = 0.5*d2phi; }
312 else if (comp_force || comp_process_dEdr) {
318 cutoff, R, &phi, &dphi);
331 if (comp_particleEnergy) {
332 particleEnergy[i] += 0.5*phi; }
334 *energy += 0.5*phi; }
337 if (comp_process_dEdr) {
338 ier = KIM_API_process_dEdr(km, &dEidr, &R, &pRij, &i, &j); }
341 if (comp_process_d2Edr2) {
342 R_pairs[0] = R_pairs[1] = R;
343 Rij_pairs[0][0] = Rij_pairs[1][0] = Rij[0];
344 Rij_pairs[0][1] = Rij_pairs[1][1] = Rij[1];
345 Rij_pairs[0][2] = Rij_pairs[1][2] = Rij[2];
346 i_pairs[0] = i_pairs[1] = i;
347 j_pairs[0] = j_pairs[1] = j;
349 ier = KIM_API_process_d2Edr2(km, &d2Eidr, &pR_pairs, &pRij_pairs,
350 &pi_pairs, &pj_pairs); }
354 for (k = 0; k <
DIM; ++k) {
355 force[i*
DIM + k] += dEidr*Rij[k]/R;
356 force[j*
DIM + k] -= dEidr*Rij[k]/R; } } }
368 intptr_t* pkim = *((intptr_t**) km);
375 ier = KIM_API_set_method(pkim,
"compute", 1, (func_ptr) &
compute);
376 if (KIM_STATUS_OK > ier) {
377 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_set_data", ier);
381 model_cutoff = (
double*) KIM_API_get_data(pkim,
"cutoff", &ier);
382 if (KIM_STATUS_OK > ier) {
383 KIM_API_report_error(__LINE__, __FILE__,
"KIM_API_get_data", ier);
387 return KIM_STATUS_OK; }
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
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)
real(c_double), parameter, public model_cutoff
static int compute(void *km)