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),
297 for (
int s = 0; s <
m; s++)
301 for (
unsigned int i = 0; i <
k_vm.size(); i++)
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);
352 std::complex<double>
one = std::complex<double>(1.0,0.0);
407 template <
typename T,
int NDIM>
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]);
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);
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++)
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++)
930 for (
unsigned int i = 0; i <
_phisa.size(); i++) ar &
_phisa[i];
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++)
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");
1080 if (
_world.
rank() == 0)
printf(
"Difference in nuclear potential: %15.8f\n\n", t1);
1081 for (
int i = 0; i < 101; i++)
1102 printf(
"Cell parameters\n");
1103 printf(
"------------------------------\n");
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)
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++)
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",
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++)
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",
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++)
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",
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");
1408 if (
_world.
rank() == 0)
print(
"Building kinetic energy matrix ...\n\n");
1422 for (
unsigned int ai = 0;
ai < ao.size();
ai++)
1424 std::vector<long> npt(3,101);
1427 std::vector<long>
npt2(3,101);
1433 for (
int i = 0; i < ao.size(); i++)
1436 sprintf(fname,
"ao_line_%d.txt",i);
1483 if (
_world.
rank() == 0)
print(
"Building potential energy matrix ...\n\n");
1486 for (
int i = 0; i < ao.size(); 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++)
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");
1567 for (
int i = 0; i <
overlap.dim(0); i++)
1569 for (
int j = 0; j <
overlap.dim(1); j++)
1732 for (
int i = 0; i <
kinetic.dim(0); i++)
1734 for (
int j = 0; j <
kinetic.dim(1); j++)
1744 for (
int i = 0; i <
kinetic.dim(0); i++)
1746 for (
int j = 0; j <
kinetic.dim(1); 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++)
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++)
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++)
1911 unsigned int sz =
k_eigsa.size();
1916 std::vector<double>
teigs;
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++)
1942 double t1 =
teigs[i];
1953 printf(
"\nSorted eigenvalues:\n");
1954 for (
unsigned int i = 0; i <
teigs.size(); i++)
1985 for (
unsigned int ik = 0;
ik < kpoints.size();
ik++)
1992 for (
unsigned int ist =
k.begin;
ist <
k.end;
ist++)
2011 std::vector<double>
occs)
2018 for (
unsigned int kp = 0;
kp < kpoints.size();
kp++)
2023 for (
unsigned int j =
kpoint.begin; j <
kpoint.end; j++)
2043 std::vector<double>
occs)
2049 std::vector<rfunctionT>
phisq(phis.size());
2050 for (
unsigned int i=0; i<phis.size(); i++) {
2060 for (
unsigned int kp = 0;
kp < kpoints.size();
kp++)
2065 for (
unsigned int j =
kpoint.begin; j <
kpoint.end; j++)
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++)
2269 std::cout.precision(8);
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);
2322 for (
unsigned int j = 0; j <
phisa.size(); j++)
2332 for (
unsigned int i = j+1; i <
phisa.size(); i++)
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++)
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++)
2572 std::vector<functionT>
pfuncsa =
2574 std::vector<functionT>
pfuncsb =
2589 std::vector<long> npt(3,101);
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 -----");
2642 for (
unsigned int ie = 0;
ie <
eo.dim(0);
ie++)
2666 if (
_world.
rank() == 0) std::cout <<
"applying BSH operator ...\n" << endl;
2677 for (
unsigned int i = 0; i <
tmpa_norms.size(); i++)
2767 std::vector<long> npt(3,101);
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();
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)
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;
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)
2950 for (
unsigned int kp = 0;
kp < kpoints.size();
kp++)
2997 unsigned int nmo =
k_wf.size();
3000 for (
long j = 0; j < nmo; j++)
3005 std::complex<T>
phase =
exp(std::complex<T>(0.0,-
ang));
3007 for (
long i = 0; i < nmo; i++)
3021 for (
long i = 0; i < nmo; i++)
3040 while (ilo < nmo-1) {
3044 if (ihi == nmo-1)
break;
3046 long nclus = ihi - ilo + 1;
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++)
3121 std::complex<T>(0.0,
k1)*
dy_wf +
3122 std::complex<T>(0.0,
k2)*
dz_wf;
3131 unsigned int eimax = -1;
3155 _eigF <<
"kpt: " <<
kp << endl;
3156 _eigF << setfill(
'-') << setw(20) <<
" " << endl;
3164 _eigF <<
"\n\n" << endl;
3173 for (
unsigned int ei = 0;
ei <
e.dim(0);
ei++)
3204 std::vector<KPoint> kpoints,
3205 std::vector<T>&
alpha,
3206 std::vector<double>& eigs)
3213 for (
unsigned int kp = 0;
kp < kpoints.size();
kp++)
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);
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++)
3262 std::complex<T>
phase =
exp(std::complex<T>(0.0,-
ang));
3264 for (
long i = 0; i < nmo; i++)
3278 for (
long i = 0; i < nmo; i++)
3297 while (ilo < nmo-1) {
3301 if (ihi == nmo-1)
break;
3303 long nclus = ihi - ilo + 1;
3349 for (
unsigned int i = 0; i <
k_wf.size(); i++)
3472 std::complex<T>(0.0,
k1)*
dy_wf +
3473 std::complex<T>(0.0,
k2)*
dz_wf;
3483 unsigned int eimax = -1;
3510 _eigF <<
"kpt: " <<
kp << endl;
3511 _eigF << setfill(
'-') << setw(20) <<
" " << endl;
3518 _eigF <<
"\n\n" << endl;
3576 if (
_world.
rank() == 0)
_outputF <<
"Building kinetic energy matrix ...\n\n" << endl;
3596 print(
"POTENTIAL:");
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));
3650 for (
int j=0; j<n; ++j) {
3651 for (
int i=0; i<n; ++i) {
3665 printf(
"orthonormalize: \n");
3666 printf(
"before matrix (after): \n");
3672 printf(
"overlap matrix (after): \n");
3692 rm.insert(rm.end(),
br.begin(),
br.end());
3705 for (
unsigned int i = 0; i <
rnvec.size(); i++)
3718 std::vector<KPoint> kpoints)
3743 for (
unsigned int kp = 0;
kp < kpoints.size();
kp++)
3756 for (
unsigned int ai = 0;
ai <
awfs.size();
ai++) {
3762 for (
unsigned int bi = 0;
bi <
bwfs.size();
bi++)
3769 for (
unsigned int ia = 0;
ia <
awfs.size();
ia++)
3784 for (
unsigned int i = 0; i <
owfs.size(); i++)
3811 std::vector<double>&
occs)
3814 double emax = eps[0];
3816 for (
int i = 0; i < eps.size(); i++)
3831 for (
int it = 0; (it < maxits)&&(!
bstop); it++)
3836 double charge = 0.0;
3838 for (
int i = 0; i < eps.size(); i++)
3840 double x = (
efermi-eps[i]) * t1;
double potential(const coord_3d &r)
Definition 3dharmonic.cc:132
double q(double t)
Definition DKops.h:18
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
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:100
static int get_k()
Returns the default wavelet order.
Definition funcdefaults.h:163
static void set_thresh(double value)
Sets the default threshold.
Definition funcdefaults.h:183
static void set_k(int value)
Sets the default wavelet order.
Definition funcdefaults.h:170
static const double & get_thresh()
Returns the default threshold.
Definition funcdefaults.h:176
static void set_bc(const BoundaryConditions< NDIM > &value)
Sets the default boundary conditions.
Definition funcdefaults.h:319
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition funcdefaults.h:369
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:362
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition funcdefaults.h:347
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:139
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:978
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition mra.h:1124
double thresh() const
Returns value of truncation threshold. No communication.
Definition mra.h:592
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:627
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition mra.h:800
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:1006
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition mra.h:746
void unaryop(T(*f)(T))
Inplace unary operation on function values.
Definition mra.h:920
Definition eigsolver.h:63
double weight()
Definition eigsolver.h:78
Definition potentialmanager.h:57
Convolutions in separated form (including Gaussian)
Definition operator.h:139
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
std::vector< KPoint > genkmesh(unsigned int ngridk0, unsigned ngridk1, unsigned int ngridk2, double koffset0, double koffset1, double koffset2, double R)
Definition solver.h:890
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
std::vector< poperatorT > make_bsh_operators(const std::vector< T > &eigs)
Definition solver.h:2081
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
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::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
A tensor is a multidimensional array.
Definition tensor.h:317
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
A simple, fixed dimension vector.
Definition vector.h:64
double ky
Definition solver.h:124
std::vector< coord_3d > specialpt
Definition solver.h:126
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition solver.h:158
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:320
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition world.h:330
WorldGopInterface & gop
Global operations.
Definition world.h:207
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.
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 rot(x, k)
Definition lookup3.c:72
#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.
constexpr double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
@ BC_PERIODIC
Definition bc.h:52
static const char * filename
Definition legendre.cc:96
std::vector< complex_function_3d > vector_complex_function_3d
Definition functypedefs.h:101
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
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:2745
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
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< 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:2003
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const array_of_bools< 3 > &lattice_sum=FunctionDefaults< 3 >::get_bc().is_periodic(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition operator.h:1818
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 truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
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
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:2078
static double r2(const coord_3d &x)
Definition smooth.h:45
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition vmra.h:133
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
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:3467
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 const Slice _(0,-1, 1)
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:1774
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
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
NDIM & f
Definition mra.h:2448
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition corepotential.cc:55
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
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:1680
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
static SeparatedConvolution< double, 3 > * BSHOperatorPtr3D(World &world, double mu, double lo, double eps, const array_of_bools< 3 > &lattice_sum=FunctionDefaults< 3 >::get_bc().is_periodic(), 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:2089
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:2431
double real(double x)
Definition complexfun.h:52
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
static SeparatedConvolution< double_complex, 3 > PeriodicHFExchangeOperator(World &world, Vector< double, 3 > args, double lo, double eps, const array_of_bools< 3 > &lattice_sum=FunctionDefaults< 3 >::get_bc().is_periodic(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition operator.h:1775
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
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:2034
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:246
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 m
Definition relops.cc:9
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
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
static const double eshift
Definition tdse1d.cc:153
AtomicInt sum
Definition test_atomicint.cc:46
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 double alpha
Definition testcosine.cc:10
#define START_TIMER
Definition testgaxpyext.cc:9
#define END_TIMER(msg)
Definition testgaxpyext.cc:10
constexpr coord_t one(1.0)
constexpr std::size_t NDIM
Definition testgconv.cc:54
double g1(const coord_t &r)
Definition testgconv.cc:122
const auto npts
Definition testgconv.cc:52
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