00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef _chemistry_qc_dft_integrator_h
00029 #define _chemistry_qc_dft_integrator_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #include <util/state/state.h>
00036 #include <util/group/thread.h>
00037 #include <chemistry/qc/dft/functional.h>
00038 #include <chemistry/qc/basis/extent.h>
00039
00040 namespace sc {
00041
00043 class DenIntegrator: virtual public SavableState {
00044 protected:
00045 Ref<Wavefunction> wfn_;
00046 Ref<ShellExtent> extent_;
00047
00048 Ref<ThreadGrp> threadgrp_;
00049 Ref<MessageGrp> messagegrp_;
00050
00051 double value_;
00052 double accuracy_;
00053
00054 double *alpha_vmat_;
00055 double *beta_vmat_;
00056
00057 double *alpha_dmat_;
00058 double *beta_dmat_;
00059 double *dmat_bound_;
00060
00061 int spin_polarized_;
00062
00063 int need_density_;
00064 double density_;
00065 int nbasis_;
00066 int nshell_;
00067 int natom_;
00068 int compute_potential_integrals_;
00069
00070 int linear_scaling_;
00071 int use_dmat_bound_;
00072
00073 void init_integration(const Ref<DenFunctional> &func,
00074 const RefSymmSCMatrix& densa,
00075 const RefSymmSCMatrix& densb,
00076 double *nuclear_gradient);
00077 void done_integration();
00078
00079 void init_object();
00080 public:
00082 DenIntegrator();
00084 DenIntegrator(const Ref<KeyVal> &);
00086 DenIntegrator(StateIn &);
00087 ~DenIntegrator();
00088 void save_data_state(StateOut &);
00089
00091 Ref<Wavefunction> wavefunction() const { return wfn_; }
00093 double value() const { return value_; }
00094
00096 void set_accuracy(double a) { accuracy_ = a; }
00097 double get_accuracy(void) {return accuracy_; }
00100 void set_compute_potential_integrals(int);
00103 const double *alpha_vmat() const { return alpha_vmat_; }
00106 const double *beta_vmat() const { return beta_vmat_; }
00107
00110 virtual void init(const Ref<Wavefunction> &);
00112 virtual void done();
00116 virtual void integrate(const Ref<DenFunctional> &,
00117 const RefSymmSCMatrix& densa =0,
00118 const RefSymmSCMatrix& densb =0,
00119 double *nuclear_grad = 0) = 0;
00120 };
00121
00122
00124 class IntegrationWeight: virtual public SavableState {
00125
00126 protected:
00127
00128 Ref<Molecule> mol_;
00129 double tol_;
00130
00131 void fd_w(int icenter, SCVector3 &point, double *fd_grad_w);
00132
00133 public:
00134 IntegrationWeight();
00135 IntegrationWeight(const Ref<KeyVal> &);
00136 IntegrationWeight(StateIn &);
00137 ~IntegrationWeight();
00138 void save_data_state(StateOut &);
00139
00140 void test(int icenter, SCVector3 &point);
00141 void test();
00142
00144 virtual void init(const Ref<Molecule> &, double tolerance);
00146 virtual void done();
00151 virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0;
00152 };
00153
00154
00156 class BeckeIntegrationWeight: public IntegrationWeight {
00157
00158 int ncenters;
00159 SCVector3 *centers;
00160 double *atomic_radius;
00161
00162 double **a_mat;
00163 double **oorab;
00164
00165 void compute_grad_p(int gc, int ic, int wc, SCVector3&r, double p,
00166 SCVector3&g);
00167 void compute_grad_nu(int gc, int bc, SCVector3 &point, SCVector3 &grad);
00168
00169 double compute_t(int ic, int jc, SCVector3 &point);
00170 double compute_p(int icenter, SCVector3&point);
00171
00172 public:
00173 BeckeIntegrationWeight();
00174 BeckeIntegrationWeight(const Ref<KeyVal> &);
00175 BeckeIntegrationWeight(StateIn &);
00176 ~BeckeIntegrationWeight();
00177 void save_data_state(StateOut &);
00178
00179 void init(const Ref<Molecule> &, double tolerance);
00180 void done();
00181 double w(int center, SCVector3 &point, double *grad_w = 0);
00182 };
00183
00185 class RadialIntegrator: virtual public SavableState {
00186 protected:
00187 int nr_;
00188 public:
00189 RadialIntegrator();
00190 RadialIntegrator(const Ref<KeyVal> &);
00191 RadialIntegrator(StateIn &);
00192 ~RadialIntegrator();
00193 void save_data_state(StateOut &);
00194
00195 virtual int nr() const = 0;
00196 virtual double radial_value(int ir, int nr, double radii,
00197 double &multiplier) = 0;
00198 };
00199
00200
00202 class AngularIntegrator: virtual public SavableState{
00203 protected:
00204 public:
00205 AngularIntegrator();
00206 AngularIntegrator(const Ref<KeyVal> &);
00207 AngularIntegrator(StateIn &);
00208 ~AngularIntegrator();
00209 void save_data_state(StateOut &);
00210
00211 virtual int nw(void) const = 0;
00212 virtual int num_angular_points(double r_value, int ir) = 0;
00213 virtual double angular_point_cartesian(int iangular, double r,
00214 SCVector3 &integration_point) const = 0;
00215 };
00216
00217
00220 class EulerMaclaurinRadialIntegrator: public RadialIntegrator {
00221 public:
00222 EulerMaclaurinRadialIntegrator();
00223 EulerMaclaurinRadialIntegrator(int i);
00224 EulerMaclaurinRadialIntegrator(const Ref<KeyVal> &);
00225 EulerMaclaurinRadialIntegrator(StateIn &);
00226 ~EulerMaclaurinRadialIntegrator();
00227 void save_data_state(StateOut &);
00228
00229 int nr() const;
00230 double radial_value(int ir, int nr, double radii, double &multiplier);
00231
00232 void print(std::ostream & =ExEnv::out0()) const;
00233 };
00234
00276 class LebedevLaikovIntegrator: public AngularIntegrator {
00277 protected:
00278 int npoint_;
00279 double *x_, *y_, *z_, *w_;
00280
00281 void init(int n);
00282 public:
00283 LebedevLaikovIntegrator();
00284 LebedevLaikovIntegrator(const Ref<KeyVal> &);
00285 LebedevLaikovIntegrator(StateIn &);
00286 LebedevLaikovIntegrator(int);
00287 ~LebedevLaikovIntegrator();
00288 void save_data_state(StateOut &);
00289
00290 int nw(void) const;
00291 int num_angular_points(double r_value, int ir);
00292 double angular_point_cartesian(int iangular, double r,
00293 SCVector3 &integration_point) const;
00294 void print(std::ostream & =ExEnv::out0()) const;
00295 };
00296
00299 class GaussLegendreAngularIntegrator: public AngularIntegrator {
00300 protected:
00301 int ntheta_;
00302 int nphi_;
00303 int Ktheta_;
00304 int ntheta_r_;
00305 int nphi_r_;
00306 int Ktheta_r_;
00307 double *theta_quad_weights_;
00308 double *theta_quad_points_;
00309
00310 int get_ntheta(void) const;
00311 void set_ntheta(int i);
00312 int get_nphi(void) const;
00313 void set_nphi(int i);
00314 int get_Ktheta(void) const;
00315 void set_Ktheta(int i);
00316 int get_ntheta_r(void) const;
00317 void set_ntheta_r(int i);
00318 int get_nphi_r(void) const;
00319 void set_nphi_r(int i);
00320 int get_Ktheta_r(void) const;
00321 void set_Ktheta_r(int i);
00322 int nw(void) const;
00323 double sin_theta(SCVector3 &point) const;
00324 void gauleg(double x1, double x2, int n);
00325 public:
00326 GaussLegendreAngularIntegrator();
00327 GaussLegendreAngularIntegrator(const Ref<KeyVal> &);
00328 GaussLegendreAngularIntegrator(StateIn &);
00329 ~GaussLegendreAngularIntegrator();
00330 void save_data_state(StateOut &);
00331
00332 int num_angular_points(double r_value, int ir);
00333 double angular_point_cartesian(int iangular, double r,
00334 SCVector3 &integration_point) const;
00335 void print(std::ostream & =ExEnv::out0()) const;
00336 };
00337
00340 class RadialAngularIntegrator: public DenIntegrator {
00341 private:
00342 int prune_grid_;
00343 double **Alpha_coeffs_;
00344 int gridtype_;
00345 int **nr_points_, *xcoarse_l_;
00346 int npruned_partitions_;
00347 double *grid_accuracy_;
00348 int dynamic_grids_;
00349 int natomic_rows_;
00350 int max_gridtype_;
00351 protected:
00352 Ref<IntegrationWeight> weight_;
00353 Ref<RadialIntegrator> radial_user_;
00354 Ref<AngularIntegrator> angular_user_;
00355 Ref<AngularIntegrator> ***angular_grid_;
00356 Ref<RadialIntegrator> **radial_grid_;
00357 public:
00358 RadialAngularIntegrator();
00359 RadialAngularIntegrator(const Ref<KeyVal> &);
00360 RadialAngularIntegrator(StateIn &);
00361 ~RadialAngularIntegrator();
00362 void save_data_state(StateOut &);
00363
00364 void integrate(const Ref<DenFunctional> &,
00365 const RefSymmSCMatrix& densa =0,
00366 const RefSymmSCMatrix& densb =0,
00367 double *nuclear_gradient = 0);
00368 void print(std::ostream & =ExEnv::out0()) const;
00369 AngularIntegrator *get_angular_grid(double radius, double atomic_radius,
00370 int charge, int deriv_order);
00371 RadialIntegrator *get_radial_grid(int charge, int deriv_order);
00372 void init_default_grids(void);
00373 int angular_grid_offset(int i);
00374 void set_grids(void);
00375 int get_atomic_row(int i);
00376 void init_parameters(void);
00377 void init_parameters(const Ref<KeyVal>& keyval);
00378 void init_pruning_coefficients(const Ref<KeyVal>& keyval);
00379 void init_pruning_coefficients(void);
00380 void init_alpha_coefficients(void);
00381 int select_dynamic_grid(void);
00382 Ref<IntegrationWeight> weight() { return weight_; }
00383 };
00384
00385 }
00386
00387 #endif
00388
00389
00390
00391
00392