MPQC 2.3.1
|
00001 // 00002 // integrator.h --- definition of the electron density integrator 00003 // 00004 // Copyright (C) 1997 Limit Point Systems, Inc. 00005 // 00006 // Author: Curtis Janssen <cljanss@limitpt.com> 00007 // Maintainer: LPS 00008 // 00009 // This file is part of the SC Toolkit. 00010 // 00011 // The SC Toolkit is free software; you can redistribute it and/or modify 00012 // it under the terms of the GNU Library General Public License as published by 00013 // the Free Software Foundation; either version 2, or (at your option) 00014 // any later version. 00015 // 00016 // The SC Toolkit is distributed in the hope that it will be useful, 00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00019 // GNU Library General Public License for more details. 00020 // 00021 // You should have received a copy of the GNU Library General Public License 00022 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to 00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 00024 // 00025 // The U.S. Government is granted a limited license as per AL 91-7. 00026 // 00027 00028 #ifndef _chemistry_qc_dft_integrator_h 00029 #define _chemistry_qc_dft_integrator_h 00030 00031 #ifdef __GNUC__ 00032 #pragma interface 00033 #endif 00034 00035 #include <util/state/state.h> 00036 #include <util/group/thread.h> 00037 #include <chemistry/qc/dft/functional.h> 00038 #include <chemistry/qc/basis/extent.h> 00039 #include <chemistry/qc/wfn/density.h> 00040 00041 namespace sc { 00042 00044 class DenIntegrator: virtual public SavableState { 00045 protected: 00046 Ref<Wavefunction> wfn_; 00047 //clj Ref<ShellExtent> extent_; 00048 Ref<BatchElectronDensity> den_; 00049 00050 Ref<ThreadGrp> threadgrp_; 00051 Ref<MessageGrp> messagegrp_; 00052 00053 double value_; 00054 double accuracy_; 00055 00056 double *alpha_vmat_; 00057 double *beta_vmat_; 00058 00059 //clj double *alpha_dmat_; 00060 //clj double *beta_dmat_; 00061 //clj double *dmat_bound_; 00062 00063 int spin_polarized_; 00064 00065 int need_density_; // specialization must set to 1 if it needs density_ 00066 double density_; 00067 int nbasis_; 00068 int nshell_; 00069 int n_integration_center_; 00070 int natom_; 00071 int compute_potential_integrals_; // 1 if potential integrals are needed 00072 00073 int linear_scaling_; 00074 int use_dmat_bound_; 00075 00076 void init_integration(const Ref<DenFunctional> &func, 00077 const RefSymmSCMatrix& densa, 00078 const RefSymmSCMatrix& densb, 00079 double *nuclear_gradient); 00080 void done_integration(); 00081 00082 void init_object(); 00083 public: 00085 DenIntegrator(); 00087 DenIntegrator(const Ref<KeyVal> &); 00089 DenIntegrator(StateIn &); 00090 ~DenIntegrator(); 00091 void save_data_state(StateOut &); 00092 00094 Ref<Wavefunction> wavefunction() const { return wfn_; } 00096 double value() const { return value_; } 00097 00099 void set_accuracy(double a); 00100 double get_accuracy(void) {return accuracy_; } 00103 void set_compute_potential_integrals(int); 00106 const double *alpha_vmat() const { return alpha_vmat_; } 00109 const double *beta_vmat() const { return beta_vmat_; } 00110 00113 virtual void init(const Ref<Wavefunction> &); 00115 virtual void done(); 00119 virtual void integrate(const Ref<DenFunctional> &, 00120 const RefSymmSCMatrix& densa =0, 00121 const RefSymmSCMatrix& densb =0, 00122 double *nuclear_grad = 0) = 0; 00123 }; 00124 00125 00127 class IntegrationWeight: virtual public SavableState { 00128 00129 protected: 00130 00131 Ref<Molecule> mol_; 00132 double tol_; 00133 00134 void fd_w(int icenter, SCVector3 &point, double *fd_grad_w); 00135 00136 public: 00137 IntegrationWeight(); 00138 IntegrationWeight(const Ref<KeyVal> &); 00139 IntegrationWeight(StateIn &); 00140 ~IntegrationWeight(); 00141 void save_data_state(StateOut &); 00142 00143 void test(int icenter, SCVector3 &point); 00144 void test(); 00145 00147 virtual void init(const Ref<Molecule> &, double tolerance); 00149 virtual void done(); 00154 virtual double w(int center, SCVector3 &point, double *grad_w = 0) = 0; 00155 }; 00156 00157 00159 class BeckeIntegrationWeight: public IntegrationWeight { 00160 00161 int n_integration_centers; 00162 SCVector3 *centers; 00163 double *atomic_radius; 00164 00165 double **a_mat; 00166 double **oorab; 00167 00168 void compute_grad_p(int gc, int ic, int wc, SCVector3&r, double p, 00169 SCVector3&g); 00170 void compute_grad_nu(int gc, int bc, SCVector3 &point, SCVector3 &grad); 00171 00172 double compute_t(int ic, int jc, SCVector3 &point); 00173 double compute_p(int icenter, SCVector3&point); 00174 00175 public: 00176 BeckeIntegrationWeight(); 00177 BeckeIntegrationWeight(const Ref<KeyVal> &); 00178 BeckeIntegrationWeight(StateIn &); 00179 ~BeckeIntegrationWeight(); 00180 void save_data_state(StateOut &); 00181 00182 void init(const Ref<Molecule> &, double tolerance); 00183 void done(); 00184 double w(int center, SCVector3 &point, double *grad_w = 0); 00185 }; 00186 00188 class RadialIntegrator: virtual public SavableState { 00189 protected: 00190 int nr_; 00191 public: 00192 RadialIntegrator(); 00193 RadialIntegrator(const Ref<KeyVal> &); 00194 RadialIntegrator(StateIn &); 00195 ~RadialIntegrator(); 00196 void save_data_state(StateOut &); 00197 00198 virtual int nr() const = 0; 00199 virtual double radial_value(int ir, int nr, double radii, 00200 double &multiplier) = 0; 00201 }; 00202 00203 00205 class AngularIntegrator: virtual public SavableState{ 00206 protected: 00207 public: 00208 AngularIntegrator(); 00209 AngularIntegrator(const Ref<KeyVal> &); 00210 AngularIntegrator(StateIn &); 00211 ~AngularIntegrator(); 00212 void save_data_state(StateOut &); 00213 00214 virtual int nw(void) const = 0; 00215 virtual int num_angular_points(double r_value, int ir) = 0; 00216 virtual double angular_point_cartesian(int iangular, double r, 00217 SCVector3 &integration_point) const = 0; 00218 }; 00219 00220 00223 class EulerMaclaurinRadialIntegrator: public RadialIntegrator { 00224 public: 00225 EulerMaclaurinRadialIntegrator(); 00226 EulerMaclaurinRadialIntegrator(int i); 00230 EulerMaclaurinRadialIntegrator(const Ref<KeyVal> &); 00231 EulerMaclaurinRadialIntegrator(StateIn &); 00232 ~EulerMaclaurinRadialIntegrator(); 00233 void save_data_state(StateOut &); 00234 00235 int nr() const; 00236 double radial_value(int ir, int nr, double radii, double &multiplier); 00237 00238 void print(std::ostream & =ExEnv::out0()) const; 00239 }; 00240 00282 class LebedevLaikovIntegrator: public AngularIntegrator { 00283 protected: 00284 int npoint_; 00285 double *x_, *y_, *z_, *w_; 00286 00287 void init(int n); 00288 public: 00289 LebedevLaikovIntegrator(); 00295 LebedevLaikovIntegrator(const Ref<KeyVal> &); 00296 LebedevLaikovIntegrator(StateIn &); 00297 LebedevLaikovIntegrator(int); 00298 ~LebedevLaikovIntegrator(); 00299 void save_data_state(StateOut &); 00300 00301 int nw(void) const; 00302 int num_angular_points(double r_value, int ir); 00303 double angular_point_cartesian(int iangular, double r, 00304 SCVector3 &integration_point) const; 00305 void print(std::ostream & =ExEnv::out0()) const; 00306 }; 00307 00310 class GaussLegendreAngularIntegrator: public AngularIntegrator { 00311 protected: 00312 int ntheta_; 00313 int nphi_; 00314 int Ktheta_; 00315 int ntheta_r_; 00316 int nphi_r_; 00317 int Ktheta_r_; 00318 double *theta_quad_weights_; 00319 double *theta_quad_points_; 00320 00321 int get_ntheta(void) const; 00322 void set_ntheta(int i); 00323 int get_nphi(void) const; 00324 void set_nphi(int i); 00325 int get_Ktheta(void) const; 00326 void set_Ktheta(int i); 00327 int get_ntheta_r(void) const; 00328 void set_ntheta_r(int i); 00329 int get_nphi_r(void) const; 00330 void set_nphi_r(int i); 00331 int get_Ktheta_r(void) const; 00332 void set_Ktheta_r(int i); 00333 int nw(void) const; 00334 double sin_theta(SCVector3 &point) const; 00335 void gauleg(double x1, double x2, int n); 00336 public: 00337 GaussLegendreAngularIntegrator(); 00344 GaussLegendreAngularIntegrator(const Ref<KeyVal> &); 00345 GaussLegendreAngularIntegrator(StateIn &); 00346 ~GaussLegendreAngularIntegrator(); 00347 void save_data_state(StateOut &); 00348 00349 int num_angular_points(double r_value, int ir); 00350 double angular_point_cartesian(int iangular, double r, 00351 SCVector3 &integration_point) const; 00352 void print(std::ostream & =ExEnv::out0()) const; 00353 }; 00354 00357 class RadialAngularIntegrator: public DenIntegrator { 00358 private: 00359 int prune_grid_; 00360 double **Alpha_coeffs_; 00361 int gridtype_; 00362 int **nr_points_, *xcoarse_l_; 00363 int npruned_partitions_; 00364 double *grid_accuracy_; 00365 int dynamic_grids_; 00366 int natomic_rows_; 00367 int max_gridtype_; 00368 protected: 00369 Ref<IntegrationWeight> weight_; 00370 Ref<RadialIntegrator> radial_user_; 00371 Ref<AngularIntegrator> angular_user_; 00372 Ref<AngularIntegrator> ***angular_grid_; 00373 Ref<RadialIntegrator> **radial_grid_; 00374 public: 00375 RadialAngularIntegrator(); 00416 RadialAngularIntegrator(const Ref<KeyVal> &); 00417 RadialAngularIntegrator(StateIn &); 00418 ~RadialAngularIntegrator(); 00419 void save_data_state(StateOut &); 00420 00421 void integrate(const Ref<DenFunctional> &, 00422 const RefSymmSCMatrix& densa =0, 00423 const RefSymmSCMatrix& densb =0, 00424 double *nuclear_gradient = 0); 00425 void print(std::ostream & =ExEnv::out0()) const; 00426 AngularIntegrator *get_angular_grid(double radius, double atomic_radius, 00427 int charge, int deriv_order); 00428 RadialIntegrator *get_radial_grid(int charge, int deriv_order); 00429 void init_default_grids(void); 00430 int angular_grid_offset(int i); 00431 void set_grids(void); 00432 int get_atomic_row(int i); 00433 void init_parameters(void); 00434 void init_parameters(const Ref<KeyVal>& keyval); 00435 void init_pruning_coefficients(const Ref<KeyVal>& keyval); 00436 void init_pruning_coefficients(void); 00437 void init_alpha_coefficients(void); 00438 int select_dynamic_grid(void); 00439 Ref<IntegrationWeight> weight() { return weight_; } 00440 }; 00441 00442 } 00443 00444 #endif 00445 00446 // Local Variables: 00447 // mode: c++ 00448 // c-file-style: "CLJ" 00449 // End: