MADNESS  0.10.1
electronic_correlation_factor.h
Go to the documentation of this file.
1 /*
2  * electronic_correlation_factor.h
3  *
4  * Created on: Jul 9, 2015
5  * Author: fbischoff
6  */
7 
8 #ifndef SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
9 #define SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
10 
11 
12 
13 #include <madness/mra/mra.h>
14 #include <madness/mra/lbdeux.h>
16 #include <iomanip>
17 
18 namespace madness {
19 /// a class holding the electronic correlation factor for R12 theory
20 class CorrelationFactor {
21 
22  World& world;
23  double _gamma; ///< the correlation factor exp(-gamma r12)
24  double dcut; ///< the cutoff for the 1/r potential
25  double lo; ///< smallest length scale to be resolved
26 
27 public:
28 
29  /// ctor, use negative gamma for linear correlation factor r12
30  CorrelationFactor(World& world) : world(world), _gamma(-1.0), dcut(1.e-10),
31  lo(1.e-10) {
32  }
33 
34  /// ctor, use negative gamma for linear correlation factor r12
35  CorrelationFactor(World& world, const double& gamma, const double dcut,
36  const Molecule& molecule) : world(world), _gamma(gamma), dcut(dcut) {
37  lo=1.e-6;//lo = molecule.smallest_length_scale();
38 // if (world.rank()==0) {
39 // if (gamma>0.0) print("constructed correlation factor with gamma=",gamma);
40 // else if (gamma==0.0) print("constructed linear correlation factor");
41 // }
42  }
43  /// ctor, use negative gamma for linear correlation factor r12
44  CorrelationFactor(World& world, const double& gamma, const double dcut,
45  const double lo) : world(world), _gamma(gamma), dcut(dcut), lo(lo) {
46 // if (world.rank()==0) {
47 // if (gamma>0.0) print("constructed correlation factor with gamma=",gamma);
48 // else if (gamma==0.0) print("constructed linear correlation factor");
49 // }
50  }
51 
52  /// copy ctor
53  CorrelationFactor(const CorrelationFactor& other) : world(other.world) {
54  _gamma=other._gamma;
55  dcut=other.dcut;
56  lo=other.lo;
57  }
58 
59  /// assignment; assume other's world is this world
60  CorrelationFactor& operator=(const CorrelationFactor& other) {
61  _gamma=other._gamma;
62  dcut=other.dcut;
63  lo=other.lo;
64  return *this;
65  }
66 
67  /// return the exponent of this correlation factor
68  double gamma() const {return _gamma;}
69 
70  /// return the value of the correlation factor
71  double operator()(const coord_6d& r) const {
72  const double rr=r12(r);
73  if (_gamma>0.0) return (1.0-exp(-_gamma*rr))/(2.0*_gamma);
74  return 0.5*rr;
75  }
76 
77  /// apply Kutzelnigg's regularized potential to an orbital product
78  real_function_6d apply_U(const real_function_3d& phi_i, const real_function_3d& phi_j,
79  const real_convolution_6d& op_mod, const bool symmetric=false) const {
80  if(not op_mod.modified()) MADNESS_EXCEPTION("ElectronicCorrelationFactor::apply_U, op_mod must be in modified_NS form",1);
81  const double thresh = FunctionDefaults<6>::get_thresh();
82  const bool debug = false;
83  if(symmetric) MADNESS_ASSERT((phi_i-phi_j).norm2() < FunctionDefaults<3>::get_thresh());
84 
85  real_function_6d result=real_factory_6d(world);
86 
87  for (int axis=0; axis<3; ++axis) {
88  //if (world.rank()==0) print("working on axis",axis);
89  real_derivative_3d D = free_space_derivative<double,3>(world, axis);
90  const real_function_3d Di=(D(phi_i)).truncate();
92  if(symmetric) Dj=madness::copy(Di);
93  else Dj=(D(phi_j)).truncate();
94 
96  real_function_6d tmp1=CompositeFactory<double,6,3>(world)
97  .g12(u).particle1(copy(Di)).particle2(copy(phi_j)).thresh(thresh);
98  tmp1.fill_cuspy_tree(op_mod).truncate();
99 
100  real_function_6d tmp2;
101  if(symmetric) tmp2 = -1.0*swap_particles(tmp1);
102  else{
103  tmp2=CompositeFactory<double,6,3>(world)
104  .g12(u).particle1(copy(phi_i)).particle2(copy(Dj)).thresh(thresh);
105  tmp2.fill_cuspy_tree(op_mod).truncate();
106  }
107 
108  result=result+(tmp1-tmp2).truncate();
109 
110 
111  tmp1.clear();
112  tmp2.clear();
113  world.gop.fence();
114  result.truncate().reduce_rank();
115  }
116 
117  // include the purely local potential that (partially) cancels 1/r12
118  if (_gamma>0.0) {
119  fg_ func(_gamma,dcut);
121  real_function_6d mul=CompositeFactory<double,6,3>(world)
122  .g12(fg3).particle1(copy(phi_i)).particle2(copy(phi_j)).thresh(thresh);;
123  mul.fill_cuspy_tree(op_mod).truncate();
124  // mul.print_size("mul");
125 
126  result=(result+mul).truncate().reduce_rank();
127  }
128  if(debug) result.print_size("Ue|ij>");
129  return result;
130  }
131 
132  /// return the U1 term of the correlation function
133  real_function_6d U1(const int axis) const {
134  U func(_gamma,axis,dcut);
135  const real_function_6d u1=real_factory_6d(world)
137  return u1;
138  }
139 
140  /// return the U1 term of the correlation function
141  real_function_6d U2() const {
142  if (world.rank()==0) print("U2 for the electronic correlation factor");
143  if (world.rank()==0) print("is expensive -- do you really need it??");
144  MADNESS_EXCEPTION("U2() not implemented, since it might be expensive",1);
145  return real_factory_6d(world);
146  }
147 
148  /// return the correlation factor as on-demand function
149  real_function_6d f() const {
150 // real_function_6d tmp=real_factory_6d(world).functor2(*this).is_on_demand();
152  real_function_6d tmp=TwoElectronFactory(world)
153  .dcut(dcut).gamma(_gamma).f12().thresh(thresh);
154  return tmp;
155  }
156 
157  /// return f^2 as on-demand function
158  real_function_6d f2() const {
159  f2_ func(_gamma);
161  return tmp;
162  }
163 
164  /// return fg+sth as on-demand function
165  real_function_6d fg() const {
166  fg_ func(_gamma,dcut);
168  return tmp;
169  }
170 
171  /// return f/r as on-demand function
172  real_function_6d f_over_r() const {
173  f_over_r_ func(_gamma,dcut);
175  return tmp;
176  }
177 
178  /// return (\nabla f)^2 as on-demand functions
179  real_function_6d nablaf2() const {
180  nablaf2_ func(_gamma);
182  return tmp;
183  }
184 
185 private:
186  /// functor for the local potential (1-f12)/r12 + sth (doubly connected term of the commutator)
187 
188  /// TODO: turn this into coeffs directly
189  struct fg_ : FunctionFunctorInterface<double,6> {
190  double gamma;
191  double dcut;
192  fg_(double gamma, double dcut) : gamma(gamma), dcut(dcut) {
193  MADNESS_ASSERT(gamma>0.0);
194  }
195  double operator()(const coord_6d& r) const {
196  const double rr=r12(r);
197  const double e=exp(-gamma*rr);
198  return (1.0-e)*u(rr,dcut) + 0.5*gamma*e;
199  }
200  };
201 
202  /// functor for the local potential (1-f12)/r12
203  struct f_over_r_ : FunctionFunctorInterface<double,6> {
204  double gamma;
205  double dcut;
206  f_over_r_(double gamma, double dcut) : gamma(gamma), dcut(dcut) {
207  MADNESS_ASSERT(gamma>0.0);
208  }
209  double operator()(const coord_6d& r) const {
210  const double rr=r12(r);
211  const double e=exp(-gamma*rr);
212  return (1.0-e)*u(rr,dcut)/(2.0*gamma);
213  }
214  };
215 
216  /// functor for the local part of the regularized potential: f12/r12*(r1-r2)(D1-D2)
217  struct U : FunctionFunctorInterface<double,6> {
218  double gamma;
219  int axis;
220  double dcut;
221  U(double gamma, int axis, double dcut) : gamma(gamma), axis(axis),
222  dcut(dcut) {
223  MADNESS_ASSERT(axis>=0 and axis<3);
224  }
225  double operator()(const coord_6d& r) const {
226  const double rr=r12(r);
227  const coord_3d vr12{r[0]-r[3],r[1]-r[4],r[2]-r[5]};
228  const coord_3d N=unitvec(vr12);
229  if (gamma>0.0) return -0.5*exp(-gamma*rr)*N[axis];
230  MADNESS_EXCEPTION("no gamma in electronic corrfac::U1",1);
231 // const double rr=r12(r);
232 // const double g12=u(rr,dcut);
233 // double a=0.5;
234 // if (gamma>0.0) a=0.5*exp(-gamma*rr);
235 // return -a*x12(r,axis) * g12;
236  }
237  };
238 
239  /// functor for the local potential (1-f12)^2
240  struct f2_ : FunctionFunctorInterface<double,6> {
241  double gamma;
242  f2_(double gamma) : gamma(gamma) {MADNESS_ASSERT(gamma>0.0);}
243  double operator()(const coord_6d& r) const {
244  const double rr=r12(r);
245  const double e=exp(-gamma*rr);
246  const double f=(1.0-e)/(2.0*gamma);
247  return f*f;
248  }
249  };
250 
251  /// functor for the local potential (\nabla f)^2
252  struct nablaf2_ : FunctionFunctorInterface<double,6> {
253  double gamma;
254  nablaf2_(double gamma) : gamma(gamma) {
255  MADNESS_ASSERT(gamma>0.0);
256  MADNESS_ASSERT(gamma==1.0);
257  }
258  double operator()(const coord_6d& r) const {
259  const double rr=r12(r);
260  const double f=exp(-2.0*gamma*rr)/(4.0*gamma*gamma);
261  return f;
262  }
263  };
264 
265  /// Smoothed 1/r potential (c is the smoothing distance)
266  static double u(double r, double c) {
267  r = r/c;
268  double r2 = r*r, pot;
269  if (r > 6.5){
270  pot = 1.0/r;
271  } else if (r > 1e-2) {
272  pot = erf(r)/r + exp(-r2)*0.56418958354775630;
273  } else{
274  pot = 1.6925687506432689-r2*(0.94031597257959381-r2*(0.39493270848342941-0.12089776790309064*r2));
275  }
276  return pot/c;
277  }
278 
279  static double r12(const coord_6d& r) {
280  const double x12=r[0]-r[3];
281  const double y12=r[1]-r[4];
282  const double z12=r[2]-r[5];
283  const double r12=sqrt(x12*x12 + y12*y12 + z12*z12);
284  return r12;
285  }
286  static double x12(const coord_6d& r, const int axis) {
287  return r[axis]-r[axis+3];
288  }
289 
290  static coord_3d smoothed_unitvec(const coord_3d& xyz, double smoothing) {
291 // if (smoothing==0.0) smoothing=molecule.get_eprec();
292  // TODO:need to test this
293  // reduce the smoothing for the unitvector
294  //if (not (this->type()==None or this->type()==Two)) smoothing=sqrt(smoothing);
295  const double r=xyz.normf();
296  const double cutoff=smoothing;
297  if (r>cutoff) {
298  return 1.0/r*xyz;
299  } else {
300  const double xi=r/cutoff;
301  const double xi2=xi*xi;
302  const double xi3=xi*xi*xi;
303 // const double nu21=0.5+1./32.*(45.*xi - 50.*xi3 + 21.*xi*xi*xi*xi*xi);
304  const double nu22=0.5 + 1./64.*(105* xi - 175 *xi3 + 147* xi2*xi3 - 45* xi3*xi3*xi);
305 // const double nu40=0.5 + 1./128.*(225 *xi - 350 *xi3 + 189*xi2*xi3);
306  const double kk=2.*nu22-1.0;
307  return kk/(r+1.e-15)*xyz;
308  }
309  }
310 };
311 
312 /// a class holding the electronic correlation factor for R12 theory
313 /// CorrelationFactor2 = (1-0.5*exp(-gamma*r12), gamma=0.5
314 /// (CorrelationFactor + 1)*2.0 = CorrelationFactor2 (currently CorrelationFactor2 is only implemented for gamma=0.5 so use this gamma also on CorrelationFactor
315 class CorrelationFactor2 {
316 
317  World& world;
318  double _gamma; ///< the correlation factor exp(-gamma r12)
319  typedef std::shared_ptr< FunctionFunctorInterface<double,6> > functorT;
320 
321 public:
322 
323  double dcut; ///< the cutoff for the 1/r potential
324  double lo; ///< smallest length scale to be resolved
325  double vtol; ///< initial projection threshold
326 
327 
328  /// ctor, use negative gamma for linear correlation factor r12
329  CorrelationFactor2(World& world) : world(world), _gamma(0.5), dcut(1.e-10),
330  lo(1.e-10), vtol(FunctionDefaults<3>::get_thresh()*0.1) {
331  MADNESS_ASSERT(_gamma==0.5);
332  }
333 
334  /// return the exponent of this correlation factor
335  double gamma() const {return _gamma;}
336 
337  real_function_6d function() const {
338  functorT R=functorT(new R_functor(_gamma,1));
339  return real_factory_6d(world).functor(R).is_on_demand();
340  }
341 
342  real_function_6d square() const {
343  functorT R2=functorT(new R_functor(_gamma,2));
344  return real_factory_6d(world).functor(R2).is_on_demand();
345  }
346 
347  real_function_6d inverse() const {
348  functorT R=functorT(new R_functor(_gamma,-1));
349  return real_factory_6d(world).functor(R).is_on_demand();
350  }
351 
352  /// return the U1 term of the correlation function
353  real_function_6d U1(const int axis) const {
354  functorT U1f=functorT(new U1_functor(_gamma,axis));
355  return real_factory_6d(world).functor(U1f).is_on_demand();
356  }
357 
358  /// return the U2 term of the correlation function
359  real_function_6d U2() const {
360  functorT U2f=functorT(new U2_functor(_gamma));
361  return real_factory_6d(world).functor(U2f).is_on_demand();
362  }
363 
364  /// apply Kutzelnigg's regularized potential to an orbital product
365  real_function_6d apply_U(const real_function_6d& psi, const double eps) const {
366  const double bsh_thresh=1.e-7;
367 
368  real_function_6d result=real_factory_6d(world);
369 
370  real_convolution_6d op_mod = BSHOperator<6>(world, sqrt(-2*eps), lo,bsh_thresh);
371  op_mod.modified()=true;
372 
373  for (int axis=0; axis<3; ++axis) {
374  //if (world.rank()==0) print("working on axis",axis);
375  real_derivative_6d D1 = free_space_derivative<double,6>(world, axis);
376  real_derivative_6d D2 = free_space_derivative<double,6>(world, axis+3);
377  const real_function_6d Drhs1=D1(psi).truncate();
378  const real_function_6d Drhs2=D2(psi).truncate();
379 
380  const real_function_6d u1=U1(axis);
381 
382  real_function_6d tmp1=CompositeFactory<double,6,3>(world)
383  .g12(u1).ket(copy(Drhs1));
384  tmp1.fill_cuspy_tree(op_mod).truncate();
385 
386  real_function_6d tmp2=CompositeFactory<double,6,3>(world)
387  .g12(u1).ket(copy(Drhs2));
388  tmp2.fill_cuspy_tree(op_mod).truncate();
389  // if (world.rank()==0) print("done with fill_tree");
390 
391  result=result+(tmp1-tmp2).truncate();
392  tmp1.clear();
393  tmp2.clear();
394  world.gop.fence();
395  result.truncate().reduce_rank();
396 
397  // if (world.rank()==0) printf("done with multiplication with U at ime %.1f\n",wall_time());
398  // result.print_size("result");
399  }
400 
401  real_function_6d u2=U2();
402  real_function_6d r2=CompositeFactory<double,6,3>(world).ket(copy(psi))
403  .g12(u2);
404  r2.fill_tree(op_mod);
405  result=(result+r2).truncate();
406  return result;
407  }
408 
409 
410 private:
411 
412  /// functor for the correlation factor R
413  class R_functor : public FunctionFunctorInterface<double,6> {
414  double gamma;
415  int exponent;
416 
417 
418  public:
419  R_functor(double gamma, int e=1) : gamma(gamma), exponent(e) {
420  MADNESS_ASSERT(gamma==0.5);
421  }
422 
423  // only valid for gamma=1
424  double operator()(const coord_6d& r) const {
425  const double rr=r12(r);
426  double val=(1.0-0.5*exp(-gamma*rr));
427  if (exponent==1) return val;
428  else if (exponent==2) return val*val;
429  else if (exponent==-1) return 1.0/val;
430  else {
431  MADNESS_EXCEPTION("fancy exponent in correlationfactor2",1);
432  }
433  }
434  };
435 
436  /// functor for the U2 local potential
437  class U2_functor : public FunctionFunctorInterface<double,6> {
438  double gamma;
439 
440  public:
441  U2_functor(double gamma) : gamma(gamma) {
442  MADNESS_ASSERT(gamma==0.5);
443  }
444 
445  // only valid for gamma=1
446  double operator()(const coord_6d& r) const {
447  const double rr=r12(r);
448  // Taylor expansion for small r
449  if (rr<1.e-4) { // valid for gamma==0.5, otherwise singular
450  return (5./4.0 - rr + (35.0* rr*rr)/48.0 - (101.0*rr*rr*rr)/192.0);
451  }
452  const double egr=exp(-gamma*rr);
453  return -(-8.*egr + 8.0 + rr*egr)/(4.0 *rr*egr - 8 *rr);
454  }
455  };
456 
457  /// functor for the U1 = -\frac{\vec\nabla_1 f_{12}}{f_{12}} potential
458 
459  /// the potential is given by
460  /// U1 = -\frac{\vec\nabla_1 f_{12}}{f_{12}}
461  /// = \frac{e^{-r12/2}{4-2e^{-r12/2}} \vec unitvec
462  /// the derivative operators are not included
463  class U1_functor : public FunctionFunctorInterface<double,6> {
464  double gamma;
465  int axis;
466 
467  public:
468  U1_functor(double gamma, int axis) : gamma(gamma), axis(axis) {
469  MADNESS_ASSERT(gamma==0.5);
470  MADNESS_ASSERT(axis<3);
471  }
472 
473  double operator()(const coord_6d& r) const {
474  const double rr=r12(r);
475  const coord_3d vr12{r[0]-r[3],r[1]-r[4],r[2]-r[5]};
476  const coord_3d N=unitvec(vr12);
477  // Taylor expansion for small r
478  double val;
479  if (rr<1.e-4) { // valid for gamma==0.5, otherwise singular
480  val = 0.5 - 0.5*rr + 0.125*(3.*rr*rr) - (13.* rr*rr*rr)/48.0;
481  } else {
482  const double egr=exp(-gamma*rr);
483  val=egr/(4.0-2.0*egr);
484  }
485  // NOTE the sign
486  return -val*N[axis];
487  }
488  };
489 
490  /// helper function
491  static double r12(const coord_6d& r) {
492  const double x12=r[0]-r[3];
493  const double y12=r[1]-r[4];
494  const double z12=r[2]-r[5];
495  const double r12=sqrt(x12*x12 + y12*y12 + z12*z12);
496  return r12;
497  }
498 
499 };
500 
501 }
502 
503 
504 
505 #endif /* SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_ */
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
virtual FunctionFactory & is_on_demand()
Definition: function_factory.h:281
FunctionFactory & functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > &f)
Definition: function_factory.h:141
Function< T, NDIM > & fill_cuspy_tree(const opT &op, const bool fence=true)
Definition: mra.h:1166
double thresh() const
Returns value of truncation threshold. No communication.
Definition: mra.h:567
Function< T, NDIM > & truncate(double tol=0.0, bool fence=true)
Truncate the function with optional fence. Compresses with fence if not compressed.
Definition: mra.h:602
T normf() const
Calculate the 2-norm of the vector elements.
Definition: vector.h:400
static const double R
Definition: csqrt.cc:46
double(* f2)(const coord_3d &)
Definition: derivatives.cc:56
static double lo
Definition: dirac-hatom.cc:23
static bool debug
Definition: dirac-hatom.cc:16
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
static const double dcut
Definition: he.cc:13
Implements (2nd generation) static load/data balancing for functions.
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition: mcpfit.cc:49
Main include file for MADNESS and defines Function interface.
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
std::enable_if_t< NDIM%2==0, Function< T, NDIM > > swap_particles(const Function< T, NDIM > &f)
swap particles 1 and 2
Definition: mra.h:2302
Function< TENSOR_RESULT_TYPE(Q, T), NDIM > mul(const Q alpha, const Function< T, NDIM > &f, bool fence=true)
Returns new function equal to alpha*f(x) with optional fence.
Definition: mra.h:1701
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
static double r2(const coord_3d &x)
Definition: smooth.h:45
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition: mra.h:2681
Vector< double, 6 > coord_6d
Definition: funcplot.h:1043
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
Function< double, 6 > real_function_6d
Definition: functypedefs.h:68
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
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::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition: corepotential.cc:55
Vector< T, N > unitvec(const Vector< T, N > &r, const double eps=1.e-6)
Construct a unit-Vector that has the same direction as r.
Definition: vector.h:729
Function< double, 3 > real_function_3d
Definition: functypedefs.h:65
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition: lapack.cc:832
SeparatedConvolution< double, 6 > real_convolution_6d
Definition: functypedefs.h:124
Vector< double, 3 > coord_3d
Definition: funcplot.h:1042
Derivative< double, 6 > real_derivative_6d
Definition: functypedefs.h:173
FunctionFactory< double, 6 > real_factory_6d
Definition: functypedefs.h:96
Derivative< double, 3 > real_derivative_3d
Definition: functypedefs.h:170
static const double c
Definition: relops.cc:10
static const double thresh
Definition: rk.cc:45
const double xi
Exponent for delta function approx.
Definition: siam_example.cc:60
Definition: test_ar.cc:204
void e()
Definition: test_sig.cc:75
#define N
Definition: testconv.cc:37
std::size_t axis
Definition: testpdiff.cc:59
double u(const double x, const double expnt)
Definition: testperiodic.cc:56
static Molecule molecule
Definition: testperiodicdft.cc:38
const double R2
Definition: vnucso.cc:84