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