Main Page Class Hierarchy Compound List File List Compound Members Related Pages

# integrator.h

    00001 //
    00002 // integrator.h --- definition of the electron density integrator
    00003 //
    00004 // Copyright (C) 1997 Limit Point Systems, Inc.
    00005 //
    00006 // Author: Curtis Janssen <cljanss@limitpt.com>
    00007 // Maintainer: LPS
    00008 //
    00009 // This file is part of the SC Toolkit.
    00010 //
    00011 // The SC Toolkit is free software; you can redistribute it and/or modify
    00012 // it under the terms of the GNU Library General Public License as published by
    00013 // the Free Software Foundation; either version 2, or (at your option)
    00014 // any later version.
    00015 //
    00016 // The SC Toolkit is distributed in the hope that it will be useful,
    00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
    00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    00019 // GNU Library General Public License for more details.
    00020 //
    00021 // You should have received a copy of the GNU Library General Public License
    00022 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
    00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
    00024 //
    00025 // The U.S. Government is granted a limited license as per AL 91-7.
    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_; // specialization must set to 1 if it needs density_
    00064     double density_;
    00065     int nbasis_;
    00066     int nshell_;
    00067     int natom_;
    00068     int compute_potential_integrals_; // 1 if potential integrals are needed
    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 // Local Variables:
    00390 // mode: c++
    00391 // c-file-style: "CLJ"
    00392 // End: