45 double exactexchange = 0.0;
48 std::vector<std::string> functional_names =
51 if (functional_names.size() > 2) {
52 throw std::runtime_error(
"Too many functional names");
53 }
else if (functional_names.size() < 1) {
54 throw std::runtime_error(
"Specify at least one functional");
57 for (
const std::string& functional_name : functional_names) {
59 int func_id = map.getID(functional_name);
65 if (xc_func_init(&func, func_id, XC_UNPOLARIZED) != 0) {
66 throw std::runtime_error(
67 (boost::format(
"Functional %s not found\n") % functional_name).str());
69 if (exactexchange > 0 && func.cam_alpha > 0) {
70 throw std::runtime_error(
71 "You have specified two functionals with exact exchange");
73 exactexchange += func.cam_alpha;
160 const Eigen::MatrixXd& density_matrix)
const {
162 assert(density_matrix.isApprox(density_matrix.transpose()) &&
163 "Density matrix has to be symmetric!");
166#pragma omp parallel for schedule(guided) reduction(+ : vxc)
167 for (
Index i = 0; i <
grid_.getBoxesSize(); ++i) {
172 double EXC_box = 0.0;
176 1.e-40 / double(density_matrix.rows()) / double(density_matrix.rows());
177 if (DMAT_here.cwiseAbs2().maxCoeff() < cutoff) {
180 Eigen::MatrixXd Vxc_here =
181 Eigen::MatrixXd::Zero(DMAT_here.rows(), DMAT_here.cols());
182 const std::vector<Eigen::Vector3d>& points = box.
getGridPoints();
188 Eigen::VectorXd temp = ao.
values.transpose() * DMAT_here;
189 double rho = 0.5 * temp.dot(ao.
values);
190 const double weight = weights[p];
191 if (rho * weight < 1.e-20) {
194 const Eigen::Vector3d rho_grad = temp.transpose() * ao.
derivatives;
198 EXC_box += weight * rho * xc.
f_xc;
202 Vxc_here.noalias() += temp * ao.
values.transpose();