Psi4
cubature.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2017 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This file is part of Psi4.
12  *
13  * Psi4 is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser General Public License as published by
15  * the Free Software Foundation, version 3.
16  *
17  * Psi4 is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License along
23  * with Psi4; if not, write to the Free Software Foundation, Inc.,
24  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25  *
26  * @END LICENSE
27  */
28 
29 #ifndef libmints_cubature_H
30 #define libmints_cubature_H
31 
32 #include "psi4/psi4-dec.h"
33 
34 #include "psi4/libmints/vector3.h"
35 #include "psi4/libmints/typedefs.h"
36 
37 #include <map>
38 #include <vector>
39 
40 namespace psi {
41 
42 class BasisSet;
43 class Matrix;
44 class Vector;
45 class BasisExtents;
46 class BlockOPoints;
47 class RadialGrid;
48 class SphericalGrid;
49 class Options;
50 
51 // This is an auxiliary structure used internally by the grid-builder class.
52 // Apparently, for performance reasons, it is not good for the final molecular grid
53 // to consist of structures like this one.
54 //
55 // RMP: You got that right. Calling nuclear weights one point at a time is another
56 // great way to get a 1000x slowdown. What an incredible smell you've discovered!
57 //
58 struct MassPoint {
59  double x,y,z,w;
60 };
61 
62 
64 protected:
65  int debug_;
66 
68  std::shared_ptr<Molecule> molecule_;
69 
70  // ==> Fast Grid Specification <== //
71 
73  int npoints_;
79  double* x_;
81  double* y_;
83  double* z_;
85  double* w_;
86 
87  // ==> Clean Grid Specification <== //
88 
90  std::shared_ptr<Matrix> orientation_;
92  std::vector<std::shared_ptr<RadialGrid> > radial_grids_;
94  std::vector<std::vector<std::shared_ptr<SphericalGrid> > > spherical_grids_;
96  int* index_;
97 
99  std::vector<std::shared_ptr<BlockOPoints> > blocks_;
100 
102  std::shared_ptr<BasisExtents> extents_;
104  std::shared_ptr<BasisSet> primary_;
105 
107  void postProcess(std::shared_ptr<BasisExtents> extents, int max_points, int min_points, double max_radius);
108  void remove_distant_points(double Rcut);
109  void block(int max_points, int min_points, double max_radius);
110 
111 public:
115  short radscheme; // Effectively an enumeration
116  short prunescheme;
117  short nucscheme;
118  short namedGrid; // -1 = None, 0 = SG-0, 1 = SG-1
119  int nradpts;
120  int nangpts;
121  };
122 protected:
125 public:
126  MolecularGrid(std::shared_ptr<Molecule> molecule);
127  virtual ~MolecularGrid();
128 
133  const std::vector<std::vector<double> >& rs, // Radial nodes, per atom
134  const std::vector<std::vector<double> >& ws, // Radial weights, per atom
135  const std::vector<std::vector<int> >& Ls); // Spherical orders, per atom
136 
138  void print(std::string out_fname = "outfile", int print = 2) const;
139  void print_details(std::string out_fname = "outfile", int print = 2) const;
140 
142  std::shared_ptr<Matrix> orientation() const { return orientation_; }
144  const std::vector<std::shared_ptr<RadialGrid> >& radial_grids() const { return radial_grids_; }
146  const std::vector<std::vector<std::shared_ptr<SphericalGrid> > >& spherical_grids() const { return spherical_grids_; }
148  int* index() const { return index_; }
149 
151  int npoints() const { return npoints_; }
153  int max_points() const { return max_points_; }
155  int max_functions() const { return max_functions_; }
156 
158  double* x() const { return x_; }
160  double* y() const { return y_; }
162  double* z() const { return z_; }
164  double* w() const { return w_; }
165 
167  std::shared_ptr<BasisExtents> extents() const { return extents_; }
169  const std::vector<std::shared_ptr<BlockOPoints> >& blocks() const { return blocks_; }
170 
171  void set_debug(int debug) { debug_ = debug; }
172 };
173 
175 
176 protected:
178  std::shared_ptr<BasisSet> primary_;
180  std::string filename_;
181 
184 
186  void buildGridFromOptions();
187 
188 public:
189 
191  PseudospectralGrid(std::shared_ptr<Molecule> molecule,
192  std::shared_ptr<BasisSet> primary,
193  Options& options);
194  virtual ~PseudospectralGrid();
195 
196 };
197 
198 class DFTGrid : public MolecularGrid {
199 
200 protected:
202  std::shared_ptr<BasisSet> primary_;
204  void buildGridFromOptions(std::map<std::string, int> int_opts_map,
205  std::map<std::string, std::string> opts_map);
208 
209 public:
210  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
211  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary,
212  std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map,
213  Options& options);
214 virtual ~DFTGrid();
215 };
216 
217 class RadialGrid {
218 
219 protected:
220 
222  std::string scheme_;
224  int npoints_;
226  double alpha_;
228  double* r_;
230  double* w_;
231 
232  // ==> Standard Radial Grids <== //
233 
235  static std::shared_ptr<RadialGrid> build_becke(int npoints, double alpha);
237  static std::shared_ptr<RadialGrid> build_treutler(int npoints, double alpha);
238  // TODO: Add more grids
239 
241  RadialGrid();
242 public:
243  // ==> Initializers <== //
244 
246  virtual ~RadialGrid();
247 
249  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double alpha);
251  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double* r, double* wr, double alpha);
252 
253  // ==> Accessors <== //
254 
256  std::string scheme() const { return scheme_; }
258  int npoints() const { return npoints_; }
260  double alpha() const { return alpha_; }
262  double* r() const { return r_; }
264  double* w() const { return w_; }
265 
267  void print(std::string out_fname = "outfile", int level = 1) const;
268 };
269 
271 
272 protected:
273 
275  std::string scheme_;
277  int npoints_;
279  int order_;
281  double* x_;
283  double* y_;
285  double* z_;
287  double* w_;
288 
290  double* phi_;
292  double* theta_;
293 
294  // ==> Unique Lebedev Grids (statically stored) <== //
295 
297  static std::map<int, int> lebedev_mapping_;
299  static void initialize_lebedev();
301  static void lebedev_error();
302 
303  // ==> Utility Routines <== //
304 
306  void build_angles();
307 
309  SphericalGrid();
310 public:
311  // ==> Initializers <== //
312 
314  virtual ~SphericalGrid();
315 
317  static std::shared_ptr<SphericalGrid> build(const std::string& scheme, int npoints, const MassPoint* points);
318 
319  // ==> Accessors <== //
320 
322  std::string scheme() const { return scheme_; }
324  int npoints() const { return npoints_; }
326  int order() const { return order_; }
328  double* x() const { return x_; }
330  double* y() const { return y_; }
332  double* z() const { return z_; }
334  double* w() const { return w_; }
335 
337  double* phi() const { return phi_; }
339  double* theta() const { return theta_; }
340 
342  void print(std::string out_fname = "outfile", int level = 1) const;
343 
344 
345 };
346 
348 
349 protected:
351  int npoints_;
352 
358 
360  double* x_;
362  double* y_;
364  double* z_;
366  double* w_;
368  std::vector<int> shells_local_to_global_;
370  std::vector<int> functions_local_to_global_;
372  std::shared_ptr<BasisExtents> extents_;
373 
377  double R_;
378 
380  void populate();
382  void bound();
383 
384 public:
386  std::shared_ptr<BasisExtents> extents);
387  BlockOPoints(int npoints, double* x, double* y, double* z, double* w,
388  std::shared_ptr<BasisExtents> extents);
389  virtual ~BlockOPoints();
390 
392  void refresh() { populate(); }
393 
395  int npoints() const { return npoints_; }
397  void print(std::string out_fname = "outfile", int print = 2);
398 
400  double* x() const { return x_; }
402  double* y() const { return y_; }
404  double* z() const { return z_; }
406  double* w() const { return w_; }
407 
409  const std::vector<int>& shells_local_to_global() const { return shells_local_to_global_; }
411  const std::vector<int>& functions_local_to_global() const { return functions_local_to_global_; }
412 };
413 
415 
416 protected:
418  std::shared_ptr<BasisSet> primary_;
420  double delta_;
422  std::shared_ptr<Vector> shell_extents_;
424  double maxR_;
425 
427  void computeExtents();
428 public:
429  BasisExtents(std::shared_ptr<BasisSet> primary, double delta);
430  virtual ~BasisExtents();
431 
433  void print(std::string out_fname = "outfile");
435  void set_delta(double delta) { delta_ = delta; computeExtents(); }
436 
438  double delta() const { return delta_; }
440  std::shared_ptr<BasisSet> basis() const { return primary_; }
442  std::shared_ptr<Vector> shell_extents() const { return shell_extents_; }
444  double maxR() const { return maxR_; }
445 };
446 
447 }
448 #endif
Options & options_
The Options object.
Definition: cubature.h:183
void print(std::string out_fname="outfile")
Print a trace of these extents.
Definition: cubature.cc:3826
virtual ~BlockOPoints()
Definition: cubature.cc:3863
Definition: cubature.h:58
void populate()
Populate significant functions given information in extents.
Definition: cubature.cc:3891
double * w_
Full weights.
Definition: cubature.h:85
double * r() const
Radial nodes (including alpha scale). You do not own this.
Definition: cubature.h:262
Definition: cubature.h:347
void block(int max_points, int min_points, double max_radius)
Definition: cubature.cc:4098
virtual ~BasisExtents()
Definition: cubature.cc:3707
std::shared_ptr< Matrix > orientation() const
Orientation matrix.
Definition: cubature.h:142
double * theta() const
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:339
int nradpts
Definition: cubature.h:119
double * theta_
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:292
const std::vector< std::vector< std::shared_ptr< SphericalGrid > > > & spherical_grids() const
Spherical grids, per atom and radial point.
Definition: cubature.h:146
double alpha_
Alpha scale (for user reference)
Definition: cubature.h:226
double * x() const
The x points. You do not own this.
Definition: cubature.h:400
std::string scheme_
Scheme.
Definition: cubature.h:222
DFTGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Definition: cubature.cc:3968
int npoints() const
Number of grid points.
Definition: cubature.h:395
double delta_
Cutoff value for basis values.
Definition: cubature.h:420
std::vector< int > shells_local_to_global_
Relevant shells, local -&gt; global.
Definition: cubature.h:368
SharedVector zvec_
Definition: cubature.h:356
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:322
static void lebedev_error()
Print valid Lebedev grids and error out (throws)
Definition: cubature.cc:4729
double * z() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:332
void computeExtents()
Recompute and shell_extents_.
Definition: cubature.cc:3710
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:256
int nangpts
Definition: cubature.h:120
double z
Definition: cubature.h:59
void buildGridFromOptions()
Master builder methods.
Definition: cubature.cc:4054
double * phi_
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:290
virtual ~SphericalGrid()
Destructor.
Definition: cubature.cc:4630
MolecularGridOptions options_
A copy of the options used, for printing purposes.
Definition: cubature.h:124
double * z() const
The z points. You do not own this.
Definition: cubature.h:404
void print_details(std::string out_fname="outfile", int print=2) const
Definition: cubature.cc:4208
BasisExtents(std::shared_ptr< BasisSet > primary, double delta)
Definition: cubature.cc:3701
virtual ~DFTGrid()
Definition: cubature.cc:3986
static std::shared_ptr< RadialGrid > build_becke(int npoints, double alpha)
Build the Becke 1988 radial grid.
Definition: cubature.cc:4575
Options & options_
The Options object.
Definition: cubature.h:207
void rs(int nm, int n, double **array, double *e_vals, int matz, double **e_vecs, double toler)
double y
Definition: cubature.h:59
double * z() const
The z points. You do not own this.
Definition: cubature.h:162
void remove_distant_points(double Rcut)
Definition: cubature.cc:4139
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4641
Definition: vector3.h:37
double * y_
Full y points.
Definition: cubature.h:81
double * phi() const
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:337
void print(std::string out_fname="outfile", int print=2)
Print a trace of this BlockOPoints.
Definition: cubature.cc:3932
Definition: cubature.h:217
SharedVector yvec_
Definition: cubature.h:355
short prunescheme
Definition: cubature.h:116
int * index_
index_[fast_index] = slow_index
Definition: cubature.h:96
static std::shared_ptr< RadialGrid > build_treutler(int npoints, double alpha)
Build the Treutler-Ahlrichs 1995 radial grid (scale power = 0.6)
Definition: cubature.cc:4596
double * w() const
The weights, normalized to 1 on R3. You do not own this.
Definition: cubature.h:164
std::vector< std::shared_ptr< BlockOPoints > > blocks_
Vector of blocks.
Definition: cubature.h:99
Definition: cubature.h:63
int max_functions() const
Maximum number of funtions in a block.
Definition: cubature.h:155
static std::shared_ptr< SphericalGrid > build(const std::string &scheme, int npoints, const MassPoint *points)
Master build routines.
Definition: cubature.cc:4668
void buildGridFromOptions(MolecularGridOptions const &opt)
Build the grid.
Definition: cubature.cc:3543
std::string filename_
The filename used to optionally build the grid.
Definition: cubature.h:180
SphericalGrid()
Protected constructor.
Definition: cubature.cc:4626
int order_
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:279
Definition: cubature.h:414
int order() const
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:326
RadialGrid()
Protected constructor.
Definition: cubature.cc:4522
const std::vector< int > & functions_local_to_global() const
Relevant functions, local -&gt; global.
Definition: cubature.h:411
int npoints() const
Number of points in radial grid.
Definition: cubature.h:324
const std::vector< std::shared_ptr< BlockOPoints > > & blocks() const
Set of spatially sieved blocks of points, generated by sieve() internally.
Definition: cubature.h:169
MolecularGrid(std::shared_ptr< Molecule > molecule)
Definition: cubature.cc:4083
Definition: cubature.h:174
void refresh()
Refresh populations (if extents_-&gt;delta() changes)
Definition: cubature.h:392
double * w() const
Spherical weights, normalized to 4pi.
Definition: cubature.h:334
std::shared_ptr< BasisSet > primary_
BasisSet from extents_.
Definition: cubature.h:104
void postProcess(std::shared_ptr< BasisExtents > extents, int max_points, int min_points, double max_radius)
Sieve and block.
Definition: cubature.cc:4177
const std::vector< int > & shells_local_to_global() const
Relevant shells, local -&gt; global.
Definition: cubature.h:409
double * w_
Spherical weights, normalized to 4pi.
Definition: cubature.h:287
virtual ~PseudospectralGrid()
Definition: cubature.cc:4050
int max_points() const
Maximum number of grid points in a block.
Definition: cubature.h:153
std::shared_ptr< Molecule > molecule_
The molecule this grid is built on.
Definition: cubature.h:68
std::shared_ptr< Vector > shell_extents() const
WCS significant extent of each shell.
Definition: cubature.h:442
double x
Definition: cubature.h:59
double * x_
Pointer to x (does not own)
Definition: cubature.h:360
void set_delta(double delta)
Reset delta and recompute extents.
Definition: cubature.h:435
int max_functions_
Maximum number of functions in a block.
Definition: cubature.h:77
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:202
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4533
int npoints_
Number of points in radial grid.
Definition: cubature.h:224
std::shared_ptr< BasisSet > primary_
Basis this corresponds to.
Definition: cubature.h:418
double * y() const
The y points. You do not own this.
Definition: cubature.h:402
double * x() const
The x points. You do not own this.
Definition: cubature.h:158
int npoints_
Total points for this molecule.
Definition: cubature.h:73
virtual ~MolecularGrid()
Definition: cubature.cc:4087
std::shared_ptr< BasisExtents > extents_
Reference to the extents object.
Definition: cubature.h:372
const std::vector< std::shared_ptr< RadialGrid > > & radial_grids() const
Radial grids, per atom.
Definition: cubature.h:144
std::shared_ptr< Molecule > molecule
Definition: dx_write.cc:58
double * y() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:330
std::shared_ptr< BasisSet > basis() const
The basis set this BasisExtents is built on.
Definition: cubature.h:440
Definition: cubature.h:198
double * w_
Pointer to w (does not own)
Definition: cubature.h:366
double * z_
Pointer to z (does not own)
Definition: cubature.h:364
double * w_
Weights (including alpha and r^2)
Definition: cubature.h:230
std::string scheme_
Scheme.
Definition: cubature.h:275
std::vector< std::shared_ptr< RadialGrid > > radial_grids_
Radial grids, per atom.
Definition: cubature.h:92
Definition: liboptions.h:355
std::shared_ptr< BasisExtents > extents_
Points to basis extents, built internally.
Definition: cubature.h:102
double * y_
Spherical nodes, on the unit sphere.
Definition: cubature.h:283
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:178
static std::shared_ptr< RadialGrid > build(const std::string &scheme, int npoints, double alpha)
Master build routine.
Definition: cubature.cc:4550
double pruning_alpha
Definition: cubature.h:114
SharedVector xvec_
Data holders if requested.
Definition: cubature.h:354
double * w() const
The weights. You do not own this.
Definition: cubature.h:406
double bs_radius_alpha
Definition: cubature.h:113
short nucscheme
Definition: cubature.h:117
BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, std::shared_ptr< BasisExtents > extents)
Definition: cubature.cc:3842
double * z_
Spherical nodes, on the unit sphere.
Definition: cubature.h:285
void set_debug(int debug)
Definition: cubature.h:171
int npoints() const
Number of grid points.
Definition: cubature.h:151
short namedGrid
Definition: cubature.h:118
double maxR() const
Maximum spatial extent over all atoms.
Definition: cubature.h:444
double R_
Bounding radius of the BlockOPoints.
Definition: cubature.h:377
int max_points_
Maximum number of points in a block.
Definition: cubature.h:75
std::shared_ptr< Vector > shell_extents_
Significant extent of shells.
Definition: cubature.h:422
void build_angles()
Build the spherical angles from &lt;x,y,z&gt;, for reference.
Definition: cubature.cc:4656
Definition: cubature.h:270
double * y() const
The y points. You do not own this.
Definition: cubature.h:160
double * z_
Full z points.
Definition: cubature.h:83
double * x_
Spherical nodes, on the unit sphere.
Definition: cubature.h:281
PseudospectralGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Constructor to use for autogeneration.
Definition: cubature.cc:4043
void buildGridFromOptions(std::map< std::string, int > int_opts_map, std::map< std::string, std::string > opts_map)
Master builder methods.
Definition: cubature.cc:3990
static void initialize_lebedev()
Initialize the above arrays with the unique Lebedev grids.
Definition: cubature.cc:4691
double * w() const
Radial weights (including alpha scale and r^2). You do not own this.
Definition: cubature.h:264
int * index() const
index_[fast_index] = slow_index. You do not own this
Definition: cubature.h:148
void print(std::string out_fname="outfile", int print=2) const
Print information about the grid.
Definition: cubature.cc:4188
short radscheme
Definition: cubature.h:115
SharedVector wvec_
Definition: cubature.h:357
std::shared_ptr< Matrix > orientation_
Orientation matrix.
Definition: cubature.h:90
double alpha() const
Alpha scale (for user reference)
Definition: cubature.h:260
int debug_
Definition: cubature.h:65
double * y_
Pointer to y (does not own)
Definition: cubature.h:362
double * x() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:328
int npoints_
Number of points in radial grid.
Definition: cubature.h:277
static std::map< int, int > lebedev_mapping_
Grid npoints to order map.
Definition: cubature.h:297
int npoints() const
Number of points in radial grid.
Definition: cubature.h:258
void bound()
Compute bounding sphere.
Definition: cubature.cc:3864
double * r_
Nodes (including alpha)
Definition: cubature.h:228
int npoints_
number of points in this block
Definition: cubature.h:351
virtual ~RadialGrid()
Destructor.
Definition: cubature.cc:4526
double delta() const
The cutoff value.
Definition: cubature.h:438
std::vector< int > functions_local_to_global_
Relevant functions, local -&gt; global.
Definition: cubature.h:370
std::shared_ptr< BasisExtents > extents() const
Pointer to basis extents.
Definition: cubature.h:167
std::shared_ptr< Vector > SharedVector
Definition: adc.h:51
Vector3 xc_
Center of this BlockOPoints.
Definition: cubature.h:375
double w
Definition: cubature.h:59
double maxR_
Maximum extent.
Definition: cubature.h:424
std::vector< std::vector< std::shared_ptr< SphericalGrid > > > spherical_grids_
Spherical grids, per atom and radial point.
Definition: cubature.h:94
double * x_
Full x points.
Definition: cubature.h:79