================================================================================================ * Used to calculate the partial power spectrum. */
| 2238 | * Used to calculate the partial power spectrum. |
| 2239 | */ |
| 2240 | void getPD( |
| 2241 | py::detail::unchecked_mutable_reference<double, 2> &descriptor_mu, |
| 2242 | py::detail::unchecked_reference<double, 4> &Cnnd_u, |
| 2243 | int Ns, |
| 2244 | int Ts, |
| 2245 | int nCenters, |
| 2246 | int lMax, |
| 2247 | bool crossover |
| 2248 | ) { |
| 2249 | |
| 2250 | // The power spectrum is multiplied by an l-dependent prefactor that comes |
| 2251 | // from the normalization of the Wigner D matrices. This prefactor is |
| 2252 | // mentioned in the arrata of the original SOAP paper: On representing |
| 2253 | // chemical environments, Phys. Rev. B 87, 184115 (2013). Here the square |
| 2254 | // root of the prefactor in the dot-product kernel is used, so that after a |
| 2255 | // possible dot-product the full prefactor is recovered. |
| 2256 | for(int i = 0; i < nCenters; i++){ |
| 2257 | int shiftAll = 0; |
| 2258 | for(int j = 0; j < Ts; j++){ |
| 2259 | int jdLimit = crossover ? Ts : j+1; |
| 2260 | for(int jd = j; jd < jdLimit; jd++){ |
| 2261 | for(int m=0; m <= lMax; m++){ |
| 2262 | double prel; |
| 2263 | if(m > 1){prel = PI*sqrt(8.0/(2.0*m+1.0))*PI3;} |
| 2264 | else{prel = PI*sqrt(8.0/(2.0*m+1.0));} |
| 2265 | if(j==jd){ |
| 2266 | for(int k = 0; k < Ns; k++){ |
| 2267 | for(int kd = k; kd < Ns; kd++){ |
| 2268 | double buffDouble = 0; |
| 2269 | for(int buffShift = m*m; buffShift < (m+1)*(m+1); buffShift++){ |
| 2270 | buffDouble += Cnnd_u(i,j,k,buffShift) * Cnnd_u(i,jd,kd,buffShift); |
| 2271 | } |
| 2272 | descriptor_mu(i, shiftAll) = prel*buffDouble; |
| 2273 | shiftAll++; |
| 2274 | } |
| 2275 | } |
| 2276 | } else { |
| 2277 | for(int k = 0; k < Ns; k++){ |
| 2278 | for(int kd = 0; kd < Ns; kd++){ |
| 2279 | double buffDouble = 0; |
| 2280 | for(int buffShift = m*m; buffShift < (m+1)*(m+1); buffShift++){ |
| 2281 | buffDouble += Cnnd_u(i,j,k,buffShift) * Cnnd_u(i,jd,kd,buffShift); |
| 2282 | } |
| 2283 | descriptor_mu(i, shiftAll) = prel*buffDouble; |
| 2284 | shiftAll++; |
| 2285 | } |
| 2286 | } |
| 2287 | } |
| 2288 | } //end ifelse |
| 2289 | } |
| 2290 | } |
| 2291 | } |
| 2292 | } |
| 2293 | |
| 2294 | |
| 2295 | //================================================================================================ |