MADNESS  0.10.1
gfit.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31 
32  $Id: key.h 2907 2012-06-14 10:15:05Z 3ru6ruWu $
33  */
34 
35 #ifndef MADNESS_MRA_GFIT_H__INCLUDED
36 #define MADNESS_MRA_GFIT_H__INCLUDED
37 
38 /// \file gfit.h
39 /// \brief fit isotropic functions to a set of Gaussians with controlled precision
40 
41 //#include <iostream>
42 
43 #include <cmath>
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"
52 
53 
54 namespace madness {
55 
56 template<typename T, std::size_t NDIM>
57 class GFit {
58 
59 public:
60 
61  /// default ctor does nothing
62  GFit() = default;
63 
65  double mu = info.mu;
66  double lo = info.lo;
67  double hi = info.hi;
68  MADNESS_CHECK_THROW(hi>0,"hi must be positive in gfit: U need to set it manually in operator.h");
69  double eps = info.thresh;
70  bool debug=info.debug;
71  OpType type = info.type;
72 
73 
74  if (type==OT_G12) {*this=CoulombFit(lo,hi,eps,debug);
75  } else if (type==OT_SLATER) {*this=SlaterFit(mu,lo,hi,eps,debug);
76  } else if (type==OT_GAUSS) {*this=GaussFit(mu,lo,hi,eps,debug);
77  } else if (type==OT_F12) {*this=F12Fit(mu,lo,hi,eps,debug);
78  } else if (type==OT_FG12) {*this=FGFit(mu,lo,hi,eps,debug);
79  } else if (type==OT_F212) {*this=F12sqFit(mu,lo,hi,eps,debug);
80  } else if (type==OT_F2G12) {*this=F2GFit(mu,lo,hi,eps,debug);
81  } else if (type==OT_BSH) {*this=BSHFit(mu,lo,hi,eps,debug);
82  } else {
83  print("Operator type not implemented: ",type);
84  MADNESS_EXCEPTION("Operator type not implemented: ",1);
85  }
86 
87  }
88 
89  /// copy constructor
90  GFit(const GFit& other) = default;
91 
92  /// assignment operator
93  GFit& operator=(const GFit& other) {
94  coeffs_ = other.coeffs_;
95  exponents_ = other.exponents_;
96  return *this;
97  }
98 
99  /// return a fit for the Coulomb function
100 
101  /// f(r) = 1/r
102  static GFit CoulombFit(double lo, double hi, double eps, bool prnt=false) {
103  GFit fit=BSHFit(0.0,lo,hi,eps/(4.0*constants::pi),prnt);
104  fit.coeffs_.scale(4.0*constants::pi); // BSHFit scales by 1/(4 pi), undo here
105  return fit;
106  }
107 
108  /// return a fit for the bound-state Helmholtz function
109 
110  /// the BSH function is defined by
111  /// f(r) = 1/ (4 pi) exp(-\mu r)/r
112  /// @param[in] mu the exponent of the BSH
113  /// @param[in] lo the smallest length scale that needs to be precisely represented
114  /// @param[in] hi the largest length scale that needs to be precisely represented
115  /// @param[in] eps the precision threshold
116  /// @parma[in] prnt print level
117  static GFit BSHFit(double mu, double lo, double hi, double eps, bool prnt=false) {
118  GFit fit;
119  bool fix_interval=false;
120  if (NDIM==3) bsh_fit(mu,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt,fix_interval);
121  else bsh_fit_ndim(NDIM,mu,lo,hi,eps,fit.coeffs_,fit.exponents_,prnt);
122 
123  if (prnt) {
124  print("bsh fit");
125  auto exact = [&mu](const double r) -> double { return 1.0/(4.0 * constants::pi) * exp(-mu * r)/r; };
126  fit.print_accuracy(exact, lo, hi);
127  }
128  return fit;
129  }
130 
131  /// return a fit for the Slater function
132 
133  /// the Slater function is defined by
134  /// f(r) = exp(-\gamma r)
135  /// @param[in] gamma the exponent of the Slater function
136  /// @param[in] lo the smallest length scale that needs to be precisely represented
137  /// @param[in] hi the largest length scale that needs to be precisely represented
138  /// @param[in] eps the precision threshold
139  /// @parma[in] prnt print level
140  static GFit SlaterFit(double gamma, double lo, double hi, double eps, bool prnt=false) {
141  GFit fit;
142  slater_fit(gamma,lo,hi,eps,fit.coeffs_,fit.exponents_,false);
143  if (prnt) {
144  print("Slater fit");
145  auto exact = [&gamma](const double r) -> double { return exp(-gamma * r); };
146  fit.print_accuracy(exact, lo, hi);
147  }
148  return fit;
149  }
150 
151  /// return a (trivial) fit for a single Gauss function
152 
153  /// the Gauss function is defined by
154  /// f(r) = exp(-\gamma r^2)
155  /// @param[in] gamma the exponent of the Gauss function
156  /// @param[in] lo the smallest length scale that needs to be precisely represented
157  /// @param[in] hi the largest length scale that needs to be precisely represented
158  /// @param[in] eps the precision threshold
159  /// @parma[in] prnt print level
160  static GFit GaussFit(double gamma, double lo, double hi, double eps, bool prnt=false) {
161  GFit fit;
162  fit.coeffs_=Tensor<T>(1);
163  fit.exponents_=Tensor<double>(1);
164  fit.coeffs_=1.0;
165  fit.exponents_=gamma;
166  if (prnt) {
167  print("Gauss fit");
168  auto exact = [&gamma](const double r) -> double { return exp(-gamma * r*r); };
169  fit.print_accuracy(exact, lo, hi);
170  }
171  return fit;
172  }
173 
174  /// return a fit for the F12 correlation factor
175 
176  /// the Slater function is defined by
177  /// f(r) = 1/(2 gamma) * (1 - exp(-\gamma r))
178  /// @param[in] gamma the exponent of the Slater function
179  /// @param[in] lo the smallest length scale that needs to be precisely represented
180  /// @param[in] hi the largest length scale that needs to be precisely represented
181  /// @param[in] eps the precision threshold
182  /// @parma[in] prnt print level
183  static GFit F12Fit(double gamma, double lo, double hi, double eps, bool prnt=false) {
184  GFit fit;
185  f12_fit(gamma,lo*0.1,hi,eps*0.01,fit.coeffs_,fit.exponents_,false);
186  fit.coeffs_*=(0.5/gamma);
187  if (prnt) {
188  print("f12 fit");
189  auto exact=[&gamma](const double r) -> double {return 0.5/gamma*(1.0-exp(-gamma*r));};
190  fit.print_accuracy(exact,lo,hi);
191  }
192  return fit;
193  }
194 
195  /// return a fit for the F12^2 correlation factor
196 
197  /// the Slater function square is defined by
198  /// f(r) = [ 1/(2 gamma) * (1 - exp(-\gamma r)) ] ^2
199  /// @param[in] gamma the exponent of the Slater function
200  /// @param[in] lo the smallest length scale that needs to be precisely represented
201  /// @param[in] hi the largest length scale that needs to be precisely represented
202  /// @param[in] eps the precision threshold
203  /// @parma[in] prnt print level
204  static GFit F12sqFit(double gamma, double lo, double hi, double eps, bool prnt=false) {
205  GFit fit;
206  f12sq_fit(gamma,lo*0.1,hi,eps*0.01,fit.coeffs_,fit.exponents_,false);
207  fit.coeffs_*=(0.25/(gamma*gamma));
208  if (prnt) {
209  print("f12sq fit");
210  auto exact=[&gamma](const double r) -> double {return std::pow(0.5/gamma*(1.0-exp(-gamma*r)),2.0);};
211  fit.print_accuracy(exact,lo,hi);
212  }
213  return fit;
214  }
215 
216 
217  /// return a fit for the FG function
218 
219  /// fg = 1/(2 mu) * (1 - exp(-gamma r12)) / r12
220  /// = 1/(2 mu) *( 1/r12 - exp(-gamma r12)/r12)
221  /// = 1/(2 mu) * (coulomb - bsh)
222  /// @param[in] gamma the exponent of the Slater function
223  /// @param[in] lo the smallest length scale that needs to be precisely represented
224  /// @param[in] hi the largest length scale that needs to be precisely represented
225  /// @param[in] eps the precision threshold
226  /// @parma[in] prnt print level
227  static GFit FGFit(double gamma, double lo, double hi, double eps, bool prnt=false) {
228  GFit bshfit,coulombfit;
229  eps*=0.1;
230  lo*=0.1;
231 // bool restrict_interval=false;
232  bool fix_interval=true;
233  bsh_fit(gamma,lo,hi,eps,bshfit.coeffs_,bshfit.exponents_,false,fix_interval);
234  bsh_fit(0.0,lo,hi,eps,coulombfit.coeffs_,coulombfit.exponents_,false,fix_interval);
235  // check the exponents are identical
236  auto diffexponents=(coulombfit.exponents() - bshfit.exponents());
237  MADNESS_CHECK(diffexponents.normf()/coulombfit.exponents().size()<1.e-12);
238  auto diffcoefficients=(coulombfit.coeffs() - bshfit.coeffs());
239  GFit fgfit;
240  fgfit.exponents_=bshfit.exponents_;
241  fgfit.coeffs_=4.0*constants::pi*0.5/gamma*diffcoefficients;
243 
244  if (prnt) {
245  print("fg fit");
246  auto exact=[&gamma](const double r) -> double {return 0.5/gamma*(1.0-exp(-gamma*r))/r;};
247  fgfit.print_accuracy(exact,lo,hi);
248  }
249  return fgfit;
250  }
251 
252  /// return a fit for the F2G function
253 
254  /// f2g = (1/(2 mu) * (1 - exp(-gamma r12)))^2 / r12
255  /// = 1/(4 mu^2) * [ 1/r12 - 2 exp(-gamma r12)/r12) + exp(-2 gamma r12)/r12 ]
256  /// @param[in] gamma the exponent of the Slater function
257  /// @param[in] lo the smallest length scale that needs to be precisely represented
258  /// @param[in] hi the largest length scale that needs to be precisely represented
259  /// @param[in] eps the precision threshold
260  /// @parma[in] prnt print level
261  static GFit F2GFit(double gamma, double lo, double hi, double eps, bool prnt=false) {
262  GFit bshfit,coulombfit,bsh2fit;
263  eps*=0.1;
264  lo*=0.1;
265 // bool restrict_interval=false;
266  bool fix_interval=true;
267  bsh_fit(gamma,lo,hi,eps,bshfit.coeffs_,bshfit.exponents_,false,fix_interval);
268  bsh_fit(2.0*gamma,lo,hi,eps,bsh2fit.coeffs_,bsh2fit.exponents_,false,fix_interval);
269  bsh_fit(0.0,lo,hi,eps,coulombfit.coeffs_,coulombfit.exponents_,false,fix_interval);
270 
271  // check the exponents are identical
272  auto diffexponents=(coulombfit.exponents() - bshfit.exponents());
273  MADNESS_CHECK(diffexponents.normf()/coulombfit.exponents().size()<1.e-12);
274  auto diffexponents1=(coulombfit.exponents() - bsh2fit.exponents());
275  MADNESS_CHECK(diffexponents1.normf()/coulombfit.exponents().size()<1.e-12);
276 
277  auto coefficients=(coulombfit.coeffs() - 2.0* bshfit.coeffs() + bsh2fit.coeffs());
278  GFit f2gfit;
279  f2gfit.exponents_=bshfit.exponents_;
280  // additional factor 4 pi due to implementation of bsh_fit
281  double fourpi=4.0*constants::pi;
282  double fourmu2=4.0*gamma*gamma;
283  f2gfit.coeffs_=fourpi/fourmu2*coefficients;
285 
286  if (prnt) {
287  print("fg fit");
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;
290  };
291  f2gfit.print_accuracy(exact,lo,hi);
292  }
293  return f2gfit;
294  }
295 
296 
297  /// return a fit for a general isotropic function
298 
299  /// note that the error is controlled over a uniform grid, the boundaries
300  /// will be poorly represented in general. Following Beylkin 2005
301  static GFit GeneralFit() {
302  MADNESS_EXCEPTION("General GFit still to be implemented",1);
303  return GFit();
304  }
305 
306  /// return the coefficients of the fit
307  Tensor<T> coeffs() const {return coeffs_;}
308 
309  /// return the exponents of the fit
310  Tensor<T> exponents() const {return exponents_;}
311 
312  void static prune_small_coefficients(const double eps, const double lo, const double hi,
313  Tensor<double>& coeff, Tensor<double>& expnt) {
314  double mid = lo + (hi-lo)*0.5;
315  long npt=coeff.size();
316  long i;
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);
323  if (std::max(std::abs(errlo),std::abs(errhi)) > 0.03*eps) break;
324  npt--;
325  coeff[i-1] = coeff[i-1] + cnew;
326  }
327  coeff = coeff(Slice(0,npt-1));
328  expnt = expnt(Slice(0,npt-1));
329  }
330 
332  double L, bool discardG0) const {
333  double tcut = 0.25/L/L;
334 
335  if (discardG0) {
336  // Relies on expnts being in decreasing order
337  for (int i=0; i<e.dim(0); ++i) {
338  if (e(i) < tcut) {
339  c = c(Slice(0,i));
340  e = e(Slice(0,i));
341  break;
342  }
343  }
344  } else {
345 // // Relies on expnts being in decreasing order
346 // int icut = -1;
347 // for (int i=0; i<e.dim(0); ++i) {
348 // if (e(i) < tcut) {
349 // icut = i;
350 // break;
351 // }
352 // }
353 // if (icut > 0) {
354 // for (int i=icut+1; i<e.dim(0); ++i) {
355 // c(icut) += c(i);
356 // }
357 // c = c(Slice(0,icut));
358 // e = e(Slice(0,icut));
359 // }
360  }
361  }
362 
363 private:
364 
365  /// ctor taking an isotropic function
366 
367  /// the function will be represented with a uniform error on a uniform grid
368  /// @param[in] f a 1d-function that implements T operator()
369  template<typename funcT>
370  GFit(const funcT& f) {
371 
372  }
373 
374  /// print coefficients and exponents, and values and errors
375 
376  /// @param[in] op the exact function, e.g. 1/r, exp(-mu r), etc
377  /// @param[in] lo lower bound for the range r
378  /// @param[in] hi higher bound for the range r
379  template<typename opT>
380  void print_accuracy(opT op, const double lo, const double hi) const {
381 
382  std::cout << "weights and roots" << std::endl;
383  for (int i=0; i<coeffs_.size(); ++i)
384  std::cout << i << " " << coeffs_[i] << " " << exponents_[i] << std::endl;
385 
386  long npt = 300;
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));
392 // double exact = exp(-mu*r)/r/4.0/constants::pi;
393  double exact = op(r);
394  double test = 0.0;
395  for (int j=0; j<coeffs_.dim(0); ++j)
396  test += coeffs_[j]*exp(-r*r*exponents_[j]);
397  double err = 0.0;
398  if (exact) err = (exact-test)/exact;
399  printf(" %.6e %8.1e %8.1e %8.1e\n",r, exact, exact-test, err);
400  }
401  }
402 
403  /// the coefficients of the expansion f(x) = \sum_m coeffs[m] exp(-exponents[m] * x^2)
405 
406  /// the exponents of the expansion f(x) = \sum_m coeffs[m] exp(-exponents[m] * x^2)
408 
409  /// fit the function exp(-mu r)/r
410 
411  /// formulas taken from
412  /// G. Beylkin and L. Monzon, On approximation of functions by exponential sums,
413  /// Appl Comput Harmon A, vol. 19, no. 1, pp. 17-48, Jul. 2005.
414  /// and
415  /// R. J. Harrison, G. I. Fann, T. Yanai, and G. Beylkin,
416  /// Multiresolution Quantum Chemistry in Multiwavelet Bases,
417  /// Lecture Notes in Computer Science, vol. 2660, p. 707, 2003.
418  static void bsh_fit(double mu, double lo, double hi, double eps,
419  Tensor<double>& pcoeff, Tensor<double>& pexpnt, bool prnt, bool fix_interval) {
420 
421  if (mu < 0.0) throw "cannot handle negative mu in bsh_fit";
422 // bool restrict_interval=(mu>0) and use_mu_for_restricting_interval;
423 
424 
425  if ((mu > 0) and (not fix_interval)) {
426 // if (restrict_interval) {
427  // Restrict hi according to the exponential decay
428  double r = -log(4*constants::pi*0.01*eps);
429  r = -log(r * 4*constants::pi*0.01*eps);
430  if (hi > r) hi = r;
431  }
432 
433  double TT;
434  double slo, shi;
435 
436  if (eps >= 1e-2) TT = 5;
437  else if (eps >= 1e-4) TT = 10;
438  else if (eps >= 1e-6) TT = 14;
439  else if (eps >= 1e-8) TT = 18;
440  else if (eps >= 1e-10) TT = 22;
441  else if (eps >= 1e-12) TT = 26;
442  else TT = 30;
443 
444  if ((mu > 0) and (not fix_interval)) {
445 // if (restrict_interval) {
446  slo = -0.5*log(4.0*TT/(mu*mu));
447  }
448  else {
449  slo = log(eps/hi) - 1.0;
450  }
451  shi = 0.5*log(TT/(lo*lo));
452  if (shi <= slo) throw "bsh_fit: logic error in slo,shi";
453 
454  // Resolution required for quadrature over s
455  double h = 1.0/(0.2-.50*log10(eps)); // was 0.5 was 0.47
456 
457  // Truncate the number of binary digits in h's mantissa
458  // so that rounding does not occur when performing
459  // manipulations to determine the quadrature points and
460  // to limit the number of distinct values in case of
461  // multiple precisions being used at the same time.
462  h = floor(64.0*h)/64.0;
463 
464 
465  // Round shi/lo up/down to an integral multiple of quadrature points
466  shi = ceil(shi/h)*h;
467  slo = floor(slo/h)*h;
468 
469  long npt = long((shi-slo)/h+0.5);
470 
471  //if (prnt)
472  //std::cout << "mu " << mu << " slo " << slo << " shi " << shi << " npt " << npt << " h " << h << " " << eps << std::endl;
473 
474  Tensor<double> coeff(npt), expnt(npt);
475 
476  for (int i=0; i<npt; ++i) {
477  double s = slo + h*(npt-i); // i+1
478  coeff[i] = h*2.0/sqrt(constants::pi)*exp(-mu*mu*exp(-2.0*s)/4.0)*exp(s);
479  coeff[i] = coeff[i]/(4.0*constants::pi);
480  expnt[i] = exp(2.0*s);
481  }
482 
483 #if ONE_TERM
484  npt=1;
485  double s=1.0;
486  coeff[0]=1.0;
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]);
491 
492 
493 #endif
494 
495  // Prune large exponents from the fit ... never necessary due to construction
496 
497  // Prune small exponents from Coulomb fit. Evaluate a gaussian at
498  // the range midpoint, and replace it there with the next most
499  // diffuse gaussian. Then examine the resulting error at the two
500  // end points ... if this error is less than the desired
501  // precision, can discard the diffuse gaussian.
502 
503  if ((mu == 0.0) and (not fix_interval)) {
504 // if (restrict_interval) {
505  GFit<double,3>::prune_small_coefficients(eps,lo,hi,coeff,expnt);
506  }
507 
508  // Modify the coeffs of the largest exponents to satisfy the moment conditions
509  //
510  // SETTING NMOM>1 TURNS OUT TO BE A BAD IDEA (AS CURRENTLY IMPLEMENTED)
511  // [It is accurate and efficient for a one-shot application but it seems to
512  // introduce fine-scale noise that amplifies during iterative solution of
513  // the SCF and DFT equations ... the symptom is negative coeffs in the fit]
514  //
515  // SET NMOM=0 or 1 (1 recommended) unless you are doing a one-shot application
516  //
517  // Determine the effective range of the four largest exponents and compute
518  // moments of the exact and remainder of the fit. Then adjust the coeffs
519  // to reproduce the exact moments in that volume.
520  //
521  // If nmom!=4 we assume that we will eliminate n=-1 which is stored first
522  // in the moment list
523  //
524  // <r^i|gj> cj = <r^i|exact-remainder>
525  const long nmom = 1;
526  if (nmom > 0) {
527  Tensor<double> q(4), qg(4);
528  double range = sqrt(-log(1e-6)/expnt[nmom-1]);
529  if (prnt) print("exponent(nmom-1)",expnt[nmom-1],"has range", range);
530 
531  bsh_spherical_moments(mu, range, q);
532  Tensor<double> M(nmom,nmom);
533  for (int i=nmom; i<npt; ++i) {
534  Tensor<double> qt(4);
535  gaussian_spherical_moments(expnt[i], range, qt);
536  qg += qt*coeff[i];
537  }
538  if (nmom != 4) {
539  q = q(Slice(1,nmom));
540  qg = qg(Slice(1,nmom));
541  }
542  if (prnt) {
543  print("moments", q);
544  print("moments", qg);
545  }
546  q = q - qg;
547  for (int j=0; j<nmom; ++j) {
548  Tensor<double> qt(4);
549  gaussian_spherical_moments(expnt[j], range, qt);
550  if (nmom != 4) qt = qt(Slice(1,nmom));
551  for (int i=0; i<nmom; ++i) {
552  M(i,j) = qt[i];
553  }
554  }
555  Tensor<double> ncoeff;
556  gesv(M, q, ncoeff);
557  if (prnt) {
558  print("M\n",M);
559  print("old coeffs", coeff(Slice(0,nmom-1)));
560  print("new coeffs", ncoeff);
561  }
562 
563  coeff(Slice(0,nmom-1)) = ncoeff;
564  }
565 
566  pcoeff = coeff;
567  pexpnt = expnt;
568  }
569 
570  /// fit a Slater function using a sum of Gaussians
571 
572  /// formula inspired by the BSH fit, with the roles of r and mu exchanged
573  /// see also Eq. (A3) in
574  /// S. Ten-no, Initiation of explicitly correlated Slater-type geminal theory,
575  /// Chem. Phys. Lett., vol. 398, no. 1, pp. 56-61, 2004.
576  static void slater_fit(double gamma, double lo, double hi, double eps,
577  Tensor<double>& pcoeff, Tensor<double>& pexpnt, bool prnt) {
578 
579  MADNESS_CHECK(gamma >0.0);
580  // empirical number TT for the upper integration limit
581  double TT;
582  if (eps >= 1e-2) TT = 5;
583  else if (eps >= 1e-4) TT = 10;
584  else if (eps >= 1e-6) TT = 14;
585  else if (eps >= 1e-8) TT = 18;
586  else if (eps >= 1e-10) TT = 22;
587  else if (eps >= 1e-12) TT = 26;
588  else TT = 30;
589 
590  // integration limits for quadrature over s: slo and shi
591  // slo and shi must not depend on gamma!!!
592  double slo=0.5 * log(eps) - 1.0;
593  double shi=log(TT/(lo*lo))*0.5;
594 
595  // Resolution required for quadrature over s
596  double h = 1.0/(0.2-.5*log10(eps)); // was 0.5 was 0.47
597 
598  // Truncate the number of binary digits in h's mantissa
599  // so that rounding does not occur when performing
600  // manipulations to determine the quadrature points and
601  // to limit the number of distinct values in case of
602  // multiple precisions being used at the same time.
603  h = floor(64.0*h)/64.0;
604 
605  // Round shi/lo up/down to an integral multiple of quadrature points
606  shi = ceil(shi/h)*h;
607  slo = floor(slo/h)*h;
608 
609  long npt = long((shi-slo)/h+0.5);
610 
611  Tensor<double> coeff(npt), expnt(npt);
612 
613  for (int i=0; i<npt; ++i) {
614  const double s = slo + h*(npt-i); // i+1
615  coeff[i] = h*exp(-gamma*gamma*exp(2.0*s) + s);
616  coeff[i]*=2.0*gamma/sqrt(constants::pi);
617  expnt[i] = 0.25*exp(-2.0*s);
618  }
619 
620  if (prnt) {
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;
624 
625  long npt = 300;
626  //double hi = 1.0;
627  //if (mu) hi = min(1.0,30.0/mu);
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);
634  double test = 0.0;
635  for (int j=0; j<coeff.dim(0); ++j)
636  test += coeff[j]*exp(-r*r*expnt[j]);
637  double err = 0.0;
638  if (exact) err = (exact-test)/exact;
639  printf(" %.6e %8.1e %8.1e %8.1e\n",r, exact, exact-test, err);
640  }
641  }
642  pcoeff = coeff;
643  pexpnt = expnt;
644  }
645 
646 
647  /// fit a correlation factor (1- exp(-mu r))
648 
649  /// use the Slater fit with an additional term: 1*exp(-0 r^2)
650  static void f12_fit(double gamma, double lo, double hi, double eps,
651  Tensor<double>& pcoeff, Tensor<double>& pexpnt, bool prnt) {
652  Tensor<double> coeff,expnt;
653  slater_fit(gamma, lo, hi, eps, coeff, expnt, prnt);
654 
655  pcoeff=Tensor<double>(coeff.size()+1);
656  pcoeff(Slice(1,-1,1))=-coeff(_);
657  pexpnt=Tensor<double>(expnt.size()+1);
658  pexpnt(Slice(1,-1,1))=expnt(_);
659 
660  pcoeff(0l)=1.0;
661  pexpnt(0l)=1.e-10;
662  }
663 
664 
665  /// fit a correlation factor f12^2 = (1- exp(-mu r))^2 = 1 - 2 exp(-mu r) + exp(-2 mu r)
666 
667  /// no factor 1/(2 mu) or square of it included!
668  /// use the Slater fit with an additional term: 1*exp(-0 r^2)
669  static void f12sq_fit(double gamma, double lo, double hi, double eps,
670  Tensor<double>& pcoeff, Tensor<double>& pexpnt, bool prnt) {
671  Tensor<double> coeff1,expnt1, coeff2, expnt2;
672  slater_fit(gamma, lo, hi, eps, coeff1, expnt1, prnt);
673  slater_fit(2.0*gamma, lo, hi, eps, coeff2, expnt2, prnt);
674 
675  // check exponents are the same
676  MADNESS_CHECK((expnt1-expnt2).normf()/expnt1.size()<1.e-12);
677 
678  // add exponential terms
679  auto coeff=coeff2-2.0*coeff1;
680  pcoeff=Tensor<double>(coeff1.size()+1);
681  pcoeff(Slice(1,-1,1))=coeff(_);
682  pexpnt=Tensor<double>(expnt1.size()+1);
683  pexpnt(Slice(1,-1,1))=expnt1(_);
684 
685  // add constant term
686  pcoeff(0l)=1.0;
687  pexpnt(0l)=1.e-10;
688  }
689 
690 
691  void static bsh_fit_ndim(int ndim, double mu, double lo, double hi, double eps,
692  Tensor<double>& pcoeff, Tensor<double>& pexpnt, bool prnt) {
693 
694  if (mu > 0) {
695  // Restrict hi according to the exponential decay
696  double r = -log(4*constants::pi*0.01*eps);
697  r = -log(r * 4*constants::pi*0.01*eps);
698  if (hi > r) hi = r;
699  }
700 
701 
702  // Determine range of quadrature by estimating when
703  // kernel drops to exp(-100)
704 
705  double slo, shi;
706  if (mu > 0) {
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));
709  }
710  else {
711  slo = log(eps/hi) - 1.0;
712  }
713  shi = 0.5*log(100.0/(lo*lo));
714 
715  // Resolution required for quadrature over s
716  double h = 1.0/(0.2-.50*log10(eps)); // was 0.5 was 0.47
717 
718  // Truncate the number of binary digits in h's mantissa
719  // so that rounding does not occur when performing
720  // manipulations to determine the quadrature points and
721  // to limit the number of distinct values in case of
722  // multiple precisions being used at the same time.
723  h = floor(64.0*h)/64.0;
724 
725 
726  // Round shi/lo up/down to an integral multiple of quadrature points
727  shi = ceil(shi/h)*h;
728  slo = floor(slo/h)*h;
729 
730  long npt = long((shi-slo)/h+0.5);
731 
732  if (prnt)
733  std::cout << "bsh: mu " << mu << " lo " << lo << " hi " << hi
734  << " eps " << eps << " slo " << slo << " shi " << shi
735  << " npt " << npt << " h " << h << std::endl;
736 
737 
738  // Compute expansion pruning small coeffs and large exponents
739  Tensor<double> coeff(npt), expnt(npt);
740  int nnpt=0;
741  for (int i=0; i<npt; ++i) {
742  double s = slo + h*(npt-i); // i+1
743  double c = exp(-0.25*mu*mu*exp(-2.0*s)+(ndim-2)*s)*0.5/pow(constants::pi,0.5*ndim);
744  double p = exp(2.0*s);
745  c = c*h;
746  if (c*exp(-p*lo*lo) > eps) {
747  coeff(nnpt) = c;
748  expnt(nnpt) = p;
749  ++nnpt;
750  }
751  }
752  npt = nnpt;
753 #if ONE_TERM
754  npt=1;
755  double s=1.0;
756  coeff[0]=1.0;
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]);
761 
762 #endif
763 
764 
765  // Prune small exponents from Coulomb fit. Evaluate a gaussian at
766  // the range midpoint, and replace it there with the next most
767  // diffuse gaussian. Then examine the resulting error at the two
768  // end points ... if this error is less than the desired
769  // precision, can discard the diffuse gaussian.
770 
771  if (mu == 0.0) {
772  double mid = lo + (hi-lo)*0.5;
773  long i;
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);
780  if (std::max(std::abs(errlo),std::abs(errhi)) > 0.03*eps) break;
781  npt--;
782  coeff[i-1] = coeff[i-1] + cnew;
783  }
784  }
785 
786  // Shrink array to correct size
787  coeff = coeff(Slice(0,npt-1));
788  expnt = expnt(Slice(0,npt-1));
789 
790 
791  if (prnt) {
792  for (int i=0; i<npt; ++i)
793  std::cout << i << " " << coeff[i] << " " << expnt[i] << std::endl;
794 
795  long npt = 300;
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));
801  double test = 0.0;
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);
805  }
806  }
807 
808  pcoeff = coeff;
809  pexpnt = expnt;
810  }
811 
812  // Returns in q[0..4] int(r^2(n+1)*exp(-alpha*r^2),r=0..R) n=-1,0,1,2
813  static void gaussian_spherical_moments(double alpha, double R, Tensor<double>& q) {
814  q[0] = -(-0.1e1 + exp(-alpha * R*R)) / alpha / 0.2e1;
815  q[1] = (-0.2e1 * R * pow(alpha, 0.3e1 / 0.2e1) + sqrt(constants::pi)
816  * erf(R * sqrt(alpha)) * alpha * exp(alpha * R*R))
817  * pow(alpha, -0.5e1 / 0.2e1) * exp(-alpha * R*R) / 0.4e1;
818  q[2] = -(-0.1e1 + exp(-alpha * R*R) + exp(-alpha * R*R) * alpha * R*R)
819  * pow(alpha, -0.2e1) / 0.2e1;
820  q[3] = -(-0.3e1 * sqrt(constants::pi) * erf(R * sqrt(alpha)) * pow(alpha, 0.2e1)
821  * exp(alpha * R*R) + 0.6e1 * R * pow(alpha, 0.5e1 / 0.2e1)
822  + 0.4e1 * pow(R, 0.3e1) * pow(alpha, 0.7e1 / 0.2e1))
823  * pow(alpha, -0.9e1 / 0.2e1) * exp(-alpha * R*R) / 0.8e1;
824  }
825 
826  // Returns in q[0..4] int(r^2(n+1)*exp(-mu*r)/(4*constants::pi*r),r=0..R) n=-1,0,1,2
827  static void bsh_spherical_moments(double mu, double R, Tensor<double>& q) {
828  if (mu == 0.0) {
829  q[0] = R / constants::pi / 0.4e1;
830  q[1] = pow(R, 0.2e1) / constants::pi / 0.8e1;
831  q[2] = pow(R, 0.3e1) / constants::pi / 0.12e2;
832  q[3] = pow(R, 0.4e1) / constants::pi / 0.16e2;
833  }
834  else {
835  q[0] = (exp(mu * R) - 0.1e1) / mu * exp(-mu * R) / constants::pi / 0.4e1;
836  q[1] = -(-exp(mu * R) + 0.1e1 + mu * R) * pow(mu, -0.2e1) / constants::pi
837  * exp(-mu * R) / 0.4e1;
838  q[2] = -(-0.2e1 * exp(mu * R) + 0.2e1 + 0.2e1 * mu * R + R*R *
839  pow(mu, 0.2e1))*pow(mu, -0.3e1) / constants::pi * exp(-mu * R) / 0.4e1;
840  q[3] = -(-0.6e1 * exp(mu * R) + 0.6e1 + 0.6e1 * mu * R + 0.3e1 * R*R
841  * pow(mu, 0.2e1) + pow(R, 0.3e1) * pow(mu, 0.3e1))
842  * pow(mu, -0.4e1) / constants::pi * exp(-mu * R) / 0.4e1;
843  }
844  }
845 
846 };
847 
848 
849 }
850 
851 #endif // MADNESS_MRA_GFIT_H__INCLUDED
852 
double q(double t)
Definition: DKops.h:18
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
Definition: gfit.h:57
void print_accuracy(opT op, const double lo, const double hi) const
print coefficients and exponents, and values and errors
Definition: gfit.h:380
GFit & operator=(const GFit &other)
assignment operator
Definition: gfit.h:93
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
Tensor< T > exponents() const
return the exponents of the fit
Definition: gfit.h:310
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
Tensor< T > coeffs() const
return the coefficients of the fit
Definition: gfit.h:307
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
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
void test(World &world, bool doloadbal=false)
Definition: dataloadbal.cc:224
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 max(a, b)
Definition: lda.h:51
#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:190
#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:210
const double pi
Mathematical constant .
Definition: constants.h:48
File holds all helper structures necessary for the CC_Operator and CC2 class.
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
const double alpha
Definition: test_coulomb.cc:54
void e()
Definition: test_sig.cc:75
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