1 // 2 // wfn.h 3 // 4 // Copyright (C) 1996 Limit Point Systems, Inc. 5 // 6 // Author: Curtis Janssen <cljanss@limitpt.com> 7 // Maintainer: LPS 8 // 9 // This file is part of the SC Toolkit. 10 // 11 // The SC Toolkit is free software; you can redistribute it and/or modify 12 // it under the terms of the GNU Library General Public License as published by 13 // the Free Software Foundation; either version 2, or (at your option) 14 // any later version. 15 // 16 // The SC Toolkit is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU Library General Public License for more details. 20 // 21 // You should have received a copy of the GNU Library General Public License 22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to 23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 24 // 25 // The U.S. Government is granted a limited license as per AL 91-7. 26 // 27 28 #ifndef _chemistry_qc_wfn_wfn_h 29 #define _chemistry_qc_wfn_wfn_h 30 31 #ifdef __GNUC__ 32 #pragma interface 33 #endif 34 35 #include <iostream> 36 37 #include <util/misc/compute.h> 38 #include <math/scmat/matrix.h> 39 #include <math/scmat/vector3.h> 40 #include <chemistry/molecule/energy.h> 41 #include <chemistry/qc/basis/basis.h> 42 #include <chemistry/qc/basis/integral.h> 43 #include <chemistry/qc/basis/orthog.h> 44 45 namespace sc { 46 47 /** A Wavefunction is a MolecularEnergy that utilizies a GaussianBasisSet. */ 48 class Wavefunction: public MolecularEnergy { 49 50 RefSCDimension aodim_; 51 RefSCDimension sodim_; 52 Ref<SCMatrixKit> basiskit_; 53 54 ResultRefSymmSCMatrix overlap_; 55 ResultRefSymmSCMatrix hcore_; 56 ResultRefSCMatrix natural_orbitals_; 57 ResultRefDiagSCMatrix natural_density_; 58 59 double * bs_values; 60 double * bsg_values; 61 62 Ref<GaussianBasisSet> gbs_; 63 Ref<Integral> integral_; 64 65 Ref<GaussianBasisSet> atom_basis_; 66 double * atom_basis_coef_; 67 68 double lindep_tol_; 69 OverlapOrthog::OrthogMethod orthog_method_; 70 Ref<OverlapOrthog> orthog_; 71 72 int print_nao_; 73 int print_npa_; 74 75 void init_orthog(); 76 77 void set_up_charge_types(std::vector<int> &q_pc, 78 std::vector<int> &q_cd, 79 std::vector<int> &n_pc, 80 std::vector<int> &n_cd); 81 82 double nuc_rep_pc_pc(const std::vector<int>&,const std::vector<int>&,bool); 83 double nuc_rep_pc_cd(const std::vector<int>&,const std::vector<int>&); 84 double nuc_rep_cd_cd(const std::vector<int>&,const std::vector<int>&,bool); 85 void scale_atom_basis_coef(); 86 87 void nuc_rep_grad_pc_pc(double **grad, 88 const std::vector<int>&c1, 89 const std::vector<int>&c2, 90 bool uniq); 91 void nuc_rep_grad_pc_cd(double **grad, 92 const std::vector<int>&c1, 93 const std::vector<int>&c2); 94 void nuc_rep_grad_cd_cd(double **grad, 95 const std::vector<int>&c1, 96 const std::vector<int>&c2, 97 bool uniq); 98 99 protected: 100 101 int debug_; 102 103 double min_orthog_res(); 104 double max_orthog_res(); 105 106 void copy_orthog_info(const Ref<Wavefunction> &); 107 108 public: 109 Wavefunction(StateIn&); 110 /** The KeyVal constructor. 111 112 <dl> 113 114 <dt><tt>basis</tt><dd> Specifies a GaussianBasisSet object. There 115 is no default. 116 117 <dt><tt>integral</tt><dd> Specifies an Integral object that 118 computes the two electron integrals. The default is a IntegralV3 119 object. 120 121 <dt><tt>orthog_method</tt><dd> This is a string that specifies the 122 orthogonalization method to be used. It can be one one canonical, 123 gramschmidt, or symmetric. The default is symmetric. 124 125 <dt><tt>lindep_tol</tt><dd> The tolerance used to detect linearly 126 dependent basis functions. The precise meaning depends on the 127 orthogonalization method. The default value is 1e-8. 128 129 <dt><tt>print_nao</tt><dd> This specifies a boolean value. If true 130 the natural atomic orbitals will be printed. Not all wavefunction 131 will be able to do this. The default is false. 132 133 <dt><tt>print_npa</tt><dd> This specifies a boolean value. If true 134 the natural population analysis will be printed. Not all 135 wavefunction will be able to do this. The default is true if 136 print_nao is true, otherwise it is false. 137 138 <dt><tt>debug</tt><dd> This integer can be used to produce output 139 for debugging. The default is 0. 140 141 </dl> */ 142 Wavefunction(const Ref<KeyVal>&); 143 virtual ~Wavefunction(); 144 145 void save_data_state(StateOut&); 146 147 double density(const SCVector3&); 148 double density_gradient(const SCVector3&,double*); 149 double natural_orbital(const SCVector3& r, int iorb); 150 double natural_orbital_density(const SCVector3& r, 151 int orb, double* orbval = 0); 152 double orbital(const SCVector3& r, int iorb, const RefSCMatrix& orbs); 153 154 double orbital_density(const SCVector3& r, 155 int iorb, 156 const RefSCMatrix& orbs, 157 double* orbval = 0); 158 159 /// Returns the charge. 160 double charge(); 161 /// Returns the number of electrons. 162 virtual int nelectron() = 0; 163 164 /// Returns the SO density. 165 virtual RefSymmSCMatrix density() = 0; 166 /// Returns the AO density. 167 virtual RefSymmSCMatrix ao_density(); 168 /// Returns the natural orbitals. 169 virtual RefSCMatrix natural_orbitals(); 170 /// Returns the natural density (a diagonal matrix). 171 virtual RefDiagSCMatrix natural_density(); 172 173 /// Return 1 if the alpha density is not equal to the beta density. 174 virtual int spin_polarized() = 0; 175 176 /// Return alpha electron densities in the SO basis. 177 virtual RefSymmSCMatrix alpha_density(); 178 /// Return beta electron densities in the SO basis. 179 virtual RefSymmSCMatrix beta_density(); 180 /// Return alpha electron densities in the AO basis. 181 virtual RefSymmSCMatrix alpha_ao_density(); 182 /// Return beta electron densities in the AO basis. 183 virtual RefSymmSCMatrix beta_ao_density(); 184 185 /// returns the ao to nao transformation matrix 186 virtual RefSCMatrix nao(double *atom_charges=0); 187 188 /// Returns the SO overlap matrix. 189 virtual RefSymmSCMatrix overlap(); 190 /// Returns the SO core Hamiltonian. 191 virtual RefSymmSCMatrix core_hamiltonian(); 192 193 /** Returns the nuclear repulsion energy. This must be used instead of 194 Molecule::nuclear_repulsion_energy() since there may be diffuse 195 atomic charges. */ 196 virtual double nuclear_repulsion_energy(); 197 /** Computes the nuclear repulsion gradient. This must be used instead 198 of Molecule::nuclear_repulsion_1der() since there may be diffuse 199 atomic charges. The gradient, g, is zeroed and set to x_0, y_0, 200 z_0, x_1, ... . */ 201 void nuclear_repulsion_energy_gradient(double *g); 202 /** Computes the nuclear repulsion gradient. This must be used instead 203 of Molecule::nuclear_repulsion_1der() since there may be diffuse 204 atomic charges. The gradient, g, is first zeroed. Its dimensions 205 are g[natom][3]. */ 206 virtual void nuclear_repulsion_energy_gradient(double **g); 207 208 /// Atomic orbital dimension. 209 RefSCDimension ao_dimension(); 210 /// Symmetry adapted orbital dimension. 211 RefSCDimension so_dimension(); 212 /// Orthogonalized symmetry adapted orbital dimension. 213 RefSCDimension oso_dimension(); 214 /// Matrix kit for AO, SO, orthogonalized SO, and MO dimensioned matrices. 215 Ref<SCMatrixKit> basis_matrixkit(); 216 /// Returns the Molecule. 217 Ref<Molecule> molecule() const; 218 /// Returns the basis set. 219 Ref<GaussianBasisSet> basis() const; 220 /// Returns the basis set describing the nuclear charge distributions 221 Ref<GaussianBasisSet> atom_basis() const; 222 /** Returns the coefficients of the nuclear charge distribution basis 223 * functions. */ 224 const double *atom_basis_coef() const; 225 /// Returns the integral evaluator. 226 Ref<Integral> integral(); 227 228 // override symmetry_changed from MolecularEnergy 229 void symmetry_changed(); 230 231 /** Returns a matrix which does the default transform from SO's to 232 orthogonal SO's. This could be either the symmetric or canonical 233 orthogonalization matrix. The row dimension is SO and the column 234 dimension is ortho SO. An operator \f$O\f$ in the ortho SO basis is 235 given by \f$X O X^T\f$ where \f$X\f$ is the return value of this 236 function. */ 237 RefSCMatrix so_to_orthog_so(); 238 239 /** Returns the inverse of the transformation returned by so_to_orthog_so. 240 */ 241 RefSCMatrix so_to_orthog_so_inverse(); 242 243 /// Returns the orthogonalization method 244 OverlapOrthog::OrthogMethod orthog_method() const; 245 /// (Re)Sets the orthogonalization method and makes this obsolete 246 void set_orthog_method(const OverlapOrthog::OrthogMethod&); 247 248 /// Returns the tolerance for linear dependencies. 249 double lindep_tol() const; 250 /// Re(Sets) the tolerance for linear dependencies. 251 void set_lindep_tol(double); 252 253 void obsolete(); 254 255 void print(std::ostream& = ExEnv::out0()) const; 256 }; 257 258 } 259 260 #endif 261 262 // Local Variables: 263 // mode: c++ 264 // c-file-style: "ETS" 265 // End: 266