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