1 /*------------------------------------------------------------------- 2 Copyright 2011 Ravishankar Sundararaman 3 Copyright 1996-2003 Sohrab Ismail-Beigi 4 5 This file is part of JDFTx. 6 7 JDFTx is free software: you can redistribute it and/or modify 8 it under the terms of the GNU General Public License as published by 9 the Free Software Foundation, either version 3 of the License, or 10 (at your option) any later version. 11 12 JDFTx is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU General Public License for more details. 16 17 You should have received a copy of the GNU General Public License 18 along with JDFTx. If not, see <http://www.gnu.org/licenses/>. 19 -------------------------------------------------------------------*/ 20 21 #ifndef JDFTX_ELECTRONIC_ELECVARS_H 22 #define JDFTX_ELECTRONIC_ELECVARS_H 23 24 #include <fluid/FluidSolverParams.h> 25 #include <core/ScalarFieldArray.h> 26 #include <string> 27 #include <memory> 28 29 struct ElecGradient; 30 31 //! @addtogroup ElecSystem 32 //! @{ 33 34 //! Electronic variables and main energy calculator 35 class ElecVars 36 { 37 public: 38 std::vector<ColumnBundle> C; //!< orthonormal electronic wavefunctions 39 std::vector<diagMatrix> Haux_eigs; //!< auxilliary hamiltonian eigenvalues 40 std::vector<diagMatrix> F; //!< the fillings (diagonal matrices) for each state 41 42 std::vector<matrix> Hsub; //!< Subspace Hamiltonian: Hsub[q]=C[q]^H*C[q] 43 std::vector<matrix> Hsub_evecs; //!< eigenvectors of Hsub[q] in columns 44 std::vector<diagMatrix> Hsub_eigs; //!< eigenvalues of Hsub[q] 45 46 std::vector< std::vector<matrix> > VdagC; //!< cached pseudopotential projections (by state and then species) 47 48 //Densities and potentials: 49 ScalarFieldArray n; //!< electron density (single ScalarField) or spin density (two ScalarFields [up,dn]) or spin density matrix (four ScalarFields [UpUp, DnDn, Re(UpDn), Im(UpDn)]) 50 ScalarFieldArray nAccum; //!< ElecVars::n accumulated over an MD trajectory 51 ScalarFieldArray get_nXC() const; //!< return the total (spin) density including core contributions get_nTot()52 ScalarField get_nTot() const { return n.size()==1 ? n[0] : n[0]+n[1]; } //!< return the total electron density (even in spin polarized situations) 53 54 ScalarFieldArray tau; //!< kinetic energy density including tauCore, if present (computed if a meta GGA is being used) 55 56 ScalarFieldTilde d_fluid; //!< electrostatic potential due to fluid 57 ScalarFieldTilde V_cavity; //!< non-electrostatic potential on electrons due to fluid 58 59 ScalarFieldArray Vscloc; //! Local part of (optionally spin-dependent) self-consistent potential 60 ScalarFieldArray Vxc; //! Exchange-correlation potential 61 ScalarFieldArray Vtau; //! Gradient w.r.t kinetic energy density (if meta-GGA) 62 63 std::vector<matrix> rhoAtom, U_rhoAtom; //!< Atomic density matrices and gradients w.r.t them (for DFT+U) 64 65 //External interactions: 66 ScalarFieldArray Vexternal; //!< external potential 67 ScalarFieldTilde rhoExternal; //!< external charge density 68 bool rhoExternalSelfEnergy; //!< whether to include self-energy of rhoExternal in output energy 69 struct BoxPotential //!< box potetial desciptor 70 { vector3<> min; //!< minimum of bounding box 71 vector3<> max; //!< maximum of bounding box 72 double Vin; //!< potential inside 73 double Vout; //!< potential outside 74 double convolve_radius; //!< smoothing radius 75 }; 76 std::vector<BoxPotential> boxPot; //!< parameters for the external box potential 77 78 //Fluid properties 79 FluidSolverParams fluidParams; 80 std::shared_ptr<struct FluidSolver> fluidSolver; 81 string fluidInitialStateFilename; 82 83 //Wavefunction initialization: 84 string wfnsFilename; //!< file to read wavefunctions from 85 std::shared_ptr<struct ElecInfo::ColumnBundleReadConversion> readConversion; //!< ColumnBundle conversion 86 bool isRandom; //!< indicates whether the electronic state is random (not yet minimized) 87 bool initLCAO; //!< initialize wave functions using linear combinations of atomic orbitals 88 bool skipWfnsInit; //!< whether to skip wavefunction initialization (used to speed up dry runs, phonon calculations) 89 90 string eigsFilename; //!< file to read eigenvalues from 91 92 //Auxiliary hamiltonian initialization 93 bool HauxInitialized; //!< whether Haux has been read in/computed 94 95 string nFilenamePattern; //!< file pattern to read electron (spin,kinetic) density from 96 string VFilenamePattern; //!< file pattern to read electron (spin,kinetic) potential from 97 98 ElecVars(); 99 void setup(const Everything &everything); 100 101 //! Compute the terms written as a functional of the electronic density, and its gradient i.e. Vscloc 102 //! If supplied, alternateExCorr replaces the main exchange and correlaton functional 103 void EdensityAndVscloc(Energies& ener, const ExCorr* alternateExCorr=0); 104 105 //! Update and return the electronic system energy. 106 //! Optionally compute the gradient, preconditioned gradient and/or the subspace hamiltonian 107 double elecEnergyAndGrad(Energies& ener, ElecGradient* grad=0, ElecGradient* Kgrad=0, bool calc_Hsub=false); 108 109 //! Compute gradient of electronic energy with respect to lattice vectors 110 matrix3<> latticeGrad() const; 111 112 //! Set C to eigenvectors of the subspace hamiltonian 113 void setEigenvectors(); 114 115 //! Compute the kinetic energy density 116 ScalarFieldArray KEdensity() const; 117 118 //! Calculate density using current orthonormal wavefunctions (C) 119 ScalarFieldArray calcDensity() const; 120 121 //! Orthonormalise wavefunctions, with an optional extra rotation 122 //! If extraRotation is present, it is applied after symmetric orthononormalization, 123 //! and on output extraRotation contains the net transformation applied to the wavefunctions. 124 void orthonormalize(int q, matrix* extraRotation=0); 125 126 //! Applies the Kohn-Sham Hamiltonian on the orthonormal wavefunctions C, and computes Hsub if necessary, for a single quantum number 127 //! Returns the Kinetic energy contribution from q, which can be used for the inverse kinetic preconditioner 128 double applyHamiltonian(int q, const diagMatrix& Fq, ColumnBundle& HCq, Energies& ener, bool need_Hsub = false); 129 130 private: 131 const Everything* e; 132 133 std::vector<string> VexternalFilename; //!< external potential filename (read in real space) 134 friend struct CommandVexternal; 135 136 string rhoExternalFilename; //!< external charge filename 137 friend struct CommandRhoExternal; 138 139 int lcaoIter; //!< number of iterations for LCAO (automatic if negative) 140 double lcaoTol; //!< tolerance for LCAO subspace minimization 141 int LCAO(); //!< Initialize LCAO wavefunctions (returns the number of bands initialized) 142 friend struct CommandWavefunction; 143 friend struct CommandLcaoParams; 144 friend class Dump; 145 }; 146 147 //! @} 148 #endif // JDFTX_ELECTRONIC_ELECVARS_H 149