diff --git a/source/module_base/mathzone_add1.cpp b/source/module_base/mathzone_add1.cpp index 459f0205b1..f22667c104 100644 --- a/source/module_base/mathzone_add1.cpp +++ b/source/module_base/mathzone_add1.cpp @@ -28,7 +28,6 @@ typedef fftw_complex FFTW_COMPLEX; namespace ModuleBase { -bool Mathzone_Add1::flag_jlx_expand_coef = false; double** Mathzone_Add1::c_ln_c = nullptr; double** Mathzone_Add1::c_ln_s = nullptr; @@ -39,605 +38,6 @@ Mathzone_Add1::~Mathzone_Add1() {} -/********************************************** - * coefficients to expand jlx using - * cos and sin, from SIESTA - * *******************************************/ -void Mathzone_Add1::expand_coef_jlx() -{ - ModuleBase::timer::tick("Mathzone_Add1","expand_coef_jlx"); - - int ir, il, in; - const int L = 20; - - //allocation - delete[] c_ln_c; - delete[] c_ln_s; - - c_ln_c = new double*[L+1]; - c_ln_s = new double*[L+1]; - - for(ir = 0; ir < L+1; ir++) - { - c_ln_c[ir] = new double[ir+1]; - c_ln_s[ir] = new double[ir+1]; - } - - //calculate initial value - c_ln_c[0][0] = 0.0; - c_ln_s[0][0] = 1.0; - - c_ln_c[1][0] = 0.0; - c_ln_c[1][1] = -1.0; - c_ln_s[1][0] = 1.0; - c_ln_s[1][1] = 0.0; - - //recursive equation - for(il = 2; il < Mathzone_Add1::sph_lmax+1; il++) - { - for(in = 0; in < il + 1; in++) - { - if(in >= 2) - { - if(in == il) - { - c_ln_c[il][in] = -c_ln_c[il-2][in-2]; - c_ln_s[il][in] = -c_ln_s[il-2][in-2]; - } - else - { - c_ln_c[il][in] = (2*il-1)*c_ln_c[il-1][in] - c_ln_c[il-2][in-2]; - c_ln_s[il][in] = (2*il-1)*c_ln_s[il-1][in] - c_ln_s[il-2][in-2]; - } - } - else - { - //in = 0 or 1, but here il = 2, so in != il - c_ln_c[il][in] = (2*il-1)*c_ln_c[il-1][in]; - c_ln_s[il][in] = (2*il-1)*c_ln_s[il-1][in]; - } - } - } - - ModuleBase::timer::tick("Mathzone_Add1","expand_coef_jlx"); - return; -} - -void Mathzone_Add1::Spherical_Bessel -( - const int &msh, //number of grid points - const double *r,//radial grid - const double &q, // - const int &l, //angular momentum - double *sj, //jl(1:msh) = j_l(q*r(i)),spherical bessel function - double *sjp -) -{ - ModuleBase::timer::tick ("Mathzone_Add1","Spherical_Bessel"); - - assert (l <= Mathzone_Add1::sph_lmax); - - //creat coefficients - if (!flag_jlx_expand_coef) - { - Mathzone_Add1::expand_coef_jlx (); - flag_jlx_expand_coef = true; - } - - //epsilon - const double eps = 1.0E-10; - - /****************************************************************** - jlx = \sum_{n=0}^{l}(c_{ln}^{s}sin(x) + c_{ln}^{c}cos(x))/x^{l-n+1} - jldx = \sum_{n=0}^{l}(c_{ln}^{s}sin(x) + c_{ln}^{c}cos(x))/x^{l-n+1} - *******************************************************************/ - for (int ir = 0; ir < msh; ir++) - { - double qr = q * r[ir]; - - //judge qr approx 0 - if (fabs(qr) < eps) - { - if (l == 0) sj[ir] = 1.0; - else sj[ir] = 0.0; - - if (l == 1) sjp[ir] = 1.0 / 3.0; - else sjp[ir] = 0.0; - } - else - { - sj[ir] = 0.0; - sjp[ir] = 0.0; - - double lqr = pow (qr, l); - - //divided fac - double xqr = 1.0; - - for (int n = 0; n <= l; n++) - { - double com1 = (c_ln_s[l][n] * sin(qr) + c_ln_c[l][n] * cos(qr)) * xqr; - double com2 = (c_ln_s[l][n] * cos(qr) - c_ln_c[l][n] * sin(qr)) * xqr * qr - + (c_ln_s[l][n] * sin(qr) + c_ln_c[l][n] * cos(qr)) * n * xqr; - - sj[ir] += com1; - sjp[ir] += (com2 - l*com1); - - xqr *= qr; - } - - sj[ir] /= lqr; - sj[ir] /= (lqr * qr); - } - } - - ModuleBase::timer::tick ("Mathzone_Add1","Spherical_Bessel"); - - return; - -} - -double Mathzone_Add1::uni_simpson -( - const double* func, - const int& mshr, - const double& dr -) -{ - ModuleBase::timer::tick("Mathzone_Add1","uni_simpson"); - - assert(mshr >= 3); - - int ir=0; - int idx=0; - int msh_left=0; - double sum = 0.0; - - //f(a) - sum += func[0]; - - //simpson 3/8 rule - if(mshr % 2 == 0) - { - //simpson 3/8 rule - sum += 9.0 / 8.0 * (func[mshr-1] + 3.0*(func[mshr-2] + func[mshr-3]) + func[mshr-4]); - msh_left = mshr - 3; - } - else - { - msh_left = mshr; - } - - //f(b) - sum += func[msh_left-1]; - - //points left - for(ir = 0; ir < msh_left / 2; ir++) - { - idx = 2*ir+1; - sum += 4.0 * func[idx]; - } - - for(ir = 1; ir < msh_left / 2; ir++) - { - idx = 2*ir; - sum += 2.0 * func[idx]; - } - - ModuleBase::timer::tick("Mathzone_Add1","uni_simpson"); - return sum * dr / 3.0; -} - -void Mathzone_Add1::uni_radfft -( - const int& n, - const int& aml, - const int& mshk, - const double* arr_k, - const int& mshr, - const double* ri, - const double& dr, - const double* phir, - double* phik -) -{ - ModuleBase::timer::tick ("Mathzone_Add1","uni_radfft"); - - //allocate memory - double* fi = new double[mshr]; - double* fi_cp = new double[mshr]; - double* jl = new double[mshr]; - - //function to be integrated: r^2 phir - for (int ir = 0; ir < mshr; ir++) - { - fi_cp[ir] = pow(ri[ir], n) * phir[ir]; - } - - //integration - for (int ik = 0; ik < mshk; ik++) - { - //calculate spherical bessel - ModuleBase::Sphbes::Spherical_Bessel(mshr, ri, arr_k[ik], aml, jl); - - //functions to be integrated - for (int ir = 0; ir < mshr; ir++) - { - fi[ir] = fi_cp[ir] * jl[ir]; - } - - phik[ik] = uni_simpson (fi, mshr, dr); - } - - //deallocate memory - delete[] fi; - delete[] fi_cp; - delete[] jl; - - ModuleBase::timer::tick("Mathzone_Add1","uni_radfft"); - - return; -} - -void Mathzone_Add1::Sph_Bes (double x, int lmax, double *sb, double *dsb) -{ - ModuleBase::timer::tick("Mathzone_Add1","Spherical_Bessel"); - - int m, n, nmax; - double j0, j1, sf, tmp, si, co, ix; // j0p, j1p, ix2 - - if (x < 0.0) - { - std::cout << "\nminus x is invalid for Spherical_Bessel" << std::endl; - exit(0); // mohan add 2021-05-06 - } - - const double xmin = 1E-10; - - /* find an appropriate nmax */ - nmax = lmax + 3*static_cast(x) + 20; - if (nmax < 100) nmax = 100; - - /* allocate tsb */ - double* tsb = new double[nmax+1]; - - /* if x is larger than xmin */ - - if (xmin < x) - { - /* initial values*/ - tsb[nmax] = 0.0; - tsb[nmax-1] = 1.0e-14; - - /* downward recurrence from nmax-2 to lmax+2 */ - - for (n = nmax-1; (lmax+2) < n; n--) - { - tsb[n-1] = (2.0*n + 1.0)/x*tsb[n] - tsb[n+1]; - - if (1.0e+250 < tsb[n-1]) - { - tmp = tsb[n-1]; - tsb[n-1] /= tmp; - tsb[n ] /= tmp; - } - } - - /* downward recurrence from lmax+1 to 0 */ - n = lmax + 3; - tmp = tsb[n-1]; - tsb[n-1] /= tmp; - tsb[n ] /= tmp; - - for (n = lmax+2; 0 < n; n--) - { - tsb[n-1] = (2.0*n + 1.0)/x*tsb[n] - tsb[n+1]; - - if (1.0e+250 < tsb[n-1]) - { - tmp = tsb[n-1]; - for (m = n-1; m <= lmax+1; m++) - { - tsb[m] /= tmp; - } - } - } - - /* normalization */ - si = sin(x); - co = cos(x); - ix = 1.0/x; - //ix2 = ix*ix; - j0 = si*ix; - j1 = si*ix*ix - co*ix; - - if (fabs(tsb[1]) < fabs(tsb[0])) sf = j0/tsb[0]; - else sf = j1/tsb[1]; - - /* tsb to sb */ -// for (n = 0; n <= lmax+1; n++) - for (n = 0; n <= lmax; n++) - { - sb[n] = tsb[n]*sf; - } - - /* derivative of sb */ - dsb[0] = co*ix - si*ix*ix; - // for (n = 1; n <= lmax; n++) - for (n = 1; n < lmax; n++) - { - dsb[n] = ( (double)n*sb[n-1] - (double)(n+1.0)*sb[n+1] )/(2.0*(double)n + 1.0); - } - - n = lmax; - dsb[n] = ( (double)n*sb[n-1] - (double)(n+1.0)*sf*tsb[n+1] )/(2.0*(double)n + 1.0); - } - - /* if x is smaller than xmin */ - else - { - /* sb */ - for (n = 0; n <= lmax; n++ ) - { - sb[n] = 0.0; - } - sb[0] = 1.0; - - /* derivative of sb */ - dsb[0] = 0.0; - dsb[1] = 1.0 / 3.0; - for (n = 2; n <= lmax; n++) - { -// dsb[n] = ( (double)n*sb[n-1] - (double)(n+1.0)*sb[n+1] )/(2.0*(double)n + 1.0); - dsb[n] = 0.0; - } - } - - /* free tsb */ - delete [] tsb; - - ModuleBase::timer::tick("Mathzone_Add1","Spherical_Bessel"); - - return; -} - -void Mathzone_Add1::Sbt_new -( - const int& polint_order, - const int& l, - const double* k, - const double& dk, - const int& mshk, - const double* r, - const double& dr, - const int& mshr, - const double* fr, - const int& rpow, - double* fk -) -{ - ModuleBase::timer::tick ("Mathzone_Add1","Sbt_new"); - - //check parameter - assert (mshr >= 1); - assert (l >= 0); - assert (rpow >= 0 && rpow <=2); - - //step 0 - //l is odd or even - bool parity_flag; - if (l % 2 == 0) - { - parity_flag = true; - } - else - { - parity_flag = false; - } - - ModuleBase::GlobalFunc::ZEROS (fk, mshk); - - if (polint_order != 3) - { - std::cout << "\nhigh order interpolation is not available!" << std::endl; - exit(0); // mohan add 2021-05-06 - //ModuleBase::QUIT(); - } - - /********************************** - function multiplied by power of r - for different polint_order - **********************************/ - double* fr2; - double* fr3; - - //polint_order == 1 - fr2 = new double[mshr]; - if (rpow == 0) for (int ir = 0; ir < mshr; ir++) fr2[ir] = fr[ir]*r[ir]*r[ir]; - else if (rpow == 1) for (int ir = 0; ir < mshr; ir++) fr2[ir] = fr[ir]*r[ir]; - else if (rpow == 2) for (int ir = 0; ir < mshr; ir++) fr2[ir] = fr[ir]; - - fr3 = new double[mshr]; - for (int ir = 0; ir < mshr; ir++) fr3[ir] = fr2[ir] * r[ir]; - -// const int polint_order = 3; - int nu_pol_coef = (polint_order+1)*(mshk-1); - double* polint_coef = new double[nu_pol_coef]; - - //step 1 - //start calc - if (parity_flag) - { - //even - const int n = l/2; - - //coef for interpolation - int ct = 0; - - double ft_save = fourier_cosine_transform (fr2, r, mshr, dr, k[0]); - double dft_save = -fourier_sine_transform (fr3, r, mshr, dr, k[0]); - - for (int ik = 0; ik < mshk-1; ik++) - { - double ft0 = ft_save; - double dft0 = dft_save; - - double ft1 = fourier_cosine_transform (fr2, r, mshr, dr, k[ik+1]); - double dft1 = -fourier_sine_transform (fr3, r, mshr, dr, k[ik+1]); - - //double d2k = dk*dk; - //double d3k = d2k*dk; - - double c0 = ft0; - double c1 = dft0; - double c2 = 3.0*(ft1-ft0)/dk/dk-(dft1+2.0*dft0)/dk; - double c3 = (-2.0*(ft1-ft0)/dk+(dft1+dft0))/dk/dk; - - double k2 = k[ik]*k[ik]; - double k3 = k2*k[ik]; - - polint_coef[ct] = c0-c1*k[ik]+c2*k2-c3*k3; - polint_coef[ct+1] = c1-2.0*c2*k[ik]+3.0*c3*k2; - polint_coef[ct+2] = c2-3.0*k[ik]*c3; - polint_coef[ct+3] = c3; - - //test - // double x = (k[ik]+k[ik+1])/2; - // double x = k[ik]; -// double tmp = polint_coef[ct]+x*polint_coef[ct+1]+x*x*polint_coef[ct+2]+x*x*x*polint_coef[ct+3]; - // double tmp_ana = fourier_cosine_transform (fr2, r, mshr, dr, x); - // std::cout << "\ninterp = " << tmp << " ana = " << tmp_ana << " diff = " << log(fabs(tmp-tmp_ana))/log(10); - - //update - ct += (polint_order+1); - ft_save = ft1; - dft_save = dft1; - } - - //store coefficients for calculation - double* coef = new double[n+1]; - double fac = dualfac (l-1) / dualfac (l); - - for (int j = 0; j < n; j++) - { - coef[j] = fac; - - //update - int twoj = 2*j; - fac *= -static_cast((l+twoj+1)*(l-twoj))/(twoj+2)/(twoj+1); - } - coef[n] = pow (-1.0, n) * dualfac (2*l-1) / factorial (l); - - //start calc - //special case k = 0; - if (n ==0) fk[0] = uni_simpson (fr2, mshr, dr); - else fk[0] = 0.0; - - //k > 0 - for (int j =0; j <=n; j++) - { - //initialize Snm - double Snm = 0.0; - for (int ik = 1; ik < mshk; ik++) - { - Snm += pol_seg_int (polint_order, polint_coef, 2*j, k, ik); - double k2j = pow (k[ik], 2*j+1); - - fk[ik] += coef[j] * Snm / k2j; - } - } - //free - delete [] coef; - } - else - { - //odd - const int n = (l-1)/2; - - //coef for interpolation - int ct = 0; - double ft_save, dft_save; - ft_save = fourier_sine_transform (fr2, r, mshr, dr, k[0]); - dft_save = fourier_cosine_transform (fr3, r, mshr, dr, k[0]); - - for (int ik = 0; ik < mshk-1; ik++) - { - double ft0 = ft_save; - double dft0 = dft_save; - - double ft1 = fourier_sine_transform (fr2, r, mshr, dr, k[ik+1]); - double dft1 = fourier_cosine_transform (fr3, r, mshr, dr, k[ik+1]); - - double c0, c1, c2, c3; - c0 = ft0; - c1 = dft0; - c2 = 3.0*(ft1-ft0)/dk/dk-(dft1+2.0*dft0)/dk; - c3 = (-2.0*(ft1-ft0)/dk+(dft1+dft0))/dk/dk; - - double k2 = k[ik]*k[ik]; - double k3 = k2*k[ik]; - - polint_coef[ct] = c0-c1*k[ik]+c2*k2-c3*k3; - polint_coef[ct+1] = c1-2.0*c2*k[ik]+3.0*c3*k2; - polint_coef[ct+2] = c2-3.0*k[ik]*c3; - polint_coef[ct+3] = c3; - -// double x = (k[ik]+k[ik+1])/2; -// double x = k[ik]; -// double tmp = polint_coef[ct]+x*polint_coef[ct+1]+x*x*polint_coef[ct+2]+x*x*x*polint_coef[ct+3]; -// std::cout << "\ninterp = " << tmp << " ana = " << fourier_sine_transform (fr2, r, mshr, dr, x); - //update - ct += (polint_order+1); - ft_save = ft1; - dft_save = dft1; - } - - //store coefficients for calculation - double* coef = new double[n+1]; - double fac = dualfac (l) / dualfac (l-1); - - for (int j = 0; j < n; j++) - { - coef[j] = fac; - - //update - int twoj = 2*j; - fac *= -static_cast((l+twoj+2)*(l-twoj-1))/(twoj+3)/(twoj+2); - - //test -// std::cout << "\ncoef[j] = " << coef[j] << std::endl; - } - coef[n] = pow (-1.0, n) * dualfac (2*l-1) / factorial (l); - - //start calc - //special case k =0 ; - fk[0] = 0.0; - - //k > 0 - for (int j = 0; j <= n; j++) - { - double Snm = 0.0; - for (int ik = 1; ik < mshk; ik++) - { - Snm += pol_seg_int (polint_order, polint_coef, 2*j+1, k, ik); - double k2j = pow(k[ik], 2*j+2); - - fk[ik] += coef[j] * Snm / k2j; - } - } - - //free - delete [] coef; - } - - delete [] fr2; - delete [] fr3; - delete [] polint_coef; - - ModuleBase::timer::tick ("Mathzone_Add1","Sbt_new"); - return; -} - double Mathzone_Add1::factorial (const int& l) { if (l == 0 || l == 1) return 1.0; @@ -650,248 +50,6 @@ double Mathzone_Add1::dualfac (const int& l) else return l * dualfac (l-2); } -double Mathzone_Add1::pol_seg_int -( - const int& polint_order, - const double* coef, - const int& n, - const double* k, - const int& ik -) -{ - double val = 0.0; - double kmf = pow (k[ik], n+1); - double kmb = pow (k[ik-1], n+1); - - int cstart = (polint_order+1)*(ik-1); - for (int i = 0; i <= polint_order; i++) - { - val += coef[cstart+i]*(kmf-kmb)/(n+i+1); - /* - if (ik == 110) - { - std::cout << "i = " << i << " coef = " << coef[cstart+i] << " df = " << kmf-kmb << std::endl; - } - */ - //update - kmf *= k[ik]; - kmb *= k[ik-1]; - } - - /* - if (ik == 110) - { - std::cout << "val = " << val << std::endl; - ModuleBase::QUIT (); - } - */ - return val; -} - -double Mathzone_Add1::fourier_sine_transform -( - const double* func, - const double* r, - const int& mshr, - const double& dr, - const double& k -) -{ - ModuleBase::timer::tick ("Mathzone_Add1","Fsin"); - double val = 0.0; - double* sinf = new double[mshr]; - for (int ir = 0; ir < mshr; ir++) - { - sinf[ir] = func[ir] * sin(k*r[ir]); - } - val = uni_simpson (sinf, mshr, dr); - delete[] sinf; - ModuleBase::timer::tick ("Mathzone_Add1","Fsin"); - return val; -} - -double Mathzone_Add1::fourier_cosine_transform -( - const double* func, - const double* r, - const int& mshr, - const double& dr, - const double& k -) -{ - ModuleBase::timer::tick ("Mathzone_Add1","Fcos"); - double val = 0.0; - double* cosf = new double[mshr]; - for (int ir = 0; ir < mshr; ir++) - { - cosf[ir] = func[ir] * cos(k*r[ir]); - } - - val = uni_simpson (cosf, mshr, dr); - delete[] cosf; - ModuleBase::timer::tick ("Mathzone_Add1","Fcos"); - return val; -} - -void Mathzone_Add1::test () -{ - int polint_order =3; - int dim = 2048; - int ci = 1; - int l = 0; - double rmax = 20; - double dr = rmax/dim; - - double* rad = new double[dim]; - double* func = new double[dim]; - double* fk = new double[dim]; - for (int ir = 0; ir < dim; ir++) - { - rad[ir] = ir * dr; - func[ir] = pow(rad[ir], l) * exp(-ci*rad[ir]*rad[ir]); - fk[ir] = 0.0; - } - - Sbt_new (polint_order, l, rad, dr, dim, rad, dr, dim, func, 0, fk); - - for (int ik = 0; ik < dim; ik++) - { - double diff = fk[ik]- sqrt(ModuleBase::PI/4/ci)/pow(2.0*ci, l+1)* std::pow(rad[ik], l) * exp(-rad[ik]*rad[ik]/4/ci); - std::cout << rad[ik] << " " << fk[ik] << " " << sqrt(ModuleBase::PI/4/ci)/pow(2.0*ci, l+1)*pow(rad[ik], l)*exp(-rad[ik]*rad[ik]/4/ci) - << " " << std::log(fabs(diff))/std::log(10.0) << std::endl; - } - - delete[] rad; - delete[] func; - delete[] fk; - return; -} - -void Mathzone_Add1::test2 () -{ - int polint_order =3; - int N = 200; - int ci = 1; - int l = 0; - double rmax = 20; - double dr = rmax/(N-1); - - double dk = ModuleBase::PI / rmax /2; -// double kmax = PI / dr; -// double dk = dr; - - double* rad = new double[N]; - double* kad = new double[N]; - double* func = new double[N]; - double* fk = new double[N]; - double* fr = new double[N]; - for (int ir = 0; ir < N; ir++) - { - rad[ir] = ir * dr; - kad[ir] = ir * dk; - func[ir] = pow(rad[ir], l) * exp(-ci*rad[ir]*rad[ir]); - fk[ir] = 0.0; - fr[ir] = 0.0; - } - - Sbt_new (polint_order, l, kad, dk, N, rad, dr, N, func, 0, fk); - -/* - for (int ik = 0; ik < N; ik++) - { - double diff = fk[ik]- sqrt(PI/4/ci)/pow(2*ci, l+1)*pow(kad[ik], l)*exp(-kad[ik]*kad[ik]/4/ci); - std::cout << kad[ik] << " " << fk[ik] << " " << sqrt(PI/4/ci)/pow(2*ci, l+1)*pow(kad[ik], l)*exp(-kad[ik]*kad[ik]/4/ci) - << " " << log(fabs(diff))/log(10) << std::endl; - } - ModuleBase::QUIT (); -*/ - Sbt_new (polint_order, l, rad, dr, N, kad, dk, N, fk, 0, fr); - - for (int ir = 0; ir < N; ir++) - { - std::cout << ir*dr << " " << func[ir] << " " << fr[ir] *2.0 / ModuleBase::PI << " " << std::log(fabs(fr[ir]*2.0/ModuleBase::PI-func[ir]))/std::log(10.0) << std::endl; - } - - - delete[] rad; - delete[] func; - delete[] fk; - delete[] fr; - delete[] kad; - return; -} - -double Mathzone_Add1::Polynomial_Interpolation -( - const double* xa, - const double* ya, - const int& n, - const double& x -) -{ - ModuleBase::timer::tick("Mathzone_Add1","Polynomial_Interpolation"); - - int i, m, ns; - double den, dif, dift, ho, hp, w, rs, drs; - - //zero offset - const double* Cxa = xa - 1; - const double* Cya = ya - 1; - - double* cn = new double[n+1]; - double* dn = new double[n+1]; - - ns = 1; - dif = fabs(x - Cxa[1]); - - for(i = 1; i <= n; i++) - { - dift = fabs(x - Cxa[i]); - if(dift < dif) - { - ns = i; - dif = dift; - } - cn[i] = Cya[i]; - dn[i] = Cya[i]; - } - - rs = Cya[ns--]; - - for(m = 1; m < n; m++) - { - for(i = 1; i <= n-m; i++) - { - ho = Cxa[i] - x; - hp = Cxa[i+m] - x; - w = cn[i+1] - dn[i]; - - den = ho - hp; - if(den == 0.0) - { - std::cout << "Two Xs are equal" << std::endl; - // ModuleBase::WARNING_QUIT("Mathzone_Add1::Polynomial_Interpolation","Two Xs are equal"); - exit(0); // mohan update 2021-05-06 - } - den = w / den; - - dn[i] = hp * den; - cn[i] = ho * den; - } - if(2 * ns < n-m) drs = cn[ns+1]; - else drs = dn[ns--]; - - rs += drs; - } - - delete[] cn; - delete[] dn; - - return rs; - ModuleBase::timer::tick("Mathzone_Add1","Polynomial_Interpolation"); - -} - void Mathzone_Add1::SplineD2 // modified by pengfei 13-8-8 add second derivative as a condition ( @@ -1024,124 +182,7 @@ void Mathzone_Add1::Cubic_Spline_Interpolation ModuleBase::timer::tick("Mathzone_Add1","Cubic_Spline_Interpolation"); } -// Interpolation for Numerical Orbitals -double Mathzone_Add1::RadialF -( - const double* rad, - const double* rad_f, - const int& msh, - const int& l, - const double& R -) -{ - ModuleBase::timer::tick("Mathzone_Add1","RadialF"); - - int mp_min, mp_max, m; - double h1, h2, h3, f1, f2, f3, f4; - double g1, g2, x1, x2, y1, y2, f; - double c, result; - - mp_min = 0; - mp_max = msh - 1; - - //assume psir behaves like r**l - if (R < rad[0]) - { - if (l == 0) - { - f = rad_f[0]; - } - else - { - c = rad_f[0] / pow(rad[0], l); - f = pow(R, l) * c; - } - } - else if (rad[mp_max] < R) - { - f = 0.0; - } - else - { - do - { - m = (mp_min + mp_max)/2; - if (rad[m] < R) mp_min = m; - else mp_max = m; - } - while((mp_max-mp_min)!=1); - m = mp_max; - - if (m < 2) - { - m = 2; - } - else if (msh <= m) - { - m = msh - 2; - } - - /**************************************************** - Spline like interpolation - ****************************************************/ - - if (m == 1) - { - h2 = rad[m] - rad[m-1]; - h3 = rad[m+1] - rad[m]; - - f2 = rad_f[m-1]; - f3 = rad_f[m]; - f4 = rad_f[m+1]; - - h1 = -(h2+h3); - f1 = f4; - } - else if (m == (msh-1)) - { - h1 = rad[m-1] - rad[m-2]; - h2 = rad[m] - rad[m-1]; - - f1 = rad_f[m-2]; - f2 = rad_f[m-1]; - f3 = rad_f[m]; - - h3 = -(h1+h2); - f4 = f1; - } - else - { - h1 = rad[m-1] - rad[m-2]; - h2 = rad[m] - rad[m-1]; - h3 = rad[m+1] - rad[m]; - - f1 = rad_f[m-2]; - f2 = rad_f[m-1]; - f3 = rad_f[m]; - f4 = rad_f[m+1]; - } - - //Calculate the value at R - - g1 = ((f3-f2)*h1/h2 + (f2-f1)*h2/h1)/(h1+h2); - g2 = ((f4-f3)*h2/h3 + (f3-f2)*h3/h2)/(h2+h3); - - x1 = R - rad[m-1]; - x2 = R - rad[m]; - y1 = x1/h2; - y2 = x2/h2; - - f = y2*y2*(3.0*f2 + h2*g1 + (2.0*f2 + h2*g1)*y2) - + y1*y1*(3.0*f3 - h2*g2 - (2.0*f3 - h2*g2)*y1); - } - - result = f; - - ModuleBase::timer::tick("Mathzone_Add1","RadialF"); - return result; -} - -// Interpolation for Numerical Orbitals +/// Interpolation for Numerical Orbitals double Mathzone_Add1::Uni_RadialF ( const double* old_phi, @@ -1368,4 +409,4 @@ void Mathzone_Add1::Uni_Deriv_Phi ModuleBase::timer::tick("Mathzone_Add1", "Uni_Deriv_Phi"); } -} \ No newline at end of file +} diff --git a/source/module_base/mathzone_add1.h b/source/module_base/mathzone_add1.h index 68d76560d4..05565ff3c2 100644 --- a/source/module_base/mathzone_add1.h +++ b/source/module_base/mathzone_add1.h @@ -1,8 +1,8 @@ #ifndef MATHZONE_ADD1_H #define MATHZONE_ADD1_H -#include #include +#include namespace ModuleBase { @@ -10,188 +10,74 @@ namespace ModuleBase /************************************************************************ LiaoChen add @ 2010/03/09 to add efficient functions in LCAO calculation ************************************************************************/ -//Only used in module_orbital +// Only used in module_orbital class Mathzone_Add1 { -public: - + public: Mathzone_Add1(); ~Mathzone_Add1(); - static void Sph_Bes (double x, int lmax, double *sb, double *dsb); - - //calculate jlx and its derivatives - static void Spherical_Bessel - ( - const int &msh, //number of grid points - const double *r,//radial grid - const double &q, // - const int &l, //angular momentum - double *sj, //jl(1:msh) = j_l(q*r(i)),spherical bessel function - double *sjp - ); - - /*********************************************** - function : integrate function within [a,b] - formula : \int func(r) dr - ***********************************************/ - static double uni_simpson (const double* func, const int& mshr, const double& dr); - - /******************************************************** - function : radial fourier transform - formula : F(k)_{n,l} = \int r^{n} dr jl(kr) \phi(r) - input : - int n (exponential index) - int aml (angular momentum) - int mshk (number of kpoints) - double* arr_k (an array of kpoints) - int mshr (number of grids) - double* ri (radial points) - double dr (delta r) - double* phir (function to be transformed) - output : - double* phik (with size of mshk) - ********************************************************/ - static void uni_radfft - ( - const int& n, - const int& aml, - const int& mshk, - const double* arr_k, - const int& mshr, - const double* ri, - const double& dr, - const double* phir, - double* phik - ); - - static void Sbt_new - ( - const int& polint_order, - const int& l, - const double* k, - const double& dk, - const int& mshk, - const double* r, - const double& dr, - const int& mshr, - const double* fr, - const int& rpow, - double* fk - ); + static double dualfac(const int& l); + static double factorial(const int& l); + /** + * @brief calculate second derivatives for cubic + * spline interpolation + * + * @param[in] rad x before interpolation + * @param[in] rad_f f(x) before interpolation + * @param[in] mesh number of x before interpolation + * @param[in] yp1 f'(0) boundary condition + * @param[in] ypn f'(n) boundary condition + * @param[out] y2 f''(x) + */ + static void SplineD2(const double* rad, + const double* rad_f, + const int& mesh, + const double& yp1, + const double& ypn, + double* y2); - static double dualfac (const int& l); - static double factorial (const int& l); - static double fourier_cosine_transform - ( - const double* func, - const double* r, - const int& mshr, - const double& dr, - const double& k - ); - - static double fourier_sine_transform - ( - const double* func, - const double* r, - const int& mshr, - const double& dr, - const double& k - ); - -/* - static void RecursiveRadfft - ( - const int& l, - const double** sb, - const double** dsb, - const double* arr_k, - const double* kpoint, - const int& kmesh, - const double& dk, - double& rs, - double& drs - ); -*/ + /** + * @brief cubic spline interpolation + * + * @param[in] rad x before interpolation + * @param[in] rad_f f(x) before inpterpolation + * @param[in] y2 f''(x) before interpolation + * @param[in] mesh number of x before interpolation + * @param[in] r x after interpolation + * @param[in] rsize number of x after interpolation + * @param[out] y f(x) after interpolation + * @param[out] dy f'(x) after interpolation + */ + static void Cubic_Spline_Interpolation(const double* const rad, + const double* const rad_f, + const double* const y2, + const int& mesh, + const double* const r, + const int& rsize, + double* const y, + double* const dy); - static double Polynomial_Interpolation - ( - const double* xa, - const double* ya, - const int& n, - const double& x - ); - - static void SplineD2 - ( - const double *rad, - const double *rad_f, - const int& mesh, - const double &yp1, - const double &ypn, - double* y2 - ); + /** + * @brief "spline like interpolation" of a uniform + * funcation of r + * + * @param[in] rad_f f(x) before interpolation + * @param[in] msh number of x known + * @param[in] dr uniform distance of x + * @param R f(R) is to be calculated + * @return double f(R) + */ + static double Uni_RadialF(const double* rad_f, const int& msh, const double& dr, const double& R); - static void Cubic_Spline_Interpolation - ( - const double * const rad, - const double * const rad_f, - const double * const y2, - const int& mesh, - const double * const r, - const int& rsize, - double * const y, - double * const dy - ); - - static double RadialF - ( - const double* rad, - const double* rad_f, - const int& msh, - const int& l, - const double& R - ); - - static double Uni_RadialF - ( - const double* rad_f, - const int& msh, - const double& dr, - const double& R - ); - - static void test (); - static void test2 (); - - - static void Uni_Deriv_Phi - ( - const double *radf, - const int &mesh, - const double &dr, - const int &nd, - double* phind - ); - - private: - const static int sph_lmax = 20; - static double** c_ln_c; - static double** c_ln_s; - static bool flag_jlx_expand_coef; + static void Uni_Deriv_Phi(const double* radf, const int& mesh, const double& dr, const int& nd, double* phind); - static void expand_coef_jlx (); - static double pol_seg_int - ( - const int& polint_order, - const double* coef, - const int& n, - const double* k, - const int& ik - ); + private: + const static int sph_lmax = 20; + static double** c_ln_c; + static double** c_ln_s; }; -} +} // namespace ModuleBase #endif diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt index ddfa687934..7fa8a6c5e2 100644 --- a/source/module_base/test/CMakeLists.txt +++ b/source/module_base/test/CMakeLists.txt @@ -72,8 +72,13 @@ AddTest( TARGET base_math_polyint SOURCES math_polyint_test.cpp ../math_polyint.cpp ../realarray.cpp ../timer.cpp ) +AddTest( + TARGET base_mathzone_add1 + LIBS ${math_libs} + SOURCES mathzone_add1_test.cpp ../mathzone_add1.cpp ../math_sphbes.cpp ../mathzone.cpp ../matrix3.cpp ../matrix.cpp ../tool_quit.cpp ../global_variable.cpp ../global_file.cpp ../memory.cpp ../timer.cpp +) AddTest( TARGET base_math_ylmreal LIBS ${math_libs} SOURCES math_ylmreal_test.cpp ../math_ylmreal.cpp ../realarray.cpp ../timer.cpp ../matrix.cpp -) \ No newline at end of file +) diff --git a/source/module_base/test/mathzone_add1_test.cpp b/source/module_base/test/mathzone_add1_test.cpp new file mode 100644 index 0000000000..c10c86c008 --- /dev/null +++ b/source/module_base/test/mathzone_add1_test.cpp @@ -0,0 +1,231 @@ +#include "../mathzone_add1.h" +#include "gtest/gtest.h" +#include "gmock/gmock.h" + +/************************************************ + * unit test of class Mathzone_Add1 + ***********************************************/ + +/** + * - Tested Functions: + * - CubicSplineBoundary1 + * - call SplineD2 and Cubic_Spline_Interpolation + * - to interpolate a function using the first + * - kind of boundary condition:f'(0) = f'(n) = 0.0 + * - CubicSplineBoundary2 + * - call SplineD2 and Cubic_Spline_Interpolation + * - to interpolate a function using the second + * - kind of boundary condition:f''(0) = f''(n) = 0.0 + * - UniRadialF + * - call Uni_RadialF to interpolate the radial part + * - of an atomic orbital function, whose discrete r + * - points are uniform. + * - Factorial + * - calculate the factorial of an integer + * - DualFac + * - calculate the double factorial or + * - semifactorial of an integer + */ + +class MathzoneAdd1Test : public testing::Test +{ +protected: + const int MaxInt = 100; + int nr_in = 17; + int nr_out = 161; + double *r_in; + double *r_out; + double *y2; + double *psi_in; + double *psi_out; + double *dpsi; + void SetUp() + { + r_in = new double[nr_in]; + r_out = new double[nr_out]; + y2 = new double[nr_in]; + psi_in = new double[nr_in]; + psi_out = new double[nr_out]; + dpsi = new double[nr_out]; + } + void TearDown() + { + delete[] r_in; + delete[] r_out; + delete[] y2; + delete[] psi_in; + delete[] psi_out; + delete[] dpsi; + } +}; + +/// first kind boundary condition: f'(0) = f'(n) = 0.0 +TEST_F(MathzoneAdd1Test, CubicSplineBoundary1) +{ + // data from abacus/tests/integrate/tools/PP_ORB/Si_gga_8au_60Ry_2s2p1d.orb + // data for d orbital of Si : L = 2, N = 0 + psi_in[0] = 0; + psi_in[1] = -2.583946346740e-01; + psi_in[2] = -4.570087269049e-01; + psi_in[3] = -4.374680500187e-01; + psi_in[4] = -3.587829079989e-01; + psi_in[5] = -2.581772323753e-01; + psi_in[6] = -1.616203660437e-01; + psi_in[7] = -9.108838081645e-02; + psi_in[8] = -5.202206559586e-02; + psi_in[9] = -3.126875315134e-02; + psi_in[10] = -1.860199873973e-02; + psi_in[11] = -8.049945178799e-03; + psi_in[12] = -1.652010824028e-03; + psi_in[13] = 1.495515249035e-03; + psi_in[14] = 3.221037475903e-03; + psi_in[15] = 3.802139894646e-03; + psi_in[16] = 0; + for (int i=0; i< nr_in; i++) + { + r_in[i] = i*0.5; + //std::cout<< r_in[i] << " " << psi_in[i] << std::endl; // for plotting + } + for (int i=0; i< nr_out; i++) + { + r_out[i] = i*0.05; + } + ModuleBase::Mathzone_Add1::SplineD2(r_in,psi_in,nr_in,0.0,0.0,y2); + //std::cout << "y2[0] "<< y2[0] << " y2[nr_in] "<< y2[nr_in-1] << std::endl; // for checking + ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(r_in,psi_in,y2,nr_in,r_out,nr_out,psi_out,dpsi); + for (int i=0; i< nr_out; i++) + { + int j = i/10; + if(i%10==0) { + EXPECT_EQ(psi_in[j],psi_out[i]); + } + //std::cout<< r_out[i] << " " << psi_out[i] << std::endl; // for plotting + } + EXPECT_NEAR(dpsi[0],0.0,1e-15); + EXPECT_NEAR(dpsi[nr_out-1],0.0,1e-15); + //std::cout<