MPQC 2.3.1

- Main Page
- Related Pages
- Classes
- Files

- File List

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_;
48 Ref<BatchElectronDensity> den_;
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:
85 DenIntegrator();
87 DenIntegrator(const Ref<KeyVal> &);
89 DenIntegrator(StateIn &);
90 ~DenIntegrator();
91 void save_data_state(StateOut &);
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_; }
103 void set_compute_potential_integrals(int);
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:
137 IntegrationWeight();
138 IntegrationWeight(const Ref<KeyVal> &);
139 IntegrationWeight(StateIn &);
140 ~IntegrationWeight();
141 void save_data_state(StateOut &);
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
159 class 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();
180 void save_data_state(StateOut &);
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:
192 RadialIntegrator();
193 RadialIntegrator(const Ref<KeyVal> &);
194 RadialIntegrator(StateIn &);
195 ~RadialIntegrator();
196 void save_data_state(StateOut &);
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:
208 AngularIntegrator();
209 AngularIntegrator(const Ref<KeyVal> &);
210 AngularIntegrator(StateIn &);
211 ~AngularIntegrator();
212 void save_data_state(StateOut &);
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
223 class EulerMaclaurinRadialIntegrator: public RadialIntegrator {
224 public:
225 EulerMaclaurinRadialIntegrator();
226 EulerMaclaurinRadialIntegrator(int i);
230 EulerMaclaurinRadialIntegrator(const Ref<KeyVal> &);
231 EulerMaclaurinRadialIntegrator(StateIn &);
232 ~EulerMaclaurinRadialIntegrator();
233 void save_data_state(StateOut &);
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
282 class LebedevLaikovIntegrator: public AngularIntegrator {
283 protected:
284 int npoint_;
285 double *x_, *y_, *z_, *w_;
286
287 void init(int n);
288 public:
289 LebedevLaikovIntegrator();
295 LebedevLaikovIntegrator(const Ref<KeyVal> &);
296 LebedevLaikovIntegrator(StateIn &);
297 LebedevLaikovIntegrator(int);
298 ~LebedevLaikovIntegrator();
299 void save_data_state(StateOut &);
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
310 class 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();
344 GaussLegendreAngularIntegrator(const Ref<KeyVal> &);
345 GaussLegendreAngularIntegrator(StateIn &);
346 ~GaussLegendreAngularIntegrator();
347 void save_data_state(StateOut &);
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
357 class 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:
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:
375 RadialAngularIntegrator();
416 RadialAngularIntegrator(const Ref<KeyVal> &);
417 RadialAngularIntegrator(StateIn &);
418 ~RadialAngularIntegrator();
419 void save_data_state(StateOut &);
420
421 void integrate(const Ref<DenFunctional> &,
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: