8#ifndef SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
9#define SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
68 double gamma()
const {
return _gamma;}
71 double operator()(
const coord_6d& r)
const {
72 const double rr=
r12(r);
82 const bool debug =
false;
114 result.truncate().reduce_rank();
128 if(
debug) result.print_size(
"Ue|ij>");
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??");
189 struct fg_ : FunctionFunctorInterface<double,6> {
192 fg_(
double gamma,
double dcut) : gamma(gamma),
dcut(
dcut) {
195 double operator()(
const coord_6d& r)
const {
196 const double rr=r12(r);
197 const double e=exp(-gamma*rr);
199 double value= 1.5*gamma - 1.*std::pow(gamma,2)*rr + 0.41666666666666663*std::pow(gamma,3)*std::pow(rr,2) -
200 0.125*std::pow(gamma,4)*std::pow(rr,3) + 0.029166666666666667*std::pow(gamma,5)*std::pow(rr,4)
201 + 0.005555555555555556*std::pow(gamma,6)*std::pow(rr,5) +
202 0.0008928571428571429*std::pow(gamma,7)*std::pow(rr,6);
206 return (1.0-
e)/rr + 0.5*gamma*
e;
212 struct f_over_r_ : FunctionFunctorInterface<double,6> {
215 f_over_r_(
double gamma,
double dcut) : gamma(gamma),
dcut(
dcut) {
218 double operator()(
const coord_6d& r)
const {
219 const double rr=r12(r);
220 const double e=exp(-gamma*rr);
221 return (1.0-
e)*
u(rr,
dcut)/(2.0*gamma);
226 struct U : FunctionFunctorInterface<double,6> {
234 double operator()(
const coord_6d& r)
const {
235 const double rr=r12(r);
236 const coord_3d vr12{r[0]-r[3],r[1]-r[4],r[2]-r[5]};
238 if (gamma>0.0)
return -0.5*exp(-gamma*rr)*
N[
axis];
249 struct f2_ : FunctionFunctorInterface<double,6> {
252 double operator()(
const coord_6d& r)
const {
253 const double rr=r12(r);
254 const double e=exp(-gamma*rr);
255 const double f=(1.0-
e)/(2.0*gamma);
261 struct nablaf2_ : FunctionFunctorInterface<double,6> {
263 nablaf2_(
double gamma) : gamma(gamma) {
267 double operator()(
const coord_6d& r)
const {
268 const double rr=r12(r);
269 const double f=exp(-2.0*gamma*rr)/(4.0*gamma*gamma);
275 static double u(
double r,
double c) {
277 double r2 = r*r, pot;
280 }
else if (r > 1
e-2) {
281 pot =
erf(r)/r +
exp(-
r2)*0.56418958354775630;
283 pot = 1.6925687506432689-
r2*(0.94031597257959381-
r2*(0.39493270848342941-0.12089776790309064*
r2));
289 const double x12=r[0]-r[3];
290 const double y12=r[1]-r[4];
291 const double z12=r[2]-r[5];
309 const double xi=r/cutoff;
315 const double kk=2.*
nu22-1.0;
316 return kk/(r+1.e-15)*
xyz;
328 typedef std::shared_ptr< FunctionFunctorInterface<double,6> >
functorT;
339 lo(1.e-10), vtol(FunctionDefaults<3>::get_thresh()*0.1) {
344 double gamma()
const {
return _gamma;}
404 result.truncate().reduce_rank();
422 class R_functor :
public FunctionFunctorInterface<double,6> {
428 R_functor(
double gamma,
int e=1) : gamma(gamma), exponent(
e) {
433 double operator()(
const coord_6d& r)
const {
434 const double rr=r12(r);
435 double val=(1.0-0.5*exp(-gamma*rr));
436 if (exponent==1)
return val;
437 else if (exponent==2)
return val*val;
438 else if (exponent==-1)
return 1.0/val;
446 class U2_functor :
public FunctionFunctorInterface<double,6> {
450 U2_functor(
double gamma) : gamma(gamma) {
455 double operator()(
const coord_6d& r)
const {
456 const double rr=r12(r);
459 return (5./4.0 - rr + (35.0* rr*rr)/48.0 - (101.0*rr*rr*rr)/192.0);
461 const double egr=exp(-gamma*rr);
462 return -(-8.*egr + 8.0 + rr*egr)/(4.0 *rr*egr - 8 *rr);
472 class U1_functor :
public FunctionFunctorInterface<double,6> {
477 U1_functor(
double gamma,
int axis) : gamma(gamma),
axis(
axis) {
482 double operator()(
const coord_6d& r)
const {
483 const double rr=r12(r);
484 const coord_3d vr12{r[0]-r[3],r[1]-r[4],r[2]-r[5]};
489 val = 0.5 - 0.5*rr + 0.125*(3.*rr*rr) - (13.* rr*rr*rr)/48.0;
491 const double egr=exp(-gamma*rr);
492 val=egr/(4.0-2.0*egr);
501 const double x12=r[0]-r[3];
502 const double y12=r[1]-r[4];
503 const double z12=r[2]-r[5];
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:176
FunctionFactory & functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > &f)
Definition function_factory.h:141
virtual FunctionFactory & is_on_demand()
Definition function_factory.h:288
Function< T, NDIM > & fill_tree(const Function< R, NDIM > &g, bool fence=true)
With this being an on-demand function, fill the MRA tree according to different criteria.
Definition mra.h:1192
Function< T, NDIM > & fill_cuspy_tree(const opT &op, const bool fence=true)
Definition mra.h:1233
double thresh() const
Returns value of truncation threshold. No communication.
Definition mra.h:607
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:642
constexpr T normf() const
Calculate the 2-norm of the vector elements.
Definition vector.h:421
static const double R
Definition csqrt.cc:46
double(* f)(const coord_3d &)
Definition derivatives.cc:54
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
static double function(const coord_3d &r)
Normalized gaussian.
Definition functionio.cc:100
double psi(const Vector< double, 3 > &r)
Definition hatom_energy.cc:78
static double u(double r, double c)
Definition he.cc:20
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.
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
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:2746
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:1765
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
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:2367
static double r2(const coord_3d &x)
Definition smooth.h:45
Vector< double, 6 > coord_6d
Definition funcplot.h:1044
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition vmra.h:871
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition lapack.cc:832
constexpr 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:768
Function< double, 6 > real_function_6d
Definition functypedefs.h:83
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:2481
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition corepotential.cc:55
Function< double, 3 > real_function_3d
Definition functypedefs.h:80
SeparatedConvolution< double, 6 > real_convolution_6d
Definition functypedefs.h:139
Vector< double, 3 > coord_3d
Definition funcplot.h:1043
Derivative< double, 6 > real_derivative_6d
Definition functypedefs.h:188
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:371
FunctionFactory< double, 6 > real_factory_6d
Definition functypedefs.h:111
Derivative< double, 3 > real_derivative_3d
Definition functypedefs.h:185
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:2066
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
static Molecule molecule
Definition testperiodicdft.cc:39
const double R2
Definition vnucso.cc:84