* Used to calculate the partial power spectrum. * * The power spectrum is multiplied by an l-dependent prefactor * PI*sqrt(8.0/(2.0*l+1.0)); that comes from the normalization of the Wigner D * matrices. This prefactor is mentioned in the errata of the original SOAP * paper: On representing chemical environments, Phys. Rev. B 87, 184115 * (2013). Here the square root of the prefactor in the d
| 1746 | * used, so that after a possible dot-product the full prefactor is recovered. |
| 1747 | */ |
| 1748 | void getP(py::detail::unchecked_mutable_reference<double, 2> &Ps, double* Cs, int Nt, int lMax, int nMax, int Hs, double rCut2, int nFeatures, bool crossover, int nCoeffs) |
| 1749 | { |
| 1750 | // The current index in the final power spectrum array. |
| 1751 | int pIdx = 0; |
| 1752 | |
| 1753 | for (int i = 0; i < Hs; i++) { |
| 1754 | pIdx = 0; |
| 1755 | for (int Z1 = 0; Z1 < Nt; Z1++) { |
| 1756 | int Z2Limit = crossover ? Nt : Z1+1; |
| 1757 | for (int Z2 = Z1; Z2 < Z2Limit; Z2++) { |
| 1758 | // If the species are identical, then there is symmetry in the |
| 1759 | // radial basis and we only loop N2 from N1 to nMax |
| 1760 | if (Z1 == Z2) { |
| 1761 | for (int l = 0; l < lMax+1; l++) { |
| 1762 | for (int N1 = 0; N1 < nMax; N1++) { |
| 1763 | for (int N2 = N1; N2 < nMax; N2++) { |
| 1764 | double sum = 0; |
| 1765 | for (int m = 0; m < l+1; m++) { |
| 1766 | if (m == 0) { |
| 1767 | sum += Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1)] // m=0 |
| 1768 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1)]; // m=0 |
| 1769 | } else { |
| 1770 | sum += 2*(Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m] |
| 1771 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m] |
| 1772 | +Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m + 1] |
| 1773 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m + 1]); |
| 1774 | } |
| 1775 | } |
| 1776 | Ps(i, pIdx) = PI*sqrt(8.0/(2.0*l+1.0))*39.478417604*rCut2*sum; // Normalization and other constants |
| 1777 | ++pIdx; |
| 1778 | } |
| 1779 | } |
| 1780 | } |
| 1781 | // If the species are different, then there is no symmetry in the |
| 1782 | // radial basis and we have to loop over all pairwise combinations. |
| 1783 | } else { |
| 1784 | for (int l = 0; l < lMax+1; l++) { |
| 1785 | for (int N1 = 0; N1 < nMax; N1++) { |
| 1786 | for (int N2 = 0; N2 < nMax; N2++) { |
| 1787 | double sum = 0; |
| 1788 | for (int m = 0; m < l+1; m++) { |
| 1789 | if (m == 0) { |
| 1790 | sum += Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1)] // m=0 |
| 1791 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1)]; // m=0 |
| 1792 | } else { |
| 1793 | sum += 2*(Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m] |
| 1794 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m] |
| 1795 | +Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m + 1] |
| 1796 | *Cs[i*nCoeffs+2*Z2*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m + 1]); |
| 1797 | } |
| 1798 | } |
| 1799 | Ps(i, pIdx) = PI*sqrt(8.0/(2.0*l+1.0))*39.478417604*rCut2*sum; // Normalization and other constants |
| 1800 | ++pIdx; |
| 1801 | } |
| 1802 | } |
| 1803 | } |
| 1804 | } |
| 1805 | } |