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 
80 namespace madness {
81 
82 /// ABC for the nuclear correlation factors
84 public:
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);
180  real_function_3d R_inverse=real_factory_3d(world).thresh(vtol)
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 
202 private:
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 
216 protected:
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 
223 private:
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 
249 public:
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 
302 public:
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]
377  coord_3d dsmoothed_unitvec(const coord_3d& xyz, const int axis,
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> {
436  int exponent;
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]
497  class U1_atomic_functor : public FunctionFunctorInterface<double,3> {
498 
500  const size_t iatom;
501  const int axis;
502 
503  public:
504  U1_atomic_functor(const NuclearCorrelationFactor* ncf, const size_t atom,
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!
532  class U1_dot_U1_functor : public FunctionFunctorInterface<double,3> {
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
618  class U2_atomic_functor : public FunctionFunctorInterface<double,3> {
619 
621  const size_t iatom;
622 
623  public:
624  U2_atomic_functor(const NuclearCorrelationFactor* ncf, const size_t atom)
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
645  class U3_atomic_functor : public FunctionFunctorInterface<double,3> {
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  }
704  const double V=-molecule.atomic_attraction_potential(
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> {
745  const Atom& thisatom;
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),
752  derivativeaxis(daxis), exponent(exponent) {
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)),
759  derivativeaxis(daxis), exponent(exponent) {
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> {
799  const Atom& thisatom;
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
957 public:
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 
977 private:
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
1089 public:
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)
1095  : NuclearCorrelationFactor(world,mol), a(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 
1110 private:
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
1233 public:
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 
1254 private:
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
1324 public:
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 
1346 private:
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 
1468 public:
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 
1511 private:
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
1780 template<std::size_t N>
1782 public:
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 
1807 private:
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 
1956 public:
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 
1998 private:
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 
2040 public:
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 
2059 private:
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 
2096 std::shared_ptr<NuclearCorrelationFactor>
2098  const Molecule& molecule,
2099  const std::shared_ptr<PotentialManager> pm,
2100  const std::string inputline);
2101 
2102 std::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
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
double operator()(const coord_3d &xyz) const
Definition: correlationfactor.h:763
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: correlationfactor.h:789
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
int exponent
Definition: correlationfactor.h:436
const NuclearCorrelationFactor * ncf
Definition: correlationfactor.h:435
double operator()(const coord_3d &xyz) const
Definition: correlationfactor.h:440
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: correlationfactor.h:456
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
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: correlationfactor.h:832
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
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
double operator()(const coord_3d &xyz) const
Definition: correlationfactor.h:572
const NuclearCorrelationFactor * ncf
Definition: correlationfactor.h:569
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: correlationfactor.h:582
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
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 Molecule & molecule
Definition: correlationfactor.h:717
const size_t iatom
Definition: correlationfactor.h:718
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
const size_t iatom
Definition: correlationfactor.h:691
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 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
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
std::vector< real_function_3d > U1vec() const
return the U1 functions in a vector
Definition: correlationfactor.h:191
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
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 double Z2(const coord_3d &r)
Definition: helium_mp2.cc:124
static double pow(const double *a, const double *b)
Definition: lda.h:74
#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
Main include file for MADNESS and defines Function interface.
const double pi
Mathematical constant .
Definition: constants.h:48
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
double smoothed_potential(double r)
Smoothed 1/r potential.
Definition: atomutil.cc:220
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:1436
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
double dsmoothed_potential(double r)
Derivative of the regularized 1/r potential.
Definition: atomutil.cc:333
double smoothing_parameter(double Z, double eprec)
Returns radius for smoothing nuclear potential with energy precision eprec.
Definition: atomutil.cc:203
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
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
NDIM const Function< R, NDIM > & g
Definition: mra.h:2416
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
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