[0b990d] | 1 | //
|
---|
| 2 | // integrator.cc
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 1997 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 | #ifdef __GNUC__
|
---|
| 29 | #pragma implementation
|
---|
| 30 | #endif
|
---|
| 31 |
|
---|
| 32 | #include <util/misc/math.h>
|
---|
| 33 |
|
---|
| 34 | #include <util/misc/timer.h>
|
---|
| 35 | #include <util/misc/formio.h>
|
---|
| 36 | #include <util/state/stateio.h>
|
---|
| 37 | #include <util/container/carray.h>
|
---|
| 38 | #include <chemistry/qc/dft/integrator.h>
|
---|
| 39 |
|
---|
| 40 | using namespace std;
|
---|
| 41 | using namespace sc;
|
---|
| 42 |
|
---|
| 43 | namespace sc {
|
---|
| 44 |
|
---|
| 45 | #ifdef EXPLICIT_TEMPLATE_INSTANTIATION
|
---|
| 46 | template void delete_c_array2<Ref<RadialIntegrator> >(Ref<RadialIntegrator>**);
|
---|
| 47 | template Ref<RadialIntegrator>**
|
---|
| 48 | new_c_array2<Ref<RadialIntegrator> >(int,int,Ref<RadialIntegrator>);
|
---|
| 49 | template void delete_c_array3<Ref<RadialIntegrator> >(Ref<RadialIntegrator>***);
|
---|
| 50 | template Ref<RadialIntegrator>***
|
---|
| 51 | new_c_array3<Ref<RadialIntegrator> >(int,int,int,Ref<RadialIntegrator>);
|
---|
| 52 |
|
---|
| 53 | template void delete_c_array2<Ref<AngularIntegrator> >(Ref<AngularIntegrator>**);
|
---|
| 54 | template Ref<AngularIntegrator>**
|
---|
| 55 | new_c_array2<Ref<AngularIntegrator> >(int,int,Ref<AngularIntegrator>);
|
---|
| 56 | template void delete_c_array3<Ref<AngularIntegrator> >(Ref<AngularIntegrator>***);
|
---|
| 57 | template Ref<AngularIntegrator>***
|
---|
| 58 | new_c_array3<Ref<AngularIntegrator> >(int,int,int,Ref<AngularIntegrator>);
|
---|
| 59 | #endif
|
---|
| 60 |
|
---|
| 61 | //#define CHECK_ALIGN(v) if(int(&v)&7)ExEnv::outn()<<"Bad Alignment: "<< ## v <<endl;
|
---|
| 62 | #define CHECK_ALIGN(v)
|
---|
| 63 |
|
---|
| 64 | ///////////////////////////////////////////////////////////////////////////
|
---|
| 65 | // utility functions
|
---|
| 66 |
|
---|
| 67 | inline static double
|
---|
| 68 | norm(double v[3])
|
---|
| 69 | {
|
---|
| 70 | double x,y,z;
|
---|
| 71 | return sqrt((x=v[0])*x + (y=v[1])*y + (z=v[2])*z);
|
---|
| 72 | }
|
---|
| 73 |
|
---|
| 74 | static double
|
---|
| 75 | get_radius(const Ref<Molecule> &mol, int iatom)
|
---|
| 76 | {
|
---|
| 77 | double r = mol->atominfo()->maxprob_radius(mol->Z(iatom));
|
---|
| 78 | if (r == 0) {
|
---|
| 79 | static bool warned = false;
|
---|
| 80 | if (!warned) {
|
---|
| 81 | ExEnv::out0()
|
---|
| 82 | << indent << "WARNING: BeckeIntegrationWeight usually uses the atomic maximum" << std::endl
|
---|
| 83 | << indent << " probability radius, however this is not available for" << std::endl
|
---|
| 84 | << indent << " one of the atoms in your system. The Bragg radius will" << std::endl
|
---|
| 85 | << indent << " be used instead." << std::endl;
|
---|
| 86 | warned = true;
|
---|
| 87 | }
|
---|
| 88 | r = mol->atominfo()->bragg_radius(mol->Z(iatom));
|
---|
| 89 | }
|
---|
| 90 | if (r == 0) {
|
---|
| 91 | ExEnv::out0()
|
---|
| 92 | << indent << "ERROR: BeckeIntegrationWeight could not find a maximum probability" << std::endl
|
---|
| 93 | << indent << " or a Bragg radius for an atom" << std::endl;
|
---|
| 94 | abort();
|
---|
| 95 | }
|
---|
| 96 | return r;
|
---|
| 97 | }
|
---|
| 98 |
|
---|
| 99 | ///////////////////////////////////////////////////////////////////////////
|
---|
| 100 | // DenIntegratorThread
|
---|
| 101 |
|
---|
| 102 | //ThreadLock *tlock;
|
---|
| 103 |
|
---|
| 104 | class DenIntegratorThread: public Thread {
|
---|
| 105 | protected:
|
---|
| 106 | // data common to all threads
|
---|
| 107 | int nthread_;
|
---|
| 108 | int nshell_;
|
---|
| 109 | int nbasis_;
|
---|
| 110 | int natom_;
|
---|
| 111 | int n_integration_center_;
|
---|
| 112 | DenIntegrator *integrator_;
|
---|
| 113 | int spin_polarized_;
|
---|
| 114 | int need_hessian_;
|
---|
| 115 | int need_gradient_;
|
---|
| 116 | GaussianBasisSet *basis_;
|
---|
| 117 | bool linear_scaling_;
|
---|
| 118 | bool use_dmat_bound_;
|
---|
| 119 | double value_;
|
---|
| 120 | DenFunctional *func_;
|
---|
| 121 |
|
---|
| 122 | double accuracy_;
|
---|
| 123 | int compute_potential_integrals_;
|
---|
| 124 |
|
---|
| 125 | // data local to thread
|
---|
| 126 | int ithread_;
|
---|
| 127 | Ref<BatchElectronDensity> den_;
|
---|
| 128 | double *alpha_vmat_; // lower triangle of xi_i(r) v(r) xi_j(r) integrals
|
---|
| 129 | double *beta_vmat_; // lower triangle of xi_i(r) v(r) xi_j(r) integrals
|
---|
| 130 | double *nuclear_gradient_;
|
---|
| 131 | double *w_gradient_;
|
---|
| 132 | double *f_gradient_;
|
---|
| 133 |
|
---|
| 134 | public:
|
---|
| 135 | DenIntegratorThread(int ithread, int nthread,
|
---|
| 136 | DenIntegrator *integrator,
|
---|
| 137 | DenFunctional *func,
|
---|
| 138 | const Ref<BatchElectronDensity> &den,
|
---|
| 139 | int linear_scaling,
|
---|
| 140 | int use_dmat_bound,
|
---|
| 141 | double accuracy,
|
---|
| 142 | int compute_potential_integrals,
|
---|
| 143 | int need_nuclear_gradient);
|
---|
| 144 | virtual ~DenIntegratorThread();
|
---|
| 145 | double do_point(int iatom, const SCVector3 &r,
|
---|
| 146 | double weight, double multiplier,
|
---|
| 147 | double *nuclear_gradient,
|
---|
| 148 | double *f_gradient, double *w_gradient);
|
---|
| 149 | double *nuclear_gradient() { return nuclear_gradient_; }
|
---|
| 150 | double *alpha_vmat() { return alpha_vmat_; }
|
---|
| 151 | double *beta_vmat() { return beta_vmat_; }
|
---|
| 152 | double value() { return value_; }
|
---|
| 153 | };
|
---|
| 154 |
|
---|
| 155 | DenIntegratorThread::DenIntegratorThread(int ithread, int nthread,
|
---|
| 156 | DenIntegrator *integrator,
|
---|
| 157 | DenFunctional *func,
|
---|
| 158 | const Ref<BatchElectronDensity> &den,
|
---|
| 159 | int linear_scaling,
|
---|
| 160 | int use_dmat_bound,
|
---|
| 161 | double accuracy,
|
---|
| 162 | int compute_potential_integrals,
|
---|
| 163 | int need_nuclear_gradient)
|
---|
| 164 | {
|
---|
| 165 | value_ = 0.0;
|
---|
| 166 |
|
---|
| 167 | den_ = den;
|
---|
| 168 | ithread_ = ithread;
|
---|
| 169 | nthread_ = nthread;
|
---|
| 170 | integrator_ = integrator;
|
---|
| 171 | spin_polarized_ = integrator->wavefunction()->spin_polarized();
|
---|
| 172 | nshell_ = integrator->wavefunction()->basis()->nshell();
|
---|
| 173 | nbasis_ = integrator->wavefunction()->basis()->nbasis();
|
---|
| 174 | natom_ = integrator->wavefunction()->molecule()->natom();
|
---|
| 175 | n_integration_center_
|
---|
| 176 | = integrator->wavefunction()->molecule()->n_non_q_atom();
|
---|
| 177 | need_gradient_ = func->need_density_gradient();
|
---|
| 178 | need_hessian_ = func->need_density_hessian();
|
---|
| 179 | basis_ = integrator->wavefunction()->basis().pointer();
|
---|
| 180 | linear_scaling_ = linear_scaling;
|
---|
| 181 | use_dmat_bound_ = use_dmat_bound;
|
---|
| 182 | func_ = func;
|
---|
| 183 | accuracy_ = accuracy;
|
---|
| 184 | compute_potential_integrals_ = compute_potential_integrals;
|
---|
| 185 | den_->set_accuracy(accuracy);
|
---|
| 186 |
|
---|
| 187 | if (need_nuclear_gradient) {
|
---|
| 188 | nuclear_gradient_ = new double[3*natom_];
|
---|
| 189 | memset(nuclear_gradient_, 0, 3*natom_*sizeof(double));
|
---|
| 190 | }
|
---|
| 191 | else nuclear_gradient_ = 0;
|
---|
| 192 |
|
---|
| 193 | alpha_vmat_ = 0;
|
---|
| 194 | beta_vmat_ = 0;
|
---|
| 195 | if (compute_potential_integrals_) {
|
---|
| 196 | int ntri = (nbasis_*(nbasis_+1))/2;
|
---|
| 197 | alpha_vmat_ = new double[ntri];
|
---|
| 198 | memset(alpha_vmat_, 0, sizeof(double)*ntri);
|
---|
| 199 | if (spin_polarized_) {
|
---|
| 200 | beta_vmat_ = new double[ntri];
|
---|
| 201 | memset(beta_vmat_, 0, sizeof(double)*ntri);
|
---|
| 202 | }
|
---|
| 203 | }
|
---|
| 204 |
|
---|
| 205 | w_gradient_ = 0;
|
---|
| 206 | f_gradient_ = 0;
|
---|
| 207 | if (nuclear_gradient_) {
|
---|
| 208 | w_gradient_ = new double[n_integration_center_*3];
|
---|
| 209 | f_gradient_ = new double[natom_*3];
|
---|
| 210 | }
|
---|
| 211 | }
|
---|
| 212 |
|
---|
| 213 | DenIntegratorThread::~DenIntegratorThread()
|
---|
| 214 | {
|
---|
| 215 | delete[] alpha_vmat_;
|
---|
| 216 | delete[] beta_vmat_;
|
---|
| 217 | delete[] nuclear_gradient_;
|
---|
| 218 | delete[] f_gradient_;
|
---|
| 219 | delete[] w_gradient_;
|
---|
| 220 | }
|
---|
| 221 |
|
---|
| 222 | ///////////////////////////////////////////////////////////////////////////
|
---|
| 223 | // DenIntegrator
|
---|
| 224 |
|
---|
| 225 | static ClassDesc DenIntegrator_cd(
|
---|
| 226 | typeid(DenIntegrator),"DenIntegrator",1,"public SavableState",
|
---|
| 227 | 0, 0, 0);
|
---|
| 228 |
|
---|
| 229 | DenIntegrator::DenIntegrator(StateIn& s):
|
---|
| 230 | SavableState(s)
|
---|
| 231 | {
|
---|
| 232 | init_object();
|
---|
| 233 | s.get(linear_scaling_);
|
---|
| 234 | s.get(use_dmat_bound_);
|
---|
| 235 | }
|
---|
| 236 |
|
---|
| 237 | DenIntegrator::DenIntegrator()
|
---|
| 238 | {
|
---|
| 239 | init_object();
|
---|
| 240 | }
|
---|
| 241 |
|
---|
| 242 | DenIntegrator::DenIntegrator(const Ref<KeyVal>& keyval)
|
---|
| 243 | {
|
---|
| 244 | init_object();
|
---|
| 245 |
|
---|
| 246 | linear_scaling_ = keyval->booleanvalue("linear_scaling",
|
---|
| 247 | KeyValValueboolean(linear_scaling_));
|
---|
| 248 | use_dmat_bound_ = keyval->booleanvalue("use_dmat_bound",
|
---|
| 249 | KeyValValueboolean(use_dmat_bound_));
|
---|
| 250 | }
|
---|
| 251 |
|
---|
| 252 | DenIntegrator::~DenIntegrator()
|
---|
| 253 | {
|
---|
| 254 | }
|
---|
| 255 |
|
---|
| 256 | void
|
---|
| 257 | DenIntegrator::save_data_state(StateOut& s)
|
---|
| 258 | {
|
---|
| 259 | s.put(linear_scaling_);
|
---|
| 260 | s.put(use_dmat_bound_);
|
---|
| 261 | }
|
---|
| 262 |
|
---|
| 263 | void
|
---|
| 264 | DenIntegrator::init_object()
|
---|
| 265 | {
|
---|
| 266 | threadgrp_ = ThreadGrp::get_default_threadgrp();
|
---|
| 267 | messagegrp_ = MessageGrp::get_default_messagegrp();
|
---|
| 268 | compute_potential_integrals_ = 0;
|
---|
| 269 | accuracy_ = DBL_EPSILON;
|
---|
| 270 | linear_scaling_ = 1;
|
---|
| 271 | use_dmat_bound_ = 1;
|
---|
| 272 | alpha_vmat_ = 0;
|
---|
| 273 | beta_vmat_ = 0;
|
---|
| 274 | }
|
---|
| 275 |
|
---|
| 276 | void
|
---|
| 277 | DenIntegrator::set_compute_potential_integrals(int i)
|
---|
| 278 | {
|
---|
| 279 | compute_potential_integrals_=i;
|
---|
| 280 | }
|
---|
| 281 |
|
---|
| 282 | void
|
---|
| 283 | DenIntegrator::init(const Ref<Wavefunction> &wfn)
|
---|
| 284 | {
|
---|
| 285 | wfn_ = wfn;
|
---|
| 286 | den_ = new BatchElectronDensity(wfn,accuracy_);
|
---|
| 287 | den_->set_linear_scaling(linear_scaling_);
|
---|
| 288 | den_->set_use_dmat_bound(use_dmat_bound_);
|
---|
| 289 | den_->init(false);
|
---|
| 290 | }
|
---|
| 291 |
|
---|
| 292 | void
|
---|
| 293 | DenIntegrator::set_accuracy(double a)
|
---|
| 294 | {
|
---|
| 295 | accuracy_ = a;
|
---|
| 296 | if (den_.nonnull()) den_->set_accuracy(a);
|
---|
| 297 | }
|
---|
| 298 |
|
---|
| 299 | void
|
---|
| 300 | DenIntegrator::done()
|
---|
| 301 | {
|
---|
| 302 | wfn_ = 0;
|
---|
| 303 | den_ = 0;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | void
|
---|
| 307 | DenIntegrator::init_integration(const Ref<DenFunctional> &func,
|
---|
| 308 | const RefSymmSCMatrix& densa,
|
---|
| 309 | const RefSymmSCMatrix& densb,
|
---|
| 310 | double *nuclear_gradient)
|
---|
| 311 | {
|
---|
| 312 | int i;
|
---|
| 313 | value_ = 0.0;
|
---|
| 314 |
|
---|
| 315 | func->set_compute_potential(
|
---|
| 316 | compute_potential_integrals_ || nuclear_gradient != 0);
|
---|
| 317 |
|
---|
| 318 | spin_polarized_ = wfn_->spin_polarized();
|
---|
| 319 | func->set_spin_polarized(spin_polarized_);
|
---|
| 320 |
|
---|
| 321 | natom_ = wfn_->molecule()->natom();
|
---|
| 322 | n_integration_center_
|
---|
| 323 | = wfn_->molecule()->n_non_q_atom();
|
---|
| 324 |
|
---|
| 325 | nshell_ = wfn_->basis()->nshell();
|
---|
| 326 | nbasis_ = wfn_->basis()->nbasis();
|
---|
| 327 |
|
---|
| 328 | den_->set_densities(densa,densb);
|
---|
| 329 |
|
---|
| 330 | delete[] alpha_vmat_;
|
---|
| 331 | delete[] beta_vmat_;
|
---|
| 332 | alpha_vmat_ = 0;
|
---|
| 333 | beta_vmat_ = 0;
|
---|
| 334 | if (compute_potential_integrals_) {
|
---|
| 335 | int ntri = (nbasis_*(nbasis_+1))/2;
|
---|
| 336 | alpha_vmat_ = new double[ntri];
|
---|
| 337 | memset(alpha_vmat_, 0, sizeof(double)*ntri);
|
---|
| 338 | if (spin_polarized_) {
|
---|
| 339 | beta_vmat_ = new double[ntri];
|
---|
| 340 | memset(beta_vmat_, 0, sizeof(double)*ntri);
|
---|
| 341 | }
|
---|
| 342 | }
|
---|
| 343 | }
|
---|
| 344 |
|
---|
| 345 | void
|
---|
| 346 | DenIntegrator::done_integration()
|
---|
| 347 | {
|
---|
| 348 | messagegrp_->sum(value_);
|
---|
| 349 | if (compute_potential_integrals_) {
|
---|
| 350 | int ntri = (nbasis_*(nbasis_+1))/2;
|
---|
| 351 | messagegrp_->sum(alpha_vmat_,ntri);
|
---|
| 352 | if (spin_polarized_) {
|
---|
| 353 | messagegrp_->sum(beta_vmat_,ntri);
|
---|
| 354 | }
|
---|
| 355 | }
|
---|
| 356 | }
|
---|
| 357 |
|
---|
| 358 | double
|
---|
| 359 | DenIntegratorThread::do_point(int iatom, const SCVector3 &r,
|
---|
| 360 | double weight, double multiplier,
|
---|
| 361 | double *nuclear_gradient,
|
---|
| 362 | double *f_gradient, double *w_gradient)
|
---|
| 363 | {
|
---|
| 364 | int i,j;
|
---|
| 365 | double w_mult = weight * multiplier;
|
---|
| 366 |
|
---|
| 367 | CHECK_ALIGN(w_mult);
|
---|
| 368 |
|
---|
| 369 | PointInputData id(r);
|
---|
| 370 |
|
---|
| 371 | den_->compute_density(r,
|
---|
| 372 | &id.a.rho,
|
---|
| 373 | (need_gradient_?id.a.del_rho:0),
|
---|
| 374 | (need_hessian_?id.a.hes_rho:0),
|
---|
| 375 | &id.b.rho,
|
---|
| 376 | (need_gradient_?id.b.del_rho:0),
|
---|
| 377 | (need_hessian_?id.b.hes_rho:0));
|
---|
| 378 |
|
---|
| 379 | id.compute_derived(spin_polarized_, need_gradient_, need_hessian_);
|
---|
| 380 |
|
---|
| 381 | int ncontrib = den_->ncontrib();
|
---|
| 382 | int *contrib = den_->contrib();
|
---|
| 383 | int ncontrib_bf = den_->ncontrib_bf();
|
---|
| 384 | int *contrib_bf = den_->contrib_bf();
|
---|
| 385 | double *bs_values = den_->bs_values();
|
---|
| 386 | double *bsg_values = den_->bsg_values();
|
---|
| 387 | double *bsh_values = den_->bsh_values();
|
---|
| 388 |
|
---|
| 389 | PointOutputData od;
|
---|
| 390 | if ( (id.a.rho + id.b.rho) > 1e2*DBL_EPSILON) {
|
---|
| 391 | if (nuclear_gradient == 0) {
|
---|
| 392 | func_->point(id, od);
|
---|
| 393 | }
|
---|
| 394 | else {
|
---|
| 395 | func_->gradient(id, od, f_gradient, iatom, basis_,
|
---|
| 396 | den_->alpha_density_matrix(),
|
---|
| 397 | den_->beta_density_matrix(),
|
---|
| 398 | ncontrib, contrib,
|
---|
| 399 | ncontrib_bf, contrib_bf,
|
---|
| 400 | bs_values, bsg_values,
|
---|
| 401 | bsh_values);
|
---|
| 402 | }
|
---|
| 403 | }
|
---|
| 404 | else {
|
---|
| 405 | return id.a.rho + id.b.rho;
|
---|
| 406 | }
|
---|
| 407 |
|
---|
| 408 | value_ += od.energy * w_mult;
|
---|
| 409 |
|
---|
| 410 | if (compute_potential_integrals_) {
|
---|
| 411 | // the contribution to the potential integrals
|
---|
| 412 | if (need_gradient_) {
|
---|
| 413 | double gradsa[3], gradsb[3];
|
---|
| 414 | gradsa[0] = w_mult*(2.0*od.df_dgamma_aa*id.a.del_rho[0] +
|
---|
| 415 | od.df_dgamma_ab*id.b.del_rho[0]);
|
---|
| 416 | gradsa[1] = w_mult*(2.0*od.df_dgamma_aa*id.a.del_rho[1] +
|
---|
| 417 | od.df_dgamma_ab*id.b.del_rho[1]);
|
---|
| 418 | gradsa[2] = w_mult*(2.0*od.df_dgamma_aa*id.a.del_rho[2] +
|
---|
| 419 | od.df_dgamma_ab*id.b.del_rho[2]);
|
---|
| 420 | double drhoa = w_mult*od.df_drho_a, drhob=0.0;
|
---|
| 421 | if (spin_polarized_) {
|
---|
| 422 | drhob = w_mult*od.df_drho_b;
|
---|
| 423 | gradsb[0] = w_mult*(2.0*od.df_dgamma_bb*id.b.del_rho[0] +
|
---|
| 424 | od.df_dgamma_ab*id.a.del_rho[0]);
|
---|
| 425 | gradsb[1] = w_mult*(2.0*od.df_dgamma_bb*id.b.del_rho[1] +
|
---|
| 426 | od.df_dgamma_ab*id.a.del_rho[1]);
|
---|
| 427 | gradsb[2] = w_mult*(2.0*od.df_dgamma_bb*id.b.del_rho[2] +
|
---|
| 428 | od.df_dgamma_ab*id.a.del_rho[2]);
|
---|
| 429 | }
|
---|
| 430 |
|
---|
| 431 | for (int j=0; j<ncontrib_bf; j++) {
|
---|
| 432 | int jt = contrib_bf[j];
|
---|
| 433 | double dfdra_phi_m = drhoa*bs_values[j];
|
---|
| 434 | double dfdga_phi_m = gradsa[0]*bsg_values[j*3+0] +
|
---|
| 435 | gradsa[1]*bsg_values[j*3+1] +
|
---|
| 436 | gradsa[2]*bsg_values[j*3+2];
|
---|
| 437 | double vamu = dfdra_phi_m + dfdga_phi_m, vbmu=0.0;
|
---|
| 438 | double dfdrb_phi_m, dfdgb_phi_m;
|
---|
| 439 | if (spin_polarized_) {
|
---|
| 440 | dfdrb_phi_m = drhob*bs_values[j];
|
---|
| 441 | dfdgb_phi_m = gradsb[0]*bsg_values[j*3+0] +
|
---|
| 442 | gradsb[1]*bsg_values[j*3+1] +
|
---|
| 443 | gradsb[2]*bsg_values[j*3+2];
|
---|
| 444 | vbmu = dfdrb_phi_m + dfdgb_phi_m;
|
---|
| 445 | }
|
---|
| 446 |
|
---|
| 447 | int jtoff = (jt*(jt+1))>>1;
|
---|
| 448 |
|
---|
| 449 | for (int k=0; k <= j; k++) {
|
---|
| 450 | int kt = contrib_bf[k];
|
---|
| 451 | int jtkt = jtoff + kt;
|
---|
| 452 |
|
---|
| 453 | double dfdga_phi_n = gradsa[0]*bsg_values[k*3+0] +
|
---|
| 454 | gradsa[1]*bsg_values[k*3+1] +
|
---|
| 455 | gradsa[2]*bsg_values[k*3+2];
|
---|
| 456 | alpha_vmat_[jtkt] += vamu * bs_values[k] +
|
---|
| 457 | dfdga_phi_n * bs_values[j];
|
---|
| 458 | if (spin_polarized_) {
|
---|
| 459 | double dfdgb_phi_n = gradsb[0]*bsg_values[k*3+0] +
|
---|
| 460 | gradsb[1]*bsg_values[k*3+1] +
|
---|
| 461 | gradsb[2]*bsg_values[k*3+2];
|
---|
| 462 | beta_vmat_[jtkt] += vbmu * bs_values[k] +
|
---|
| 463 | dfdgb_phi_n * bs_values[j];
|
---|
| 464 | }
|
---|
| 465 | }
|
---|
| 466 | }
|
---|
| 467 | }
|
---|
| 468 | else {
|
---|
| 469 | double drhoa = w_mult*od.df_drho_a;
|
---|
| 470 | double drhob = w_mult*od.df_drho_b;
|
---|
| 471 | for (int j=0; j<ncontrib_bf; j++) {
|
---|
| 472 | int jt = contrib_bf[j];
|
---|
| 473 | double bsj = bs_values[j];
|
---|
| 474 | double dfa_phi_m = drhoa * bsj;
|
---|
| 475 | double dfb_phi_m = drhob * bsj;
|
---|
| 476 | int jtoff = (jt*(jt+1))>>1;
|
---|
| 477 | for (int k=0; k <= j; k++) {
|
---|
| 478 | int kt = contrib_bf[k];
|
---|
| 479 | int jtkt = jtoff + kt;
|
---|
| 480 | double bsk = bs_values[k];
|
---|
| 481 | alpha_vmat_[jtkt] += dfa_phi_m * bsk;
|
---|
| 482 | if (spin_polarized_)
|
---|
| 483 | beta_vmat_[jtkt] += dfb_phi_m * bsk;
|
---|
| 484 | }
|
---|
| 485 | }
|
---|
| 486 | }
|
---|
| 487 | }
|
---|
| 488 |
|
---|
| 489 | if (nuclear_gradient != 0) {
|
---|
| 490 | // the contribution from f dw/dx
|
---|
| 491 | if (w_gradient) {
|
---|
| 492 | for (int icenter = 0; icenter<n_integration_center_; icenter++) {
|
---|
| 493 | int iatom = basis_->molecule()->non_q_atom(icenter);
|
---|
| 494 | for (int ixyz=0; ixyz<3; ixyz++) {
|
---|
| 495 | nuclear_gradient[iatom*3+ixyz]
|
---|
| 496 | += w_gradient[icenter*3+ixyz] * od.energy * multiplier;
|
---|
| 497 | }
|
---|
| 498 | }
|
---|
| 499 | }
|
---|
| 500 | // the contribution from (df/dx) w
|
---|
| 501 | for (i=0; i<natom_*3; i++) {
|
---|
| 502 | nuclear_gradient[i] += f_gradient[i] * w_mult;
|
---|
| 503 | }
|
---|
| 504 | }
|
---|
| 505 |
|
---|
| 506 | return id.a.rho + id.b.rho;
|
---|
| 507 | }
|
---|
| 508 |
|
---|
| 509 | ///////////////////////////////////////////////////////////////////////////
|
---|
| 510 | // IntegrationWeight
|
---|
| 511 |
|
---|
| 512 | static ClassDesc IntegrationWeight_cd(
|
---|
| 513 | typeid(IntegrationWeight),"IntegrationWeight",1,"public SavableState",
|
---|
| 514 | 0, 0, 0);
|
---|
| 515 |
|
---|
| 516 | IntegrationWeight::IntegrationWeight(StateIn& s):
|
---|
| 517 | SavableState(s)
|
---|
| 518 | {
|
---|
| 519 | }
|
---|
| 520 |
|
---|
| 521 | IntegrationWeight::IntegrationWeight()
|
---|
| 522 | {
|
---|
| 523 | }
|
---|
| 524 |
|
---|
| 525 | IntegrationWeight::IntegrationWeight(const Ref<KeyVal>& keyval)
|
---|
| 526 | {
|
---|
| 527 | }
|
---|
| 528 |
|
---|
| 529 | IntegrationWeight::~IntegrationWeight()
|
---|
| 530 | {
|
---|
| 531 | }
|
---|
| 532 |
|
---|
| 533 | void
|
---|
| 534 | IntegrationWeight::save_data_state(StateOut& s)
|
---|
| 535 | {
|
---|
| 536 | }
|
---|
| 537 |
|
---|
| 538 | void
|
---|
| 539 | IntegrationWeight::init(const Ref<Molecule> &mol, double tolerance)
|
---|
| 540 | {
|
---|
| 541 | mol_ = mol;
|
---|
| 542 | tol_ = tolerance;
|
---|
| 543 | }
|
---|
| 544 |
|
---|
| 545 | void
|
---|
| 546 | IntegrationWeight::done()
|
---|
| 547 | {
|
---|
| 548 | }
|
---|
| 549 |
|
---|
| 550 | void
|
---|
| 551 | IntegrationWeight::fd_w(int icenter, SCVector3 &point,
|
---|
| 552 | double *fd_grad_w)
|
---|
| 553 | {
|
---|
| 554 | if (!fd_grad_w) return;
|
---|
| 555 | double delta = 0.001;
|
---|
| 556 | int natom = mol_->natom();
|
---|
| 557 | Ref<Molecule> molsav = mol_;
|
---|
| 558 | Ref<Molecule> dmol = new Molecule(*mol_.pointer());
|
---|
| 559 | for (int i=0; i<natom; i++) {
|
---|
| 560 | for (int j=0; j<3; j++) {
|
---|
| 561 | dmol->r(i,j) += delta;
|
---|
| 562 | if (icenter == i) point[j] += delta;
|
---|
| 563 | init(dmol,tol_);
|
---|
| 564 | double w_plus = w(icenter, point);
|
---|
| 565 | dmol->r(i,j) -= 2*delta;
|
---|
| 566 | if (icenter == i) point[j] -= 2*delta;
|
---|
| 567 | init(dmol,tol_);
|
---|
| 568 | double w_minus = w(icenter, point);
|
---|
| 569 | dmol->r(i,j) += delta;
|
---|
| 570 | if (icenter == i) point[j] += delta;
|
---|
| 571 | fd_grad_w[i*3+j] = (w_plus-w_minus)/(2.0*delta);
|
---|
| 572 | // ExEnv::outn() << scprintf("%d,%d %12.10f %12.10f %12.10f",
|
---|
| 573 | // i,j,w_plus,w_minus,fd_grad_w[i*3+j])
|
---|
| 574 | // << endl;
|
---|
| 575 | }
|
---|
| 576 | }
|
---|
| 577 | init(molsav, tol_);
|
---|
| 578 | }
|
---|
| 579 |
|
---|
| 580 | void
|
---|
| 581 | IntegrationWeight::test(int icenter, SCVector3 &point)
|
---|
| 582 | {
|
---|
| 583 | int natom = mol_->natom();
|
---|
| 584 | int natom3 = natom*3;
|
---|
| 585 |
|
---|
| 586 | // tests over sums of weights
|
---|
| 587 | int i;
|
---|
| 588 | double sum_weight = 0.0;
|
---|
| 589 | for (i=0; i<natom; i++) {
|
---|
| 590 | double weight = w(i,point);
|
---|
| 591 | sum_weight += weight;
|
---|
| 592 | }
|
---|
| 593 | if (fabs(1.0 - sum_weight) > DBL_EPSILON) {
|
---|
| 594 | ExEnv::out0() << "IntegrationWeight::test: failed on weight" << endl;
|
---|
| 595 | ExEnv::out0() << "sum_w = " << sum_weight << endl;
|
---|
| 596 | }
|
---|
| 597 |
|
---|
| 598 | // finite displacement tests of weight gradients
|
---|
| 599 | double *fd_grad_w = new double[natom3];
|
---|
| 600 | double *an_grad_w = new double[natom3];
|
---|
| 601 | w(icenter, point, an_grad_w);
|
---|
| 602 | fd_w(icenter, point, fd_grad_w);
|
---|
| 603 | for (i=0; i<natom3; i++) {
|
---|
| 604 | double mag = fabs(fd_grad_w[i]);
|
---|
| 605 | double err = fabs(fd_grad_w[i]-an_grad_w[i]);
|
---|
| 606 | int bad = 0;
|
---|
| 607 | if (mag > 0.00001 && err/mag > 0.01) bad = 1;
|
---|
| 608 | else if (err > 0.00001) bad = 1;
|
---|
| 609 | if (bad) {
|
---|
| 610 | ExEnv::out0() << "iatom = " << i/3
|
---|
| 611 | << " ixyx = " << i%3
|
---|
| 612 | << " icenter = " << icenter << " point = " << point << endl;
|
---|
| 613 | ExEnv::out0() << scprintf("dw/dx bad: fd_val=%16.13f an_val=%16.13f err=%16.13f",
|
---|
| 614 | fd_grad_w[i], an_grad_w[i],
|
---|
| 615 | fd_grad_w[i]-an_grad_w[i])
|
---|
| 616 | << endl;
|
---|
| 617 | }
|
---|
| 618 | }
|
---|
| 619 | delete[] fd_grad_w;
|
---|
| 620 | delete[] an_grad_w;
|
---|
| 621 | }
|
---|
| 622 |
|
---|
| 623 | void
|
---|
| 624 | IntegrationWeight::test()
|
---|
| 625 | {
|
---|
| 626 | SCVector3 point;
|
---|
| 627 | for (int icenter=0; icenter<mol_->natom(); icenter++) {
|
---|
| 628 | for (point[0]=-1; point[0]<=1; point[0]++) {
|
---|
| 629 | for (point[1]=-1; point[1]<=1; point[1]++) {
|
---|
| 630 | for (point[2]=-1; point[2]<=1; point[2]++) {
|
---|
| 631 | test(icenter, point);
|
---|
| 632 | }
|
---|
| 633 | }
|
---|
| 634 | }
|
---|
| 635 | }
|
---|
| 636 | }
|
---|
| 637 |
|
---|
| 638 | ///////////////////////////////////////////////////////////////////////////
|
---|
| 639 | // BeckeIntegrationWeight
|
---|
| 640 |
|
---|
| 641 | // utility functions
|
---|
| 642 |
|
---|
| 643 | inline static double
|
---|
| 644 | calc_s(double m)
|
---|
| 645 | {
|
---|
| 646 | double m1 = 1.5*m - 0.5*m*m*m;
|
---|
| 647 | double m2 = 1.5*m1 - 0.5*m1*m1*m1;
|
---|
| 648 | double m3 = 1.5*m2 - 0.5*m2*m2*m2;
|
---|
| 649 | return 0.5*(1.0-m3);
|
---|
| 650 | }
|
---|
| 651 |
|
---|
| 652 | inline static double
|
---|
| 653 | calc_f3_prime(double m)
|
---|
| 654 | {
|
---|
| 655 | double m1 = 1.5*m - 0.5*m*m*m;
|
---|
| 656 | double m2 = 1.5*m1 - 0.5*m1*m1*m1;
|
---|
| 657 | double m3 = 1.5 *(1.0 - m2*m2);
|
---|
| 658 | double n2 = 1.5 *(1.0 - m1*m1);
|
---|
| 659 | double o1 = 1.5 *(1.0 - m*m);
|
---|
| 660 | return m3*n2*o1;
|
---|
| 661 | }
|
---|
| 662 |
|
---|
| 663 | static ClassDesc BeckeIntegrationWeight_cd(
|
---|
| 664 | typeid(BeckeIntegrationWeight),"BeckeIntegrationWeight",1,"public IntegrationWeight",
|
---|
| 665 | 0, create<BeckeIntegrationWeight>, create<BeckeIntegrationWeight>);
|
---|
| 666 |
|
---|
| 667 | BeckeIntegrationWeight::BeckeIntegrationWeight(StateIn& s):
|
---|
| 668 | SavableState(s),
|
---|
| 669 | IntegrationWeight(s)
|
---|
| 670 | {
|
---|
| 671 | n_integration_centers = 0;
|
---|
| 672 | atomic_radius = 0;
|
---|
| 673 | a_mat = 0;
|
---|
| 674 | oorab = 0;
|
---|
| 675 | centers = 0;
|
---|
| 676 | }
|
---|
| 677 |
|
---|
| 678 | BeckeIntegrationWeight::BeckeIntegrationWeight()
|
---|
| 679 | {
|
---|
| 680 | n_integration_centers = 0;
|
---|
| 681 | centers = 0;
|
---|
| 682 | atomic_radius = 0;
|
---|
| 683 | a_mat = 0;
|
---|
| 684 | oorab = 0;
|
---|
| 685 | }
|
---|
| 686 |
|
---|
| 687 | BeckeIntegrationWeight::BeckeIntegrationWeight(const Ref<KeyVal>& keyval):
|
---|
| 688 | IntegrationWeight(keyval)
|
---|
| 689 | {
|
---|
| 690 | n_integration_centers = 0;
|
---|
| 691 | centers = 0;
|
---|
| 692 | atomic_radius = 0;
|
---|
| 693 | a_mat = 0;
|
---|
| 694 | oorab = 0;
|
---|
| 695 | }
|
---|
| 696 |
|
---|
| 697 | BeckeIntegrationWeight::~BeckeIntegrationWeight()
|
---|
| 698 | {
|
---|
| 699 | done();
|
---|
| 700 | }
|
---|
| 701 |
|
---|
| 702 | void
|
---|
| 703 | BeckeIntegrationWeight::save_data_state(StateOut& s)
|
---|
| 704 | {
|
---|
| 705 | IntegrationWeight::save_data_state(s);
|
---|
| 706 | }
|
---|
| 707 |
|
---|
| 708 | void
|
---|
| 709 | BeckeIntegrationWeight::init(const Ref<Molecule> &mol, double tolerance)
|
---|
| 710 | {
|
---|
| 711 | done();
|
---|
| 712 | IntegrationWeight::init(mol, tolerance);
|
---|
| 713 |
|
---|
| 714 | // We only want to include to include "atoms" that correspond to nuclei,
|
---|
| 715 | // not charges.
|
---|
| 716 | n_integration_centers = mol->n_non_q_atom();
|
---|
| 717 |
|
---|
| 718 | double *atomic_radius = new double[n_integration_centers];
|
---|
| 719 | centers = new SCVector3[n_integration_centers];
|
---|
| 720 |
|
---|
| 721 | for (int icenter=0; icenter<n_integration_centers; icenter++) {
|
---|
| 722 | int iatom = mol->non_q_atom(icenter);
|
---|
| 723 |
|
---|
| 724 | atomic_radius[icenter] = get_radius(mol, iatom);
|
---|
| 725 |
|
---|
| 726 | centers[icenter].x() = mol->r(iatom,0);
|
---|
| 727 | centers[icenter].y() = mol->r(iatom,1);
|
---|
| 728 | centers[icenter].z() = mol->r(iatom,2);
|
---|
| 729 | }
|
---|
| 730 |
|
---|
| 731 | a_mat = new double*[n_integration_centers];
|
---|
| 732 | a_mat[0] = new double[n_integration_centers*n_integration_centers];
|
---|
| 733 | oorab = new double*[n_integration_centers];
|
---|
| 734 | oorab[0] = new double[n_integration_centers*n_integration_centers];
|
---|
| 735 |
|
---|
| 736 | for (int icenter=0; icenter<n_integration_centers; icenter++) {
|
---|
| 737 | if (icenter) {
|
---|
| 738 | a_mat[icenter] = &a_mat[icenter-1][n_integration_centers];
|
---|
| 739 | oorab[icenter] = &oorab[icenter-1][n_integration_centers];
|
---|
| 740 | }
|
---|
| 741 |
|
---|
| 742 | double atomic_radius_a = atomic_radius[icenter];
|
---|
| 743 | for (int jcenter=0; jcenter < n_integration_centers; jcenter++) {
|
---|
| 744 | double chi=atomic_radius_a/atomic_radius[jcenter];
|
---|
| 745 | double uab=(chi-1.)/(chi+1.);
|
---|
| 746 | a_mat[icenter][jcenter] = uab/(uab*uab-1.);
|
---|
| 747 | if (icenter!=jcenter) {
|
---|
| 748 | oorab[icenter][jcenter]
|
---|
| 749 | = 1./centers[icenter].dist(centers[jcenter]);
|
---|
| 750 | }
|
---|
| 751 | else {
|
---|
| 752 | oorab[icenter][jcenter] = 0.0;
|
---|
| 753 | }
|
---|
| 754 | }
|
---|
| 755 | }
|
---|
| 756 |
|
---|
| 757 | }
|
---|
| 758 |
|
---|
| 759 | void
|
---|
| 760 | BeckeIntegrationWeight::done()
|
---|
| 761 | {
|
---|
| 762 | delete[] atomic_radius;
|
---|
| 763 | atomic_radius = 0;
|
---|
| 764 |
|
---|
| 765 | delete[] centers;
|
---|
| 766 | centers = 0;
|
---|
| 767 |
|
---|
| 768 | if (a_mat) {
|
---|
| 769 | delete[] a_mat[0];
|
---|
| 770 | delete[] a_mat;
|
---|
| 771 | a_mat = 0;
|
---|
| 772 | }
|
---|
| 773 |
|
---|
| 774 | if (oorab) {
|
---|
| 775 | delete[] oorab[0];
|
---|
| 776 | delete[] oorab;
|
---|
| 777 | oorab = 0;
|
---|
| 778 | }
|
---|
| 779 |
|
---|
| 780 | n_integration_centers = 0;
|
---|
| 781 | }
|
---|
| 782 |
|
---|
| 783 | double
|
---|
| 784 | BeckeIntegrationWeight::compute_p(int icenter, SCVector3&point)
|
---|
| 785 | {
|
---|
| 786 | double ra = point.dist(centers[icenter]);
|
---|
| 787 | double *ooraba = oorab[icenter];
|
---|
| 788 | double *aa = a_mat[icenter];
|
---|
| 789 |
|
---|
| 790 | double p = 1.0;
|
---|
| 791 | for (int jcenter=0; jcenter < n_integration_centers; jcenter++) {
|
---|
| 792 | if (icenter != jcenter) {
|
---|
| 793 | double mu = (ra-point.dist(centers[jcenter]))*ooraba[jcenter];
|
---|
| 794 |
|
---|
| 795 | if (mu <= -1.)
|
---|
| 796 | continue; // s(-1) == 1.0
|
---|
| 797 | else if (mu >= 1.) {
|
---|
| 798 | return 0.0; // s(1) == 0.0
|
---|
| 799 | }
|
---|
| 800 | else
|
---|
| 801 | p *= calc_s(mu + aa[jcenter]*(1.-mu*mu));
|
---|
| 802 | }
|
---|
| 803 | }
|
---|
| 804 |
|
---|
| 805 | return p;
|
---|
| 806 | }
|
---|
| 807 |
|
---|
| 808 | // compute derivative of mu(grad_center,bcenter) wrt grad_center;
|
---|
| 809 | // NB: the derivative is independent of the (implicit) wcenter
|
---|
| 810 | // provided that wcenter!=grad_center
|
---|
| 811 | void
|
---|
| 812 | BeckeIntegrationWeight::compute_grad_nu(int grad_center, int bcenter,
|
---|
| 813 | SCVector3 &point, SCVector3 &grad)
|
---|
| 814 | {
|
---|
| 815 | SCVector3 r_g = point - centers[grad_center];
|
---|
| 816 | SCVector3 r_b = point - centers[bcenter];
|
---|
| 817 | SCVector3 r_gb = centers[grad_center] - centers[bcenter];
|
---|
| 818 | double mag_r_g = r_g.norm();
|
---|
| 819 | double mag_r_b = r_b.norm();
|
---|
| 820 | double oorgb = oorab[grad_center][bcenter];
|
---|
| 821 | double mu = (mag_r_g-mag_r_b)*oorgb;
|
---|
| 822 | double a_gb = a_mat[grad_center][bcenter];
|
---|
| 823 | double coef = 1.0-2.0*a_gb*mu;
|
---|
| 824 | double r_g_coef;
|
---|
| 825 |
|
---|
| 826 | if (mag_r_g < 10.0 * DBL_EPSILON) r_g_coef = 0.0;
|
---|
| 827 | else r_g_coef = -coef*oorgb/mag_r_g;
|
---|
| 828 | int ixyz;
|
---|
| 829 | for (ixyz=0; ixyz<3; ixyz++) grad[ixyz] = r_g_coef * r_g[ixyz];
|
---|
| 830 | double r_gb_coef = coef*(mag_r_b - mag_r_g)*oorgb*oorgb*oorgb;
|
---|
| 831 | for (ixyz=0; ixyz<3; ixyz++) grad[ixyz] += r_gb_coef * r_gb[ixyz];
|
---|
| 832 | }
|
---|
| 833 |
|
---|
| 834 | // compute t(nu_ij)
|
---|
| 835 | double
|
---|
| 836 | BeckeIntegrationWeight::compute_t(int icenter, int jcenter, SCVector3 &point)
|
---|
| 837 | {
|
---|
| 838 | // Cf. Johnson et al., JCP v. 98, p. 5612 (1993) (Appendix B)
|
---|
| 839 | // NB: t is zero if s is zero
|
---|
| 840 |
|
---|
| 841 | SCVector3 r_i = point - centers[icenter];
|
---|
| 842 | SCVector3 r_j = point - centers[jcenter];
|
---|
| 843 | SCVector3 r_ij = centers[icenter] - centers[jcenter];
|
---|
| 844 | double t;
|
---|
| 845 | double mag_r_j = r_j.norm();
|
---|
| 846 | double mag_r_i = r_i.norm();
|
---|
| 847 | double mu = (mag_r_i-mag_r_j)*oorab[icenter][jcenter];
|
---|
| 848 | if (mu >= 1.0-100*DBL_EPSILON) {
|
---|
| 849 | t = 0.0;
|
---|
| 850 | return t;
|
---|
| 851 | }
|
---|
| 852 |
|
---|
| 853 | double a_ij = a_mat[icenter][jcenter];
|
---|
| 854 | double nu = mu + a_ij*(1.-mu*mu);
|
---|
| 855 | double s;
|
---|
| 856 | if (mu <= -1.0) s = 1.0;
|
---|
| 857 | else s = calc_s(nu);
|
---|
| 858 | if (fabs(s) < 10*DBL_EPSILON) {
|
---|
| 859 | t = 0.0;
|
---|
| 860 | return t;
|
---|
| 861 | }
|
---|
| 862 | double p1 = 1.5*nu - 0.5*nu*nu*nu;
|
---|
| 863 | double p2 = 1.5*p1 - 0.5*p1*p1*p1;
|
---|
| 864 |
|
---|
| 865 | t = -(27.0/16.0) * (1 - p2*p2) * (1 - p1*p1) * (1 - nu*nu) / s;
|
---|
| 866 | return t;
|
---|
| 867 | }
|
---|
| 868 |
|
---|
| 869 | void
|
---|
| 870 | BeckeIntegrationWeight::compute_grad_p(int grad_center, int bcenter,
|
---|
| 871 | int wcenter, SCVector3&point,
|
---|
| 872 | double p, SCVector3&grad)
|
---|
| 873 | {
|
---|
| 874 | // the gradient of p is computed using the formulae from
|
---|
| 875 | // Johnson et al., JCP v. 98, p. 5612 (1993) (Appendix B)
|
---|
| 876 |
|
---|
| 877 | if (grad_center == bcenter) {
|
---|
| 878 | grad = 0.0;
|
---|
| 879 | for (int dcenter=0; dcenter<n_integration_centers; dcenter++) {
|
---|
| 880 | if (dcenter == bcenter) continue;
|
---|
| 881 | SCVector3 grad_nu;
|
---|
| 882 | compute_grad_nu(grad_center, dcenter, point, grad_nu);
|
---|
| 883 | double t = compute_t(grad_center,dcenter,point);
|
---|
| 884 | for (int ixyz=0; ixyz<3; ixyz++) grad[ixyz] += t * grad_nu[ixyz];
|
---|
| 885 | }
|
---|
| 886 | }
|
---|
| 887 | else {
|
---|
| 888 | SCVector3 grad_nu;
|
---|
| 889 | compute_grad_nu(grad_center, bcenter, point, grad_nu);
|
---|
| 890 | double t = compute_t(bcenter,grad_center,point);
|
---|
| 891 | for (int ixyz=0; ixyz<3; ixyz++) grad[ixyz] = -t * grad_nu[ixyz];
|
---|
| 892 | }
|
---|
| 893 | grad *= p;
|
---|
| 894 | }
|
---|
| 895 |
|
---|
| 896 | double
|
---|
| 897 | BeckeIntegrationWeight::w(int acenter, SCVector3 &point,
|
---|
| 898 | double *w_gradient)
|
---|
| 899 | {
|
---|
| 900 | int icenter;
|
---|
| 901 | double p_sum=0.0, p_a=0.0;
|
---|
| 902 |
|
---|
| 903 | for (icenter=0; icenter<n_integration_centers; icenter++) {
|
---|
| 904 | double p_tmp = compute_p(icenter, point);
|
---|
| 905 | if (icenter==acenter) p_a=p_tmp;
|
---|
| 906 | p_sum += p_tmp;
|
---|
| 907 | }
|
---|
| 908 | double w_a = p_a/p_sum;
|
---|
| 909 |
|
---|
| 910 | if (w_gradient) {
|
---|
| 911 | // w_gradient is computed using the formulae from
|
---|
| 912 | // Johnson et al., JCP v. 98, p. 5612 (1993) (Appendix B)
|
---|
| 913 | int i,j;
|
---|
| 914 | for (i=0; i<n_integration_centers*3; i++ ) w_gradient[i] = 0.0;
|
---|
| 915 | // fd_w(acenter, point, w_gradient); // imbn commented out for debug
|
---|
| 916 | // ExEnv::outn() << point << " ";
|
---|
| 917 | // for (int i=0; i<n_integration_centers*3; i++) {
|
---|
| 918 | // ExEnv::outn() << scprintf(" %10.6f", w_gradient[i]);
|
---|
| 919 | // }
|
---|
| 920 | // ExEnv::outn() << endl;
|
---|
| 921 | // return w_a; // imbn commented out for debug
|
---|
| 922 | for (int ccenter = 0; ccenter < n_integration_centers; ccenter++) {
|
---|
| 923 | // NB: for ccenter==acenter, use translational invariance
|
---|
| 924 | // to get the corresponding component of the gradient
|
---|
| 925 | if (ccenter != acenter) {
|
---|
| 926 | SCVector3 grad_c_w_a;
|
---|
| 927 | SCVector3 grad_c_p_a;
|
---|
| 928 | compute_grad_p(ccenter, acenter, acenter, point, p_a, grad_c_p_a);
|
---|
| 929 | for (i=0; i<3; i++) grad_c_w_a[i] = grad_c_p_a[i]/p_sum;
|
---|
| 930 | for (int bcenter=0; bcenter<n_integration_centers; bcenter++) {
|
---|
| 931 | SCVector3 grad_c_p_b;
|
---|
| 932 | double p_b = compute_p(bcenter,point);
|
---|
| 933 | compute_grad_p(ccenter, bcenter, acenter, point, p_b,
|
---|
| 934 | grad_c_p_b);
|
---|
| 935 | for (i=0; i<3; i++) grad_c_w_a[i] -= w_a*grad_c_p_b[i]/p_sum;
|
---|
| 936 | }
|
---|
| 937 | for (i=0; i<3; i++) w_gradient[ccenter*3+i] = grad_c_w_a[i];
|
---|
| 938 | }
|
---|
| 939 | }
|
---|
| 940 | // fill in w_gradient for ccenter==acenter
|
---|
| 941 | for (j=0; j<3; j++) {
|
---|
| 942 | for (i=0; i<n_integration_centers; i++) {
|
---|
| 943 | if (i != acenter) {
|
---|
| 944 | w_gradient[acenter*3+j] -= w_gradient[i*3+j];
|
---|
| 945 | }
|
---|
| 946 | }
|
---|
| 947 | }
|
---|
| 948 | }
|
---|
| 949 |
|
---|
| 950 | return w_a;
|
---|
| 951 | }
|
---|
| 952 |
|
---|
| 953 | ///////////////////////////////////////////////////
|
---|
| 954 | // RadialIntegrator
|
---|
| 955 | static ClassDesc RadialIntegrator_cd(
|
---|
| 956 | typeid(RadialIntegrator),"RadialIntegrator",1,"public SavableState",
|
---|
| 957 | 0, 0, 0);
|
---|
| 958 |
|
---|
| 959 | RadialIntegrator::RadialIntegrator(StateIn& s):
|
---|
| 960 | SavableState(s)
|
---|
| 961 | {
|
---|
| 962 | }
|
---|
| 963 |
|
---|
| 964 | RadialIntegrator::RadialIntegrator()
|
---|
| 965 | {
|
---|
| 966 | }
|
---|
| 967 |
|
---|
| 968 | RadialIntegrator::RadialIntegrator(const Ref<KeyVal>& keyval)
|
---|
| 969 | {
|
---|
| 970 | }
|
---|
| 971 |
|
---|
| 972 | RadialIntegrator::~RadialIntegrator()
|
---|
| 973 | {
|
---|
| 974 | }
|
---|
| 975 |
|
---|
| 976 | void
|
---|
| 977 | RadialIntegrator::save_data_state(StateOut& s)
|
---|
| 978 | {
|
---|
| 979 | }
|
---|
| 980 |
|
---|
| 981 | ///////////////////////////////////////
|
---|
| 982 | // AngularIntegrator
|
---|
| 983 |
|
---|
| 984 | static ClassDesc AngularIntegrator_cd(
|
---|
| 985 | typeid(AngularIntegrator),"AngularIntegrator",1,"public SavableState",
|
---|
| 986 | 0, 0, 0);
|
---|
| 987 |
|
---|
| 988 | AngularIntegrator::AngularIntegrator(StateIn& s):
|
---|
| 989 | SavableState(s)
|
---|
| 990 | {
|
---|
| 991 | }
|
---|
| 992 |
|
---|
| 993 | AngularIntegrator::AngularIntegrator()
|
---|
| 994 | {
|
---|
| 995 | }
|
---|
| 996 |
|
---|
| 997 | AngularIntegrator::AngularIntegrator(const Ref<KeyVal>& keyval)
|
---|
| 998 | {
|
---|
| 999 | }
|
---|
| 1000 |
|
---|
| 1001 | AngularIntegrator::~AngularIntegrator()
|
---|
| 1002 | {
|
---|
| 1003 | }
|
---|
| 1004 |
|
---|
| 1005 | void
|
---|
| 1006 | AngularIntegrator::save_data_state(StateOut& s)
|
---|
| 1007 | {
|
---|
| 1008 | }
|
---|
| 1009 |
|
---|
| 1010 | ///////////////////////////////////////
|
---|
| 1011 | // EulerMaclaurinRadialIntegrator
|
---|
| 1012 |
|
---|
| 1013 | static ClassDesc EulerMaclaurinRadialIntegrator_cd(
|
---|
| 1014 | typeid(EulerMaclaurinRadialIntegrator),"EulerMaclaurinRadialIntegrator",1,"public RadialIntegrator",
|
---|
| 1015 | 0, create<EulerMaclaurinRadialIntegrator>, create<EulerMaclaurinRadialIntegrator>);
|
---|
| 1016 |
|
---|
| 1017 | EulerMaclaurinRadialIntegrator::EulerMaclaurinRadialIntegrator(StateIn& s):
|
---|
| 1018 | SavableState(s),
|
---|
| 1019 | RadialIntegrator(s)
|
---|
| 1020 | {
|
---|
| 1021 | s.get(nr_);
|
---|
| 1022 | }
|
---|
| 1023 |
|
---|
| 1024 | EulerMaclaurinRadialIntegrator::EulerMaclaurinRadialIntegrator()
|
---|
| 1025 | {
|
---|
| 1026 | nr_ = 75;
|
---|
| 1027 | }
|
---|
| 1028 |
|
---|
| 1029 | EulerMaclaurinRadialIntegrator::EulerMaclaurinRadialIntegrator(int nr_points)
|
---|
| 1030 | {
|
---|
| 1031 | nr_ = nr_points;
|
---|
| 1032 | }
|
---|
| 1033 |
|
---|
| 1034 | EulerMaclaurinRadialIntegrator::EulerMaclaurinRadialIntegrator(const Ref<KeyVal>& keyval):
|
---|
| 1035 | RadialIntegrator(keyval)
|
---|
| 1036 | {
|
---|
| 1037 | nr_ = keyval->intvalue("nr", KeyValValueint(75));
|
---|
| 1038 | }
|
---|
| 1039 |
|
---|
| 1040 | EulerMaclaurinRadialIntegrator::~EulerMaclaurinRadialIntegrator()
|
---|
| 1041 | {
|
---|
| 1042 | }
|
---|
| 1043 |
|
---|
| 1044 | void
|
---|
| 1045 | EulerMaclaurinRadialIntegrator::save_data_state(StateOut& s)
|
---|
| 1046 | {
|
---|
| 1047 | RadialIntegrator::save_data_state(s);
|
---|
| 1048 | s.put(nr_);
|
---|
| 1049 | }
|
---|
| 1050 |
|
---|
| 1051 | int
|
---|
| 1052 | EulerMaclaurinRadialIntegrator::nr() const
|
---|
| 1053 | {
|
---|
| 1054 | return nr_;
|
---|
| 1055 | }
|
---|
| 1056 |
|
---|
| 1057 | double
|
---|
| 1058 | EulerMaclaurinRadialIntegrator::radial_value(int ir, int nr, double radii,
|
---|
| 1059 | double &multiplier)
|
---|
| 1060 | {
|
---|
| 1061 | double q = (double)ir/(double)nr;
|
---|
| 1062 | double value = q/(1.-q);
|
---|
| 1063 | double r = radii*value*value;
|
---|
| 1064 | double dr_dq = 2.*radii*q*pow(1.-q,-3.);
|
---|
| 1065 | double dr_dqr2 = dr_dq*r*r;
|
---|
| 1066 | multiplier = dr_dqr2/nr;
|
---|
| 1067 | return r;
|
---|
| 1068 | }
|
---|
| 1069 |
|
---|
| 1070 | void
|
---|
| 1071 | EulerMaclaurinRadialIntegrator::print(ostream &o) const
|
---|
| 1072 | {
|
---|
| 1073 | o << indent
|
---|
| 1074 | << scprintf("%s: nr = %d", class_name(), nr()) << endl;
|
---|
| 1075 | }
|
---|
| 1076 |
|
---|
| 1077 | //////////////////////////////////////////////////////////////////////////
|
---|
| 1078 | // LebedevLaikovIntegrator
|
---|
| 1079 |
|
---|
| 1080 | static ClassDesc LebedevLaikovIntegrator_cd(
|
---|
| 1081 | typeid(LebedevLaikovIntegrator),"LebedevLaikovIntegrator",1,"public AngularIntegrator",
|
---|
| 1082 | 0, create<LebedevLaikovIntegrator>, create<LebedevLaikovIntegrator>);
|
---|
| 1083 |
|
---|
| 1084 | LebedevLaikovIntegrator::LebedevLaikovIntegrator(StateIn& s):
|
---|
| 1085 | SavableState(s),
|
---|
| 1086 | AngularIntegrator(s)
|
---|
| 1087 | {
|
---|
| 1088 | s.get(npoint_);
|
---|
| 1089 | init(npoint_);
|
---|
| 1090 | }
|
---|
| 1091 |
|
---|
| 1092 | LebedevLaikovIntegrator::LebedevLaikovIntegrator()
|
---|
| 1093 | {
|
---|
| 1094 | init(302);
|
---|
| 1095 | }
|
---|
| 1096 |
|
---|
| 1097 | LebedevLaikovIntegrator::LebedevLaikovIntegrator(int npoint)
|
---|
| 1098 | {
|
---|
| 1099 | init(npoint);
|
---|
| 1100 | }
|
---|
| 1101 |
|
---|
| 1102 | LebedevLaikovIntegrator::LebedevLaikovIntegrator(const Ref<KeyVal>& keyval)
|
---|
| 1103 | {
|
---|
| 1104 | KeyValValueint defnpoint(302);
|
---|
| 1105 |
|
---|
| 1106 | init(keyval->intvalue("n", defnpoint));
|
---|
| 1107 | }
|
---|
| 1108 |
|
---|
| 1109 | LebedevLaikovIntegrator::~LebedevLaikovIntegrator()
|
---|
| 1110 | {
|
---|
| 1111 | delete [] x_;
|
---|
| 1112 | delete [] y_;
|
---|
| 1113 | delete [] z_;
|
---|
| 1114 | delete [] w_;
|
---|
| 1115 | }
|
---|
| 1116 |
|
---|
| 1117 | void
|
---|
| 1118 | LebedevLaikovIntegrator::save_data_state(StateOut& s)
|
---|
| 1119 | {
|
---|
| 1120 | AngularIntegrator::save_data_state(s);
|
---|
| 1121 | s.put(npoint_);
|
---|
| 1122 | }
|
---|
| 1123 |
|
---|
| 1124 | extern "C" {
|
---|
| 1125 | int Lebedev_Laikov_sphere (int N, double *X, double *Y, double *Z,
|
---|
| 1126 | double *W);
|
---|
| 1127 | int Lebedev_Laikov_npoint (int lvalue);
|
---|
| 1128 | }
|
---|
| 1129 |
|
---|
| 1130 | int
|
---|
| 1131 | LebedevLaikovIntegrator::nw(void) const
|
---|
| 1132 | {
|
---|
| 1133 | return npoint_;
|
---|
| 1134 | }
|
---|
| 1135 |
|
---|
| 1136 | void
|
---|
| 1137 | LebedevLaikovIntegrator::init(int n)
|
---|
| 1138 | {
|
---|
| 1139 | // ExEnv::outn() << " LebedevLaikovIntegrator::init -> before x_, y_, z_, and w_ malloc's " << endl;
|
---|
| 1140 | // ExEnv::outn() << " n = " << n << endl;
|
---|
| 1141 |
|
---|
| 1142 | x_ = new double[n];
|
---|
| 1143 | y_ = new double[n];
|
---|
| 1144 | z_ = new double[n];
|
---|
| 1145 | w_ = new double[n];
|
---|
| 1146 |
|
---|
| 1147 | // ExEnv::outn() << " LebedevLaikovIntegrator::init -> nw_points = " << n << endl;
|
---|
| 1148 |
|
---|
| 1149 | npoint_ = Lebedev_Laikov_sphere(n, x_, y_, z_, w_);
|
---|
| 1150 | if (npoint_ != n) {
|
---|
| 1151 | ExEnv::outn() << class_name() << ": bad number of points given: " << n << endl;
|
---|
| 1152 | abort();
|
---|
| 1153 | }
|
---|
| 1154 | }
|
---|
| 1155 |
|
---|
| 1156 | int
|
---|
| 1157 | LebedevLaikovIntegrator::num_angular_points(double r_value, int ir)
|
---|
| 1158 | {
|
---|
| 1159 | if (ir == 0) return 1;
|
---|
| 1160 | return npoint_;
|
---|
| 1161 | }
|
---|
| 1162 |
|
---|
| 1163 | double
|
---|
| 1164 | LebedevLaikovIntegrator
|
---|
| 1165 | ::angular_point_cartesian(int iangular, double r,
|
---|
| 1166 | SCVector3 &integration_point) const
|
---|
| 1167 | {
|
---|
| 1168 | integration_point.x() = r*x_[iangular];
|
---|
| 1169 | integration_point.y() = r*y_[iangular];
|
---|
| 1170 | integration_point.z() = r*z_[iangular];
|
---|
| 1171 |
|
---|
| 1172 | return 4.0*M_PI*w_[iangular];
|
---|
| 1173 | }
|
---|
| 1174 |
|
---|
| 1175 | void
|
---|
| 1176 | LebedevLaikovIntegrator::print(ostream &o) const
|
---|
| 1177 | {
|
---|
| 1178 | o << indent
|
---|
| 1179 | << scprintf("%s: n = %d", class_name(), npoint_) << endl;
|
---|
| 1180 | }
|
---|
| 1181 |
|
---|
| 1182 | /////////////////////////////////
|
---|
| 1183 | // GaussLegendreAngularIntegrator
|
---|
| 1184 |
|
---|
| 1185 | static ClassDesc GaussLegendreAngularIntegrator_cd(
|
---|
| 1186 | typeid(GaussLegendreAngularIntegrator),"GaussLegendreAngularIntegrator",1,"public AngularIntegrator",
|
---|
| 1187 | 0, create<GaussLegendreAngularIntegrator>, create<GaussLegendreAngularIntegrator>);
|
---|
| 1188 |
|
---|
| 1189 | GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator(StateIn& s):
|
---|
| 1190 | SavableState(s),
|
---|
| 1191 | AngularIntegrator(s)
|
---|
| 1192 | {
|
---|
| 1193 | s.get(ntheta_);
|
---|
| 1194 | s.get(nphi_);
|
---|
| 1195 | s.get(Ktheta_);
|
---|
| 1196 | theta_quad_weights_ = new double[ntheta_];
|
---|
| 1197 | theta_quad_points_ = new double[ntheta_];
|
---|
| 1198 | }
|
---|
| 1199 |
|
---|
| 1200 | GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator()
|
---|
| 1201 | {
|
---|
| 1202 | set_ntheta(16);
|
---|
| 1203 | set_nphi(32);
|
---|
| 1204 | set_Ktheta(5);
|
---|
| 1205 | int ntheta = get_ntheta();
|
---|
| 1206 | theta_quad_weights_ = new double [ntheta];
|
---|
| 1207 | theta_quad_points_ = new double [ntheta];
|
---|
| 1208 | }
|
---|
| 1209 |
|
---|
| 1210 | GaussLegendreAngularIntegrator::GaussLegendreAngularIntegrator(const Ref<KeyVal>& keyval)
|
---|
| 1211 | {
|
---|
| 1212 | set_ntheta( keyval->intvalue("ntheta") );
|
---|
| 1213 | if (keyval->error() != KeyVal::OK) set_ntheta(16);
|
---|
| 1214 | set_nphi( keyval->intvalue("nphi") );
|
---|
| 1215 | if (keyval->error() != KeyVal::OK) set_nphi(2*get_ntheta());
|
---|
| 1216 | set_Ktheta( keyval->intvalue("Ktheta") );
|
---|
| 1217 | if (keyval->error() != KeyVal::OK) set_Ktheta(5);
|
---|
| 1218 |
|
---|
| 1219 | int ntheta = get_ntheta();
|
---|
| 1220 | theta_quad_weights_ = new double [ntheta];
|
---|
| 1221 | theta_quad_points_ = new double [ntheta];
|
---|
| 1222 | }
|
---|
| 1223 |
|
---|
| 1224 | GaussLegendreAngularIntegrator::~GaussLegendreAngularIntegrator()
|
---|
| 1225 | {
|
---|
| 1226 | delete [] theta_quad_points_;
|
---|
| 1227 | delete [] theta_quad_weights_;
|
---|
| 1228 | }
|
---|
| 1229 |
|
---|
| 1230 | void
|
---|
| 1231 | GaussLegendreAngularIntegrator::save_data_state(StateOut& s)
|
---|
| 1232 | {
|
---|
| 1233 | AngularIntegrator::save_data_state(s);
|
---|
| 1234 | s.put(ntheta_);
|
---|
| 1235 | s.put(nphi_);
|
---|
| 1236 | s.put(Ktheta_);
|
---|
| 1237 | }
|
---|
| 1238 |
|
---|
| 1239 | int
|
---|
| 1240 | GaussLegendreAngularIntegrator::get_ntheta(void) const
|
---|
| 1241 | {
|
---|
| 1242 | return ntheta_;
|
---|
| 1243 | }
|
---|
| 1244 |
|
---|
| 1245 | void
|
---|
| 1246 | GaussLegendreAngularIntegrator::set_ntheta(int i)
|
---|
| 1247 | {
|
---|
| 1248 | ntheta_ = i;
|
---|
| 1249 | }
|
---|
| 1250 |
|
---|
| 1251 | int
|
---|
| 1252 | GaussLegendreAngularIntegrator::get_nphi(void) const
|
---|
| 1253 | {
|
---|
| 1254 | return nphi_;
|
---|
| 1255 | }
|
---|
| 1256 |
|
---|
| 1257 | void
|
---|
| 1258 | GaussLegendreAngularIntegrator::set_nphi(int i)
|
---|
| 1259 | {
|
---|
| 1260 | nphi_ = i;
|
---|
| 1261 | }
|
---|
| 1262 |
|
---|
| 1263 | int
|
---|
| 1264 | GaussLegendreAngularIntegrator::get_Ktheta(void) const
|
---|
| 1265 | {
|
---|
| 1266 | return Ktheta_;
|
---|
| 1267 | }
|
---|
| 1268 |
|
---|
| 1269 | void
|
---|
| 1270 | GaussLegendreAngularIntegrator::set_Ktheta(int i)
|
---|
| 1271 | {
|
---|
| 1272 | Ktheta_ = i;
|
---|
| 1273 | }
|
---|
| 1274 |
|
---|
| 1275 | int
|
---|
| 1276 | GaussLegendreAngularIntegrator::get_ntheta_r(void) const
|
---|
| 1277 | {
|
---|
| 1278 | return ntheta_r_;
|
---|
| 1279 | }
|
---|
| 1280 |
|
---|
| 1281 | void
|
---|
| 1282 | GaussLegendreAngularIntegrator::set_ntheta_r(int i)
|
---|
| 1283 | {
|
---|
| 1284 | ntheta_r_ = i;
|
---|
| 1285 | }
|
---|
| 1286 |
|
---|
| 1287 | int
|
---|
| 1288 | GaussLegendreAngularIntegrator::get_nphi_r(void) const
|
---|
| 1289 | {
|
---|
| 1290 | return nphi_r_;
|
---|
| 1291 | }
|
---|
| 1292 |
|
---|
| 1293 | void
|
---|
| 1294 | GaussLegendreAngularIntegrator::set_nphi_r(int i)
|
---|
| 1295 | {
|
---|
| 1296 | nphi_r_ = i;
|
---|
| 1297 | }
|
---|
| 1298 |
|
---|
| 1299 | int
|
---|
| 1300 | GaussLegendreAngularIntegrator::get_Ktheta_r(void) const
|
---|
| 1301 | {
|
---|
| 1302 | return Ktheta_r_;
|
---|
| 1303 | }
|
---|
| 1304 |
|
---|
| 1305 | void
|
---|
| 1306 | GaussLegendreAngularIntegrator::set_Ktheta_r(int i)
|
---|
| 1307 | {
|
---|
| 1308 | Ktheta_r_ = i;
|
---|
| 1309 | }
|
---|
| 1310 |
|
---|
| 1311 | int
|
---|
| 1312 | GaussLegendreAngularIntegrator::nw(void) const
|
---|
| 1313 | {
|
---|
| 1314 | return nphi_*ntheta_;
|
---|
| 1315 | }
|
---|
| 1316 |
|
---|
| 1317 | double
|
---|
| 1318 | GaussLegendreAngularIntegrator::sin_theta(SCVector3 &point) const
|
---|
| 1319 | {
|
---|
| 1320 | return sin(point.theta());
|
---|
| 1321 | }
|
---|
| 1322 |
|
---|
| 1323 | int
|
---|
| 1324 | GaussLegendreAngularIntegrator::num_angular_points(double r_value,
|
---|
| 1325 | int ir)
|
---|
| 1326 | {
|
---|
| 1327 | int Ktheta, ntheta, ntheta_r;
|
---|
| 1328 |
|
---|
| 1329 | if (ir == 0) {
|
---|
| 1330 | set_ntheta_r(1);
|
---|
| 1331 | set_nphi_r(1);
|
---|
| 1332 | }
|
---|
| 1333 | else {
|
---|
| 1334 | Ktheta = get_Ktheta();
|
---|
| 1335 | ntheta = get_ntheta();
|
---|
| 1336 | ntheta_r= (int) (r_value*Ktheta*ntheta);
|
---|
| 1337 | set_ntheta_r(ntheta_r);
|
---|
| 1338 | if (ntheta_r > ntheta) set_ntheta_r(ntheta);
|
---|
| 1339 | if (ntheta_r < 6) set_ntheta_r(6);
|
---|
| 1340 | set_nphi_r(2*get_ntheta_r());
|
---|
| 1341 | }
|
---|
| 1342 |
|
---|
| 1343 | gauleg(0.0, M_PI, get_ntheta_r());
|
---|
| 1344 |
|
---|
| 1345 | return get_ntheta_r()*get_nphi_r();
|
---|
| 1346 | }
|
---|
| 1347 |
|
---|
| 1348 | void
|
---|
| 1349 | GaussLegendreAngularIntegrator::gauleg(double x1, double x2, int n)
|
---|
| 1350 | {
|
---|
| 1351 | int m,j,i;
|
---|
| 1352 | double z1,z,xm,xl,pp,p3,p2,p1;
|
---|
| 1353 | const double EPS = 10.0 * DBL_EPSILON;
|
---|
| 1354 |
|
---|
| 1355 | m=(n+1)/2;
|
---|
| 1356 | xm=0.5*(x2+x1);
|
---|
| 1357 | xl=0.5*(x2-x1);
|
---|
| 1358 | for (i=1;i<=m;i++) {
|
---|
| 1359 | z=cos(M_PI*(i-0.25)/(n+0.5));
|
---|
| 1360 | do {
|
---|
| 1361 | p1=1.0;
|
---|
| 1362 | p2=0.0;
|
---|
| 1363 | for (j=1;j<=n;j++) {
|
---|
| 1364 | p3=p2;
|
---|
| 1365 | p2=p1;
|
---|
| 1366 | p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
|
---|
| 1367 | }
|
---|
| 1368 | pp=n*(z*p1-p2)/(z*z-1.0);
|
---|
| 1369 | z1=z;
|
---|
| 1370 | z=z1-p1/pp;
|
---|
| 1371 | } while (fabs(z-z1) > EPS);
|
---|
| 1372 | theta_quad_points_[i-1]=xm-xl*z;
|
---|
| 1373 | theta_quad_points_[n-i]=xm+xl*z;
|
---|
| 1374 | theta_quad_weights_[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
|
---|
| 1375 | theta_quad_weights_[n-i]=theta_quad_weights_[i-1];
|
---|
| 1376 | }
|
---|
| 1377 | }
|
---|
| 1378 |
|
---|
| 1379 | double
|
---|
| 1380 | GaussLegendreAngularIntegrator
|
---|
| 1381 | ::angular_point_cartesian(int iangular, double r,
|
---|
| 1382 | SCVector3 &integration_point) const
|
---|
| 1383 | {
|
---|
| 1384 | int itheta, iphi, nphi_r;
|
---|
| 1385 |
|
---|
| 1386 | nphi_r = get_nphi_r();
|
---|
| 1387 | itheta = iangular/nphi_r;
|
---|
| 1388 | iphi = iangular - itheta*nphi_r;
|
---|
| 1389 | SCVector3 point;
|
---|
| 1390 | point.theta() = theta_quad_points_[itheta];
|
---|
| 1391 | point.phi() = (double) iphi/ (double) nphi_r * 2.0 * M_PI;
|
---|
| 1392 | point.r() = r;
|
---|
| 1393 | point.spherical_to_cartesian(integration_point);
|
---|
| 1394 | return ( sin_theta(point)*theta_quad_weights_[itheta]*2.0*M_PI/(double)nphi_r );
|
---|
| 1395 | }
|
---|
| 1396 |
|
---|
| 1397 | void
|
---|
| 1398 | GaussLegendreAngularIntegrator::print(ostream &o) const
|
---|
| 1399 | {
|
---|
| 1400 | o << indent << class_name() << ":" << endl;
|
---|
| 1401 | o << incindent;
|
---|
| 1402 | o << indent << scprintf("ntheta = %5d", get_ntheta()) << endl;
|
---|
| 1403 | o << indent << scprintf("nphi = %5d", get_nphi()) << endl;
|
---|
| 1404 | o << indent << scprintf("Ktheta = %5d", get_Ktheta()) << endl;
|
---|
| 1405 | o << decindent;
|
---|
| 1406 | }
|
---|
| 1407 |
|
---|
| 1408 | //////////////////////////////////////////////
|
---|
| 1409 | // RadialAngularIntegratorThread
|
---|
| 1410 |
|
---|
| 1411 | class RadialAngularIntegratorThread: public DenIntegratorThread {
|
---|
| 1412 | protected:
|
---|
| 1413 | SCVector3 *centers_;
|
---|
| 1414 | int *nr_;
|
---|
| 1415 | double *atomic_radius_;
|
---|
| 1416 | Molecule *mol_;
|
---|
| 1417 | RadialAngularIntegrator *ra_integrator_;
|
---|
| 1418 | IntegrationWeight *weight_;
|
---|
| 1419 | int point_count_total_;
|
---|
| 1420 | double total_density_;
|
---|
| 1421 | public:
|
---|
| 1422 | RadialAngularIntegratorThread(int ithread, int nthread,
|
---|
| 1423 | RadialAngularIntegrator *integrator,
|
---|
| 1424 | DenFunctional *func,
|
---|
| 1425 | const Ref<BatchElectronDensity> &den,
|
---|
| 1426 | int linear_scaling, int use_dmat_bound,
|
---|
| 1427 | double accuracy,
|
---|
| 1428 | int compute_potential_integrals,
|
---|
| 1429 | int need_nuclear_gradient);
|
---|
| 1430 | ~RadialAngularIntegratorThread();
|
---|
| 1431 | void run();
|
---|
| 1432 | double total_density() { return total_density_; }
|
---|
| 1433 | int point_count() { return point_count_total_; }
|
---|
| 1434 | };
|
---|
| 1435 |
|
---|
| 1436 | RadialAngularIntegratorThread
|
---|
| 1437 | ::RadialAngularIntegratorThread(int ithread, int nthread,
|
---|
| 1438 | RadialAngularIntegrator *integrator,
|
---|
| 1439 | DenFunctional *func,
|
---|
| 1440 | const Ref<BatchElectronDensity> &den,
|
---|
| 1441 | int linear_scaling, int use_dmat_bound,
|
---|
| 1442 | double accuracy,
|
---|
| 1443 | int compute_potential_integrals,
|
---|
| 1444 | int need_nuclear_gradient):
|
---|
| 1445 | DenIntegratorThread(ithread,nthread,
|
---|
| 1446 | integrator, func,
|
---|
| 1447 | den,
|
---|
| 1448 | linear_scaling, use_dmat_bound,
|
---|
| 1449 | accuracy,
|
---|
| 1450 | compute_potential_integrals,
|
---|
| 1451 | need_nuclear_gradient)
|
---|
| 1452 | {
|
---|
| 1453 | int icenter;
|
---|
| 1454 | ra_integrator_ = integrator;
|
---|
| 1455 | int deriv_order = (need_nuclear_gradient==0?0:1);
|
---|
| 1456 |
|
---|
| 1457 | mol_ = integrator_->wavefunction()->molecule().pointer();
|
---|
| 1458 |
|
---|
| 1459 | weight_ = ra_integrator_->weight().pointer();
|
---|
| 1460 |
|
---|
| 1461 | nr_ = new int[n_integration_center_];
|
---|
| 1462 |
|
---|
| 1463 | for (icenter=0; icenter<n_integration_center_; icenter++) {
|
---|
| 1464 | int iatom = mol_->non_q_atom(icenter);
|
---|
| 1465 | nr_[icenter]
|
---|
| 1466 | = ra_integrator_->get_radial_grid(mol_->Z(iatom),deriv_order)->nr();
|
---|
| 1467 | }
|
---|
| 1468 |
|
---|
| 1469 | centers_ = new SCVector3[n_integration_center_];
|
---|
| 1470 | for (icenter=0; icenter<n_integration_center_; icenter++) {
|
---|
| 1471 | int iatom = mol_->non_q_atom(icenter);
|
---|
| 1472 | centers_[icenter].x() = mol_->r(iatom,0);
|
---|
| 1473 | centers_[icenter].y() = mol_->r(iatom,1);
|
---|
| 1474 | centers_[icenter].z() = mol_->r(iatom,2);
|
---|
| 1475 | }
|
---|
| 1476 |
|
---|
| 1477 | atomic_radius_ = new double[n_integration_center_];
|
---|
| 1478 | for (icenter=0; icenter<n_integration_center_; icenter++) {
|
---|
| 1479 | int iatom = mol_->non_q_atom(icenter);
|
---|
| 1480 | atomic_radius_[icenter] = get_radius(mol_, iatom);
|
---|
| 1481 | }
|
---|
| 1482 |
|
---|
| 1483 | point_count_total_ = 0;
|
---|
| 1484 | total_density_ = 0.0;
|
---|
| 1485 | }
|
---|
| 1486 |
|
---|
| 1487 | RadialAngularIntegratorThread::~RadialAngularIntegratorThread()
|
---|
| 1488 | {
|
---|
| 1489 | delete[] centers_;
|
---|
| 1490 | delete[] atomic_radius_;
|
---|
| 1491 | delete[] nr_;
|
---|
| 1492 | }
|
---|
| 1493 |
|
---|
| 1494 | void
|
---|
| 1495 | RadialAngularIntegratorThread::run()
|
---|
| 1496 | {
|
---|
| 1497 | int icenter;
|
---|
| 1498 | int nangular;
|
---|
| 1499 | int ir, iangular; // Loop indices for diff. integration dim
|
---|
| 1500 | int point_count; // Counter for # integration points per center
|
---|
| 1501 | int nr;
|
---|
| 1502 |
|
---|
| 1503 | SCVector3 center; // Cartesian position of center
|
---|
| 1504 | SCVector3 integration_point;
|
---|
| 1505 |
|
---|
| 1506 | double w,radial_multiplier,angular_multiplier;
|
---|
| 1507 | int deriv_order = (nuclear_gradient_==0?0:1);
|
---|
| 1508 |
|
---|
| 1509 | int parallel_counter = 0;
|
---|
| 1510 |
|
---|
| 1511 | for (icenter=0; icenter < n_integration_center_; icenter++) {
|
---|
| 1512 | int iatom = mol_->non_q_atom(icenter);
|
---|
| 1513 | point_count=0;
|
---|
| 1514 | center = centers_[icenter];
|
---|
| 1515 | // get current radial grid: depends on convergence threshold
|
---|
| 1516 | RadialIntegrator *radial
|
---|
| 1517 | = ra_integrator_->get_radial_grid(mol_->Z(iatom), deriv_order);
|
---|
| 1518 | nr = radial->nr();
|
---|
| 1519 | for (ir=0; ir < nr; ir++) {
|
---|
| 1520 | if (! (parallel_counter++%nthread_ == ithread_)) continue;
|
---|
| 1521 | double r = radial->radial_value(ir, nr, atomic_radius_[icenter],
|
---|
| 1522 | radial_multiplier);
|
---|
| 1523 | // get current angular grid: depends on radial point and threshold
|
---|
| 1524 | AngularIntegrator *angular
|
---|
| 1525 | = ra_integrator_->get_angular_grid(r, atomic_radius_[icenter],
|
---|
| 1526 | mol_->Z(iatom),
|
---|
| 1527 | deriv_order);
|
---|
| 1528 | nangular = angular->num_angular_points(r/atomic_radius_[icenter],ir);
|
---|
| 1529 | for (iangular=0; iangular<nangular; iangular++) {
|
---|
| 1530 | angular_multiplier
|
---|
| 1531 | = angular->angular_point_cartesian(iangular,r,
|
---|
| 1532 | integration_point);
|
---|
| 1533 | integration_point += center;
|
---|
| 1534 | w=weight_->w(icenter, integration_point, w_gradient_);
|
---|
| 1535 | point_count++;
|
---|
| 1536 | double multiplier = angular_multiplier * radial_multiplier;
|
---|
| 1537 | total_density_
|
---|
| 1538 | += w * multiplier
|
---|
| 1539 | * do_point(iatom, integration_point,
|
---|
| 1540 | w, multiplier,
|
---|
| 1541 | nuclear_gradient_, f_gradient_, w_gradient_);
|
---|
| 1542 | }
|
---|
| 1543 | }
|
---|
| 1544 | point_count_total_ += point_count;
|
---|
| 1545 | }
|
---|
| 1546 | }
|
---|
| 1547 |
|
---|
| 1548 | //////////////////////////////////////////////
|
---|
| 1549 | // RadialAngularIntegrator
|
---|
| 1550 |
|
---|
| 1551 | static ClassDesc RadialAngularIntegrator_cd(
|
---|
| 1552 | typeid(RadialAngularIntegrator),"RadialAngularIntegrator",1,"public DenIntegrator",
|
---|
| 1553 | 0, create<RadialAngularIntegrator>, create<RadialAngularIntegrator>);
|
---|
| 1554 |
|
---|
| 1555 | RadialAngularIntegrator::RadialAngularIntegrator(StateIn& s):
|
---|
| 1556 | SavableState(s),
|
---|
| 1557 | DenIntegrator(s)
|
---|
| 1558 | {
|
---|
| 1559 | s.get(natomic_rows_);
|
---|
| 1560 | s.get(max_gridtype_);
|
---|
| 1561 | s.get(prune_grid_);
|
---|
| 1562 | s.get(gridtype_);
|
---|
| 1563 | s.get(npruned_partitions_);
|
---|
| 1564 | s.get(dynamic_grids_);
|
---|
| 1565 |
|
---|
| 1566 | // ExEnv::outn() << "natomic_rows_ = " << natomic_rows_ << endl;
|
---|
| 1567 | // ExEnv::outn() << "max_gridtype_ = " << max_gridtype_ << endl;
|
---|
| 1568 | // ExEnv::outn() << "prune_grid_ = " << prune_grid_ << endl;
|
---|
| 1569 | // ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
|
---|
| 1570 | // ExEnv::outn() << "npruned_partitions_ = " << npruned_partitions_ << endl;
|
---|
| 1571 | // ExEnv::outn() << "dynamic_grids_ = " << dynamic_grids_ << endl;
|
---|
| 1572 |
|
---|
| 1573 |
|
---|
| 1574 | // ExEnv::outn() << "In StateIn Constructor!" << endl;
|
---|
| 1575 | weight_ = new BeckeIntegrationWeight;
|
---|
| 1576 |
|
---|
| 1577 | int i;
|
---|
| 1578 | grid_accuracy_ = new double[max_gridtype_];
|
---|
| 1579 | grid_accuracy_[0] = 1e-4;
|
---|
| 1580 | for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;
|
---|
| 1581 |
|
---|
| 1582 | Alpha_coeffs_ = new_c_array2(natomic_rows_,npruned_partitions_-1,
|
---|
| 1583 | double(0));
|
---|
| 1584 | s.get_array_double(Alpha_coeffs_[0], natomic_rows_*(npruned_partitions_-1));
|
---|
| 1585 |
|
---|
| 1586 | radial_user_ << SavableState::restore_state(s);
|
---|
| 1587 | angular_user_ << SavableState::restore_state(s);
|
---|
| 1588 |
|
---|
| 1589 | init_default_grids();
|
---|
| 1590 | set_grids();
|
---|
| 1591 |
|
---|
| 1592 | }
|
---|
| 1593 |
|
---|
| 1594 | RadialAngularIntegrator::RadialAngularIntegrator()
|
---|
| 1595 | {
|
---|
| 1596 | weight_ = new BeckeIntegrationWeight;
|
---|
| 1597 |
|
---|
| 1598 | init_parameters();
|
---|
| 1599 | init_default_grids();
|
---|
| 1600 | set_grids();
|
---|
| 1601 |
|
---|
| 1602 | }
|
---|
| 1603 |
|
---|
| 1604 | RadialAngularIntegrator::RadialAngularIntegrator(const Ref<KeyVal>& keyval):
|
---|
| 1605 | DenIntegrator(keyval)
|
---|
| 1606 | {
|
---|
| 1607 |
|
---|
| 1608 | radial_user_ << keyval->describedclassvalue("radial");
|
---|
| 1609 | angular_user_ << keyval->describedclassvalue("angular");
|
---|
| 1610 |
|
---|
| 1611 | weight_ << keyval->describedclassvalue("weight");
|
---|
| 1612 | if (weight_.null()) weight_ = new BeckeIntegrationWeight;
|
---|
| 1613 | // ExEnv::outn() << "In Ref<KeyVal> Constructor" << endl;
|
---|
| 1614 |
|
---|
| 1615 | init_parameters(keyval);
|
---|
| 1616 | init_default_grids();
|
---|
| 1617 | set_grids();
|
---|
| 1618 |
|
---|
| 1619 | }
|
---|
| 1620 |
|
---|
| 1621 | RadialAngularIntegrator::~RadialAngularIntegrator()
|
---|
| 1622 | {
|
---|
| 1623 | delete_c_array2(Alpha_coeffs_);
|
---|
| 1624 | delete_c_array2(radial_grid_);
|
---|
| 1625 | delete_c_array3(angular_grid_);
|
---|
| 1626 | delete_c_array2(nr_points_);
|
---|
| 1627 | delete[] xcoarse_l_;
|
---|
| 1628 | delete[] grid_accuracy_;
|
---|
| 1629 | }
|
---|
| 1630 |
|
---|
| 1631 | void
|
---|
| 1632 | RadialAngularIntegrator::save_data_state(StateOut& s)
|
---|
| 1633 | {
|
---|
| 1634 | DenIntegrator::save_data_state(s);
|
---|
| 1635 | s.put(natomic_rows_);
|
---|
| 1636 | s.put(max_gridtype_);
|
---|
| 1637 | s.put(prune_grid_);
|
---|
| 1638 | s.put(gridtype_);
|
---|
| 1639 | // s.put(nr_points_[0], natomic_rows_*max_gridtype_);
|
---|
| 1640 | // s.put(xcoarse_l_, natomic_rows_);
|
---|
| 1641 | s.put(npruned_partitions_);
|
---|
| 1642 | s.put(dynamic_grids_);
|
---|
| 1643 | s.put_array_double(Alpha_coeffs_[0],natomic_rows_*(npruned_partitions_-1));
|
---|
| 1644 |
|
---|
| 1645 | // ExEnv::outn() << "natomic_rows_ = " << natomic_rows_ << endl;
|
---|
| 1646 | // ExEnv::outn() << "max_gridtype_ = " << max_gridtype_ << endl;
|
---|
| 1647 | // ExEnv::outn() << "prune_grid_ = " << prune_grid_ << endl;
|
---|
| 1648 | // ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
|
---|
| 1649 | // ExEnv::outn() << "npruned_partitions_ = " << npruned_partitions_ << endl;
|
---|
| 1650 | // ExEnv::outn() << "dynamic_grids_ = " << dynamic_grids_ << endl;
|
---|
| 1651 |
|
---|
| 1652 | SavableState::save_state(radial_user_.pointer(),s);
|
---|
| 1653 | SavableState::save_state(angular_user_.pointer(),s);
|
---|
| 1654 | }
|
---|
| 1655 |
|
---|
| 1656 | void
|
---|
| 1657 | RadialAngularIntegrator::init_parameters(void)
|
---|
| 1658 | {
|
---|
| 1659 |
|
---|
| 1660 | prune_grid_ = 1;
|
---|
| 1661 | gridtype_ = 3;
|
---|
| 1662 | npruned_partitions_ = 5;
|
---|
| 1663 | dynamic_grids_ = 1;
|
---|
| 1664 | max_gridtype_ = 6;
|
---|
| 1665 | natomic_rows_ = 5;
|
---|
| 1666 | grid_accuracy_ = new double[max_gridtype_];
|
---|
| 1667 |
|
---|
| 1668 | int i;
|
---|
| 1669 | grid_accuracy_[0] = 1e-4;
|
---|
| 1670 | for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;
|
---|
| 1671 |
|
---|
| 1672 | init_pruning_coefficients();
|
---|
| 1673 |
|
---|
| 1674 | // ExEnv::outn() << "gridtype_ = " << gridtype_ << endl;
|
---|
| 1675 |
|
---|
| 1676 | }
|
---|
| 1677 |
|
---|
| 1678 | void
|
---|
| 1679 | RadialAngularIntegrator::init_parameters(const Ref<KeyVal>& keyval)
|
---|
| 1680 | {
|
---|
| 1681 | char *grid = 0;
|
---|
| 1682 |
|
---|
| 1683 | max_gridtype_ = 6;
|
---|
| 1684 | natomic_rows_ = 5;
|
---|
| 1685 |
|
---|
| 1686 | grid = keyval->pcharvalue("grid");
|
---|
| 1687 |
|
---|
| 1688 | if (grid) {
|
---|
| 1689 | if (!strcmp(grid,"xcoarse")) gridtype_ = 0;
|
---|
| 1690 | else if (!strcmp(grid,"coarse")) gridtype_ = 1;
|
---|
| 1691 | else if (!strcmp(grid,"medium")) gridtype_ = 2;
|
---|
| 1692 | else if (!strcmp(grid,"fine")) gridtype_ = 3;
|
---|
| 1693 | else if (!strcmp(grid,"xfine")) gridtype_ = 4;
|
---|
| 1694 | else if (!strcmp(grid,"ultrafine")) gridtype_ = 5;
|
---|
| 1695 | else {
|
---|
| 1696 | ExEnv::out0()
|
---|
| 1697 | << indent
|
---|
| 1698 | << "ERROR: grid = \"" << grid << "\" not recognized."
|
---|
| 1699 | << endl
|
---|
| 1700 | << indent
|
---|
| 1701 | << "Choices are: xcoarse, coarse, medium, fine, xfine."
|
---|
| 1702 | << endl
|
---|
| 1703 | << indent
|
---|
| 1704 | << "The default is grid = fine."
|
---|
| 1705 | << endl;
|
---|
| 1706 | abort();
|
---|
| 1707 | }
|
---|
| 1708 |
|
---|
| 1709 | }
|
---|
| 1710 | else {
|
---|
| 1711 | gridtype_ = 3;
|
---|
| 1712 | }
|
---|
| 1713 |
|
---|
| 1714 |
|
---|
| 1715 |
|
---|
| 1716 | //ExEnv::outn() << " gridtype = " << gridtype_ << endl;
|
---|
| 1717 | //ExEnv::outn() << " max_gridtype = " << max_gridtype_ << endl;
|
---|
| 1718 | dynamic_grids_ = keyval->intvalue("dynamic");
|
---|
| 1719 | if (keyval->error() != KeyVal::OK) dynamic_grids_ = 1;
|
---|
| 1720 | grid_accuracy_ = new double[max_gridtype_];
|
---|
| 1721 | //ExEnv::outn() << "init_parameters:: max_gridtype_ = " << max_gridtype_;
|
---|
| 1722 |
|
---|
| 1723 | int i;
|
---|
| 1724 | grid_accuracy_[0] = 1e-4;
|
---|
| 1725 | for (i=1; i<max_gridtype_; i++) grid_accuracy_[i] = grid_accuracy_[i-1]*1e-1;
|
---|
| 1726 |
|
---|
| 1727 | init_pruning_coefficients(keyval);
|
---|
| 1728 |
|
---|
| 1729 | delete[] grid;
|
---|
| 1730 | }
|
---|
| 1731 |
|
---|
| 1732 |
|
---|
| 1733 | void
|
---|
| 1734 | RadialAngularIntegrator::set_grids(void)
|
---|
| 1735 | {
|
---|
| 1736 | int i, j, k;
|
---|
| 1737 |
|
---|
| 1738 | radial_grid_ = new_c_array2(natomic_rows_,gridtype_+1,
|
---|
| 1739 | Ref<RadialIntegrator>());
|
---|
| 1740 | angular_grid_ = new_c_array3(natomic_rows_, npruned_partitions_,
|
---|
| 1741 | gridtype_+1, Ref<AngularIntegrator>());
|
---|
| 1742 |
|
---|
| 1743 | int prune_formula_1[5] = {26, 18, 12, 0, 12}; // for H to Ne
|
---|
| 1744 | int prune_formula_2[5] = {36, 24, 12, 0, 12}; // for Na and up
|
---|
| 1745 |
|
---|
| 1746 | double prune_factor[5] = {0.1, 0.4, 0.6, 1.0, 0.6};
|
---|
| 1747 |
|
---|
| 1748 | if (npruned_partitions_ == 1) {
|
---|
| 1749 | prune_factor[0] = 1.0;
|
---|
| 1750 | prune_formula_1[0] = 0;
|
---|
| 1751 | prune_formula_2[0] = 0;
|
---|
| 1752 | }
|
---|
| 1753 | else if (npruned_partitions_ != 5) {
|
---|
| 1754 | ExEnv::errn() << "RadialAngularIntegrator::set_grids: "
|
---|
| 1755 | << "npruned_partitations must be 1 or 5" << endl;
|
---|
| 1756 | abort();
|
---|
| 1757 | }
|
---|
| 1758 |
|
---|
| 1759 | for (i=0; i<natomic_rows_; i++) {
|
---|
| 1760 | for (j=0; j<=gridtype_; j++)
|
---|
| 1761 | radial_grid_[i][j]
|
---|
| 1762 | = new EulerMaclaurinRadialIntegrator(nr_points_[i][j]);
|
---|
| 1763 |
|
---|
| 1764 | for (j=0; j<npruned_partitions_; j++) {
|
---|
| 1765 | for (k=0; k<=gridtype_; k++) {
|
---|
| 1766 | int grid_l = xcoarse_l_[i] + angular_grid_offset(k);
|
---|
| 1767 | int use_l;
|
---|
| 1768 |
|
---|
| 1769 | // the constant shift depends on the row and the partition
|
---|
| 1770 | int prune_formula;
|
---|
| 1771 | if (i<=1) prune_formula = prune_formula_1[j]; // H to Ne
|
---|
| 1772 | else prune_formula = prune_formula_2[j]; // Na and up
|
---|
| 1773 |
|
---|
| 1774 | // Compute the l to be used from both the constant shift and
|
---|
| 1775 | // the multiplicative factor method. Use the method giving
|
---|
| 1776 | // the highest l.
|
---|
| 1777 | int use_l_formula = grid_l - prune_formula;
|
---|
| 1778 | int use_l_factor = int(grid_l*prune_factor[j] + 0.5);
|
---|
| 1779 | if (use_l_formula > use_l_factor) use_l = use_l_formula;
|
---|
| 1780 | else use_l = use_l_factor;
|
---|
| 1781 |
|
---|
| 1782 | angular_grid_[i][j][k]
|
---|
| 1783 | = new LebedevLaikovIntegrator(Lebedev_Laikov_npoint(use_l));
|
---|
| 1784 | // ExEnv::outn() << " angular_grid_["
|
---|
| 1785 | // << i << "]["
|
---|
| 1786 | // << j << "]["
|
---|
| 1787 | // << k
|
---|
| 1788 | // << "]->nw = " << angular_grid_[i][j][k]->nw()
|
---|
| 1789 | // << " xc_l = " << xcoarse_l_[i]
|
---|
| 1790 | // << " off = " << angular_grid_offset(k)
|
---|
| 1791 | // << " grid_l = " << grid_l
|
---|
| 1792 | // << " use_l = " << use_l
|
---|
| 1793 | // << endl;
|
---|
| 1794 | }
|
---|
| 1795 | }
|
---|
| 1796 | }
|
---|
| 1797 | }
|
---|
| 1798 |
|
---|
| 1799 | void
|
---|
| 1800 | RadialAngularIntegrator::init_pruning_coefficients(void)
|
---|
| 1801 | {
|
---|
| 1802 | // Set up Alpha arrays for pruning
|
---|
| 1803 | //ExEnv::outn() << "npruned_partitions = " << npruned_partitions_ << endl;
|
---|
| 1804 | //ExEnv::outn() << "natomic_rows = " << natomic_rows_ << endl;
|
---|
| 1805 | int num_boundaries = npruned_partitions_-1;
|
---|
| 1806 | Alpha_coeffs_ = new_zero_c_array2(natomic_rows_, num_boundaries,
|
---|
| 1807 | double(0));
|
---|
| 1808 |
|
---|
| 1809 | // set pruning cutoff variables - Alpha -> radial shell cutoffs
|
---|
| 1810 | init_alpha_coefficients();
|
---|
| 1811 |
|
---|
| 1812 | }
|
---|
| 1813 |
|
---|
| 1814 | void
|
---|
| 1815 | RadialAngularIntegrator::init_pruning_coefficients(const Ref<KeyVal>& keyval)
|
---|
| 1816 | {
|
---|
| 1817 | int i, j;
|
---|
| 1818 |
|
---|
| 1819 | prune_grid_ = keyval->booleanvalue("prune_grid");
|
---|
| 1820 | if (keyval->error() != KeyVal::OK) prune_grid_ = 1;
|
---|
| 1821 |
|
---|
| 1822 | // ExEnv::outn() << "prune_grid = " << prune_grid_ << endl;
|
---|
| 1823 |
|
---|
| 1824 | // Need to generalize this to parse input for any number of grids
|
---|
| 1825 | if (prune_grid_) {
|
---|
| 1826 | npruned_partitions_ = keyval->count("angular_points");
|
---|
| 1827 | if (keyval->error() != KeyVal::OK) npruned_partitions_ = 5;
|
---|
| 1828 |
|
---|
| 1829 | // set pruning cutoff variables - Alpha -> radial shell cutoffs
|
---|
| 1830 | int num_boundaries = npruned_partitions_-1;
|
---|
| 1831 | Alpha_coeffs_ = new_zero_c_array2(natomic_rows_, num_boundaries,
|
---|
| 1832 | double(0));
|
---|
| 1833 | int alpha_rows = keyval->count("alpha_coeffs");
|
---|
| 1834 | if (keyval->error() != KeyVal::OK) {
|
---|
| 1835 | if (npruned_partitions_ != 5) {
|
---|
| 1836 | ExEnv::outn() << " RadialAngularIntegrator:: Need to supply alpha coefficients "
|
---|
| 1837 | << "for the " << num_boundaries << " partition boundaries " << endl;
|
---|
| 1838 | abort();
|
---|
| 1839 | }
|
---|
| 1840 | init_alpha_coefficients();
|
---|
| 1841 | }
|
---|
| 1842 | else { // alpha coefficients defined in input
|
---|
| 1843 | int check;
|
---|
| 1844 | for (i=0; i<alpha_rows; i++) {
|
---|
| 1845 | check = keyval->count("alpha_coeffs", i);
|
---|
| 1846 | if (check != num_boundaries) {
|
---|
| 1847 | ExEnv::outn() << "RadialAngularIntegrator:: Number of alpha coefficients does "
|
---|
| 1848 | << "not match the number of boundaries (" << check << " != "
|
---|
| 1849 | << num_boundaries << ")" << endl;
|
---|
| 1850 | abort();
|
---|
| 1851 | }
|
---|
| 1852 | for (j=0; j<num_boundaries; j++)
|
---|
| 1853 | Alpha_coeffs_[i][j] = keyval->doublevalue("alpha_coeffs", i, j);
|
---|
| 1854 | }
|
---|
| 1855 | }
|
---|
| 1856 | }
|
---|
| 1857 | else {
|
---|
| 1858 | npruned_partitions_ = 1;
|
---|
| 1859 | Alpha_coeffs_ = new_zero_c_array2(natomic_rows_,0,
|
---|
| 1860 | double(0));
|
---|
| 1861 | }
|
---|
| 1862 | }
|
---|
| 1863 |
|
---|
| 1864 | void
|
---|
| 1865 | RadialAngularIntegrator::init_alpha_coefficients(void)
|
---|
| 1866 | {
|
---|
| 1867 | // assumes Alpha_coeffs_ is allocated and zeroed.
|
---|
| 1868 |
|
---|
| 1869 | Alpha_coeffs_[0][0] = 0.25; Alpha_coeffs_[0][1] = 0.5;
|
---|
| 1870 | Alpha_coeffs_[0][2] = 0.9; Alpha_coeffs_[0][3] = 4.5;
|
---|
| 1871 | Alpha_coeffs_[1][0] = 0.1667; Alpha_coeffs_[1][1] = 0.5;
|
---|
| 1872 | Alpha_coeffs_[1][2] = 0.8; Alpha_coeffs_[1][3] = 4.5;
|
---|
| 1873 | Alpha_coeffs_[2][0] = 0.1; Alpha_coeffs_[2][1] = 0.4;
|
---|
| 1874 | Alpha_coeffs_[2][2] = 0.7; Alpha_coeffs_[2][3] = 2.5;
|
---|
| 1875 |
|
---|
| 1876 | // No pruning for atoms past second row
|
---|
| 1877 | int i;
|
---|
| 1878 | for (i=3; i<natomic_rows_; i++) Alpha_coeffs_[i][2] = DBL_MAX;
|
---|
| 1879 |
|
---|
| 1880 | }
|
---|
| 1881 |
|
---|
| 1882 | void
|
---|
| 1883 | RadialAngularIntegrator::init_default_grids(void)
|
---|
| 1884 | {
|
---|
| 1885 | xcoarse_l_ = new int[natomic_rows_];
|
---|
| 1886 |
|
---|
| 1887 | nr_points_ = new_c_array2(natomic_rows_,max_gridtype_,int(0));
|
---|
| 1888 |
|
---|
| 1889 | // Set angular momentum level of reference xcoarse grids for each atomic row
|
---|
| 1890 | xcoarse_l_[0] = 17; xcoarse_l_[1] = 17; xcoarse_l_[2] = 21;
|
---|
| 1891 | xcoarse_l_[3] = 25; xcoarse_l_[4] = 31;
|
---|
| 1892 |
|
---|
| 1893 | // Set number of radial points for each atomic row and grid type
|
---|
| 1894 | nr_points_[0][0] = 30; nr_points_[0][1] = 50; nr_points_[0][2] = 75;
|
---|
| 1895 | nr_points_[0][3] = 85; nr_points_[0][4] = 100; nr_points_[0][5] = 140;
|
---|
| 1896 |
|
---|
| 1897 | nr_points_[1][0] = 30; nr_points_[1][1] = 50; nr_points_[1][2] = 75;
|
---|
| 1898 | nr_points_[1][3] = 85; nr_points_[1][4] = 100; nr_points_[1][5] = 140;
|
---|
| 1899 |
|
---|
| 1900 | nr_points_[2][0] = 45; nr_points_[2][1] = 75; nr_points_[2][2] = 95;
|
---|
| 1901 | nr_points_[2][3] = 110; nr_points_[2][4] = 125; nr_points_[2][5] = 175;
|
---|
| 1902 |
|
---|
| 1903 | nr_points_[3][0] = 75; nr_points_[3][1] = 95; nr_points_[3][2] = 110;
|
---|
| 1904 | nr_points_[3][3] = 130; nr_points_[3][4] = 160; nr_points_[3][5] = 210;
|
---|
| 1905 |
|
---|
| 1906 | nr_points_[4][0] = 105; nr_points_[4][1] = 130; nr_points_[4][2] = 155;
|
---|
| 1907 | nr_points_[4][3] = 180; nr_points_[4][4] = 205; nr_points_[4][5] = 235;
|
---|
| 1908 |
|
---|
| 1909 | // prune_grid_ = 1; npruned_partitions_ = 5; gridtype_ = 2;
|
---|
| 1910 | }
|
---|
| 1911 |
|
---|
| 1912 | int
|
---|
| 1913 | RadialAngularIntegrator::angular_grid_offset(int gridtype)
|
---|
| 1914 | {
|
---|
| 1915 | switch (gridtype) {
|
---|
| 1916 | case 0: return 0;
|
---|
| 1917 | case 1: return 6;
|
---|
| 1918 | case 2: return 12;
|
---|
| 1919 | case 3: return 18;
|
---|
| 1920 | case 4: return 30;
|
---|
| 1921 | case 5: return 42;
|
---|
| 1922 | default: abort();
|
---|
| 1923 | }
|
---|
| 1924 | return 0;
|
---|
| 1925 | }
|
---|
| 1926 |
|
---|
| 1927 | RadialIntegrator *
|
---|
| 1928 | RadialAngularIntegrator::get_radial_grid(int charge, int deriv_order)
|
---|
| 1929 | {
|
---|
| 1930 | if (radial_user_.null()) {
|
---|
| 1931 |
|
---|
| 1932 | int select_grid;
|
---|
| 1933 |
|
---|
| 1934 | if (dynamic_grids_ && deriv_order == 0)
|
---|
| 1935 | select_grid = select_dynamic_grid();
|
---|
| 1936 | else select_grid = gridtype_;
|
---|
| 1937 | //ExEnv::out << "RAI::get_radial_grid -> select_grid = " << select_grid;
|
---|
| 1938 |
|
---|
| 1939 | if (charge<3) return radial_grid_[0][select_grid].pointer();
|
---|
| 1940 | else if (charge<11) return radial_grid_[1][select_grid].pointer();
|
---|
| 1941 | else if (charge<19) return radial_grid_[2][select_grid].pointer();
|
---|
| 1942 | else if (charge<37) return radial_grid_[3][select_grid].pointer();
|
---|
| 1943 | else if (charge<55) return radial_grid_[4][select_grid].pointer();
|
---|
| 1944 | else {
|
---|
| 1945 | ExEnv::outn() << " No default radial grids for atomic charge " << charge << endl;
|
---|
| 1946 | abort();
|
---|
| 1947 | }
|
---|
| 1948 | }
|
---|
| 1949 |
|
---|
| 1950 | return radial_user_.pointer();
|
---|
| 1951 |
|
---|
| 1952 | }
|
---|
| 1953 |
|
---|
| 1954 | int
|
---|
| 1955 | RadialAngularIntegrator::select_dynamic_grid(void)
|
---|
| 1956 | {
|
---|
| 1957 | double accuracy = get_accuracy();
|
---|
| 1958 | // accurate_grid = gridtype_ to get original non-dynamic version
|
---|
| 1959 | // ExEnv::out << " accuracy = " << accuracy << endl;
|
---|
| 1960 | int select_grid;
|
---|
| 1961 | int i;
|
---|
| 1962 |
|
---|
| 1963 | if (accuracy >= grid_accuracy_[0]) select_grid=0;
|
---|
| 1964 | else if (accuracy <= grid_accuracy_[gridtype_]) select_grid=gridtype_;
|
---|
| 1965 | else {
|
---|
| 1966 | for (i=gridtype_; i>=0; i--)
|
---|
| 1967 | if (accuracy >= grid_accuracy_[i]) select_grid=i;
|
---|
| 1968 | }
|
---|
| 1969 |
|
---|
| 1970 | // ExEnv::out << " select_grid = " << select_grid << endl;
|
---|
| 1971 | return select_grid;
|
---|
| 1972 | }
|
---|
| 1973 |
|
---|
| 1974 | int
|
---|
| 1975 | RadialAngularIntegrator::get_atomic_row(int i)
|
---|
| 1976 | {
|
---|
| 1977 | if (i<3) return 0;
|
---|
| 1978 | else if (i<11) return 1;
|
---|
| 1979 | else if (i<19) return 2;
|
---|
| 1980 | else if (i<37) return 3;
|
---|
| 1981 | else if (i<55) return 4;
|
---|
| 1982 | else if (i<87) return 5;
|
---|
| 1983 |
|
---|
| 1984 | ExEnv::outn() << " RadialAngularIntegrator::get_atomic_row: Z too large: "
|
---|
| 1985 | << i << endl;
|
---|
| 1986 | abort();
|
---|
| 1987 | return 0;
|
---|
| 1988 | }
|
---|
| 1989 |
|
---|
| 1990 | AngularIntegrator *
|
---|
| 1991 | RadialAngularIntegrator::get_angular_grid(double radius, double atomic_radius,
|
---|
| 1992 | int Z, int deriv_order)
|
---|
| 1993 | {
|
---|
| 1994 | int atomic_row, i;
|
---|
| 1995 |
|
---|
| 1996 | int select_grid;
|
---|
| 1997 | if (dynamic_grids_ && deriv_order == 0) select_grid = select_dynamic_grid();
|
---|
| 1998 | else select_grid = gridtype_;
|
---|
| 1999 |
|
---|
| 2000 | //ExEnv::out << "RAI::get_angular_grid -> select_grid = " << select_grid;
|
---|
| 2001 | //ExEnv::out << " prune_grid_ = " << prune_grid_ << endl;
|
---|
| 2002 | atomic_row = get_atomic_row(Z);
|
---|
| 2003 | if (angular_user_.null()) {
|
---|
| 2004 | // Seleted Alpha's
|
---|
| 2005 | double *Alpha = Alpha_coeffs_[atomic_row];
|
---|
| 2006 | // gridtype_ will need to be adjusted for dynamic grids
|
---|
| 2007 | for (i=0; i<npruned_partitions_-1; i++) {
|
---|
| 2008 | if (radius/atomic_radius < Alpha[i]) {
|
---|
| 2009 | return angular_grid_[atomic_row][i][select_grid].pointer();
|
---|
| 2010 | }
|
---|
| 2011 | }
|
---|
| 2012 | return angular_grid_[atomic_row][npruned_partitions_-1][select_grid]
|
---|
| 2013 | .pointer();
|
---|
| 2014 | }
|
---|
| 2015 | else {
|
---|
| 2016 | return angular_user_.pointer();
|
---|
| 2017 | }
|
---|
| 2018 | }
|
---|
| 2019 |
|
---|
| 2020 | void
|
---|
| 2021 | RadialAngularIntegrator::integrate(const Ref<DenFunctional> &denfunc,
|
---|
| 2022 | const RefSymmSCMatrix& densa,
|
---|
| 2023 | const RefSymmSCMatrix& densb,
|
---|
| 2024 | double *nuclear_gradient)
|
---|
| 2025 | {
|
---|
| 2026 | int i;
|
---|
| 2027 |
|
---|
| 2028 | tim_enter("integrate");
|
---|
| 2029 |
|
---|
| 2030 | init_integration(denfunc, densa, densb, nuclear_gradient);
|
---|
| 2031 |
|
---|
| 2032 | weight_->init(wavefunction()->molecule(), DBL_EPSILON);
|
---|
| 2033 |
|
---|
| 2034 | int me = messagegrp_->me();
|
---|
| 2035 | int nthread = threadgrp_->nthread();
|
---|
| 2036 | int nthread_overall = nthread;
|
---|
| 2037 | messagegrp_->sum(nthread_overall);
|
---|
| 2038 | int ithread_overall = 0;
|
---|
| 2039 | if (me > 0) {
|
---|
| 2040 | messagegrp_->recv(me - 1,ithread_overall);
|
---|
| 2041 | }
|
---|
| 2042 | if (me < messagegrp_->n() - 1) {
|
---|
| 2043 | int ithread_overall_next = ithread_overall + nthread;
|
---|
| 2044 | messagegrp_->send(me + 1, ithread_overall_next);
|
---|
| 2045 | }
|
---|
| 2046 |
|
---|
| 2047 | // create threads
|
---|
| 2048 | //cout << "creating test lock" << endl;
|
---|
| 2049 | //Ref<ThreadLock> reflock = threadgrp_->new_lock();
|
---|
| 2050 | //tlock = reflock.pointer();
|
---|
| 2051 | RadialAngularIntegratorThread **threads =
|
---|
| 2052 | new RadialAngularIntegratorThread*[nthread];
|
---|
| 2053 | for (i=0; i<nthread; i++) {
|
---|
| 2054 | Ref<BatchElectronDensity> bed = new BatchElectronDensity(den_, true);
|
---|
| 2055 | if (nuclear_gradient != 0) {
|
---|
| 2056 | bed->set_need_basis_gradient(true);
|
---|
| 2057 | if (denfunc->need_density_gradient()) bed->set_need_basis_hessian(true);
|
---|
| 2058 | }
|
---|
| 2059 | threads[i] = new RadialAngularIntegratorThread(
|
---|
| 2060 | i + ithread_overall, nthread_overall,
|
---|
| 2061 | this, denfunc.pointer(),
|
---|
| 2062 | bed,
|
---|
| 2063 | linear_scaling_, use_dmat_bound_,
|
---|
| 2064 | accuracy_, compute_potential_integrals_,
|
---|
| 2065 | nuclear_gradient != 0);
|
---|
| 2066 | threadgrp_->add_thread(i, threads[i]);
|
---|
| 2067 | }
|
---|
| 2068 |
|
---|
| 2069 | //for (i=0; i<nthread; i++) {
|
---|
| 2070 | // ExEnv::outn() << "Intial Thread " << i
|
---|
| 2071 | // << ": point count = " << threads[i]->point_count()
|
---|
| 2072 | // << " total density = " << threads[i]->total_density()
|
---|
| 2073 | // << " value = " << threads[i]->value()
|
---|
| 2074 | // << endl;
|
---|
| 2075 | //}
|
---|
| 2076 |
|
---|
| 2077 | // run threads
|
---|
| 2078 | threadgrp_->start_threads();
|
---|
| 2079 | threadgrp_->wait_threads();
|
---|
| 2080 | //for (i=0; i<nthread; i++) {
|
---|
| 2081 | // cout << "Running thread " << i << endl;
|
---|
| 2082 | // threads[i]->run();
|
---|
| 2083 | // }
|
---|
| 2084 |
|
---|
| 2085 | // sum results
|
---|
| 2086 | int point_count_total = 0;
|
---|
| 2087 | double total_density = 0.0;
|
---|
| 2088 | value_ = 0.0;
|
---|
| 2089 | for (i=0; i<nthread; i++) {
|
---|
| 2090 | //ExEnv::outn() << "Thread " << i
|
---|
| 2091 | // << ": point count = " << threads[i]->point_count()
|
---|
| 2092 | // << " total density = " << threads[i]->total_density()
|
---|
| 2093 | // << " value = " << threads[i]->value()
|
---|
| 2094 | // << endl;
|
---|
| 2095 | point_count_total += threads[i]->point_count();
|
---|
| 2096 | total_density += threads[i]->total_density();
|
---|
| 2097 | value_ += threads[i]->value();
|
---|
| 2098 | if (compute_potential_integrals_) {
|
---|
| 2099 | int ntri = (nbasis_*(nbasis_+1))/2;
|
---|
| 2100 | double *alpha_vmat_i = threads[i]->alpha_vmat();
|
---|
| 2101 | for (int j=0; j<ntri; j++) alpha_vmat_[j] += alpha_vmat_i[j];
|
---|
| 2102 | if (spin_polarized_) {
|
---|
| 2103 | double *beta_vmat_i = threads[i]->beta_vmat();
|
---|
| 2104 | for (int j=0; j<ntri; j++) beta_vmat_[j] += beta_vmat_i[j];
|
---|
| 2105 | }
|
---|
| 2106 | }
|
---|
| 2107 | if (nuclear_gradient != 0) {
|
---|
| 2108 | int natom3 = 3 * wavefunction()->molecule()->natom();
|
---|
| 2109 | double *th_nuclear_gradient = threads[i]->nuclear_gradient();
|
---|
| 2110 | for (int j=0; j<natom3; j++) {
|
---|
| 2111 | nuclear_gradient[j] += th_nuclear_gradient[j];
|
---|
| 2112 | }
|
---|
| 2113 | }
|
---|
| 2114 | }
|
---|
| 2115 |
|
---|
| 2116 | threadgrp_->delete_threads();
|
---|
| 2117 | delete[] threads;
|
---|
| 2118 |
|
---|
| 2119 | messagegrp_->sum(point_count_total);
|
---|
| 2120 | messagegrp_->sum(total_density);
|
---|
| 2121 | done_integration();
|
---|
| 2122 | weight_->done();
|
---|
| 2123 |
|
---|
| 2124 | ExEnv::out0() << indent
|
---|
| 2125 | << "Total integration points = " << point_count_total << endl;
|
---|
| 2126 | ExEnv::out0() << indent
|
---|
| 2127 | << "Integrated electron density error = "
|
---|
| 2128 | << scprintf("%14.12f", total_density-wfn_->nelectron())
|
---|
| 2129 | << endl;
|
---|
| 2130 |
|
---|
| 2131 | tim_exit("integrate");
|
---|
| 2132 | }
|
---|
| 2133 |
|
---|
| 2134 | void
|
---|
| 2135 | RadialAngularIntegrator::print(ostream &o) const
|
---|
| 2136 | {
|
---|
| 2137 | o << indent << class_name() << ":" << endl;
|
---|
| 2138 | o << incindent;
|
---|
| 2139 | if (radial_user_.nonnull()) {
|
---|
| 2140 | o << indent << "User defined radial grid:" << endl;
|
---|
| 2141 | o << incindent;
|
---|
| 2142 | radial_user_->print(o);
|
---|
| 2143 | o << decindent;
|
---|
| 2144 | }
|
---|
| 2145 | if (angular_user_.nonnull()) {
|
---|
| 2146 | o << indent << "User defined angular grid:" << endl;
|
---|
| 2147 | o << incindent;
|
---|
| 2148 | angular_user_->print(o);
|
---|
| 2149 | o << decindent;
|
---|
| 2150 | }
|
---|
| 2151 | if (angular_user_.null() || radial_user_.null()) {
|
---|
| 2152 | if (prune_grid_) o << indent << "Pruned ";
|
---|
| 2153 | switch (gridtype_) {
|
---|
| 2154 | case 0: o << "xcoarse"; break;
|
---|
| 2155 | case 1: o << "coarse"; break;
|
---|
| 2156 | case 2: o << "medium"; break;
|
---|
| 2157 | case 3: o << "fine"; break;
|
---|
| 2158 | case 4: o << "xfine"; break;
|
---|
| 2159 | case 5: o << "ultrafine"; break;
|
---|
| 2160 | default: o << "unknown"; break;
|
---|
| 2161 | }
|
---|
| 2162 | o << " grid employed" << endl;
|
---|
| 2163 | }
|
---|
| 2164 |
|
---|
| 2165 | o << decindent;
|
---|
| 2166 | }
|
---|
| 2167 |
|
---|
| 2168 | /////////////////////////////////////////////////////////////////////////////
|
---|
| 2169 |
|
---|
| 2170 | }
|
---|
| 2171 |
|
---|
| 2172 | // Local Variables:
|
---|
| 2173 | // mode: c++
|
---|
| 2174 | // c-file-style: "CLJ"
|
---|
| 2175 | // End:
|
---|
| 2176 |
|
---|