Google

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:

Generated at Fri Jan 10 08:14:09 2003 for MPQC 2.1.3 using the documentation package Doxygen 1.2.14.