8#ifndef SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
9#define SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
20class CorrelationFactor {
30 CorrelationFactor(World& world) : world(world), _gamma(-1.0),
dcut(1.e-10),
35 CorrelationFactor(World& world,
const double& gamma,
const double dcut,
44 CorrelationFactor(World& world,
const double& gamma,
const double dcut,
45 const double lo) : world(world), _gamma(gamma),
dcut(
dcut),
lo(
lo) {
53 CorrelationFactor(
const CorrelationFactor& other) : world(other.world) {
60 CorrelationFactor& operator=(
const CorrelationFactor& other) {
68 double gamma()
const {
return _gamma;}
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);
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;
98 tmp1.fill_cuspy_tree(op_mod).truncate();
103 tmp2=CompositeFactory<double,6,3>(world)
105 tmp2.fill_cuspy_tree(op_mod).truncate();
108 result=result+(tmp1-tmp2).
truncate();
114 result.truncate().reduce_rank();
122 .g12(fg3).particle1(
copy(phi_i)).particle2(
copy(phi_j)).thresh(
thresh);;
123 mul.fill_cuspy_tree(op_mod).truncate();
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??");
151 double thresh=FunctionDefaults<3>::get_thresh();
180 nablaf2_
func(_gamma);
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));
288 static double r12(
const coord_6d& r) {
289 const double x12=r[0]-r[3];
290 const double y12=r[1]-r[4];
291 const double z12=r[2]-r[5];
292 const double r12=sqrt(x12*x12 + y12*y12 + z12*z12);
295 static double x12(
const coord_6d& r,
const int axis) {
299 static coord_3d smoothed_unitvec(
const coord_3d& xyz,
double smoothing) {
304 const double r=xyz.
normf();
305 const double cutoff=smoothing;
309 const double xi=r/cutoff;
310 const double xi2=
xi*
xi;
311 const double xi3=
xi*
xi*
xi;
313 const double nu22=0.5 + 1./64.*(105*
xi - 175 *xi3 + 147* xi2*xi3 - 45* xi3*xi3*
xi);
315 const double kk=2.*nu22-1.0;
316 return kk/(r+1.e-15)*xyz;
324class CorrelationFactor2 {
328 typedef std::shared_ptr< FunctionFunctorInterface<double,6> >
functorT;
338 CorrelationFactor2(World& world) : world(world), _gamma(0.5),
dcut(1.e-10),
339 lo(1.e-10), vtol(FunctionDefaults<3>::get_thresh()*0.1) {
344 double gamma()
const {
return _gamma;}
375 const double bsh_thresh=1.e-7;
380 op_mod.modified()=
true;
392 .g12(u1).ket(
copy(Drhs1));
396 .g12(u1).ket(
copy(Drhs2));
400 result=result+(tmp1-tmp2).
truncate();
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);
500 static double r12(
const coord_6d& r) {
501 const double x12=r[0]-r[3];
502 const double y12=r[1]-r[4];
503 const double z12=r[2]-r[5];
504 const double r12=sqrt(x12*x12 + y12*y12 + z12*z12);
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
std::vector< Spinor > truncate(std::vector< Spinor > arg)
Definition dirac-hatom.cc:499
Fcwf copy(Fcwf psi)
Definition fcwf.cc:338
static double function(const coord_3d &r)
Normalized gaussian.
Definition functionio.cc:100
double square(double x)
Definition gfit.cc:59
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
double func(double x)
A simple program for testing the CubicInterpolationTable class.
Definition interp3.cc:43
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
void print(const tensorT &t)
Definition mcpfit.cc:140
Main include file for MADNESS and defines Function interface.
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
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
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
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition vmra.h:871
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
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
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