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
118  short prunescheme;
119  short nucscheme;
120  short namedGrid; // -1 = None, 0 = SG-0, 1 = SG-1
121  int nradpts;
122  int nangpts;
123  };
124 
125  protected:
128 
129  public:
130  MolecularGrid(std::shared_ptr<Molecule> molecule);
131  virtual ~MolecularGrid();
132 
137  const std::vector<std::vector<double> >& rs, // Radial nodes, per atom
138  const std::vector<std::vector<double> >& ws, // Radial weights, per atom
139  const std::vector<std::vector<int> >& Ls); // Spherical orders, per atom
140 
142  void print(std::string out_fname = "outfile", int print = 2) const;
143  void print_details(std::string out_fname = "outfile", int print = 2) const;
144 
146  std::shared_ptr<Matrix> orientation() const { return orientation_; }
148  const std::vector<std::shared_ptr<RadialGrid> >& radial_grids() const { return radial_grids_; }
150  const std::vector<std::vector<std::shared_ptr<SphericalGrid> > >& spherical_grids() const {
151  return spherical_grids_;
152  }
154  int* index() const { return index_; }
155 
157  int npoints() const { return npoints_; }
159  int max_points() const { return max_points_; }
161  int max_functions() const { return max_functions_; }
163  size_t collocation_size() { return collocation_size_; }
164 
166  double* x() const { return x_; }
168  double* y() const { return y_; }
170  double* z() const { return z_; }
172  double* w() const { return w_; }
173 
175  std::shared_ptr<BasisExtents> extents() const { return extents_; }
177  const std::vector<std::shared_ptr<BlockOPoints> >& blocks() const { return blocks_; }
178 
179  void set_debug(int debug) { debug_ = debug; }
180 };
181 
183  protected:
185  std::shared_ptr<BasisSet> primary_;
187  std::string filename_;
188 
191 
193  void buildGridFromOptions();
194 
195  public:
197  PseudospectralGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
198  ~PseudospectralGrid() override;
199 };
200 
201 class DFTGrid : public MolecularGrid {
202  protected:
204  std::shared_ptr<BasisSet> primary_;
206  void buildGridFromOptions(std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map);
209 
210  public:
211  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
212  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary,
213  std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map, Options& options);
214  ~DFTGrid() override;
215 };
216 
217 class RadialGrid {
218  protected:
220  std::string scheme_;
222  int npoints_;
224  double alpha_;
226  double* r_;
228  double* w_;
229 
230  // ==> Standard Radial Grids <== //
231 
233  static std::shared_ptr<RadialGrid> build_becke(int npoints, double alpha);
235  static std::shared_ptr<RadialGrid> build_treutler(int npoints, double alpha);
236  // TODO: Add more grids
237 
239  RadialGrid();
240 
241  public:
242  // ==> Initializers <== //
243 
245  virtual ~RadialGrid();
246 
248  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double alpha);
250  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double* r, double* wr,
251  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  protected:
273  std::string scheme_;
275  int npoints_;
277  int order_;
279  double* x_;
281  double* y_;
283  double* z_;
285  double* w_;
286 
288  double* phi_;
290  double* theta_;
291 
292  // ==> Unique Lebedev Grids (statically stored) <== //
293 
295  static std::map<int, int> lebedev_mapping_;
297  static void initialize_lebedev();
299  static void lebedev_error();
300 
301  // ==> Utility Routines <== //
302 
304  void build_angles();
305 
307  SphericalGrid();
308 
309  public:
310  // ==> Initializers <== //
311 
313  virtual ~SphericalGrid();
314 
316  static std::shared_ptr<SphericalGrid> build(const std::string& scheme, int npoints, const MassPoint* points);
317 
318  // ==> Accessors <== //
319 
321  std::string scheme() const { return scheme_; }
323  int npoints() const { return npoints_; }
325  int order() const { return order_; }
327  double* x() const { return x_; }
329  double* y() const { return y_; }
331  double* z() const { return z_; }
333  double* w() const { return w_; }
334 
336  double* phi() const { return phi_; }
338  double* theta() const { return theta_; }
339 
341  void print(std::string out_fname = "outfile", int level = 1) const;
342 };
343 
345  protected:
347  size_t index_;
348  size_t npoints_;
349  size_t local_nbf_;
350 
356 
358  double* x_;
360  double* y_;
362  double* z_;
364  double* w_;
366  std::vector<int> shells_local_to_global_;
368  std::vector<int> functions_local_to_global_;
370  std::shared_ptr<BasisExtents> extents_;
371 
375  double R_;
376 
378  void populate();
380  void bound();
381 
382  public:
383  BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, std::shared_ptr<BasisExtents> extents);
384  BlockOPoints(size_t index, size_t npoints, double* x, double* y, double* z, double* w,
385  std::shared_ptr<BasisExtents> extents);
386  virtual ~BlockOPoints();
387 
389  void refresh() { populate(); }
390 
392  size_t npoints() const { return npoints_; }
394  size_t local_nbf() const { return local_nbf_; }
396  size_t index() const { return index_; }
398  void print(std::string out_fname = "outfile", int print = 2);
399 
401  double* x() const { return x_; }
403  double* y() const { return y_; }
405  double* z() const { return z_; }
407  double* w() const { return w_; }
408 
410  const std::vector<int>& shells_local_to_global() const { return shells_local_to_global_; }
412  const std::vector<int>& functions_local_to_global() const { return functions_local_to_global_; }
413 };
414 
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 
429  public:
430  BasisExtents(std::shared_ptr<BasisSet> primary, double delta);
431  virtual ~BasisExtents();
432 
434  void print(std::string out_fname = "outfile");
436  void set_delta(double delta) {
437  delta_ = delta;
438  computeExtents();
439  }
440 
442  double delta() const { return delta_; }
444  std::shared_ptr<BasisSet> basis() const { return primary_; }
446  std::shared_ptr<Vector> shell_extents() const { return shell_extents_; }
448  double maxR() const { return maxR_; }
449 };
450 }
451 #endif
Options & options_
The Options object.
Definition: cubature.h:190
void print(std::string out_fname="outfile")
Print a trace of these extents.
Definition: cubature.cc:3799
virtual ~BlockOPoints()
Definition: cubature.cc:3834
Definition: cubature.h:59
void populate()
Populate significant functions given information in extents.
Definition: cubature.cc:3859
double * w_
Full weights.
Definition: cubature.h:87
double * r() const
Radial nodes (including alpha scale). You do not own this.
Definition: cubature.h:262
Definition: cubature.h:344
~DFTGrid() override
Definition: cubature.cc:3938
void block(int max_points, int min_points, double max_radius)
Definition: cubature.cc:4039
virtual ~BasisExtents()
Definition: cubature.cc:3685
std::shared_ptr< Matrix > orientation() const
Orientation matrix.
Definition: cubature.h:146
double * theta() const
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:338
int nradpts
Definition: cubature.h:121
double * theta_
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:290
const std::vector< std::vector< std::shared_ptr< SphericalGrid > > > & spherical_grids() const
Spherical grids, per atom and radial point.
Definition: cubature.h:150
double alpha_
Alpha scale (for user reference)
Definition: cubature.h:224
size_t local_nbf_
Definition: cubature.h:349
double * x() const
The x points. You do not own this.
Definition: cubature.h:401
std::string scheme_
Scheme.
Definition: cubature.h:220
DFTGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Definition: cubature.cc:3927
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:366
SharedVector zvec_
Definition: cubature.h:354
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:321
static void lebedev_error()
Print valid Lebedev grids and error out (throws)
Definition: cubature.cc:4652
double * z() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:331
size_t npoints() const
Number of grid points.
Definition: cubature.h:392
void computeExtents()
Recompute and shell_extents_.
Definition: cubature.cc:3686
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:256
int nangpts
Definition: cubature.h:122
double z
Definition: cubature.h:60
void buildGridFromOptions()
Master builder methods.
Definition: cubature.cc:3999
double * phi_
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:288
virtual ~SphericalGrid()
Destructor.
Definition: cubature.cc:4558
MolecularGridOptions options_
A copy of the options used, for printing purposes.
Definition: cubature.h:127
double * z() const
The z points. You do not own this.
Definition: cubature.h:405
void print_details(std::string out_fname="outfile", int print=2) const
Definition: cubature.cc:4143
BasisExtents(std::shared_ptr< BasisSet > primary, double delta)
Definition: cubature.cc:3681
static std::shared_ptr< RadialGrid > build_becke(int npoints, double alpha)
Build the Becke 1988 radial grid.
Definition: cubature.cc:4507
Options & options_
The Options object.
Definition: cubature.h:208
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:170
void remove_distant_points(double Rcut)
Definition: cubature.cc:4082
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4568
Definition: vector3.h:40
double * y_
Full y points.
Definition: cubature.h:83
double * phi() const
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:336
void print(std::string out_fname="outfile", int print=2)
Print a trace of this BlockOPoints.
Definition: cubature.cc:3896
size_t index() const
Index of the currently owned block.
Definition: cubature.h:396
Definition: cubature.h:217
SharedVector yvec_
Definition: cubature.h:353
short prunescheme
Definition: cubature.h:118
int * index_
index_[fast_index] = slow_index
Definition: cubature.h:98
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:4528
~PseudospectralGrid() override
Definition: cubature.cc:3997
double * w() const
The weights, normalized to 1 on R3. You do not own this.
Definition: cubature.h:172
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:394
int max_functions() const
Maximum number of funtions in a block.
Definition: cubature.h:161
size_t npoints_
Definition: cubature.h:348
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:4593
void buildGridFromOptions(MolecularGridOptions const &opt)
Build the grid.
Definition: cubature.cc:3500
std::string filename_
The filename used to optionally build the grid.
Definition: cubature.h:187
SphericalGrid()
Protected constructor.
Definition: cubature.cc:4557
size_t collocation_size()
Total collocation size of all blocks.
Definition: cubature.h:163
int order_
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:277
Definition: cubature.h:415
int order() const
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:325
RadialGrid()
Protected constructor.
Definition: cubature.cc:4461
const std::vector< int > & functions_local_to_global() const
Relevant functions, local -&gt; global.
Definition: cubature.h:412
int npoints() const
Number of points in radial grid.
Definition: cubature.h:323
const std::vector< std::shared_ptr< BlockOPoints > > & blocks() const
Set of spatially sieved blocks of points, generated by sieve() internally.
Definition: cubature.h:177
MolecularGrid(std::shared_ptr< Molecule > molecule)
Definition: cubature.cc:4027
Definition: cubature.h:182
void refresh()
Refresh populations (if extents_-&gt;delta() changes)
Definition: cubature.h:389
double * w() const
Spherical weights, normalized to 4pi.
Definition: cubature.h:333
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:4114
const std::vector< int > & shells_local_to_global() const
Relevant shells, local -&gt; global.
Definition: cubature.h:410
double * w_
Spherical weights, normalized to 4pi.
Definition: cubature.h:285
int max_points() const
Maximum number of grid points in a block.
Definition: cubature.h:159
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:446
double x
Definition: cubature.h:60
double * x_
Pointer to x (does not own)
Definition: cubature.h:358
void set_delta(double delta)
Reset delta and recompute extents.
Definition: cubature.h:436
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:204
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4468
int npoints_
Number of points in radial grid.
Definition: cubature.h:222
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:403
double * x() const
The x points. You do not own this.
Definition: cubature.h:166
int npoints_
Total points for this molecule.
Definition: cubature.h:73
virtual ~MolecularGrid()
Definition: cubature.cc:4029
std::shared_ptr< BasisExtents > extents_
Reference to the extents object.
Definition: cubature.h:370
const std::vector< std::shared_ptr< RadialGrid > > & radial_grids() const
Radial grids, per atom.
Definition: cubature.h:148
std::shared_ptr< Molecule > molecule
Definition: dx_write.cc:58
double * y() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:329
std::shared_ptr< BasisSet > basis() const
The basis set this BasisExtents is built on.
Definition: cubature.h:444
Definition: cubature.h:201
double * w_
Pointer to w (does not own)
Definition: cubature.h:364
double * z_
Pointer to z (does not own)
Definition: cubature.h:362
double * w_
Weights (including alpha and r^2)
Definition: cubature.h:228
std::string scheme_
Scheme.
Definition: cubature.h:273
std::vector< std::shared_ptr< RadialGrid > > radial_grids_
Radial grids, per atom.
Definition: cubature.h:94
Definition: liboptions.h:352
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:281
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:185
static std::shared_ptr< RadialGrid > build(const std::string &scheme, int npoints, double alpha)
Master build routine.
Definition: cubature.cc:4483
double pruning_alpha
Definition: cubature.h:116
SharedVector xvec_
Data holders if requested.
Definition: cubature.h:352
double * w() const
The weights. You do not own this.
Definition: cubature.h:407
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:3812
double * z_
Spherical nodes, on the unit sphere.
Definition: cubature.h:283
void set_debug(int debug)
Definition: cubature.h:179
int npoints() const
Number of grid points.
Definition: cubature.h:157
short namedGrid
Definition: cubature.h:120
double maxR() const
Maximum spatial extent over all atoms.
Definition: cubature.h:448
double R_
Bounding radius of the BlockOPoints.
Definition: cubature.h:375
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:4582
Definition: cubature.h:270
double * y() const
The y points. You do not own this.
Definition: cubature.h:168
double * z_
Full z points.
Definition: cubature.h:85
double * x_
Spherical nodes, on the unit sphere.
Definition: cubature.h:279
PseudospectralGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Constructor to use for autogeneration.
Definition: cubature.cc:3992
void buildGridFromOptions(std::map< std::string, int > int_opts_map, std::map< std::string, std::string > opts_map)
Master builder methods.
Definition: cubature.cc:3940
static void initialize_lebedev()
Initialize the above arrays with the unique Lebedev grids.
Definition: cubature.cc:4615
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:154
void print(std::string out_fname="outfile", int print=2) const
Print information about the grid.
Definition: cubature.cc:4125
short radscheme
Definition: cubature.h:117
SharedVector wvec_
Definition: cubature.h:355
std::shared_ptr< Matrix > orientation_
Orientation matrix.
Definition: cubature.h:92
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:360
double * x() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:327
int npoints_
Number of points in radial grid.
Definition: cubature.h:275
static std::map< int, int > lebedev_mapping_
Grid npoints to order map.
Definition: cubature.h:295
int npoints() const
Number of points in radial grid.
Definition: cubature.h:258
void bound()
Compute bounding sphere.
Definition: cubature.cc:3835
double * r_
Nodes (including alpha)
Definition: cubature.h:226
virtual ~RadialGrid()
Destructor.
Definition: cubature.cc:4462
double delta() const
The cutoff value.
Definition: cubature.h:442
std::vector< int > functions_local_to_global_
Relevant functions, local -&gt; global.
Definition: cubature.h:368
std::shared_ptr< BasisExtents > extents() const
Pointer to basis extents.
Definition: cubature.h:175
std::shared_ptr< Vector > SharedVector
Definition: adc.h:51
Vector3 xc_
Center of this BlockOPoints.
Definition: cubature.h:373
double w
Definition: cubature.h:60
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:96
size_t index_
number of points in this block
Definition: cubature.h:347
double * x_
Full x points.
Definition: cubature.h:81