MADNESS  0.10.1
dft.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  $Id$
32 */
33 #ifndef DFT_H_
34 #define DFT_H_
35 
36 #include <madness/mra/mra.h>
37 #include <madness/world/MADworld.h>
38 #include <vector>
39 
40 #include "eigsolver.h"
42 
43 class ElectronicStructureAppParams;
44 
45 namespace madness
46 {
47 
48 ////*****************************************************************************
49 //void xc_rks_generic_lda(Tensor<double> rho_alpha, ///< Alpha-spin density at each grid point
50 // Tensor<double> f, ///< Value of functional at each grid point
51 // Tensor<double> df_drho); ///< Derivative of functional w.r.t. rho_alpha
52 ////*****************************************************************************
53 
54  //***************************************************************************
55  template <typename T, int NDIM>
56  class DFTNuclearPotentialOp : public EigSolverOp<T,NDIM>
57  {
58  public:
60  //*************************************************************************
61  // Constructor
62  DFTNuclearPotentialOp(World& world, funcT V, double coeff, double thresh);
63  //*************************************************************************
64 
65  //*************************************************************************
66  // Is there an orbitally-dependent term?
67  virtual bool is_od() {return false;}
68  //*************************************************************************
69 
70  //*************************************************************************
71  // Is there a density-dependent term?
72  virtual bool is_rd() {return true;}
73  //*************************************************************************
74 
75  //*************************************************************************
76  virtual funcT op_r(const funcT& rho, const funcT& psi);
77  //*************************************************************************
78 
79  private:
80  //*************************************************************************
82  //*************************************************************************
83  };
84  //***************************************************************************
85 
86  //***************************************************************************
87  template <typename T, int NDIM>
88  class DFTCoulombOp : public EigSolverOp<T,NDIM>
89  {
90  // Typedef's
92  public:
93  //*************************************************************************
94  // Constructor
95  DFTCoulombOp(World& world, double coeff, double thresh);
96  //*************************************************************************
97 
98  //*************************************************************************
99  // Is there an orbitally-dependent term?
100  virtual bool is_od() {return false;}
101  //*************************************************************************
102 
103  //*************************************************************************
104  // Is there a density-dependent term?
105  virtual bool is_rd() {return true;}
106  //*************************************************************************
107 
108  //*************************************************************************
109  virtual funcT op_r(const funcT& rho, const funcT& psi);
110  //*************************************************************************
111 
112  //*************************************************************************
113  // Build the potential from a density if a density-dependent operator.
114  virtual void prepare_op(Function<double,NDIM> rho);
115  //*************************************************************************
116 
117  //*************************************************************************
119  //*************************************************************************
120 
121  private:
122  //*************************************************************************
124  //*************************************************************************
125 
126  //*************************************************************************
127  bool _spinpol;
128  //*************************************************************************
129  };
130  //***************************************************************************
131 
132  //***************************************************************************
133  template <typename T, int NDIM>
134  class DFTCoulombPeriodicOp : public EigSolverOp<T,NDIM>
135  {
136  // Typedef's
138  public:
139  //*************************************************************************
140  // Constructor
141  DFTCoulombPeriodicOp(World& world, funcT rhon, double coeff, double thresh);
142  //*************************************************************************
143 
144  //*************************************************************************
145  // Is there an orbitally-dependent term?
146  virtual bool is_od() {return false;}
147  //*************************************************************************
148 
149  //*************************************************************************
150  // Is there a density-dependent term?
151  virtual bool is_rd() {return true;}
152  //*************************************************************************
153 
154  //*************************************************************************
155  virtual funcT op_r(const funcT& rho, const funcT& psi);
156  //*************************************************************************
157 
158  //*************************************************************************
159  // Build the potential from a density if a density-dependent operator.
160  virtual void prepare_op(Function<double,NDIM> rho);
161  //*************************************************************************
162 
163  //*************************************************************************
165  //*************************************************************************
166 
167  private:
168  //*************************************************************************
170  //*************************************************************************
171 
172  //*************************************************************************
173  bool _spinpol;
174  //*************************************************************************
175 
176  //*************************************************************************
178  //*************************************************************************
179 
180  };
181  //***************************************************************************
182 
183  //***************************************************************************
184  template <typename T, int NDIM>
185  class XCFunctionalLDA : public EigSolverOp<T,NDIM>
186  {
187  // Typedef's
189  public:
190  //*************************************************************************
191  // Constructor
192  XCFunctionalLDA(World& world, double coeff, double thresh);
193  //*************************************************************************
194 
195  //*************************************************************************
196  // Is there an orbitally-dependent term?
197  virtual bool is_od() {return false;}
198  //*************************************************************************
199 
200  //*************************************************************************
201  // Is there a density-dependent term?
202  virtual bool is_rd() {return true;}
203  //*************************************************************************
204 
205  //*************************************************************************
206  virtual funcT op_r(const funcT& rho, const funcT& psi);
207  //*************************************************************************
208  };
209  //***************************************************************************
210 
211  //***************************************************************************
212  template <typename T, int NDIM>
213  class DFTNuclearChargeDensityOp : public EigSolverOp<T,NDIM>
214  {
215  public:
217  //*************************************************************************
218  // Constructor
219  DFTNuclearChargeDensityOp(World& world, funcT rhon, double coeff,
220  double thresh, bool periodic);
221  //*************************************************************************
222 
223  //*************************************************************************
225  {
226  }
227  //*************************************************************************
228 
229  //*************************************************************************
230  // Is there an orbitally-dependent term?
231  virtual bool is_od() {return false;}
232  //*************************************************************************
233 
234  //*************************************************************************
235  // Is there a density-dependent term?
236  virtual bool is_rd() {return true;}
237  //*************************************************************************
238 
239  //*************************************************************************
241  //*************************************************************************
242 
243  //*************************************************************************
244  virtual funcT op_r(const funcT& rho, const funcT& psi);
245  //*************************************************************************
246 
247  private:
248  //*************************************************************************
250  //*************************************************************************
251 
252  //*************************************************************************
254  //*************************************************************************
255  };
256  //***************************************************************************
257 
258  //***************************************************************************
259  template <typename T, int NDIM>
260  class DFT : public IEigSolverObserver<T,NDIM>
261  {
262  // Typedef's
265  public:
266  //*************************************************************************
267  // Constructor
268  DFT(World& world, funcT vnucrhon, std::vector<funcT> phis,
269  std::vector<double> eigs, ElectronicStructureParams params);
270  //*************************************************************************
271 
272  //*************************************************************************
273  DFT();
274  //*************************************************************************
275 
276  //*************************************************************************
277  virtual ~DFT();
278  //*************************************************************************
279 
280  //*************************************************************************
281  void solve(int maxits);
282  //*************************************************************************
283 
284  //***************************************************************************
285  static double calculate_ke_sp(funcT psi, bool periodic = false);
286  //***************************************************************************
287 
288  //***************************************************************************
289  static double calculate_tot_ke_sp(const std::vector<funcT>& phis,
290  bool spinpol, bool periodic = false);
291  //***************************************************************************
292 
293  //***************************************************************************
294  static double calculate_tot_pe_sp(const World& world,
295  const Function<double,NDIM>& rho, const Function<double,NDIM>& vnucrhon,
296  bool spinpol, const double thresh, bool periodic, bool ispotential);
297  //***************************************************************************
298 
299  //***************************************************************************
300  static double calculate_tot_coulomb_energy(const World& world,
301  const Function<double,NDIM>& rho, bool spinpol, const double thresh,
302  bool periodic = false);
303  //***************************************************************************
304 
305  //***************************************************************************
306  static double calculate_tot_xc_energy(const Function<double,NDIM>& rho);
307  //***************************************************************************
308 
309  //*************************************************************************
310  T matrix_element(const funcT& phii, const funcT& phij)
311  {
312  return _solver->matrix_element(phii, phij);
313  }
314  //*************************************************************************
315 
316  //*************************************************************************
317  void print_matrix_elements(const funcT& phii, const funcT& phij)
318  {
319  _solver->print_matrix_elements(phii, phij);
320  }
321  //*************************************************************************
322 
323  //*************************************************************************
324  virtual void iterateOutput(const std::vector<funcT>& phis,
325  const std::vector<double>& eigs, const Function<double,NDIM>& rho,
326  const int& iter, bool periodic = false);
327  //*************************************************************************
328 
329  //*************************************************************************
330  double get_eig(int indx)
331  {
332  return _solver->get_eig(indx);
333  }
334  //*************************************************************************
335 
336  //*************************************************************************
337  funcT get_phi(int indx)
338  {
339  return _solver->get_phi(indx);
340  }
341  //*************************************************************************
342 
343  //*************************************************************************
344  const std::vector<double>& eigs()
345  {
346  return _solver->eigs();
347  }
348  //*************************************************************************
349 
350  //*************************************************************************
351  const std::vector<funcT>& phis()
352  {
353  return _solver->phis();
354  }
355  //*************************************************************************
356 
357  private:
358 
359  //*************************************************************************
360  // Eigenvalue solver
362  //*************************************************************************
363 
364  //*************************************************************************
366  //*************************************************************************
367 
368  //*************************************************************************
369  // This variable could either be a nuclear potiential or a nuclear charge
370  // density depending on the "ispotential" variable in the
371  // ElectronicStructureParams class.
373  //*************************************************************************
374 
375  //*************************************************************************
376  // Exchange-correlation functional. Needed to compute the energy Exc[rho]
377  // Gets deleted my the EigSolver class during the EigSolver destructor
379  //*************************************************************************
380 
381  //*************************************************************************
383  //*************************************************************************
384 
385  //*************************************************************************
386  World& world() {return _world;}
387  //*************************************************************************
388 
389  };
390  //***************************************************************************
391 
392 }
393 #endif /*DFT_H_*/
This header should include pretty much everything needed for the parallel runtime.
Definition: dft.h:89
Function< T, NDIM > funcT
Definition: dft.h:91
DFTCoulombOp(World &world, double coeff, double thresh)
Definition: dft.cc:84
funcT _Vc
Definition: dft.h:123
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: dft.cc:153
SeparatedConvolution< T, NDIM > * _cop
Definition: dft.h:118
virtual void prepare_op(Function< double, NDIM > rho)
Definition: dft.cc:119
bool _spinpol
Definition: dft.h:127
virtual bool is_od()
Is there an orbitally-dependent term?
Definition: dft.h:100
virtual bool is_rd()
Is there a density-dependent term?
Definition: dft.h:105
Definition: dft.h:135
Function< T, NDIM > funcT
Definition: dft.h:137
SeparatedConvolution< T, NDIM > * _cop
Definition: dft.h:164
virtual void prepare_op(Function< double, NDIM > rho)
Definition: dft.cc:127
DFTCoulombPeriodicOp(World &world, funcT rhon, double coeff, double thresh)
Definition: dft.cc:101
bool _spinpol
Definition: dft.h:173
funcT _rhon
Definition: dft.h:177
virtual bool is_rd()
Is there a density-dependent term?
Definition: dft.h:151
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: dft.cc:162
funcT _Vc
Definition: dft.h:169
virtual bool is_od()
Is there an orbitally-dependent term?
Definition: dft.h:146
funcT _rhon
Definition: dft.h:249
Function< T, NDIM > funcT
Definition: dft.h:216
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: dft.cc:135
~DFTNuclearChargeDensityOp()
Definition: dft.h:224
DFTNuclearChargeDensityOp(World &world, funcT rhon, double coeff, double thresh, bool periodic)
Definition: dft.cc:46
void prepare_op(Function< double, NDIM > rho)
Definition: dft.h:240
funcT _Vnuc
Definition: dft.h:253
virtual bool is_od()
Is there an orbitally-dependent term?
Definition: dft.h:231
virtual bool is_rd()
Is there a density-dependent term?
Definition: dft.h:236
Definition: dft.h:57
Function< T, NDIM > funcT
Definition: dft.h:59
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: dft.cc:144
virtual bool is_rd()
Is there a density-dependent term?
Definition: dft.h:72
funcT _V
Definition: dft.h:81
virtual bool is_od()
Is there an orbitally-dependent term?
Definition: dft.h:67
DFTNuclearPotentialOp(World &world, funcT V, double coeff, double thresh)
Definition: dft.cc:73
Definition: dft.h:261
static double calculate_tot_coulomb_energy(const World &world, const Function< double, NDIM > &rho, bool spinpol, const double thresh, bool periodic=false)
Definition: dft.cc:325
double get_eig(int indx)
Definition: dft.h:330
static double calculate_tot_ke_sp(const std::vector< funcT > &phis, bool spinpol, bool periodic=false)
Definition: dft.cc:276
EigSolverOp< T, NDIM > * _xcfunc
Definition: dft.h:378
Function< double, NDIM > _vnucrhon
Definition: dft.h:372
Vector< double, NDIM > kvecT
Definition: dft.h:264
void print_matrix_elements(const funcT &phii, const funcT &phij)
Definition: dft.h:317
static double calculate_tot_pe_sp(const World &world, const Function< double, NDIM > &rho, const Function< double, NDIM > &vnucrhon, bool spinpol, const double thresh, bool periodic, bool ispotential)
Definition: dft.cc:294
funcT get_phi(int indx)
Definition: dft.h:337
static double calculate_tot_xc_energy(const Function< double, NDIM > &rho)
Definition: dft.cc:352
EigSolver< T, NDIM > * _solver
Definition: dft.h:361
const std::vector< double > & eigs()
Definition: dft.h:344
World & _world
Definition: dft.h:365
virtual ~DFT()
Definition: dft.cc:245
virtual void iterateOutput(const std::vector< funcT > &phis, const std::vector< double > &eigs, const Function< double, NDIM > &rho, const int &iter, bool periodic=false)
Definition: dft.cc:363
ElectronicStructureParams _params
Definition: dft.h:382
T matrix_element(const funcT &phii, const funcT &phij)
Definition: dft.h:310
const std::vector< funcT > & phis()
Definition: dft.h:351
Function< T, NDIM > funcT
Definition: dft.h:263
static double calculate_ke_sp(funcT psi, bool periodic=false)
Definition: dft.cc:261
World & world()
Definition: dft.h:386
void solve(int maxits)
Definition: dft.cc:253
Definition: eigsolver.h:99
double thresh()
Definition: eigsolver.h:193
double coeff()
Definition: eigsolver.h:177
Definition: eigsolver.h:226
A multiresolution adaptive numerical function.
Definition: mra.h:122
Definition: eigsolver.h:52
A parallel world class.
Definition: world.h:132
Definition: dft.h:186
virtual funcT op_r(const funcT &rho, const funcT &psi)
Density-dependent portion of operator.
Definition: dft.cc:181
XCFunctionalLDA(World &world, double coeff, double thresh)
Definition: dft.cc:171
virtual bool is_rd()
Is there a density-dependent term?
Definition: dft.h:202
virtual bool is_od()
Is there an orbitally-dependent term?
Definition: dft.h:197
Function< T, NDIM > funcT
Definition: dft.h:188
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
static const double thresh
Definition: rk.cc:45
Definition: electronicstructureparams.h:46
static double V(const coordT &r)
Definition: tdse.cc:288