density.h

//
// density.h
//
// Copyright (C) 1996 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_wfn_density_h
#define _chemistry_qc_wfn_density_h
 
#include <math/isosurf/volume.h>
#include <chemistry/qc/wfn/wfn.h>
#include <chemistry/qc/basis/extent.h>
#include <chemistry/molecule/molrender.h>
#include <chemistry/molecule/molecule.h>
#include <math/mmisc/grid.h>
 
namespace sc {
 
class ElectronDensity: public Volume {
  protected:
    Ref<Wavefunction> wfn_;
    virtual void compute();
  public:
    ElectronDensity(const Ref<KeyVal>&);
    ElectronDensity(const Ref<Wavefunction>&);
    ~ElectronDensity();
    virtual void boundingbox(double valuemin,
                             double valuemax,
                             SCVector3& p1, SCVector3& p2);
};
 
class BatchElectronDensity: public Volume {
    void zero_pointers();
  protected:
    // wfn_ might be null in which case basis_ and integral_
    // must be initialized and set_densities must be used.
    Ref<Wavefunction> wfn_;
 
    Ref<GaussianBasisSet> basis_;
    Ref<Integral> integral_;
 
    bool initialized_;
 
    // shared between threads
    double *alpha_dmat_;
    double *beta_dmat_;
    double *dmat_bound_;
    ShellExtent *extent_;
 
    // private data
    GaussianBasisSet::ValueData *valdat_;
    int ncontrib_;
    int *contrib_;
    int ncontrib_bf_;
    int *contrib_bf_;
    double *bs_values_;
    double *bsg_values_;
    double *bsh_values_;
 
    int nshell_;
    int nbasis_;
    bool spin_polarized_;
    int linear_scaling_;
    int use_dmat_bound_;
 
    bool need_hessian_, need_gradient_;
    bool need_basis_hessian_, need_basis_gradient_;
 
    bool using_shared_data_;
 
    double accuracy_;
    virtual void init_common_data();
    // this must be called after common data is initialized,
    // either with init_common_data or by copying
    virtual void init_scratch_data();
    void compute_basis_values(const SCVector3&r);
    void compute_spin_density(const double *RESTRICT dmat,
                              double *RESTRICT rho,
                              double *RESTRICT grad,
                              double *RESTRICT hess);
 
    virtual void compute();
  public:
 
    enum {X=0, Y=1, Z=2};
    enum {XX=0, YX=1, YY=2, ZX=3, ZY=4, ZZ=5};
 
    BatchElectronDensity(const Ref<KeyVal>&);
    BatchElectronDensity(const Ref<GaussianBasisSet> &basis,
                         const Ref<Integral> &integral,
                         double accuracy=DBL_EPSILON);
    BatchElectronDensity(const Ref<Wavefunction> &wfn,
                         double accuracy=DBL_EPSILON);
    BatchElectronDensity(const Ref<BatchElectronDensity>& d,
                         bool reference_parent_data=false);
    ~BatchElectronDensity();
    virtual void boundingbox(double valuemin,
                             double valuemax,
                             SCVector3& p1, SCVector3& p2);
 
    void clear();
 
    void compute_density(const SCVector3 &r,
                         double *alpha_density,
                         double *alpha_density_grad,
                         double *alpha_density_hessian,
                         double *beta_density,
                         double *beta_density_grad,
                         double *beta_density_hessian);
 
    virtual void init();
 
    virtual void set_densities(const RefSymmSCMatrix &aden,
                               const RefSymmSCMatrix &bden = 0);
 
    virtual void set_densities(const Ref<Wavefunction> &wfn);
 
    void set_linear_scaling(bool b) { linear_scaling_ = b; }
 
    void set_accuracy(double a) { accuracy_ = a; }
 
    void set_use_dmat_bound(bool b) { use_dmat_bound_ = b; }
 
    bool spin_polarized() const { return spin_polarized_; }
 
    double *alpha_density_matrix() { return alpha_dmat_; }
    double *beta_density_matrix()
        { return (spin_polarized_?beta_dmat_:alpha_dmat_); }
    int ncontrib() { return ncontrib_; }
    int *contrib() { return contrib_; }
    int ncontrib_bf() { return ncontrib_bf_; }
    int *contrib_bf() { return contrib_bf_; }
    double *bs_values() { return bs_values_; }
    double *bsg_values() { return bsg_values_; }
    double *bsh_values() { return bsh_values_; }
    void set_need_basis_gradient(bool b) { need_basis_gradient_ = b; }
    void set_need_basis_hessian(bool b) { need_basis_hessian_ = b; }
};
 
class WriteElectronDensity: public WriteGrid {
  private:
    double df_alpha(double alpha, double beta);
    double df_beta(double alpha, double beta);
    double df_sum(double alpha, double beta);
    double df_spin(double alpha, double beta);
  protected:
    Ref<Wavefunction> wfn_;
    Ref<BatchElectronDensity> bed_;
    double accuracy_;
    std::string type_;
    double (WriteElectronDensity::*density_function_)(double, double);
    
    void initialize();
    void label(char* buffer);
    Ref<Molecule> get_molecule();
    double calculate_value(SCVector3 point);
  public:
    WriteElectronDensity(const Ref<KeyVal> &);
};
 
class DensityColorizer: public MoleculeColorizer {
  protected:
    Ref<Wavefunction> wfn_;
    double scale_;
    double reference_;
    int have_scale_;
    int have_reference_;
  public:
    DensityColorizer(const Ref<KeyVal>&);
    ~DensityColorizer();
 
    void colorize(const Ref<RenderedPolygons> &);
};
 
class GradDensityColorizer: public MoleculeColorizer {
  protected:
    Ref<Wavefunction> wfn_;
    double scale_;
    double reference_;
    int have_scale_;
    int have_reference_;
  public:
    GradDensityColorizer(const Ref<KeyVal>&);
    ~GradDensityColorizer();
 
    void colorize(const Ref<RenderedPolygons> &);
};
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: