MPQC  2.3.1
integrator.h
1 //
2 // integrator.h --- definition of the electron density integrator
3 //
4 // Copyright (C) 1997 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _chemistry_qc_dft_integrator_h
29 #define _chemistry_qc_dft_integrator_h
30 
31 #ifdef __GNUC__
32 #pragma interface
33 #endif
34 
35 #include <util/state/state.h>
36 #include <util/group/thread.h>
37 #include <chemistry/qc/dft/functional.h>
38 #include <chemistry/qc/basis/extent.h>
39 #include <chemistry/qc/wfn/density.h>
40 
41 namespace sc {
42 
44 class DenIntegrator: virtual public SavableState {
45  protected:
46  Ref<Wavefunction> wfn_;
47 //clj Ref<ShellExtent> extent_;
49 
50  Ref<ThreadGrp> threadgrp_;
51  Ref<MessageGrp> messagegrp_;
52 
53  double value_;
54  double accuracy_;
55 
56  double *alpha_vmat_;
57  double *beta_vmat_;
58 
59 //clj double *alpha_dmat_;
60 //clj double *beta_dmat_;
61 //clj double *dmat_bound_;
62 
63  int spin_polarized_;
64 
65  int need_density_; // specialization must set to 1 if it needs density_
66  double density_;
67  int nbasis_;
68  int nshell_;
69  int n_integration_center_;
70  int natom_;
71  int compute_potential_integrals_; // 1 if potential integrals are needed
72 
73  int linear_scaling_;
74  int use_dmat_bound_;
75 
76  void init_integration(const Ref<DenFunctional> &func,
77  const RefSymmSCMatrix& densa,
78  const RefSymmSCMatrix& densb,
79  double *nuclear_gradient);
80  void done_integration();
81 
82  void init_object();
83  public:
90  ~DenIntegrator();
92 
94  Ref<Wavefunction> wavefunction() const { return wfn_; }
96  double value() const { return value_; }
97 
99  void set_accuracy(double a);
100  double get_accuracy(void) {return accuracy_; }
106  const double *alpha_vmat() const { return alpha_vmat_; }
109  const double *beta_vmat() const { return beta_vmat_; }
110 
113  virtual void init(const Ref<Wavefunction> &);
115  virtual void done();
119  virtual void integrate(const Ref<DenFunctional> &,
120  const RefSymmSCMatrix& densa =0,
121  const RefSymmSCMatrix& densb =0,
122  double *nuclear_grad = 0) = 0;
123 };
124 
125 
127 class IntegrationWeight: virtual public SavableState {
128 
129  protected:
130 
131  Ref<Molecule> mol_;
132  double tol_;
133 
134  void fd_w(int icenter, SCVector3 &point, double *fd_grad_w);
135 
136  public:
142 
143  void test(int icenter, SCVector3 &point);
144  void test();
145 
147  virtual void init(const Ref<Molecule> &, double tolerance);
149  virtual void done();
154  virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0;
155 };
156 
157 
160 
161  int n_integration_centers;
162  SCVector3 *centers;
163  double *atomic_radius;
164 
165  double **a_mat;
166  double **oorab;
167 
168  void compute_grad_p(int gc, int ic, int wc, SCVector3&r, double p,
169  SCVector3&g);
170  void compute_grad_nu(int gc, int bc, SCVector3 &point, SCVector3 &grad);
171 
172  double compute_t(int ic, int jc, SCVector3 &point);
173  double compute_p(int icenter, SCVector3&point);
174 
175  public:
181 
182  void init(const Ref<Molecule> &, double tolerance);
183  void done();
184  double w(int center, SCVector3 &point, double *grad_w = 0);
185 };
186 
188 class RadialIntegrator: virtual public SavableState {
189  protected:
190  int nr_;
191  public:
193  RadialIntegrator(const Ref<KeyVal> &);
195  ~RadialIntegrator();
197 
198  virtual int nr() const = 0;
199  virtual double radial_value(int ir, int nr, double radii,
200  double &multiplier) = 0;
201 };
202 
203 
205 class AngularIntegrator: virtual public SavableState{
206  protected:
207  public:
213 
214  virtual int nw(void) const = 0;
215  virtual int num_angular_points(double r_value, int ir) = 0;
216  virtual double angular_point_cartesian(int iangular, double r,
217  SCVector3 &integration_point) const = 0;
218 };
219 
220 
224  public:
234 
235  int nr() const;
236  double radial_value(int ir, int nr, double radii, double &multiplier);
237 
238  void print(std::ostream & =ExEnv::out0()) const;
239 };
240 
283  protected:
284  int npoint_;
285  double *x_, *y_, *z_, *w_;
286 
287  void init(int n);
288  public:
300 
301  int nw(void) const;
302  int num_angular_points(double r_value, int ir);
303  double angular_point_cartesian(int iangular, double r,
304  SCVector3 &integration_point) const;
305  void print(std::ostream & =ExEnv::out0()) const;
306 };
307 
311  protected:
312  int ntheta_;
313  int nphi_;
314  int Ktheta_;
315  int ntheta_r_;
316  int nphi_r_;
317  int Ktheta_r_;
318  double *theta_quad_weights_;
319  double *theta_quad_points_;
320 
321  int get_ntheta(void) const;
322  void set_ntheta(int i);
323  int get_nphi(void) const;
324  void set_nphi(int i);
325  int get_Ktheta(void) const;
326  void set_Ktheta(int i);
327  int get_ntheta_r(void) const;
328  void set_ntheta_r(int i);
329  int get_nphi_r(void) const;
330  void set_nphi_r(int i);
331  int get_Ktheta_r(void) const;
332  void set_Ktheta_r(int i);
333  int nw(void) const;
334  double sin_theta(SCVector3 &point) const;
335  void gauleg(double x1, double x2, int n);
336  public:
348 
349  int num_angular_points(double r_value, int ir);
350  double angular_point_cartesian(int iangular, double r,
351  SCVector3 &integration_point) const;
352  void print(std::ostream & =ExEnv::out0()) const;
353 };
354 
358  private:
359  int prune_grid_;
360  double **Alpha_coeffs_;
361  int gridtype_;
362  int **nr_points_, *xcoarse_l_;
363  int npruned_partitions_;
364  double *grid_accuracy_;
365  int dynamic_grids_;
366  int natomic_rows_;
367  int max_gridtype_;
368  protected:
369  Ref<IntegrationWeight> weight_;
370  Ref<RadialIntegrator> radial_user_;
371  Ref<AngularIntegrator> angular_user_;
372  Ref<AngularIntegrator> ***angular_grid_;
373  Ref<RadialIntegrator> **radial_grid_;
374  public:
420 
422  const RefSymmSCMatrix& densa =0,
423  const RefSymmSCMatrix& densb =0,
424  double *nuclear_gradient = 0);
425  void print(std::ostream & =ExEnv::out0()) const;
426  AngularIntegrator *get_angular_grid(double radius, double atomic_radius,
427  int charge, int deriv_order);
428  RadialIntegrator *get_radial_grid(int charge, int deriv_order);
429  void init_default_grids(void);
430  int angular_grid_offset(int i);
431  void set_grids(void);
432  int get_atomic_row(int i);
433  void init_parameters(void);
434  void init_parameters(const Ref<KeyVal>& keyval);
435  void init_pruning_coefficients(const Ref<KeyVal>& keyval);
436  void init_pruning_coefficients(void);
437  void init_alpha_coefficients(void);
438  int select_dynamic_grid(void);
439  Ref<IntegrationWeight> weight() { return weight_; }
440 };
441 
442 }
443 
444 #endif
445 
446 // Local Variables:
447 // mode: c++
448 // c-file-style: "CLJ"
449 // End:
An abstract base class for angular integrators.
Definition: integrator.h:205
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
Implements Becke's integration weight scheme.
Definition: integrator.h:159
double w(int center, SCVector3 &point, double *grad_w=0)
Computes the weight for a given center at a given point in space.
void init(const Ref< Molecule > &, double tolerance)
Initialize the integration weight object.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
void done()
Called when finished with the integration weight object.
An abstract base class for integrating the electron density.
Definition: integrator.h:44
void set_accuracy(double a)
Sets the accuracy to use in the integration.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
DenIntegrator(const Ref< KeyVal > &)
Construct a new DenIntegrator given the KeyVal input.
virtual void init(const Ref< Wavefunction > &)
Called before integrate.
virtual void integrate(const Ref< DenFunctional > &, const RefSymmSCMatrix &densa=0, const RefSymmSCMatrix &densb=0, double *nuclear_grad=0)=0
Performs the integration of the given functional using the given alpha and beta density matrices.
DenIntegrator(StateIn &)
Construct a new DenIntegrator given the StateIn data.
DenIntegrator()
Construct a new DenIntegrator.
Ref< Wavefunction > wavefunction() const
Returns the wavefunction used for the integration.
Definition: integrator.h:94
const double * beta_vmat() const
Returns the beta potential integrals.
Definition: integrator.h:109
virtual void done()
Must be called between calls to init.
void set_compute_potential_integrals(int)
Call with non zero if the potential integrals are to be computed.
const double * alpha_vmat() const
Returns the alpha potential integrals.
Definition: integrator.h:106
double value() const
Returns the result of the integration.
Definition: integrator.h:96
An implementation of a radial integrator using the Euler-Maclaurin weights and grid points.
Definition: integrator.h:223
void print(std::ostream &=ExEnv::out0()) const
Print the object.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
EulerMaclaurinRadialIntegrator(const Ref< KeyVal > &)
Constructs a EulerMaclaurinRadialIntegrator from KeyVal input.
static std::ostream & out0()
Return an ostream that writes from node 0.
An implementation of an angular integrator using the Gauss-Legendre weights and grid points.
Definition: integrator.h:310
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
void print(std::ostream &=ExEnv::out0()) const
Print the object.
GaussLegendreAngularIntegrator(const Ref< KeyVal > &)
Contract a GaussLegendreAngularIntegrator from KeyVal input.
An abstract base class for computing grid weights.
Definition: integrator.h:127
virtual double w(int center, SCVector3 &point, double *grad_w=0)=0
Computes the weight for a given center at a given point in space.
virtual void done()
Called when finished with the integration weight object.
virtual void init(const Ref< Molecule > &, double tolerance)
Initialize the integration weight object.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
An implementation of a Lebedev angular integrator.
Definition: integrator.h:282
LebedevLaikovIntegrator(const Ref< KeyVal > &)
Construct a LebedevLaikovIntegrator using the given KeyVal input.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
void print(std::ostream &=ExEnv::out0()) const
Print the object.
An implementation of an integrator using any combination of a RadialIntegrator and an AngularIntegrat...
Definition: integrator.h:357
void integrate(const Ref< DenFunctional > &, const RefSymmSCMatrix &densa=0, const RefSymmSCMatrix &densb=0, double *nuclear_gradient=0)
Performs the integration of the given functional using the given alpha and beta density matrices.
void print(std::ostream &=ExEnv::out0()) const
Print the object.
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
RadialAngularIntegrator(const Ref< KeyVal > &)
Construct a RadialAngularIntegrator from KeyVal input.
An abstract base class for radial integrators.
Definition: integrator.h:188
void save_data_state(StateOut &)
Save the base classes (with save_data_state) and the members in the same order that the StateIn CTOR ...
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:261
A template class that maintains references counts.
Definition: ref.h:332
Definition: vector3.h:46
Base class for objects that can save/restore state.
Definition: state.h:46
Restores objects that derive from SavableState.
Definition: statein.h:70
Serializes objects that derive from SavableState.
Definition: stateout.h:61
Definition: implicit.h:5

Generated at Thu Jan 20 2022 00:00:00 for MPQC 2.3.1 using the documentation package Doxygen 1.9.1.