72 #ifndef MADNESS_CHEM_NUCLEARCORRELATIONFACTOR_H_
73 #define MADNESS_CHEM_NUCLEARCORRELATIONFACTOR_H_
87 typedef std::shared_ptr< FunctionFunctorInterface<double,3> >
functorT;
111 .functor(U1f).truncate_on_project());
118 .functor(U2f).truncate_on_project();
124 .functor(U3f).truncate_on_project();
154 .functor(Rf).truncate_on_project();
162 .functor(r).truncate_on_project();
173 .functor(
func).truncate_on_project();
181 .functor(r).truncate_on_project();
191 std::vector<real_function_3d>
U1vec()
const {
192 std::vector<real_function_3d> uvec(3);
229 virtual double S(
const double& r,
const double&
Z)
const = 0;
247 virtual double Spp_div_S(
const double& r,
const double&
Z)
const = 0;
259 virtual double Sr_div_S(
const double& r,
const double&
Z)
const = 0;
269 virtual double Srr_div_S(
const double& r,
const double&
Z)
const = 0;
279 virtual double Srrr_div_S(
const double& r,
const double&
Z)
const = 0;
296 print(
"you can't compute the Hessian matrix");
297 print(
"U2X_spherical is not implemented for the nuclear correlation factor");
323 smoothing=sqrt(smoothing);
324 const double r=xyz.
normf();
325 const double rs=r/smoothing;
328 double erfrs_div_r=2.0/(smoothing*sqrtpi)-2.0/3.0*rs*rs/(sqrtpi*smoothing);
329 return erfrs_div_r*xyz;
331 return erf(rs)/r*xyz;
339 if (smoothing==0.0) smoothing=
eprec;
343 const double r=xyz.
normf();
344 const double cutoff=smoothing;
348 const double xi=r/cutoff;
349 const double xi2=
xi*
xi;
350 const double xi3=
xi*
xi*
xi;
352 const double nu22=0.5 + 1./64.*(105*
xi - 175 *xi3 + 147* xi2*xi3 - 45* xi3*xi3*
xi);
354 const double kk=2.*nu22-1.0;
378 double smoothing=0.0)
const {
380 const double r=xyz.
normf();
382 if (smoothing==0.0) smoothing=
eprec;
388 smoothing=sqrt(smoothing);
390 const double rs=r/smoothing;
392 const double sqrtpis3=sqrtpi*smoothing*smoothing*smoothing;
396 double p=-4.0/(3.0*sqrtpis3) + 4.0*rs*rs/(5.0*sqrtpis3);
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;
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;
410 double s2=smoothing*smoothing;
411 double s7=s2*s2*s2*smoothing;
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));
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;
427 result=xyz*(-xyz[
axis]/(r*r*r));
442 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
443 const Atom& atom=
ncf->molecule.get_atom(i);
445 const double r=vr1A.
normf();
446 result*=
ncf->S(r,atom.
q);
448 if (
exponent==-1)
return 1.0/result;
449 else if (
exponent==2)
return result*result;
450 else if (
exponent==1)
return result;
457 return ncf->molecule.get_all_coords_vec();
475 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
476 const Atom& atom=
ncf->molecule.get_atom(i);
478 const double r=vr1A.
normf();
479 const double&
Z=atom.
q;
481 result-=
ncf->Sr_div_S(r,
Z)*
ncf->smoothed_unitvec(vr1A)[
axis];
486 return ncf->molecule.get_all_coords_vec();
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];
516 std::vector< madness::Vector<double,3> >
c(1);
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);
545 const double r=vr1A.
normf();
551 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
552 for (
size_t j=0; j<
ncf->molecule.natom(); ++j) {
563 return ncf->molecule.get_all_coords_vec();
574 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
575 const Atom& atom=
ncf->molecule.get_atom(i);
577 const double r=vr1A.
normf();
578 result+=
ncf->Spp_div_S(r,atom.
q);
583 return ncf->molecule.get_all_coords_vec();
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);
596 const double r=vr1A.
normf();
598 all_terms[i]=
ncf->Sr_div_S(r,atom.
q)*
ncf->smoothed_unitvec(vr1A);
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];
613 return ncf->molecule.get_all_coords_vec();
630 const double r=vr1A.
normf();
631 return ncf->Spp_div_S(r,atom.
q);
635 std::vector< madness::Vector<double,3> >
c(1);
657 const double rA=vr1A.
normf();
659 double Sr_div_SA=
ncf->Sr_div_S(rA,atomA.
q);
663 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
664 if (i==
iatom)
continue;
666 const Atom& atomB=
ncf->molecule.get_atom(i);
668 const double rB=vr1B.
normf();
670 double Sr_div_SB=
ncf->Sr_div_S(rB,atomB.
q);
672 double dot=nA[0]*nB[0] + nA[1]*nB[1] + nA[2]*nB[2];
673 result+=Sr_div_SB*Sr_div_SA*
dot;
679 std::vector< madness::Vector<double,3> >
c(1);
694 const Molecule& mol,
const size_t iatom1)
698 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
699 const Atom& atom=
ncf->molecule.get_atom(i);
701 const double r=vr1A.
normf();
702 result*=
ncf->S(r,atom.
q);
705 iatom, xyz[0], xyz[1], xyz[2]);
706 return result*result*
V;
710 return ncf->molecule.get_all_coords_vec();
722 const Molecule& molecule1,
const size_t atom1,
const int axis1)
726 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
727 const Atom& atom=
ncf->molecule.get_atom(i);
729 const double r=vr1A.
normf();
730 result*=
ncf->S(r,atom.
q);
734 return result*result*Vprime;
738 return ncf->molecule.get_all_coords_vec();
768 for (
size_t i=0; i<
ncf->molecule.natom(); ++i) {
769 const Atom& atom=
ncf->molecule.get_atom(i);
771 const double r=vr1A.
normf();
772 result*=
ncf->S(r,atom.
q);
774 if (
exponent==2) result=result*result;
780 const double r=vr1A.
normf();
782 const double S1=-
ncf->Sr_div_S(r,
Z)
790 return ncf->molecule.get_all_coords_vec();
820 const double r=vr1A.
normf();
822 const double S1=
ncf->Sr_div_S(r,
Z);
823 const double S2=
ncf->Srr_div_S(r,
Z);
828 return drhodx*(S2-S1*S1)*
ncf->smoothed_unitvec(vr1A)[
U1axis]
833 std::vector< madness::Vector<double,3> >
c(1);
852 double lo=1.0/atom.
q;
859 const double r=vr1A.
normf();
860 const double&
Z=atom.
q;
865 const double drhodx=-
ncf->smoothed_unitvec(vr1A)[
axis];
866 return drhodx*
ncf->U2X_spherical(r,
Z,
rcut);
870 std::vector< madness::Vector<double,3> >
c(1);
903 const double r1A=vr1A.
normf();
904 const double& ZA=atomA.
q;
906 double S1A=
ncf->Sr_div_S(r1A,ZA);
907 double S2A=
ncf->Srr_div_S(r1A,ZA);
908 double termA=S2A-S1A*S1A;
915 const double drhodx=-nA[
axis];
918 for (
size_t jatom=0; jatom<
ncf->molecule.natom(); ++jatom) {
919 if (
iatom==jatom)
continue;
921 const Atom& atomB=
ncf->molecule.get_atom(jatom);
923 const double r1B=vr1B.
normf();
924 const double& ZB=atomB.
q;
926 double S1B=
ncf->Sr_div_S(r1B,ZB);
931 for (
int i=0; i<3; ++i) {
935 term+=(+drhodx*termA*S1B*
dot + S1A*S1B*ddot);
943 return ncf->molecule.get_all_coords_vec();
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))");
970 print(
"which is of Gaussian-Slater type\n");
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)));
988 const double r=sqrt(vr1A[0]*vr1A[0] +
989 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
991 const double eA=exp(-
Z*r);
992 const double gA=exp(-
Z*
Z*r*r);
1001 const double rho=
Z*r;
1003 return Z*
Z*(-3.5 - 4.0*rho + 6.0*rho*rho + 12.0*rho*rho*rho);
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;
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;
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;
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;
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);
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;
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))");
1103 print(
"which is of Gradiental Gaussian-Slater type\n");
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)));
1123 const double r=sqrt(vr1A[0]*vr1A[0] +
1124 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
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;
1136 const double rho=
Z*r;
1137 const double sqrtz=sqrt(
Z);
1139 const double zfivehalf=
Z*
Z*sqrtz;
1140 const double a2=
a*
a;
1141 const double a4=
a2*
a2;
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;
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);
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;
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);
1177 const double denom=1.0-
g+
e/sqrtz;
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);
1188 const double denom=
e+sqrtz-
g*sqrtz;
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;
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);
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;
1244 print(
"constructed nuclear correlation factor of the form");
1245 print(
" S_A = -Z_A r_{1A} exp(-Z_A r_{1A}) + 1");
1248 print(
"which is of linear Slater type\n");
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;
1271 const double b=a_param();
1272 const double r=sqrt(vr1A[0]*vr1A[0] +
1273 vr1A[1]*vr1A[1] + vr1A[2]*vr1A[2]);
1275 const double ebrz=exp(-
b*r*
Z);
1285 const double b=a_param();
1286 const double rho=
Z*r;
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);
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);
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);
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);
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);
1336 print(
"\nconstructed nuclear correlation factor of the form");
1337 print(
" S_A = 1/(a-1) exp(-a Z_A r_{1A}) + 1");
1339 print(
"with eprec ",eprec_);
1340 print(
"which is of Slater type\n");
1361 const double&
a=a_param();
1362 return -
a*
Z/(1.0+(
a-1.0)*exp(
a*r*
Z));
1371 const double&
a=a_param();
1372 const double aZ=
a*
Z;
1373 return aZ*aZ/(1.0+(
a-1.0)*exp(r*aZ));
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));
1388 double S(
const double& r,
const double&
Z)
const {
1389 const double a=a_param();
1391 return 1.0+1.0/(
a-1.0) * exp(-
a*
Z*r);
1400 const double a=a_param();
1401 const double r=vr1A.
normf();
1407 const double a=a_param();
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);
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);
1437 const double a=a_param();
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);
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;
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)");
1483 print(
"with eprec ",eprec_);
1484 print(
"which is of poly4erfc type\n");
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;
1504 print(
"invalid parameter a for poly4erfc: only 0.5 and 1.0 implemented");
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))));
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)*
1552 (16.532205048243702951212522516 + x*
1553 (-11.89273747187279634240347945 +
1554 x*(15.157537549656745468895369276 + x*(-4.9102960292797655978798640519 + 1.*x))));
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))
1564 (886.4859678528423041797649741 + x*(269.17746130370931387996124706 +
1565 x*(-130.21383548057958115685397713 + x*(37.644499985765056273193347388 + 1.*x))));
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))));
1580 }
else if (
a==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
1590 (1.6046957999975126909719683196 + x*
1591 (-0.8234878506688316215458988304 +
1592 x*(3.4607641859010551501639903314 + x*(-1.8955210085531557309609670978 + 1.*x))));
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
1602 (5.4509691924562038998044993827 + x*
1603 (-2.2732931206867811068721699796 +
1604 x*(4.1530634219989859344450742259 + x*(-2.0183662125874366044951391259 + 1.*x))));
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))));
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;
1648 return 1.0 + (a0 +
a1* arZ+
a2 * arZ2 + a3*arZ*arZ2 + + a4*arZ2*arZ2) *erfc(
a*r*
Z);
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))));
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))));
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))));
1705 }
else if (
a==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))));
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))));
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))));
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))));
1780 template<std::
size_t N>
1791 a_=(2. + (-2. + sqrt(-1. +
N))*
N)/(-2. +
N);
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 ");
1801 print(
"which is of polynomial type with exponent N = ",
N);
1818 double S(
const double& r,
const double&
Z)
const {
1820 const double rho=r*
Z;
1825 const double arg=-1.0 + rho/
b;
1826 return 1.0 + power<N>(-1.0) *
a*power<N>(
arg);
1836 const double r=vr1A.
normf();
1837 const double rho=r*
Z;
1852 const double rho=r*
Z;
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 +
1862 return Z*
Z*(c0 + c1*r + c2*r*r);
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 +
1870 const double denom=2.* (r + power<N>(-1) *
a* r* power<N>(-1 + rho/
b));
1879 const double rho=r*
Z;
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)));
1895 const double rho=r*
Z;
1900 const double negn= power<N>(-1.0);
1902 (
a*
N*(1 + negn*
a*power<N>(-1 + ((1 +
a)*r*
Z)/(
a*
N))));
1909 const double rho=r*
Z;
1914 const double negn= power<N>(-1.0);
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) {
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))/
1935 const double r2=((-4*(17 + 9*rn) +
N*(92 + 80*rn +
1938 result=(r0 + r*r1 + r*r*
r2);
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;
1962 const std::shared_ptr<PotentialManager> pot,
const double fac)
1964 eprec(mol.get_eprec()), fac(fac) {
1967 print(
"constructed nuclear correlation factor of the form");
1970 print(
"which means it's (nearly) a conventional calculation\n");
1987 return potentialmanager->vnuclear();
2007 double Sr_div_S(
const double& r,
const double&
Z)
const {
return 0.0;}
2009 double Srr_div_S(
const double& r,
const double&
Z)
const {
return 0.0;}
2011 double Srrr_div_S(
const double& r,
const double&
Z)
const {
return 0.0;}
2014 double S(
const double& r,
const double&
Z)
const {
2046 const std::vector<real_function_3d>&
U1)
2053 print(
"constructed ad hoc nuclear correlation factor");
2077 double S(
const double& r,
const double&
Z)
const {
2096 std::shared_ptr<NuclearCorrelationFactor>
2099 const std::shared_ptr<PotentialManager> pm,
2100 const std::string inputline);
2102 std::shared_ptr<NuclearCorrelationFactor>
2105 const std::shared_ptr<PotentialManager> pm,
2106 const std::pair<std::string,double>&
ncf);
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
Definition: correlationfactor.h:715
const int axis
Definition: correlationfactor.h:719
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
Definition: correlationfactor.h:688
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