1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth Weaver
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_ELECINFO_H
22 #define JDFTX_ELECTRONIC_ELECINFO_H
23
24 #include <core/vector3.h>
25 #include <core/MPIUtil.h>
26
27 class matrix;
28 class diagMatrix;
29 class Energies;
30 class Everything;
31
32 //! @addtogroup ElecSystem
33 //! @{
34 //! @file ElecInfo.h
35
36 //! Spin polarization options
37 enum SpinType {
38 SpinNone, //!< unpolarized
39 SpinZ, //!< spin-polarized
40 SpinVector, //!< noncollinear magnetism (supports spin-orbit)
41 SpinOrbit //!< noncollinear but unpolarized (spin-orbit in nonmagnetic systems)
42 };
43
44 //! Quantum number of Bloch state: k-point and spin
45 class QuantumNumber
46 {
47 public:
48 vector3<> k; //!< k-point wave vector
49 int spin; //!< possible spin orientation. up=1, down=-1, none=0
50 double weight; //!< state weight (= 1x or 2x k-point weight depending on spintype)
51
QuantumNumber()52 QuantumNumber() : spin(0) {}
index()53 int index() const { return spin<0 ? 1 : 0; } //!< return the appropriate index into electron (spin) density/potential arrays
54 };
55
56 //! Conversion function needed for PeriodicLookup<QuantumNumber> used for k-point reduction
getCoord(const QuantumNumber & qnum)57 inline vector3<> getCoord(const QuantumNumber& qnum) { return qnum.k; }
58
59 //! Electronic system auxiliary information, fillings and related functions
60 class ElecInfo
61 {
62 public:
63 int nBands, nStates; //!< Number of bands and total number of states
64 int nDensities, spinWeight, qWeightSum; //!< number of density components, spin weight factor (= max occupation per state) and sum of k-point weights
65 int qStart, qStop; //!< Range of states handled by current process (= 0 and nStates for non-MPI jobs)
isMine(int q)66 bool isMine(int q) const { return qDivision.isMine(q); } //!< check if state index is local
whose(int q)67 int whose(int q) const { return qDivision.whose(q); } //!< find out which process this state index belongs to
qStartOther(int iProc)68 int qStartOther(int iProc) const { return qDivision.start(iProc); } //!< find out qStart for another process
qStopOther(int iProc)69 int qStopOther(int iProc) const { return qDivision.stop(iProc); } //!< find out qStop for another process
70
71 SpinType spinType; //!< type of spin treatment
72 double nElectrons; //!< the number of electrons = Sum w Tr[F]
73 std::vector<QuantumNumber> qnums; //!< k-points, spins and weights for each state
74
isNoncollinear()75 bool isNoncollinear() const { return spinType==SpinVector || spinType==SpinOrbit; }
spinorLength()76 int spinorLength() const { return isNoncollinear() ? 2 : 1; }
nSpins()77 int nSpins() const { return spinType==SpinZ ? 2 : 1; }
kFoldingCount()78 vector3<int> kFoldingCount() const { return kfold; }
79
80 //! Fillings update mode
81 enum FillingsUpdate
82 { FillingsConst, //!< constant fillings (T=0)
83 FillingsHsub //!< fillings are a function of subspace Hamiltonian
84 }
85 fillingsUpdate; //!< fillings update mode
86 bool scalarFillings; //!< whether fillings are scalar (equal for all bands) at all quantum numbers
87
88 //! Smearing options
89 enum SmearingType
90 { SmearingFermi, //!< Fermi-Dirac smearing
91 SmearingGauss, //!< Gaussian smearing
92 SmearingMP1, //!< Methfessel-Paxton order 1
93 SmearingCold //!< Cold smearing
94 }
95 smearingType; //!< smearing option
96 double smearingWidth; //!< Smearing width (temperature in the Fermi case)
97 double mu; //!< If NaN, fix nElectrons, otherwise fix/target chemical potential to this
98 double Bz; //!< If NaN, fix magnetization, otherwise fix/target magnetic field to this value
99 bool muLoop; //!< Whether to optimize mu in an outer loop over fixed charge calculations
100
101 bool hasU; //! Flag to check whether the calculation has a DFT+U self-interaction correction
102
103 string initialFillingsFilename; //!< filename for initial fillings (zero-length if none)
104
105 ElecInfo();
106 void setup(const Everything &e, std::vector<diagMatrix>& F, Energies& ener); //!< setup bands and initial fillings
107 void printFillings(FILE* fp) const;
108 void smearReport(const double* muOverride=0) const; //Smearing report (compute mu from eigenvalues in eVars if muOverride not provided)
109 void updateFillingsEnergies(const std::vector<diagMatrix>& eps, Energies&) const; //!< Calculate variable fillings Legendre multipliers (TS/muN)
110
111 //Smearing function utilities:
muEff(double mu,double Bz,int q)112 inline double muEff(double mu, double Bz, int q) const { return mu + Bz*qnums[q].spin; } //!< effective mu for each spin
113 double smear(double mu, double eps) const; //!< smearing function
114 double smearPrime(double mu, double eps) const; //!< derivative of smearing function
115 double smearEntropy(double mu, double eps) const; //!< entropy associated with smearing function
116 diagMatrix smear(double mu, const diagMatrix& eps) const; //!< elementwise smearing function
117 diagMatrix smearPrime(double mu, const diagMatrix& eps) const; //!< elementwise smearing function derivative
118 diagMatrix smearEntropy(double mu, const diagMatrix& eps) const; //!< element-wise smearing function entropy
119
120 //! Propagate matrix gradient w.r.t F to gradient w.r.t. eps (in the basis where fillings are diagonal)
121 matrix smearGrad(double mu, const diagMatrix& eps, const matrix& gradF) const;
122
123 //! Compute number of electrons for the smearing function with specified eigenvalues
124 //! If magnetization is constrained, bisect on corresponding Lagrange multiplier Bz, and retrieve it too
125 double nElectronsCalc(double mu, const std::vector<diagMatrix>& eps, double& Bz) const;
126
127 //! Find the chemical potential for which the smearing function with specified eigenvalues adds up to nElectrons
128 //! If magnetization is constrained, retrieve corresponding Lagrange multiplier Bz as well
129 double findMu(const std::vector<diagMatrix>& eps, double nElectrons, double& Bz) const;
130
131 void kpointsPrint(FILE* fp, bool printSpin=false) const; //!< Output k-points, weights and optionally spins
132 void kpointPrint(FILE* fp, int q, bool printSpin=false) const; //!< Output k-points, weights and optionally spins
133
134 int findHOMO(int q) const; //!< Returns the band index of the Highest Occupied Kohn-Sham Orbital
135
136 //Parallel I/O utilities for diagMatrix/matrix array (one-per-kpoint, with nBands rows and columns unless overridden):
137 void read(std::vector<diagMatrix>&, const char *fname, int nRowsOverride=0) const; //!< parallel read array of diagonal matrices
138 void read(std::vector<matrix>&, const char *fname, int nRowsOverride=0, int nColsOverride=0) const; //!< parallel read array of matrices
139 void write(const std::vector<diagMatrix>&, const char *fname, int nRowsOverride=0) const; //!< parallel write array of diagonal matrices
140 void write(const std::vector<matrix>&, const char *fname, int nRowsOverride=0, int nColsOverride=0) const; //!< parallel write array of matrices
141
142 //Parallel I/O utilities for ColumnBundle array (defined in COlumnBUndle.cpp):
143 struct ColumnBundleReadConversion //!< Utility to convert columnbundle basis / bands
144 { bool realSpace; //!< whether to read realspace wavefunctions
145 int nBandsOld; //!< nBands for the input wavefunction
146 double Ecut, EcutOld; //!< Ecut for the current calculation and input wavefunction in fourier space
147 vector3<int> S_old; //!< fftbox size for the input wavefunction in double space
148 ColumnBundleReadConversion();
149 };
150 void read(std::vector<class ColumnBundle>&, const char *fname, const ColumnBundleReadConversion* conversion=0) const; //!< Read array of columnbundles, optionally with conversion
151 void write(const std::vector<class ColumnBundle>&, const char *fname) const; //!< write an array of columnbundles to file
152
153 private:
154 const Everything* e;
155 TaskDivision qDivision; //!< MPI division of k-points
156
157 //Initial fillings:
158 int nBandsOld; //!<number of bands in file being read
159 double Qinitial, Minitial; //!< net excess electrons and initial magnetization
160 friend struct CommandElecInitialFillings;
161 friend struct CommandElecInitialCharge;
162 friend struct CommandElecInitialMagnetization;
163 friend struct CommandInitialState;
164 friend class ElecVars;
165 friend struct LCAOminimizer;
166 friend void dumpFCI(const Everything& e, const char* filename);
167
168 //!Calculate nElectrons and return magnetization at given mu, Bz and eigenvalues eps
169 double magnetizationCalc(double mu, double Bz, const std::vector<diagMatrix>& eps, double& nElectrons) const;
170
171 //k-points:
172 vector3<int> kfold; //!< kpoint fold vector
173 friend struct CommandKpointFolding;
174 friend class Everything;
175 friend class Phonon;
176 void kpointsFold(); //!< Fold k-points by kfold
177 void kpointsReduce(); //!< Reduce folded k-points under symmetries
178 };
179
180 //! @}
181 #endif // JDFTX_ELECTRONIC_ELECINFO_H
182