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

#ifdef __GNUC__
#pragma interface
#endif

#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:
    Ref<Wavefunction> wfn_;

    Ref<GaussianBasisSet> basis_;

    // 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_;
    int 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(bool initialize_density_matrices);
    // 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 *restrictxx dmat,
                              double *restrictxx rho,
                              double *restrictxx grad,
                              double *restrictxx 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<Wavefunction>&, 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(bool initialize_density_matrices = true);

    virtual void set_densities(const RefSymmSCMatrix &aden,
                               const RefSymmSCMatrix &bden);

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

    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_;
    char *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:

---