| 1513 | } |
| 1514 | } |
| 1515 | pair<int, int> getDeltas(double* dx, double* dy, double* dz, double* ri, double* rw, double rCut, double* oOri, double* oO4arri, double* minExp, double* pluExp, double eta, const py::array_t<double> &positions, const double ix, const double iy, const double iz, const vector<int> &indices, int rsize, int Ihpos, int Itype) |
| 1516 | { |
| 1517 | int iNeighbour = 0; |
| 1518 | int iCenter = 0; |
| 1519 | double ri2; |
| 1520 | double oOa = 1/eta; |
| 1521 | double Xi; double Yi; double Zi; |
| 1522 | int nNeighbours = indices.size(); |
| 1523 | double* oO4ari = (double*) malloc(sd*nNeighbours); |
| 1524 | |
| 1525 | auto pos = positions.unchecked<2>(); |
| 1526 | for (const int &i : indices) { |
| 1527 | Xi = pos(i, 0) - ix; |
| 1528 | Yi = pos(i, 1) - iy; |
| 1529 | Zi = pos(i, 2) - iz; |
| 1530 | ri2 = Xi*Xi + Yi*Yi + Zi*Zi; |
| 1531 | |
| 1532 | // When an atom is very close to the center (=approximately on top of |
| 1533 | // it), we do not add it to the calculations, as the numerical |
| 1534 | // integration cannot handle these cases. Instead, we gather the number |
| 1535 | // of such centered atoms and report them back for later correction. |
| 1536 | if (ri2<=1e-12) { |
| 1537 | iCenter++; |
| 1538 | } else { |
| 1539 | ri[iNeighbour] = sqrt(ri2); |
| 1540 | dx[iNeighbour] = Xi; |
| 1541 | dy[iNeighbour] = Yi; |
| 1542 | dz[iNeighbour] = Zi; |
| 1543 | oOri[iNeighbour] = 1/ri[iNeighbour]; |
| 1544 | oO4ari[iNeighbour] = 0.25*oOa*oOri[iNeighbour]; |
| 1545 | iNeighbour++; |
| 1546 | } |
| 1547 | } |
| 1548 | |
| 1549 | // If there is at least one atom at the center, we add a zero element to |
| 1550 | // the end of ris so that the weights can be calculated. This way they do |
| 1551 | // not interfere with the calculations for non-centered atoms. |
| 1552 | if (iCenter > 0) { |
| 1553 | ri[iNeighbour] = 0; |
| 1554 | } |
| 1555 | |
| 1556 | double* oOr = getoOr(rw, rsize); |
| 1557 | for (int i = 0; i < iNeighbour; i++) { |
| 1558 | for (int w = 0; w < rsize; w++) { |
| 1559 | oO4arri[rsize*i + w] = oO4ari[i]*oOr[w]; |
| 1560 | } |
| 1561 | } |
| 1562 | expMs(minExp, eta, rw, ri, iNeighbour, rsize); |
| 1563 | expPs(pluExp, eta, rw, ri, iNeighbour, rsize); |
| 1564 | |
| 1565 | free(oO4ari); |
| 1566 | return make_pair(iNeighbour, iCenter); |
| 1567 | } |
| 1568 | double* getFlir(double* oO4arri,double* ri, double* minExp, double* pluExp, int icount, int rsize, int lMax) |
| 1569 | { |
| 1570 | double* Flir = (double*) malloc(sd*(lMax+1)*icount*rsize); |
no test coverage detected