1 /* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * HF/DFT from Hel 7 * 8 * Copyright © 2015 The Regents of the University of California 9 * All Rights Reserved 10 * 11 * Written by Susi Lehtola, Lawrence Berkeley National Laboratory 12 * 13 * This program is free software; you can redistribute it and/or 14 * modify it under the terms of the GNU General Public License 15 * as published by the Free Software Foundation; either version 2 16 * of the License, or (at your option) any later version. 17 */ 18 19 #ifndef ERKALE_PZSTAB 20 #define ERKALE_PZSTAB 21 22 #include "scf.h" 23 #include "dftgrid.h" 24 class Timer; 25 26 class FDHessian { 27 protected: 28 /// Verbose operation? 29 bool verbose; 30 /// Finite difference derivative step size 31 double ss_fd; 32 /// Line search step size 33 double ss_ls; 34 35 /// Print optimization status 36 virtual void print_status(size_t iiter, const arma::vec & g, const Timer & t) const; 37 38 public: 39 /// Constructor 40 FDHessian(bool verbose=true); 41 /// Destructor 42 virtual ~FDHessian(); 43 44 /// Get amount of parameters 45 virtual size_t count_params() const=0; 46 /// Evaluate function 47 virtual double eval(const arma::vec & x)=0; 48 /// Update solution 49 virtual void update(const arma::vec & x); 50 51 /// Evaluate finite difference gradient 52 virtual arma::vec gradient(); 53 /// Evaluate finite difference gradient at point x 54 virtual arma::vec gradient(const arma::vec & x); 55 /// Evaluate finite difference Hessian 56 virtual arma::mat hessian(); 57 58 /// Run optimization 59 virtual double optimize(size_t maxiter=1000, double gthr=1e-4, bool max=false); 60 }; 61 62 /// Classify parameters 63 typedef struct { 64 /// Name of the block 65 std::string name; 66 /// Degrees of freedom in block 67 arma::uvec idx; 68 } pz_rot_par_t; 69 70 /// Orbital depedent scaling 71 typedef enum { 72 /// Constant scaling 73 PZ_SCALE_CONSTANT, 74 /// Density based scaling 75 PZ_SCALE_DENSITY, 76 /// Kinetic energy based scaling 77 PZ_SCALE_KINETIC 78 } pz_scaling_t; 79 80 /// PZ optimizer and stability analysis 81 class PZStability: public FDHessian { 82 protected: 83 /// SCF solver, used for energy calculations 84 SCF * solverp; 85 /// Basis set 86 BasisSet basis; 87 /// DFT grid 88 DFTGrid grid; 89 /// NL grid 90 DFTGrid nlgrid; 91 92 /// OV method 93 dft_t ovmethod; 94 /// OO method 95 dft_t oomethod; 96 /// Weight for PZ correction 97 double pzw; 98 /// or scaling method 99 pz_scaling_t scale; 100 /// and scaling exponent 101 double scaleexp; 102 103 /// Reference solution. Spin-restricted 104 rscf_t rsol; 105 /// or unrestricted 106 uscf_t usol; 107 /// Reference self-interaction energies 108 arma::vec ref_Eorb, ref_Eorba, ref_Eorbb; 109 /// Reference orbital Fock matrices 110 std::vector<arma::cx_mat> ref_Forb, ref_Forba, ref_Forbb; 111 /// Reference weighting factors 112 arma::vec ref_worb, ref_worba, ref_worbb; 113 114 /// Real part of transformations? 115 bool real; 116 /// Imaginary part of transformations? 117 bool imag; 118 /// Check stability of canonical orbitals? 119 bool cancheck; 120 /// Check stability of oo block 121 bool oocheck; 122 123 /// Spin-restricted? 124 bool restr; 125 /// Amount of occupied orbitals 126 size_t oa, ob; 127 /// Amount of virtual orbitals 128 size_t va, vb; 129 130 /// Maximum step size 131 double Tmu; 132 133 /// Count amount of parameters for rotations 134 size_t count_ov_params(size_t o, size_t v) const; 135 /// Count amount of parameters for rotations 136 size_t count_oo_params(size_t o) const; 137 /// Count amount of parameters for rotations 138 size_t count_params(size_t o, size_t v) const; 139 /// Count amount of parameters 140 size_t count_params() const; 141 142 /// Classify parameters 143 std::vector<pz_rot_par_t> classify() const; 144 145 /// Calculate rotation matrix 146 arma::cx_mat rotation(const arma::vec & x, bool spin=false) const; 147 /// Form rotation parameter matrix 148 arma::cx_mat rotation_pars(const arma::vec & x, bool spin=false) const; 149 /// Calculate matrix exponential 150 arma::cx_mat matexp(const arma::cx_mat & X) const; 151 152 /// Construct unified Hamiltonian 153 arma::cx_mat unified_H(const arma::cx_mat & CO, const arma::cx_mat & CV, const std::vector<arma::cx_mat> & Forb, const arma::vec & worb, const arma::cx_mat & H0) const; 154 155 /// Evaluate analytic gradient 156 arma::vec gradient(); 157 /// Evaluate analytic gradient at point x 158 arma::vec gradient(const arma::vec & x, bool ref); 159 /// Evaluate semi-analytic Hessian 160 arma::mat hessian(); 161 162 /// Parallel transport 163 void parallel_transport(arma::vec & gold, const arma::vec & sd, double step) const; 164 165 /// Update step size 166 void update_step(const arma::vec & g); 167 /// Perform quasicanonical diagonalisation 168 void diagonalize(); 169 170 /// Get orbital weights 171 arma::vec compute_worb(const arma::cx_mat & C); 172 /// Put in the scaling part of the OO gradient 173 void scaling_gradient_oo(arma::cx_mat & gOO, const arma::cx_mat & CO, const arma::vec & Eorb); 174 /// Put in the scaling part of the OV gradient 175 void scaling_gradient_ov(arma::cx_mat & gOV, const arma::cx_mat & CO, const arma::vec & Eorb, const arma::cx_mat & CV); 176 177 /// Evaluate function at x, and possibly orbital Fock matrices 178 double eval(const arma::vec & x, rscf_t & sol, std::vector<arma::cx_mat> & Forb, arma::vec & Eorb, arma::vec & worb, bool ks, bool fock, bool useref); 179 /// Evaluate function at x, and possibly orbital Fock matrices 180 double eval(const arma::vec & x, uscf_t & sol, std::vector<arma::cx_mat> & Forba, arma::vec & Eorba, arma::vec & worba, std::vector<arma::cx_mat> & Forbb, arma::vec & Eorbb, arma::vec & worbb, bool ks, bool fock, bool useref); 181 /// Evaluate function at x 182 double eval(const arma::vec & x); 183 184 /// Update solution 185 void update(const arma::vec & x); 186 /// Update reference 187 void update_reference(bool sort); 188 /// Update (adaptive) integration grid. If init=true, initialization is done for a static grid 189 void update_grid(bool init); 190 191 /// Print status of optimization 192 void print_status(size_t iiter, const arma::vec & g, const Timer & t) const; 193 194 /// Print information on solution 195 void print_info(const arma::cx_mat & CO, const arma::cx_mat & CV, const std::vector<arma::cx_mat> & Forb, const arma::cx_mat & H0, const arma::vec & Eorb, const arma::vec & worb); 196 197 /// Get the full Fock matrix 198 arma::cx_mat get_H(const rscf_t & sol) const; 199 /// Get the full Fock matrix 200 arma::cx_mat get_H(const uscf_t & sol, bool spin) const; 201 202 /// Precondition gradient vector with unified Hamiltonian 203 arma::vec precondition_unified(const arma::vec & g) const; 204 /// Precondition gradient vector with orbital Hamiltonian 205 arma::vec precondition_orbital(const arma::vec & g) const; 206 207 /// Get occupied orbitals (restricted) 208 arma::cx_mat get_CO(const rscf_t & sol) const; 209 arma::cx_mat get_CO() const; 210 /// Get occupied orbitals (unrestricted) 211 arma::cx_mat get_CO(bool spin, const uscf_t & sol) const; 212 arma::cx_mat get_CO(bool spin) const; 213 /// Get virtual orbitals (restricted) 214 arma::cx_mat get_CV(const rscf_t & sol) const; 215 arma::cx_mat get_CV() const; 216 /// Get virtual orbitals (unrestricted) 217 arma::cx_mat get_CV(bool spin, const uscf_t & sol) const; 218 arma::cx_mat get_CV(bool spin) const; 219 220 public: 221 /// Constructor 222 PZStability(SCF *solver, bool verbose=true); 223 /// Destructor 224 ~PZStability(); 225 226 /// Set method and weight 227 void set_method(const dft_t & ovmethod, const dft_t & oomethod, double pzw, pz_scaling_t scale, double scaleexp); 228 /// Set parameters. real: real rotations? imag: imaginary rotations? ov: ov rotations? oo: oo rotations? 229 void set_params(bool real, bool imag, bool ov, bool oo); 230 231 /// Set reference 232 void set(const rscf_t & sol); 233 /// Set reference 234 void set(const uscf_t & sol); 235 236 /// Evaluate energy 237 double get_E(); 238 239 /// Get updated solution 240 rscf_t get_rsol() const; 241 /// Get updated solution 242 uscf_t get_usol() const; 243 244 /// Check stability of solution. 245 bool check(bool stability=false, double cutoff=-1e-3, double dEthr=-1e-7); 246 /// Print out a line search 247 void linesearch(const std::string & fname="pz_ls.dat", int prec=1, int Np=100); 248 249 /// Print information 250 void print_info(); 251 252 /// Add in a small random perturbation to the solution 253 void perturb(double h=1e-6); 254 255 /// Run optimization 256 virtual double optimize(size_t maxiter=1000, double gthr=1e-4, double nrthr=1e-4, double dEthr=1e-9, int preconditioning=1); 257 }; 258 259 #endif 260