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 #ifndef DFTGRID 18 #define DFTGRID 19 20 #include "basis.h" 21 22 namespace helfem { 23 namespace diatomic { 24 namespace dftgrid { 25 26 /// Worker class 27 class DFTGridWorker { 28 protected: 29 /// Basis set 30 const helfem::diatomic::basis::TwoDBasis *basp; 31 32 /// Angular grid 33 arma::vec cth, phi, wang; 34 /// Total quadrature weight 35 arma::rowvec wtot; 36 37 /// Scale factors 38 arma::rowvec scale_r, scale_theta, scale_phi; 39 40 /// List of basis functions in element 41 arma::uvec bf_ind; 42 /// Values of important functions in grid points, Nbf * Ngrid 43 arma::cx_mat bf; 44 /// Radial gradient 45 arma::cx_mat bf_rho; 46 /// Theta gradient 47 arma::cx_mat bf_theta; 48 /// Phi gradient 49 arma::cx_mat bf_phi; 50 /// Values of laplacians in grid points, (3*Nbf) * Ngrid 51 arma::cx_mat bf_lapl; 52 53 /// Density helper matrices: P_{uv} chi_v, and P_{uv} nabla(chi_v) 54 arma::cx_mat Pv, Pv_rho, Pv_theta, Pv_phi; 55 /// Same for spin-polarized 56 arma::cx_mat Pav, Pav_rho, Pav_theta, Pav_phi; 57 arma::cx_mat Pbv, Pbv_rho, Pbv_theta, Pbv_phi; 58 59 /// Is gradient needed? 60 bool do_grad; 61 /// Is kinetic energy density needed? 62 bool do_tau; 63 /// Is laplacian needed? 64 bool do_lapl; 65 66 /// Spin-polarized calculation? 67 bool polarized; 68 69 /// GGA functional used? (Set in compute_xc, only affects eval_Fxc) 70 bool do_gga; 71 /// Meta-GGA tau used? (Set in compute_xc, only affects eval_Fxc) 72 bool do_mgga_t; 73 /// Meta-GGA lapl used? (Set in compute_xc, only affects eval_Fxc) 74 bool do_mgga_l; 75 76 // LDA stuff: 77 78 /// Density, Nrho x Npts 79 arma::mat rho; 80 /// Energy density, Npts 81 arma::rowvec exc; 82 /// Functional derivative of energy wrt electron density, Nrho x Npts 83 arma::mat vxc; 84 85 // GGA stuff 86 87 /// Gradient of electron density, (3 x Nrho) x Npts 88 arma::mat grho; 89 /// Dot products of gradient of electron density, N x Npts; N=1 for closed-shell and 3 for open-shell 90 arma::mat sigma; 91 /// Functional derivative of energy wrt gradient of electron density 92 arma::mat vsigma; 93 94 // Meta-GGA stuff 95 96 /// Laplacian of electron density 97 arma::mat lapl; 98 /// Kinetic energy density 99 arma::mat tau; 100 101 /// Functional derivative of energy wrt laplacian of electron density 102 arma::mat vlapl; 103 /// Functional derivative of energy wrt kinetic energy density 104 arma::mat vtau; 105 106 public: 107 /// Dummy constructor 108 DFTGridWorker(); 109 /// Constructor 110 DFTGridWorker(const helfem::diatomic::basis::TwoDBasis * basp, int lang, int mang); 111 /// Destructor 112 ~DFTGridWorker(); 113 114 /// Check necessity of computing gradient and laplacians, necessary for compute_bf! 115 void check_grad_tau_lapl(int x_func, int c_func); 116 /// Get necessity of computing gradient and laplacians 117 void get_grad_tau_lapl(bool & grad, bool & tau, bool & lapl) const; 118 /// Set necessity of computing gradient and laplacians, necessary for compute_bf! 119 void set_grad_tau_lapl(bool grad, bool tau, bool lapl); 120 121 /// Compute basis functions on grid points 122 void compute_bf(size_t iel, size_t irad); 123 /// Free memory 124 void free(); 125 /// Save data 126 void save(const std::string & info) const; 127 128 /// Update values of density, restricted calculation 129 void update_density(const arma::mat & P); 130 /// Update values of density, unrestricted calculation 131 void update_density(const arma::mat & Pa, const arma::mat & Pb); 132 /// Screen out small densities 133 void screen_density(double thr); 134 135 /// Compute number of electrons 136 double compute_Nel() const; 137 /// Compute kinetic energy 138 double compute_Ekin() const; 139 140 /// Initialize XC arrays 141 void init_xc(); 142 /// Compute XC functional from density and add to total XC 143 /// array. Pot toggles evaluation of potential 144 void compute_xc(int func_id, const arma::vec & params, bool pot=true); 145 /// Evaluate exchange/correlation energy 146 double eval_Exc() const; 147 /// Zero out energy 148 void zero_Exc(); 149 /// Numerical clean up of xc 150 void check_xc(); 151 152 /// Evaluate overlap matrix 153 void eval_overlap(arma::mat & S) const; 154 /// Evaluate kinetic energy matrix 155 void eval_kinetic(arma::mat & T) const; 156 157 /// Evaluate Fock matrix, restricted calculation 158 void eval_Fxc(arma::mat & H) const; 159 /// Evaluate Fock matrix, unrestricted calculation 160 void eval_Fxc(arma::mat & Ha, arma::mat & Hb, bool beta=true) const; 161 }; 162 163 /// Wrapper routine 164 class DFTGrid { 165 private: 166 /// Pointer to basis set 167 const helfem::diatomic::basis::TwoDBasis * basp; 168 /// Angular rule 169 int lang, mang; 170 171 public: 172 /// Dummy constructor 173 DFTGrid(); 174 /// Constructor 175 DFTGrid(const helfem::diatomic::basis::TwoDBasis * basp, int lang, int mang); 176 /// Destructor 177 ~DFTGrid(); 178 179 /// Compute Fock matrix, exchange-correlation energy and integrated electron density, restricted case 180 void 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); 181 /// Compute Fock matrix, exchange-correlation energy and integrated electron density, unrestricted case 182 void 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); 183 184 /// Evaluate overlap 185 arma::mat eval_overlap(); 186 /// Evaluate kinetic energy matrix 187 arma::mat eval_kinetic(); 188 }; 189 190 /// BLAS routine for LDA-type quadrature increment_lda(arma::mat & H,const arma::rowvec & vxc,const arma::Mat<T> & f)191 template<typename T> void increment_lda(arma::mat & H, const arma::rowvec & vxc, const arma::Mat<T> & f) { 192 if(f.n_cols != vxc.n_elem) { 193 std::ostringstream oss; 194 oss << "Number of functions " << f.n_cols << " and potential values " << vxc.n_elem << " do not match!\n"; 195 throw std::runtime_error(oss.str()); 196 } 197 if(H.n_rows != f.n_rows || H.n_cols != f.n_rows) { 198 std::ostringstream oss; 199 oss << "Size of basis function (" << f.n_rows << "," << f.n_cols << ") and Fock matrix (" << H.n_rows << "," << H.n_cols << ") doesn't match!\n"; 200 throw std::runtime_error(oss.str()); 201 } 202 203 // Form helper matrix 204 arma::Mat<T> fhlp(f); 205 for(size_t i=0;i<fhlp.n_rows;i++) 206 for(size_t j=0;j<fhlp.n_cols;j++) 207 fhlp(i,j)*=vxc(j); 208 H+=arma::real(fhlp*arma::trans(f)); 209 } 210 211 /// BLAS routine for GGA-type quadrature increment_gga(arma::mat & H,const arma::mat & gn,const arma::Mat<T> & f,arma::Mat<T> f_x,arma::Mat<T> f_y,arma::Mat<T> f_z)212 template<typename T> void increment_gga(arma::mat & H, const arma::mat & gn, const arma::Mat<T> & f, arma::Mat<T> f_x, arma::Mat<T> f_y, arma::Mat<T> f_z) { 213 if(gn.n_cols!=3) { 214 throw std::runtime_error("Grad rho must have three columns!\n"); 215 } 216 if(f.n_rows != f_x.n_rows || f.n_cols != f_x.n_cols || f.n_rows != f_y.n_rows || f.n_cols != f_y.n_cols || f.n_rows != f_z.n_rows || f.n_cols != f_z.n_cols) { 217 throw std::runtime_error("Sizes of basis function and derivative matrices doesn't match!\n"); 218 } 219 if(H.n_rows != f.n_rows || H.n_cols != f.n_rows) { 220 throw std::runtime_error("Sizes of basis function and Fock matrices doesn't match!\n"); 221 } 222 223 // Compute helper: gamma_{ip} = \sum_c \chi_{ip;c} gr_{p;c} 224 // (N, Np) = (N Np; c) (Np, 3) 225 arma::Mat<T> gamma(f.n_rows,f.n_cols); 226 gamma.zeros(); 227 { 228 // Helper 229 arma::rowvec gc; 230 231 // x gradient 232 gc=arma::strans(gn.col(0)); 233 for(size_t j=0;j<f_x.n_cols;j++) 234 for(size_t i=0;i<f_x.n_rows;i++) 235 f_x(i,j)*=gc(j); 236 gamma+=f_x; 237 238 // x gradient 239 gc=arma::strans(gn.col(1)); 240 for(size_t j=0;j<f_y.n_cols;j++) 241 for(size_t i=0;i<f_y.n_rows;i++) 242 f_y(i,j)*=gc(j); 243 gamma+=f_y; 244 245 // z gradient 246 gc=arma::strans(gn.col(2)); 247 for(size_t j=0;j<f_z.n_cols;j++) 248 for(size_t i=0;i<f_z.n_rows;i++) 249 f_z(i,j)*=gc(j); 250 gamma+=f_z; 251 } 252 253 // Form Fock matrix 254 H+=arma::real(gamma*arma::trans(f) + f*arma::trans(gamma)); 255 } 256 } 257 } 258 } 259 260 #endif 261