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:

Generated at Mon Dec 3 23:23:37 2007 for MPQC 2.3.1 using the documentation package Doxygen 1.5.2.