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
41namespace sc {
42
44class DenIntegrator: virtual public SavableState {
45 protected:
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:
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
127class 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:
137 IntegrationWeight();
138 IntegrationWeight(const Ref<KeyVal> &);
139 IntegrationWeight(StateIn &);
140 ~IntegrationWeight();
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
159class BeckeIntegrationWeight: public IntegrationWeight {
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:
176 BeckeIntegrationWeight();
177 BeckeIntegrationWeight(const Ref<KeyVal> &);
178 BeckeIntegrationWeight(StateIn &);
179 ~BeckeIntegrationWeight();
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
188class RadialIntegrator: virtual public SavableState {
189 protected:
190 int nr_;
191 public:
192 RadialIntegrator();
193 RadialIntegrator(const Ref<KeyVal> &);
194 RadialIntegrator(StateIn &);
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
205class AngularIntegrator: virtual public SavableState{
206 protected:
207 public:
208 AngularIntegrator();
209 AngularIntegrator(const Ref<KeyVal> &);
210 AngularIntegrator(StateIn &);
211 ~AngularIntegrator();
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
223class EulerMaclaurinRadialIntegrator: public RadialIntegrator {
224 public:
225 EulerMaclaurinRadialIntegrator();
226 EulerMaclaurinRadialIntegrator(int i);
231 EulerMaclaurinRadialIntegrator(StateIn &);
232 ~EulerMaclaurinRadialIntegrator();
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
282class LebedevLaikovIntegrator: public AngularIntegrator {
283 protected:
284 int npoint_;
285 double *x_, *y_, *z_, *w_;
286
287 void init(int n);
288 public:
289 LebedevLaikovIntegrator();
296 LebedevLaikovIntegrator(StateIn &);
297 LebedevLaikovIntegrator(int);
298 ~LebedevLaikovIntegrator();
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
310class GaussLegendreAngularIntegrator: public AngularIntegrator {
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:
337 GaussLegendreAngularIntegrator();
345 GaussLegendreAngularIntegrator(StateIn &);
346 ~GaussLegendreAngularIntegrator();
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
357class RadialAngularIntegrator: public DenIntegrator {
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:
370 Ref<RadialIntegrator> radial_user_;
371 Ref<AngularIntegrator> angular_user_;
372 Ref<AngularIntegrator> ***angular_grid_;
373 Ref<RadialIntegrator> **radial_grid_;
374 public:
375 RadialAngularIntegrator();
417 RadialAngularIntegrator(StateIn &);
418 ~RadialAngularIntegrator();
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 ...
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.
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.
Ref< Wavefunction > wavefunction() const
Returns the wavefunction used for the integration.
Definition integrator.h:94
DenIntegrator(StateIn &)
Construct a new DenIntegrator given the StateIn data.
DenIntegrator()
Construct a new DenIntegrator.
const double * alpha_vmat() const
Returns the alpha potential integrals.
Definition integrator.h:106
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.
double value() const
Returns the result of the integration.
Definition integrator.h:96
const double * beta_vmat() const
Returns the beta potential integrals.
Definition integrator.h:109
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.
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.
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 ...
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.
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
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 Mon Apr 28 2025 00:00:00 for MPQC 2.3.1 using the documentation package Doxygen 1.13.2.