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);
198 return (1.0-
e)*
u(rr,
dcut) + 0.5*gamma*
e;
203 struct f_over_r_ : FunctionFunctorInterface<double,6> {
206 f_over_r_(
double gamma,
double dcut) : gamma(gamma),
dcut(
dcut) {
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);
217 struct U : FunctionFunctorInterface<double,6> {
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]};
229 if (gamma>0.0)
return -0.5*exp(-gamma*rr)*
N[
axis];
240 struct f2_ : FunctionFunctorInterface<double,6> {
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);
252 struct nablaf2_ : FunctionFunctorInterface<double,6> {
254 nablaf2_(
double gamma) : gamma(gamma) {
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);
266 static double u(
double r,
double c) {
268 double r2 = r*r, pot;
271 }
else if (r > 1
e-2) {
272 pot = erf(r)/r + exp(-r2)*0.56418958354775630;
274 pot = 1.6925687506432689-r2*(0.94031597257959381-r2*(0.39493270848342941-0.12089776790309064*r2));
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);
286 static double x12(
const coord_6d& r,
const int axis) {
290 static coord_3d smoothed_unitvec(
const coord_3d& xyz,
double smoothing) {
295 const double r=xyz.
normf();
296 const double cutoff=smoothing;
300 const double xi=r/cutoff;
301 const double xi2=
xi*
xi;
302 const double xi3=
xi*
xi*
xi;
304 const double nu22=0.5 + 1./64.*(105*
xi - 175 *xi3 + 147* xi2*xi3 - 45* xi3*xi3*
xi);
306 const double kk=2.*nu22-1.0;
307 return kk/(r+1.e-15)*xyz;
315class CorrelationFactor2 {
319 typedef std::shared_ptr< FunctionFunctorInterface<double,6> >
functorT;
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) {
335 double gamma()
const {
return _gamma;}
366 const double bsh_thresh=1.e-7;
371 op_mod.modified()=
true;
383 .g12(u1).ket(
copy(Drhs1));
387 .g12(u1).ket(
copy(Drhs2));
391 result=result+(tmp1-tmp2).
truncate();
395 result.truncate().reduce_rank();
413 class R_functor :
public FunctionFunctorInterface<double,6> {
419 R_functor(
double gamma,
int e=1) : gamma(gamma), exponent(
e) {
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;
437 class U2_functor :
public FunctionFunctorInterface<double,6> {
441 U2_functor(
double gamma) : gamma(gamma) {
446 double operator()(
const coord_6d& r)
const {
447 const double rr=r12(r);
450 return (5./4.0 - rr + (35.0* rr*rr)/48.0 - (101.0*rr*rr*rr)/192.0);
452 const double egr=exp(-gamma*rr);
453 return -(-8.*egr + 8.0 + rr*egr)/(4.0 *rr*egr - 8 *rr);
463 class U1_functor :
public FunctionFunctorInterface<double,6> {
468 U1_functor(
double gamma,
int axis) : gamma(gamma),
axis(
axis) {
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]};
480 val = 0.5 - 0.5*rr + 0.125*(3.*rr*rr) - (13.* rr*rr*rr)/48.0;
482 const double egr=exp(-gamma*rr);
483 val=egr/(4.0-2.0*egr);
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);
FunctionFactory & functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > &f)
Definition function_factory.h:141
virtual FunctionFactory & is_on_demand()
Definition function_factory.h:281
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:1125
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(* 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:1701
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
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
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition vmra.h:851
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition lapack.cc:832
Function< double, 6 > real_function_6d
Definition functypedefs.h:68
Function< double, 3 > real_function_3d
Definition functypedefs.h:65
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
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
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:38
const double R2
Definition vnucso.cc:84