MADNESS  0.10.1
eigsolver.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 */
31 
32 #ifndef EIGSOLVER_H_
33 #define EIGSOLVER_H_
34 #include <madness/mra/mra.h>
35 #include <madness/world/MADworld.h>
36 #include <vector>
38 
39 namespace madness
40 {
41 //***************************************************************************
42 /// This is the interface the an observer wishing to receive output must
43 /// implement. This call back gives the current eigenfunctions, eigenvalues,
44 /// and the density.
45 /// This is a test LaTeX formula
46 /// The Pythagorean theorem is
47 /// \f[
48 /// c^2 = a^2 + b^2
49 /// \f]
50 template <typename T, int NDIM>
52 {
54 public:
55  virtual void iterateOutput(const std::vector<funcT>& phis,
56  const std::vector<double>& eigs, const Function<double, NDIM>& rho, const int& iter, bool periodic) = 0;
57 
58  virtual ~IEigSolverObserver() {};
59 };
60 //***************************************************************************
61 
62 class KPoint
63 {
64 public:
65  KPoint(double kx, double ky, double kz, double weight)
66  {
67  _kx = kx; _ky = ky; _kz = kz;
68  _weight = weight;
69  }
70 
71  //*************************************************************************
72  double kx() {return _kx;}
73  double ky() {return _ky;}
74  double kz() {return _kz;}
75  //*************************************************************************
76 
77  //*************************************************************************
78  double weight() {return _weight;}
79  //*************************************************************************
80 
81 private:
82  //*************************************************************************
83  // the actual k-point
84  double _kx;
85  double _ky;
86  double _kz;
87  //*************************************************************************
88 
89  //*************************************************************************
90  // weight
91  double _weight;
92  //*************************************************************************
93 
94 };
95 
96 //***************************************************************************
97 template <typename T, int NDIM>
99 {
100  // Typedef's
102 public:
103  //*************************************************************************
104  // Constructor
105  EigSolverOp(World& world, double coeff, double thresh)
106  : _world(world), _coeff(coeff), _thresh(thresh) {}
107  //*************************************************************************
108 
109  //*************************************************************************
110  // Destructor
111  virtual ~EigSolverOp() {}
112  //*************************************************************************
113 
114  //*************************************************************************
115  /// Is there an orbitally-dependent term?
116  virtual bool is_od() = 0;
117  //*************************************************************************
118 
119  //*************************************************************************
120  /// Is there a density-dependent term?
121  virtual bool is_rd() = 0;
122  //*************************************************************************
123 
124  //*************************************************************************
125  /// Build the potential from a density if a density-dependent operator.
126  virtual void prepare_op(funcT rho) {}
127  //*************************************************************************
128 
129  //*************************************************************************
130  /// Orbital-dependent portion of operator
131  virtual funcT op_o(const std::vector<funcT>& phis, const funcT& psi)
132  {
134  return func;
135  }
136  //*************************************************************************
137 
138  //*************************************************************************
139  /// Density-dependent portion of operator
140  virtual funcT op_r(const funcT& rho, const funcT& psi)
141  {
143  return func;
144  }
145  //*************************************************************************
146 
147  //*************************************************************************
148  /// Orbital-dependent portion of operator
149  virtual std::vector<funcT> multi_op_o(const std::vector<funcT>& phis)
150  {
151  // Collection of empty functions
152  std::vector<funcT> newphis(phis.size(), FunctionFactory<T,NDIM>(_world));
153  for (unsigned int pi = 0; pi < phis.size(); pi++)
154  {
155  newphis[pi] = op_o(phis, phis[pi]);
156  }
157  _world.gop.fence();
158  return newphis;
159  }
160  //*************************************************************************
161 
162  //*************************************************************************
163  /// Density-dependent portion of operator
164  virtual std::vector<funcT> multi_op_r(const funcT& rho, const std::vector<funcT>& phis)
165  {
166  std::vector<funcT> newphis(phis.size(), FunctionFactory<T,NDIM>(_world));
167  for (unsigned int pi = 0; pi < phis.size(); pi++)
168  {
169  newphis[pi] = op_r(rho, phis[pi]);
170  }
171  _world.gop.fence();
172  return newphis;
173  }
174  //*************************************************************************
175 
176  //*************************************************************************
177  double coeff() {return _coeff;}
178  //*************************************************************************
179 
180  //*************************************************************************
181  std::string messsageME()
182  {
183  return _messageME;
184  }
185  //*************************************************************************
186 
187  //*************************************************************************
189  //*************************************************************************
190 
191 protected:
192  //*************************************************************************
193  double thresh() {return _thresh;}
194  //*************************************************************************
195 
196  //*************************************************************************
197  void messageME(std::string messageME)
198  {
200  }
201  //*************************************************************************
202 
203 private:
204  //*************************************************************************
205  double _coeff;
206  //*************************************************************************
207 
208  //*************************************************************************
209  double _thresh;
210  //*************************************************************************
211 
212  //*************************************************************************
213  std::string _messageME;
214  //*************************************************************************
215 
216 };
217 //***************************************************************************
218 
219 //***************************************************************************
220 /// The EigSolver class is the class that is the workhorse of both the Hartree
221 /// Fock and the DFT algorithms. This class relies on the wrapper class to
222 /// give it a list of operators to implement as its potential. This should
223 /// allow for much more reuse.
224 template <typename T, int NDIM>
226 {
227 public:
228  //*************************************************************************
229  // Typedef's
231 // typedef KPoint<NDIM> kvecT;
234  typedef std::shared_ptr<operatorT> poperatorT;
235  //*************************************************************************
236 
237  //*************************************************************************
238  /// Constructor for periodic system
239  EigSolver(World& world, std::vector<funcT> phis, std::vector<double> eigs,
240  std::vector<EigSolverOp<T,NDIM>*> ops, std::vector<kvecT> kpoints,
242  //*************************************************************************
243 
244  //*************************************************************************
245  /// Constructor for non-periodic system
246  EigSolver(World& world, std::vector<funcT> phis, std::vector<double> eigs,
247  std::vector<EigSolverOp<T,NDIM>*> ops, ElectronicStructureParams params);
248  //*************************************************************************
249 
250  //*************************************************************************
251  /// Destructor
252  virtual ~EigSolver();
253  //*************************************************************************
254 
255  //*************************************************************************
256  /// This solver has not been optimized for usage in parallel. This solver
257  /// processes each eigenfunction in a serial fashion.
258  void solve(int maxits);
259  //*************************************************************************
260 
261  //*************************************************************************
262  /// This solver has been optimized for usage in parallel. This solver
263  /// processes each eigenfunction in a parallel fashion.
264  void multi_solve(int maxits);
265  //*************************************************************************
266 
267  //*************************************************************************
268  double get_eig(int indx)
269  {
270  return _eigs[indx];
271  }
272  //*************************************************************************
273 
274  //*************************************************************************
275  funcT get_phi(int indx)
276  {
277  return _phis[indx];
278  }
279  //*************************************************************************
280 
281  //*************************************************************************
282  const std::vector<funcT>& phis()
283  {
284  return _phis;
285  }
286  //*************************************************************************
287 
288  //*************************************************************************
289  const std::vector<double>& eigs()
290  {
291  return _eigs;
292  }
293  //*************************************************************************
294 
295  //*************************************************************************
297  {
298  _obs.push_back(obs);
299  }
300  //*************************************************************************
301 
302  //*************************************************************************
303  /// Computes a matrix element given the left and right functions.
304  T matrix_element(const funcT& phii, const funcT& phij);
305  //*************************************************************************
306 
307  //*************************************************************************
308  /// Prints a matrix element given the left and right functions.
309  void print_matrix_elements(const funcT& phii, const funcT& phij);
310  //*************************************************************************
311 
312  //*************************************************************************
313  /// Preprocesses the operators for doing an iteration of "eigensolving".
314  void prepare_ops();
315  //*************************************************************************
316 
317  //*************************************************************************
318  /// Makes the BSH Green's functions for the parallel solver (multi_solve()).
319  void make_bsh_operators();
320  //*************************************************************************
321 
322  //*************************************************************************
323  void update_occupation();
324  //*************************************************************************
325 
326  //*************************************************************************
327  /// Computes the electronic density
328  static funcT compute_rho(std::vector<funcT> phis, std::vector<double> occs,
329  const World& world);
330  //*************************************************************************
331 
332 private:
333  //*************************************************************************
334  /// List of the functions
335  std::vector<funcT> _phis;
336  //*************************************************************************
337 
338  //*************************************************************************
339  /// List of the eigenvalues
340  std::vector<double> _eigs;
341  //*************************************************************************
342 
343  //*************************************************************************
344  /// List of the ops
345  std::vector< EigSolverOp<T,NDIM>* > _ops;
346  //*************************************************************************
347 
348  //*************************************************************************
349  /// List of the ops
350  std::vector<kvecT> _kpoints;
351  //*************************************************************************
352 
353  //*************************************************************************
355  //*************************************************************************
356 
357  //*************************************************************************
358  // List of the obs
359  std::vector<IEigSolverObserver<T,NDIM>*> _obs;
360  //*************************************************************************
361 
362  // Electronic charge density
363  //*************************************************************************
365  //*************************************************************************
366 
367  //*************************************************************************
368  // List of the ops
369  std::vector<poperatorT> _bops;
370  //*************************************************************************
371 
372  //*************************************************************************
373  // List of the occupation numbers
374  std::vector<double> _occs;
375  //*************************************************************************
376 
377  //*************************************************************************
379  //*************************************************************************
380 
381 };
382 //***************************************************************************
383 
384 }
385 
386 #endif /*EIGSOLVER_H_*/
387 
This header should include pretty much everything needed for the parallel runtime.
Definition: eigsolver.h:99
virtual std::vector< funcT > multi_op_o(const std::vector< funcT > &phis)
Orbital-dependent portion of operator.
Definition: eigsolver.h:149
World & _world
Definition: eigsolver.h:188
virtual funcT op_o(const std::vector< funcT > &phis, const funcT &psi)
Orbital-dependent portion of operator.
Definition: eigsolver.h:131
virtual std::vector< funcT > multi_op_r(const funcT &rho, const std::vector< funcT > &phis)
Density-dependent portion of operator.
Definition: eigsolver.h:164
std::string messsageME()
Definition: eigsolver.h:181
virtual bool is_rd()=0
Is there a density-dependent term?
double _thresh
Definition: eigsolver.h:209
double thresh()
Definition: eigsolver.h:193
std::string _messageME
Definition: eigsolver.h:213
double _coeff
Definition: eigsolver.h:205
virtual void prepare_op(funcT rho)
Build the potential from a density if a density-dependent operator.
Definition: eigsolver.h:126
EigSolverOp(World &world, double coeff, double thresh)
Definition: eigsolver.h:105
virtual ~EigSolverOp()
Definition: eigsolver.h:111
Function< T, NDIM > funcT
Definition: eigsolver.h:101
void messageME(std::string messageME)
Definition: eigsolver.h:197
double coeff()
Definition: eigsolver.h:177
virtual bool is_od()=0
Is there an orbitally-dependent term?
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: eigsolver.h:140
Definition: eigsolver.h:226
World & _world
Definition: eigsolver.h:354
std::vector< double > _occs
Definition: eigsolver.h:374
void print_matrix_elements(const funcT &phii, const funcT &phij)
Prints a matrix element given the left and right functions.
Definition: eigsolver.cc:242
void make_bsh_operators()
Makes the BSH Green's functions for the parallel solver (multi_solve()).
Definition: eigsolver.cc:214
void prepare_ops()
Preprocesses the operators for doing an iteration of "eigensolving".
Definition: eigsolver.cc:128
std::vector< double > _eigs
List of the eigenvalues.
Definition: eigsolver.h:340
std::vector< IEigSolverObserver< T, NDIM > * > _obs
Definition: eigsolver.h:359
std::vector< funcT > _phis
List of the functions.
Definition: eigsolver.h:335
std::vector< kvecT > _kpoints
List of the ops.
Definition: eigsolver.h:350
static funcT compute_rho(std::vector< funcT > phis, std::vector< double > occs, const World &world)
Computes the electronic density.
Definition: eigsolver.cc:107
Function< T, NDIM > funcT
Definition: eigsolver.h:230
ElectronicStructureParams _params
Definition: eigsolver.h:378
std::vector< poperatorT > _bops
Definition: eigsolver.h:369
void multi_solve(int maxits)
Definition: eigsolver.cc:413
std::vector< EigSolverOp< T, NDIM > * > _ops
List of the ops.
Definition: eigsolver.h:345
T matrix_element(const funcT &phii, const funcT &phij)
Computes a matrix element given the left and right functions.
Definition: eigsolver.cc:143
Vector< double, NDIM > kvecT
Definition: eigsolver.h:232
void update_occupation()
Definition: eigsolver.cc:168
EigSolver(World &world, std::vector< funcT > phis, std::vector< double > eigs, std::vector< EigSolverOp< T, NDIM > * > ops, std::vector< kvecT > kpoints, ElectronicStructureParams params)
Constructor for periodic system.
Definition: eigsolver.cc:55
double get_eig(int indx)
Definition: eigsolver.h:268
void solve(int maxits)
Definition: eigsolver.cc:279
funcT get_phi(int indx)
Definition: eigsolver.h:275
void addObserver(IEigSolverObserver< T, NDIM > *obs)
Definition: eigsolver.h:296
Function< double, NDIM > _rho
Definition: eigsolver.h:364
SeparatedConvolution< double, NDIM > operatorT
Definition: eigsolver.h:233
const std::vector< funcT > & phis()
Definition: eigsolver.h:282
virtual ~EigSolver()
Destructor.
Definition: eigsolver.cc:88
const std::vector< double > & eigs()
Definition: eigsolver.h:289
std::shared_ptr< operatorT> poperatorT
Definition: eigsolver.h:234
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
A multiresolution adaptive numerical function.
Definition: mra.h:122
Definition: eigsolver.h:52
virtual ~IEigSolverObserver()
Definition: eigsolver.h:58
Function< T, NDIM > funcT
Definition: eigsolver.h:53
virtual void iterateOutput(const std::vector< funcT > &phis, const std::vector< double > &eigs, const Function< double, NDIM > &rho, const int &iter, bool periodic)=0
Definition: eigsolver.h:63
double kz()
Definition: eigsolver.h:74
double ky()
Definition: eigsolver.h:73
double _kz
Definition: eigsolver.h:86
double _kx
Definition: eigsolver.h:84
double weight()
Definition: eigsolver.h:78
double _weight
Definition: eigsolver.h:91
double kx()
Definition: eigsolver.h:72
KPoint(double kx, double ky, double kz, double weight)
Definition: eigsolver.h:65
double _ky
Definition: eigsolver.h:85
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
A parallel world class.
Definition: world.h:132
WorldGopInterface & gop
Global operations.
Definition: world.h:205
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
Main include file for MADNESS and defines Function interface.
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
Definition: electronicstructureparams.h:46
static const double pi
Definition: testcosine.cc:6