MADNESS 0.10.1
correlationfactor.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30*/
31
32/*!
33 \file apps/chem/correlationfactor.h
34 \brief class for regularizing singular potentials in the molecular
35 Hamilton operator
36
37 \par Introduction
38
39 The correlation factors are intended to represent the cusps (nuclear or
40 electronic) in the molecular wave function. Their commutator over the
41 kinetic energy operator give a potential operator that cancels the
42 singular potential
43
44 [T, f12] = U + 1/r12
45 R^-1[T, R] = U_nuc - sum_A Z_A/r1A
46
47 The regularized potentials U and U_nuc contain two terms each, a local
48 potential, (denoted U2), and a potential that is multiplied with the
49 derivative operator \nabla (denoted U1)
50
51 U = U1 . (\nabla_1 - \nabla_2) + U2
52 U_nuc = U1 . \nabla + U2
53
54 with
55
56 U2 = (\nabla^2 f12)
57 U2 = R^{-1}(\nabla^2 R)
58
59 \vec U1 = (\vec \nabla f12)
60 \vec U1 = R^{-1}(\vec \nabla R)
61
62 To construct a nuclear correlation factor write:
63
64 std::shared_ptr<NuclearCorrelationFactor> nuclear_correlation
65 = create_nuclear_correlation_factor(world,*calc);
66
67 where calc is an SCF calculation which holds the molecule and the
68 nuclear_corrfac parameter name.
69*/
70
71
72#ifndef MADNESS_CHEM_NUCLEARCORRELATIONFACTOR_H_
73#define MADNESS_CHEM_NUCLEARCORRELATIONFACTOR_H_
74
75#include <madness/mra/mra.h>
79
80namespace madness {
81
82/// ABC for the nuclear correlation factors
84public:
87 typedef std::shared_ptr< FunctionFunctorInterface<double,3> > functorT;
88
89 /// ctor
90
91 /// @param[in] world the world
92 /// @param[in] mol molecule with the sites of the nuclei
94 : world(world), vtol(FunctionDefaults<3>::get_thresh()*0.1)
95 , eprec(mol.get_eprec()), molecule(mol) {}
96
97 /// virtual destructor
99
100 /// initialize the regularized potentials U1 and U2
101 void initialize(const double vtol1) {
102
103 // set threshold for projections
104 vtol=vtol1;
105
106 // construct the potential functions
107 // keep tighter threshold for orthogonalization
108 for (int axis=0; axis<3; ++axis) {
109 functorT U1f=functorT(new U1_functor(this,axis));
111 .functor(U1f).truncate_on_project());
112 U1_function.back().set_thresh(FunctionDefaults<3>::get_thresh());
113 }
114
115 // U2 is the term -S"/S - Z/r
116 functorT U2f=functorT(new U2_functor(this));
118 .functor(U2f).truncate_on_project();
120
121 // U3 is the term SA'/SA . SB'/SB
122 functorT U3f=functorT(new U3_functor(this));
124 .functor(U3f).truncate_on_project();
126 U2_function+=tmp;
128 }
129
130 virtual corrfactype type() const = 0;
131
132 /// apply the regularized potential U_nuc on a given function rhs
133 virtual real_function_3d apply_U(const real_function_3d& rhs) const {
134
135 // the purely local part
136 real_function_3d result=(U2()*rhs).truncate();
137
138 // the part with the derivative operators
139 result.compress();
140 for (int axis=0; axis<3; ++axis) {
141 real_derivative_3d D = free_space_derivative<double,3>(world, axis);
142 const real_function_3d Drhs=D(rhs).truncate();
143 result+=U1(axis)*Drhs;
144 }
145
146 result.truncate();
147 return result;
148 }
149
150 /// return the nuclear correlation factor
151 virtual real_function_3d function() const {
152 functorT Rf=functorT(new R_functor(this,1));
154 .functor(Rf).truncate_on_project();
155 return r;
156 }
157
158 /// return the square of the nuclear correlation factor
159 virtual real_function_3d square() const {
160 R_functor r(this,2);
162 .functor(r).truncate_on_project();
163 return R2;
164 }
165
166 /// return the square of the nuclear correlation factor multiplied with
167 /// the derivative of the nuclear potential for the specified atom
168
169 /// @return R^2 * \frac{\partial Z_A/r_{1A}}{\partial X_A}
170 virtual real_function_3d square_times_V_derivative(const int iatom, const int axis) const {
173 .functor(func).truncate_on_project();
174 return R2;
175 }
176
177 /// return the inverse nuclear correlation factor
178 virtual real_function_3d inverse() const {
179 R_functor r(this,-1);
181 .functor(r).truncate_on_project();
182 return R_inverse;
183 }
184
185 /// return the U1 term of the correlation function
186 virtual const real_function_3d U1(const int axis) const {
187 return U1_function[axis];
188 }
189
190 /// return the U1 functions in a vector
191 std::vector<real_function_3d> U1vec() const {
192 std::vector<real_function_3d> uvec(3);
193 uvec[0]=U1_function[0];
194 uvec[1]=U1_function[1];
195 uvec[2]=U1_function[2];
196 return uvec;
197 }
198
199 /// return the U2 term of the correlation function
200 virtual const real_function_3d U2() const {return U2_function;}
201
202private:
203
204 /// the world
206
207 /// the threshold for initial projection
208 double vtol;
209
210 /// smoothing of the potential/step function
211 double eprec;
212
213 /// the molecule
215
216protected:
217 /// the three components of the U1 potential
218 std::vector<real_function_3d> U1_function;
219
220 /// the purely local U2 potential, having absorbed the nuclear pot V_nuc
222
223private:
224 /// the correlation factor S wrt a given atom
225
226 /// @param[in] r the distance of the req'd coord to the nucleus
227 /// @param[in] Z the nuclear charge
228 /// @return the nuclear correlation factor S_A(r_1A)
229 virtual double S(const double& r, const double& Z) const = 0;
230
231 /// the partial derivative of correlation factor S' wrt the cartesian coordinates
232
233 /// @param[in] vr1A the vector of the req'd coord to the nucleus
234 /// @param[in] Z the nuclear charge
235 /// @return the gradient of the nuclear correlation factor S'_A(r_1A)
236 virtual coord_3d Sp(const coord_3d& vr1A, const double& Z) const = 0;
237
238 /// the regularized potential wrt a given atom wrt the cartesian coordinate
239
240 /// S" is the Cartesian Laplacian applied on the NCF. Note the difference
241 /// to Srr_div_S, which is the second derivative wrt the distance rho.
242 /// this is: -S"/S - Z/r
243 /// @param[in] r the distance of the req'd coord to the nucleus
244 /// @param[in] Z the nuclear charge
245 /// @return the Laplacian of the nuclear correlation factor divided
246 /// by the correlation factor minus the nuclear potential
247 virtual double Spp_div_S(const double& r, const double& Z) const = 0;
248
249public:
250
251 /// first derivative of the NCF with respect to the relative distance rho
252 /// \f[
253 /// \frac{\partial S(\rho)}{\partial \rho} \frac{1}{S(\rho)}
254 /// \f]
255 /// where the distance of the electron to the nucleus A is given by
256 /// \f[
257 /// \rho = |\vec r - \vec R_A |
258 /// \f]
259 virtual double Sr_div_S(const double& r, const double& Z) const = 0;
260
261 /// second derivative of the NCF with respect to the relative distance rho
262 /// \f[
263 /// \frac{\partial^2 S(\rho)}{\partial \rho^2} \frac{1}{S(\rho)}
264 /// \f]
265 /// where the distance of the electron to the nucleus A is given by
266 /// \f[
267 /// \rho = |\vec r - \vec R_A |
268 /// \f]
269 virtual double Srr_div_S(const double& r, const double& Z) const = 0;
270
271 /// third derivative of the NCF with respect to the relative distance rho
272 /// \f[
273 /// \frac{\partial^3 S(\rho)}{\partial \rho^3} \frac{1}{S(\rho)}
274 /// \f]
275 /// where the distance of the electron to the nucleus A is given by
276 /// \f[
277 /// \rho = |\vec r - \vec R_A |
278 /// \f]
279 virtual double Srrr_div_S(const double& r, const double& Z) const = 0;
280
281 /// derivative of the U2 potential wrt nuclear coordinate X (spherical part)
282
283 /// need to reimplement this for all derived classes due to the
284 /// range for r -> 0, where the singular terms cancel. With
285 /// \f[
286 /// \rho = \left| \vec r- \vec R_A \right|
287 /// \f]
288 /// returns the term in the parenthesis without the the derivative of rho
289 /// \f[
290 /// \frac{\partial U_2}{\partial X_A} = \frac{\partial \rho}{\partial X}
291 /// \left(-\frac{1}{2}\frac{S''' S - S'' S'}{S^2} + \frac{1}{\rho^2}\frac{S'}{S}
292 /// - \frac{1}{\rho} \frac{S''S - S'^2}{S^2} + \frac{Z_A}{\rho^2}\right)
293 /// \f]
294 virtual double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
295 if (world.rank()==0) {
296 print("you can't compute the Hessian matrix");
297 print("U2X_spherical is not implemented for the nuclear correlation factor");
298 }
299 MADNESS_EXCEPTION("do more implementation work",1);
300 }
301
302public:
303
304 /// smoothed unit vector for the computation of the U1 potential
305
306 /// note the identity for exchanging nuclear and electronic coordinates
307 /// (there is a sign change, unlike for the smoothed potential)
308 /// \f[
309 /// \vec n = \frac{\partial \rho}{\partial x} = -\frac{\partial \rho}{\partial X}
310 /// \f]
311 /// \f[
312 /// \vec n = \left\{\frac{x \mathrm{erf}\left(\frac{r}{s}\right)}{r},
313 /// \frac{y \mathrm{erf}\left(\frac{r}{s}\right)}{r},
314 /// \frac{z \mathrm{erf}\left(\frac{r}{s}\right)}{r}\right\}
315 /// \f]
316 coord_3d smoothed_unitvec(const coord_3d& xyz, double smoothing=0.0) const {
317#if 0
318
319 if (smoothing==0.0) smoothing=molecule.get_eprec();
320 // TODO:need to test this
321 // reduce the smoothing for the unitvector
322 //if (not (this->type()==None or this->type()==Two)) smoothing=sqrt(smoothing);
323 smoothing=sqrt(smoothing);
324 const double r=xyz.normf();
325 const double rs=r/smoothing;
326 if (r<1.e-4) {
327 const double sqrtpi=sqrt(constants::pi);
328 double erfrs_div_r=2.0/(smoothing*sqrtpi)-2.0/3.0*rs*rs/(sqrtpi*smoothing);
329 return erfrs_div_r*xyz;
330 } else if (r<6.0) {
331 return erf(rs)/r*xyz;
332 } else {
333 return 1.0/r*xyz;
334 }
335
336
337
338#else
339 if (smoothing==0.0) smoothing=eprec;
340 // TODO:need to test this
341 // reduce the smoothing for the unitvector
342 //if (not (this->type()==None or this->type()==Two)) smoothing=sqrt(smoothing);
343 const double r=xyz.normf();
344 const double cutoff=smoothing;
345 if (r>cutoff) {
346 return 1.0/r*xyz;
347 } else {
348 const double xi=r/cutoff;
349 const double xi2=xi*xi;
350 const double xi3=xi*xi*xi;
351// const double nu21=0.5+1./32.*(45.*xi - 50.*xi3 + 21.*xi*xi*xi*xi*xi);
352 const double nu22=0.5 + 1./64.*(105* xi - 175 *xi3 + 147* xi2*xi3 - 45* xi3*xi3*xi);
353// const double nu40=0.5 + 1./128.*(225 *xi - 350 *xi3 + 189*xi2*xi3);
354 const double kk=2.*nu22-1.0;
355 return kk/r*xyz;
356 }
357
358#endif
359 }
360
361 /// derivative of smoothed unit vector wrt the *electronic* coordinate
362
363 /// note the sign change for exchanging nuclear and electronic coordinates
364 /// \f[
365 /// \frac{\partial \vec n}{\partial x} = -\frac{\partial \vec n}{\partial X}
366 /// \f]
367 /// the derivative wrt x is given by
368 /// \f[
369 /// \frac{\partial\vec n}{\partial x} =
370 /// \left\{\frac{\left(r^2-x^2\right) \mathrm{erf}\left(\frac{r}{s}\right)}{r^3}
371 /// +\frac{2 x^2 e^{-\frac{r^2}{s^2}}}{\sqrt{\pi } r^2 s},
372 /// x y \left(\frac{2 e^{-\frac{r^2}{s^2}}}{\sqrt{\pi } r^2 s}
373 /// -\frac{\mathrm{erf}\left(\frac{r}{s}\right)}{r^3}\right),
374 /// x z \left(\frac{2 e^{-\frac{r^2}{s^2}}}{\sqrt{\pi } r^2 s}
375 /// -\frac{\mathrm{erf}\left(\frac{r}{s}\right)}{r^3}\right)\right\}
376 /// \f]
378 double smoothing=0.0) const {
379
380 const double r=xyz.normf();
381 coord_3d result;
382 if (smoothing==0.0) smoothing=eprec;
383
384#if 1
385 // TODO:need to test this
386 // reduce the smoothing for the unitvector
387 //if (not (this->type()==None or this->type()==Two)) smoothing=sqrt(smoothing);
388 smoothing=sqrt(smoothing);
389
390 const double rs=r/smoothing;
391 const static double sqrtpi=sqrt(constants::pi);
392 const double sqrtpis3=sqrtpi*smoothing*smoothing*smoothing;
393
394 if (r<1.e-4) {
395 // series expansion
396 double p=-4.0/(3.0*sqrtpis3) + 4.0*rs*rs/(5.0*sqrtpis3);
397
398 double erfrs_div_r=2.0/(smoothing*sqrtpi)-2.0/3.0*rs*rs/(sqrtpi*smoothing);
399 result=xyz*xyz[axis]*p;
400 result[axis]+=erfrs_div_r;
401
402 } else if (r<6.0) {
403 const double erfrs_div_r=erf(rs)/r;
404 const double term1=2.0*exp(-rs*rs)/(sqrtpi*r*r*smoothing);
405 result=xyz*xyz[axis]*(term1-erfrs_div_r/(r*r));
406 result[axis]+=erfrs_div_r;
407#else
408 if (r<smoothing) {
409 double r2=r*r;
410 double s2=smoothing*smoothing;
411 double s7=s2*s2*s2*smoothing;
412 double x2=xyz[axis]*xyz[axis];
413
414 double fac_offdiag=-(((135. *r2*r2 - 294.* r2 *s2
415 + 175.*s2*s2))/(16.* s7));
416 double fac_diag=-((45.* r2*r2*r2 - 147.* r2*r2* s2
417 + 175.* r2*s2*s2 - 105.* s2*s2*s2 + 270.* r2*r2* x2
418 - 588.* r2* s2* x2 + 350.*s2* s2 *x2)/(32.* s7));
419
420 result[0]=fac_offdiag*xyz[0]*xyz[axis];
421 result[1]=fac_offdiag*xyz[1]*xyz[axis];
422 result[2]=fac_offdiag*xyz[2]*xyz[axis];
423 result[axis]=fac_diag;
424
425#endif
426 } else {
427 result=xyz*(-xyz[axis]/(r*r*r));
428 result[axis]+=1/r;
429 }
430 return result;
431 }
432
433
434 class R_functor : public FunctionFunctorInterface<double,3> {
437 public:
439 : ncf(ncf), exponent(e) {}
440 double operator()(const coord_3d& xyz) const {
441 double result=1.0;
442 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
443 const Atom& atom=ncf->molecule.get_atom(i);
444 const coord_3d vr1A=xyz-atom.get_coords();
445 const double r=vr1A.normf();
446 result*=ncf->S(r,atom.q);
447 }
448 if (exponent==-1) return 1.0/result;
449 else if (exponent==2) return result*result;
450 else if (exponent==1) return result;
451 else {
452 return std::pow(result,double(exponent));
453 }
454
455 }
456 std::vector<coord_3d> special_points() const {
457 return ncf->molecule.get_all_coords_vec();
458 }
459 };
460
461 /// functor for the local part of the U1 potential -- NOTE THE SIGN
462
463 /// U1 = -S'/S
464 class U1_functor : public FunctionFunctorInterface<double,3> {
465
467 const int axis;
468
469 public:
471 : ncf(ncf), axis(axis) {}
472
473 double operator()(const coord_3d& xyz) const {
474 double result=0.0;
475 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
476 const Atom& atom=ncf->molecule.get_atom(i);
477 const coord_3d vr1A=xyz-atom.get_coords();
478 const double r=vr1A.normf();
479 const double& Z=atom.q;
480// result-=(ncf->Sp(vr1A,Z)[axis]/ncf->S(r,Z));
481 result-=ncf->Sr_div_S(r,Z)*ncf->smoothed_unitvec(vr1A)[axis];
482 }
483 return result;
484 }
485 std::vector<coord_3d> special_points() const {
486 return ncf->molecule.get_all_coords_vec();
487 }
488 };
489
490 /// U1 functor for a specific atom
491
492 /// NOTE THE SIGN !!
493 /// this is
494 /// \f[
495 /// -\frac{\partial \rho}{\partial X_A}\frac{\partial S}{\partial \rho}\frac{1}{S}
496 /// \f]
498
500 const size_t iatom;
501 const int axis;
502
503 public:
505 const int axis) : ncf(ncf), iatom(atom), axis(axis) {}
506
507 double operator()(const coord_3d& xyz) const {
508 const Atom& atom=ncf->molecule.get_atom(iatom);
509 const coord_3d vr1A=xyz-atom.get_coords();
510 const double r=vr1A.normf();
511 const double& Z=atom.q;
512 return ncf->Sr_div_S(r,Z)*ncf->smoothed_unitvec(vr1A)[axis];
513 }
514
515 std::vector<coord_3d> special_points() const {
516 std::vector< madness::Vector<double,3> > c(1);
517 const Atom& atom=ncf->molecule.get_atom(iatom);
518 c[0][0]=atom.x;
519 c[0][1]=atom.y;
520 c[0][2]=atom.z;
521 return c;
522 }
523 };
524
525
526 /// functor for a local U1 dot U1 potential
527
528 /// \f[
529 /// U1\dot U1 = \frac{\left(S^r_A S^r_B\right)}{S_A S_B} n_A \cdot n_B
530 /// \f]
531 /// with positive sign!
533
535
536 public:
538
539 double operator()(const coord_3d& xyz) const {
540 std::vector<double> Sr_div_S(ncf->molecule.natom());
541 std::vector<coord_3d> unitvec(ncf->molecule.natom());
542 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
543 const Atom& atom=ncf->molecule.get_atom(i);
544 const coord_3d vr1A=xyz-atom.get_coords();
545 const double r=vr1A.normf();
546 Sr_div_S[i]=ncf->Sr_div_S(r,atom.q);
547 unitvec[i]=ncf->smoothed_unitvec(vr1A);
548 }
549
550 double result=0.0;
551 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
552 for (size_t j=0; j<ncf->molecule.natom(); ++j) {
553 double tmp=Sr_div_S[i]*Sr_div_S[j];
554 if (i!=j) tmp*=inner(unitvec[i],unitvec[j]);
555 result+=tmp;
556 }
557 }
558
559
560 return result;
561 }
562 std::vector<coord_3d> special_points() const {
563 return ncf->molecule.get_all_coords_vec();
564 }
565 };
566
567
568 class U2_functor : public FunctionFunctorInterface<double,3> {
570 public:
572 double operator()(const coord_3d& xyz) const {
573 double result=0.0;
574 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
575 const Atom& atom=ncf->molecule.get_atom(i);
576 const coord_3d vr1A=xyz-atom.get_coords();
577 const double r=vr1A.normf();
578 result+=ncf->Spp_div_S(r,atom.q);
579 }
580 return result;
581 }
582 std::vector<coord_3d> special_points() const {
583 return ncf->molecule.get_all_coords_vec();
584 }
585 };
586
587 class U3_functor : public FunctionFunctorInterface<double,3> {
589 public:
591 double operator()(const coord_3d& xyz) const {
592 std::vector<coord_3d> all_terms(ncf->molecule.natom());
593 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
594 const Atom& atom=ncf->molecule.get_atom(i);
595 const coord_3d vr1A=xyz-atom.get_coords();
596 const double r=vr1A.normf();
597// all_terms[i]=ncf->Sp(vr1A,atom.q)*(1.0/ncf->S(r,atom.q));
598 all_terms[i]=ncf->Sr_div_S(r,atom.q)*ncf->smoothed_unitvec(vr1A);
599 }
600
601 double result=0.0;
602 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
603 for (size_t j=0; j<i; ++j) {
604 result+=all_terms[i][0]*all_terms[j][0]
605 +all_terms[i][1]*all_terms[j][1]
606 +all_terms[i][2]*all_terms[j][2];
607 }
608 }
609
610 return -1.0*result;
611 }
612 std::vector<coord_3d> special_points() const {
613 return ncf->molecule.get_all_coords_vec();
614 }
615 };
616
617 /// U2 functor for a specific atom
619
621 const size_t iatom;
622
623 public:
625 : ncf(ncf), iatom(atom) {}
626
627 double operator()(const coord_3d& xyz) const {
628 const Atom& atom=ncf->molecule.get_atom(iatom);
629 const coord_3d vr1A=xyz-atom.get_coords();
630 const double r=vr1A.normf();
631 return ncf->Spp_div_S(r,atom.q);
632 }
633
634 std::vector<coord_3d> special_points() const {
635 std::vector< madness::Vector<double,3> > c(1);
636 const Atom& atom=ncf->molecule.get_atom(iatom);
637 c[0][0]=atom.x;
638 c[0][1]=atom.y;
639 c[0][2]=atom.z;
640 return c;
641 }
642 };
643
644 /// U3 functor for a specific atom
646
648 const size_t iatom;
649
650 public:
652 : ncf(ncf), iatom(atom) {}
653
654 double operator()(const coord_3d& xyz) const {
655 const Atom& atomA=ncf->molecule.get_atom(iatom);
656 const coord_3d vr1A=xyz-atomA.get_coords();
657 const double rA=vr1A.normf();
658 const coord_3d nA=ncf->smoothed_unitvec(vr1A);
659 double Sr_div_SA=ncf->Sr_div_S(rA,atomA.q);
660
661 double result=0.0;
662 // sum over B
663 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
664 if (i==iatom) continue; // restricted sum
665
666 const Atom& atomB=ncf->molecule.get_atom(i);
667 const coord_3d vr1B=xyz-atomB.get_coords();
668 const double rB=vr1B.normf();
669 const coord_3d nB=ncf->smoothed_unitvec(vr1B);
670 double Sr_div_SB=ncf->Sr_div_S(rB,atomB.q);
671
672 double dot=nA[0]*nB[0] + nA[1]*nB[1] + nA[2]*nB[2];
673 result+=Sr_div_SB*Sr_div_SA*dot;
674 }
675 return -0.5*result;
676 }
677
678 std::vector<coord_3d> special_points() const {
679 std::vector< madness::Vector<double,3> > c(1);
680 const Atom& atom=ncf->molecule.get_atom(iatom);
681 c[0][0]=atom.x;
682 c[0][1]=atom.y;
683 c[0][2]=atom.z;
684 return c;
685 }
686 };
687
691 const size_t iatom;
692 public:
694 const Molecule& mol, const size_t iatom1)
695 : ncf(ncf), molecule(mol), iatom(iatom1) {}
696 double operator()(const coord_3d& xyz) const {
697 double result=1.0;
698 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
699 const Atom& atom=ncf->molecule.get_atom(i);
700 const coord_3d vr1A=xyz-atom.get_coords();
701 const double r=vr1A.normf();
702 result*=ncf->S(r,atom.q);
703 }
705 iatom, xyz[0], xyz[1], xyz[2]);
706 return result*result*V;
707
708 }
709 std::vector<coord_3d> special_points() const {
710 return ncf->molecule.get_all_coords_vec();
711 }
712 };
713
714
718 const size_t iatom;
719 const int axis;
720 public:
722 const Molecule& molecule1, const size_t atom1, const int axis1)
723 : ncf(ncf), molecule(molecule1), iatom(atom1), axis(axis1) {}
724 double operator()(const coord_3d& xyz) const {
725 double result=1.0;
726 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
727 const Atom& atom=ncf->molecule.get_atom(i);
728 const coord_3d vr1A=xyz-atom.get_coords();
729 const double r=vr1A.normf();
730 result*=ncf->S(r,atom.q);
731 }
733 iatom, axis, xyz[0], xyz[1], xyz[2]);
734 return result*result*Vprime;
735
736 }
737 std::vector<coord_3d> special_points() const {
738 return ncf->molecule.get_all_coords_vec();
739 }
740 };
741
742 /// compute the derivative of R wrt the displacement of atom A, coord axis
743 class RX_functor : public FunctionFunctorInterface<double,3> {
746 const int derivativeaxis; /// direction of the derivative operator
747 const int exponent; /// 1 or 2 -> R^X or R^X R
748
749 public:
751 const int daxis, const int exponent) : ncf(ncf), thisatom(atom1),
753 MADNESS_ASSERT((exponent==1) or (exponent==2) or (exponent==-1));
754 }
755
756 RX_functor(const NuclearCorrelationFactor* ncf, const int iatom,
757 const int daxis, const int exponent) : ncf(ncf),
758 thisatom(ncf->molecule.get_atom(iatom)),
760 MADNESS_ASSERT((exponent==1) or (exponent==2) or (exponent==-1));
761 }
762
763 double operator()(const coord_3d& xyz) const {
764
765 // compute the R term
766 double result=1.0;
767 if ((exponent==1) or (exponent==2)) {
768 for (size_t i=0; i<ncf->molecule.natom(); ++i) {
769 const Atom& atom=ncf->molecule.get_atom(i);
770 const coord_3d vr1A=xyz-atom.get_coords();
771 const double r=vr1A.normf();
772 result*=ncf->S(r,atom.q);
773 }
774 if (exponent==2) result=result*result;
775 }
776
777 // compute the derivative term
778 {
779 const coord_3d vr1A=xyz-thisatom.get_coords();
780 const double r=vr1A.normf();
781 const double& Z=thisatom.q;
782 const double S1=-ncf->Sr_div_S(r,Z) // note the sign
783 *ncf->smoothed_unitvec(vr1A)[derivativeaxis];
784 result*=S1;
785 }
786 return result;
787 }
788
789 std::vector<coord_3d> special_points() const {
790 return ncf->molecule.get_all_coords_vec();
791 }
792
793 };
794
795
796 /// compute the derivative of U1 wrt the displacement of atom A, coord axis
797 class U1X_functor : public FunctionFunctorInterface<double,3> {
800 const int U1axis; /// U1x/U1y/U1z potential?
801 const int derivativeaxis; /// direction of the derivative operator
802 public:
804 const int U1axis, const int daxis) : ncf(ncf), thisatom(atom1),
805 U1axis(U1axis), derivativeaxis(daxis) {
806 double lo=1.0/thisatom.q;
808 }
809
810 U1X_functor(const NuclearCorrelationFactor* ncf, const int iatom,
811 const int U1axis, const int daxis) : ncf(ncf),
812 thisatom(ncf->molecule.get_atom(iatom)),
813 U1axis(U1axis), derivativeaxis(daxis) {
814 double lo=1.0/thisatom.q;
816 }
817
818 double operator()(const coord_3d& xyz) const {
819 const coord_3d vr1A=xyz-thisatom.get_coords();
820 const double r=vr1A.normf();
821 const double& Z=thisatom.q;
822 const double S1=ncf->Sr_div_S(r,Z);
823 const double S2=ncf->Srr_div_S(r,Z);
824
825 // note the sign change smoothed_unitvec due to the
826 // change in the derivative variable x: electronic -> nuclear
827 const double drhodx=-ncf->smoothed_unitvec(vr1A)[derivativeaxis];
828 return drhodx*(S2-S1*S1)*ncf->smoothed_unitvec(vr1A)[U1axis]
829 -S1*(ncf->dsmoothed_unitvec(vr1A,derivativeaxis)[U1axis]);
830 }
831
832 std::vector<coord_3d> special_points() const {
833 std::vector< madness::Vector<double,3> > c(1);
834 c[0][0]=thisatom.x;
835 c[0][1]=thisatom.y;
836 c[0][2]=thisatom.z;
837 return c;
838 }
839
840 };
841
842
843 /// compute the derivative of U2 wrt the displacement of atom A
844 class U2X_functor : public FunctionFunctorInterface<double,3> {
846 const int iatom;
847 const int axis;
848 public:
849 U2X_functor(const NuclearCorrelationFactor* ncf, const int& atom1,
850 const int axis) : ncf(ncf), iatom(atom1), axis(axis) {
851 const Atom& atom=ncf->molecule.get_atom(iatom);
852 double lo=1.0/atom.q;
854 }
855
856 double operator()(const coord_3d& xyz) const {
857 const Atom& atom=ncf->molecule.get_atom(iatom);
858 const coord_3d vr1A=xyz-atom.get_coords();
859 const double r=vr1A.normf();
860 const double& Z=atom.q;
861 const double rcut=ncf->molecule.get_rcut()[iatom];
862
863 // note the sign change due to the change in the derivative
864 // variable x: electronic -> nuclear in drho/dx
865 const double drhodx=-ncf->smoothed_unitvec(vr1A)[axis];
866 return drhodx*ncf->U2X_spherical(r,Z,rcut);
867 }
868
869 std::vector<coord_3d> special_points() const {
870 std::vector< madness::Vector<double,3> > c(1);
871 const Atom& atom=ncf->molecule.get_atom(iatom);
872 c[0][0]=atom.x;
873 c[0][1]=atom.y;
874 c[0][2]=atom.z;
875 return c;
876 }
877 };
878
879
880 /// compute the derivative of U3 wrt the displacement of atom A, coord axis
881
882 /// \f[
883 /// U_3^{X_A} = -\sum_{B\neq A}\left(\frac{\vec S_A'}{S_A}\right)^X\cdot\left(\frac{\vec S_B'}{S_B}\right)
884 /// \f]
885 /// with
886 /// \f[
887 /// \left(\frac{\vec S_A'}{S_A}\right)^X =
888 /// \frac{\partial \rho}{\partial X}\left(\frac{S''_A}{S_A}
889 /// -\left(\frac{S'_A}{S_A}\right)^2\right)\vec n_{1A}
890 /// + \left(\frac{S'_A}{S_A}\right)\frac{\partial \vec n_{1A}}{\partial X}
891 /// \f]
892 class U3X_functor : public FunctionFunctorInterface<double,3> {
894 const size_t iatom;
895 const int axis;
896 public:
898 const int axis) : ncf(ncf), iatom(iatom), axis(axis) {}
899
900 double operator()(const coord_3d& xyz) const {
901 const Atom& atomA=ncf->molecule.get_atom(iatom);
902 const coord_3d vr1A=xyz-atomA.get_coords();
903 const double r1A=vr1A.normf();
904 const double& ZA=atomA.q;
905
906 double S1A=ncf->Sr_div_S(r1A,ZA);
907 double S2A=ncf->Srr_div_S(r1A,ZA);
908 double termA=S2A-S1A*S1A;
909
910 // unit vector \vec n_A = \vec r_{1A}/r_{1A}
911 const coord_3d nA=ncf->smoothed_unitvec(vr1A);
912 // derivative of the unit vector \frac{\partial \vec n_A}{\partial X}
913 const coord_3d dnA=ncf->dsmoothed_unitvec(vr1A,axis)*(-1.0);
914 // \frac{\partial \rho}{\partial X}
915 const double drhodx=-nA[axis];
916
917 double term=0.0;
918 for (size_t jatom=0; jatom<ncf->molecule.natom(); ++jatom) {
919 if (iatom==jatom) continue; // restricted sum B \neq A
920
921 const Atom& atomB=ncf->molecule.get_atom(jatom);
922 const coord_3d vr1B=xyz-atomB.get_coords();
923 const double r1B=vr1B.normf();
924 const double& ZB=atomB.q;
925
926 double S1B=ncf->Sr_div_S(r1B,ZB);
927 const coord_3d nB=ncf->smoothed_unitvec(vr1B);
928
929 double dot=0.0; // n_A.n_B
930 double ddot=0.0; // n'_A.n_B
931 for (int i=0; i<3; ++i) {
932 ddot+=dnA[i]*nB[i];
933 dot+=nA[i]*nB[i];
934 }
935 term+=(+drhodx*termA*S1B*dot + S1A*S1B*ddot);
936
937 }
938
939 return term;
940 }
941
942 std::vector<coord_3d> special_points() const {
943 return ncf->molecule.get_all_coords_vec();
944 }
945 };
946
947};
948
949
950/// A nuclear correlation factor class
951
952/// The nuclear correlation factor is given by
953/// \[f
954/// R = \prod S_A ; S_A=exp(-Z_A r_{1A}) + ( 1 - exp(-r_{1A}^2) )
955/// \]f
957public:
958 /// ctor
959
960 /// @param[in] world the world
961 /// @param[in] mol molecule with the sites of the nuclei
964
965 if (world.rank()==0) {
966 print("constructed nuclear correlation factor of the form");
967 print(" R = Prod_A S_A");
968 print(" S_A = exp(-Z_A r_{1A}) + (1 - exp(-Z_A^2*r_{1A}^2))");
969 print("with eprec ",mol.get_eprec());
970 print("which is of Gaussian-Slater type\n");
971 }
972
973 }
974
976
977private:
978
979 /// the nuclear correlation factor
980 double S(const double& r, const double& Z) const {
981 const double rho=r*Z;
982 return exp(-rho)+(1.0-exp(-(rho*rho)));
983 }
984
985 /// radial part first derivative of the nuclear correlation factor
986 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
987
988 const double r=sqrt(vr1A[0]*vr1A[0] +
989 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
990
991 const double eA=exp(-Z*r);
992 const double gA=exp(-Z*Z*r*r);
993 coord_3d term=(2.0*gA*Z*Z*vr1A-Z*eA*smoothed_unitvec(vr1A));
994 return term;
995 }
996
997 /// second derivative of the nuclear correlation factor
998
999 /// -1/2 S"/S - Z/r
1000 double Spp_div_S(const double& r, const double& Z) const {
1001 const double rho=Z*r;
1002 if (rho<1.e-4) {
1003 return Z*Z*(-3.5 - 4.0*rho + 6.0*rho*rho + 12.0*rho*rho*rho);
1004 } else {
1005 const double e=exp(-rho);
1006 const double g=exp(-rho*rho);
1007 const double term1=-Z/r*(1.0-g);
1008 const double term2=-g*Z*Z*(3.0-2.0*Z*Z*r*r) - Z*Z/2.0*e;
1009 const double S_inv=exp(-rho)+(1.0-exp(-(rho*rho)));
1010 return (term1+term2)/S_inv;
1011 }
1012 }
1013
1014 double Sr_div_S(const double& r, const double& Z) const {
1015 const double Zr=r*Z;
1016 const double eA=exp(-Zr);
1017 const double gA=exp(-Zr*Zr);
1018 const double num=Z*(2.0*Zr*gA-eA);
1019 const double denom=1.0+eA-gA;
1020 return num/denom;
1021 }
1022
1023 double Srr_div_S(const double& r, const double& Z) const {
1024 const double Zr=r*Z;
1025 const double eA=exp(-Zr);
1026 const double gA=exp(-Zr*Zr);
1027 const double num=Z*Z*(eA+gA*(2.0-4.0*Zr*Zr));
1028 const double denom=1.0+eA-gA;
1029 return num/denom;
1030 }
1031
1032 double Srrr_div_S(const double& r, const double& Z) const {
1033 const double Zr=r*Z;
1034 const double eA=exp(-Zr);
1035 const double gA=exp(-Zr*Zr);
1036 const double num=Z*Z*Z*(-eA - 12.0*gA*Zr + 8.0*gA*Zr*Zr*Zr);
1037 const double denom=1.0+eA-gA;
1038 return num/denom;
1039
1040 }
1041
1042 /// derivative of the U2 potential wrt X (scalar part)
1043
1044 /// with
1045 /// \f[
1046 /// \rho = \left| \vec r- \vec R_A \right|
1047 /// \f]
1048 /// returns the term in the parenthesis without the the derivative of rho
1049 /// \f[
1050 /// \frac{\partial U_2}{\partial X_A} = \frac{\partial \rho}{\partial X}
1051 /// \left(-\frac{1}{2}\frac{S''' S - S'' S'}{S^2} + \frac{1}{\rho^2}\frac{S'}{S}
1052 /// - \frac{1}{\rho} \frac{S''S - S'^2}{S^2} + \frac{Z_A}{\rho^2}\right)
1053 /// \f]
1054 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
1055
1056 double result=0.0;
1057 if (r*Z<1.e-4) {
1058 const double ZZ=Z*Z;
1059 const double ZZZ=ZZ*Z;
1060 const double Z4=ZZ*ZZ;
1061 const double r0=-4.0*ZZZ;
1062 const double r1=12.0*Z4;
1063 const double r2=36*Z4*Z;
1064 const double r3=-67.0/6.0*Z4*ZZ;
1065 result=(r0 + r*r1 + r*r*r2 + r*r*r*r3);
1066
1067 } else {
1068 const double S1=Sr_div_S(r,Z);
1069 const double S2=Srr_div_S(r,Z);
1070 const double S3=Srrr_div_S(r,Z);
1071 const double term1=-0.5*(S3-S1*S2);
1072 const double term2=(S1+Z)/(r*r);
1073 const double term3=(S2-S1*S1)/r;
1074 result=term1+term2-term3;
1075 }
1076 return result;
1077 }
1078
1079
1080};
1081
1082/// A nuclear correlation factor class
1083
1084/// The nuclear correlation factor is given by
1085/// \[f
1086/// R = \prod S_A ; S_A=exp(-Z_A r_{1A}) + ( 1 - exp(-r_{1A}^2) )
1087/// \]f
1089public:
1090 /// ctor
1091
1092 /// @param[in] world the world
1093 /// @param[in] mol molecule with the sites of the nuclei
1094 GradientalGaussSlater(World& world, const Molecule& mol, const double a)
1096
1097 if (world.rank()==0) {
1098 print("constructed nuclear correlation factor of the form");
1099 print(" R = Prod_A S_A");
1100 print(" S_A = 1/sqrt{Z} exp(-Z_A r_{1A}) + (1 - exp(-a^2*Z_A^2*r_{1A}^2))");
1101 print(" a = ",a);
1102 print("with eprec ",mol.get_eprec());
1103 print("which is of Gradiental Gaussian-Slater type\n");
1104 }
1105
1106 }
1107
1109
1110private:
1111
1112 const double a;
1113
1114 /// the nuclear correlation factor
1115 double S(const double& r, const double& Z) const {
1116 const double rho=r*Z;
1117 return 1/sqrt(Z) * exp(-rho)+(1.0-exp(-(a*a*rho*rho)));
1118 }
1119
1120 /// radial part first derivative of the nuclear correlation factor
1121 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
1122
1123 const double r=sqrt(vr1A[0]*vr1A[0] +
1124 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
1125
1126 const double rho=Z*r;
1127 const double sqrtz=sqrt(Z);
1128 const double term=-exp(-rho)*sqrtz + 2.0*a*a*exp(-a*a*rho*rho)*Z*rho;
1129 return term*smoothed_unitvec(vr1A);
1130 }
1131
1132 /// second derivative of the nuclear correlation factor
1133
1134 /// -1/2 S"/S - Z/r
1135 double Spp_div_S(const double& r, const double& Z) const {
1136 const double rho=Z*r;
1137 const double sqrtz=sqrt(Z);
1138 if (rho<1.e-4) {
1139 const double zfivehalf=Z*Z*sqrtz;
1140 const double a2=a*a;
1141 const double a4=a2*a2;
1142 return -0.5*Z*Z
1143 - 3. *a2 * zfivehalf
1144 - 4.* a2 *rho* zfivehalf
1145 - 2. *a2 * rho*rho*zfivehalf
1146 + 5. *a4 *rho*rho*zfivehalf
1147 + 3. *a4 *rho*rho*Z*Z*Z
1148 -0.5 *a2 *rho*rho*rho*zfivehalf
1149 +5.5 *a4 *rho*rho*rho*zfivehalf
1150 +7. *a4 *rho*rho*rho*Z*Z*Z;
1151 } else {
1152 const double e=exp(-rho);
1153 const double g=exp(-a*a*rho*rho);
1154 const double poly=(2.0-6.0*a*a*rho + 4.0*a*a*a*a*rho*rho*rho);
1155 const double num=Z*(-2.0 - e*r*sqrtz + g*poly);
1156 const double denom=2.0*r*(1.0-g+e/sqrtz);
1157 return num/denom;
1158 }
1159 }
1160
1161 double Sr_div_S(const double& r, const double& Z) const {
1162 const double rZ=r*Z;
1163 const double e=exp(-rZ);
1164 const double g=exp(-a*a*rZ*rZ);
1165 const double sqrtz=sqrt(Z);
1166 const double num=-sqrtz*e + 2.0*a*a*g*Z*rZ;
1167 const double denom=1.0-g+e/sqrtz;
1168 return num/denom;
1169 }
1170
1171 double Srr_div_S(const double& r, const double& Z) const {
1172 const double rZ=r*Z;
1173 const double e=exp(-rZ);
1174 const double g=exp(-a*a*rZ*rZ);
1175 const double sqrtz=sqrt(Z);
1176 const double num=e*Z*sqrtz + g*(2.0*a*a - 4.0*power<4>(a)*rZ*rZ)*Z*Z;
1177 const double denom=1.0-g+e/sqrtz;
1178 return num/denom;
1179 }
1180
1181 double Srrr_div_S(const double& r, const double& Z) const {
1182 const double rZ=r*Z;
1183 const double e=exp(-rZ);
1184 const double g=exp(-a*a*rZ*rZ);
1185 const double sqrtz=sqrt(Z);
1186 const double num=e*power<3>(Z) + (12.0*power<4>(a)*g*rZ
1187 -8.0*power<6>(a)*g*power<3>(rZ))*sqrtz*power<3>(Z);
1188 const double denom=e+sqrtz-g*sqrtz;
1189 return -num/denom;
1190 }
1191
1192 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
1193
1194 double result=0.0;
1195 if (r*Z<1.e-4) {
1196 const double sqrtz=sqrt(Z);
1197 const double Z2=Z*Z;
1198 const double Z4=Z2*Z2;
1199 const double Z5=Z4*Z;
1200 const double Z6=Z5*Z;
1201 const double Z7=Z6*Z;
1202 const double a2=a*a;
1203 const double a4=a2*a2;
1204
1205 const double r0=-4.* a2* sqrt(Z7);
1206 const double r1=2.* (-2.* a2* Z*sqrt(Z7)+ 5.* a4* Z*sqrt(Z7) + 3.* a4 *Z5) *r;
1207 const double r2=1.5 * (-a2* sqrtz*Z5 + 11.* a4* sqrtz*Z5 + 14.*a4* Z6)* r*r;
1208 const double r3=1./6.* (-a2* sqrtz*Z6 + 66.* a4*sqrtz*Z6 - 84.* a2*a4* sqrtz*Z6 +
1209 180. *a4* Z7 - 156.*a2*a4* Z7 - 72.* a2*a4*sqrtz*Z7) *r*r*r;
1210 result=(r0 + r1 + r2 + r3);
1211
1212 } else {
1213 const double S1=Sr_div_S(r,Z);
1214 const double S2=Srr_div_S(r,Z);
1215 const double S3=Srrr_div_S(r,Z);
1216 const double term1=-0.5*(S3-S1*S2);
1217 const double term2=(S1+Z)/(r*r);
1218 const double term3=(S2-S1*S1)/r;
1219 result=term1+term2-term3;
1220 }
1221 return result;
1222 }
1223};
1224
1225
1226/// A nuclear correlation factor class
1227
1228/// The nuclear correlation factor is given by
1229/// \[f
1230/// R = \prod S_A ; S_A= -Z_A r_{1A} exp(-Z_A r_{1A}) + 1
1231/// \]f
1233public:
1234 /// ctor
1235
1236 /// @param[in] world the world
1237 /// @param[in] mol molecule with the sites of the nuclei
1238 LinearSlater(World& world, const Molecule& mol, const double a)
1239 : NuclearCorrelationFactor(world,mol), a_(1.0) {
1240
1241 if (a!=0.0) a_=a;
1242
1243 if (world.rank()==0) {
1244 print("constructed nuclear correlation factor of the form");
1245 print(" S_A = -Z_A r_{1A} exp(-Z_A r_{1A}) + 1");
1246 print(" a = ",a_);
1247 print("with eprec ",mol.get_eprec());
1248 print("which is of linear Slater type\n");
1249 }
1250 }
1251
1253
1254private:
1255
1256 /// the length scale parameter a
1257 double a_;
1258
1259 double a_param() const {return 1.0;}
1260
1261 /// the nuclear correlation factor
1262 double S(const double& r, const double& Z) const {
1263 const double rho=r*Z;
1264 const double b=a_param();
1265 return (-rho)*exp(-b*rho)+1.0;
1266 }
1267
1268 /// radial part first derivative of the nuclear correlation factor
1269 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
1270
1271 const double b=a_param();
1272 const double r=sqrt(vr1A[0]*vr1A[0] +
1273 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
1274
1275 const double ebrz=exp(-b*r*Z);
1276 const coord_3d term=Z*ebrz*(b*Z*(vr1A) - smoothed_unitvec(vr1A));
1277 return term;
1278 }
1279
1280 /// second derivative of the nuclear correlation factor
1281
1282 /// -1/2 S"/S - Z/r
1283 double Spp_div_S(const double& r, const double& Z) const {
1284
1285 const double b=a_param();
1286 const double rho=Z*r;
1287 if (rho<1.e-4) {
1288 const double O0=1.0- 3.0* b;
1289 const double O1=Z - 4.0*b*Z + 3.0*b*b*Z;
1290 const double O2=Z*Z - 5.0*b*Z*Z + 6.5*b*b*Z*Z - 5.0/3.0*b*b*b*Z*Z;
1291 return Z*Z*(O0 + O1*r + O2*r*r);
1292
1293 } else {
1294 const double ebrz=exp(-b*rho);
1295 const double num=Z* (ebrz - 1.0 + 0.5*ebrz*rho* (2.0 + b*(b*rho-4.0)));
1296 const double denom=r*(rho*ebrz-1.0);
1297 return -num/denom;
1298 }
1299 }
1300
1301 double Sr_div_S(const double& r, const double& Z) const {
1302 const double& a=a_param();
1303 const double earz=exp(-a*r*Z);
1304 return Z*earz*(a*r*Z-1.0)/(1.0-r*Z*earz);
1305 }
1306
1307 double Srr_div_S(const double& r, const double& Z) const {
1308 const double& a=a_param();
1309 const double earz=exp(-a*r*Z);
1310 return a*Z*Z*earz*(a*r*Z-2.0)/(-1.0+r*Z*earz);
1311 }
1312
1313 double Srrr_div_S(const double& r, const double& Z) const {
1314 const double& a=a_param();
1315 const double earz=exp(-a*r*Z);
1316 return a*a*Z*Z*Z*earz*(a*r*Z-3.0)/(1.0-r*Z*earz);
1317 }
1318
1319};
1320
1321
1322/// A nuclear correlation factor class
1324public:
1325 /// ctor
1326
1327 /// @param[in] world the world
1328 /// @param[in] mol molecule with the sites of the nuclei
1329 Slater(World& world, const Molecule& mol, const double a)
1330 : NuclearCorrelationFactor(world,mol), a_(1.5) {
1331
1332 if (a!=0.0) a_=a;
1333 eprec_=mol.get_eprec();
1334
1335 if (world.rank()==0) {
1336 print("\nconstructed nuclear correlation factor of the form");
1337 print(" S_A = 1/(a-1) exp(-a Z_A r_{1A}) + 1");
1338 print(" a = ",a_);
1339 print("with eprec ",eprec_);
1340 print("which is of Slater type\n");
1341 }
1342 }
1343
1345
1346private:
1347
1348 /// the length scale parameter
1349 double a_;
1350 double eprec_;
1351
1352 double a_param() const {return a_;}
1353 double eprec_param() const {return eprec_;}
1354
1355 /// first derivative of the correlation factor wrt (r-R_A)
1356
1357 /// \f[
1358 /// Sr_div_S = \frac{1}{S(r)}\frac{\partial S(r)}{\partial r}
1359 /// \f]
1360 double Sr_div_S(const double& r, const double& Z) const {
1361 const double& a=a_param();
1362 return -a*Z/(1.0+(a-1.0)*exp(a*r*Z));
1363 }
1364
1365 /// second derivative of the correlation factor wrt (r-R_A)
1366
1367 /// \f[
1368 /// result = \frac{1}{S(r)}\frac{\partial^2 S(r)}{\partial r^2}
1369 /// \f]
1370 double Srr_div_S(const double& r, const double& Z) const {
1371 const double& a=a_param();
1372 const double aZ=a*Z;
1373 return aZ*aZ/(1.0+(a-1.0)*exp(r*aZ));
1374 }
1375
1376 /// third derivative of the correlation factor wrt (r-R_A)
1377
1378 /// \f[
1379 /// result = \frac{1}{S(r)}\frac{\partial^3 S(r)}{\partial r^3}
1380 /// \f]
1381 double Srrr_div_S(const double& r, const double& Z) const {
1382 const double& a=a_param();
1383 const double aZ=a*Z;
1384 return -aZ*aZ*aZ/(1.0+(a-1.0)*exp(r*aZ));
1385 }
1386
1387 /// the nuclear correlation factor
1388 double S(const double& r, const double& Z) const {
1389 const double a=a_param();
1390 //const double eprec=eprec_param();
1391 return 1.0+1.0/(a-1.0) * exp(-a*Z*r);
1392
1393// return 1.0 + 0.5/(a-1.0) *
1394// (exp(-a*r*Z + 0.5*a*a*Z*Z*eprec) * erfc((-r+a*eprec*Z)/sqrt(2*eprec))
1395// + exp(-a*r*Z + 0.5*a*Z*(4.0*r+a*eprec*Z)) * erfc((r+a*eprec*Z)/sqrt(2*eprec)));
1396 }
1397
1398 /// radial part first derivative of the nuclear correlation factor
1399 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
1400 const double a=a_param();
1401 const double r=vr1A.normf();
1402 return -(a*exp(-a*Z*r)*Z)/(a-1.0)*smoothed_unitvec(vr1A);
1403 }
1404
1405 /// second derivative of the nuclear correlation factor
1406 double Spp_div_S(const double& r, const double& Z) const {
1407 const double a=a_param();
1408
1409 if (r*Z<1.e-4) {
1410 const double O0=1.0-(1.5*a);
1411 const double O1=(a-1.0)*(a-1.0)*Z;
1412 const double O2=(1.0/12.0 * (a-1.0)*(12.0+a*(5*a-18.0)))*Z*Z;
1413 return Z*Z*(O0 + O1*r + O2*r*r);
1414
1415 } else {
1416 const double earz=exp(-a*r*Z);
1417 const double num=Z*(-earz + a*earz - (a-1.0) - 0.5*a*a*r*Z*earz);
1418 const double denom=(r*earz + (a-1.0) * r);
1419 return num/denom;
1420 }
1421 }
1422
1423
1424 /// derivative of the U2 potential wrt X (scalar part)
1425
1426 /// with
1427 /// \f[
1428 /// \rho = \left| \vec r- \vec R_A \right|
1429 /// \f]
1430 /// returns the term in the parenthesis without the the derivative of rho
1431 /// \f[
1432 /// \frac{\partial U_2}{\partial X_A} = \frac{\partial \rho}{\partial X}
1433 /// \left(-\frac{1}{2}\frac{S''' S - S'' S'}{S^2} + \frac{1}{\rho^2}\frac{S'}{S}
1434 /// - \frac{1}{\rho} \frac{S''S - S'^2}{S^2} + \frac{Z_A}{\rho^2}\right)
1435 /// \f]
1436 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
1437 const double a=a_param();
1438
1439 double result=0.0;
1440 if (r*Z<1.e-4) {
1441 const double ZZ=Z*Z;
1442 const double ZZZ=ZZ*Z;
1443 const double a2=a*a;
1444 const double a4=a2*a2;
1445 const double r0=ZZZ*(1. - 2.* a + a2);
1446 const double r1=ZZ*ZZ/6.* (12.0 - 30.* a + 23. *a2 - 5.*a*a2);
1447 const double r2=1./8.*ZZ*ZZZ* (24. - 72.*a + 74.*a2 - 29.*a2*a + 3.*a4);
1448 const double r3=1./60.*ZZZ*ZZZ* (240. - 840.*a + 1080.*a2 - 610.*a2*a
1449 + 137.*a2*a2 - 7.*a4*a);
1450 result=(r0 + r*r1 + r*r*r2 + r*r*r*r3);
1451
1452 } else {
1453 const double S1=Sr_div_S(r,Z);
1454 const double S2=Srr_div_S(r,Z);
1455 const double S3=Srrr_div_S(r,Z);
1456 const double term1=-0.5*(S3-S1*S2);
1457 const double term2=(S1+Z)/(r*r);
1458 const double term3=(S2-S1*S1)/r;
1459 result=term1+term2-term3;
1460 }
1461 return result;
1462 }
1463
1464};
1465
1466
1468public:
1469 /// ctor
1470
1471 /// @param[in] world the world
1472 /// @param[in] mol molecule with the sites of the nuclei
1473 poly4erfc(World& world, const Molecule& mol, const double aa)
1474 : NuclearCorrelationFactor(world,mol), a(1.0) {
1475
1476 if (aa!=0.0) a=aa;
1477 eprec_=mol.get_eprec();
1478
1479 if (world.rank()==0) {
1480 print("\nconstructed nuclear correlation factor of the form");
1481 print(" S_A = 1 + (a0 + a1 arZ + a2 (arZ)^2 + a3 (arZ)^3 + a4 (arZ)^4) erfc(arZ)");
1482 print(" a = ",a);
1483 print("with eprec ",eprec_);
1484 print("which is of poly4erfc type\n");
1485 }
1486 //const double pi32=std::pow(constants::pi,1.5);
1487 //const double sqrtpi=sqrt(constants::pi);
1488 //const double Pi=constants::pi;
1489
1490 if (a==0.5) {
1491 a0=0.5083397721116242769;
1492 a1=-2.4430795355664112811;
1493 a2=3.569312300653802680;
1494 a3=-1.9812471972342746507;
1495 a4=0.3641705622093696564;
1496 } else if (a==1.0) {
1497 a0=0.20265985404508529127;
1498 a1=-0.9739826967339938056;
1499 a2=1.4229779953809877198;
1500 a3=-0.7898639647077711196;
1501 a4=0.14518390461225107425;
1502
1503 } else {
1504 print("invalid parameter a for poly4erfc: only 0.5 and 1.0 implemented");
1505 MADNESS_EXCEPTION("stupid you",1);
1506 }
1507 }
1508
1510
1511private:
1512
1513 /// the length scale parameter
1514 double a;
1515 double a0, a1, a2, a3, a4;
1516 double eprec_;
1517
1518 double a_param() const {return a;}
1519 double eprec_param() const {return eprec_;}
1520
1521 /// first derivative of the correlation factor wrt (r-R_A)
1522
1523 /// \f[
1524 /// Sr_div_S = \frac{1}{S(r)}\frac{\partial S(r)}{\partial r}
1525 /// \f]
1526 double Sr_div_S(const double& r, const double& Z) const {
1527 const double x=r*Z;
1528
1529 double result=0.0;
1530 if (a==0.5) {
1531 if (x<1.0) {
1532 result=(-17.97663543396820624361586474 + x*(32.78290319470982346841067868 +
1533 x*(-18.158574783628271659713638233 +
1534 x*(2.472138984374094343735335913 +
1535 x*(0.5516975358315341628276502285 +
1536 x*(-0.008573693952875097234391220137 +
1537 x*(-0.05596791202351071993748992739 +
1538 (0.002673799219133696315436690424 +
1539 0.0013386538660557369902632289083*x)*x)))))))/
1540 (17.976635433967702922140233083 + x*
1541 (-13.062204300852085089323266568 +
1542 x*(16.397871971437618641239835027 + x*(-5.383337491559214163188757918 + 1.*x))));
1543 } else if (x<2.0) {
1544 result=(-16.53370050883888159389958126 + x*(30.04151304875517461538549269 +
1545 x*(-16.692529697855029750948871179 +
1546 x*(2.50341008323651011875249567 +
1547 x*(0.3106921665634860719234742532 +
1548 x*(0.08721948207311506458903445571 +
1549 x*(-0.10041387133168708232852057948 +
1550 (0.02000987266876476192949541524 - 0.0012508983745483161308604975792*x)*
1551 x)))))))/
1552 (16.532205048243702951212522516 + x*
1553 (-11.89273747187279634240347945 +
1554 x*(15.157537549656745468895369276 + x*(-4.9102960292797655978798640519 + 1.*x))));
1555 } else if (x<5.0) {
1556 result=(-2352.191894900273554810118278 + x*(5782.846962269399174183661793 +
1557 x*(-5653.246084369776756298278851 +
1558 x*(2948.18046377483570925427449 +
1559 x*(-913.4583247839311453090142452 +
1560 x*(174.39391722588915386106331206 +
1561 x*(-20.22035127074332315930567933 +
1562 (1.3107321165711966663791114988 - 0.03655666729452579098523876463*x)*x))
1563 )))))/
1564 (886.4859678528423041797649741 + x*(269.17746130370931387996124706 +
1565 x*(-130.21383548057958115685397713 + x*(37.644499985765056273193347388 + 1.*x))));
1566 } else if (x<10) {
1567 result=(2.2759176275121988686860433041 + x*(-1.8014283464827425541637211503 +
1568 x*(0.60570955276433317373251991152 +
1569 x*(-0.11235368819003308411943926239 +
1570 x*(0.01243211635244600976892077538 +
1571 x*(-0.0008211800260491381826149891865 +
1572 x*(0.000029973534470203049782417744015 +
1573 (-4.6423722605763162431293872646e-7 -
1574 1.3224615425412157194986675329e-10*x)*x)))))))/
1575 (1039.800013929971888016478838 + x*(-702.5378531848183775210948787 +
1576 x*(183.17476380259879599459789974 + x*(-21.74106003575315304073197254 + 1.*x))));
1577 } else {
1578 result=0.0;
1579 }
1580 } else if (a==1.0) {
1581 if (x<1.0) {
1582 result=(-1.6046958001953006847027538457 + x*(5.948945186367159977879486279 +
1583 x*(-6.884321742840291285733040882 +
1584 x*(2.296896506418919905405783368 +
1585 x*(0.616939354810622973212914089 +
1586 x*(-0.13679830198890803519235207564 +
1587 x*(-0.3356576872066501398893403439 +
1588 (0.14876798925798674488426727928 - 0.016049886728185297028755535226*x)*x
1589 )))))))/
1590 (1.6046957999975126909719683196 + x*
1591 (-0.8234878506688316215458988304 +
1592 x*(3.4607641859010551501639903314 + x*(-1.8955210085531557309609670978 + 1.*x))));
1593 } else if (x<2.0) {
1594 result=(-7.143856421301985019580778813 + x*(33.35568129248075086686087865 +
1595 x*(-60.0412343766569246898209957 +
1596 x*(55.46407913315151830247939138 +
1597 x*(-28.7874770749840240158264326 +
1598 x*(8.379317837934083469035061852 +
1599 x*(-1.2103278392957399092107317741 +
1600 (0.04275186003977071121074860478 + 0.005247730112126140726731063205*x)*x
1601 )))))))/
1602 (5.4509691924562038998044993827 + x*
1603 (-2.2732931206867811068721699796 +
1604 x*(4.1530634219989859344450742259 + x*(-2.0183662125874366044951391259 + 1.*x))));
1605 } else if (x<5.0) {
1606 result=(-0.4869290414611847276899694883 + x*(1.0513417728218375522016338562 +
1607 x*(-0.9694851629317255156038942437 +
1608 x*(0.5007460889402078102011673774 +
1609 x*(-0.15897436623286639052954575417 +
1610 x*(0.03185042477631079356985872518 +
1611 x*(-0.003940899912654543218183049969 +
1612 (0.0002757919409481032696079686825 -
1613 8.369138363906178282041501561e-6*x)*x)))))))/
1614 (30.81121613134246115634780413 + x*(-49.989974505724725146933397436 +
1615 x*(31.455643953462635691568128729 + x*(-8.992097794824270871044305786 + 1.*x))));
1616 } else {
1617 result=0.0;
1618 }
1619 }
1620 return result*Z;
1621 }
1622
1623 /// second derivative of the correlation factor wrt (r-R_A)
1624
1625 /// \f[
1626 /// result = \frac{1}{S(r)}\frac{\partial^2 S(r)}{\partial r^2}
1627 /// \f]
1628 double Srr_div_S(const double& r, const double& Z) const {
1629 MADNESS_EXCEPTION("no Srr_div_S in Slater2 yet",0);
1630 return 0.0;
1631 }
1632
1633 /// third derivative of the correlation factor wrt (r-R_A)
1634
1635 /// \f[
1636 /// result = \frac{1}{S(r)}\frac{\partial^3 S(r)}{\partial r^3}
1637 /// \f]
1638 double Srrr_div_S(const double& r, const double& Z) const {
1639 MADNESS_EXCEPTION("no Srrr_div_S in Slater2 yet",0);
1640 return 0.0;
1641 }
1642
1643 /// the nuclear correlation factor
1644 double S(const double& r, const double& Z) const {
1645 const double arZ=a*r*Z;
1646 const double arZ2=a*a*r*r*Z*Z;
1647
1648 return 1.0 + (a0 +a1* arZ+ a2 * arZ2 + a3*arZ*arZ2 + + a4*arZ2*arZ2) *erfc( a*r*Z);
1649 }
1650
1651 /// radial part first derivative of the nuclear correlation factor
1652 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
1653 MADNESS_EXCEPTION("no Sp in Slater2 yet",0);
1654 return smoothed_unitvec(vr1A);;
1655 }
1656
1657 /// second derivative of the nuclear correlation factor
1658 double Spp_div_S(const double& r, const double& Z) const {
1659 const double x=r*Z;
1660 double result=0.0;
1661 if (a==0.5) {
1662 if (x<1.0) {
1663 result=(-37.75186842343823465059 + x*(21.3476988467348903615 +
1664 x*(0.014608707333424946750026 +
1665 x*(-3.704945314726273312722 +
1666 x*(0.08808104845382944292252 +
1667 x*(0.3362428305409206967066 +
1668 x*(-0.02288039625102549092766 +
1669 (-0.017240056622850307571001 + 0.002412740490618117536527*x)*x)))))))/
1670 (17.595611361287293183447 + x*(-12.421065066789808273295 +
1671 x*(15.957564552156320207938 + x*(-5.0239100389132464317907 + 1.*x))));
1672 } else if (x<2.0) {
1673 result=(-35.46082798344982189328 + x*(20.3565414350879797493 +
1674 x*(-0.3046488975234456455871 +
1675 x*(-3.298067641731523613442 +
1676 x*(-0.10712268902574947466554 +
1677 x*(0.4829111116912435121877 +
1678 x*(-0.10089023589818206760275 +
1679 (0.0013899828749998182176948 + 0.0008305150868988335610678*x)*x)))))))/
1680 (16.526499668833810286111 + x*(-11.79820182874024593006 +
1681 x*(15.115407083582433322342 + x*(-4.8449266197426825174319 + 1.*x))));
1682 } else if (x<5.0) {
1683 result=(-414.5311559516023264104 + x*(615.8720414166440747539 +
1684 x*(-422.5938440932793094888 + x*
1685 (159.72497494584873155352 +
1686 x*(-35.6790348104188081907 +
1687 x*(4.658777521872728594702 +
1688 x*(-0.328305094433759490678 +
1689 (0.009162754689309596905172 + 0.00005047926659010755662873*x)*x)))))))/
1690 (112.7044543272820830484 + x*(-56.072518762714727894479 +
1691 x*(28.188843903409059322224 + x*(-6.519554545057610040741 + 1.*x))));
1692 } else if (x<10.0) {
1693 result=(-146.68256559112012287314 + x*(188.52807059353309478385 +
1694 x*(-82.07176992590032524431 + x*
1695 (18.107697718347802322776 +
1696 x*(-2.3221393933622638466979 +
1697 x*(0.18681223946803275939642 +
1698 x*(-0.009601011427143648072501 +
1699 (0.00028623295732553770894583 - 3.77364531155112782235e-6*x)*x)))))))/
1700 (633.1319105785227043552 + x*(-574.29230199331494406798 +
1701 x*(171.872612865808639376 + x*(-21.906495121260483674385 + 1.*x))));
1702 } else {
1703 result=-1.0/x;
1704 }
1705 } else if (a==1.0) {
1706 if (x<1.0) {
1707 result=(-8.85288955131414420808 + x*(11.359434597010419928515 +
1708 x*(-2.982094072176256165405 + x*
1709 (-4.924201512880445076103 +
1710 x*(0.6773928043289287907009 +
1711 x*(2.061680233090141287213 +
1712 x*(-0.5528612000728412913713 +
1713 (-0.3024276764350212595121 + 0.10765958646570631264003*x)*x)))))))/
1714 (1.6731803922436094206275 + x*(-0.8242052614481793487834 +
1715 x*(3.5912910648514747558597 + x*(-1.9146644266104277222142 + 1.*x))));
1716 } else if (x<2.0) {
1717 result=(-8.538839434207355205544 + x*(-37.85611260277589742674 +
1718 x*(155.08466711228234382211 + x*
1719 (-241.1247881821171613212 +
1720 x*(199.33293375918859716585 +
1721 x*(-96.80750216221383928113 +
1722 x*(27.71704080319043943288 +
1723 (-4.354251431065185902435 + 0.2906781051124817624589*x)*x)))))))/
1724 (2.2233877901091612794108 + x*(3.2084406367698851844258 +
1725 x*(1.021615116328133792993 + x*(-0.9819667717342250528772 + 1.*x))));
1726 } else if (x<5.0) {
1727 result=(-7.338671998425412491091 + x*(18.593867285988629508481 +
1728 x*(-19.15344657249203844577 + x*
1729 (10.072359942912659732126 +
1730 x*(-2.99771628023376895151 +
1731 x*(0.5324416109232007541639 +
1732 x*(-0.05993976172902101376555 +
1733 (0.003887110757665478374745 - 0.00011079212872583089652945*x)*x)))))))/
1734 (13.273604111590967002999 + x*(-28.012330064250556933961 +
1735 x*(22.161347591665727407914 + x*(-7.617582259533338164664 + 1.*x))));
1736 } else if (x<10) {
1737 result=(132.78397650676876308801 + x*(-78.01898251977413923271 +
1738 x*(15.292039515801252996504 + x*
1739 (-1.000000154745351328125 +
1740 x*(1.950300262619412492847e-8 +
1741 x*(-1.6337541591423364318692e-9 +
1742 x*(8.771557330523968791519e-11 +
1743 (-2.7388800346031121159372e-12 + 3.7894572289998844980568e-14*x)*x))))))
1744 )/(-4.724326520908917985827e-6 + x*
1745 (-132.78397108375581258096 + x*(78.01897976124964688489 +
1746 x*(-15.292038699709498552356 + 1.*x))));
1747 } else {
1748 result=-1.0/x;
1749 }
1750 }
1751 return result*Z*Z;
1752 }
1753
1754
1755 /// derivative of the U2 potential wrt X (scalar part)
1756
1757 /// with
1758 /// \f[
1759 /// \rho = \left| \vec r- \vec R_A \right|
1760 /// \f]
1761 /// returns the term in the parenthesis without the the derivative of rho
1762 /// \f[
1763 /// \frac{\partial U_2}{\partial X_A} = \frac{\partial \rho}{\partial X}
1764 /// \left(-\frac{1}{2}\frac{S''' S - S'' S'}{S^2} + \frac{1}{\rho^2}\frac{S'}{S}
1765 /// - \frac{1}{\rho} \frac{S''S - S'^2}{S^2} + \frac{Z_A}{\rho^2}\right)
1766 /// \f]
1767 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
1768 MADNESS_EXCEPTION("no U2X_spherical in Slater2 yet",0);
1769 return 0.0;
1770 }
1771
1772};
1773
1774
1775
1776/// A nuclear correlation factor class
1777
1778/// should reduce to quartic for N=4
1779/// @tparam N the exponent of the polynomial
1780template<std::size_t N>
1782public:
1783 /// ctor
1784
1785 /// @param[in] world the world
1786 /// @param[in] mol molecule with the sites of the nuclei
1787 Polynomial(World& world, const Molecule& mol, const double a)
1789
1790 /// length scale parameter a, default chosen that linear terms in U2 vanish
1791 a_=(2. + (-2. + sqrt(-1. + N))*N)/(-2. + N);
1792
1793 if (a!=0.0) a_=a;
1794
1795 if (world.rank()==0) {
1796 print("constructed nuclear correlation factor of the form");
1797 print(" R = Prod_A S_A");
1798 print(" S_A = 1 + a (r/b -1)^N if r<b, with b= (N*a)/((1+a) Z)");
1799 print(" = 1 else ");
1800 print("with eprec ",mol.get_eprec());
1801 print("which is of polynomial type with exponent N = ",N);
1802 }
1803 }
1804
1806
1807private:
1808
1809 /// length scale parameter a, default chosen that linear terms in U2 vanish
1810 double a_;
1811
1812 double a_param() const {return a_;}
1813
1814 /// the cutoff
1815 static double b_param(const double& a) {return N*a/(1.0+a);}
1816
1817 /// the nuclear correlation factor
1818 double S(const double& r, const double& Z) const {
1819
1820 const double rho=r*Z;
1821 const double a=Polynomial<N>::a_param();
1822 const double b=Polynomial<N>::b_param(a);
1823
1824 if (rho<b) {
1825 const double arg=-1.0 + rho/b;
1826 return 1.0 + power<N>(-1.0) * a*power<N>(arg);
1827 } else {
1828 return 1.0;
1829 }
1830
1831 }
1832
1833 /// radial part first derivative of the nuclear correlation factor
1834 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
1835
1836 const double r=vr1A.normf();
1837 const double rho=r*Z;
1838 const double a=Polynomial<N>::a_param();
1839 const double b=Polynomial<N>::b_param(a);
1840
1841 if (rho<b) {
1842 return power<N>(-1.)*(1.+a)* Z* power<N-1>(-1.+rho/b)*smoothed_unitvec(vr1A);
1843 }
1844 return coord_3d(0.0);
1845 }
1846
1847 /// second derivative of the nuclear correlation factor
1848
1849 /// -1/2 S"/S - Z/r
1850 double Spp_div_S(const double& r, const double& Z) const {
1851
1852 const double rho=r*Z;
1853 const double a=Polynomial<N>::a_param();
1854 const double b=Polynomial<N>::b_param(a);
1855
1856 if (rho<1.e-6) {
1857 const double ap1=1.0+a;
1858 const double c0=((3. *(1. + a) - (3. + a) * N))/(2.* a*N);
1859 const double c1=((2.* ap1*ap1 - ap1* (3. + a)*N + N*N)*Z)/(a*a*N*N);
1860 const double c2=((30.*ap1*ap1*ap1- ap1*ap1* (55 + 18* a)*N +
1861 30 *ap1 *N*N + (-5 + a* (8 + a)) *N*N*N)* Z*Z)/(12 *a*a*a*N*N*N);
1862 return Z*Z*(c0 + c1*r + c2*r*r);
1863
1864 } else if (rho<b) {
1865
1866 const double num=Z* (2 + (power<N>(-1)* a* power<N>(-1 + rho/b)
1867 * (-2 *a*N*N + (1 + a) *N* (1 + a *(-3 + N) + N)* rho +
1868 2 *(1 + a)*(1+a)* rho*rho))/power<2>(a* N - (1 + a)*rho));
1869
1870 const double denom=2.* (r + power<N>(-1) *a* r* power<N>(-1 + rho/b));
1871 return -num/denom;
1872
1873 } else {
1874 return -Z*Z/rho;
1875 }
1876 }
1877
1878 double Sr_div_S(const double& r, const double& Z) const {
1879 const double rho=r*Z;
1880 const double a=Polynomial<N>::a_param();
1881 const double b=Polynomial<N>::b_param(a);
1882
1883 if (rho<b) {
1884 const double negn= power<N>(-1.0);
1885 const double num=(negn*(1 + a)*Z*power<N-1>(-1 + ((1 + a)*r*Z)/(a*N)));
1886 const double denom=(1 + negn*a*power<N>(-1 + ((1 + a)*r*Z)/(a*N)));
1887 return num/denom;
1888 } else {
1889 return 0.0;
1890 }
1891
1892 }
1893
1894 double Srr_div_S(const double& r, const double& Z) const {
1895 const double rho=r*Z;
1896 const double a=Polynomial<N>::a_param();
1897 const double b=Polynomial<N>::b_param(a);
1898
1899 if (rho<b) {
1900 const double negn= power<N>(-1.0);
1901 return (negn*power<2>(1 + a)*(-1 + N)*power<2>(Z)*power<N-2>(-1 + ((1 + a)*r*Z)/(a*N)))/
1902 (a*N*(1 + negn*a*power<N>(-1 + ((1 + a)*r*Z)/(a*N))));
1903 } else {
1904 return 0.0;
1905 }
1906 }
1907
1908 double Srrr_div_S(const double& r, const double& Z) const {
1909 const double rho=r*Z;
1910 const double a=Polynomial<N>::a_param();
1911 const double b=Polynomial<N>::b_param(a);
1912
1913 if (rho<b) {
1914 const double negn= power<N>(-1.0);
1915 return (negn*power<3>(1 + a)*(-2 + N)*(-1 + N)*power<3>(Z)*power<N-3>(-1 + ((1 + a)*r*Z)/(a*N)))/
1916 (power<2>(a*N)*(1 + negn*a*power<N>(-1 + ((1 + a)*r*Z)/(a*N))));
1917 } else {
1918 return 0.0;
1919 }
1920 }
1921
1922 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
1923 const double a=a_param();
1924 const double aopt=(2. + (-2. + sqrt(-1. + N))*N)/(-2. + N);
1925 if (fabs(a-aopt)>1.e-10) {
1926 MADNESS_EXCEPTION("U2X_spherical for polynomial ncf only with aopt",1);
1927 }
1928
1929 double result=0.0;
1930 if (r*Z<1.e-4) {
1931 const double rn=sqrt(N-1);
1932 const double r0=0.0;
1933 const double r1=((2.*(-8. + 9.*rn) + N*(25. + 10.*rn + N))*r*power<4>(Z))/
1934 (6.*power<2>(-2 + N)*rn);
1935 const double r2=((-4*(17 + 9*rn) + N*(92 + 80*rn +
1936 N*(-29 - 33*rn + N*(4 + 7*rn + N))))*power<5>(Z))/
1937 (8.*power<3>(-2 + N)*(-1 + N)*rn);
1938 result=(r0 + r*r1 + r*r*r2);
1939
1940 } else {
1941 const double S1=Sr_div_S(r,Z);
1942 const double S2=Srr_div_S(r,Z);
1943 const double S3=Srrr_div_S(r,Z);
1944 const double term1=-0.5*(S3-S1*S2);
1945 const double term2=(S1+Z)/(r*r);
1946 const double term3=(S2-S1*S1)/r;
1947 result=term1+term2-term3;
1948 }
1949 return result;
1950 }
1951
1952};
1953
1955
1956public:
1957 /// ctor
1958
1959 /// @param[in] world the world
1960 /// @param[in] mol molecule with the sites of the nuclei
1962 const std::shared_ptr<PotentialManager> pot, const double fac)
1963 : NuclearCorrelationFactor(world,mol), potentialmanager(pot),
1964 eprec(mol.get_eprec()), fac(fac) {
1965
1966 if (world.rank()==0) {
1967 print("constructed nuclear correlation factor of the form");
1968 print(" R = ",fac);
1969 print("with eprec ",mol.get_eprec());
1970 print("which means it's (nearly) a conventional calculation\n");
1971 }
1972
1973 // add the missing -Z/r part to U2!
1974 }
1975
1976 corrfactype type() const {return None;}
1977
1978 /// return the U2 term of the correlation function
1979
1980 /// overloading to avoid inconsistent state of U2, which needs the
1981 /// nuclear potential
1982 const real_function_3d U2() const {
1983
1984// if (not U2_function.is_initialized()) {
1985 MADNESS_ASSERT(potentialmanager->vnuclear().is_initialized());
1986// }
1987 return potentialmanager->vnuclear();
1988 }
1989
1990 /// apply the regularized potential U_nuc on a given function rhs
1991
1992 /// overload the base class method for efficiency
1994 return (U2()*rhs).truncate();
1995 }
1996
1997
1998private:
1999
2000 /// underlying potential (=molecule)
2001 std::shared_ptr<PotentialManager> potentialmanager;
2002 double eprec;
2003
2004 /// the factor of the correlation factor: R=fac;
2005 const double fac;
2006
2007 double Sr_div_S(const double& r, const double& Z) const {return 0.0;}
2008
2009 double Srr_div_S(const double& r, const double& Z) const {return 0.0;}
2010
2011 double Srrr_div_S(const double& r, const double& Z) const {return 0.0;}
2012
2013 /// the nuclear correlation factor
2014 double S(const double& r, const double& Z) const {
2015 return fac;
2016 }
2017
2018 /// radial part first derivative of the nuclear correlation factor
2019 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
2020 return coord_3d(0.0);
2021 }
2022
2023 /// second derivative of the nuclear correlation factor
2024 double Spp_div_S(const double& r, const double& Z) const {
2025 double rcut= 1.0 / smoothing_parameter(Z, eprec);
2026 return - Z * smoothed_potential(r*rcut)*rcut;
2027 }
2028
2029 double U2X_spherical(const double& r, const double& Z, const double& rcut) const {
2030 // factor -1 from the definition of the dsmoothed_potential as -1/r^2
2031 return -Z*dsmoothed_potential(r * rcut) * (rcut * rcut);
2032 }
2033
2034};
2035
2036
2037/// this ncf has no information about itself, only U2 and U1 assigned
2039
2040public:
2041 /// ctor
2042
2043 /// @param[in] world the world
2044 /// @param[in] mol molecule with the sites of the nuclei
2046 const std::vector<real_function_3d>& U1)
2048
2049 U2_function=U2;
2050 U1_function=U1;
2051
2052 if (world.rank()==0) {
2053 print("constructed ad hoc nuclear correlation factor");
2054 }
2055 }
2056
2057 corrfactype type() const {return Adhoc;}
2058
2059private:
2060
2061 double Sr_div_S(const double& r, const double& Z) const {
2062 MADNESS_EXCEPTION("no Sr_div_S() in AdhocNuclearCorrelationFactor",0);
2063 return 0.0;
2064 }
2065
2066 double Srr_div_S(const double& r, const double& Z) const {
2067 MADNESS_EXCEPTION("no Srr_div_S() in AdhocNuclearCorrelationFactor",0);
2068 return 0.0;
2069 }
2070
2071 double Srrr_div_S(const double& r, const double& Z) const {
2072 MADNESS_EXCEPTION("no Srrr_div_S() in AdhocNuclearCorrelationFactor",0);
2073 return 0.0;
2074 }
2075
2076 /// the nuclear correlation factor
2077 double S(const double& r, const double& Z) const {
2078 MADNESS_EXCEPTION("no S() in AdhocNuclearCorrelationFactor",0);
2079 return 0.0;
2080 }
2081
2082 /// radial part first derivative of the nuclear correlation factor
2083 coord_3d Sp(const coord_3d& vr1A, const double& Z) const {
2084 MADNESS_EXCEPTION("no Sp() in AdhocNuclearCorrelationFactor",0);
2085 return coord_3d(0.0);
2086 }
2087
2088 /// second derivative of the nuclear correlation factor
2089 double Spp_div_S(const double& r, const double& Z) const {
2090 MADNESS_EXCEPTION("no Spp_div_S() in AdhocNuclearCorrelationFactor",0);
2091 return 0.0;
2092 }
2093};
2094
2095
2096std::shared_ptr<NuclearCorrelationFactor>
2098 const Molecule& molecule,
2099 const std::shared_ptr<PotentialManager> pm,
2100 const std::string inputline);
2101
2102std::shared_ptr<NuclearCorrelationFactor>
2104 const Molecule& molecule,
2105 const std::shared_ptr<PotentialManager> pm,
2106 const std::pair<std::string,double>& ncf);
2107
2108}
2109
2110
2111#endif /* NUCLEARCORRELATIONFACTOR_H_ */
Declaration of utility class and functions for atom.
this ncf has no information about itself, only U2 and U1 assigned
Definition correlationfactor.h:2038
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:2089
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2061
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:2083
corrfactype type() const
Definition correlationfactor.h:2057
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2066
AdhocNuclearCorrelationFactor(World &world, const real_function_3d U2, const std::vector< real_function_3d > &U1)
ctor
Definition correlationfactor.h:2045
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:2077
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2071
Definition molecule.h:58
double y
Definition molecule.h:60
double x
Definition molecule.h:60
double z
Definition molecule.h:60
madness::Vector< double, 3 > get_coords() const
Definition molecule.h:99
double q
Coordinates and charge in atomic units.
Definition molecule.h:60
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition derivative.h:266
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
void set_length_scale(double lo)
adapt the special level to resolve the smallest length scale
Definition function_interface.h:80
double thresh() const
Returns value of truncation threshold. No communication.
Definition mra.h:567
void set_thresh(double value, bool fence=true)
Sets the value of the truncation threshold. Optional global fence.
Definition mra.h:577
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
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition mra.h:721
A nuclear correlation factor class.
Definition correlationfactor.h:956
GaussSlater(World &world, const Molecule &mol)
ctor
Definition correlationfactor.h:962
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1023
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:980
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt X (scalar part)
Definition correlationfactor.h:1054
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1032
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1000
corrfactype type() const
Definition correlationfactor.h:975
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:986
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1014
A nuclear correlation factor class.
Definition correlationfactor.h:1088
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1171
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:1115
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:1121
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt nuclear coordinate X (spherical part)
Definition correlationfactor.h:1192
const double a
Definition correlationfactor.h:1112
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1181
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1135
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1161
GradientalGaussSlater(World &world, const Molecule &mol, const double a)
ctor
Definition correlationfactor.h:1094
corrfactype type() const
Definition correlationfactor.h:1108
A nuclear correlation factor class.
Definition correlationfactor.h:1232
LinearSlater(World &world, const Molecule &mol, const double a)
ctor
Definition correlationfactor.h:1238
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1307
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:1269
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:1262
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1301
double a_param() const
Definition correlationfactor.h:1259
double a_
the length scale parameter a
Definition correlationfactor.h:1257
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1313
corrfactype type() const
Definition correlationfactor.h:1252
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1283
Definition molecule.h:124
double atomic_attraction_potential(int iatom, double x, double y, double z) const
nuclear attraction potential for a specific atom in the molecule
Definition molecule.cc:995
double get_eprec() const
Definition molecule.h:425
double nuclear_attraction_potential_derivative(int atom, int axis, double x, double y, double z) const
Definition molecule.cc:1008
compute the derivative of R wrt the displacement of atom A, coord axis
Definition correlationfactor.h:743
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:744
const Atom & thisatom
Definition correlationfactor.h:745
const int derivativeaxis
Definition correlationfactor.h:746
RX_functor(const NuclearCorrelationFactor *ncf, const Atom &atom1, const int daxis, const int exponent)
1 or 2 -> R^X or R^X R
Definition correlationfactor.h:750
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:789
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:763
const int exponent
direction of the derivative operator
Definition correlationfactor.h:747
RX_functor(const NuclearCorrelationFactor *ncf, const int iatom, const int daxis, const int exponent)
Definition correlationfactor.h:756
Definition correlationfactor.h:434
R_functor(const NuclearCorrelationFactor *ncf, const int e=1)
Definition correlationfactor.h:438
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:456
int exponent
Definition correlationfactor.h:436
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:435
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:440
compute the derivative of U1 wrt the displacement of atom A, coord axis
Definition correlationfactor.h:797
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:798
const int U1axis
Definition correlationfactor.h:800
U1X_functor(const NuclearCorrelationFactor *ncf, const int iatom, const int U1axis, const int daxis)
Definition correlationfactor.h:810
U1X_functor(const NuclearCorrelationFactor *ncf, const Atom &atom1, const int U1axis, const int daxis)
direction of the derivative operator
Definition correlationfactor.h:803
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:818
const Atom & thisatom
Definition correlationfactor.h:799
const int derivativeaxis
U1x/U1y/U1z potential?
Definition correlationfactor.h:801
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:832
U1 functor for a specific atom.
Definition correlationfactor.h:497
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:507
const size_t iatom
Definition correlationfactor.h:500
U1_atomic_functor(const NuclearCorrelationFactor *ncf, const size_t atom, const int axis)
Definition correlationfactor.h:504
const int axis
Definition correlationfactor.h:501
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:515
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:499
functor for a local U1 dot U1 potential
Definition correlationfactor.h:532
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:562
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:534
U1_dot_U1_functor(const NuclearCorrelationFactor *ncf)
Definition correlationfactor.h:537
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:539
functor for the local part of the U1 potential – NOTE THE SIGN
Definition correlationfactor.h:464
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:473
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:485
const int axis
Definition correlationfactor.h:467
U1_functor(const NuclearCorrelationFactor *ncf, const int axis)
Definition correlationfactor.h:470
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:466
compute the derivative of U2 wrt the displacement of atom A
Definition correlationfactor.h:844
const int iatom
Definition correlationfactor.h:846
const int axis
Definition correlationfactor.h:847
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:869
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:845
U2X_functor(const NuclearCorrelationFactor *ncf, const int &atom1, const int axis)
Definition correlationfactor.h:849
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:856
U2 functor for a specific atom.
Definition correlationfactor.h:618
const size_t iatom
Definition correlationfactor.h:621
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:627
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:634
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:620
U2_atomic_functor(const NuclearCorrelationFactor *ncf, const size_t atom)
Definition correlationfactor.h:624
Definition correlationfactor.h:568
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:582
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:572
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:569
U2_functor(const NuclearCorrelationFactor *ncf)
Definition correlationfactor.h:571
compute the derivative of U3 wrt the displacement of atom A, coord axis
Definition correlationfactor.h:892
U3X_functor(const NuclearCorrelationFactor *ncf, const size_t iatom, const int axis)
Definition correlationfactor.h:897
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:893
const int axis
Definition correlationfactor.h:895
const size_t iatom
Definition correlationfactor.h:894
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:900
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:942
U3 functor for a specific atom.
Definition correlationfactor.h:645
const size_t iatom
Definition correlationfactor.h:648
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:647
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:654
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:678
U3_atomic_functor(const NuclearCorrelationFactor *ncf, const int atom)
Definition correlationfactor.h:651
Definition correlationfactor.h:587
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:612
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:591
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:588
U3_functor(const NuclearCorrelationFactor *ncf)
Definition correlationfactor.h:590
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:724
const Molecule & molecule
Definition correlationfactor.h:717
const size_t iatom
Definition correlationfactor.h:718
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:737
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:716
square_times_V_derivative_functor(const NuclearCorrelationFactor *ncf, const Molecule &molecule1, const size_t atom1, const int axis1)
Definition correlationfactor.h:721
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition correlationfactor.h:709
const size_t iatom
Definition correlationfactor.h:691
const NuclearCorrelationFactor * ncf
Definition correlationfactor.h:689
square_times_V_functor(const NuclearCorrelationFactor *ncf, const Molecule &mol, const size_t iatom1)
Definition correlationfactor.h:693
const Molecule & molecule
Definition correlationfactor.h:690
double operator()(const coord_3d &xyz) const
Definition correlationfactor.h:696
ABC for the nuclear correlation factors.
Definition correlationfactor.h:83
double eprec
smoothing of the potential/step function
Definition correlationfactor.h:211
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition correlationfactor.h:87
virtual const real_function_3d U2() const
return the U2 term of the correlation function
Definition correlationfactor.h:200
virtual real_function_3d square() const
return the square of the nuclear correlation factor
Definition correlationfactor.h:159
double vtol
the threshold for initial projection
Definition correlationfactor.h:208
virtual double Srrr_div_S(const double &r, const double &Z) const =0
virtual real_function_3d function() const
return the nuclear correlation factor
Definition correlationfactor.h:151
coord_3d smoothed_unitvec(const coord_3d &xyz, double smoothing=0.0) const
smoothed unit vector for the computation of the U1 potential
Definition correlationfactor.h:316
const Molecule & molecule
the molecule
Definition correlationfactor.h:214
virtual real_function_3d apply_U(const real_function_3d &rhs) const
apply the regularized potential U_nuc on a given function rhs
Definition correlationfactor.h:133
virtual real_function_3d inverse() const
return the inverse nuclear correlation factor
Definition correlationfactor.h:178
virtual double Spp_div_S(const double &r, const double &Z) const =0
the regularized potential wrt a given atom wrt the cartesian coordinate
World & world
the world
Definition correlationfactor.h:205
real_function_3d U2_function
the purely local U2 potential, having absorbed the nuclear pot V_nuc
Definition correlationfactor.h:221
virtual ~NuclearCorrelationFactor()
virtual destructor
Definition correlationfactor.h:98
virtual double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt nuclear coordinate X (spherical part)
Definition correlationfactor.h:294
virtual double S(const double &r, const double &Z) const =0
the correlation factor S wrt a given atom
virtual real_function_3d square_times_V_derivative(const int iatom, const int axis) const
Definition correlationfactor.h:170
void initialize(const double vtol1)
initialize the regularized potentials U1 and U2
Definition correlationfactor.h:101
virtual coord_3d Sp(const coord_3d &vr1A, const double &Z) const =0
the partial derivative of correlation factor S' wrt the cartesian coordinates
virtual const real_function_3d U1(const int axis) const
return the U1 term of the correlation function
Definition correlationfactor.h:186
NuclearCorrelationFactor(World &world, const Molecule &mol)
ctor
Definition correlationfactor.h:93
std::vector< real_function_3d > U1_function
the three components of the U1 potential
Definition correlationfactor.h:218
virtual corrfactype type() const =0
std::vector< real_function_3d > U1vec() const
return the U1 functions in a vector
Definition correlationfactor.h:191
coord_3d dsmoothed_unitvec(const coord_3d &xyz, const int axis, double smoothing=0.0) const
derivative of smoothed unit vector wrt the electronic coordinate
Definition correlationfactor.h:377
virtual double Srr_div_S(const double &r, const double &Z) const =0
corrfactype
Definition correlationfactor.h:85
@ Two
Definition correlationfactor.h:86
@ LinearSlater
Definition correlationfactor.h:85
@ GradientalGaussSlater
Definition correlationfactor.h:85
@ poly4erfc
Definition correlationfactor.h:86
@ Slater
Definition correlationfactor.h:86
@ Adhoc
Definition correlationfactor.h:86
@ GaussSlater
Definition correlationfactor.h:85
@ None
Definition correlationfactor.h:85
@ Polynomial
Definition correlationfactor.h:86
virtual double Sr_div_S(const double &r, const double &Z) const =0
A nuclear correlation factor class.
Definition correlationfactor.h:1781
double a_
length scale parameter a, default chosen that linear terms in U2 vanish
Definition correlationfactor.h:1810
Polynomial(World &world, const Molecule &mol, const double a)
ctor
Definition correlationfactor.h:1787
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:1834
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1850
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:1818
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1908
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt nuclear coordinate X (spherical part)
Definition correlationfactor.h:1922
corrfactype type() const
Definition correlationfactor.h:1805
double a_param() const
Definition correlationfactor.h:1812
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1894
static double b_param(const double &a)
the cutoff
Definition correlationfactor.h:1815
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:1878
Definition correlationfactor.h:1954
double Srrr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2011
double eprec
Definition correlationfactor.h:2002
corrfactype type() const
Definition correlationfactor.h:1976
std::shared_ptr< PotentialManager > potentialmanager
underlying potential (=molecule)
Definition correlationfactor.h:2001
real_function_3d apply_U(const real_function_3d &rhs) const
apply the regularized potential U_nuc on a given function rhs
Definition correlationfactor.h:1993
double Srr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2009
const double fac
the factor of the correlation factor: R=fac;
Definition correlationfactor.h:2005
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt nuclear coordinate X (spherical part)
Definition correlationfactor.h:2029
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:2014
double Sr_div_S(const double &r, const double &Z) const
Definition correlationfactor.h:2007
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:2019
PseudoNuclearCorrelationFactor(World &world, const Molecule &mol, const std::shared_ptr< PotentialManager > pot, const double fac)
ctor
Definition correlationfactor.h:1961
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:2024
const real_function_3d U2() const
return the U2 term of the correlation function
Definition correlationfactor.h:1982
A nuclear correlation factor class.
Definition correlationfactor.h:1323
double eprec_param() const
Definition correlationfactor.h:1353
corrfactype type() const
Definition correlationfactor.h:1344
double a_param() const
Definition correlationfactor.h:1352
double Srr_div_S(const double &r, const double &Z) const
second derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1370
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:1399
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1406
Slater(World &world, const Molecule &mol, const double a)
ctor
Definition correlationfactor.h:1329
double eprec_
Definition correlationfactor.h:1350
double Sr_div_S(const double &r, const double &Z) const
first derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1360
double Srrr_div_S(const double &r, const double &Z) const
third derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1381
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt X (scalar part)
Definition correlationfactor.h:1436
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:1388
double a_
the length scale parameter
Definition correlationfactor.h:1349
T normf() const
Calculate the 2-norm of the vector elements.
Definition vector.h:400
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
Definition correlationfactor.h:1467
double Srrr_div_S(const double &r, const double &Z) const
third derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1638
double a
the length scale parameter
Definition correlationfactor.h:1514
double Spp_div_S(const double &r, const double &Z) const
second derivative of the nuclear correlation factor
Definition correlationfactor.h:1658
poly4erfc(World &world, const Molecule &mol, const double aa)
ctor
Definition correlationfactor.h:1473
double a0
Definition correlationfactor.h:1515
double eprec_
Definition correlationfactor.h:1516
coord_3d Sp(const coord_3d &vr1A, const double &Z) const
radial part first derivative of the nuclear correlation factor
Definition correlationfactor.h:1652
corrfactype type() const
Definition correlationfactor.h:1509
double U2X_spherical(const double &r, const double &Z, const double &rcut) const
derivative of the U2 potential wrt X (scalar part)
Definition correlationfactor.h:1767
double a_param() const
Definition correlationfactor.h:1518
double eprec_param() const
Definition correlationfactor.h:1519
double Srr_div_S(const double &r, const double &Z) const
second derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1628
double Sr_div_S(const double &r, const double &Z) const
first derivative of the correlation factor wrt (r-R_A)
Definition correlationfactor.h:1526
double S(const double &r, const double &Z) const
the nuclear correlation factor
Definition correlationfactor.h:1644
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
static double lo
Definition dirac-hatom.cc:23
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only)
Definition tensor.h:2502
static const double rcut
Definition hatom_sf_dirac.cc:19
static const double eprec
Definition hatom_sf_dirac.cc:18
static double Z2(const coord_3d &r)
Definition helium_mp2.cc:124
#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
void print(const tensorT &t)
Definition mcpfit.cc:140
static double dsmoothed_potential(double r)
Derivative of the regularized 1/r potential.
Definition mentity.cc:223
static double smoothed_potential(double r)
Regularized 1/r potential.
Definition mentity.cc:206
static double smoothing_parameter(double Z, double eprec)
Returns radius for smoothing nuclear potential with energy precision eprec.
Definition mentity.cc:192
Main include file for MADNESS and defines Function interface.
const double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
std::shared_ptr< NuclearCorrelationFactor > create_nuclear_correlation_factor(World &world, const Molecule &molecule, const std::shared_ptr< PotentialManager > potentialmanager, const std::string inputline)
create and return a new nuclear correlation factor
Definition correlationfactor.cc:45
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
int power< 4 >(int base)
Definition power.h:68
static double r2(const coord_3d &x)
Definition smooth.h:45
int power< 5 >(int base)
Definition power.h:73
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
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
FunctionFactory< double, 3 > real_factory_3d
Definition functypedefs.h:93
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
Function< TENSOR_RESULT_TYPE(T, R), NDIM > dot(World &world, const std::vector< Function< T, NDIM > > &a, const std::vector< Function< R, NDIM > > &b, bool fence=true)
Multiplies and sums two vectors of functions r = \sum_i a[i] * b[i].
Definition vmra.h:1493
NDIM const Function< R, NDIM > & g
Definition mra.h:2416
int power< 2 >(int base)
Definition power.h:58
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
int power< 3 >(int base)
Definition power.h:63
Vector< double, 3 > coord_3d
Definition funcplot.h:1042
int power< 6 >(int base)
Definition power.h:78
static const double b
Definition nonlinschro.cc:119
static const double a
Definition nonlinschro.cc:118
static const double c
Definition relops.cc:10
static double Z
Definition rk.cc:35
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
Definition dirac-hatom.cc:108
static double V(const coordT &r)
Definition tdse.cc:288
void e()
Definition test_sig.cc:75
double aa
Definition testbsh.cc:68
#define N
Definition testconv.cc:37
std::size_t axis
Definition testpdiff.cc:59
const double a2
Definition vnucso.cc:86
const double R2
Definition vnucso.cc:84
const double a1
Definition vnucso.cc:85