1 /* 2 * This source code is part of 3 * 4 * HelFEM 5 * - 6 * Finite element methods for electronic structure calculations on small systems 7 * 8 * Written by Susi Lehtola, 2018- 9 * Copyright (c) 2018- Susi Lehtola 10 * 11 * This program is free software; you can redistribute it and/or 12 * modify it under the terms of the GNU General Public License 13 * as published by the Free Software Foundation; either version 2 14 * of the License, or (at your option) any later version. 15 */ 16 17 #include <cfloat> 18 #include <cmath> 19 #include <cstdio> 20 // LibXC 21 #include <xc.h> 22 23 #include "dftgrid.h" 24 #include "../general/dftfuncs.h" 25 // Angular quadrature 26 #include "../general/angular.h" 27 28 // OpenMP parallellization for XC calculations 29 #ifdef _OPENMP 30 #include <omp.h> 31 #endif 32 33 namespace helfem { 34 namespace diatomic { 35 namespace dftgrid { DFTGridWorker()36 DFTGridWorker::DFTGridWorker() { 37 } 38 DFTGridWorker(const helfem::diatomic::basis::TwoDBasis * basp_,int lang,int mang)39 DFTGridWorker::DFTGridWorker(const helfem::diatomic::basis::TwoDBasis * basp_, int lang, int mang) : basp(basp_) { 40 do_grad=false; 41 do_tau=false; 42 do_lapl=false; 43 44 // Get angular grid 45 helfem::angular::angular_chebyshev(lang,mang,cth,phi,wang); 46 } 47 ~DFTGridWorker()48 DFTGridWorker::~DFTGridWorker() { 49 } 50 update_density(const arma::mat & P0)51 void DFTGridWorker::update_density(const arma::mat & P0) { 52 // Update values of density 53 if(!P0.n_elem) { 54 throw std::runtime_error("Error - density matrix is empty!\n"); 55 } 56 arma::mat P(basp->expand_boundaries(P0)(bf_ind,bf_ind)); 57 58 // Non-polarized calculation. 59 polarized=false; 60 61 // Update density vector 62 Pv=P*arma::conj(bf); 63 64 // Calculate density 65 rho.zeros(1,wtot.n_elem); 66 #ifdef _OPENMP 67 #pragma omp parallel for 68 #endif 69 for(size_t ip=0;ip<wtot.n_elem;ip++) 70 rho(0,ip)=std::real(arma::dot(Pv.col(ip),bf.col(ip))); 71 72 // Calculate gradient 73 if(do_grad) { 74 grho.zeros(3,wtot.n_elem); 75 sigma.zeros(1,wtot.n_elem); 76 #ifdef _OPENMP 77 #pragma omp parallel for 78 #endif 79 for(size_t ip=0;ip<wtot.n_elem;ip++) { 80 // Calculate values 81 double g_rad=grho(0,ip)=2.0*std::real(arma::dot(Pv.col(ip),bf_rho.col(ip)))/scale_r(ip); 82 double g_th=grho(1,ip)=2.0*std::real(arma::dot(Pv.col(ip),bf_theta.col(ip)))/scale_theta(ip); 83 double g_phi=grho(2,ip)=2.0*std::real(arma::dot(Pv.col(ip),bf_phi.col(ip)))/scale_phi(ip); 84 // Compute sigma as well 85 sigma(0,ip)=g_rad*g_rad + g_th*g_th + g_phi*g_phi; 86 } 87 } 88 89 // Calculate laplacian and kinetic energy density 90 if(do_tau) { 91 // Adjust size of grid 92 tau.zeros(1,wtot.n_elem); 93 94 // Update helpers 95 Pv_rho=P*arma::conj(bf_rho); 96 Pv_theta=P*arma::conj(bf_theta); 97 Pv_phi=P*arma::conj(bf_phi); 98 99 // Calculate values 100 #ifdef _OPENMP 101 #pragma omp parallel for 102 #endif 103 for(size_t ip=0;ip<wtot.n_elem;ip++) { 104 // Gradient term 105 double kinrho(std::real(arma::dot(Pv_rho.col(ip),bf_rho.col(ip)))/std::pow(scale_r(ip),2)); 106 double kintheta(std::real(arma::dot(Pv_theta.col(ip),bf_theta.col(ip)))/std::pow(scale_theta(ip),2)); 107 double kinphi(std::real(arma::dot(Pv_phi.col(ip),bf_phi.col(ip)))/std::pow(scale_phi(ip),2)); 108 double kin(kinrho + kintheta + kinphi); 109 110 // Store values 111 tau(0,ip)=0.5*kin; 112 } 113 } 114 115 if(do_lapl) 116 throw std::logic_error("Laplacian not implemented!\n"); 117 } 118 update_density(const arma::mat & Pa0,const arma::mat & Pb0)119 void DFTGridWorker::update_density(const arma::mat & Pa0, const arma::mat & Pb0) { 120 if(!Pa0.n_elem || !Pb0.n_elem) { 121 throw std::runtime_error("Error - density matrix is empty!\n"); 122 } 123 124 // Polarized calculation. 125 polarized=true; 126 127 // Update density vector 128 arma::mat Pa(basp->expand_boundaries(Pa0)(bf_ind,bf_ind)); 129 arma::mat Pb(basp->expand_boundaries(Pb0)(bf_ind,bf_ind)); 130 131 Pav=Pa*arma::conj(bf); 132 Pbv=Pb*arma::conj(bf); 133 134 // Calculate density 135 rho.zeros(2,wtot.n_elem); 136 #ifdef _OPENMP 137 #pragma omp parallel for 138 #endif 139 for(size_t ip=0;ip<wtot.n_elem;ip++) { 140 rho(0,ip)=std::real(arma::dot(Pav.col(ip),bf.col(ip))); 141 rho(1,ip)=std::real(arma::dot(Pbv.col(ip),bf.col(ip))); 142 143 /* 144 double na=compute_density(Pa0,*basp,grid[ip].r); 145 double nb=compute_density(Pb0,*basp,grid[ip].r); 146 if(fabs(da-na)>1e-6 || fabs(db-nb)>1e-6) 147 printf("Density at point % .3f % .3f % .3f: %e vs %e, %e vs %e\n",grid[ip].r.x,grid[ip].r.y,grid[ip].r.z,da,na,db,nb); 148 */ 149 } 150 151 // Calculate gradient 152 153 if(do_grad) { 154 grho.zeros(6,wtot.n_elem); 155 sigma.zeros(3,wtot.n_elem); 156 #ifdef _OPENMP 157 #pragma omp parallel for 158 #endif 159 for(size_t ip=0;ip<wtot.n_elem;ip++) { 160 double ga_rad=grho(0,ip)=2.0*std::real(arma::dot(Pav.col(ip),bf_rho.col(ip)))/scale_r(ip); 161 double ga_th=grho(1,ip)=2.0*std::real(arma::dot(Pav.col(ip),bf_theta.col(ip)))/scale_theta(ip); 162 double ga_phi=grho(2,ip)=2.0*std::real(arma::dot(Pav.col(ip),bf_phi.col(ip)))/scale_phi(ip); 163 164 double gb_rad=grho(3,ip)=2.0*std::real(arma::dot(Pbv.col(ip),bf_rho.col(ip)))/scale_r(ip); 165 double gb_th=grho(4,ip)=2.0*std::real(arma::dot(Pbv.col(ip),bf_theta.col(ip)))/scale_theta(ip); 166 double gb_phi=grho(5,ip)=2.0*std::real(arma::dot(Pbv.col(ip),bf_phi.col(ip)))/scale_phi(ip); 167 168 // Compute sigma as well 169 sigma(0,ip)=ga_rad*ga_rad + ga_th*ga_th + ga_phi*ga_phi; 170 sigma(1,ip)=ga_rad*gb_rad + ga_th*gb_th + ga_phi*gb_phi; 171 sigma(2,ip)=gb_rad*gb_rad + gb_th*gb_th + gb_phi*gb_phi; 172 } 173 } 174 175 // Calculate kinetic energy density 176 if(do_tau) { 177 // Adjust size of grid 178 tau.resize(2,wtot.n_elem); 179 180 // Update helpers 181 Pav_rho=Pa*arma::conj(bf_rho); 182 Pav_theta=Pa*arma::conj(bf_theta); 183 Pav_phi=Pa*arma::conj(bf_phi); 184 185 Pbv_rho=Pb*arma::conj(bf_rho); 186 Pbv_theta=Pb*arma::conj(bf_theta); 187 Pbv_phi=Pb*arma::conj(bf_phi); 188 189 // Calculate values 190 #ifdef _OPENMP 191 #pragma omp parallel for 192 #endif 193 for(size_t ip=0;ip<wtot.n_elem;ip++) { 194 // Gradient term 195 double kinar=std::real(arma::dot(Pav_rho.col(ip),bf_rho.col(ip)))/std::pow(scale_r(ip),2); 196 double kinath=std::real(arma::dot(Pav_theta.col(ip),bf_theta.col(ip)))/std::pow(scale_theta(ip),2); 197 double kinaphi=std::real(arma::dot(Pav_phi.col(ip),bf_phi.col(ip)))/std::pow(scale_phi(ip),2); 198 double kina(kinar + kinath + kinaphi); 199 200 double kinbr=std::real(arma::dot(Pbv_rho.col(ip),bf_rho.col(ip)))/std::pow(scale_r(ip),2); 201 double kinbth=std::real(arma::dot(Pbv_theta.col(ip),bf_theta.col(ip)))/std::pow(scale_theta(ip),2); 202 double kinbphi=std::real(arma::dot(Pbv_phi.col(ip),bf_phi.col(ip)))/std::pow(scale_phi(ip),2); 203 double kinb(kinbr + kinbth + kinbphi); 204 205 // Store values 206 tau(0,ip)=0.5*kina; 207 tau(1,ip)=0.5*kinb; 208 } 209 if(do_lapl) 210 throw std::logic_error("Laplacian not implemented!\n"); 211 } 212 } 213 screen_density(double thr)214 void DFTGridWorker::screen_density(double thr) { 215 if(polarized) { 216 #ifdef _OPENMP 217 #pragma omp parallel for 218 #endif 219 for(size_t ip=0;ip<wtot.n_elem;ip++) { 220 if(rho(0,ip)+rho(1,ip) <= thr) { 221 rho(0,ip)=0.0; 222 rho(1,ip)=0.0; 223 224 if(do_grad) { 225 sigma(0,ip)=0.0; 226 sigma(1,ip)=0.0; 227 sigma(2,ip)=0.0; 228 } 229 230 if(do_tau) { 231 tau(0,ip)=0.0; 232 tau(1,ip)=0.0; 233 } 234 } 235 } 236 } else { 237 #ifdef _OPENMP 238 #pragma omp parallel for 239 #endif 240 for(size_t ip=0;ip<wtot.n_elem;ip++) { 241 if(rho(0,ip) <= thr) { 242 rho(0,ip)=0.0; 243 if(do_grad) { 244 sigma(0,ip)=0.0; 245 } 246 if(do_tau) { 247 tau(0,ip)=0.0; 248 } 249 } 250 } 251 } 252 } 253 compute_Nel() const254 double DFTGridWorker::compute_Nel() const { 255 double nel=0.0; 256 if(!polarized) { 257 for(size_t ip=0;ip<wtot.n_elem;ip++) 258 nel+=wtot(ip)*rho(0,ip); 259 } else { 260 for(size_t ip=0;ip<wtot.n_elem;ip++) 261 nel+=wtot(ip)*(rho(0,ip)+rho(1,ip)); 262 } 263 264 return nel; 265 } 266 compute_Ekin() const267 double DFTGridWorker::compute_Ekin() const { 268 double ekin=0.0; 269 270 if(do_tau) { 271 if(!polarized) { 272 for(size_t ip=0;ip<wtot.n_elem;ip++) 273 ekin+=wtot(ip)*tau(0,ip); 274 } else { 275 for(size_t ip=0;ip<wtot.n_elem;ip++) 276 ekin+=wtot(ip)*(tau(0,ip)+tau(1,ip)); 277 } 278 } 279 return ekin; 280 } 281 init_xc()282 void DFTGridWorker::init_xc() { 283 // Size of grid. 284 const size_t N=wtot.n_elem; 285 286 // Zero energy 287 zero_Exc(); 288 289 if(!polarized) { 290 // Restricted case 291 vxc.zeros(1,N); 292 if(do_grad) 293 vsigma.zeros(1,N); 294 if(do_tau) 295 vtau.zeros(1,N); 296 if(do_lapl) 297 vlapl.zeros(1,N); 298 } else { 299 // Unrestricted case 300 vxc.zeros(2,N); 301 if(do_grad) 302 vsigma.zeros(3,N); 303 if(do_tau) 304 vtau.zeros(2,N); 305 if(do_lapl) { 306 vlapl.zeros(2,N); 307 } 308 } 309 310 // Initial values 311 do_gga=false; 312 do_mgga_l=false; 313 do_mgga_t=false; 314 } 315 zero_Exc()316 void DFTGridWorker::zero_Exc() { 317 exc.zeros(wtot.n_elem); 318 } 319 check_xc()320 void DFTGridWorker::check_xc() { 321 size_t inf=0; 322 323 for(arma::uword i=0;i<exc.n_elem;i++) 324 if(!std::isfinite(exc[i])) { 325 inf++; 326 exc[i]=0.0; 327 } 328 329 for(arma::uword i=0;i<vxc.n_elem;i++) 330 if(!std::isfinite(vxc[i])) { 331 inf++; 332 vxc[i]=0.0; 333 } 334 335 for(arma::uword i=0;i<vsigma.n_elem;i++) 336 if(!std::isfinite(vsigma[i])) { 337 inf++; 338 vsigma[i]=0.0; 339 } 340 341 for(arma::uword i=0;i<vlapl.n_elem;i++) 342 if(!std::isfinite(vlapl[i])) { 343 inf++; 344 vlapl[i]=0.0; 345 } 346 347 348 for(arma::uword i=0;i<vtau.n_elem;i++) 349 if(!std::isfinite(vtau[i])) { 350 inf++; 351 vtau[i]=0.0; 352 } 353 354 if(inf) { 355 printf("Warning - %i non-finite entries found in xc energy / potential.\n",(int) inf); 356 } 357 } 358 check_array(const std::vector<double> & x,size_t n,std::vector<size_t> & idx)359 void check_array(const std::vector<double> & x, size_t n, std::vector<size_t> & idx) { 360 if(x.size()%n!=0) { 361 std::ostringstream oss; 362 oss << "Size of array " << x.size() << " is not divisible by " << n << "!\n"; 363 throw std::runtime_error(oss.str()); 364 } 365 366 for(size_t i=0;i<x.size()/n;i++) { 367 // Check for failed entry 368 bool fail=false; 369 for(size_t j=0;j<n;j++) 370 if(!std::isfinite(x[i*n+j])) 371 fail=true; 372 373 // If failed i is not in the list, add it 374 if(fail) { 375 if (!std::binary_search (idx.begin(), idx.end(), i)) { 376 idx.push_back(i); 377 std::sort(idx.begin(),idx.end()); 378 } 379 } 380 } 381 } 382 compute_xc(int func_id,const arma::vec & p,bool pot)383 void DFTGridWorker::compute_xc(int func_id, const arma::vec & p, bool pot) { 384 // Compute exchange-correlation functional 385 386 // Which functional is in question? 387 bool gga, mgga_t, mgga_l; 388 is_gga_mgga(func_id,gga,mgga_t,mgga_l); 389 390 // Update controlling flags for eval_Fxc (exchange and correlation 391 // parts might be of different type) 392 do_gga=do_gga || gga || mgga_t || mgga_l; 393 do_mgga_t=do_mgga_t || mgga_t; 394 do_mgga_l=do_mgga_l || mgga_l; 395 396 // Amount of grid points 397 const size_t N=wtot.n_elem; 398 399 // Work arrays - exchange and correlation are computed separately 400 arma::rowvec exc_wrk; 401 arma::mat vxc_wrk; 402 arma::mat vsigma_wrk; 403 arma::mat vlapl_wrk; 404 arma::mat vtau_wrk; 405 406 if(has_exc(func_id)) 407 exc_wrk.zeros(exc.n_elem); 408 if(pot) { 409 vxc_wrk.zeros(vxc.n_rows,vxc.n_cols); 410 if(gga || mgga_t || mgga_l) 411 vsigma_wrk.zeros(vsigma.n_rows,vsigma.n_cols); 412 if(mgga_t) 413 vtau_wrk.zeros(vtau.n_rows,vtau.n_cols); 414 if(mgga_l) 415 vlapl_wrk.zeros(vlapl.n_rows,vlapl.n_cols); 416 } 417 418 // Spin variable for libxc 419 int nspin; 420 if(!polarized) 421 nspin=XC_UNPOLARIZED; 422 else 423 nspin=XC_POLARIZED; 424 425 // Initialize libxc worker 426 xc_func_type func; 427 if(xc_func_init(&func, func_id, nspin) != 0) { 428 std::ostringstream oss; 429 oss << "Functional "<<func_id<<" not found!"; 430 throw std::runtime_error(oss.str()); 431 } 432 433 // Set parameters 434 if(p.n_elem) { 435 // Check sanity 436 if(p.n_elem != (arma::uword) xc_func_info_get_n_ext_params((xc_func_info_type *)func.info)) 437 throw std::logic_error("Incompatible number of parameters!\n"); 438 arma::vec phlp(p); 439 xc_func_set_ext_params(&func, phlp.memptr()); 440 } 441 442 // Evaluate functionals. 443 if(has_exc(func_id)) { 444 if(pot) { 445 if(mgga_t || mgga_l) {// meta-GGA 446 double * laplp = mgga_l ? lapl.memptr() : NULL; 447 double * taup = mgga_t ? tau.memptr() : NULL; 448 double * vlaplp = mgga_t ? vlapl_wrk.memptr() : NULL; 449 double * vtaup = mgga_t ? vtau_wrk.memptr() : NULL; 450 xc_mgga_exc_vxc(&func, N, rho.memptr(), sigma.memptr(), laplp, taup, exc_wrk.memptr(), vxc_wrk.memptr(), vsigma_wrk.memptr(), vlaplp, vtaup); 451 } else if(gga) // GGA 452 xc_gga_exc_vxc(&func, N, rho.memptr(), sigma.memptr(), exc_wrk.memptr(), vxc_wrk.memptr(), vsigma_wrk.memptr()); 453 else // LDA 454 xc_lda_exc_vxc(&func, N, rho.memptr(), exc_wrk.memptr(), vxc_wrk.memptr()); 455 } else { 456 if(mgga_t || mgga_l) { // meta-GGA 457 double * laplp = mgga_l ? lapl.memptr() : NULL; 458 double * taup = mgga_t ? tau.memptr() : NULL; 459 xc_mgga_exc(&func, N, rho.memptr(), sigma.memptr(), laplp, taup, exc_wrk.memptr()); 460 } else if(gga) // GGA 461 xc_gga_exc(&func, N, rho.memptr(), sigma.memptr(), exc_wrk.memptr()); 462 else // LDA 463 xc_lda_exc(&func, N, rho.memptr(), exc_wrk.memptr()); 464 } 465 466 } else { 467 if(pot) { 468 if(mgga_t || mgga_l) { // meta-GGA 469 double * laplp = mgga_l ? lapl.memptr() : NULL; 470 double * taup = mgga_t ? tau.memptr() : NULL; 471 double * vlaplp = mgga_l ? vlapl_wrk.memptr() : NULL; 472 double * vtaup = mgga_t ? vtau_wrk.memptr() : NULL; 473 xc_mgga_vxc(&func, N, rho.memptr(), sigma.memptr(), laplp, taup, vxc_wrk.memptr(), vsigma_wrk.memptr(), vlaplp, vtaup); 474 } else if(gga) // GGA 475 xc_gga_vxc(&func, N, rho.memptr(), sigma.memptr(), vxc_wrk.memptr(), vsigma_wrk.memptr()); 476 else // LDA 477 xc_lda_vxc(&func, N, rho.memptr(), vxc_wrk.memptr()); 478 } 479 } 480 481 // Free functional 482 xc_func_end(&func); 483 484 // Sum to total arrays containing both exchange and correlation 485 if(has_exc(func_id)) 486 exc+=exc_wrk; 487 if(pot) { 488 if(mgga_l) 489 vlapl+=vlapl_wrk; 490 if(mgga_t) 491 vtau+=vtau_wrk; 492 if(mgga_t || mgga_l || gga) 493 vsigma+=vsigma_wrk; 494 vxc+=vxc_wrk; 495 } 496 } 497 save(const std::string & info) const498 void DFTGridWorker::save(const std::string & info) const { 499 rho.save("rho"+info+".dat",arma::raw_ascii); 500 sigma.save("sigma"+info+".dat",arma::raw_ascii); 501 tau.save("tau"+info+".dat",arma::raw_ascii); 502 503 vxc.save("vxc"+info+".dat",arma::raw_ascii); 504 vsigma.save("vsigma"+info+".dat",arma::raw_ascii); 505 vtau.save("vtau"+info+".dat",arma::raw_ascii); 506 exc.save("exc"+info+".dat",arma::raw_ascii); 507 } 508 eval_Exc() const509 double DFTGridWorker::eval_Exc() const { 510 arma::rowvec dens(rho.row(0)); 511 if(polarized) 512 dens+=rho.row(1); 513 514 return arma::sum(wtot%exc%dens); 515 } 516 eval_overlap(arma::mat & So) const517 void DFTGridWorker::eval_overlap(arma::mat & So) const { 518 // Calculate in subspace 519 arma::mat S(bf_ind.n_elem,bf_ind.n_elem); 520 S.zeros(); 521 increment_lda< std::complex<double> >(S,wtot,bf); 522 // Increment 523 So.submat(bf_ind,bf_ind)+=S; 524 } 525 eval_kinetic(arma::mat & To) const526 void DFTGridWorker::eval_kinetic(arma::mat & To) const { 527 // Calculate in subspace 528 arma::mat T(bf_ind.n_elem,bf_ind.n_elem); 529 T.zeros(); 530 increment_lda< std::complex<double> >(T,wtot/(scale_r%scale_r),bf_rho); 531 increment_lda< std::complex<double> >(T,wtot/(scale_theta%scale_theta),bf_theta); 532 increment_lda< std::complex<double> >(T,wtot/(scale_phi%scale_phi),bf_phi); 533 // Increment 534 To.submat(bf_ind,bf_ind)+=0.5*T; 535 } 536 eval_Fxc(arma::mat & Ho) const537 void DFTGridWorker::eval_Fxc(arma::mat & Ho) const { 538 if(polarized) { 539 throw std::runtime_error("Refusing to compute restricted Fock matrix with unrestricted density.\n"); 540 } 541 542 // Work matrix 543 arma::mat H(bf_ind.n_elem,bf_ind.n_elem); 544 H.zeros(); 545 546 { 547 // LDA potential 548 arma::rowvec vrho(vxc.row(0)); 549 // Multiply weights into potential 550 vrho%=wtot; 551 // Increment matrix 552 increment_lda< std::complex<double> >(H,vrho,bf); 553 } 554 555 if(do_gga) { 556 // Get vsigma 557 arma::rowvec vs(vsigma.row(0)); 558 // Get grad rho 559 arma::uvec idx(arma::linspace<arma::uvec>(0,2,3)); 560 arma::mat gr(arma::trans(grho.rows(idx))); 561 // Multiply grad rho by vsigma and the weights 562 for(size_t i=0;i<gr.n_rows;i++) { 563 gr(i,0)*=2.0*wtot(i)*vs(i)/scale_r(i); 564 gr(i,1)*=2.0*wtot(i)*vs(i)/scale_theta(i); 565 gr(i,2)*=2.0*wtot(i)*vs(i)/scale_phi(i); 566 } 567 // Increment matrix 568 increment_gga< std::complex<double> >(H,gr,bf,bf_rho,bf_theta,bf_phi); 569 } 570 571 if(do_mgga_t) { 572 arma::rowvec vt(vtau.row(0)); 573 vt%=0.5*wtot; 574 575 increment_lda< std::complex<double> >(H,vt/arma::square(scale_r),bf_rho); 576 increment_lda< std::complex<double> >(H,vt/arma::square(scale_theta),bf_theta); 577 increment_lda< std::complex<double> >(H,vt/arma::square(scale_phi),bf_phi); 578 } 579 if(do_mgga_l) 580 throw std::logic_error("Laplacian not implemented!\n"); 581 582 Ho(bf_ind,bf_ind)+=H; 583 } 584 eval_Fxc(arma::mat & Hao,arma::mat & Hbo,bool beta) const585 void DFTGridWorker::eval_Fxc(arma::mat & Hao, arma::mat & Hbo, bool beta) const { 586 if(!polarized) { 587 throw std::runtime_error("Refusing to compute unrestricted Fock matrix with restricted density.\n"); 588 } 589 590 arma::mat Ha, Hb; 591 Ha.zeros(bf_ind.n_elem,bf_ind.n_elem); 592 if(beta) 593 Hb.zeros(bf_ind.n_elem,bf_ind.n_elem); 594 595 { 596 // LDA potential 597 arma::rowvec vrhoa(vxc.row(0)); 598 // Multiply weights into potential 599 vrhoa%=wtot; 600 // Increment matrix 601 increment_lda< std::complex<double> >(Ha,vrhoa,bf); 602 603 if(beta) { 604 arma::rowvec vrhob(vxc.row(1)); 605 vrhob%=wtot; 606 increment_lda< std::complex<double> >(Hb,vrhob,bf); 607 } 608 } 609 if(Ha.has_nan() || (beta && Hb.has_nan())) 610 //throw std::logic_error("NaN encountered!\n"); 611 fprintf(stderr,"NaN in Hamiltonian!\n"); 612 613 if(do_gga) { 614 // Get vsigma 615 arma::rowvec vs_aa(vsigma.row(0)); 616 arma::rowvec vs_ab(vsigma.row(1)); 617 618 // Get grad rho 619 arma::uvec idxa(arma::linspace<arma::uvec>(0,2,3)); 620 arma::uvec idxb(arma::linspace<arma::uvec>(3,5,3)); 621 arma::mat gr_a0(arma::trans(grho.rows(idxa))); 622 arma::mat gr_b0(arma::trans(grho.rows(idxb))); 623 624 // Multiply grad rho by vsigma and the weights 625 arma::mat gr_a(gr_a0); 626 for(size_t i=0;i<gr_a.n_rows;i++) { 627 gr_a(i,0)=wtot(i)*(2.0*vs_aa(i)*gr_a0(i,0) + vs_ab(i)*gr_b0(i,0))/scale_r(i); 628 gr_a(i,1)=wtot(i)*(2.0*vs_aa(i)*gr_a0(i,1) + vs_ab(i)*gr_b0(i,1))/scale_theta(i); 629 gr_a(i,2)=wtot(i)*(2.0*vs_aa(i)*gr_a0(i,2) + vs_ab(i)*gr_b0(i,2))/scale_phi(i); 630 } 631 // Increment matrix 632 increment_gga< std::complex<double> >(Ha,gr_a,bf,bf_rho,bf_theta,bf_phi); 633 634 if(beta) { 635 arma::rowvec vs_bb(vsigma.row(2)); 636 arma::mat gr_b(gr_b0); 637 for(size_t i=0;i<gr_b.n_rows;i++) { 638 gr_b(i,0)=wtot(i)*(2.0*vs_bb(i)*gr_b0(i,0) + vs_ab(i)*gr_a0(i,0))/scale_r(i); 639 gr_b(i,1)=wtot(i)*(2.0*vs_bb(i)*gr_b0(i,1) + vs_ab(i)*gr_a0(i,1))/scale_theta(i); 640 gr_b(i,2)=wtot(i)*(2.0*vs_bb(i)*gr_b0(i,2) + vs_ab(i)*gr_a0(i,2))/scale_phi(i); 641 } 642 increment_gga< std::complex<double> >(Hb,gr_b,bf,bf_rho,bf_theta,bf_phi); 643 } 644 } 645 646 647 if(do_mgga_t) { 648 arma::rowvec vt_a(vtau.row(0)); 649 vt_a%=0.5*wtot; 650 651 increment_lda< std::complex<double> >(Ha,vt_a/arma::square(scale_r),bf_rho); 652 increment_lda< std::complex<double> >(Ha,vt_a/arma::square(scale_theta),bf_theta); 653 increment_lda< std::complex<double> >(Ha,vt_a/arma::square(scale_phi),bf_phi); 654 if(beta) { 655 arma::rowvec vt_b(vtau.row(1)); 656 vt_b%=0.5*wtot; 657 658 increment_lda< std::complex<double> >(Hb,vt_b/arma::square(scale_r),bf_rho); 659 increment_lda< std::complex<double> >(Hb,vt_b/arma::square(scale_theta),bf_theta); 660 increment_lda< std::complex<double> >(Hb,vt_b/arma::square(scale_phi),bf_phi); 661 } 662 } 663 if(do_mgga_l) { 664 throw std::logic_error("Laplacian not implemented!\n"); 665 } 666 667 Hao(bf_ind,bf_ind)+=Ha; 668 if(beta) 669 Hbo(bf_ind,bf_ind)+=Hb; 670 } 671 check_grad_tau_lapl(int x_func,int c_func)672 void DFTGridWorker::check_grad_tau_lapl(int x_func, int c_func) { 673 // Do we need gradients? 674 do_grad=false; 675 if(x_func>0) 676 do_grad=do_grad || gradient_needed(x_func); 677 if(c_func>0) 678 do_grad=do_grad || gradient_needed(c_func); 679 680 // Do we need laplacians? 681 do_tau=false; 682 if(x_func>0) 683 do_tau=do_tau || tau_needed(x_func); 684 if(c_func>0) 685 do_tau=do_tau || tau_needed(c_func); 686 687 // Do we need laplacians? 688 do_lapl=false; 689 if(x_func>0) 690 do_lapl=do_lapl || laplacian_needed(x_func); 691 if(c_func>0) 692 do_lapl=do_lapl || laplacian_needed(c_func); 693 } 694 get_grad_tau_lapl(bool & grad_,bool & tau_,bool & lap_) const695 void DFTGridWorker::get_grad_tau_lapl(bool & grad_, bool & tau_, bool & lap_) const { 696 grad_=do_grad; 697 tau_=do_tau; 698 lap_=do_lapl; 699 } 700 set_grad_tau_lapl(bool grad_,bool tau_,bool lap_)701 void DFTGridWorker::set_grad_tau_lapl(bool grad_, bool tau_, bool lap_) { 702 do_grad=grad_; 703 do_tau=tau_; 704 do_lapl=lap_; 705 } 706 compute_bf(size_t iel,size_t irad)707 void DFTGridWorker::compute_bf(size_t iel, size_t irad) { 708 // Update function list 709 bf_ind=basp->bf_list(iel); 710 711 // Get radial weights. Only do one radial quadrature point at a 712 // time, since this is an easy way to save a lot of memory. 713 arma::vec wrad(1), r(1); 714 wrad(0)=basp->get_wrad(iel)(irad); 715 r(0)=basp->get_r(iel)(irad); 716 717 double Rhalf(basp->get_Rhalf()); 718 719 // Calculate helpers 720 arma::vec shmu(arma::sinh(r)); 721 722 arma::vec sth(cth.n_elem); 723 for(size_t ia=0;ia<cth.n_elem;ia++) 724 sth(ia)=sqrt(1.0 - cth(ia)*cth(ia)); 725 726 // Radial is 727 scale_r.resize(wrad.n_elem*wang.n_elem); 728 for(size_t ia=0;ia<wang.n_elem;ia++) 729 for(size_t ir=0;ir<wrad.n_elem;ir++) 730 // h_mu = R_{h}\sqrt{\sinh^{2}\mu+\sin^{2}\nu} 731 scale_r(ia*wrad.n_elem+ir)=Rhalf*sqrt(std::pow(shmu(ir),2) + std::pow(sth(ia),2)); 732 // Theta is same as radial 733 scale_theta=scale_r; 734 // phi is simple 735 scale_phi.resize(wrad.n_elem*wang.n_elem); 736 for(size_t ia=0;ia<wang.n_elem;ia++) 737 for(size_t ir=0;ir<wrad.n_elem;ir++) 738 scale_phi(ia*wrad.n_elem+ir)=Rhalf*shmu(ir)*sth(ia); 739 // Update total weights 740 wtot.zeros(wrad.n_elem*wang.n_elem); 741 for(size_t ia=0;ia<wang.n_elem;ia++) 742 for(size_t ir=0;ir<wrad.n_elem;ir++) { 743 size_t idx=ia*wrad.n_elem+ir; 744 // sin(th) is already contained within wang, but we don't want to divide by it since it may be zero. 745 wtot(idx)=wang(ia)*wrad(ir)*std::pow(Rhalf,3)*shmu(ir)*(std::pow(shmu(ir),2)+std::pow(sth(ia),2)); 746 } 747 748 // Compute basis function values 749 bf.zeros(bf_ind.n_elem,wtot.n_elem); 750 // Loop over angular grid 751 #ifdef _OPENMP 752 #pragma omp parallel for 753 #endif 754 for(size_t ia=0;ia<cth.n_elem;ia++) { 755 // Evaluate basis functions at angular point 756 arma::cx_mat abf(basp->eval_bf(iel, irad, cth(ia), phi(ia))); 757 if(abf.n_cols != bf_ind.n_elem) { 758 std::ostringstream oss; 759 oss << "Mismatch! Have " << bf_ind.n_elem << " basis function indices but " << abf.n_cols << " basis functions!\n"; 760 throw std::logic_error(oss.str()); 761 } 762 // Store functions 763 bf.cols(ia*wrad.n_elem,(ia+1)*wrad.n_elem-1)=arma::trans(abf); 764 } 765 766 if(do_grad) { 767 bf_rho.zeros(bf_ind.n_elem,wtot.n_elem); 768 bf_theta.zeros(bf_ind.n_elem,wtot.n_elem); 769 bf_phi.zeros(bf_ind.n_elem,wtot.n_elem); 770 771 #ifdef _OPENMP 772 #pragma omp parallel for 773 #endif 774 for(size_t ia=0;ia<cth.n_elem;ia++) { 775 // Evaluate basis functions at angular point 776 arma::cx_mat dr, dth, dphi; 777 basp->eval_df(iel, irad, cth(ia), phi(ia), dr, dth, dphi); 778 if(dr.n_cols != bf_ind.n_elem) { 779 std::ostringstream oss; 780 oss << "Mismatch! Have " << bf_ind.n_elem << " basis function indices but " << dr.n_cols << " basis functions!\n"; 781 throw std::logic_error(oss.str()); 782 } 783 // Store functions 784 bf_rho.cols(ia*wrad.n_elem,(ia+1)*wrad.n_elem-1)=arma::trans(dr); 785 bf_theta.cols(ia*wrad.n_elem,(ia+1)*wrad.n_elem-1)=arma::trans(dth); 786 bf_phi.cols(ia*wrad.n_elem,(ia+1)*wrad.n_elem-1)=arma::trans(dphi); 787 } 788 } 789 790 if(do_lapl) { 791 throw std::logic_error("Laplacian not implemented.\n"); 792 } 793 } 794 DFTGrid()795 DFTGrid::DFTGrid() { 796 } 797 DFTGrid(const helfem::diatomic::basis::TwoDBasis * basp_,int lang_,int mang_)798 DFTGrid::DFTGrid(const helfem::diatomic::basis::TwoDBasis * basp_, int lang_, int mang_) : basp(basp_), lang(lang_), mang(mang_) { 799 arma::vec cth, phi, wang; 800 helfem::angular::angular_chebyshev(lang,mang,cth,phi,wang); 801 printf("DFT angular grid of order l=%i m=%i has %i points\n",lang,mang,(int) wang.n_elem); 802 } 803 ~DFTGrid()804 DFTGrid::~DFTGrid() { 805 } 806 eval_Fxc(int x_func,const arma::vec & x_pars,int c_func,const arma::vec & c_pars,const arma::mat & P,arma::mat & H,double & Exc,double & Nel,double & Ekin,double thr)807 void DFTGrid::eval_Fxc(int x_func, const arma::vec & x_pars, int c_func, const arma::vec & c_pars, const arma::mat & P, arma::mat & H, double & Exc, double & Nel, double & Ekin, double thr) { 808 H.zeros(basp->Ndummy(),basp->Ndummy()); 809 810 double exc=0.0; 811 double ekin=0.0; 812 double nel=0.0; 813 { 814 DFTGridWorker grid(basp,lang,mang); 815 grid.check_grad_tau_lapl(x_func,c_func); 816 817 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 818 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 819 grid.compute_bf(iel,irad); 820 grid.update_density(P); 821 nel+=grid.compute_Nel(); 822 ekin+=grid.compute_Ekin(); 823 824 grid.init_xc(); 825 if(thr>0.0) 826 grid.screen_density(thr); 827 if(x_func>0) 828 grid.compute_xc(x_func, x_pars); 829 if(c_func>0) 830 grid.compute_xc(c_func, c_pars); 831 832 exc+=grid.eval_Exc(); 833 grid.eval_Fxc(H); 834 835 #if 0 836 std::ostringstream oss; 837 oss << "_" << iel << "_" << irad; 838 grid.save(oss.str()); 839 #endif 840 } 841 } 842 } 843 844 // Save outputs 845 Exc=exc; 846 Ekin=ekin; 847 Nel=nel; 848 849 H=basp->remove_boundaries(H); 850 } 851 eval_Fxc(int x_func,const arma::vec & x_pars,int c_func,const arma::vec & c_pars,const arma::mat & Pa,const arma::mat & Pb,arma::mat & Ha,arma::mat & Hb,double & Exc,double & Nel,double & Ekin,bool beta,double thr)852 void DFTGrid::eval_Fxc(int x_func, const arma::vec & x_pars, int c_func, const arma::vec & c_pars, const arma::mat & Pa, const arma::mat & Pb, arma::mat & Ha, arma::mat & Hb, double & Exc, double & Nel, double & Ekin, bool beta, double thr) { 853 Ha.zeros(basp->Ndummy(),basp->Ndummy()); 854 Hb.zeros(basp->Ndummy(),basp->Ndummy()); 855 856 double exc=0.0; 857 double nel=0.0; 858 double ekin=0.0; 859 { 860 DFTGridWorker grid(basp,lang,mang); 861 grid.check_grad_tau_lapl(x_func,c_func); 862 863 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 864 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 865 grid.compute_bf(iel,irad); 866 grid.update_density(Pa,Pb); 867 nel+=grid.compute_Nel(); 868 ekin+=grid.compute_Ekin(); 869 870 grid.init_xc(); 871 if(thr>0.0) 872 grid.screen_density(thr); 873 if(x_func>0) 874 grid.compute_xc(x_func, x_pars); 875 if(c_func>0) 876 grid.compute_xc(c_func, c_pars); 877 878 exc+=grid.eval_Exc(); 879 grid.eval_Fxc(Ha,Hb,beta); 880 881 #if 0 882 std::ostringstream oss; 883 oss << "_" << iel << "_" << irad; 884 grid.save(oss.str()); 885 #endif 886 } 887 } 888 } 889 890 // Save outputs 891 Exc=exc; 892 Ekin=ekin; 893 Nel=nel; 894 895 // Clean up matrices 896 Ha=basp->remove_boundaries(Ha); 897 Hb=basp->remove_boundaries(Hb); 898 } 899 eval_overlap()900 arma::mat DFTGrid::eval_overlap() { 901 arma::mat S(basp->Ndummy(),basp->Ndummy()); 902 S.zeros(); 903 904 { 905 DFTGridWorker grid(basp,lang,mang); 906 grid.set_grad_tau_lapl(false,false,false); 907 908 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 909 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 910 grid.compute_bf(iel,irad); 911 grid.eval_overlap(S); 912 } 913 } 914 } 915 916 // Clean up matrices 917 return basp->remove_boundaries(S); 918 } 919 eval_kinetic()920 arma::mat DFTGrid::eval_kinetic() { 921 arma::mat T(basp->Ndummy(),basp->Ndummy()); 922 T.zeros(); 923 924 { 925 DFTGridWorker grid(basp,lang,mang); 926 grid.set_grad_tau_lapl(true,false,false); 927 928 for(size_t iel=0;iel<basp->get_rad_Nel();iel++) { 929 for(size_t irad=0;irad<basp->get_r(iel).n_elem;irad++) { 930 grid.compute_bf(iel,irad); 931 grid.eval_kinetic(T); 932 } 933 } 934 } 935 936 // Clean up matrices 937 return basp->remove_boundaries(T); 938 } 939 } 940 } 941 } 942