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
#ifdef __GNUC__
#pragma interface
#endif
#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<Wavefunction> wfn_;
//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 &);
Ref<Wavefunction> wavefunction() const { return wfn_; }
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 done();
virtual void integrate(const Ref<DenFunctional> &,
const RefSymmSCMatrix& densa =0,
const RefSymmSCMatrix& densb =0,
double *nuclear_grad = 0) = 0;
};
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: