35#ifndef MADNESS_MRA_GFIT_H__INCLUDED
36#define MADNESS_MRA_GFIT_H__INCLUDED
44#include "../constants.h"
45#include "../tensor/basetensor.h"
46#include "../tensor/slice.h"
47#include "../tensor/tensor.h"
48#include "../tensor/tensor_lapack.h"
49#include "../world/madness_exception.h"
50#include "../world/print.h"
56template<
typename T, std::
size_t NDIM>
68 MADNESS_CHECK_THROW(hi>0,
"hi must be positive in gfit: U need to set it manually in operator.h");
83 print(
"Operator type not implemented: ",
type);
145 auto exact = [&gamma](
const double r) ->
double {
return exp(-gamma * r); };
165 fit.exponents_=gamma;
168 auto exact = [&gamma](
const double r) ->
double {
return exp(-gamma * r*r); };
186 fit.coeffs_*=(0.5/gamma);
189 auto exact=[&gamma](
const double r) ->
double {
return 0.5/gamma*(1.0-
exp(-gamma*r));};
207 fit.coeffs_*=(0.25/(gamma*gamma));
210 auto exact=[&gamma](
const double r) ->
double {
return std::pow(0.5/gamma*(1.0-
exp(-gamma*r)),2.0);};
246 auto exact=[&gamma](
const double r) ->
double {
return 0.5/gamma*(1.0-
exp(-gamma*r))/r;};
282 double fourmu2=4.0*gamma*gamma;
288 auto exact=[&gamma](
const double r) ->
double {
289 return 0.25/(gamma*gamma)*(1.0-2.0*
exp(-gamma*r)+
exp(-2.0*gamma*r))/r;
315 long npt=coeff.
size();
317 for (i=npt-1; i>0; --i) {
321 double errhi = coeff[i]*
exp(-expnt[i]*hi*hi) -
325 coeff[i-1] = coeff[i-1] +
cnew;
327 coeff = coeff(
Slice(0,npt-1));
328 expnt = expnt(
Slice(0,npt-1));
337 for (
int i=0; i<
e.dim(0); ++i) {
369 template<
typename funcT>
379 template<
typename opT>
382 std::cout <<
"weights and roots" << std::endl;
383 for (
int i=0; i<
coeffs_.size(); ++i)
387 std::cout <<
" x value abserr relerr" << std::endl;
388 std::cout <<
" ------------ ------- -------- -------- " << std::endl;
389 double step =
exp(log(hi/
lo)/(npt+1));
390 for (
int i=0; i<=npt; ++i) {
391 double r =
lo*(
pow(step,i+0.5));
395 for (
int j=0; j<
coeffs_.dim(0); ++j)
421 if (
mu < 0.0)
throw "cannot handle negative mu in bsh_fit";
436 if (eps >= 1
e-2)
TT = 5;
437 else if (eps >= 1
e-4)
TT = 10;
438 else if (eps >= 1
e-6)
TT = 14;
439 else if (eps >= 1
e-8)
TT = 18;
440 else if (eps >= 1
e-10)
TT = 22;
441 else if (eps >= 1
e-12)
TT = 26;
449 slo = log(eps/hi) - 1.0;
452 if (
shi <=
slo)
throw "bsh_fit: logic error in slo,shi";
455 double h = 1.0/(0.2-.50*
log10(eps));
476 for (
int i=0; i<npt; ++i) {
477 double s =
slo +
h*(npt-i);
480 expnt[i] =
exp(2.0*s);
487 expnt[0] =
exp(2.0*s);
488 coeff=coeff(
Slice(0,0));
489 expnt=expnt(
Slice(0,0));
490 print(
"only one term in gfit",s,coeff[0],expnt[0]);
528 double range =
sqrt(-log(1
e-6)/expnt[
nmom-1]);
529 if (
prnt)
print(
"exponent(nmom-1)",expnt[
nmom-1],
"has range", range);
533 for (
int i=
nmom; i<npt; ++i) {
547 for (
int j=0; j<
nmom; ++j) {
551 for (
int i=0; i<
nmom; ++i) {
582 if (eps >= 1
e-2)
TT = 5;
583 else if (eps >= 1
e-4)
TT = 10;
584 else if (eps >= 1
e-6)
TT = 14;
585 else if (eps >= 1
e-8)
TT = 18;
586 else if (eps >= 1
e-10)
TT = 22;
587 else if (eps >= 1
e-12)
TT = 26;
592 double slo=0.5 * log(eps) - 1.0;
596 double h = 1.0/(0.2-.5*
log10(eps));
613 for (
int i=0; i<npt; ++i) {
614 const double s =
slo +
h*(npt-i);
615 coeff[i] =
h*
exp(-gamma*gamma*
exp(2.0*s) + s);
617 expnt[i] = 0.25*
exp(-2.0*s);
621 std::cout <<
"weights and roots for a Slater function with gamma=" << gamma << std::endl;
622 for (
int i=0; i<npt; ++i)
623 std::cout << i <<
" " << coeff[i] <<
" " << expnt[i] << std::endl;
628 std::cout <<
" x value abserr relerr" << std::endl;
629 std::cout <<
" ------------ ------- -------- -------- " << std::endl;
630 double step =
exp(log(hi/
lo)/(npt+1));
631 for (
int i=0; i<=npt; ++i) {
632 double r =
lo*(
pow(step,i+0.5));
635 for (
int j=0; j<coeff.dim(0); ++j)
636 test += coeff[j]*
exp(-r*r*expnt[j]);
650 static void f12_fit(
double gamma,
double lo,
double hi,
double eps,
669 static void f12sq_fit(
double gamma,
double lo,
double hi,
double eps,
707 slo = -0.5*log(4.0*100.0/(
mu*
mu));
711 slo = log(eps/hi) - 1.0;
716 double h = 1.0/(0.2-.50*
log10(eps));
733 std::cout <<
"bsh: mu " <<
mu <<
" lo " <<
lo <<
" hi " << hi
734 <<
" eps " << eps <<
" slo " <<
slo <<
" shi " <<
shi
735 <<
" npt " << npt <<
" h " <<
h << std::endl;
741 for (
int i=0; i<npt; ++i) {
742 double s =
slo +
h*(npt-i);
744 double p =
exp(2.0*s);
757 expnt[0] =
exp(2.0*s);
758 coeff=coeff(
Slice(0,0));
759 expnt=expnt(
Slice(0,0));
760 print(
"only one term in gfit",s,coeff[0],expnt[0]);
774 for (i=npt-1; i>0; --i) {
778 double errhi = coeff[i]*
exp(-expnt[i]*hi*hi) -
782 coeff[i-1] = coeff[i-1] +
cnew;
787 coeff = coeff(
Slice(0,npt-1));
788 expnt = expnt(
Slice(0,npt-1));
792 for (
int i=0; i<npt; ++i)
793 std::cout << i <<
" " << coeff[i] <<
" " << expnt[i] << std::endl;
796 std::cout <<
" x value" << std::endl;
797 std::cout <<
" ------------ ---------------------" << std::endl;
798 double step =
exp(log(hi/
lo)/(npt+1));
799 for (
int i=0; i<=npt; ++i) {
800 double r =
lo*(
pow(step,i+0.5));
802 for (
int j=0; j<coeff.dim(0); ++j)
803 test += coeff[j]*
exp(-r*r*expnt[j]);
838 q[2] = -(-0.2e1 *
exp(
mu *
R) + 0.2e1 + 0.2e1 *
mu *
R +
R*
R *
840 q[3] = -(-0.6e1 *
exp(
mu *
R) + 0.6e1 + 0.6e1 *
mu *
R + 0.3e1 *
R*
R
double q(double t)
Definition DKops.h:18
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
Tensor< T > coeffs() const
return the coefficients of the fit
Definition gfit.h:307
void print_accuracy(opT op, const double lo, const double hi) const
print coefficients and exponents, and values and errors
Definition gfit.h:380
Tensor< T > exponents() const
return the exponents of the fit
Definition gfit.h:310
static GFit GaussFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a (trivial) fit for a single Gauss function
Definition gfit.h:160
GFit()=default
default ctor does nothing
static GFit F2GFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the F2G function
Definition gfit.h:261
GFit(const funcT &f)
ctor taking an isotropic function
Definition gfit.h:370
static void prune_small_coefficients(const double eps, const double lo, const double hi, Tensor< double > &coeff, Tensor< double > &expnt)
Definition gfit.h:312
static void f12_fit(double gamma, double lo, double hi, double eps, Tensor< double > &pcoeff, Tensor< double > &pexpnt, bool prnt)
fit a correlation factor (1- exp(-mu r))
Definition gfit.h:650
static void slater_fit(double gamma, double lo, double hi, double eps, Tensor< double > &pcoeff, Tensor< double > &pexpnt, bool prnt)
fit a Slater function using a sum of Gaussians
Definition gfit.h:576
static void bsh_spherical_moments(double mu, double R, Tensor< double > &q)
Definition gfit.h:827
static GFit BSHFit(double mu, double lo, double hi, double eps, bool prnt=false)
return a fit for the bound-state Helmholtz function
Definition gfit.h:117
static GFit GeneralFit()
return a fit for a general isotropic function
Definition gfit.h:301
Tensor< T > exponents_
the exponents of the expansion f(x) = \sum_m coeffs[m] exp(-exponents[m] * x^2)
Definition gfit.h:407
static GFit FGFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the FG function
Definition gfit.h:227
static void f12sq_fit(double gamma, double lo, double hi, double eps, Tensor< double > &pcoeff, Tensor< double > &pexpnt, bool prnt)
fit a correlation factor f12^2 = (1- exp(-mu r))^2 = 1 - 2 exp(-mu r) + exp(-2 mu r)
Definition gfit.h:669
static GFit F12sqFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the F12^2 correlation factor
Definition gfit.h:204
static void bsh_fit(double mu, double lo, double hi, double eps, Tensor< double > &pcoeff, Tensor< double > &pexpnt, bool prnt, bool fix_interval)
fit the function exp(-mu r)/r
Definition gfit.h:418
static void gaussian_spherical_moments(double alpha, double R, Tensor< double > &q)
Definition gfit.h:813
static GFit F12Fit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the F12 correlation factor
Definition gfit.h:183
GFit & operator=(const GFit &other)
assignment operator
Definition gfit.h:93
static GFit SlaterFit(double gamma, double lo, double hi, double eps, bool prnt=false)
return a fit for the Slater function
Definition gfit.h:140
static GFit CoulombFit(double lo, double hi, double eps, bool prnt=false)
return a fit for the Coulomb function
Definition gfit.h:102
static void bsh_fit_ndim(int ndim, double mu, double lo, double hi, double eps, Tensor< double > &pcoeff, Tensor< double > &pexpnt, bool prnt)
Definition gfit.h:691
GFit(const GFit &other)=default
copy constructor
void truncate_periodic_expansion(Tensor< double > &c, Tensor< double > &e, double L, bool discardG0) const
Definition gfit.h:331
GFit(OperatorInfo info)
Definition gfit.h:64
Tensor< T > coeffs_
the coefficients of the expansion f(x) = \sum_m coeffs[m] exp(-exponents[m] * x^2)
Definition gfit.h:404
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
A tensor is a multidimensional array.
Definition tensor.h:317
static const double R
Definition csqrt.cc:46
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
static double lo
Definition dirac-hatom.cc:23
static bool debug
Definition dirac-hatom.cc:16
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
static double pow(const double *a, const double *b)
Definition lda.h:74
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:207
constexpr double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static const Slice _(0,-1, 1)
OpType
operator types
Definition operatorinfo.h:11
@ OT_FG12
1-exp(-r)
Definition operatorinfo.h:18
@ OT_SLATER
1/r
Definition operatorinfo.h:15
@ OT_GAUSS
exp(-r)
Definition operatorinfo.h:16
@ OT_BSH
(1-exp(-r))^2/r = 1/r + exp(-2r)/r - 2 exp(-r)/r
Definition operatorinfo.h:21
@ OT_F12
exp(-r2)
Definition operatorinfo.h:17
@ OT_F212
(1-exp(-r))/r
Definition operatorinfo.h:19
@ OT_G12
indicates the identity
Definition operatorinfo.h:14
@ OT_F2G12
(1-exp(-r))^2
Definition operatorinfo.h:20
void gesv(const Tensor< T > &a, const Tensor< T > &b, Tensor< T > &x)
Solve Ax = b for general A using the LAPACK *gesv routines.
Definition lapack.cc:804
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
NDIM & f
Definition mra.h:2448
std::string type(const PairType &n)
Definition PNOParameters.h:18
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
static long abs(long a)
Definition tensor.h:218
const double mu
Definition navstokes_cosines.cc:95
static const double c
Definition relops.cc:10
static const double L
Definition rk.cc:46
Definition operatorinfo.h:58
double hi
Definition operatorinfo.h:67
double thresh
Definition operatorinfo.h:65
bool debug
Definition operatorinfo.h:69
OpType type
introspection
Definition operatorinfo.h:66
double mu
some introspection
Definition operatorinfo.h:63
double lo
Definition operatorinfo.h:64
void e()
Definition test_sig.cc:75
static const double alpha
Definition testcosine.cc:10
double(* exact)(double, double, double)
Definition testfuns.cc:6
std::vector< double > fit(size_t m, size_t n, const std::vector< double > N, const std::vector< double > &f)
Definition testfuns.cc:36
constexpr std::size_t NDIM
Definition testgconv.cc:54
double h(const coord_1d &r)
Definition testgconv.cc:175
void test()
Definition y.cc:696