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
 
#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>
#include <chemistry/qc/wfn/orbitalspace.h>
#ifdef MPQC_NEW_FEATURES
#include <util/misc/xml.h>
#endif // MPQC_NEW_FEATURES
 
 
namespace sc {
 
 
class Wavefunction: public MolecularEnergy {
 
    RefSCDimension aodim_;
    RefSCDimension sodim_;
    Ref<SCMatrixKit> basiskit_;
 
    ResultRefSymmSCMatrix overlap_;
    ResultRefSymmSCMatrix hcore_;
 
  protected:
    ResultRefSCMatrix natural_orbitals_;
    ResultRefDiagSCMatrix natural_density_;
    Ref<GaussianBasisSet> gbs_;
 
  private:
    double * bs_values;
    double * bsg_values;
 
    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_;
 
    mutable double magnetic_moment_; 
 
    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);
 
    bool nonzero_efield_supported() const;
 
  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&);
#ifdef MPQC_NEW_FEATURES
    virtual boost::property_tree::ptree& write_xml(boost::property_tree::ptree& parent, const XMLWriter& writer);
#endif
 
    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);
    // same as above, except computes all orbital AND orbs must be transposed (i.e. mo rows, ao columns)
    void orbitals(const SCVector3& r, const RefSCMatrix& orbs, RefSCVector& values);
    // same as above except it evaluates on a set of values, and is static
    static void orbitals(const Ref<OrbitalSpace>& orbs,
                         const std::vector<SCVector3>& r,
                         std::vector<double>& values);
 
    double orbital_density(const SCVector3& r,
                           int iorb,
                           const RefSCMatrix& orbs,
                           double* orbval = 0);
 
    double total_charge() const;
    virtual int nelectron() = 0;
    virtual double magnetic_moment() const;
 
    virtual RefSymmSCMatrix density() = 0;
    virtual RefSymmSCMatrix ao_density();
    virtual RefSCMatrix natural_orbitals();
    virtual RefDiagSCMatrix natural_density();
 
    int spin_polarized() { return magnetic_moment() != 0.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;
    virtual 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;
 
 
    void writeorbitals();
};
 
// end of addtogroup ChemistryElectronicStructure
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: