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 #include <chemistry/qc/wfn/density.h>
00040
00041 namespace sc {
00042
00044 class DenIntegrator: virtual public SavableState {
00045 protected:
00046 Ref<Wavefunction> wfn_;
00047
00048 Ref<BatchElectronDensity> den_;
00049
00050 Ref<ThreadGrp> threadgrp_;
00051 Ref<MessageGrp> messagegrp_;
00052
00053 double value_;
00054 double accuracy_;
00055
00056 double *alpha_vmat_;
00057 double *beta_vmat_;
00058
00059
00060
00061
00062
00063 int spin_polarized_;
00064
00065 int need_density_;
00066 double density_;
00067 int nbasis_;
00068 int nshell_;
00069 int n_integration_center_;
00070 int natom_;
00071 int compute_potential_integrals_;
00072
00073 int linear_scaling_;
00074 int use_dmat_bound_;
00075
00076 void init_integration(const Ref<DenFunctional> &func,
00077 const RefSymmSCMatrix& densa,
00078 const RefSymmSCMatrix& densb,
00079 double *nuclear_gradient);
00080 void done_integration();
00081
00082 void init_object();
00083 public:
00085 DenIntegrator();
00087 DenIntegrator(const Ref<KeyVal> &);
00089 DenIntegrator(StateIn &);
00090 ~DenIntegrator();
00091 void save_data_state(StateOut &);
00092
00094 Ref<Wavefunction> wavefunction() const { return wfn_; }
00096 double value() const { return value_; }
00097
00099 void set_accuracy(double a);
00100 double get_accuracy(void) {return accuracy_; }
00103 void set_compute_potential_integrals(int);
00106 const double *alpha_vmat() const { return alpha_vmat_; }
00109 const double *beta_vmat() const { return beta_vmat_; }
00110
00113 virtual void init(const Ref<Wavefunction> &);
00115 virtual void done();
00119 virtual void integrate(const Ref<DenFunctional> &,
00120 const RefSymmSCMatrix& densa =0,
00121 const RefSymmSCMatrix& densb =0,
00122 double *nuclear_grad = 0) = 0;
00123 };
00124
00125
00127 class IntegrationWeight: virtual public SavableState {
00128
00129 protected:
00130
00131 Ref<Molecule> mol_;
00132 double tol_;
00133
00134 void fd_w(int icenter, SCVector3 &point, double *fd_grad_w);
00135
00136 public:
00137 IntegrationWeight();
00138 IntegrationWeight(const Ref<KeyVal> &);
00139 IntegrationWeight(StateIn &);
00140 ~IntegrationWeight();
00141 void save_data_state(StateOut &);
00142
00143 void test(int icenter, SCVector3 &point);
00144 void test();
00145
00147 virtual void init(const Ref<Molecule> &, double tolerance);
00149 virtual void done();
00154 virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0;
00155 };
00156
00157
00159 class BeckeIntegrationWeight: public IntegrationWeight {
00160
00161 int n_integration_centers;
00162 SCVector3 *centers;
00163 double *atomic_radius;
00164
00165 double **a_mat;
00166 double **oorab;
00167
00168 void compute_grad_p(int gc, int ic, int wc, SCVector3&r, double p,
00169 SCVector3&g);
00170 void compute_grad_nu(int gc, int bc, SCVector3 &point, SCVector3 &grad);
00171
00172 double compute_t(int ic, int jc, SCVector3 &point);
00173 double compute_p(int icenter, SCVector3&point);
00174
00175 public:
00176 BeckeIntegrationWeight();
00177 BeckeIntegrationWeight(const Ref<KeyVal> &);
00178 BeckeIntegrationWeight(StateIn &);
00179 ~BeckeIntegrationWeight();
00180 void save_data_state(StateOut &);
00181
00182 void init(const Ref<Molecule> &, double tolerance);
00183 void done();
00184 double w(int center, SCVector3 &point, double *grad_w = 0);
00185 };
00186
00188 class RadialIntegrator: virtual public SavableState {
00189 protected:
00190 int nr_;
00191 public:
00192 RadialIntegrator();
00193 RadialIntegrator(const Ref<KeyVal> &);
00194 RadialIntegrator(StateIn &);
00195 ~RadialIntegrator();
00196 void save_data_state(StateOut &);
00197
00198 virtual int nr() const = 0;
00199 virtual double radial_value(int ir, int nr, double radii,
00200 double &multiplier) = 0;
00201 };
00202
00203
00205 class AngularIntegrator: virtual public SavableState{
00206 protected:
00207 public:
00208 AngularIntegrator();
00209 AngularIntegrator(const Ref<KeyVal> &);
00210 AngularIntegrator(StateIn &);
00211 ~AngularIntegrator();
00212 void save_data_state(StateOut &);
00213
00214 virtual int nw(void) const = 0;
00215 virtual int num_angular_points(double r_value, int ir) = 0;
00216 virtual double angular_point_cartesian(int iangular, double r,
00217 SCVector3 &integration_point) const = 0;
00218 };
00219
00220
00223 class EulerMaclaurinRadialIntegrator: public RadialIntegrator {
00224 public:
00225 EulerMaclaurinRadialIntegrator();
00226 EulerMaclaurinRadialIntegrator(int i);
00230 EulerMaclaurinRadialIntegrator(const Ref<KeyVal> &);
00231 EulerMaclaurinRadialIntegrator(StateIn &);
00232 ~EulerMaclaurinRadialIntegrator();
00233 void save_data_state(StateOut &);
00234
00235 int nr() const;
00236 double radial_value(int ir, int nr, double radii, double &multiplier);
00237
00238 void print(std::ostream & =ExEnv::out0()) const;
00239 };
00240
00282 class LebedevLaikovIntegrator: public AngularIntegrator {
00283 protected:
00284 int npoint_;
00285 double *x_, *y_, *z_, *w_;
00286
00287 void init(int n);
00288 public:
00289 LebedevLaikovIntegrator();
00295 LebedevLaikovIntegrator(const Ref<KeyVal> &);
00296 LebedevLaikovIntegrator(StateIn &);
00297 LebedevLaikovIntegrator(int);
00298 ~LebedevLaikovIntegrator();
00299 void save_data_state(StateOut &);
00300
00301 int nw(void) const;
00302 int num_angular_points(double r_value, int ir);
00303 double angular_point_cartesian(int iangular, double r,
00304 SCVector3 &integration_point) const;
00305 void print(std::ostream & =ExEnv::out0()) const;
00306 };
00307
00310 class GaussLegendreAngularIntegrator: public AngularIntegrator {
00311 protected:
00312 int ntheta_;
00313 int nphi_;
00314 int Ktheta_;
00315 int ntheta_r_;
00316 int nphi_r_;
00317 int Ktheta_r_;
00318 double *theta_quad_weights_;
00319 double *theta_quad_points_;
00320
00321 int get_ntheta(void) const;
00322 void set_ntheta(int i);
00323 int get_nphi(void) const;
00324 void set_nphi(int i);
00325 int get_Ktheta(void) const;
00326 void set_Ktheta(int i);
00327 int get_ntheta_r(void) const;
00328 void set_ntheta_r(int i);
00329 int get_nphi_r(void) const;
00330 void set_nphi_r(int i);
00331 int get_Ktheta_r(void) const;
00332 void set_Ktheta_r(int i);
00333 int nw(void) const;
00334 double sin_theta(SCVector3 &point) const;
00335 void gauleg(double x1, double x2, int n);
00336 public:
00337 GaussLegendreAngularIntegrator();
00344 GaussLegendreAngularIntegrator(const Ref<KeyVal> &);
00345 GaussLegendreAngularIntegrator(StateIn &);
00346 ~GaussLegendreAngularIntegrator();
00347 void save_data_state(StateOut &);
00348
00349 int num_angular_points(double r_value, int ir);
00350 double angular_point_cartesian(int iangular, double r,
00351 SCVector3 &integration_point) const;
00352 void print(std::ostream & =ExEnv::out0()) const;
00353 };
00354
00357 class RadialAngularIntegrator: public DenIntegrator {
00358 private:
00359 int prune_grid_;
00360 double **Alpha_coeffs_;
00361 int gridtype_;
00362 int **nr_points_, *xcoarse_l_;
00363 int npruned_partitions_;
00364 double *grid_accuracy_;
00365 int dynamic_grids_;
00366 int natomic_rows_;
00367 int max_gridtype_;
00368 protected:
00369 Ref<IntegrationWeight> weight_;
00370 Ref<RadialIntegrator> radial_user_;
00371 Ref<AngularIntegrator> angular_user_;
00372 Ref<AngularIntegrator> ***angular_grid_;
00373 Ref<RadialIntegrator> **radial_grid_;
00374 public:
00375 RadialAngularIntegrator();
00416 RadialAngularIntegrator(const Ref<KeyVal> &);
00417 RadialAngularIntegrator(StateIn &);
00418 ~RadialAngularIntegrator();
00419 void save_data_state(StateOut &);
00420
00421 void integrate(const Ref<DenFunctional> &,
00422 const RefSymmSCMatrix& densa =0,
00423 const RefSymmSCMatrix& densb =0,
00424 double *nuclear_gradient = 0);
00425 void print(std::ostream & =ExEnv::out0()) const;
00426 AngularIntegrator *get_angular_grid(double radius, double atomic_radius,
00427 int charge, int deriv_order);
00428 RadialIntegrator *get_radial_grid(int charge, int deriv_order);
00429 void init_default_grids(void);
00430 int angular_grid_offset(int i);
00431 void set_grids(void);
00432 int get_atomic_row(int i);
00433 void init_parameters(void);
00434 void init_parameters(const Ref<KeyVal>& keyval);
00435 void init_pruning_coefficients(const Ref<KeyVal>& keyval);
00436 void init_pruning_coefficients(void);
00437 void init_alpha_coefficients(void);
00438 int select_dynamic_grid(void);
00439 Ref<IntegrationWeight> weight() { return weight_; }
00440 };
00441
00442 }
00443
00444 #endif
00445
00446
00447
00448
00449