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-2018 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 
63 
65 protected:
66  int debug_;
67 
69  std::shared_ptr<Molecule> molecule_;
70 
71  // ==> Fast Grid Specification <== //
72 
74  int npoints_;
79  // The total collocation size
82  double* x_;
84  double* y_;
86  double* z_;
88  double* w_;
89 
90  // ==> Clean Grid Specification <== //
91 
93  std::shared_ptr<Matrix> orientation_;
95  std::vector<std::shared_ptr<RadialGrid> > radial_grids_;
97  std::vector<std::vector<std::shared_ptr<SphericalGrid> > > spherical_grids_;
99  int* index_;
100 
102  std::vector<std::shared_ptr<BlockOPoints> > blocks_;
103 
105  std::shared_ptr<BasisExtents> extents_;
107  std::shared_ptr<BasisSet> primary_;
108 
110  void postProcess(std::shared_ptr<BasisExtents> extents, int max_points, int min_points, double max_radius);
111  void remove_distant_points(double Rcut);
112  void block(int max_points, int min_points, double max_radius);
113 
114 public:
118  short radscheme; // Effectively an enumeration
119  short prunescheme;
120  short nucscheme;
121  short namedGrid; // -1 = None, 0 = SG-0, 1 = SG-1
122  int nradpts;
123  int nangpts;
124  };
125 protected:
128 public:
129  MolecularGrid(std::shared_ptr<Molecule> molecule);
130  virtual ~MolecularGrid();
131 
136  const std::vector<std::vector<double> >& rs, // Radial nodes, per atom
137  const std::vector<std::vector<double> >& ws, // Radial weights, per atom
138  const std::vector<std::vector<int> >& Ls); // Spherical orders, per atom
139 
141  void print(std::string out_fname = "outfile", int print = 2) const;
142  void print_details(std::string out_fname = "outfile", int print = 2) const;
143 
145  std::shared_ptr<Matrix> orientation() const { return orientation_; }
147  const std::vector<std::shared_ptr<RadialGrid> >& radial_grids() const { return radial_grids_; }
149  const std::vector<std::vector<std::shared_ptr<SphericalGrid> > >& spherical_grids() const { return spherical_grids_; }
151  int* index() const { return index_; }
152 
154  int npoints() const { return npoints_; }
156  int max_points() const { return max_points_; }
158  int max_functions() const { return max_functions_; }
160  size_t collocation_size() { return collocation_size_; }
161 
163  double* x() const { return x_; }
165  double* y() const { return y_; }
167  double* z() const { return z_; }
169  double* w() const { return w_; }
170 
172  std::shared_ptr<BasisExtents> extents() const { return extents_; }
174  const std::vector<std::shared_ptr<BlockOPoints> >& blocks() const { return blocks_; }
175 
176  void set_debug(int debug) { debug_ = debug; }
177 };
178 
180 
181 protected:
183  std::shared_ptr<BasisSet> primary_;
185  std::string filename_;
186 
189 
191  void buildGridFromOptions();
192 
193 public:
194 
196  PseudospectralGrid(std::shared_ptr<Molecule> molecule,
197  std::shared_ptr<BasisSet> primary,
198  Options& options);
199  virtual ~PseudospectralGrid();
200 
201 };
202 
203 class DFTGrid : public MolecularGrid {
204 
205 protected:
207  std::shared_ptr<BasisSet> primary_;
209  void buildGridFromOptions(std::map<std::string, int> int_opts_map,
210  std::map<std::string, std::string> opts_map);
213 
214 public:
215  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary, Options& options);
216  DFTGrid(std::shared_ptr<Molecule> molecule, std::shared_ptr<BasisSet> primary,
217  std::map<std::string, int> int_opts_map, std::map<std::string, std::string> opts_map,
218  Options& options);
219 virtual ~DFTGrid();
220 };
221 
222 class RadialGrid {
223 
224 protected:
225 
227  std::string scheme_;
229  int npoints_;
231  double alpha_;
233  double* r_;
235  double* w_;
236 
237  // ==> Standard Radial Grids <== //
238 
240  static std::shared_ptr<RadialGrid> build_becke(int npoints, double alpha);
242  static std::shared_ptr<RadialGrid> build_treutler(int npoints, double alpha);
243  // TODO: Add more grids
244 
246  RadialGrid();
247 public:
248  // ==> Initializers <== //
249 
251  virtual ~RadialGrid();
252 
254  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double alpha);
256  static std::shared_ptr<RadialGrid> build(const std::string& scheme, int npoints, double* r, double* wr, double alpha);
257 
258  // ==> Accessors <== //
259 
261  std::string scheme() const { return scheme_; }
263  int npoints() const { return npoints_; }
265  double alpha() const { return alpha_; }
267  double* r() const { return r_; }
269  double* w() const { return w_; }
270 
272  void print(std::string out_fname = "outfile", int level = 1) const;
273 };
274 
276 
277 protected:
278 
280  std::string scheme_;
282  int npoints_;
284  int order_;
286  double* x_;
288  double* y_;
290  double* z_;
292  double* w_;
293 
295  double* phi_;
297  double* theta_;
298 
299  // ==> Unique Lebedev Grids (statically stored) <== //
300 
302  static std::map<int, int> lebedev_mapping_;
304  static void initialize_lebedev();
306  static void lebedev_error();
307 
308  // ==> Utility Routines <== //
309 
311  void build_angles();
312 
314  SphericalGrid();
315 public:
316  // ==> Initializers <== //
317 
319  virtual ~SphericalGrid();
320 
322  static std::shared_ptr<SphericalGrid> build(const std::string& scheme, int npoints, const MassPoint* points);
323 
324  // ==> Accessors <== //
325 
327  std::string scheme() const { return scheme_; }
329  int npoints() const { return npoints_; }
331  int order() const { return order_; }
333  double* x() const { return x_; }
335  double* y() const { return y_; }
337  double* z() const { return z_; }
339  double* w() const { return w_; }
340 
342  double* phi() const { return phi_; }
344  double* theta() const { return theta_; }
345 
347  void print(std::string out_fname = "outfile", int level = 1) const;
348 
349 
350 };
351 
353 
354 protected:
356  size_t index_;
357  size_t npoints_;
358  size_t local_nbf_;
359 
365 
367  double* x_;
369  double* y_;
371  double* z_;
373  double* w_;
375  std::vector<int> shells_local_to_global_;
377  std::vector<int> functions_local_to_global_;
379  std::shared_ptr<BasisExtents> extents_;
380 
384  double R_;
385 
387  void populate();
389  void bound();
390 
391 public:
393  std::shared_ptr<BasisExtents> extents);
394  BlockOPoints(size_t index, size_t npoints, double* x, double* y, double* z, double* w,
395  std::shared_ptr<BasisExtents> extents);
396  virtual ~BlockOPoints();
397 
399  void refresh() { populate(); }
400 
402  size_t npoints() const { return npoints_; }
404  size_t local_nbf() const { return local_nbf_; }
406  size_t index() const { return index_; }
408  void print(std::string out_fname = "outfile", int print = 2);
409 
411  double* x() const { return x_; }
413  double* y() const { return y_; }
415  double* z() const { return z_; }
417  double* w() const { return w_; }
418 
420  const std::vector<int>& shells_local_to_global() const { return shells_local_to_global_; }
422  const std::vector<int>& functions_local_to_global() const { return functions_local_to_global_; }
423 };
424 
426 
427 protected:
429  std::shared_ptr<BasisSet> primary_;
431  double delta_;
433  std::shared_ptr<Vector> shell_extents_;
435  double maxR_;
436 
438  void computeExtents();
439 public:
440  BasisExtents(std::shared_ptr<BasisSet> primary, double delta);
441  virtual ~BasisExtents();
442 
444  void print(std::string out_fname = "outfile");
446  void set_delta(double delta) { delta_ = delta; computeExtents(); }
447 
449  double delta() const { return delta_; }
451  std::shared_ptr<BasisSet> basis() const { return primary_; }
453  std::shared_ptr<Vector> shell_extents() const { return shell_extents_; }
455  double maxR() const { return maxR_; }
456 };
457 
458 }
459 #endif
Options & options_
The Options object.
Definition: cubature.h:188
void print(std::string out_fname="outfile")
Print a trace of these extents.
Definition: cubature.cc:3855
virtual ~BlockOPoints()
Definition: cubature.cc:3893
Definition: cubature.h:59
void populate()
Populate significant functions given information in extents.
Definition: cubature.cc:3921
double * w_
Full weights.
Definition: cubature.h:88
double * r() const
Radial nodes (including alpha scale). You do not own this.
Definition: cubature.h:267
Definition: cubature.h:352
void block(int max_points, int min_points, double max_radius)
Definition: cubature.cc:4129
virtual ~BasisExtents()
Definition: cubature.cc:3736
std::shared_ptr< Matrix > orientation() const
Orientation matrix.
Definition: cubature.h:145
double * theta() const
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:344
int nradpts
Definition: cubature.h:122
double * theta_
Spherical nodes, in spherical coordinates (inclination)
Definition: cubature.h:297
const std::vector< std::vector< std::shared_ptr< SphericalGrid > > > & spherical_grids() const
Spherical grids, per atom and radial point.
Definition: cubature.h:149
double alpha_
Alpha scale (for user reference)
Definition: cubature.h:231
size_t local_nbf_
Definition: cubature.h:358
double * x() const
The x points. You do not own this.
Definition: cubature.h:411
std::string scheme_
Scheme.
Definition: cubature.h:227
DFTGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Definition: cubature.cc:3999
double delta_
Cutoff value for basis values.
Definition: cubature.h:431
std::vector< int > shells_local_to_global_
Relevant shells, local -&gt; global.
Definition: cubature.h:375
SharedVector zvec_
Definition: cubature.h:363
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:327
static void lebedev_error()
Print valid Lebedev grids and error out (throws)
Definition: cubature.cc:4773
double * z() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:337
size_t npoints() const
Number of grid points.
Definition: cubature.h:402
void computeExtents()
Recompute and shell_extents_.
Definition: cubature.cc:3739
std::string scheme() const
Scheme of this radial grid.
Definition: cubature.h:261
int nangpts
Definition: cubature.h:123
double z
Definition: cubature.h:60
void buildGridFromOptions()
Master builder methods.
Definition: cubature.cc:4085
double * phi_
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:295
virtual ~SphericalGrid()
Destructor.
Definition: cubature.cc:4674
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:415
void print_details(std::string out_fname="outfile", int print=2) const
Definition: cubature.cc:4243
BasisExtents(std::shared_ptr< BasisSet > primary, double delta)
Definition: cubature.cc:3730
virtual ~DFTGrid()
Definition: cubature.cc:4017
static std::shared_ptr< RadialGrid > build_becke(int npoints, double alpha)
Build the Becke 1988 radial grid.
Definition: cubature.cc:4619
Options & options_
The Options object.
Definition: cubature.h:212
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:167
void remove_distant_points(double Rcut)
Definition: cubature.cc:4173
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4685
Definition: vector3.h:38
double * y_
Full y points.
Definition: cubature.h:84
double * phi() const
Spherical nodes, in spherical coordinates (azimuth)
Definition: cubature.h:342
void print(std::string out_fname="outfile", int print=2)
Print a trace of this BlockOPoints.
Definition: cubature.cc:3963
size_t index() const
Index of the currently owned block.
Definition: cubature.h:406
Definition: cubature.h:222
SharedVector yvec_
Definition: cubature.h:362
short prunescheme
Definition: cubature.h:119
int * index_
index_[fast_index] = slow_index
Definition: cubature.h:99
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:4640
double * w() const
The weights, normalized to 1 on R3. You do not own this.
Definition: cubature.h:169
std::vector< std::shared_ptr< BlockOPoints > > blocks_
Vector of blocks.
Definition: cubature.h:102
Definition: cubature.h:64
size_t local_nbf() const
Number of basis functions in the block.
Definition: cubature.h:404
int max_functions() const
Maximum number of funtions in a block.
Definition: cubature.h:158
size_t npoints_
Definition: cubature.h:357
size_t collocation_size_
Definition: cubature.h:80
static std::shared_ptr< SphericalGrid > build(const std::string &scheme, int npoints, const MassPoint *points)
Master build routines.
Definition: cubature.cc:4712
void buildGridFromOptions(MolecularGridOptions const &opt)
Build the grid.
Definition: cubature.cc:3547
std::string filename_
The filename used to optionally build the grid.
Definition: cubature.h:185
SphericalGrid()
Protected constructor.
Definition: cubature.cc:4670
size_t collocation_size()
Total collocation size of all blocks.
Definition: cubature.h:160
int order_
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:284
Definition: cubature.h:425
int order() const
Order of spherical harmonics in spherical grid (integrates products up to L_tot = 2 * order_ + 1) ...
Definition: cubature.h:331
RadialGrid()
Protected constructor.
Definition: cubature.cc:4566
const std::vector< int > & functions_local_to_global() const
Relevant functions, local -&gt; global.
Definition: cubature.h:422
int npoints() const
Number of points in radial grid.
Definition: cubature.h:329
const std::vector< std::shared_ptr< BlockOPoints > > & blocks() const
Set of spatially sieved blocks of points, generated by sieve() internally.
Definition: cubature.h:174
MolecularGrid(std::shared_ptr< Molecule > molecule)
Definition: cubature.cc:4114
Definition: cubature.h:179
void refresh()
Refresh populations (if extents_-&gt;delta() changes)
Definition: cubature.h:399
double * w() const
Spherical weights, normalized to 4pi.
Definition: cubature.h:339
std::shared_ptr< BasisSet > primary_
BasisSet from extents_.
Definition: cubature.h:107
void postProcess(std::shared_ptr< BasisExtents > extents, int max_points, int min_points, double max_radius)
Sieve and block.
Definition: cubature.cc:4211
const std::vector< int > & shells_local_to_global() const
Relevant shells, local -&gt; global.
Definition: cubature.h:420
double * w_
Spherical weights, normalized to 4pi.
Definition: cubature.h:292
virtual ~PseudospectralGrid()
Definition: cubature.cc:4081
int max_points() const
Maximum number of grid points in a block.
Definition: cubature.h:156
std::shared_ptr< Molecule > molecule_
The molecule this grid is built on.
Definition: cubature.h:69
std::shared_ptr< Vector > shell_extents() const
WCS significant extent of each shell.
Definition: cubature.h:453
double x
Definition: cubature.h:60
double * x_
Pointer to x (does not own)
Definition: cubature.h:367
void set_delta(double delta)
Reset delta and recompute extents.
Definition: cubature.h:446
int max_functions_
Maximum number of functions in a block.
Definition: cubature.h:78
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:207
void print(std::string out_fname="outfile", int level=1) const
Reflection.
Definition: cubature.cc:4577
int npoints_
Number of points in radial grid.
Definition: cubature.h:229
std::shared_ptr< BasisSet > primary_
Basis this corresponds to.
Definition: cubature.h:429
double * y() const
The y points. You do not own this.
Definition: cubature.h:413
double * x() const
The x points. You do not own this.
Definition: cubature.h:163
int npoints_
Total points for this molecule.
Definition: cubature.h:74
virtual ~MolecularGrid()
Definition: cubature.cc:4118
std::shared_ptr< BasisExtents > extents_
Reference to the extents object.
Definition: cubature.h:379
const std::vector< std::shared_ptr< RadialGrid > > & radial_grids() const
Radial grids, per atom.
Definition: cubature.h:147
std::shared_ptr< Molecule > molecule
Definition: dx_write.cc:58
double * y() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:335
std::shared_ptr< BasisSet > basis() const
The basis set this BasisExtents is built on.
Definition: cubature.h:451
Definition: cubature.h:203
double * w_
Pointer to w (does not own)
Definition: cubature.h:373
double * z_
Pointer to z (does not own)
Definition: cubature.h:371
double * w_
Weights (including alpha and r^2)
Definition: cubature.h:235
std::string scheme_
Scheme.
Definition: cubature.h:280
std::vector< std::shared_ptr< RadialGrid > > radial_grids_
Radial grids, per atom.
Definition: cubature.h:95
Definition: liboptions.h:353
std::shared_ptr< BasisExtents > extents_
Points to basis extents, built internally.
Definition: cubature.h:105
double * y_
Spherical nodes, on the unit sphere.
Definition: cubature.h:288
std::shared_ptr< BasisSet > primary_
The primary basis.
Definition: cubature.h:183
static std::shared_ptr< RadialGrid > build(const std::string &scheme, int npoints, double alpha)
Master build routine.
Definition: cubature.cc:4594
double pruning_alpha
Definition: cubature.h:117
SharedVector xvec_
Data holders if requested.
Definition: cubature.h:361
double * w() const
The weights. You do not own this.
Definition: cubature.h:417
double bs_radius_alpha
Definition: cubature.h:116
short nucscheme
Definition: cubature.h:120
BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, std::shared_ptr< BasisExtents > extents)
Definition: cubature.cc:3871
double * z_
Spherical nodes, on the unit sphere.
Definition: cubature.h:290
void set_debug(int debug)
Definition: cubature.h:176
int npoints() const
Number of grid points.
Definition: cubature.h:154
short namedGrid
Definition: cubature.h:121
double maxR() const
Maximum spatial extent over all atoms.
Definition: cubature.h:455
double R_
Bounding radius of the BlockOPoints.
Definition: cubature.h:384
int max_points_
Maximum number of points in a block.
Definition: cubature.h:76
std::shared_ptr< Vector > shell_extents_
Significant extent of shells.
Definition: cubature.h:433
void build_angles()
Build the spherical angles from &lt;x,y,z&gt;, for reference.
Definition: cubature.cc:4700
Definition: cubature.h:275
double * y() const
The y points. You do not own this.
Definition: cubature.h:165
double * z_
Full z points.
Definition: cubature.h:86
double * x_
Spherical nodes, on the unit sphere.
Definition: cubature.h:286
PseudospectralGrid(std::shared_ptr< Molecule > molecule, std::shared_ptr< BasisSet > primary, Options &options)
Constructor to use for autogeneration.
Definition: cubature.cc:4074
void buildGridFromOptions(std::map< std::string, int > int_opts_map, std::map< std::string, std::string > opts_map)
Master builder methods.
Definition: cubature.cc:4021
static void initialize_lebedev()
Initialize the above arrays with the unique Lebedev grids.
Definition: cubature.cc:4735
double * w() const
Radial weights (including alpha scale and r^2). You do not own this.
Definition: cubature.h:269
int * index() const
index_[fast_index] = slow_index. You do not own this
Definition: cubature.h:151
void print(std::string out_fname="outfile", int print=2) const
Print information about the grid.
Definition: cubature.cc:4222
short radscheme
Definition: cubature.h:118
SharedVector wvec_
Definition: cubature.h:364
std::shared_ptr< Matrix > orientation_
Orientation matrix.
Definition: cubature.h:93
double alpha() const
Alpha scale (for user reference)
Definition: cubature.h:265
int debug_
Definition: cubature.h:66
double * y_
Pointer to y (does not own)
Definition: cubature.h:369
double * x() const
Spherical nodes, on the unit sphere.
Definition: cubature.h:333
int npoints_
Number of points in radial grid.
Definition: cubature.h:282
static std::map< int, int > lebedev_mapping_
Grid npoints to order map.
Definition: cubature.h:302
int npoints() const
Number of points in radial grid.
Definition: cubature.h:263
void bound()
Compute bounding sphere.
Definition: cubature.cc:3894
double * r_
Nodes (including alpha)
Definition: cubature.h:233
virtual ~RadialGrid()
Destructor.
Definition: cubature.cc:4570
double delta() const
The cutoff value.
Definition: cubature.h:449
std::vector< int > functions_local_to_global_
Relevant functions, local -&gt; global.
Definition: cubature.h:377
std::shared_ptr< BasisExtents > extents() const
Pointer to basis extents.
Definition: cubature.h:172
std::shared_ptr< Vector > SharedVector
Definition: adc.h:51
Vector3 xc_
Center of this BlockOPoints.
Definition: cubature.h:382
double w
Definition: cubature.h:60
double maxR_
Maximum extent.
Definition: cubature.h:435
std::vector< std::vector< std::shared_ptr< SphericalGrid > > > spherical_grids_
Spherical grids, per atom and radial point.
Definition: cubature.h:97
size_t index_
number of points in this block
Definition: cubature.h:356
double * x_
Full x points.
Definition: cubature.h:82