8 #ifndef SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
9 #define SRC_APPS_CHEM_ELECTRONIC_CORRELATION_FACTOR_H_
20 class 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);
80 if(not op_mod.modified())
MADNESS_EXCEPTION(
"ElectronicCorrelationFactor::apply_U, op_mod must be in modified_NS form",1);
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??");
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);
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;
315 class 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();
404 r2.fill_tree(op_mod);
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);
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