integrator.h

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// integrator.h --- definition of the electron density integrator
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// 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>

namespace sc {

/** An abstract base class for integrating the electron density. */
class DenIntegrator: virtual public SavableState {
  protected:
    Ref<Wavefunction> wfn_;
    Ref<ShellExtent> extent_;

    Ref<ThreadGrp> threadgrp_;
    Ref<MessageGrp> messagegrp_;

    double value_;
    double accuracy_;

    double *alpha_vmat_;
    double *beta_vmat_;

    double *alpha_dmat_;
    double *beta_dmat_;
    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 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:
    /// Construct a new DenIntegrator.
    DenIntegrator();
    /// Construct a new DenIntegrator given the KeyVal input.
    DenIntegrator(const Ref<KeyVal> &);
    /// Construct a new DenIntegrator given the StateIn data.
    DenIntegrator(StateIn &);
    ~DenIntegrator();
    void save_data_state(StateOut &);

    /// Returns the wavefunction used for the integration.
    Ref<Wavefunction> wavefunction() const { return wfn_; }
    /// Returns the result of the integration.
    double value() const { return value_; }

    /// Sets the accuracy to use in the integration.
    void set_accuracy(double a) { accuracy_ = a; }
    double get_accuracy(void) {return accuracy_; }
    /** Call with non zero if the potential integrals are to be computed.
        They can be returned with the vmat() member. */
    void set_compute_potential_integrals(int);
    /** Returns the alpha potential integrals. Stored as
        the lower triangular, row-major format. */
    const double *alpha_vmat() const { return alpha_vmat_; }
    /** Returns the beta potential integrals. Stored as
        the lower triangular, row-major format. */
    const double *beta_vmat() const { return beta_vmat_; }

    /** Called before integrate.  Does not need to be called again
        unless the geometry changes or done is called. */
    virtual void init(const Ref<Wavefunction> &);
    /// Must be called between calls to init.
    virtual void done();
    /** Performs the integration of the given functional using the given
        alpha and beta density matrices.  The nuclear derivative
        contribution is placed in nuclear_grad, if it is non-null. */
    virtual void integrate(const Ref<DenFunctional> &,
                           const RefSymmSCMatrix& densa =0,
                           const RefSymmSCMatrix& densb =0,
                           double *nuclear_grad = 0) = 0;
};


/** An abstract base class for computing grid weights. */
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();

    /// Initialize the integration weight object.
    virtual void init(const Ref<Molecule> &, double tolerance);
    /// Called when finished with the integration weight object.
    virtual void done();
    /** Computes the weight for a given center at a given point in space.
        Derivatives of the weigth with respect to nuclear coordinates are
        optionally returned in grad_w.  This must be called after init, but
        before done. It must also be thread-safe. */
    virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0;
};


/** Implements Becke's integration weight scheme. */
class BeckeIntegrationWeight: public IntegrationWeight {

    int ncenters;
    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);
};

/** An abstract base class for radial integrators. */
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;
};


/** An abstract base class for angular integrators. */
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;
};


/** An implementation of a radial integrator using the Euler-Maclaurin
    weights and grid points. */
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;
};

/** An implementation of a Lebedev angular integrator.  It uses code
    written by Dr. Dmitri N. Laikov.

    This can generate grids with the following numbers of points: 6, 14,
    26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 386, 434,
    482, 530, 590, 650, 698, 770, 830, 890, 974, 1046, 1118, 1202, 1274,
    1358, 1454, 1538, 1622, 1730, 1814, 1910, 2030, 2126, 2222, 2354, 2450,
    2558, 2702, 2810, 2930, 3074, 3182, 3314, 3470, 3590, 3722, 3890, 4010,
    4154, 4334, 4466, 4610, 4802, 4934, 5090, 5294, 5438, 5606, and 5810.

    V.I. Lebedev, and D.N. Laikov
    "A quadrature formula for the sphere of the 131st
     algebraic order of accuracy"
    Doklady Mathematics, Vol. 59, No. 3, 1999, pp. 477-481.
   
    V.I. Lebedev
    "A quadrature formula for the sphere of 59th algebraic
     order of accuracy"
    Russian Acad. Sci. Dokl. Math., Vol. 50, 1995, pp. 283-286.
   
    V.I. Lebedev, and A.L. Skorokhodov
    "Quadrature formulas of orders 41, 47, and 53 for the sphere"
    Russian Acad. Sci. Dokl. Math., Vol. 45, 1992, pp. 587-592.
   
    V.I. Lebedev
    "Spherical quadrature formulas exact to orders 25-29"
    Siberian Mathematical Journal, Vol. 18, 1977, pp. 99-107.
   
    V.I. Lebedev
    "Quadratures on a sphere"
    Computational Mathematics and Mathematical Physics, Vol. 16,
    1976, pp. 10-24.
   
    V.I. Lebedev
    "Values of the nodes and weights of ninth to seventeenth
     order Gauss-Markov quadrature formulae invariant under the
     octahedron group with inversion"
    Computational Mathematics and Mathematical Physics, Vol. 15,
       1975, pp. 44-51.

 */
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;
};

/** An implementation of an angular integrator using the Gauss-Legendre
    weights and grid points. */
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;
};

/** An implementation of an integrator using any combination of
    a RadianIntegrator and an AngularIntegrator. */
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: