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);
117 static GFit BSHFit(
double mu,
double lo,
double hi,
double eps,
bool prnt=
false) {
119 bool fix_interval=
false;
140 static GFit SlaterFit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
145 auto exact = [&gamma](
const double r) ->
double {
return exp(-gamma * r); };
160 static GFit GaussFit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
165 fit.exponents_=gamma;
168 auto exact = [&gamma](
const double r) ->
double {
return exp(-gamma * r*r); };
183 static GFit F12Fit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
186 fit.coeffs_*=(0.5/gamma);
189 auto exact=[&gamma](
const double r) ->
double {
return 0.5/gamma*(1.0-exp(-gamma*r));};
204 static GFit F12sqFit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
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);};
227 static GFit FGFit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
228 GFit bshfit,coulombfit;
232 bool fix_interval=
true;
238 auto diffcoefficients=(coulombfit.
coeffs() - bshfit.
coeffs());
246 auto exact=[&gamma](
const double r) ->
double {
return 0.5/gamma*(1.0-exp(-gamma*r))/r;};
261 static GFit F2GFit(
double gamma,
double lo,
double hi,
double eps,
bool prnt=
false) {
262 GFit bshfit,coulombfit,bsh2fit;
266 bool fix_interval=
true;
282 double fourmu2=4.0*gamma*gamma;
283 f2gfit.
coeffs_=fourpi/fourmu2*coefficients;
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;
314 double mid =
lo + (hi-
lo)*0.5;
315 long npt=coeff.
size();
317 for (i=npt-1; i>0; --i) {
318 double cnew = coeff[i]*exp(-(expnt[i]-expnt[i-1])*mid*mid);
319 double errlo = coeff[i]*exp(-expnt[i]*
lo*
lo) -
320 cnew*exp(-expnt[i-1]*
lo*
lo);
321 double errhi = coeff[i]*exp(-expnt[i]*hi*hi) -
322 cnew*exp(-expnt[i-1]*hi*hi);
325 coeff[i-1] = coeff[i-1] + cnew;
327 coeff = coeff(
Slice(0,npt-1));
328 expnt = expnt(
Slice(0,npt-1));
332 double L,
bool discardG0)
const {
333 double tcut = 0.25/
L/
L;
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";
425 if ((
mu > 0) and (not fix_interval)) {
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;
444 if ((
mu > 0) and (not fix_interval)) {
446 slo = -0.5*log(4.0*TT/(
mu*
mu));
449 slo = log(eps/hi) - 1.0;
451 shi = 0.5*log(TT/(
lo*
lo));
452 if (shi <= slo)
throw "bsh_fit: logic error in slo,shi";
455 double h = 1.0/(0.2-.50*log10(eps));
462 h = floor(64.0*
h)/64.0;
467 slo = floor(slo/
h)*
h;
469 long npt = long((shi-slo)/
h+0.5);
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]);
503 if ((
mu == 0.0) and (not fix_interval)) {
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) {
540 qg = qg(
Slice(1,nmom));
544 print(
"moments", qg);
547 for (
int j=0; j<nmom; ++j) {
550 if (nmom != 4) qt = qt(
Slice(1,nmom));
551 for (
int i=0; i<nmom; ++i) {
560 print(
"new coeffs", ncoeff);
563 coeff(
Slice(0,nmom-1)) = ncoeff;
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;
593 double shi=log(TT/(
lo*
lo))*0.5;
596 double h = 1.0/(0.2-.5*log10(eps));
603 h = floor(64.0*
h)/64.0;
607 slo = floor(slo/
h)*
h;
609 long npt = long((shi-slo)/
h+0.5);
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));
633 double exact = exp(-gamma*r);
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,
656 pcoeff(
Slice(1,-1,1))=-coeff(
_);
658 pexpnt(
Slice(1,-1,1))=expnt(
_);
669 static void f12sq_fit(
double gamma,
double lo,
double hi,
double eps,
673 slater_fit(2.0*gamma,
lo, hi, eps, coeff2, expnt2, prnt);
679 auto coeff=coeff2-2.0*coeff1;
681 pcoeff(
Slice(1,-1,1))=coeff(
_);
683 pexpnt(
Slice(1,-1,1))=expnt1(
_);
707 slo = -0.5*log(4.0*100.0/(
mu*
mu));
708 slo = -0.5*log(4.0*(slo*ndim - 2.0*slo + 100.0)/(
mu*
mu));
711 slo = log(eps/hi) - 1.0;
713 shi = 0.5*log(100.0/(
lo*
lo));
716 double h = 1.0/(0.2-.50*log10(eps));
723 h = floor(64.0*
h)/64.0;
728 slo = floor(slo/
h)*
h;
730 long npt = long((shi-slo)/
h+0.5);
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);
746 if (
c*exp(-
p*
lo*
lo) > eps) {
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]);
772 double mid =
lo + (hi-
lo)*0.5;
774 for (i=npt-1; i>0; --i) {
775 double cnew = coeff[i]*exp(-(expnt[i]-expnt[i-1])*mid*mid);
776 double errlo = coeff[i]*exp(-expnt[i]*
lo*
lo) -
777 cnew*exp(-expnt[i-1]*
lo*
lo);
778 double errhi = coeff[i]*exp(-expnt[i]*hi*hi) -
779 cnew*exp(-expnt[i-1]*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]);
804 printf(
" %.6e %20.10e\n",r,
test);
837 * exp(-
mu *
R) / 0.4e1;
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 multidimension 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
const 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:2416
std::string type(const PairType &n)
Definition PNOParameters.h:18
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:65
double thresh
Definition operatorinfo.h:63
bool debug
Definition operatorinfo.h:66
OpType type
introspection
Definition operatorinfo.h:64
double mu
some introspection
Definition operatorinfo.h:61
double lo
Definition operatorinfo.h:62
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
double h(const coord_1d &r)
Definition testgconv.cc:68
static const std::size_t NDIM
Definition testpdiff.cc:42
void test()
Definition y.cc:696