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 SAD_SOLVER_H 18 #define SAD_SOLVER_H 19 20 #include <armadillo> 21 #include "dftgrid.h" 22 #include "../atomic/basis.h" 23 #include "../atomic/dftgrid.h" 24 25 namespace helfem { 26 namespace sadatom { 27 namespace solver { 28 29 /// Helper for printing out configurations 30 typedef struct { 31 int n; 32 int l; 33 double E; 34 int nocc; 35 } shell_occupation_t; 36 /// Sort in energy 37 bool operator<(const shell_occupation_t & lh, const shell_occupation_t & rh); 38 39 /// Helper for defining orbital channels 40 class OrbitalChannel { 41 protected: 42 /// Orbital coefficients: (nbf,nmo,lmax+1) 43 arma::cube C; 44 /// Orbitals energies: (nmo,lmax+1) 45 arma::mat E; 46 /// Orbital occupations per l channel 47 arma::ivec occs; 48 /// Restricted occupations? 49 bool restr; 50 /// Maximum angular channel 51 int lmax; 52 53 /// Get capacity of shell l 54 arma::sword ShellCapacity(arma::sword l) const; 55 /// Count number of occupied orbitals 56 size_t CountOccupied(int l) const; 57 /// Get occupied orbitals 58 std::vector<shell_occupation_t> GetOccupied() const; 59 60 public: 61 /// Dummy constructor 62 OrbitalChannel(); 63 /// Constructor 64 OrbitalChannel(bool restr_); 65 /// Destructor 66 ~OrbitalChannel(); 67 68 /// Is this a doubly occupied configuration? 69 bool Restricted() const; 70 /// Changes restricted setting 71 void SetRestricted(bool restr_); 72 /// Checks if orbitals have been initialized 73 bool OrbitalsInitialized() const; 74 /// Checks if occupations have been initialized 75 bool OccupationsInitialized() const; 76 77 /// Get lmax 78 int Lmax() const; 79 /// Set lmax 80 void SetLmax(int lmax); 81 /// Get coefficients 82 arma::cube Coeffs() const; 83 84 /// Counts the number of electrons 85 arma::sword Nel() const; 86 /// Gives the occupations per l channel 87 arma::ivec Occs() const; 88 /// Sets the occupations 89 void SetOccs(const arma::ivec & occs_); 90 91 /// Get HOMO-LUMO gaps 92 arma::vec GetGap() const; 93 94 /// Characterizes the configuration 95 std::string Characterize() const; 96 /// Print out orbital info 97 void Print(const basis::TwoDBasis & basp) const; 98 /// Save radial part to disk 99 void Save(const basis::TwoDBasis & basis, const std::string & symbol) const; 100 101 /// Checks if the occupations are the same 102 bool operator==(const OrbitalChannel & rh) const; 103 104 /// Updates the orbitals by diagonalization 105 void UpdateOrbitals(const arma::cube & Fl, const arma::mat & Sinvh); 106 /// Updates the orbitals by a damped diagonalization (ov and vo blocks scaled) 107 void UpdateOrbitalsDamped(const arma::cube & Fl, const arma::mat & Sinvh, const arma::mat & S, double dampov); 108 /// Updates the orbitals by diagonalization with a level shift 109 void UpdateOrbitalsShifted(const arma::cube & Fl, const arma::mat & Sinvh, const arma::mat & S, double shift); 110 /// Computes a new density matrix 111 void UpdateDensity(arma::cube & Pl) const; 112 /// Computes a full atomic density matrix 113 arma::mat FullDensity() const; 114 /// Computes an angular density matrix 115 arma::cube AngularDensity() const; 116 117 /// Determines new occupations 118 void AufbauOccupations(arma::sword numel); 119 /// Gives new trial configurations by moving electrons around 120 std::vector<OrbitalChannel> MoveElectrons() const; 121 }; 122 123 /// Restricted configuration 124 typedef struct { 125 /// Orbitals 126 OrbitalChannel orbs; 127 /// Densities 128 arma::cube Pl; 129 /// Fock matrix 130 arma::cube Fl; 131 /// Total energy of configuration 132 double Econf; 133 /// Kinetic energy 134 double Ekin; 135 /// Potential energy 136 double Epot; 137 /// Coulomb repulsion energy 138 double Ecoul; 139 /// Exchange-correlation energy 140 double Exc; 141 /// Converged? 142 bool converged; 143 } rconf_t; 144 /// Checks if orbital occupations match 145 bool operator==(const rconf_t & lh, const rconf_t & rh); 146 /// Orders configurations in energy 147 bool operator<(const rconf_t & lh, const rconf_t & rh); 148 149 typedef struct { 150 /// Orbitals 151 OrbitalChannel orbsa, orbsb; 152 /// Densities 153 arma::cube Pal, Pbl; 154 /// Fock matrices 155 arma::cube Fal, Fbl; 156 /// Total energy of configuration 157 double Econf; 158 /// Kinetic energy 159 double Ekin; 160 /// Potential energy 161 double Epot; 162 /// Coulomb repulsion energy 163 double Ecoul; 164 /// Exchange-correlation energy 165 double Exc; 166 /// Converged? 167 bool converged; 168 } uconf_t; 169 /// Checks if orbital occupations match 170 bool operator==(const uconf_t & lh, const uconf_t & rh); 171 /// Orders configurations in energy 172 bool operator<(const uconf_t & lh, const uconf_t & rh); 173 174 /// SCF Solver 175 class SCFSolver { 176 protected: 177 /// Maximum l value 178 int lmax; 179 180 /// Basis set to use 181 basis::TwoDBasis basis; 182 /// Integration grid 183 dftgrid::DFTGrid grid; 184 185 /// Full atomic basis set for meta-GGAs 186 atomic::basis::TwoDBasis atbasis; 187 /// Full atomic intgeration grid for meta-GGAs 188 atomic::dftgrid::DFTGrid atgrid; 189 190 /// Exchange functional 191 int x_func; 192 /// Exchange functional parameters 193 arma::vec x_pars; 194 /// Correlation functional 195 int c_func; 196 /// Correlation functional parameters 197 arma::vec c_pars; 198 199 /// Overlap matrix 200 arma::mat S; 201 /// Half-inverse overlap 202 arma::mat Sinvh; 203 204 /// Kinetic energy, l-independent part 205 arma::mat T; 206 /// Kinetic energy, l-dependent part 207 arma::mat Tl; 208 /// Nuclear attraction 209 arma::mat Vnuc; 210 /// Core Hamiltonian 211 arma::mat H0; 212 213 /// Maximum number of SCF iterations 214 int maxit; 215 /// Level shift 216 double shift; 217 218 /// SCF convergence threshold (DIIS error) 219 double convthr; 220 /// DFT density threshold 221 double dftthr; 222 /// Start mixing in DIIS when error smaller than 223 double diiseps; 224 /// Use DIIS exclusively when error smaller than 225 double diisthr; 226 /// Number of matrices to keep in memory 227 int diisorder; 228 229 /// Verbose operation? 230 bool verbose; 231 232 /// Form supermatrix 233 arma::mat SuperMat(const arma::mat & M) const; 234 /// Form supermatrix from cube 235 arma::mat SuperCube(const arma::cube & M) const; 236 /// Form cube from supermatrix 237 arma::cube MiniMat(const arma::mat & Msuper) const; 238 /// Form l(l+1)/r^2 cube 239 arma::cube KineticCube() const; 240 /// Replicate matrix into a cube 241 arma::cube ReplicateCube(const arma::mat & M) const; 242 243 public: 244 /// Constructor 245 SCFSolver(int Z, int finitenuc, double Rrms, int lmax, const polynomial_basis::PolynomialBasis * poly, int Nquad, const arma::vec & bval, int x_func_, int c_func_, int maxit_, double shift_, double convthr_, double dftthr_, double diiseps_, double diisthr_, int diisorder_); 246 /// Destructor 247 ~SCFSolver(); 248 249 /// Set functional 250 void set_func(int x_func_, int c_func_); 251 /// Set parameters 252 void set_params(const arma::vec & px, const arma::vec & pc); 253 254 /// Build total density 255 arma::mat TotalDensity(const arma::cube & Pl) const; 256 257 /// Initialize orbitals 258 void Initialize(OrbitalChannel & orbs) const; 259 /// Build the Fock operator, return the energy 260 double FockBuild(rconf_t & conf); 261 /// Build the Fock operator, return the energy 262 double FockBuild(uconf_t & conf); 263 264 /// Solve the SCF problem, return the energy 265 double Solve(rconf_t & conf); 266 /// Solve the SCF problem, return the energy 267 double Solve(uconf_t & conf); 268 269 /// Compute the spin-restricted effective potential 270 arma::mat RestrictedPotential(rconf_t & conf); 271 /// Compute the effective potential as the mean of spin-unrestricted potentials 272 arma::mat UnrestrictedPotential(uconf_t & conf); 273 /// Compute the effective potential as the spin-restricted potential of the average density 274 arma::mat AveragePotential(uconf_t & conf); 275 /// Compute the effective potential as the density weighted average of spin-unrestricted potentials 276 arma::mat WeightedPotential(uconf_t & conf); 277 /// Compute the effective potential for the high-spin case i.e. majority spin 278 arma::mat HighSpinPotential(uconf_t & conf); 279 /// Compute the effective potential for the low-spin case i.e. minority spin 280 arma::mat LowSpinPotential(uconf_t & conf); 281 282 /// Get the basis 283 const basis::TwoDBasis & Basis() const; 284 /// Compute the nuclear density 285 double nuclear_density(const rconf_t & conf) const; 286 /// Compute the nuclear density 287 double nuclear_density(const uconf_t & conf) const; 288 /// Compute the nuclear density gradient 289 double nuclear_density_gradient(const rconf_t & conf) const; 290 /// Compute the nuclear density gradient 291 double nuclear_density_gradient(const uconf_t & conf) const; 292 }; 293 } 294 } 295 } 296 297 #endif 298