[0b990d] | 1 | //
|
---|
| 2 | // integral.cc --- implementation of the Integral class
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 1996 Limit Point Systems, Inc.
|
---|
| 5 | //
|
---|
| 6 | // Author: Edward Seidl <seidl@janed.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 <stdexcept>
|
---|
| 33 | #include <sstream>
|
---|
| 34 |
|
---|
| 35 | #include <util/state/stateio.h>
|
---|
| 36 | #include <chemistry/qc/basis/integral.h>
|
---|
| 37 | #include <chemistry/qc/basis/shellrot.h>
|
---|
| 38 | #include <chemistry/qc/basis/petite.h>
|
---|
| 39 | #include <chemistry/qc/basis/obint.h>
|
---|
| 40 |
|
---|
| 41 | using namespace std;
|
---|
| 42 | using namespace sc;
|
---|
| 43 |
|
---|
| 44 | static ClassDesc Integral_cd(
|
---|
| 45 | typeid(Integral),"Integral",2,"public SavableState",
|
---|
| 46 | 0, 0, 0);
|
---|
| 47 |
|
---|
| 48 | Integral::Integral(const Ref<GaussianBasisSet> &b1,
|
---|
| 49 | const Ref<GaussianBasisSet> &b2,
|
---|
| 50 | const Ref<GaussianBasisSet> &b3,
|
---|
| 51 | const Ref<GaussianBasisSet> &b4)
|
---|
| 52 | {
|
---|
| 53 | storage_ = 0;
|
---|
| 54 | storage_used_ = 0;
|
---|
| 55 | grp_ = MessageGrp::get_default_messagegrp();
|
---|
| 56 | set_basis(b1,b2,b3,b4);
|
---|
| 57 | }
|
---|
| 58 |
|
---|
| 59 | Integral::~Integral()
|
---|
| 60 | {
|
---|
| 61 | }
|
---|
| 62 |
|
---|
| 63 | Integral::Integral(StateIn& s) :
|
---|
| 64 | SavableState(s)
|
---|
| 65 | {
|
---|
| 66 | storage_used_ = 0;
|
---|
| 67 | bs1_ << SavableState::restore_state(s);
|
---|
| 68 | bs2_ << SavableState::restore_state(s);
|
---|
| 69 | bs3_ << SavableState::restore_state(s);
|
---|
| 70 | bs4_ << SavableState::restore_state(s);
|
---|
| 71 | if (s.version(::class_desc<Integral>()) >= 2) {
|
---|
| 72 | double dstorage;
|
---|
| 73 | s.get(dstorage);
|
---|
| 74 | storage_ = size_t(dstorage);
|
---|
| 75 | }
|
---|
| 76 | else {
|
---|
| 77 | unsigned int istorage;
|
---|
| 78 | s.get(istorage);
|
---|
| 79 | storage_ = istorage;
|
---|
| 80 | }
|
---|
| 81 | grp_ = MessageGrp::get_default_messagegrp();
|
---|
| 82 | }
|
---|
| 83 |
|
---|
| 84 | Integral::Integral(const Ref<KeyVal>&)
|
---|
| 85 | {
|
---|
| 86 | storage_used_ = 0;
|
---|
| 87 | storage_ = 0;
|
---|
| 88 | grp_ = MessageGrp::get_default_messagegrp();
|
---|
| 89 | }
|
---|
| 90 |
|
---|
| 91 | void
|
---|
| 92 | Integral::save_data_state(StateOut&o)
|
---|
| 93 | {
|
---|
| 94 | SavableState::save_state(bs1_.pointer(),o);
|
---|
| 95 | SavableState::save_state(bs2_.pointer(),o);
|
---|
| 96 | SavableState::save_state(bs3_.pointer(),o);
|
---|
| 97 | SavableState::save_state(bs4_.pointer(),o);
|
---|
| 98 | double dstorage = storage_;
|
---|
| 99 | o.put(dstorage);
|
---|
| 100 | }
|
---|
| 101 |
|
---|
| 102 | Ref<Integral> default_integral;
|
---|
| 103 |
|
---|
| 104 | void
|
---|
| 105 | Integral::set_default_integral(const Ref<Integral>& intf)
|
---|
| 106 | {
|
---|
| 107 | default_integral = intf;
|
---|
| 108 | }
|
---|
| 109 |
|
---|
| 110 | // Liberally borrowed from ThreadGrp
|
---|
| 111 | Integral*
|
---|
| 112 | Integral::initial_integral(int& argc, char ** argv)
|
---|
| 113 | {
|
---|
| 114 | Integral *intf = 0;
|
---|
| 115 | char * keyval_string = 0;
|
---|
| 116 |
|
---|
| 117 | // see if an integral factory is given on the command line
|
---|
| 118 | if (argc && argv) {
|
---|
| 119 | for (int i=0; i < argc; i++) {
|
---|
| 120 | if (argv[i] && !strcmp(argv[i], "-integral")) {
|
---|
| 121 | char *integral_string = argv[i];
|
---|
| 122 | i++;
|
---|
| 123 | if (i >= argc) {
|
---|
| 124 | throw runtime_error("-integral must be followed by an argument");
|
---|
| 125 | }
|
---|
| 126 | keyval_string = argv[i];
|
---|
| 127 | // move the integral arguments to the end of argv
|
---|
| 128 | int j;
|
---|
| 129 | for (j=i+1; j<argc; j++) {
|
---|
| 130 | argv[j-2] = argv[j];
|
---|
| 131 | }
|
---|
| 132 | argv[j++] = integral_string;
|
---|
| 133 | argv[j++] = keyval_string;
|
---|
| 134 | // decrement argc to hide the last two arguments
|
---|
| 135 | argc -= 2;
|
---|
| 136 | break;
|
---|
| 137 | }
|
---|
| 138 | }
|
---|
| 139 | }
|
---|
| 140 |
|
---|
| 141 | if (!keyval_string) {
|
---|
| 142 | // find out if the environment gives the containing integral
|
---|
| 143 | keyval_string = getenv("INTEGRAL");
|
---|
| 144 | if (keyval_string) {
|
---|
| 145 | if (!strncmp("INTEGRAL=", keyval_string, 11)) {
|
---|
| 146 | keyval_string = strchr(keyval_string, '=');
|
---|
| 147 | }
|
---|
| 148 | if (*keyval_string == '=') keyval_string++;
|
---|
| 149 | }
|
---|
| 150 | }
|
---|
| 151 |
|
---|
| 152 | // if keyval input for a integral was found, then
|
---|
| 153 | // create it.
|
---|
| 154 | if (keyval_string) {
|
---|
| 155 | if (keyval_string[0] == '\0') return 0;
|
---|
| 156 | Ref<ParsedKeyVal> strkv = new ParsedKeyVal();
|
---|
| 157 | strkv->parse_string(keyval_string);
|
---|
| 158 | Ref<DescribedClass> dc = strkv->describedclassvalue();
|
---|
| 159 | intf = dynamic_cast<Integral*>(dc.pointer());
|
---|
| 160 | if (dc.null()) {
|
---|
| 161 | ostringstream errmsg;
|
---|
| 162 | errmsg << "Integral::initial_integral: couldn't find a Integral in " << keyval_string << ends;
|
---|
| 163 | throw runtime_error(errmsg.str());
|
---|
| 164 | } else if (!intf) {
|
---|
| 165 | ostringstream errmsg;
|
---|
| 166 | errmsg << "initial_integral: wanted Integral but got " << dc->class_name() << ends;
|
---|
| 167 | throw runtime_error(errmsg.str());
|
---|
| 168 | }
|
---|
| 169 | // prevent an accidental delete
|
---|
| 170 | intf->reference();
|
---|
| 171 | strkv = 0;
|
---|
| 172 | dc = 0;
|
---|
| 173 | // accidental delete not a problem anymore since all smart pointers
|
---|
| 174 | // to intf are dead
|
---|
| 175 | intf->dereference();
|
---|
| 176 | return intf;
|
---|
| 177 | }
|
---|
| 178 |
|
---|
| 179 | return 0;
|
---|
| 180 | }
|
---|
| 181 |
|
---|
| 182 | int
|
---|
| 183 | Integral::equiv(const Ref<Integral> &integral)
|
---|
| 184 | {
|
---|
| 185 | return eq(class_desc(),integral->class_desc());
|
---|
| 186 | }
|
---|
| 187 |
|
---|
| 188 | Ref<PetiteList>
|
---|
| 189 | Integral::petite_list()
|
---|
| 190 | {
|
---|
| 191 | return new PetiteList(bs1_, this);
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | Ref<PetiteList>
|
---|
| 195 | Integral::petite_list(const Ref<GaussianBasisSet>& gbs)
|
---|
| 196 | {
|
---|
| 197 | return new PetiteList(gbs, this);
|
---|
| 198 | }
|
---|
| 199 |
|
---|
| 200 | ShellRotation
|
---|
| 201 | Integral::shell_rotation(int am, SymmetryOperation& so, int pure)
|
---|
| 202 | {
|
---|
| 203 | this->reference();
|
---|
| 204 | ShellRotation r(am, so, this, pure);
|
---|
| 205 | this->dereference();
|
---|
| 206 | return r;
|
---|
| 207 | }
|
---|
| 208 |
|
---|
| 209 | void
|
---|
| 210 | Integral::set_basis(const Ref<GaussianBasisSet> &b1,
|
---|
| 211 | const Ref<GaussianBasisSet> &b2,
|
---|
| 212 | const Ref<GaussianBasisSet> &b3,
|
---|
| 213 | const Ref<GaussianBasisSet> &b4)
|
---|
| 214 | {
|
---|
| 215 | bs1_ = b1;
|
---|
| 216 | bs2_ = b2;
|
---|
| 217 | bs3_ = b3;
|
---|
| 218 | bs4_ = b4;
|
---|
| 219 | if (bs2_.null()) bs2_ = bs1_;
|
---|
| 220 | if (bs3_.null()) bs3_ = bs2_;
|
---|
| 221 | if (bs4_.null()) bs4_ = bs3_;
|
---|
| 222 | }
|
---|
| 223 |
|
---|
| 224 | size_t
|
---|
| 225 | Integral::storage_required_eri(const Ref<GaussianBasisSet> &b1,
|
---|
| 226 | const Ref<GaussianBasisSet> &b2,
|
---|
| 227 | const Ref<GaussianBasisSet> &b3,
|
---|
| 228 | const Ref<GaussianBasisSet> &b4)
|
---|
| 229 | {
|
---|
| 230 | // By default, generated ERI evaluator will not need
|
---|
| 231 | // any significant amount of memory
|
---|
| 232 | return 0;
|
---|
| 233 | }
|
---|
| 234 |
|
---|
| 235 | size_t
|
---|
| 236 | Integral::storage_required_eri_deriv(const Ref<GaussianBasisSet> &b1,
|
---|
| 237 | const Ref<GaussianBasisSet> &b2,
|
---|
| 238 | const Ref<GaussianBasisSet> &b3,
|
---|
| 239 | const Ref<GaussianBasisSet> &b4)
|
---|
| 240 | {
|
---|
| 241 | // By default, generated derivative ERI evaluator will not need
|
---|
| 242 | // any significant amount of memory
|
---|
| 243 | return 0;
|
---|
| 244 | }
|
---|
| 245 |
|
---|
| 246 | size_t
|
---|
| 247 | Integral::storage_required_grt(const Ref<GaussianBasisSet> &b1,
|
---|
| 248 | const Ref<GaussianBasisSet> &b2,
|
---|
| 249 | const Ref<GaussianBasisSet> &b3,
|
---|
| 250 | const Ref<GaussianBasisSet> &b4)
|
---|
| 251 | {
|
---|
| 252 | // By default, generated GRT evaluator will not need
|
---|
| 253 | // any significant amount of memory
|
---|
| 254 | return 0;
|
---|
| 255 | }
|
---|
| 256 |
|
---|
| 257 | size_t
|
---|
| 258 | Integral::storage_unused()
|
---|
| 259 | {
|
---|
| 260 | ptrdiff_t tmp=storage_-storage_used_;
|
---|
| 261 | return (tmp<0?0:tmp);
|
---|
| 262 | }
|
---|
| 263 |
|
---|
| 264 | Ref<TwoBodyInt>
|
---|
| 265 | Integral::grt()
|
---|
| 266 | {
|
---|
| 267 | throw std::runtime_error("Integral::grt(): not implemented in this particular integrals factory.");
|
---|
| 268 | }
|
---|
| 269 |
|
---|
| 270 | Ref<OneBodyOneCenterInt>
|
---|
| 271 | Integral::point_charge1(const Ref<PointChargeData>&)
|
---|
| 272 | {
|
---|
| 273 | throw std::runtime_error("Integral::point_charge1(): not implemented in this particular integrals factory.");
|
---|
| 274 | }
|
---|
| 275 |
|
---|
| 276 | Ref<TwoBodyThreeCenterInt>
|
---|
| 277 | Integral::electron_repulsion3()
|
---|
| 278 | {
|
---|
| 279 | throw std::runtime_error("Integral::electron_repulsion3(): not implemented in this particular integrals factory.");
|
---|
| 280 | }
|
---|
| 281 |
|
---|
| 282 | Ref<TwoBodyThreeCenterDerivInt>
|
---|
| 283 | Integral::electron_repulsion3_deriv()
|
---|
| 284 | {
|
---|
| 285 | throw std::runtime_error("Integral::electron_repulsion3_deriv(): not implemented in this particular integrals factory.");
|
---|
| 286 | }
|
---|
| 287 |
|
---|
| 288 | Ref<TwoBodyTwoCenterInt>
|
---|
| 289 | Integral::electron_repulsion2()
|
---|
| 290 | {
|
---|
| 291 | throw std::runtime_error("Integral::electron_repulsion2(): not implemented in this particular integrals factory.");
|
---|
| 292 | }
|
---|
| 293 |
|
---|
| 294 | Ref<TwoBodyTwoCenterDerivInt>
|
---|
| 295 | Integral::electron_repulsion2_deriv()
|
---|
| 296 | {
|
---|
| 297 | throw std::runtime_error("Integral::electron_repulsion2_deriv(): not implemented in this particular integrals factory.");
|
---|
| 298 | }
|
---|
| 299 |
|
---|
| 300 | /////////////////////////////////////////////////////////////////////////////
|
---|
| 301 |
|
---|
| 302 | // Local Variables:
|
---|
| 303 | // mode: c++
|
---|
| 304 | // c-file-style: "ETS"
|
---|
| 305 | // End:
|
---|