From 84f4a5bb814e1da5b6872b3008b743900447cdbf Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Wed, 27 Dec 2023 14:05:10 -0500 Subject: [PATCH] updates?? --- src/madness/bspline/CMakeLists.txt | 4 +- src/madness/bspline/bspline.cc | 26 ++++---- src/madness/bspline/myqd.h | 87 +++++++++++++++++++++++++-- src/madness/bspline/solidharmonics.cc | 2 +- 4 files changed, 100 insertions(+), 19 deletions(-) diff --git a/src/madness/bspline/CMakeLists.txt b/src/madness/bspline/CMakeLists.txt index 25577a66556..36492c428c9 100644 --- a/src/madness/bspline/CMakeLists.txt +++ b/src/madness/bspline/CMakeLists.txt @@ -1,5 +1,5 @@ set (MADBSPLINE_SOURCES - bspline.cc + bspline.cc madbessel.cc solidharmonics.cc ) set (MADBSPLINE_HEADERS @@ -8,6 +8,8 @@ set (MADBSPLINE_HEADERS ) add_mad_executable(bspline bspline.cc "MADchem;-lqd") +add_mad_executable(solidharmonics solidharmonics.cc "MADchem;-lqd") +add_mad_executable(madbessel madbessel.cc "MADchem;-lqd") ##add_dependencies(applications-madness bspline) diff --git a/src/madness/bspline/bspline.cc b/src/madness/bspline/bspline.cc index 8a973369ba0..a5742cd7d1e 100644 --- a/src/madness/bspline/bspline.cc +++ b/src/madness/bspline/bspline.cc @@ -1074,24 +1074,24 @@ class BsplineFunction { } } - static void test_fit() { - const T Z = 10.0; - size_t nknots = 61; - size_t order = 20; - T xlo = 0.1; // a for KnotsRational //.01 for 103 //1e-7; - T xhi = 26.0; + // static void test_fit() { + // const T Z = 10.0; + // size_t nknots = 61; + // size_t order = 20; + // T xlo = 0.1; // a for KnotsRational //.01 for 103 //1e-7; + // T xhi = 26.0; - bdataT::init(order, nknots, xlo, xhi); + // bdataT::init(order, nknots, xlo, xhi); - print("knots", bdataT::knots()); + // print("knots", bdataT::knots()); - auto A = bdataT::basis_at_GL_points(); - auto [X, W] = bdataT::quadrature(); - size_t ltest = 1; - ????????????????????????????? + // auto A = bdataT::basis_at_GL_points(); + // auto [X, W] = bdataT::quadrature(); + // size_t ltest = 1; + // ????????????????????????????? - } + // } static void test() { const T Z = 10.0; diff --git a/src/madness/bspline/myqd.h b/src/madness/bspline/myqd.h index c7f46696bdb..d10f61cfae6 100644 --- a/src/madness/bspline/myqd.h +++ b/src/madness/bspline/myqd.h @@ -14,6 +14,8 @@ // Inject math routines into std name space to enable templated types namespace std { qd_real abs(const qd_real& x) {return ::abs(x);} + qd_real exp(const qd_real& x) {return ::exp(x);} + qd_real ldexp(const qd_real &a, int exp) {return ::ldexp(a,exp);} qd_real cos(const qd_real& x) {return ::cos(x);} qd_real sin(const qd_real& x) {return ::sin(x);} qd_real pow(const qd_real& x, int n) {return ::pow(x,n);} @@ -22,6 +24,8 @@ namespace std { qd_real sqrt(const qd_real& x) {return ::sqrt(x);} dd_real abs(const dd_real& x) {return ::abs(x);} + dd_real exp(const dd_real& x) {return ::exp(x);} + dd_real ldexp(const dd_real &a, int exp) {return ::ldexp(a,exp);} dd_real cos(const dd_real& x) {return ::cos(x);} dd_real sin(const dd_real& x) {return ::sin(x);} dd_real pow(const dd_real& x, int n) {return ::pow(x,n);} @@ -30,16 +34,87 @@ namespace std { dd_real sqrt(const dd_real& x) {return ::sqrt(x);} } +template +T myread(const std::string& s) { + T x = 0; + int sx=1, sexpnt=1, expnt=0; + int n = 0; + int point = 0; + bool doingexpt = false; + for (char c : s) { + if (doingexpt) { + switch (c) { + case '-': + sexpnt = -1; + break; + + case '+': + break; + + default: + expnt = 10*expnt + int(c-'0'); + } + + } + else { + switch (c) { + case '-': + sx = -1; + break; + + case '.': + point = n; + break; + + case 'e': + case 'E': + doingexpt = true; + break; + + default: + x = T(10)*x + int(c-'0'); + n++; + } + } + } + //std::cout << n << " " << expnt << " " << point << std::endl; + n += -sexpnt*expnt - point; + n = -n; + + while (abs(n)) { + int m; + if (n<0) { + m = std::max(n, -20); + } + else { + m = std::min(n, 20); + } + n -= m; + x *= pow(T(10), m); + } + + if (sx<0) { + x = -x; + } + + return x; +} + +template +T myread(const char* s) { + return myread(std::string(s)); +} + template T from_str(const char* s); template <> qd_real from_str(const char* s) { - return qd_real(s); + return myread(s); } template <> dd_real from_str(const char* s) { - return dd_real(s); + return myread(s); } template <> @@ -60,12 +135,16 @@ float from_str(const char* s) { return d; } +template T from_str(const std::string& s) { + return from_str(s.c_str()); +} + std::string to_str(const qd_real& t) {return t.to_string();} std::string to_str(const dd_real& t) {return t.to_string();} -std::string to_str(double t) {char buf[256]; sprintf(buf,"%.17f",t); return buf;} +std::string to_str(double t) {char buf[256]; sprintf(buf,"%.17e",t); return buf;} +std::string to_str(float t) {char buf[256]; sprintf(buf,"%.8e",t); return buf;} double to_double(double d) {return d;} - double to_double(float f) {return f;} // Needed since double(qd_real) not implemented and I fear what would break if I added it. diff --git a/src/madness/bspline/solidharmonics.cc b/src/madness/bspline/solidharmonics.cc index e091b1fff17..3c0445ca636 100644 --- a/src/madness/bspline/solidharmonics.cc +++ b/src/madness/bspline/solidharmonics.cc @@ -546,7 +546,7 @@ auto threeJ() { //std::cout << std::setw(3) << l1 << " " << std::setw(3) << m1 << " " << std::setw(3) << l2 << " " << std::setw(3) << m2 << " --> " << std::setw(3) << l3 << " " << std::setw(3) << m3 << " " << std::setw(15) << j3value << " " << test << std::endl; j3.push_back({l3,m3,j3value}); - if (l2==1 and m2==0) { + if ((l2==1 and m2==1) or (l1==1 and m2==1)) { T unnorm = j3value*Nnorm.value(l3,m3)/(Nnorm.value(l1,m1)*Nnorm.value(l2,m2)); std::cout << std::setw(3) << l1 << " " << std::setw(3) << m1 << " " << std::setw(3) << l2 << " " << std::setw(3) << m2 << " --> " << std::setw(3) << l3 << " " << std::setw(3) << m3 << " " << std::setw(15) << j3value << " " << unnorm << std::endl; }