87 return 1.0/(1.0 + exp(-x));
97 template <std::
size_t NDIM>
110 for (std::size_t i=0; i<
NDIM; ++i) {
134 r[0]=x; r[1]=y; r[2]=z;
144 for (
int i=-
R; i<=+
R; i++) {
145 double xx = x[0]+i*
L;
146 for (
int j=-
R; j<=+
R; j++) {
147 double yy = x[1]+j*
L;
148 for (
int k=-
R;
k<=+
R;
k++) {
149 double zz = x[2]+
k*
L;
191 template <
typename T,
int NDIM>
232 const std::vector<KPoint>& kpoints) :
_world(world),
_kpoints(kpoints),
236 for (
unsigned int kp = 0; kp <
_kpoints.size(); kp++)
255 t_rm.insert(t_rm.end(), br.begin(), br.end());
257 std::vector<double> rnvec = norm2<valueT,NDIM>(
_world, t_rm);
261 for (
unsigned int i = 0; i < rnvec.size(); i++) rnorm += rnvec[i];
267 for (
unsigned int kp = 0; kp <
_kpoints.size(); kp++)
271 vecfuncT k_phisa(awfs_old.begin() + kpoint.begin, awfs_old.begin() + kpoint.end);
272 vecfuncT k_phisb(bwfs_old.begin() + kpoint.begin, bwfs_old.begin() + kpoint.end);
273 vecfuncT k_awfs(awfs_new.begin() + kpoint.begin, awfs_new.begin() + kpoint.end);
274 vecfuncT k_bwfs(bwfs_new.begin() + kpoint.begin, bwfs_new.begin() + kpoint.end);
283 k_rm.insert(k_rm.end(), k_br.begin(), k_br.end());
284 k_vm.insert(k_vm.end(), k_phisb.begin(), k_phisb.end());
294 int m = k_subspace.size();
297 for (
int s = 0; s <
m; s++)
299 const vecfuncT& k_vs = k_subspace[s].first;
300 const vecfuncT& k_rs = k_subspace[s].second;
301 for (
unsigned int i = 0; i < k_vm.size(); i++)
303 ms[s] += k_vm[i].inner_local(k_rs[i]);
304 sm[s] += k_vs[i].inner_local(k_rm[i]);
321 double rcond = 1
e-12;
324 if (
abs(
c[
m-1]) < 3.0) {
327 else if (rcond < 0.01) {
329 print(
"Increasing subspace singular value threshold ",
c[
m-1], rcond);
334 print(
"Forcing full step due to subspace malfunction");
346 print(
"Subspace solution",
c);
350 vecfuncT k_phisa_new = zero_functions_compressed<valueT,NDIM>(
_world, k_phisa.
size());
351 vecfuncT k_phisb_new = zero_functions_compressed<valueT,NDIM>(
_world, k_phisb.
size());
352 std::complex<double> one = std::complex<double>(1.0,0.0);
353 unsigned int norbitals = awfs_old.size() /
_kpoints.size();
354 for (
unsigned int m = 0;
m < k_subspace.size();
m++)
356 const vecfuncT& k_vm = k_subspace[
m].first;
357 const vecfuncT& k_rm = k_subspace[
m].second;
359 const vecfuncT vma(k_vm.begin(),k_vm.begin() + norbitals);
360 const vecfuncT rma(k_rm.begin(),k_rm.begin() + norbitals);
361 const vecfuncT vmb(k_vm.end() - norbitals, k_vm.end());
362 const vecfuncT rmb(k_rm.end() - norbitals, k_rm.end());
377 k_subspace.erase(k_subspace.begin());
383 for (
unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
386 awfs_new[wi] = k_phisa_new[fi];
387 bwfs_new[wi] = k_phisb_new[fi];
407 template <
typename T,
int NDIM>
458 vm.insert(vm.end(), bwfs_old.begin(), bwfs_old.end());
470 for (
int s=0; s<
m; s++)
474 for (
unsigned int i=0; i<vm.size(); i++)
476 ms[s] += vm[i].inner_local(rs[i]);
477 sm[s] += vs[i].inner_local(rm[i]);
494 double rcond = 1
e-12;
497 if (
abs(
c[
m-1]) < 3.0) {
500 else if (rcond < 0.01) {
502 print(
"Increasing subspace singular value threshold ",
c[
m-1], rcond);
507 print(
"Forcing full step due to subspace malfunction");
519 print(
"Subspace solution",
c);
525 std::complex<double> one = std::complex<double>(1.0,0.0);
529 const vecfuncT vma(vm.begin(),vm.begin()+awfs_old.size());
530 const vecfuncT rma(rm.begin(),rm.begin()+awfs_old.size());
531 const vecfuncT vmb(vm.end()-bwfs_old.size(), vm.end());
532 const vecfuncT rmb(rm.end()-bwfs_old.size(), rm.end());
550 awfs_new = phisa_new;
551 bwfs_new = phisb_new;
617 template <
typename T,
int NDIM>
827 _matF.open(
"matrix.txt");
828 _eigF.open(
"eigs.txt");
845 coordT c2 {0.5*t1,0.5*t1,0.5*t1};
865 _kmeshF <<
"ik" << setw(10) <<
"kpt" << setw(30) <<
"weight" << endl;
866 _kmeshF <<
"--" << setw(10) <<
"---" << setw(30) <<
"------" << endl;
871 for (
unsigned int i = 0; i <
_kpoints.size(); i++)
874 _kmeshF << i << setw(10) << kpoint.k[0];
875 _kmeshF << setw(10) << kpoint.k[1];
876 _kmeshF << setw(10) << kpoint.k[2];
890 std::vector<KPoint>
genkmesh(
unsigned int ngridk0,
unsigned ngridk1,
unsigned int ngridk2,
891 double koffset0,
double koffset1,
double koffset2,
double R)
893 std::vector<KPoint> kmesh;
894 double step0 = 1.0/ngridk0;
895 double step1 = 1.0/ngridk1;
896 double step2 = 1.0/ngridk2;
897 double weight = 1.0/(ngridk0*ngridk1*ngridk2);
899 for (
unsigned int i = 0; i < ngridk0; i++)
901 for (
unsigned int j = 0; j < ngridk1; j++)
903 for (
unsigned int k = 0;
k < ngridk2;
k++)
908 double k0 = ((i*step0)+koffset0) * TWO_PI/
R;
909 double k1 = ((j*step1)+koffset1) * TWO_PI/
R;
910 double k2 = ((
k*step2)+koffset2) * TWO_PI/
R;
912 kmesh.push_back(kpoint);
925 ar & (
unsigned int)(
_kpoints.size());
929 ar & (
unsigned int)(
_phisa.size());
930 for (
unsigned int i = 0; i <
_phisa.size(); i++) ar &
_phisa[i];
935 ar & (
unsigned int)(
_phisb.size());
936 for (
unsigned int i = 0; i <
_phisb.size(); i++) ar &
_phisb[i];
956 for (
unsigned int i = 0; i < nkpts; i++)
977 for (
unsigned int i = 0; i < norbs; i++)
988 for (
unsigned int i = 0; i <
_phisa.size(); i++)
993 for (
unsigned int i = 0; i <
_kpoints.size(); i++)
1000 for (
unsigned int i = 0; i < norbs; i++)
1011 for (
unsigned int i = 0; i <
_phisb.size(); i++)
1016 for (
unsigned int i = 0; i <
_kpoints.size(); i++)
1021 for (
unsigned int i = 0; i < norbs; i++)
1028 _occsa = std::vector<double>(norbs, 0.0);
1029 _occsb = std::vector<double>(norbs, 0.0);
1044 std::vector<coordT> specialpts;
1049 pt[0] = atom.
x; pt[1] = atom.
y; pt[2] = atom.
z;
1050 specialpts.push_back(pt);
1071 if (
_world.
rank() == 0)
print(
"Done creating nuclear potential ..\n");
1079 double t1 = vnucdiff.
trace();
1080 if (
_world.
rank() == 0) printf(
"Difference in nuclear potential: %15.8f\n\n", t1);
1081 for (
int i = 0; i < 101; i++)
1086 double val = vnucdiff(cpt);
1087 if (
_world.
rank() == 0) printf(
"%10.5f %15.8f\n",pt2,val);
1102 printf(
"Cell parameters\n");
1103 printf(
"------------------------------\n");
1104 print(
"cell(x) is ",csize(0,0), csize(0,1));
1105 print(
"cell(y) is ",csize(1,0), csize(1,1));
1106 print(
"cell(z) is ",csize(2,0), csize(2,1));
1132 for (
int xr = -2; xr <= 2; xr += 1)
1134 for (
int yr = -2; yr <= 2; yr += 1)
1136 for (
int zr = -2; zr <= 2; zr += 1)
1139 x[0]+xr*
R, x[1]+yr*
R, x[2]+zr*
R);
1158 template <
typename Q>
1169 printf(
"periodicity along x-axis: \n");
1170 printf(
"-------------------------------------------------------------\n\n");
1171 for (
int i = 0; i < npts; i++)
1173 printf(
"\n-------------------- edge --------------------\n");
1174 for (
int j = 0; j < npts; j++)
1177 r1[0] = begin+eps; r1[1] = (i*
delta)+begin; r1[2] = (j*
delta)+begin;
1179 dr1[0] = begin+2*eps; dr1[1] = (i*
delta)+begin; dr1[2] = (j*
delta)+begin;
1180 dr2[0] = end-2*eps; dr2[1] = (i*
delta)+begin; dr2[2] = (j*
delta)+begin;
1182 std::string success = val < tol ?
"PASS!" :
"FAIL!";
1183 printf(
"%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1184 r1[0],r1[1],r1[2],
r2[0],
r2[1],
r2[2],val,success.c_str());
1188 printf(
"\nperiodicity along y-axis: \n");
1189 printf(
"-------------------------------------------------------------\n\n");
1190 for (
int i = 0; i < npts; i++)
1192 printf(
"\n-------------------- edge --------------------\n");
1193 for (
int j = 0; j < npts; j++)
1196 r1[0] = (i*
delta)+begin; r1[1] = begin+eps; r1[2] = (j*
delta)+begin;
1198 dr1[0] = (i*
delta)+begin; dr1[1] = begin+2*eps; dr1[2] = (j*
delta)+begin;
1199 dr2[0] = (i*
delta)+begin; dr2[1] = end-2*eps; dr2[2] = (j*
delta)+begin;
1201 std::string success = val < tol ?
"PASS!" :
"FAIL!";
1202 printf(
"%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1203 r1[0],r1[1],r1[2],
r2[0],
r2[1],
r2[2],val,success.c_str());
1207 printf(
"\nperiodicity along z-axis: \n");
1208 printf(
"-------------------------------------------------------------\n\n");
1209 for (
int i = 0; i < npts; i++)
1211 printf(
"\n-------------------- edge --------------------\n");
1212 for (
int j = 0; j < npts; j++)
1215 r1[0] = (i*
delta)+begin; r1[1] = (j*
delta)+begin; r1[2] = begin+eps;
1217 dr1[0] = (i*
delta)+begin; dr1[1] = (j*
delta)+begin; dr1[2] = begin+2*eps;
1218 dr2[0] = (i*
delta)+begin; dr2[1] = (j*
delta)+begin; dr2[2] = end-2*eps;
1220 std::string success = val < tol ?
"PASS!" :
"FAIL!";
1221 printf(
"%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1222 r1[0],r1[1],r1[2],
r2[0],
r2[1],
r2[2],val,success.c_str());
1239 vlda.
unaryop(&::ldaop<double>);
1263 _params.
L, kpt.k[0], kpt.k[1], kpt.k[2]));
1264 ao[i] =
factoryT(world).functor(aofunc).truncate_on_project().truncate_mode(0);
1293 std::vector<T> eigsa,
1294 std::vector<T> eigsb,
1370 if (
_world.
rank() == 0)
print(
"Building effective potential ...\n\n");
1372 vlocal =
_vnuc + vc;
1379 vlocal = vlocal + vlda;
1404 _eigsa = std::vector<double>(norbs, 0.0);
1405 _eigsb = std::vector<double>(norbs, 0.0);
1406 _occsa = std::vector<double>(norbs, 0.0);
1407 _occsb = std::vector<double>(norbs, 0.0);
1408 if (
_world.
rank() == 0)
print(
"Building kinetic energy matrix ...\n\n");
1410 for (
int ki = 0; ki < nkpts; ki++)
1422 for (
unsigned int ai = 0; ai < ao.size(); ai++)
1424 std::vector<long> npt(3,101);
1426 sprintf(fnamedx,
"aofunc_k_%2.2d__%2.2d__.dx",ki,ai);
1427 std::vector<long> npt2(3,101);
1433 for (
int i = 0; i < ao.size(); i++)
1436 sprintf(fname,
"ao_line_%d.txt",i);
1440 sprintf(fname2,
"vnuc_line.txt");
1443 sprintf(fname3,
"vlocal_line.txt");
1483 if (
_world.
rank() == 0)
print(
"Building potential energy matrix ...\n\n");
1486 for (
int i = 0; i < ao.size(); i++)
1487 vpsi.push_back(vlocal*ao[i]);
1506 print(
"Potential:");
1508 print(
"Fock: (pre-symmetrized)");
1510 print(
"FockZero: (should be zero)");
1516 for (
unsigned int i = 0; i < fock.
dim(0); i++)
1522 sygv(fock, overlap, 1,
c,
e);
1526 sygv(kinetic, overlap, 1, ck, ek);
1532 syev(overlap, co, eo);
1536 print(
"fock eigenvectors dims:",
c.dim(0),
c.dim(1));
1537 print(
"fock eigenvectors:");
1543 print(
"kinetic eigenvalues");
1549 print(
"potential eigenvalues");
1555 print(
"overlap eigenvalues");
1560 print(
"fock eigenvalues");
1566 printf(
"Overlap: \n");
1567 for (
int i = 0; i < overlap.
dim(0); i++)
1569 for (
int j = 0; j < overlap.
dim(1); j++)
1571 printf(
"%10.5f",
real(overlap(i,j)));
1721 if (
_world.
rank() == 0) printf(
"(%8.4f,%8.4f,%8.4f)\n",kpoint.k[0], kpoint.k[1], kpoint.k[2]);
1731 printf(
"Overlap: \n");
1732 for (
int i = 0; i < kinetic.
dim(0); i++)
1734 for (
int j = 0; j < kinetic.
dim(1); j++)
1736 printf(
"%10.5f",
real(overlap(i,j)));
1743 printf(
"Kinetic: \n");
1744 for (
int i = 0; i < kinetic.
dim(0); i++)
1746 for (
int j = 0; j < kinetic.
dim(1); j++)
1748 printf(
"%10.5f",
real(kinetic(i,j)));
1756 for (
int i = 0; i <
potential.dim(0); i++)
1758 for (
int j = 0; j <
potential.dim(1); j++)
1768 for (
int i = 0; i < fock.
dim(0); i++)
1770 for (
int j = 0; j < fock.
dim(1); j++)
1772 printf(
"%10.5f",
real(fock(i,j)));
1779 printf(
"New overlap: \n");
1780 for (
int i = 0; i < overlap2.
dim(0); i++)
1782 for (
int j = 0; j < overlap2.
dim(1); j++)
1784 printf(
"%10.5f",
real(overlap2(i,j)));
1795 for (
unsigned int oi = kpoint.begin, ti = 0;
oi < kpoint.end;
oi++, ti++)
1800 _phisa.push_back(tmp_orbitals[ti]);
1801 _phisb.push_back(tmp_orbitals[ti]);
1815 const std::vector<T>& eigs,
1850 std::vector<T> eigs,
1851 std::vector<KPoint> kpoints,
1852 std::vector<double> occs,
1894 const std::vector<double>& eigsa,
1895 const std::vector<double>& eigsb,
1896 std::vector<double>& occsa,
1897 std::vector<double>& occsb)
1899 for (
unsigned int ik = 0; ik < kpoints.size(); ik++)
1905 const std::vector<double> k_eigsa(eigsa.begin() +
k.begin, eigsa.begin() +
k.end);
1906 const std::vector<double> k_eigsb(eigsb.begin() +
k.begin, eigsb.begin() +
k.end);
1907 std::vector<double> k_occsa(occsa.begin() +
k.begin, occsa.begin() +
k.end);
1908 std::vector<double> k_occsb(occsb.begin() +
k.begin, occsb.begin() +
k.end);
1911 unsigned int sz = k_eigsa.size();
1916 std::vector<double> teigs;
1920 for (
unsigned int ist = 0; ist < k_eigsa.size(); ist++) teigs.push_back(k_eigsa[ist]);
1921 for (
unsigned int ist = 0; ist < k_eigsb.size(); ist++) teigs.push_back(k_eigsb[ist]);
1923 if (
_world.
rank() == 0) printf(
"setting occs ....\n\n");
1924 for (
unsigned int ist = 0; ist < teigs.size(); ist++)
1928 printf(
"%5d %15.8f\n", ist, teigs[ist]);
1933 std::vector<unsigned int> inds(2*sz);
1934 for (
unsigned int i = 0; i < 2*sz; i++) inds[i] = i;
1936 for (
unsigned int i = 0; i < 2*sz; i++)
1938 for (
unsigned int j = i+1; j < 2*sz; j++)
1940 if (teigs[j] < teigs[i])
1942 double t1 = teigs[i];
1943 teigs[i] = teigs[j];
1953 printf(
"\nSorted eigenvalues:\n");
1954 for (
unsigned int i = 0; i < teigs.size(); i++)
1957 printf(
"%10d%10d %15.8f\n",i,inds[i],teigs[i]);
1962 for (
unsigned int i = 0; (i < 2*sz) && (availablecharge > 0.0) ; i++)
1964 unsigned int current = inds[i];
1968 k_occsb[current] = 1.0;
1969 availablecharge -= 1.0;
1973 k_occsa[current] = 1.0;
1974 availablecharge -= 1.0;
1978 for (
unsigned int ik1 =
k.begin, ik2 = 0; ik1 <
k.end; ik1++,ik2++)
1980 occsa[ik1] = k_occsa[ik2];
1981 occsb[ik1] = k_occsb[ik2];
1985 for (
unsigned int ik = 0; ik < kpoints.size(); ik++)
1990 printf(
"k-point is: %d: \n",ik);
1992 for (
unsigned int ist =
k.begin; ist <
k.end; ist++)
1996 printf(
"occa: %12.5f occb: %12.5f \n",occsa[ist],occsb[ist]);
2011 std::vector<double> occs)
2018 for (
unsigned int kp = 0; kp < kpoints.size(); kp++)
2021 KPoint kpoint = kpoints[kp];
2023 for (
unsigned int j = kpoint.begin; j < kpoint.end; j++)
2030 double rnrm = prod.
trace();
2033 rho += occs[j] * kpoint.
weight * prod;
2043 std::vector<double> occs)
2049 std::vector<rfunctionT> phisq(phis.size());
2050 for (
unsigned int i=0; i<phis.size(); i++) {
2051 phisq[i] =
abssq(phis[i],
false);
2060 for (
unsigned int kp = 0; kp < kpoints.size(); kp++)
2063 KPoint kpoint = kpoints[kp];
2065 for (
unsigned int j = kpoint.begin; j < kpoint.end; j++)
2067 rho.
gaxpy(1.0, phisq[j], occs[j] * kpoint.
weight / (phinorm[j]*phinorm[j]),
false);
2084 std::vector<poperatorT> bops;
2089 int sz = eigs.size();
2090 for (
int i = 0; i < sz; i++)
2097 std::cout <<
"bsh: warning: positive eigenvalue" << i << eps << endl;
2143 for (
unsigned int i = 0; i <
_phisa.size(); i++)
2153 for (
unsigned int i = 0; i <
_phisb.size(); i++)
2188 double rhotrace = rho.
trace();
2189 double rho2trace = rho2.
trace();
2190 if (
_world.
rank() == 0) printf(
"_vnucrhon trace: %10e\n", vnucrhontrace);
2191 if (
_world.
rank() == 0) printf(
"rho trace: %10e\n", rhotrace);
2192 if (
_world.
rank() == 0) printf(
"rho2 trace: %10e\n", rho2trace);
2195 double vlerr = (vlocal-vlocal2).
norm2();
2196 if (
_world.
rank() == 0) printf(
"vlerr trace: %10e\n\n", vlerr);
2200 double ce = 0.5*
inner(vc,rho);
2227 vxc.
unaryop(&::ldaop<double>);
2248 fc.
unaryop(&::ldaeop<double>);
2269 std::cout.precision(8);
2272 printf(
"Energies:\n");
2273 printf(
"Kinetic energy: %20.10f\n", ke);
2274 printf(
"Potential energy: %20.10f\n", pe);
2275 printf(
"Coulomb energy: %20.10f\n", ce);
2276 printf(
"Exchange energy: %20.10f\n",
xc);
2277 printf(
"Total energy: %20.10f\n\n", ke + pe + ce +
xc);
2322 for (
unsigned int j = 0; j < phisa.size(); j++)
2330 funcsa[j] -= tf_jjj;
2332 for (
unsigned int i = j+1; i < phisa.size(); i++)
2340 funcsa[i] -= tf_jij;
2345 funcsa[j] -= tf_iji;
2355 for (
unsigned int i = 0; i <
_kpoints.size(); i++)
2358 if (
k1.is_orb_in_kpoint(idx))
return k1;
2370 for (
unsigned int i = 0; i < phisa.size(); i++)
2374 for (
unsigned int j = 0; j < phisa.size(); j++)
2409 for (
unsigned int ink1 = 0, ik1 = 0; ink1 <
_phisa.size(); ++ink1)
2411 for (
unsigned int ink2 = 0, ik2 = 0; ink2 <= ink1; ink2++)
2416 if (ink1 ==
k1.end) ik1++;
2417 if (ink2 ==
k2.end) ik2++;
2439 funcsa[ink1] -= funcsa[ink2]*fr;
2440 funcsa[ink2] -= funcsa[ink1]*
conj(fr);
2447 printf(
"norm diff: %20.8f\n",
abs(ff.
norm2()-
f.norm2()));
2483 for(
unsigned int i = 0; i <
_phisa.size(); i++)
2494 for(
unsigned int i = 0; i <
_phisb.size(); i++)
2550 double drhoa = (
_rhoa-rhoa).trace();
2551 double drhob = (
_rhob-rhob).trace();
2554 printf(
"diff of alpha rho: %15.7f\n",drhoa);
2555 printf(
"diff of beta rho: %15.7f\n",drhob);
2572 std::vector<functionT> pfuncsa =
2574 std::vector<functionT> pfuncsb =
2584 std::vector<double>
alpha(pfuncsa.size(), 0.0);
2589 std::vector<long> npt(3,101);
2590 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2594 for (
unsigned int kst = kpoint.begin; kst < kpoint.end; kst++, ist++)
2596 std::ostringstream strm;
2597 strm <<
"pre_unk_" << ik <<
"_" << ist <<
".dx";
2598 std::string fname = strm.str();
2606 printf(
"\n\n\n\n------ Debugging BSH operator -----");
2607 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2610 std::vector<double> k_alpha(
alpha.begin() + kpoint.begin,
alpha.begin() + kpoint.end);
2612 printf(
"alpha: (%6.4f,%6.4f,%6.4f)\n",kpoint.k[0],kpoint.k[1],kpoint.k[2]);
2613 for (
unsigned int ia = 0; ia < k_alpha.size(); ia++)
2615 if (
_world.
rank() == 0) printf(
"%15.8f\n", k_alpha[ia]);
2621 if (
_world.
rank() == 0) printf(
"before BSH application ...\n\n");
2622 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2625 double k0 = kpoint.k[0];
2626 double k1 = kpoint.k[1];
2627 double k2 = kpoint.k[2];
2629 printf(
"(%6.5f, %6.5f, %6.5f)\n",
k0,
k1,
k2);
2630 std::vector<functionT> k_phisa(
_phisa.begin() + kpoint.begin,
_phisa.begin() + kpoint.end);
2631 std::vector<functionT> k_pfuncsa(pfuncsa.begin() + kpoint.begin, pfuncsa.begin() + kpoint.end);
2638 syev(overlap, co, eo);
2640 if (
_world.
rank() == 0) printf(
"Overlap eigenvalues: \n");
2642 for (
unsigned int ie = 0; ie < eo.
dim(0); ie++)
2644 if (
_world.
rank() == 0) printf(
"%d %15.8e\n", ie, eo(ie,ie));
2662 std::vector<T> sfactor(pfuncsa.size(), -2.0);
2666 if (
_world.
rank() == 0) std::cout <<
"applying BSH operator ...\n" << endl;
2667 truncate<valueT,NDIM>(
_world, pfuncsa);
2669 std::vector<functionT> tmpa =
apply(
_world, bopsa, pfuncsa);
2675 if (
_world.
rank() == 0) printf(
"norms of tmpa\n");
2677 for (
unsigned int i = 0; i < tmpa_norms.size(); i++)
2679 if (
_world.
rank() == 0) printf(
"%10d %15.8f\n", i, tmpa_norms[i]);
2691 truncate<valueT,NDIM>(
_world, pfuncsb);
2701 std::vector<functionT> pfuncsa2=
2703 std::vector<functionT> pfuncsb2=
2709 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2712 double k0 = kpoint.k[0];
2713 double k1 = kpoint.k[1];
2714 double k2 = kpoint.k[2];
2716 printf(
"(%6.5f, %6.5f, %6.5f)\n",
k0,
k1,
k2);
2717 std::vector<functionT> k_tmpa(tmpa.begin() + kpoint.begin, tmpa.begin() + kpoint.end);
2718 std::vector<functionT> k_pfuncsa2(pfuncsa2.begin() + kpoint.begin, pfuncsa2.begin() + kpoint.end);
2746 pfuncsa2=zero_functions<valueT,NDIM>(
_world,
_phisa.size());
2747 pfuncsb2=zero_functions<valueT,NDIM>(
_world,
_phisa.size());
2749 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2752 double k0 = kpoint.k[0];
2753 double k1 = kpoint.k[1];
2754 double k2 = kpoint.k[2];
2756 printf(
"(%6.5f, %6.5f, %6.5f)\n",
k0,
k1,
k2);
2757 std::vector<functionT> k_phisa(
_phisa.begin() + kpoint.begin,
_phisa.begin() + kpoint.end);
2758 std::vector<functionT> k_pfuncsa2(pfuncsa2.begin() + kpoint.begin, pfuncsa2.begin() + kpoint.end);
2767 std::vector<long> npt(3,101);
2768 for (
unsigned int ik = 0; ik <
_kpoints.size(); ik++)
2772 for (
unsigned int kst = kpoint.begin; kst < kpoint.end; kst++, ist++)
2774 std::ostringstream strm;
2775 strm <<
"unk_" << ik <<
"_" << ist <<
".dx";
2776 std::string fname = strm.str();
2787 const double tol = 1
e-13;
2791 double anorm =
A.normf();
2794 while (anorm*
scale > 0.1)
2803 for (
int i = 0; i < expB.
dim(0); i++) expB(i,i) = std::complex<T>(1.0,0.0);
2807 while (term.
normf() > tol)
2818 expB =
inner(expB,expB);
2826 template<
typename Q>
2830 for (
int i = 0; i < t.
dim(0); i++)
2832 for (
int j = 0; j < t.
dim(1); j++)
2834 os << t(i,j) << setw(12);
2859 print(
"potential eigenvectors dims:",cp.
dim(0),cp.
dim(1));
2860 print(
"potential eigenvectors:");
2863 print(
"potential eigenvalues:");
2882 if (
_world.
rank() == 0)
_outputF <<
"Building kinetic energy matrix ...\n\n" << endl;
2886 tensorT kinetic = ::kinetic_energy_matrix<T,NDIM>(
_world, wf,
2904 sygv(kinetic, overlap, 1, ck, ek);
2910 syev(overlap, co, eo);
2912 sygv(fock, overlap, 1,
c,
e);
2916 print(
"kinetic matrix:");
2918 print(
"\nkinetic eigenvalues:");
2922 print(
"potential matrix:");
2924 print(
"\npotential eigenvalues:");
2928 print(
"fock matrix:");
2930 print(
"\nfock eigenvalues:");
2941 std::vector<KPoint> kpoints,
2942 std::vector<T>&
alpha,
2943 std::vector<double>& eigs)
2946 double trantol = 0.1*
_params.
thresh/std::min(30.0,
double(wf.size()));
2950 for (
unsigned int kp = 0; kp < kpoints.size(); kp++)
2953 KPoint kpoint = kpoints[kp];
2954 double k0 = kpoint.k[0];
2955 double k1 = kpoint.k[1];
2956 double k2 = kpoint.k[2];
2959 vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
2960 vecfuncT k_vwf(vwf.begin() + kpoint.begin, vwf.begin() + kpoint.end);
2995 sygv(fock, overlap, 1, U,
e);
2997 unsigned int nmo = k_wf.size();
3000 for (
long j = 0; j < nmo; j++)
3004 T ang =
arg(U(imax,j));
3005 std::complex<T>
phase = exp(std::complex<T>(0.0,-ang));
3007 for (
long i = 0; i < nmo; i++)
3018 for (
int pass = 0; pass < maxpass; pass++)
3021 for (
long i = 0; i < nmo; i++)
3032 e[i] = tj;
e[j] = ti;
3040 while (ilo < nmo-1) {
3044 if (ihi == nmo-1)
break;
3046 long nclus = ihi - ilo + 1;
3079 printf(
"(%10.5f, %10.5f, %10.5f)\n",
k0,
k1,
k2);
3080 print(
"Overlap matrix:");
3083 print(
"Fock matrix:");
3086 print(
"U matrix: (eigenvectors)");
3089 print(
"Fock matrix eigenvalues:");
3114 for (
unsigned int i = 0; i < k_wf.size(); i++)
3120 d_wf[i] = std::complex<T>(0.0,
k0)*dx_wf +
3121 std::complex<T>(0.0,
k1)*dy_wf +
3122 std::complex<T>(0.0,
k2)*dz_wf;
3125 k_vwf[i] += 0.5 * ksq * k_wf[i];
3126 k_vwf[i] -= d_wf[i];
3131 unsigned int eimax = -1;
3132 double eimaxval = -1e10;
3133 for (
unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3136 if ((
real(
e(fi,fi)) > 0.0) && (
real(
e(fi,fi)) > eimaxval))
3139 eimaxval =
real(
e(fi,fi));
3143 double eshift = (eimaxval > 0.0) ? eimaxval + 0.1 : 0.0;
3144 for (
unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3148 eigs[ei] =
real(
e(fi,fi));
3150 k_vwf[fi] += (
alpha[ei]-eigs[ei])*k_wf[fi];
3155 _eigF <<
"kpt: " << kp << endl;
3156 _eigF << setfill(
'-') << setw(20) <<
" " << endl;
3157 for (
unsigned int ei = kpoint.begin; ei < kpoint.end; ei++)
3160 printf(
"%3d%15.10f",ei,
real(eigs[ei]));
3162 _eigF << eigstr << endl;
3164 _eigF <<
"\n\n" << endl;
3172 sygv(fock, overlap, 1,
c,
e);
3173 for (
unsigned int ei = 0; ei <
e.dim(0); ei++)
3175 double diffe = (ei == 0) ? 0.0 :
real(
e(ei,ei))-
real(
e(ei-1,ei-1));
3177 print(
"kpoint ", kp,
"ei ", ei,
"eps ",
real(
e(ei,ei)),
"\tdiff\t", diffe);
3180 for (
unsigned int ei = kpoint.begin, fi = 0;
3181 ei < kpoint.end; ei++, fi++)
3183 alpha[ei] = std::min(-0.1,
real(fock(fi,fi)));
3184 fock(fi,fi) -= std::complex<T>(
alpha[ei], 0.0);
3191 for (
unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3195 vwf[wi] = k_vwf[fi];
3204 std::vector<KPoint> kpoints,
3205 std::vector<T>&
alpha,
3206 std::vector<double>& eigs)
3209 double trantol = 0.1*
_params.
thresh/std::min(30.0,
double(wf.size()));
3213 for (
unsigned int kp = 0; kp < kpoints.size(); kp++)
3216 KPoint kpoint = kpoints[kp];
3217 double k0 = kpoint.k[0];
3218 double k1 = kpoint.k[1];
3219 double k2 = kpoint.k[2];
3222 vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
3223 vecfuncT k_vwf(vwf.begin() + kpoint.begin, vwf.begin() + kpoint.end);
3228 for (
unsigned int i = 0; i < k_wf.size(); i++)
3229 fock(i,i) += std::complex<double>(i*
_params.
thresh*0.1,0.0);
3232 sygv(fock, overlap, 1, U,
e);
3236 printf(
"(%10.5f, %10.5f, %10.5f)\n",
k0,
k1,
k2);
3237 print(
"Overlap matrix:");
3240 print(
"Fock matrix:");
3243 print(
"U matrix: (eigenvectors)");
3246 print(
"Fock matrix eigenvalues:");
3254 unsigned int nmo = k_wf.size();
3257 for (
long j = 0; j < nmo; j++)
3261 T ang =
arg(U(imax,j));
3262 std::complex<T>
phase = exp(std::complex<T>(0.0,-ang));
3264 for (
long i = 0; i < nmo; i++)
3275 for (
int pass = 0; pass < maxpass; pass++)
3278 for (
long i = 0; i < nmo; i++)
3289 e[i] = tj;
e[j] = ti;
3297 while (ilo < nmo-1) {
3301 if (ihi == nmo-1)
break;
3303 long nclus = ihi - ilo + 1;
3335 k_vwf =
transform(
_world, k_vwf, U, 1
e-5 / std::min(30.0,
double(k_wf.size())),
false);
3349 for (
unsigned int i = 0; i < k_wf.size(); i++)
3471 d_wf[i] = std::complex<T>(0.0,
k0)*dx_wf +
3472 std::complex<T>(0.0,
k1)*dy_wf +
3473 std::complex<T>(0.0,
k2)*dz_wf;
3476 k_vwf[i] += 0.5 * ksq * k_wf[i];
3477 k_vwf[i] -= d_wf[i];
3483 unsigned int eimax = -1;
3484 double eimaxval = -1e10;
3485 for (
unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3488 if ((
real(
e(fi,fi)) > 0.0) && (
real(
e(fi,fi)) > eimaxval))
3491 eimaxval =
real(
e(fi,fi));
3495 double eshift = (eimaxval > 0.0) ? eimaxval + 0.1 : 0.0;
3496 for (
unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3500 eigs[ei] =
real(
e(fi,fi));
3502 k_vwf[fi] += (
alpha[ei]-eigs[ei])*k_wf[fi];
3510 _eigF <<
"kpt: " << kp << endl;
3511 _eigF << setfill(
'-') << setw(20) <<
" " << endl;
3512 for (
unsigned int ei = kpoint.begin; ei < kpoint.end; ei++)
3515 sprintf(eigstr,
"%3d%15.10f",ei,
real(eigs[ei]));
3516 _eigF << eigstr << endl;
3518 _eigF <<
"\n\n" << endl;
3520 for (
unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3524 vwf[wi] = k_vwf[fi];
3576 if (
_world.
rank() == 0)
_outputF <<
"Building kinetic energy matrix ...\n\n" << endl;
3596 print(
"POTENTIAL:");
3609 for (
unsigned int fi = kpoint.begin; fi < kpoint.end; ++fi)
3612 for (
unsigned int fj = kpoint.begin; fj < fi; ++fj)
3615 f[fi] -= overlap*
f[fj];
3617 f[fi].scale(1.0/
f[fi].
norm2());
3626 Q.gaxpy(0.2,s,-2.0/3.0);
3627 for (
int i=0; i<s.dim(0); ++i)
Q(i,i) += 1.0;
3628 return Q.scale(15.0/8.0);
3635 int n=s.dim(0),
m=s.dim(1);
3640 for (
int i=0; i<n; ++i) {
3644 else if (
e(i) < tol) {
3645 print(
"Matrix square root: Warning: small eigenvalue ", i,
e(i));
3648 e(i) = 1.0/sqrt(
e(i));
3650 for (
int j=0; j<n; ++j) {
3651 for (
int i=0; i<n; ++i) {
3663 vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
3665 printf(
"orthonormalize: \n");
3666 printf(
"before matrix (after): \n");
3672 printf(
"overlap matrix (after): \n");
3674 for (
unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3692 rm.insert(rm.end(), br.begin(), br.end());
3695 std::vector<double> rnvec = norm2s<valueT,NDIM>(
_world, rm);
3697 for (
unsigned int i = 0; i < rnvec.size(); i++) rnorm += rnvec[i];
3705 for (
unsigned int i = 0; i < rnvec.size(); i++)
3707 _outputF <<
"residual" << i <<
"\t" << rnvec[i] << endl;
3718 std::vector<KPoint> kpoints)
3721 truncate<valueT,NDIM> (
_world, awfs);
3725 truncate<valueT,NDIM> (
_world, bwfs);
3743 for (
unsigned int kp = 0; kp < kpoints.size(); kp++)
3755 truncate<valueT,NDIM>(
_world, awfs);
3756 for (
unsigned int ai = 0; ai < awfs.size(); ai++) {
3757 _phisa[ai] = awfs[ai].scale(1.0 / awfs[ai].
norm2());
3761 truncate<valueT,NDIM>(
_world, bwfs);
3762 for (
unsigned int bi = 0; bi < bwfs.size(); bi++)
3764 _phisb[bi] = bwfs[bi].scale(1.0 / bwfs[bi].
norm2());
3769 for (
unsigned int ia = 0; ia < awfs.size(); ia++)
3784 for (
unsigned int i = 0; i < owfs.size(); i++)
3785 nwfs[i].
gaxpy(1.0 - s, owfs[i], s,
false);
3811 std::vector<double>& occs)
3814 double emax = eps[0];
3816 for (
int i = 0; i < eps.size(); i++)
3818 emax = (eps[i] > emax) ? eps[i] : emax;
3819 emin = (eps[i] < emin) ? eps[i] : emin;
3824 double occmax = 2.0;
3826 double efermi = 0.0;
3831 for (
int it = 0; (it < maxits)&&(!bstop); it++)
3834 efermi = 0.5 * (emax + emin);
3836 double charge = 0.0;
3838 for (
int i = 0; i < eps.size(); i++)
3840 double x = (efermi-eps[i]) * t1;
3843 charge +=
_kpoints[i].weight() * occs[i];
3845 if (fabs(emax-emin) < 1
e-5)
double potential(const coord_3d &r)
Definition: 3dharmonic.cc:132
double q(double t)
Definition: DKops.h:18
real_convolution_3d A(World &world)
Definition: DKops.h:230
This header should include pretty much everything needed for the parallel runtime.
std::complex< double > double_complex
Definition: cfft.h:14
Definition: test_ar.cc:118
Definition: test_ar.cc:141
void read_file(const std::string &filename, bool fractional)
Definition: mentity.cc:289
int natom() const
Definition: mentity.h:112
double total_nuclear_charge() const
Definition: mentity.cc:444
void center()
Moves the center of nuclear charge to the origin.
Definition: mentity.cc:413
const Atom & get_atom(unsigned int i) const
Definition: mentity.cc:369
Definition: electronicstructureapp.h:85
Definition: molecule.h:58
double y
Definition: molecule.h:60
double x
Definition: molecule.h:60
double z
Definition: molecule.h:60
Used to represent one basis function from a shell on a specific center.
Definition: madness/chem/molecularbasis.h:405
void get_coords(double &x, double &y, double &z) const
Definition: madness/chem/molecularbasis.h:450
Contracted Gaussian basis.
Definition: madness/chem/molecularbasis.h:465
double eval_guess_density(const Molecule &molecule, double x, double y, double z) const
Evaluates the guess density.
Definition: madness/chem/molecularbasis.h:642
AtomicBasisFunction get_atomic_basis_function(const Molecule &molecule, size_t ibf) const
Returns the ibf'th atomic basis function.
Definition: madness/chem/molecularbasis.h:596
void read_file(std::string filename)
read the atomic basis set from file
Definition: molecularbasis.cc:118
int nbf(const Molecule &molecule) const
Given a molecule count the number of basis functions.
Definition: madness/chem/molecularbasis.h:620
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
const vec3dT exponent
Definition: solver.h:103
Vector< double, NDIM > coordT
Definition: solver.h:100
double_complex operator()(const coordT &x) const
You should implement this to return f(x)
Definition: solver.h:108
ComplexExp(vec3dT exponent, double_complex coeff)
Definition: solver.h:105
Vector< double, NDIM > vec3dT
Definition: solver.h:101
const double_complex coeff
Definition: solver.h:102
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
A MADNESS functor to compute either x, y, or z.
Definition: solver.h:165
virtual ~DipoleFunctor()
Definition: solver.h:173
DipoleFunctor(int axis)
Definition: solver.h:169
double operator()(const coordT &x) const
Definition: solver.h:170
const int axis
Definition: solver.h:167
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static int get_k()
Returns the default wavelet order.
Definition: funcdefaults.h:266
static void set_thresh(double value)
Sets the default threshold.
Definition: funcdefaults.h:286
static void set_k(int value)
Sets the default wavelet order.
Definition: funcdefaults.h:273
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
static void set_bc(const BoundaryConditions< NDIM > &value)
Sets the default boundary conditions.
Definition: funcdefaults.h:418
static void set_cubic_cell(double lo, double hi)
Sets the user cell to be cubic with each dimension having range [lo,hi].
Definition: funcdefaults.h:461
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition: funcdefaults.h:446
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
A multiresolution adaptive numerical function.
Definition: mra.h:122
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition: mra.h:1099
double norm2() const
Returns the 2-norm of the function ... global sum ... works in either basis.
Definition: mra.h:688
Function< T, NDIM > & scale(const Q q, bool fence=true)
Inplace, scale the function by a constant. No communication except for optional fence.
Definition: mra.h:953
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 > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition: mra.h:775
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition: mra.h:721
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
void unaryop(T(*f)(T))
Inplace unary operation on function values.
Definition: mra.h:895
Definition: eigsolver.h:63
double weight()
Definition: eigsolver.h:78
Definition: potentialmanager.h:55
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
The main class of the periodic DFT solver.
Definition: solver.h:619
Tensor< valueT > tensorT
Definition: solver.h:632
rfunctionT make_lda_potential(World &world, const rfunctionT &arho, const rfunctionT &brho, const rfunctionT &adelrhosq, const rfunctionT &bdelrhosq)
Definition: solver.h:1230
vecfuncT compute_residual(const vecfuncT &awfs, const vecfuncT &bwfs)
Definition: solver.h:3683
std::vector< functionT > vecfuncT
Definition: solver.h:625
int _it
Definition: solver.h:746
void apply_hf_exchange3(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb, double &xc)
Definition: solver.h:2318
void apply_hf_exchange4(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb, double &xc)
Definition: solver.h:2366
World & _world
Definition: solver.h:639
void do_rhs_simple(vecfuncT &wf, vecfuncT &vwf, std::vector< KPoint > kpoints, std::vector< T > &alpha, std::vector< double > &eigs)
Definition: solver.h:3202
MolecularEntity _mentity
Definition: solver.h:714
std::vector< KPoint > genkmesh(unsigned int ngridk0, unsigned ngridk1, unsigned int ngridk2, double koffset0, double koffset1, double koffset2, double R)
Definition: solver.h:890
void init(const std::string &filename)
Definition: solver.h:789
cvecfuncT _aobasisf
Definition: solver.h:750
ctensorT csqrt(const ctensorT &s, double tol=1e-8)
Computes matrix square root (not used any more?)
Definition: solver.h:3634
std::vector< subspaceT > vecsubspaceT
Definition: solver.h:636
vecfuncT _phisa
Definition: solver.h:650
Function< valueT, NDIM > functionT
Definition: solver.h:624
ctensorT matrix_exponential(const ctensorT &A)
Definition: solver.h:2786
Solver(World &world, rfunctionT vnucrhon, vecfuncT phisa, vecfuncT phisb, std::vector< T > eigsa, std::vector< T > eigsb, ElectronicStructureParams params, MolecularEntity mentity)
Definition: solver.h:1289
AtomicBasisSet _aobasis
Definition: solver.h:722
int _nao
Definition: solver.h:754
tensor_complex make_kinetic_matrix(World &world, const vector_complex_function_3d &v, const KPoint &k)
Definition: solver.h:3531
rfunctionT _vnucrhon
Definition: solver.h:646
SeparatedConvolution< T, NDIM > * _cop
Definition: solver.h:698
std::shared_ptr< operatorT> poperatorT
Definition: solver.h:629
void set_occs2(const std::vector< KPoint > &kpoints, const std::vector< double > &eigsa, const std::vector< double > &eigsb, std::vector< double > &occsa, std::vector< double > &occsb)
Definition: solver.h:1893
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:633
std::ofstream _matF
Definition: solver.h:734
std::complex< T > valueT
Definition: solver.h:621
std::vector< T > _eigsa
Definition: solver.h:658
std::vector< KPoint > _kpoints
Definition: solver.h:670
vecfuncT project_ao_basis(World &world, KPoint kpt)
Definition: solver.h:1243
std::vector< tensorT > vectensorT
Definition: solver.h:635
void print_fock_matrix_eigs(const vecfuncT &wf, const vecfuncT &vwf, KPoint kpoint)
Definition: solver.h:2872
std::vector< double > _occsb
Definition: solver.h:678
double calculate_kinetic_energy()
Definition: solver.h:2135
void END_TIMER(World &world, const char *msg)
Definition: solver.h:765
FunctionFactory< T, NDIM > rfactoryT
Definition: solver.h:623
KPoint find_kpt_from_orb(unsigned int idx)
Definition: solver.h:2353
rfunctionT _vnuc
Definition: solver.h:694
ElectronicStructureParams _params
Definition: solver.h:666
Vector< double, NDIM > kvecT
Definition: solver.h:627
std::vector< double > _occsa
Definition: solver.h:674
void print_potential_matrix_eigs(const vecfuncT &wf, const vecfuncT &vwf)
Definition: solver.h:2843
void save_orbitals()
Definition: solver.h:921
void step_restriction(vecfuncT &owfs, vecfuncT &nwfs, int aorb)
Definition: solver.h:3778
FunctionFactory< valueT, NDIM > factoryT
Definition: solver.h:626
void reproject()
Definition: solver.h:2474
void START_TIMER(World &world)
Definition: solver.h:761
std::vector< poperatorT > make_bsh_operators(const std::vector< T > &eigs)
Definition: solver.h:2081
std::ofstream _outputF
Definition: solver.h:730
double _residual
Definition: solver.h:718
void initial_guess()
Initializes alpha and beta mos, occupation numbers, eigenvalues.
Definition: solver.h:1316
virtual ~Solver()
Definition: solver.h:1882
void test_periodicity(const Function< Q, 3 > &f)
Definition: solver.h:1159
void make_nuclear_potential_impl()
Definition: solver.h:1034
Solver(World &world, rfunctionT vnucrhon, vecfuncT phis, std::vector< T > eigs, std::vector< KPoint > kpoints, std::vector< double > occs, ElectronicStructureParams params, MolecularEntity mentity)
Definition: solver.h:1847
rfunctionT _rhob
Definition: solver.h:686
std::vector< T > _eigsb
Definition: solver.h:662
tensorT build_fock_matrix(vecfuncT &psi, vecfuncT &vpsi, KPoint kpoint)
Definition: solver.h:3564
void orthonormalize(vecfuncT &wf, KPoint kpoint)
Definition: solver.h:3660
Solver(World &world, const rfunctionT &vnucrhon, const vecfuncT &phis, const std::vector< T > &eigs, const ElectronicStructureParams ¶ms, MolecularEntity mentity)
Definition: solver.h:1812
rfunctionT compute_rho(const vecfuncT &phis, std::vector< KPoint > kpoints, std::vector< double > occs)
Definition: solver.h:2042
Tensor< double > rtensorT
Definition: solver.h:630
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:634
SeparatedConvolution< T, 3 > operatorT
Definition: solver.h:628
std::ofstream _eigF
Definition: solver.h:738
Tensor< std::complex< double > > ctensorT
Definition: solver.h:631
void apply_hf_exchange(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb)
Definition: solver.h:2406
Subspace< T, NDIM > * _subspace
Definition: solver.h:710
rfunctionT _rho
Definition: solver.h:690
double sss
Definition: solver.h:760
void make_nuclear_potential()
Definition: solver.h:1094
void make_nuclear_charge_density_impl()
Definition: solver.h:1042
rfunctionT _rhoa
Definition: solver.h:682
void gram_schmidt(vecfuncT &f, KPoint kpoint)
Definition: solver.h:3607
vecfuncT _phisb
Definition: solver.h:654
double _maxthresh
Definition: solver.h:726
void solve()
Definition: solver.h:2513
Function< T, NDIM > rfunctionT
Definition: solver.h:622
void load_orbitals()
Definition: solver.h:942
void do_rhs(vecfuncT &wf, vecfuncT &vwf, std::vector< KPoint > kpoints, std::vector< T > &alpha, std::vector< double > &eigs)
Definition: solver.h:2939
tensorT Q3(const tensorT &s)
Given overlap matrix, return rotation with 3rd order error to orthonormalize the vectors.
Definition: solver.h:3624
double ttt
Definition: solver.h:760
void fix_occupations(const std::vector< T > &eps, std::vector< double > &occs)
Definition: solver.h:3810
void update_orbitals(vecfuncT &awfs, vecfuncT &bwfs, std::vector< KPoint > kpoints)
Definition: solver.h:3716
void print_tensor2d(ostream &os, Tensor< Q > t)
Definition: solver.h:2827
std::ofstream _kmeshF
Definition: solver.h:742
Solver(World &world, const std::string &filename)
Definition: solver.h:771
The SubspaceK class is a container class holding previous orbitals and residuals.
Definition: solver.h:193
Tensor< valueT > tensorT
Definition: solver.h:200
World & _world
Definition: solver.h:205
vectensorT _Q
Definition: solver.h:209
Function< valueT, NDIM > functionT
Definition: solver.h:196
void update_subspace(vecfuncT &awfs_new, vecfuncT &bwfs_new, const vecfuncT &awfs_old, const vecfuncT &bwfs_old)
Definition: solver.h:245
std::vector< KPoint > _kpoints
Definition: solver.h:217
std::complex< T > valueT
Definition: solver.h:195
vecsubspaceT _subspace
Definition: solver.h:213
std::vector< tensorT > vectensorT
Definition: solver.h:201
std::vector< subspaceT > vecsubspaceT
Definition: solver.h:202
SubspaceK(World &world, const ElectronicStructureParams ¶ms, const std::vector< KPoint > &kpoints)
Definition: solver.h:231
std::vector< functionT > vecfuncT
Definition: solver.h:197
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:199
ElectronicStructureParams _params
Definition: solver.h:221
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:198
double _residual
Definition: solver.h:225
The SubspaceK class is a container class holding previous orbitals and residuals.
Definition: solver.h:409
tensorT _Q
Definition: solver.h:423
std::vector< KPoint > _kpoints
Definition: solver.h:431
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:414
ElectronicStructureParams _params
Definition: solver.h:435
Subspace(World &world, const ElectronicStructureParams ¶ms)
Definition: solver.h:441
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:415
subspaceT _subspace
Definition: solver.h:427
World & _world
Definition: solver.h:419
std::vector< functionT > vecfuncT
Definition: solver.h:413
Tensor< valueT > tensorT
Definition: solver.h:416
void update_subspace(vecfuncT &awfs_new, vecfuncT &bwfs_new, const vecfuncT &awfs_old, const vecfuncT &bwfs_old, const vecfuncT &rm)
Definition: solver.h:448
Function< valueT, NDIM > functionT
Definition: solver.h:412
std::complex< T > valueT
Definition: solver.h:411
void reproject()
Definition: solver.h:556
scalar_type absmax(long *ind=0) const
Return the absolute maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1754
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
double ky
Definition: solver.h:124
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: solver.h:158
std::vector< coord_3d > specialpt
Definition: solver.h:126
double kx
Definition: solver.h:124
double_complex operator()(const coord_3d &x) const
Definition: solver.h:140
double L
Definition: solver.h:123
WSTAtomicBasisFunctor(const AtomicBasisFunction &aofunc, double L, double kx, double ky, double kz)
Definition: solver.h:128
const AtomicBasisFunction aofunc
Definition: solver.h:122
double kz
Definition: solver.h:124
~WSTAtomicBasisFunctor()
Definition: solver.h:138
void broadcast_serializable(objT &obj, ProcessID root)
Broadcast a serializable object.
Definition: worldgop.h:754
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
void broadcast(void *buf, size_t nbyte, ProcessID root, bool dowork=true, Tag bcast_tag=-1)
Broadcasts bytes from process root while still processing AM & tasks.
Definition: worldgop.cc:173
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:870
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
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition: world.h:328
WorldGopInterface & gop
Global operations.
Definition: world.h:205
An archive for storing local or parallel data wrapping a BinaryFstreamOutputArchive.
Definition: parallel_archive.h:321
static const double R
Definition: csqrt.cc:46
double(* f1)(const coord_3d &)
Definition: derivatives.cc:55
double(* f3)(const coord_3d &)
Definition: derivatives.cc:57
double(* f2)(const coord_3d &)
Definition: derivatives.cc:56
const double delta
Definition: dielectric_external_field.cc:119
bool is_equal(double val1, double val2, double eps)
Definition: esolver.h:190
std::shared_ptr< FunctionFunctorInterface< double, 3 > > rfunctorT
Definition: esolver.h:46
Correspondence between C++ and Fortran types.
const double m
Definition: gfit.cc:199
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
rfunctionT compute_rho_slow(const vecfuncT &phis, std::vector< KPoint > kpoints, std::vector< double > occs)
Compute the electronic density for either a molecular or periodic system.
Definition: solver.h:2010
T stheta_fd(const T &x)
Definition: solver.h:75
void apply_potential(vecfuncT &pfuncsa, vecfuncT &pfuncsb, const vecfuncT &phisa, const vecfuncT &phisb, const rfunctionT &rhoa, const rfunctionT &rhob, const rfunctionT &rho)
Applies the LDA effective potential to each orbital. Currently only lda and spin-polarized is not imp...
Definition: solver.h:2177
Tensor< T > conj_transpose(const Tensor< T > &t)
Returns a new deep copy of the complex conjugate transpose of the input tensor.
Definition: tensor.h:2027
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
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
static const double v
Definition: hatom_sf_dirac.cc:20
static const long oi
Definition: he.cc:16
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
#define max(a, b)
Definition: lda.h:51
#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
@ BC_PERIODIC
Definition: funcdefaults.h:56
Function< TENSOR_RESULT_TYPE(L, R), NDIM > sub(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator- but with optional fence and no automatic compression.
Definition: mra.h:1971
static const char * filename
Definition: legendre.cc:96
std::vector< complex_function_3d > vector_complex_function_3d
Definition: functypedefs.h:86
double abs(double x)
Definition: complexfun.h:48
static double cpu_time()
Returns the cpu time in seconds relative to an arbitrary origin.
Definition: timers.h:127
Vector< double, 3 > coordT
Definition: corepotential.cc:54
void sygv(const Tensor< T > &A, const Tensor< T > &B, int itype, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Generalized real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:1052
response_space scale(response_space a, double b)
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
Function< double, NDIM > abssq(const Function< double_complex, NDIM > &z, bool fence=true)
Returns a new function that is the square of the absolute value of the input.
Definition: mra.h:2713
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
Function< T, NDIM > project(const Function< T, NDIM > &other, int k=FunctionDefaults< NDIM >::get_k(), double thresh=FunctionDefaults< NDIM >::get_thresh(), bool fence=true)
Definition: mra.h:2399
Tensor< T > KAIN(const Tensor< T > &Q, double rcond=1e-12)
Solves non-linear equation using KAIN (returns coefficients to compute next vector)
Definition: solvers.h:98
static double r2(const coord_3d &x)
Definition: smooth.h:45
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:133
void plotdx(const Function< T, NDIM > &f, const char *filename, const Tensor< double > &cell=FunctionDefaults< NDIM >::get_cell(), const std::vector< long > &npt=std::vector< long >(NDIM, 201L), bool binary=true)
Writes an OpenDX format file with a cube/slice of points on a uniform grid.
Definition: mraimpl.h:3448
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
static SeparatedConvolution< double_complex, 3 > PeriodicHFExchangeOperator(World &world, Vector< double, 3 > args, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition: operator.h:1721
static const Slice _(0,-1, 1)
std::vector< Function< T, NDIM > > rot(const std::vector< Function< T, NDIM > > &v, bool do_refine=false, bool fence=true)
shorthand rot operator
Definition: vmra.h:1976
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
NDIM & f
Definition: mra.h:2416
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition: corepotential.cc:55
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition: mra.h:1742
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition: timers.cc:48
std::vector< complex_functionT > cvecfuncT
Definition: SCF.h:80
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition: operator.h:1762
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition: vmra.h:1614
void plot_line(World &world, const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const opT &op)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition: funcplot.h:438
double real(double x)
Definition: complexfun.h:52
static SeparatedConvolution< double, 3 > * BSHOperatorPtr3D(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:2023
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition: vmra.h:689
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition: distpm.cc:46
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition: vmra.h:827
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:969
static double phase(long i)
Definition: twoscale.cc:85
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
static long abs(long a)
Definition: tensor.h:218
XCfunctional xc
Definition: newsolver_lda.cc:53
double Q(double a)
Definition: relops.cc:20
static const double c
Definition: relops.cc:10
static const double thresh
Definition: rk.cc:45
static const long k
Definition: rk.cc:44
Defines interfaces for optimization and non-linear equation solvers.
Definition: electronicstructureparams.h:46
unsigned int maxsub
Definition: electronicstructureparams.h:86
bool spinpol
Definition: electronicstructureparams.h:56
int nio
Definition: electronicstructureparams.h:98
bool fractional
Definition: electronicstructureparams.h:84
double ncharge
Definition: electronicstructureparams.h:103
double koffset2
Definition: electronicstructureparams.h:94
int ngridk1
Definition: electronicstructureparams.h:78
bool ispotential
Definition: electronicstructureparams.h:62
int restart
Definition: electronicstructureparams.h:101
bool canon
Definition: electronicstructureparams.h:90
int maxits
Definition: electronicstructureparams.h:60
double rcriterion
Definition: electronicstructureparams.h:111
int nempty
Definition: electronicstructureparams.h:72
int nelec
Definition: electronicstructureparams.h:50
double koffset1
Definition: electronicstructureparams.h:94
int nbands
Definition: electronicstructureparams.h:76
std::string basis
Definition: electronicstructureparams.h:96
void read_file(const std::string &filename)
Definition: electronicstructureparams.h:170
int waveorder
Definition: electronicstructureparams.h:66
double swidth
Definition: electronicstructureparams.h:105
double sd
Definition: electronicstructureparams.h:115
double lo
Definition: electronicstructureparams.h:54
double thresh
Definition: electronicstructureparams.h:64
bool plotorbs
Definition: electronicstructureparams.h:109
int functional
Definition: electronicstructureparams.h:52
int ngridk0
Definition: electronicstructureparams.h:78
double L
Definition: electronicstructureparams.h:48
bool centered
Definition: electronicstructureparams.h:113
bool print_matrices
Definition: electronicstructureparams.h:107
int solver
Definition: electronicstructureparams.h:92
bool periodic
Definition: electronicstructureparams.h:58
int ngridk2
Definition: electronicstructureparams.h:78
double koffset0
Definition: electronicstructureparams.h:94
Definition: solver.h:1121
double operator()(const coordT &x) const
Definition: solver.h:1127
const bool periodic
Definition: solver.h:1125
const MolecularEntity & mentity
Definition: solver.h:1122
double R
Definition: solver.h:1124
GuessDensity(const MolecularEntity &mentity, const AtomicBasisSet &aobasis, const double &R, const bool &periodic)
Definition: solver.h:1151
const AtomicBasisSet & aobasis
Definition: solver.h:1123
static const double_complex I
Definition: tdse1d.cc:164
Tensor< double > B
Definition: tdse1d.cc:166
static const double eshift
Definition: tdse1d.cc:153
const double alpha
Definition: test_coulomb.cc:54
void e()
Definition: test_sig.cc:75
static const double ky
Definition: testcosine.cc:16
static const double kz
Definition: testcosine.cc:16
static const double kx
Definition: testcosine.cc:16
static const std::size_t NDIM
Definition: testpdiff.cc:42
double k0
Definition: testperiodic.cc:66
double k2
Definition: testperiodic.cc:68
double k1
Definition: testperiodic.cc:67
FLOAT weight(const FLOAT &x)
Definition: y.cc:309