src/lib/chemistry/qc/wfn/wfn.h

//
// wfn.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_wfn_h
#define _chemistry_qc_wfn_wfn_h

#ifdef __GNUC__
#pragma interface
#endif

#include <iostream>

#include <util/misc/compute.h>
#include <math/scmat/matrix.h>
#include <math/scmat/vector3.h>
#include <chemistry/molecule/energy.h>
#include <chemistry/qc/basis/basis.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/orthog.h>

namespace sc {

class Wavefunction: public MolecularEnergy {

    RefSCDimension aodim_;
    RefSCDimension sodim_;
    Ref<SCMatrixKit> basiskit_;
    
    ResultRefSymmSCMatrix overlap_;
    ResultRefSymmSCMatrix hcore_;
    ResultRefSCMatrix natural_orbitals_;
    ResultRefDiagSCMatrix natural_density_;

    double * bs_values;
    double * bsg_values;

    Ref<GaussianBasisSet> gbs_;
    Ref<Integral> integral_;

    Ref<GaussianBasisSet> atom_basis_;
    double * atom_basis_coef_;

    double lindep_tol_;
    OverlapOrthog::OrthogMethod orthog_method_;
    Ref<OverlapOrthog> orthog_;

    int print_nao_;
    int print_npa_;

    int dk_; // 0, 1, 2, or 3
    Ref<GaussianBasisSet> momentum_basis_;

    void init_orthog();

    void set_up_charge_types(std::vector<int> &q_pc,
                             std::vector<int> &q_cd,
                             std::vector<int> &n_pc,
                             std::vector<int> &n_cd);

    double nuc_rep_pc_pc(const std::vector<int>&,const std::vector<int>&,bool);
    double nuc_rep_pc_cd(const std::vector<int>&,const std::vector<int>&);
    double nuc_rep_cd_cd(const std::vector<int>&,const std::vector<int>&,bool);
    void scale_atom_basis_coef();

    void nuc_rep_grad_pc_pc(double **grad,
                            const std::vector<int>&c1,
                            const std::vector<int>&c2,
                            bool uniq);
    void nuc_rep_grad_pc_cd(double **grad,
                            const std::vector<int>&c1,
                            const std::vector<int>&c2);
    void nuc_rep_grad_cd_cd(double **grad,
                            const std::vector<int>&c1,
                            const std::vector<int>&c2,
                            bool uniq);

    // Compute the dk relativistic core hamiltonian.
    RefSymmSCMatrix core_hamiltonian_dk(
      int dk,
      const Ref<GaussianBasisSet> &bas,
      const Ref<GaussianBasisSet> &pbas = 0);
    void core_hamiltonian_dk2_contrib(const RefSymmSCMatrix &h_pbas,
                                      const RefDiagSCMatrix &E,
                                      const RefDiagSCMatrix &K,
                                      const RefDiagSCMatrix &p2,
                                      const RefDiagSCMatrix &p2K2,
                                      const RefDiagSCMatrix &p2K2_inv,
                                      const RefSymmSCMatrix &AVA_pbas,
                                      const RefSymmSCMatrix &BpVpB_pbas);
    void core_hamiltonian_dk3_contrib(const RefSymmSCMatrix &h_pbas,
                                      const RefDiagSCMatrix &E,
                                      const RefDiagSCMatrix &B,
                                      const RefDiagSCMatrix &p2K2_inv,
                                      const RefSCMatrix &so_to_p,
                                      const RefSymmSCMatrix &pxVp);


  protected:

    int debug_;

    double min_orthog_res();
    double max_orthog_res();

    void copy_orthog_info(const Ref<Wavefunction> &);
    
  public:
    Wavefunction(StateIn&);
    Wavefunction(const Ref<KeyVal>&);
    virtual ~Wavefunction();

    void save_data_state(StateOut&);

    double density(const SCVector3&);
    double density_gradient(const SCVector3&,double*);
    double natural_orbital(const SCVector3& r, int iorb);
    double natural_orbital_density(const SCVector3& r,
                                   int orb, double* orbval = 0);
    double orbital(const SCVector3& r, int iorb, const RefSCMatrix& orbs);

    double orbital_density(const SCVector3& r,
                           int iorb,
                           const RefSCMatrix& orbs,
                           double* orbval = 0);

    double charge();
    virtual int nelectron() = 0;

    virtual RefSymmSCMatrix density() = 0;
    virtual RefSymmSCMatrix ao_density();
    virtual RefSCMatrix natural_orbitals();
    virtual RefDiagSCMatrix natural_density();

    virtual int spin_polarized() = 0;

    int dk() const { return dk_; }

    virtual RefSymmSCMatrix alpha_density();
    virtual RefSymmSCMatrix beta_density();
    virtual RefSymmSCMatrix alpha_ao_density();
    virtual RefSymmSCMatrix beta_ao_density();

    virtual RefSCMatrix nao(double *atom_charges=0);

    virtual RefSymmSCMatrix overlap();
    virtual RefSymmSCMatrix core_hamiltonian_for_basis(
      const Ref<GaussianBasisSet> &bas,
      const Ref<GaussianBasisSet> &pbas = 0);
    virtual RefSymmSCMatrix core_hamiltonian();
    // Compute the non-relativistic core hamiltonian.
    RefSymmSCMatrix core_hamiltonian_nr(
      const Ref<GaussianBasisSet> &bas);


    virtual double nuclear_repulsion_energy();
    void nuclear_repulsion_energy_gradient(double *g);
    virtual void nuclear_repulsion_energy_gradient(double **g);

    RefSCDimension ao_dimension();
    RefSCDimension so_dimension();
    RefSCDimension oso_dimension();
    Ref<SCMatrixKit> basis_matrixkit();
    Ref<Molecule> molecule() const;
    Ref<GaussianBasisSet> basis() const;
    Ref<GaussianBasisSet> momentum_basis() const;
    Ref<GaussianBasisSet> atom_basis() const;
    const double *atom_basis_coef() const;
    Ref<Integral> integral();

    // override symmetry_changed from MolecularEnergy
    void symmetry_changed();

    RefSCMatrix so_to_orthog_so();

    RefSCMatrix so_to_orthog_so_inverse();

    OverlapOrthog::OrthogMethod orthog_method() const;
    void set_orthog_method(const OverlapOrthog::OrthogMethod&);

    double lindep_tol() const;
    void set_lindep_tol(double);

    void obsolete();

    void print(std::ostream& = ExEnv::out0()) const;
};

}

#endif

// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End:

---