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-2019 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 #include "psi4/pragma.h"
34 
35 #include "psi4/libmints/vector3.h"
36 #include "psi4/libmints/typedefs.h"
37 
38 #include <map>
39 #include <vector>
40 
41 namespace psi {
42 
43 class BasisSet;
44 class Matrix;
45 class Vector;
46 class BasisExtents;
47 class BlockOPoints;
48 class RadialGrid;
49 class SphericalGrid;
50 class Options;
51 
52 // This is an auxiliary structure used internally by the grid-builder class.
53 // Apparently, for performance reasons, it is not good for the final molecular grid
54 // to consist of structures like this one.
55 //
56 // RMP: You got that right. Calling nuclear weights one point at a time is another
57 // great way to get a 1000x slowdown. What an incredible smell you've discovered!
58 //
59 struct MassPoint {
60  double x, y, z, w;
61 };
62 
64  protected:
65  int debug_;
66 
68  std::shared_ptr<Molecule> molecule_;
69 
70  // ==> Fast Grid Specification <== //
71 
73  int npoints_;
78  // The total collocation size
81  double* x_;
83  double* y_;
85  double* z_;
87  double* w_;
88 
89  // ==> Clean Grid Specification <== //
90 
92  std::shared_ptr<Matrix> orientation_;
94  std::vector<std::shared_ptr<RadialGrid> > radial_grids_;
96  std::vector<std::vector<std::shared_ptr<SphericalGrid> > > spherical_grids_;
98  int* index_;
99 
101  std::vector<std::shared_ptr<BlockOPoints> > blocks_;
102 
104  std::shared_ptr<BasisExtents> extents_;
106  std::shared_ptr<BasisSet> primary_;
107 
109  void postProcess(std::shared_ptr<BasisExtents> extents, int max_points, int min_points, double max_radius);
110  void remove_distant_points(double Rcut);
111  void block(int max_points, int min_points, double max_radius);
112 
113  public:
117  short radscheme; // Effectively an enumeration
119  short nucscheme;
120  short namedGrid; // -1 = None, 0 = SG-0, 1 = SG-1
121  int nradpts;
122  int nangpts;
123  std::string prunescheme;
124  std::string prunetype;
125  };
126 
127  protected:
130 
131  public:
132  MolecularGrid(std::shared_ptr<Molecule> molecule);
133  virtual ~MolecularGrid();
134 
139  const std::vector<std::vector<double> >& rs, // Radial nodes, per atom
140  const std::vector<std::vector<double> >& ws, // Radial weights, per atom
141  const std::vector<std::vector<int> >& Ls); // Spherical orders, per atom
142 
144  void print(std::string out_fname = "outfile", int print = 2) const;
145  void print_details(std::string out_fname = "outfile", int print = 2) const;
146 
148  std::shared_ptr<Matrix> orientation() const { return orientation_; }
150  const std::vector<std::shared_ptr<RadialGrid> >& radial_grids() const { return radial_grids_; }
152  const std::vector<std::vector<std::shared_ptr<SphericalGrid> > >& spherical_grids() const {
153  return spherical_grids_;
154  }
156  int* index() const { return index_; }
157 
159  int npoints() const { return npoints_; }
161  int max_points() const { return max_points_; }
163  int max_functions() const { return max_functions_; }
165  size_t collocation_size() { return collocation_size_; }
166 
168  double* x() const { return x_; }
170  double* y() const { return y_; }
172  double* z() const { return z_; }
174  double* w() const { return w_; }
175 
177  std::shared_ptr<BasisExtents> extents() const { return extents_; }
179  const std::vector<std::shared_ptr<BlockOPoints> >& blocks() const { return blocks_; }
180 
181  void set_debug(int debug) { debug_ = debug; }
182 };
183 
185  protected:
187  std::shared_ptr<BasisSet> primary_;
189  std::string filename_;
190 
193 
195  void buildGridFromOptions();
196 
197  public:
199  PseudospectralGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
200  ~PseudospectralGrid() override;
201 };
202 
203 class DFTGrid : public MolecularGrid {
204  protected:
206  std::shared_ptr<BasisSet> primary_;
208  void buildGridFromOptions(std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map);
211 
212  public:
213  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
214  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary,
215  std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map, Options& options);
216  ~DFTGrid() override;
217 };
218 
219 class RadialGrid {
220  protected:
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, int Z);
237  static std::shared_ptr<RadialGrid> build_treutler(int npoints, double alpha, int Z);
238  // TODO: Add more grids
239 
241  RadialGrid();
242 
243  public:
244  // ==> Initializers <== //
245 
247  virtual ~RadialGrid();
248 
250  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double alpha, int Z);
252  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double* r, double* wr,
253  double alpha, int Z);
254 
255  // ==> Accessors <== //
256 
258  std::string scheme() const { return scheme_; }
260  int npoints() const { return npoints_; }
262  double alpha() const { return alpha_; }
264  double* r() const { return r_; }
266  double* w() const { return w_; }
267 
269  void print(std::string out_fname = "outfile", int level = 1) const;
270 };
271 
273  protected:
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 
311  public:
312  // ==> Initializers <== //
313 
315  virtual ~SphericalGrid();
316 
318  static std::shared_ptr<SphericalGrid> build(const std::string& scheme, int npoints, const MassPoint* points);
319 
320  // ==> Accessors <== //
321 
323  std::string scheme() const { return scheme_; }
325  int npoints() const { return npoints_; }
327  int order() const { return order_; }
329  double* x() const { return x_; }
331  double* y() const { return y_; }
333  double* z() const { return z_; }
335  double* w() const { return w_; }
336 
338  double* phi() const { return phi_; }
340  double* theta() const { return theta_; }
341 
343  void print(std::string out_fname = "outfile", int level = 1) const;
344 };
345 
347  protected:
349  size_t index_;
350  size_t npoints_;
351  size_t local_nbf_;
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:
385  BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, std::shared_ptr<BasisExtents> extents);
386  BlockOPoints(size_t index, size_t npoints, double* x, double* y, double* z, double* w,
387  std::shared_ptr<BasisExtents> extents);
388  virtual ~BlockOPoints();
389 
391  void refresh() { populate(); }
392 
394  size_t npoints() const { return npoints_; }
396  size_t local_nbf() const { return local_nbf_; }
398  size_t index() const { return index_; }
400  void print(std::string out_fname = "outfile", int print = 2);
401 
403  double* x() const { return x_; }
405  double* y() const { return y_; }
407  double* z() const { return z_; }
409  double* w() const { return w_; }
410 
412  const std::vector<int>& shells_local_to_global() const { return shells_local_to_global_; }
414  const std::vector<int>& functions_local_to_global() const { return functions_local_to_global_; }
415 };
416 
418  protected:
420  std::shared_ptr<BasisSet> primary_;
422  double delta_;
424  std::shared_ptr<Vector> shell_extents_;
426  double maxR_;
427 
429  void computeExtents();
430 
431  public:
432  BasisExtents(std::shared_ptr<BasisSet> primary, double delta);
433  virtual ~BasisExtents();
434 
436  void print(std::string out_fname = "outfile");
438  void set_delta(double delta) {
439  delta_ = delta;
440  computeExtents();
441  }
442 
444  double delta() const { return delta_; }
446  std::shared_ptr<BasisSet> basis() const { return primary_; }
448  std::shared_ptr<Vector> shell_extents() const { return shell_extents_; }
450  double maxR() const { return maxR_; }
451 };
452 }
453 #endif
Options & options_
The Options object.
Definition: cubature.h:192
void print(std::string out_fname="outfile")
Print a trace of these extents.
Definition: cubature.cc:3924
virtual ~BlockOPoints()
Definition: cubature.cc:3959
Definition: cubature.h:59
void populate()
Populate significant functions given information in extents.
Definition: cubature.cc:3984
double * w_
Full weights.
Definition: cubature.h:87
std::string prunescheme
Definition: cubature.h:123
double * r() const
Radial nodes (including alpha scale). You do not own this.
Definition: cubature.h:264
Definition: cubature.h:346
~DFTGrid() override
Definition: cubature.cc:4063
void block(int max_points, int min_points, double max_radius)
Definition: cubature.cc:4187
virtual ~BasisExtents()
Definition: cubature.cc:3810
std::shared_ptr< Matrix > orientation() const
Orientation matrix.
Definition: cubature.h:148
double * theta() const
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:340
int nradpts
Definition: cubature.h:121
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:152
double alpha_
Alpha scale (for user reference)
Definition: cubature.h:226
size_t local_nbf_
Definition: cubature.h:351
double * x() const
The x points. You do not own this.
Definition: cubature.h:403
std::string scheme_
Scheme.
Definition: cubature.h:222
DFTGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Definition: cubature.cc:4052
double delta_
Cutoff value for basis values.
Definition: cubature.h:422
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:323
static void lebedev_error()
Print valid Lebedev grids and error out (throws)
Definition: cubature.cc:4823
double * z() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:333
size_t npoints() const
Number of grid points.
Definition: cubature.h:394
void computeExtents()
Recompute and shell_extents_.
Definition: cubature.cc:3811
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:258
int nangpts
Definition: cubature.h:122
double z
Definition: cubature.h:60
void buildGridFromOptions()
Master builder methods.
Definition: cubature.cc:4147
double * phi_
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:290
virtual ~SphericalGrid()
Destructor.
Definition: cubature.cc:4729
MolecularGridOptions options_
A copy of the options used, for printing purposes.
Definition: cubature.h:129
double * z() const
The z points. You do not own this.
Definition: cubature.h:407
void print_details(std::string out_fname="outfile", int print=2) const
Definition: cubature.cc:4297
BasisExtents(std::shared_ptr< BasisSet > primary, double delta)
Definition: cubature.cc:3806
Options & options_
The Options object.
Definition: cubature.h:210
void rs(int nm, int n, double **array, double *e_vals, int matz, double **e_vecs, double toler)
double y
Definition: cubature.h:60
double * z() const
The z points. You do not own this.
Definition: cubature.h:172
void remove_distant_points(double Rcut)
Definition: cubature.cc:4230
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4739
Definition: vector3.h:40
static std::shared_ptr< RadialGrid > build_treutler(int npoints, double alpha, int Z)
Build the Treutler-Ahlrichs 1995 radial grid (scale power = 0.6)
Definition: cubature.cc:4683
double * y_
Full y points.
Definition: cubature.h:83
double * phi() const
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:338
void print(std::string out_fname="outfile", int print=2)
Print a trace of this BlockOPoints.
Definition: cubature.cc:4021
size_t index() const
Index of the currently owned block.
Definition: cubature.h:398
Definition: cubature.h:219
SharedVector yvec_
Definition: cubature.h:355
int * index_
index_[fast_index] = slow_index
Definition: cubature.h:98
~PseudospectralGrid() override
Definition: cubature.cc:4145
double * w() const
The weights, normalized to 1 on R3. You do not own this.
Definition: cubature.h:174
std::vector< std::shared_ptr< BlockOPoints > > blocks_
Vector of blocks.
Definition: cubature.h:101
Definition: cubature.h:63
size_t local_nbf() const
Number of basis functions in the block.
Definition: cubature.h:396
int max_functions() const
Maximum number of funtions in a block.
Definition: cubature.h:163
size_t npoints_
Definition: cubature.h:350
size_t collocation_size_
Definition: cubature.h:79
static std::shared_ptr< SphericalGrid > build(const std::string &scheme, int npoints, const MassPoint *points)
Master build routines.
Definition: cubature.cc:4764
void buildGridFromOptions(MolecularGridOptions const &opt)
Build the grid.
Definition: cubature.cc:3615
std::string filename_
The filename used to optionally build the grid.
Definition: cubature.h:189
SphericalGrid()
Protected constructor.
Definition: cubature.cc:4728
size_t collocation_size()
Total collocation size of all blocks.
Definition: cubature.h:165
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:417
int order() const
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:327
RadialGrid()
Protected constructor.
Definition: cubature.cc:4616
const std::vector< int > & functions_local_to_global() const
Relevant functions, local -&gt; global.
Definition: cubature.h:414
int npoints() const
Number of points in radial grid.
Definition: cubature.h:325
static std::shared_ptr< RadialGrid > build_becke(int npoints, double alpha, int Z)
Build the Becke 1988 radial grid.
Definition: cubature.cc:4662
static std::shared_ptr< RadialGrid > build(const std::string &scheme, int npoints, double alpha, int Z)
Master build routine.
Definition: cubature.cc:4638
const std::vector< std::shared_ptr< BlockOPoints > > & blocks() const
Set of spatially sieved blocks of points, generated by sieve() internally.
Definition: cubature.h:179
MolecularGrid(std::shared_ptr< Molecule > molecule)
Definition: cubature.cc:4175
Definition: cubature.h:184
void refresh()
Refresh populations (if extents_-&gt;delta() changes)
Definition: cubature.h:391
double * w() const
Spherical weights, normalized to 4pi.
Definition: cubature.h:335
std::shared_ptr< BasisSet > primary_
BasisSet from extents_.
Definition: cubature.h:106
void postProcess(std::shared_ptr< BasisExtents > extents, int max_points, int min_points, double max_radius)
Sieve and block.
Definition: cubature.cc:4262
const std::vector< int > & shells_local_to_global() const
Relevant shells, local -&gt; global.
Definition: cubature.h:412
double * w_
Spherical weights, normalized to 4pi.
Definition: cubature.h:287
int max_points() const
Maximum number of grid points in a block.
Definition: cubature.h:161
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:448
double x
Definition: cubature.h:60
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:438
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:206
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4623
int npoints_
Number of points in radial grid.
Definition: cubature.h:224
std::shared_ptr< BasisSet > primary_
Basis this corresponds to.
Definition: cubature.h:420
double * y() const
The y points. You do not own this.
Definition: cubature.h:405
double * x() const
The x points. You do not own this.
Definition: cubature.h:168
int npoints_
Total points for this molecule.
Definition: cubature.h:73
virtual ~MolecularGrid()
Definition: cubature.cc:4177
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:150
std::shared_ptr< Molecule > molecule
Definition: dx_write.cc:58
double * y() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:331
std::shared_ptr< BasisSet > basis() const
The basis set this BasisExtents is built on.
Definition: cubature.h:446
Definition: cubature.h:203
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:94
Definition: liboptions.h:353
std::shared_ptr< BasisExtents > extents_
Points to basis extents, built internally.
Definition: cubature.h:104
double * y_
Spherical nodes, on the unit sphere.
Definition: cubature.h:283
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:187
double pruning_alpha
Definition: cubature.h:116
SharedVector xvec_
Data holders if requested.
Definition: cubature.h:354
double * w() const
The weights. You do not own this.
Definition: cubature.h:409
double bs_radius_alpha
Definition: cubature.h:115
short nucscheme
Definition: cubature.h:119
BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, std::shared_ptr< BasisExtents > extents)
Definition: cubature.cc:3937
double * z_
Spherical nodes, on the unit sphere.
Definition: cubature.h:285
void set_debug(int debug)
Definition: cubature.h:181
int npoints() const
Number of grid points.
Definition: cubature.h:159
std::string prunetype
Definition: cubature.h:124
short namedGrid
Definition: cubature.h:120
double maxR() const
Maximum spatial extent over all atoms.
Definition: cubature.h:450
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:424
void build_angles()
Build the spherical angles from &lt;x,y,z&gt;, for reference.
Definition: cubature.cc:4753
Definition: cubature.h:272
double * y() const
The y points. You do not own this.
Definition: cubature.h:170
double * z_
Full z points.
Definition: cubature.h:85
short prunefunction
Definition: cubature.h:118
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:4140
void buildGridFromOptions(std::map< std::string, int > int_opts_map, std::map< std::string, std::string > opts_map)
Master builder methods.
Definition: cubature.cc:4065
static void initialize_lebedev()
Initialize the above arrays with the unique Lebedev grids.
Definition: cubature.cc:4786
double * w() const
Radial weights (including alpha scale and r^2). You do not own this.
Definition: cubature.h:266
int * index() const
index_[fast_index] = slow_index. You do not own this
Definition: cubature.h:156
void print(std::string out_fname="outfile", int print=2) const
Print information about the grid.
Definition: cubature.cc:4273
short radscheme
Definition: cubature.h:117
SharedVector wvec_
Definition: cubature.h:357
std::shared_ptr< Matrix > orientation_
Orientation matrix.
Definition: cubature.h:92
double alpha() const
Alpha scale (for user reference)
Definition: cubature.h:262
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:329
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:260
void bound()
Compute bounding sphere.
Definition: cubature.cc:3960
double * r_
Nodes (including alpha)
Definition: cubature.h:228
virtual ~RadialGrid()
Destructor.
Definition: cubature.cc:4617
double delta() const
The cutoff value.
Definition: cubature.h:444
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:177
std::shared_ptr< Vector > SharedVector
Definition: adc.h:51
Vector3 xc_
Center of this BlockOPoints.
Definition: cubature.h:375
double w
Definition: cubature.h:60
double maxR_
Maximum extent.
Definition: cubature.h:426
std::vector< std::vector< std::shared_ptr< SphericalGrid > > > spherical_grids_
Spherical grids, per atom and radial point.
Definition: cubature.h:96
size_t index_
number of points in this block
Definition: cubature.h:349
double * x_
Full x points.
Definition: cubature.h:81