integrator.h

//
// integrator.h --- definition of the electron density integrator
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
 
#ifndef _chemistry_qc_dft_integrator_h
#define _chemistry_qc_dft_integrator_h
 
#include <util/state/state.h>
#include <util/group/thread.h>
#include <chemistry/qc/dft/functional.h>
#include <chemistry/qc/basis/extent.h>
#include <chemistry/qc/wfn/density.h>
 
namespace sc {
 
class DenIntegrator: virtual public SavableState {
  protected:
    Ref<GaussianBasisSet> basis_;
    Ref<Integral> integral_;
 
//clj    Ref<ShellExtent> extent_;
    Ref<BatchElectronDensity> den_;
 
    Ref<ThreadGrp> threadgrp_;
    Ref<MessageGrp> messagegrp_;
 
    double value_;
    double accuracy_;
 
    double *alpha_vmat_;
    double *beta_vmat_;
 
//clj    double *alpha_dmat_;
//clj    double *beta_dmat_;
//clj    double *dmat_bound_;
 
    int spin_polarized_;
 
    int need_density_; // specialization must set to 1 if it needs density_
    double density_;
    int nbasis_;
    int nshell_;
    int n_integration_center_;
    int natom_;
    int compute_potential_integrals_; // 1 if potential integrals are needed
 
    int linear_scaling_;
    int use_dmat_bound_;
 
    void init_integration(const Ref<DenFunctional> &func,
                          const RefSymmSCMatrix& densa,
                          const RefSymmSCMatrix& densb,
                          double *nuclear_gradient);
    void done_integration();
 
    void init_object();
  public:
    DenIntegrator();
    DenIntegrator(const Ref<KeyVal> &);
    DenIntegrator(StateIn &);
    ~DenIntegrator();
    void save_data_state(StateOut &);
 
    const Ref<GaussianBasisSet> &basis() const { return basis_; }
    const Ref<Integral> &integral() const { return integral_; }
    double value() const { return value_; }
 
    void set_accuracy(double a);
    double get_accuracy(void) {return accuracy_; }
    void set_compute_potential_integrals(int);
    const double *alpha_vmat() const { return alpha_vmat_; }
    const double *beta_vmat() const { return beta_vmat_; }
 
    virtual void init(const Ref<Wavefunction> &);
    virtual void init(const Ref<GaussianBasisSet> &,
                      const Ref<Integral> &);
    virtual void done();
    virtual void integrate(const Ref<DenFunctional> &,
                           const RefSymmSCMatrix& densa =0,
                           const RefSymmSCMatrix& densb =0,
                           double *nuclear_grad = 0) = 0;
 
    bool spin_polarized() const { return spin_polarized_; }
};
 
 
class IntegrationWeight: virtual public SavableState {
 
  protected:
 
    Ref<Molecule> mol_;
    double tol_;
 
    void fd_w(int icenter, SCVector3 &point, double *fd_grad_w);
 
  public:
    IntegrationWeight();
    IntegrationWeight(const Ref<KeyVal> &);
    IntegrationWeight(StateIn &);
    ~IntegrationWeight();
    void save_data_state(StateOut &);
 
    void test(int icenter, SCVector3 &point);
    void test();
 
    virtual void init(const Ref<Molecule> &, double tolerance);
    virtual void done();
    virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0;
};
 
 
class BeckeIntegrationWeight: public IntegrationWeight {
 
    int n_integration_centers;
    SCVector3 *centers;
    double *atomic_radius;
 
    double **a_mat;
    double **oorab;
 
    void compute_grad_p(int gc, int ic, int wc, SCVector3&r, double p,
                           SCVector3&g);
    void compute_grad_nu(int gc, int bc, SCVector3 &point, SCVector3 &grad);
 
    double compute_t(int ic, int jc, SCVector3 &point);
    double compute_p(int icenter, SCVector3&point);
 
  public:
    BeckeIntegrationWeight();
    BeckeIntegrationWeight(const Ref<KeyVal> &);
    BeckeIntegrationWeight(StateIn &);
    ~BeckeIntegrationWeight();
    void save_data_state(StateOut &);
 
    void init(const Ref<Molecule> &, double tolerance);
    void done();
    double w(int center, SCVector3 &point, double *grad_w = 0);
};
 
class RadialIntegrator: virtual public SavableState {
  protected:
    int nr_;
  public:
    RadialIntegrator();
    RadialIntegrator(const Ref<KeyVal> &);
    RadialIntegrator(StateIn &);
    ~RadialIntegrator();
    void save_data_state(StateOut &);
 
    virtual int nr() const = 0;
    virtual double radial_value(int ir, int nr, double radii,
                                double &multiplier) = 0;
};
 
 
class AngularIntegrator: virtual public SavableState{
  protected:
  public:
    AngularIntegrator();
    AngularIntegrator(const Ref<KeyVal> &);
    AngularIntegrator(StateIn &);
    ~AngularIntegrator();
    void save_data_state(StateOut &);
 
    virtual int nw(void) const = 0;
    virtual int num_angular_points(double r_value, int ir) = 0;
    virtual double angular_point_cartesian(int iangular, double r,
        SCVector3 &integration_point) const = 0;
};
 
 
class EulerMaclaurinRadialIntegrator: public RadialIntegrator {
  public:
    EulerMaclaurinRadialIntegrator();
    EulerMaclaurinRadialIntegrator(int i);
    EulerMaclaurinRadialIntegrator(const Ref<KeyVal> &);
    EulerMaclaurinRadialIntegrator(StateIn &);
    ~EulerMaclaurinRadialIntegrator();
    void save_data_state(StateOut &);
 
    int nr() const;
    double radial_value(int ir, int nr, double radii, double &multiplier);
 
    void print(std::ostream & =ExEnv::out0()) const;
};
 
class LebedevLaikovIntegrator: public AngularIntegrator {
  protected:
    int npoint_;
    double *x_, *y_, *z_, *w_;
    
    void init(int n);
  public:
    LebedevLaikovIntegrator();
    LebedevLaikovIntegrator(const Ref<KeyVal> &);
    LebedevLaikovIntegrator(StateIn &);
    LebedevLaikovIntegrator(int);
    ~LebedevLaikovIntegrator();
    void save_data_state(StateOut &);
 
    int nw(void) const;
    int num_angular_points(double r_value, int ir);
    double angular_point_cartesian(int iangular, double r,
                                   SCVector3 &integration_point) const;
    void print(std::ostream & =ExEnv::out0()) const;
};
 
class GaussLegendreAngularIntegrator: public AngularIntegrator {
  protected:
    int ntheta_;
    int nphi_;
    int Ktheta_;
    int ntheta_r_;
    int nphi_r_;
    int Ktheta_r_;
    double *theta_quad_weights_;
    double *theta_quad_points_;
 
    int get_ntheta(void) const;
    void set_ntheta(int i);
    int get_nphi(void) const;
    void set_nphi(int i);
    int get_Ktheta(void) const;
    void set_Ktheta(int i);
    int get_ntheta_r(void) const;
    void set_ntheta_r(int i);
    int get_nphi_r(void) const;
    void set_nphi_r(int i);
    int get_Ktheta_r(void) const;
    void set_Ktheta_r(int i);
    int nw(void) const;
    double sin_theta(SCVector3 &point) const;
    void gauleg(double x1, double x2, int n);    
  public:
    GaussLegendreAngularIntegrator();
    GaussLegendreAngularIntegrator(const Ref<KeyVal> &);
    GaussLegendreAngularIntegrator(StateIn &);
    ~GaussLegendreAngularIntegrator();
    void save_data_state(StateOut &);
    
    int num_angular_points(double r_value, int ir);
    double angular_point_cartesian(int iangular, double r,
        SCVector3 &integration_point) const;
    void print(std::ostream & =ExEnv::out0()) const;
};
 
class RadialAngularIntegrator: public DenIntegrator {
  private:
    int prune_grid_;
    double **Alpha_coeffs_;
    int gridtype_;
    int **nr_points_, *xcoarse_l_;
    int npruned_partitions_;
    double *grid_accuracy_;
    int dynamic_grids_;
    int natomic_rows_;
    int max_gridtype_;
  protected:
    Ref<IntegrationWeight> weight_;
    Ref<RadialIntegrator> radial_user_;
    Ref<AngularIntegrator> angular_user_;
    Ref<AngularIntegrator> ***angular_grid_;
    Ref<RadialIntegrator> **radial_grid_;
  public:
    RadialAngularIntegrator();
    RadialAngularIntegrator(const Ref<KeyVal> &);
    RadialAngularIntegrator(StateIn &);
    ~RadialAngularIntegrator();
    void save_data_state(StateOut &);
 
    void integrate(const Ref<DenFunctional> &,
                   const RefSymmSCMatrix& densa =0,
                   const RefSymmSCMatrix& densb =0,
                   double *nuclear_gradient = 0);
    void print(std::ostream & =ExEnv::out0()) const;
    AngularIntegrator *get_angular_grid(double radius, double atomic_radius,
                                        int charge, int deriv_order);
    RadialIntegrator *get_radial_grid(int charge, int deriv_order);
    void init_default_grids(void);
    int angular_grid_offset(int i);
    void set_grids(void);
    int get_atomic_row(int i);
    void init_parameters(void);
    void init_parameters(const Ref<KeyVal>& keyval);
    void init_pruning_coefficients(const Ref<KeyVal>& keyval);
    void init_pruning_coefficients(void);
    void init_alpha_coefficients(void);
    int select_dynamic_grid(void);
    Ref<IntegrationWeight> weight() { return weight_; }
};
 
}
    
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: