1 /* 2 * Copyright (c) 2020 Robert Shaw 3 * This file is a part of Libecpint. 4 * 5 * Permission is hereby granted, free of charge, to any person obtaining 6 * a copy of this software and associated documentation files (the 7 * "Software"), to deal in the Software without restriction, including 8 * without limitation the rights to use, copy, modify, merge, publish, 9 * distribute, sublicense, and/or sell copies of the Software, and to 10 * permit persons to whom the Software is furnished to do so, subject to 11 * the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be 14 * included in all copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 17 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 18 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 19 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 20 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 21 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 22 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 23 */ 24 25 #include "mathutil.hpp" 26 #include <iostream> 27 28 namespace libecpint { 29 30 double FAC[MAX_FAC]; 31 double DFAC[MAX_DFAC]; 32 pow_m2(const double z)33 double pow_m2(const double z) { return 1.0/(z*z); } pow_m1(const double z)34 double pow_m1(const double z) { return 1.0/z; } pow_0(const double z)35 double pow_0(const double z) { return 1.0; } pow_1(const double z)36 double pow_1(const double z) { return z; } pow_2(const double z)37 double pow_2(const double z) { return z*z; } pow_3(const double z)38 double pow_3(const double z) { return z*z*z; } pow_4(const double z)39 double pow_4(const double z) { double z2 = z*z; return z2*z2; } pow_5(const double z)40 double pow_5(const double z) { double z2 = z*z; double z3 = z2*z; return z2*z3; }; pow_6(const double z)41 double pow_6(const double z) { double z2 = z*z; return z2*z2*z2; } pow_7(const double z)42 double pow_7(const double z) { double z2 = z*z; double z3 = z*z2; return z3*z2*z2; } pow_8(const double z)43 double pow_8(const double z) { double z2 = z*z; double z4 = z2*z2; return z4*z4; } pow_9(const double z)44 double pow_9(const double z) { double z3 = z*z*z; return z3*z3*z3; } pow_10(const double z)45 double pow_10(const double z) { double z2 = z*z; double z3 = z*z2; double z5= z2*z3; return z5*z5; } pow_11(const double z)46 double pow_11(const double z) { double z2 = z*z; double z3 = z*z2; return z3*z3*z3*z2; } pow_12(const double z)47 double pow_12(const double z) { double z3 = z*z*z; double z6 = z3*z3; return z6*z6; } pow_13(const double z)48 double pow_13(const double z) { double z3 = z*z*z; double z6 = z3*z3; return z6*z6*z; } pow_14(const double z)49 double pow_14(const double z) { double z2 = z*z; double z3 = z*z2; double z7 = z2*z2*z3; return z7*z7; } pow_15(const double z)50 double pow_15(const double z) { double z2 = z*z; double z3 = z*z2; double z5 = z2*z3; return z5*z5*z5; } pow_16(const double z)51 double pow_16(const double z) { double z2 = z*z; double z4 = z2*z2; double z8 = z4*z4; return z8*z8; } pow_17(const double z)52 double pow_17(const double z) { double z2 = z*z; double z4 = z2*z2; double z8 = z4*z4; return z8*z8*z;} pow_18(const double z)53 double pow_18(const double z) { double z3 = z*z*z; double z9 = z3*z3*z3; return z9*z9; } pow_19(const double z)54 double pow_19(const double z) { double z3 = z*z*z; double z9 = z3*z3*z3; return z9*z9*z; } pow_20(const double z)55 double pow_20(const double z) { double z2 = z*z; double z4 = z2*z2; double z8 = z4*z4; return z8*z8*z4; } 56 initFactorials()57 void initFactorials() { 58 #ifndef FAC_INIT 59 #define FAC_INIT 60 FAC[0] = 1.0; 61 DFAC[0] = 1.0; 62 DFAC[1] = 1.0; 63 64 for (int i = 1; i < MAX_FAC; i++) FAC[i] = double(i) * FAC[i-1]; 65 for (int i = 2; i < MAX_DFAC; i++) DFAC[i] = double(i) * DFAC[i-2]; 66 #endif 67 } 68 69 // Compute all the real spherical harmonics Slm(theta, phi) for l,m up to lmax 70 // x = std::cos (theta) realSphericalHarmonics(const int lmax,const double x,const double phi)71 TwoIndex<double> realSphericalHarmonics(const int lmax, const double x, const double phi){ 72 TwoIndex<double> rshValues(lmax+1, 2*lmax+1, 0.0); 73 74 if (lmax > 0) { 75 // First calculate the associated Legendre polynomials, Plm(std::cos theta), ustd::sing the recursion relation 76 // (l-m)Plm = x(2l - 1)P{l-1}m - (l+m-1)P{l-2}m 77 // along with the zeroth order term 78 // Pmm = (-1)^m (2m-1)!!(1-x^2)^{m/2} 79 double x2 = x * x; 80 double Plm[lmax+1][lmax+1]; 81 // First get all Pmm terms 82 Plm[0][0] = 1.0; 83 double sox2 = std::sqrt(std::max(0.0, 1.0 - x2)); 84 double ox2m = 1.0; 85 for (int m = 1; m <= lmax; m++) { 86 ox2m *= -sox2; 87 Plm[m][m] = ox2m * DFAC[2*m-1]; 88 } 89 90 // Then increment l for each m 91 Plm[1][0] = x; 92 Plm[0][1] = 0.0; 93 for (int l = 2; l <= lmax; l++) { 94 ox2m = x * (2*l - 1); 95 for (int m = 0; m < l; m++) { 96 Plm[l][m] = ox2m * Plm[l-1][m] - (l + m - 1)*Plm[l-2][m]; 97 Plm[l][m] /= ((double) (l -m)); 98 } 99 Plm[l-1][l] = 0.0; 100 } 101 102 // Now we compute the spherical harmonics via 103 // Slm(theta, phi) = Clm * Plm(std::cos(theta)) * std::cos(m * phi), m > 0 104 // Sl{-m}(theta, phi) = Clm * Plm(std::cos(theta)) * std::sin(m * phi) 105 // Sl0(theta, phi) = std::sqrt(2) * Cl0 * Pl0(std::cos(theta)) 106 // where Clm^2 = (2l + 1)*(l - m)! / (8*pi * (l+m)!) 107 double osq4pi = 1.0 / std::sqrt(4.0 * M_PI); 108 int sign; 109 for (int l = 0; l <= lmax; l++) { 110 rshValues(l, l) = osq4pi * std::sqrt(2.0 * l + 1.0) * Plm[l][0]; 111 sign = -1; 112 for (int m = 1; m <= l; m++) { 113 ox2m = (2.0 * l + 1.0) * FAC[l-m] / FAC[l+m]; 114 ox2m = sign * osq4pi * std::sqrt(2.0 * ox2m) * Plm[l][m]; 115 rshValues(l, l+m) = ox2m * std::cos(m * phi); 116 rshValues(l, l-m) = ox2m * std::sin(m * phi); 117 sign *= -1; 118 } 119 } 120 121 } else { 122 rshValues(0, 0) = 1.0 / std::sqrt(4.0 * M_PI); 123 } 124 125 return rshValues; 126 } 127 frobenius_norm(const TwoIndex<double> & mat)128 double frobenius_norm(const TwoIndex<double>& mat) { 129 return std::sqrt(std::inner_product(mat.data.begin(), mat.data.end(), mat.data.begin(), 0.0)); 130 } 131 132 } 133