# 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 #include <chemistry/qc/wfn/density.h>
    00040
    00041 namespace sc {
    00042
    00044 class DenIntegrator: virtual public SavableState {
    00045   protected:
    00046     Ref<Wavefunction> wfn_;
    00047 //clj    Ref<ShellExtent> extent_;
    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 //clj    double *alpha_dmat_;
    00060 //clj    double *beta_dmat_;
    00061 //clj    double *dmat_bound_;
    00062
    00063     int spin_polarized_;
    00064
    00065     int need_density_; // specialization must set to 1 if it needs density_
    00066     double density_;
    00067     int nbasis_;
    00068     int nshell_;
    00069     int n_integration_center_;
    00070     int natom_;
    00071     int compute_potential_integrals_; // 1 if potential integrals are needed
    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 // Local Variables:
    00447 // mode: c++
    00448 // c-file-style: "CLJ"
    00449 // End: