Function getPWithCompression

dscribe/ext/soapGeneral.cpp:1827–1859

* Used to calculate the partial power spectrum with feature compression * as described in Darby et al. Currently only for the mu = 1, nu = 1 * compression described in Darby et al. (by far the best performing * alternative that they evaluated). It could be adapted however to * provide other compression alternatives. This compression scheme * scales linearly with the number of elements and is

 * used, so that after a possible dot-product the full prefactor is recovered.
 */
void getPWithCompression(py::detail::unchecked_mutable_reference<double, 2> &Ps, double* Cs,
    double* CsSummed, int Nt, int lMax, int nMax, int Hs,
    double rCut2, int nFeatures, int nCoeffs, int nCompressionCoeffs)
{
   // The current index in the final power spectrum array.
    int pIdx = 0;

    for (int i = 0; i < Hs; i++) {
        pIdx = 0;
        for (int Z1 = 0; Z1 < Nt; Z1++) {
            for (int l = 0; l < lMax+1; l++) {
                for (int N1 = 0; N1 < nMax; N1++) {
                    for (int N2 = 0; N2 < nMax; N2++) {
                        double sum = 0;
                        for (int m = 0; m < l+1; m++) {
                            if (m == 0) {
                                sum += Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1)] // m=0
                                            *CsSummed[i*nCompressionCoeffs + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1)]; // m=0
                            } else {
                                sum += 2*(Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m]
                                                *CsSummed[i*nCompressionCoeffs + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m]
                                                +Cs[i*nCoeffs+2*Z1*(lMax+1)*(lMax+1)*nMax + 2*(lMax+1)*(lMax+1)*N1 + l*2*(lMax+1) + 2*m + 1]
                                                *CsSummed[i*nCompressionCoeffs + 2*(lMax+1)*(lMax+1)*N2 + l*2*(lMax+1) + 2*m + 1]);
                            }
                        }
                        Ps(i, pIdx) = PI*sqrt(8.0/(2.0*l+1.0))*39.478417604*rCut2*sum;  // Normalization and other constants
                        ++pIdx;
                    }
                }
            }
        }
    }
}