diff --git a/AMSS_NCKU_source/surface_integral.C b/AMSS_NCKU_source/surface_integral.C index 243c865..1b92167 100644 --- a/AMSS_NCKU_source/surface_integral.C +++ b/AMSS_NCKU_source/surface_integral.C @@ -36,8 +36,15 @@ using namespace std; //| Constructor //|============================================================================ -surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry) -{ +surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry), + wave_cache_spinw(-1), + wave_cache_maxl(-1), + wave_cache_modes(0), + wave_theta_pos(0), + wave_theta_neg(0), + wave_phi_cos(0), + wave_phi_sin(0) +{ MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_size(MPI_COMM_WORLD, &cpusize); int N = 40; @@ -180,20 +187,80 @@ surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry) //|============================================================================ //| Destructor //|============================================================================ -surface_integral::~surface_integral() -{ - delete[] nx_g; - delete[] ny_g; - delete[] nz_g; - delete[] arcostheta; -#ifdef GaussInt - delete[] wtcostheta; -#endif -} -//|---------------------------------------------------------------- -// spin weighted spinw component of psi4, general routine -// l takes from spinw to maxl; m takes from -l to l -//|---------------------------------------------------------------- +surface_integral::~surface_integral() +{ + clear_wave_cache(); + delete[] nx_g; + delete[] ny_g; + delete[] nz_g; + delete[] arcostheta; +#ifdef GaussInt + delete[] wtcostheta; +#endif +} +void surface_integral::clear_wave_cache() +{ + delete[] wave_theta_pos; + delete[] wave_theta_neg; + delete[] wave_phi_cos; + delete[] wave_phi_sin; + wave_theta_pos = 0; + wave_theta_neg = 0; + wave_phi_cos = 0; + wave_phi_sin = 0; + wave_cache_spinw = -1; + wave_cache_maxl = -1; + wave_cache_modes = 0; +} +void surface_integral::build_wave_cache(int spinw, int maxl) +{ + if (wave_cache_spinw == spinw && wave_cache_maxl == maxl && wave_theta_pos && wave_theta_neg && wave_phi_cos && wave_phi_sin) + return; + + clear_wave_cache(); + + int modes = 0; + for (int pl = spinw; pl < maxl + 1; pl++) + for (int pm = -pl; pm < pl + 1; pm++) + modes++; + + wave_theta_pos = new double[modes * N_theta]; + wave_theta_neg = new double[modes * N_theta]; + wave_phi_cos = new double[modes * N_phi]; + wave_phi_sin = new double[modes * N_phi]; + + int countlm = 0; + for (int pl = spinw; pl < maxl + 1; pl++) + for (int pm = -pl; pm < pl + 1; pm++) + { + const double prefactor = sqrt((2 * pl + 1.0) / 4.0 / PI); + for (int i = 0; i < N_theta; i++) + { +#ifdef GaussInt + const double weight = wtcostheta[i]; +#else + const double weight = 1.0; +#endif + wave_theta_pos[countlm * N_theta + i] = prefactor * misc::Wigner_d_function(pl, pm, spinw, arcostheta[i]) * weight; + wave_theta_neg[countlm * N_theta + i] = prefactor * misc::Wigner_d_function(pl, pm, spinw, -arcostheta[i]) * weight; + } + for (int j = 0; j < N_phi; j++) + { + const double phase = pm * (j + 0.5) * dphi; + wave_phi_cos[countlm * N_phi + j] = cos(phase); + wave_phi_sin[countlm * N_phi + j] = sin(phase); + } + countlm++; + } + + wave_cache_spinw = spinw; + wave_cache_maxl = maxl; + wave_cache_modes = modes; +} +//|---------------------------------------------------------------- +// spin weighted spinw component of psi4, general routine +// l takes from spinw to maxl; m takes from -l to l +//|---------------------------------------------------------------- void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP, monitor *Monitor) // NN is the length of RP and IP @@ -254,100 +321,134 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * double costheta, thetap; double cosmphi, sinmphi; - int i, j; - int lpsy = 0; - if (Symmetry == 0) - lpsy = 1; - else if (Symmetry == 1) - lpsy = 2; - else if (Symmetry == 2) - lpsy = 8; - - double psi4RR, psi4II; - for (n = Nmin; n <= Nmax; n++) - { - // need round off always - i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 - j = n - i * N_phi; - - int countlm = 0; - for (int pl = spinw; pl < maxl + 1; pl++) - for (int pm = -pl; pm < pl + 1; pm++) - { - for (int lp = 0; lp < lpsy; lp++) - { - switch (lp) - { - case 0: //+++ (theta, phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = sin(pm * (j + 0.5) * dphi); - psi4RR = shellf[InList * n]; - psi4II = shellf[InList * n + 1]; - break; - case 1: //++- (pi-theta, phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; - break; - case 2: //+-+ (theta, 2*pi-phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = -sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; - psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; - break; - case 3: //+-- (pi-theta, 2*pi-phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = -sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; - break; - case 4: //-++ (theta, pi-phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 5: //-+- (pi-theta, pi-phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 6: //--+ (theta, pi+phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 7: //--- (pi-theta, pi+phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - } - - thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 -#ifdef GaussInt - // wtcostheta is even function respect costheta - RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; - IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; -#else - RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); - IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); -#endif - } - countlm++; // no sanity check for countlm and NN which should be noted in the input parameters - } - } + int i, j; + int lpsy = 0; + if (Symmetry == 0) + lpsy = 1; + else if (Symmetry == 1) + lpsy = 2; + else if (Symmetry == 2) + lpsy = 8; + + double psi4RR, psi4II; + if (Symmetry == 0 || Symmetry == 1) + { + build_wave_cache(spinw, maxl); + for (n = Nmin; n <= Nmax; n++) + { + i = int(n / N_phi); + j = n - i * N_phi; + + const double psi4RR0 = shellf[InList * n]; + const double psi4II0 = shellf[InList * n + 1]; + const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0; + const double psi4II1 = Ipsi4->SoA[2] * psi4II0; + for (int countlm = 0; countlm < wave_cache_modes; countlm++) + { + const int theta_idx = countlm * N_theta + i; + const int phi_idx = countlm * N_phi + j; + const double theta_pos = wave_theta_pos[theta_idx]; + const double cosmphi_here = wave_phi_cos[phi_idx]; + const double sinmphi_here = wave_phi_sin[phi_idx]; + + RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here); + IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here); + if (Symmetry == 1) + { + const double theta_neg = wave_theta_neg[theta_idx]; + RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here); + IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here); + } + } + } + } + else + { + for (n = Nmin; n <= Nmax; n++) + { + // need round off always + i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 + j = n - i * N_phi; + + int countlm = 0; + for (int pl = spinw; pl < maxl + 1; pl++) + for (int pm = -pl; pm < pl + 1; pm++) + { + for (int lp = 0; lp < lpsy; lp++) + { + switch (lp) + { + case 0: //+++ (theta, phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = sin(pm * (j + 0.5) * dphi); + psi4RR = shellf[InList * n]; + psi4II = shellf[InList * n + 1]; + break; + case 1: //++- (pi-theta, phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; + break; + case 2: //+-+ (theta, 2*pi-phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = -sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; + psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; + break; + case 3: //+-- (pi-theta, 2*pi-phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = -sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; + break; + case 4: //-++ (theta, pi-phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 5: //-+- (pi-theta, pi-phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 6: //--+ (theta, pi+phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 7: //--- (pi-theta, pi+phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + } + + thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 +#ifdef GaussInt + // wtcostheta is even function respect costheta + RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; + IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; +#else + RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); + IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); +#endif + } + countlm++; // no sanity check for countlm and NN which should be noted in the input parameters + } + } + } for (int ii = 0; ii < NN; ii++) { @@ -456,100 +557,134 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var * double costheta, thetap; double cosmphi, sinmphi; - int i, j; - int lpsy = 0; - if (Symmetry == 0) - lpsy = 1; - else if (Symmetry == 1) - lpsy = 2; - else if (Symmetry == 2) - lpsy = 8; - - double psi4RR, psi4II; - for (n = Nmin; n <= Nmax; n++) - { - // need round off always - i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 - j = n - i * N_phi; - - int countlm = 0; - for (int pl = spinw; pl < maxl + 1; pl++) - for (int pm = -pl; pm < pl + 1; pm++) - { - for (int lp = 0; lp < lpsy; lp++) - { - switch (lp) - { - case 0: //+++ (theta, phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = sin(pm * (j + 0.5) * dphi); - psi4RR = shellf[InList * n]; - psi4II = shellf[InList * n + 1]; - break; - case 1: //++- (pi-theta, phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; - break; - case 2: //+-+ (theta, 2*pi-phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = -sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; - psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; - break; - case 3: //+-- (pi-theta, 2*pi-phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (j + 0.5) * dphi); - sinmphi = -sin(pm * (j + 0.5) * dphi); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; - break; - case 4: //-++ (theta, pi-phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 5: //-+- (pi-theta, pi-phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 6: //--+ (theta, pi+phi) - costheta = arcostheta[i]; - cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - break; - case 7: //--- (pi-theta, pi+phi) - costheta = -arcostheta[i]; - cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); - sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); - psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; - psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; - } - - thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 -#ifdef GaussInt - // wtcostheta is even function respect costheta - RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; - IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; -#else - RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); - IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); -#endif - } - countlm++; // no sanity check for countlm and NN which should be noted in the input parameters - } - } + int i, j; + int lpsy = 0; + if (Symmetry == 0) + lpsy = 1; + else if (Symmetry == 1) + lpsy = 2; + else if (Symmetry == 2) + lpsy = 8; + + double psi4RR, psi4II; + if (Symmetry == 0 || Symmetry == 1) + { + build_wave_cache(spinw, maxl); + for (n = Nmin; n <= Nmax; n++) + { + i = int(n / N_phi); + j = n - i * N_phi; + + const double psi4RR0 = shellf[InList * n]; + const double psi4II0 = shellf[InList * n + 1]; + const double psi4RR1 = Rpsi4->SoA[2] * psi4RR0; + const double psi4II1 = Ipsi4->SoA[2] * psi4II0; + for (int countlm = 0; countlm < wave_cache_modes; countlm++) + { + const int theta_idx = countlm * N_theta + i; + const int phi_idx = countlm * N_phi + j; + const double theta_pos = wave_theta_pos[theta_idx]; + const double cosmphi_here = wave_phi_cos[phi_idx]; + const double sinmphi_here = wave_phi_sin[phi_idx]; + + RP_out[countlm] += theta_pos * (psi4RR0 * cosmphi_here + psi4II0 * sinmphi_here); + IP_out[countlm] += theta_pos * (psi4II0 * cosmphi_here - psi4RR0 * sinmphi_here); + if (Symmetry == 1) + { + const double theta_neg = wave_theta_neg[theta_idx]; + RP_out[countlm] += theta_neg * (psi4RR1 * cosmphi_here + psi4II1 * sinmphi_here); + IP_out[countlm] += theta_neg * (psi4II1 * cosmphi_here - psi4RR1 * sinmphi_here); + } + } + } + } + else + { + for (n = Nmin; n <= Nmax; n++) + { + // need round off always + i = int(n / N_phi); // int(1.723) = 1, int(-1.732) = -1 + j = n - i * N_phi; + + int countlm = 0; + for (int pl = spinw; pl < maxl + 1; pl++) + for (int pm = -pl; pm < pl + 1; pm++) + { + for (int lp = 0; lp < lpsy; lp++) + { + switch (lp) + { + case 0: //+++ (theta, phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = sin(pm * (j + 0.5) * dphi); + psi4RR = shellf[InList * n]; + psi4II = shellf[InList * n + 1]; + break; + case 1: //++- (pi-theta, phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[2] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * shellf[InList * n + 1]; + break; + case 2: //+-+ (theta, 2*pi-phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = -sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[1] * shellf[InList * n]; + psi4II = Ipsi4->SoA[1] * shellf[InList * n + 1]; + break; + case 3: //+-- (pi-theta, 2*pi-phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (j + 0.5) * dphi); + sinmphi = -sin(pm * (j + 0.5) * dphi); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * shellf[InList * n + 1]; + break; + case 4: //-++ (theta, pi-phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 5: //-+- (pi-theta, pi-phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (PI - (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI - (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 6: //--+ (theta, pi+phi) + costheta = arcostheta[i]; + cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + break; + case 7: //--- (pi-theta, pi+phi) + costheta = -arcostheta[i]; + cosmphi = cos(pm * (PI + (j + 0.5) * dphi)); + sinmphi = sin(pm * (PI + (j + 0.5) * dphi)); + psi4RR = Rpsi4->SoA[2] * Rpsi4->SoA[1] * Rpsi4->SoA[0] * shellf[InList * n]; + psi4II = Ipsi4->SoA[2] * Ipsi4->SoA[1] * Ipsi4->SoA[0] * shellf[InList * n + 1]; + } + + thetap = sqrt((2 * pl + 1.0) / 4.0 / PI) * misc::Wigner_d_function(pl, pm, spinw, costheta); // note the variation from -2 to 2 +#ifdef GaussInt + // wtcostheta is even function respect costheta + RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi) * wtcostheta[i]; + IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi) * wtcostheta[i]; +#else + RP_out[countlm] = RP_out[countlm] + thetap * (psi4RR * cosmphi + psi4II * sinmphi); + IP_out[countlm] = IP_out[countlm] + thetap * (psi4II * cosmphi - psi4RR * sinmphi); +#endif + } + countlm++; // no sanity check for countlm and NN which should be noted in the input parameters + } + } + } for (int ii = 0; ii < NN; ii++) { diff --git a/AMSS_NCKU_source/surface_integral.h b/AMSS_NCKU_source/surface_integral.h index 7a24b81..e1e05e8 100644 --- a/AMSS_NCKU_source/surface_integral.h +++ b/AMSS_NCKU_source/surface_integral.h @@ -27,19 +27,24 @@ using namespace std; class surface_integral { -private: - int Symmetry, factor; - int N_theta, N_phi; // Number of points in Theta & Phi directions - double dphi, dcostheta; - double *arcostheta, *wtcostheta; - int n_tot; // size of arrays - - double *nx_g, *ny_g, *nz_g; // global list of unit normals - int myrank, cpusize; - -public: - surface_integral(int iSymmetry); - ~surface_integral(); +private: + int Symmetry, factor; + int N_theta, N_phi; // Number of points in Theta & Phi directions + double dphi, dcostheta; + double *arcostheta, *wtcostheta; + int n_tot; // size of arrays + + double *nx_g, *ny_g, *nz_g; // global list of unit normals + int myrank, cpusize; + int wave_cache_spinw, wave_cache_maxl, wave_cache_modes; + double *wave_theta_pos, *wave_theta_neg; + double *wave_phi_cos, *wave_phi_sin; + void clear_wave_cache(); + void build_wave_cache(int spinw, int maxl); + +public: + surface_integral(int iSymmetry); + ~surface_integral(); void surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4, int spinw, int maxl, int NN, double *RP, double *IP,