| [0b990d] | 1 | //
|
|---|
| 2 | // int2e.cc
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (C) 2004 Sandia National Laboratories.
|
|---|
| 5 | //
|
|---|
| 6 | // Author: Joseph Kenny <jpkenny@sandia.gov>
|
|---|
| 7 | // Maintainer: Joseph Kenny
|
|---|
| 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 __GNUG__
|
|---|
| 29 | #pragma implementation
|
|---|
| 30 | #endif
|
|---|
| 31 |
|
|---|
| 32 | #include <chemistry/qc/intcca/int2e.h>
|
|---|
| 33 | #include <Chemistry_Chemistry_QC_GaussianBasis_DerivCenters.hh>
|
|---|
| 34 | #include <util/class/scexception.h>
|
|---|
| 35 |
|
|---|
| 36 | using namespace std;
|
|---|
| 37 | using namespace sc;
|
|---|
| 38 | using namespace MPQC;
|
|---|
| 39 | using namespace Chemistry;
|
|---|
| 40 | using namespace Chemistry::QC::GaussianBasis;
|
|---|
| 41 |
|
|---|
| 42 | Int2eCCA::Int2eCCA(Integral *integral,
|
|---|
| 43 | const Ref<GaussianBasisSet>&b1,
|
|---|
| 44 | const Ref<GaussianBasisSet>&b2,
|
|---|
| 45 | const Ref<GaussianBasisSet>&b3,
|
|---|
| 46 | const Ref<GaussianBasisSet>&b4,
|
|---|
| 47 | int order, size_t storage,
|
|---|
| 48 | IntegralEvaluatorFactory eval_factory,
|
|---|
| 49 | bool use_opaque, string eval_type ):
|
|---|
| 50 | bs1_(b1), bs2_(b2), bs3_(b3), bs4_(b4),
|
|---|
| 51 | erep_ptr_(0), integral_(integral), eval_factory_(eval_factory),
|
|---|
| 52 | use_opaque_(use_opaque), buffer_(0)
|
|---|
| 53 | {
|
|---|
| 54 |
|
|---|
| 55 | /* Allocate storage for the integral buffer. */
|
|---|
| 56 | int maxsize = bs1_->max_ncartesian_in_shell()
|
|---|
| 57 | *bs2_->max_ncartesian_in_shell()
|
|---|
| 58 | *bs3_->max_ncartesian_in_shell()
|
|---|
| 59 | *bs4_->max_ncartesian_in_shell();
|
|---|
| 60 | if( order == 1 )
|
|---|
| 61 | maxsize *= 9;
|
|---|
| 62 | else if( order != 0 )
|
|---|
| 63 | throw FeatureNotImplemented("only first order derivatives are available",
|
|---|
| 64 | __FILE__,__LINE__);
|
|---|
| 65 | if( !use_opaque ) buffer_ = new double[maxsize];
|
|---|
| 66 |
|
|---|
| 67 | cca_bs1_ = GaussianBasis_Molecular::_create();
|
|---|
| 68 | cca_bs2_ = GaussianBasis_Molecular::_create();
|
|---|
| 69 | cca_bs3_ = GaussianBasis_Molecular::_create();
|
|---|
| 70 | cca_bs4_ = GaussianBasis_Molecular::_create();
|
|---|
| 71 | cca_bs1_.initialize( bs1_.pointer(), bs1_->name() );
|
|---|
| 72 | cca_bs2_.initialize( bs2_.pointer(), bs2_->name() );
|
|---|
| 73 | cca_bs3_.initialize( bs3_.pointer(), bs3_->name() );
|
|---|
| 74 | cca_bs4_.initialize( bs4_.pointer(), bs4_->name() );
|
|---|
| 75 |
|
|---|
| 76 | cca_dc_ = Chemistry_QC_GaussianBasis_DerivCenters::_create();
|
|---|
| 77 |
|
|---|
| 78 | if( eval_type == "eri" ) {
|
|---|
| 79 | erep_ = eval_factory_.get_integral_evaluator4( "eri2", 0,
|
|---|
| 80 | cca_bs1_, cca_bs2_,
|
|---|
| 81 | cca_bs3_, cca_bs4_ );
|
|---|
| 82 | erep_ptr_ = &erep_;
|
|---|
| 83 | if( use_opaque_ )
|
|---|
| 84 | buffer_ = static_cast<double*>( erep_ptr_->get_buffer() );
|
|---|
| 85 | }
|
|---|
| 86 | else if( eval_type == "eri_1der") {
|
|---|
| 87 | erep_1der_ = eval_factory_.get_integral_evaluator4( "eri2", 1,
|
|---|
| 88 | cca_bs1_, cca_bs2_,
|
|---|
| 89 | cca_bs3_, cca_bs4_ );
|
|---|
| 90 | erep_1der_ptr_ = &erep_1der_;
|
|---|
| 91 | if( use_opaque_ )
|
|---|
| 92 | buffer_ = static_cast<double*>( erep_1der_ptr_->get_buffer() );
|
|---|
| 93 | }
|
|---|
| 94 | else {
|
|---|
| 95 | std::cerr << "integral type: " << eval_type << std::endl;
|
|---|
| 96 | throw InputError("unrecognized integral type",
|
|---|
| 97 | __FILE__,__LINE__);
|
|---|
| 98 | }
|
|---|
| 99 | if (!buffer_)
|
|---|
| 100 | throw ProgrammingError("buffer not assigned",
|
|---|
| 101 | __FILE__,__LINE__);
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 | void
|
|---|
| 105 | Int2eCCA::compute_erep( int is, int js, int ks, int ls )
|
|---|
| 106 | {
|
|---|
| 107 | cca_dc_.clear();
|
|---|
| 108 | if( use_opaque_ )
|
|---|
| 109 | erep_ptr_->compute( is, js, ks, ls, 0, cca_dc_ );
|
|---|
| 110 | else {
|
|---|
| 111 | sidl_buffer_ = erep_ptr_->compute_array( is, js, ks, ls, 0, cca_dc_ );
|
|---|
| 112 | int nelem = bs1_->shell(is).nfunction() * bs2_->shell(js).nfunction() *
|
|---|
| 113 | bs3_->shell(ks).nfunction() * bs4_->shell(ls).nfunction();
|
|---|
| 114 | copy_buffer(nelem);
|
|---|
| 115 | }
|
|---|
| 116 |
|
|---|
| 117 | if(!redundant_) {
|
|---|
| 118 | remove_redundant(is,js,ks,ls);
|
|---|
| 119 | }
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | void
|
|---|
| 123 | Int2eCCA::compute_erep_1der( int is, int js, int ks, int ls,
|
|---|
| 124 | Chemistry::QC::GaussianBasis::DerivCenters &dc )
|
|---|
| 125 | {
|
|---|
| 126 |
|
|---|
| 127 | if( use_opaque_ )
|
|---|
| 128 | erep_1der_ptr_->compute( is, js, ks, ls, 1, dc );
|
|---|
| 129 | else {
|
|---|
| 130 | sidl_buffer_ = erep_ptr_->compute_array( is, js, ks, ls, 1, dc );
|
|---|
| 131 | int nelem = bs1_->shell(is).nfunction() * bs2_->shell(js).nfunction() *
|
|---|
| 132 | bs3_->shell(ks).nfunction() * bs4_->shell(ls).nfunction();
|
|---|
| 133 | copy_buffer(nelem);
|
|---|
| 134 | }
|
|---|
| 135 |
|
|---|
| 136 | if(!redundant_) {
|
|---|
| 137 | remove_redundant(is,js,ks,ls);
|
|---|
| 138 | }
|
|---|
| 139 | }
|
|---|
| 140 |
|
|---|
| 141 | void
|
|---|
| 142 | Int2eCCA::copy_buffer( int n )
|
|---|
| 143 | {
|
|---|
| 144 | for( int i=0; i<n; ++i)
|
|---|
| 145 | buffer_[i] = sidl_buffer_.get(i);
|
|---|
| 146 | }
|
|---|
| 147 |
|
|---|
| 148 | #ifndef INTV3_ORDER
|
|---|
| 149 |
|
|---|
| 150 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 151 | // Code for removing redundant integrals
|
|---|
| 152 | // copied liberally from cints
|
|---|
| 153 |
|
|---|
| 154 | void get_nonredundant_ints_(double *source, double *target,
|
|---|
| 155 | int e13e24, int e12, int e34,
|
|---|
| 156 | GaussianShell* int_shell1,GaussianShell* int_shell2,
|
|---|
| 157 | GaussianShell* int_shell3,GaussianShell* int_shell4)
|
|---|
| 158 | {
|
|---|
| 159 |
|
|---|
| 160 | //std::cout << "\nremoving redundant integrals";
|
|---|
| 161 |
|
|---|
| 162 | int i,j,k,l;
|
|---|
| 163 |
|
|---|
| 164 | int nbf1 = int_shell1->nfunction();
|
|---|
| 165 | int nbf2 = int_shell2->nfunction();
|
|---|
| 166 | int nbf3 = int_shell3->nfunction();
|
|---|
| 167 | int nbf4 = int_shell4->nfunction();
|
|---|
| 168 |
|
|---|
| 169 | double *redundant_ptr = source;
|
|---|
| 170 | double *nonredundant_ptr = target;
|
|---|
| 171 |
|
|---|
| 172 | int nbf34 = nbf3*nbf4;
|
|---|
| 173 | for (i=0; i<nbf1; i++) {
|
|---|
| 174 | int jmax = e12 ? i : nbf2-1;
|
|---|
| 175 | for (j=0; j<=jmax; j++) {
|
|---|
| 176 | int kmax = e13e24 ? i : nbf3-1;
|
|---|
| 177 | for (k=0; k<=kmax; k++) {
|
|---|
| 178 | int lmax = e34 ? ( (e13e24&&(i==k)) ? j : k) :
|
|---|
| 179 | ( (e13e24&&(i==k)) ? j : nbf4-1);
|
|---|
| 180 | for (l=0; l<=lmax; l++) {
|
|---|
| 181 | *(nonredundant_ptr++) = redundant_ptr[l];
|
|---|
| 182 | }
|
|---|
| 183 | redundant_ptr += nbf4;
|
|---|
| 184 | }
|
|---|
| 185 | redundant_ptr += (nbf3-(kmax+1))*nbf4;
|
|---|
| 186 | }
|
|---|
| 187 | redundant_ptr += (nbf2-(jmax+1))*nbf34;
|
|---|
| 188 | }
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | void
|
|---|
| 192 | Int2eCCA::remove_redundant(int sh1, int sh2, int sh3, int sh4) {
|
|---|
| 193 |
|
|---|
| 194 | GaussianShell* int_shell1(&bs1_->shell(sh1));
|
|---|
| 195 | GaussianShell* int_shell2(&bs2_->shell(sh2));
|
|---|
| 196 | GaussianShell* int_shell3(&bs3_->shell(sh3));
|
|---|
| 197 | GaussianShell* int_shell4(&bs4_->shell(sh4));
|
|---|
| 198 |
|
|---|
| 199 | bool need_unique_ints_only = false;
|
|---|
| 200 | int e12,e34,e13e24;
|
|---|
| 201 | e12 = 0;
|
|---|
| 202 | if (int_shell1 == int_shell2 && int_shell1->nfunction()>1)
|
|---|
| 203 | e12 = 1;
|
|---|
| 204 | e34 = 0;
|
|---|
| 205 | if (int_shell3 == int_shell4 && int_shell3->nfunction()>1)
|
|---|
| 206 | e34 = 1;
|
|---|
| 207 | e13e24 = 0;
|
|---|
| 208 | if (int_shell1 == int_shell3 && int_shell2 ==
|
|---|
| 209 | int_shell4 && int_shell1->nfunction()*int_shell2->nfunction()>1)
|
|---|
| 210 | e13e24 = 1;
|
|---|
| 211 |
|
|---|
| 212 | if ( e12 || e34 || e13e24 )
|
|---|
| 213 | need_unique_ints_only = true;
|
|---|
| 214 |
|
|---|
| 215 | if (need_unique_ints_only) {
|
|---|
| 216 | std::cout.flush();
|
|---|
| 217 | get_nonredundant_ints_( buffer_, buffer_, e13e24, e12, e34,
|
|---|
| 218 | int_shell1, int_shell2, int_shell3, int_shell4 );
|
|---|
| 219 | }
|
|---|
| 220 | }
|
|---|
| 221 |
|
|---|
| 222 | #else
|
|---|
| 223 |
|
|---|
| 224 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 225 | // Code for removing redundant integrals
|
|---|
| 226 | // copied liberally from intV3
|
|---|
| 227 |
|
|---|
| 228 | static int
|
|---|
| 229 | shell_offset(Ref<GaussianBasisSet> cs, int off)
|
|---|
| 230 | {
|
|---|
| 231 | return off + cs->nshell();
|
|---|
| 232 | }
|
|---|
| 233 |
|
|---|
| 234 | #define INT_MAX1(n1) ((n1)-1)
|
|---|
| 235 | #define INT_MAX2(e12,i,n2) ((e12)?(i):((n2)-1))
|
|---|
| 236 | #define INT_MAX3(e13e24,i,n3) ((e13e24)?(i):((n3)-1))
|
|---|
| 237 | #define INT_MAX4(e13e24,e34,i,j,k,n4) \
|
|---|
| 238 | ((e34)?(((e13e24)&&((k)==(i)))?(j):(k)) \
|
|---|
| 239 | :((e13e24)&&((k)==(i)))?(j):(n4)-1)
|
|---|
| 240 |
|
|---|
| 241 | void
|
|---|
| 242 | nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
|
|---|
| 243 | int n1, int n2, int n3, int n4,
|
|---|
| 244 | int *red_off, int *nonred_off)
|
|---|
| 245 | {
|
|---|
| 246 | int nonredundant_index;
|
|---|
| 247 | int i,j,k,l;
|
|---|
| 248 |
|
|---|
| 249 | double *redundant_ptr = &buffer[*red_off];
|
|---|
| 250 | double *nonredundant_ptr = &buffer[*nonred_off];
|
|---|
| 251 |
|
|---|
| 252 | nonredundant_index = 0;
|
|---|
| 253 | int n34 = n3*n4;
|
|---|
| 254 | for (i=0; i<n1; i++) {
|
|---|
| 255 | int jmax = INT_MAX2(e12,i,n2);
|
|---|
| 256 | for (j=0; j<=jmax; j++) {
|
|---|
| 257 | int kmax = INT_MAX3(e13e24,i,n3);
|
|---|
| 258 | for (k=0; k<=kmax; k++) {
|
|---|
| 259 | int lmax = INT_MAX4(e13e24,e34,i,j,k,n4);
|
|---|
| 260 | for (l=0; l<=lmax; l++) {
|
|---|
| 261 | nonredundant_ptr[l] = redundant_ptr[l];
|
|---|
| 262 | }
|
|---|
| 263 | redundant_ptr += n4;
|
|---|
| 264 | nonredundant_index += lmax+1;
|
|---|
| 265 | nonredundant_ptr += lmax+1;
|
|---|
| 266 | }
|
|---|
| 267 | redundant_ptr += (n3-(kmax+1))*n4;
|
|---|
| 268 | }
|
|---|
| 269 | redundant_ptr += (n2-(jmax+1))*n34;
|
|---|
| 270 | }
|
|---|
| 271 | *red_off += n1*n2*n34;
|
|---|
| 272 | *nonred_off += nonredundant_index;
|
|---|
| 273 | }
|
|---|
| 274 |
|
|---|
| 275 | void
|
|---|
| 276 | Int2eCCA::remove_redundant(int is, int js, int ks, int ls) {
|
|---|
| 277 |
|
|---|
| 278 | int bs1_shell_offset = 0;
|
|---|
| 279 | int bs2_shell_offset, bs3_shell_offset, bs4_shell_offset;
|
|---|
| 280 |
|
|---|
| 281 | int shell_offset1 = shell_offset(bs1_,0);
|
|---|
| 282 | int shell_offset2, shell_offset3, shell_offset4;
|
|---|
| 283 | if (bs2_ == bs1_) {
|
|---|
| 284 | shell_offset2 = shell_offset1;
|
|---|
| 285 | bs2_shell_offset = bs1_shell_offset;
|
|---|
| 286 | }
|
|---|
| 287 | else {
|
|---|
| 288 | shell_offset2 = shell_offset(bs2_,shell_offset1);
|
|---|
| 289 | bs2_shell_offset = shell_offset1;
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | if (bs3_ == bs1_) {
|
|---|
| 293 | shell_offset3 = shell_offset2;
|
|---|
| 294 | bs3_shell_offset = bs1_shell_offset;
|
|---|
| 295 | }
|
|---|
| 296 | else if (bs3_ == bs2_) {
|
|---|
| 297 | shell_offset3 = shell_offset2;
|
|---|
| 298 | bs3_shell_offset = bs2_shell_offset;
|
|---|
| 299 | }
|
|---|
| 300 | else {
|
|---|
| 301 | shell_offset3 = shell_offset(bs3_,shell_offset2);
|
|---|
| 302 | bs3_shell_offset = shell_offset2;
|
|---|
| 303 | }
|
|---|
| 304 |
|
|---|
| 305 | if (bs4_ == bs1_) {
|
|---|
| 306 | bs4_shell_offset = bs1_shell_offset;
|
|---|
| 307 | }
|
|---|
| 308 | else if (bs4_ == bs2_) {
|
|---|
| 309 | bs4_shell_offset = bs2_shell_offset;
|
|---|
| 310 | }
|
|---|
| 311 | else if (bs4_ == bs3_) {
|
|---|
| 312 | bs4_shell_offset = bs3_shell_offset;
|
|---|
| 313 | }
|
|---|
| 314 | else {
|
|---|
| 315 | bs4_shell_offset = shell_offset3;
|
|---|
| 316 | }
|
|---|
| 317 |
|
|---|
| 318 | int int_unit2 = 0;
|
|---|
| 319 | int int_unit4 = 0;
|
|---|
| 320 | int osh1 = is + bs1_shell_offset;
|
|---|
| 321 | int osh2;
|
|---|
| 322 | if (!int_unit2) osh2 = js + bs2_shell_offset;
|
|---|
| 323 | int osh3 = ks + bs3_shell_offset;
|
|---|
| 324 | int osh4;
|
|---|
| 325 | if (!int_unit4) osh4 = ls + bs4_shell_offset;
|
|---|
| 326 |
|
|---|
| 327 | int sh1 = is;
|
|---|
| 328 | int sh2;
|
|---|
| 329 | if (!int_unit2) sh2 = js;
|
|---|
| 330 | int sh3 = ks;
|
|---|
| 331 | int sh4;
|
|---|
| 332 | if (!int_unit4) sh4 = ls;
|
|---|
| 333 |
|
|---|
| 334 | GaussianShell* int_unit_shell = 0;
|
|---|
| 335 | GaussianShell* int_shell1 = &bs1_->shell(sh1);
|
|---|
| 336 | GaussianShell* int_shell2;
|
|---|
| 337 | if (!int_unit2) int_shell2 = &bs2_->shell(sh2);
|
|---|
| 338 | else int_shell2 = int_unit_shell;
|
|---|
| 339 | GaussianShell* int_shell3 = &bs3_->shell(sh3);
|
|---|
| 340 | GaussianShell* int_shell4;
|
|---|
| 341 | if (!int_unit4) int_shell4 = &bs4_->shell(sh4);
|
|---|
| 342 | else int_shell4 = int_unit_shell;
|
|---|
| 343 |
|
|---|
| 344 | int redundant_offset = 0;
|
|---|
| 345 | int nonredundant_offset = 0;
|
|---|
| 346 |
|
|---|
| 347 | if ((osh1 == osh4)&&(osh2 == osh3)&&(osh1 != osh2)) {
|
|---|
| 348 | throw ProgrammingError("nonredundant integrals cannot be generated",
|
|---|
| 349 | __FILE__,__LINE__);
|
|---|
| 350 | }
|
|---|
| 351 |
|
|---|
| 352 | int e12 = (int_unit2?0:(osh1 == osh2));
|
|---|
| 353 | int e13e24 = ((osh1 == osh3)
|
|---|
| 354 | && ((int_unit2 && int_unit4)
|
|---|
| 355 | || ((int_unit2||int_unit4)?0:(osh2 == osh4))));
|
|---|
| 356 | int e34 = (int_unit4?0:(osh3 == osh4));
|
|---|
| 357 | if (e12||e34||e13e24) {
|
|---|
| 358 | nonredundant_erep(buffer_,e12,e34,e13e24,
|
|---|
| 359 | int_shell1->nfunction(),
|
|---|
| 360 | int_shell2->nfunction(),
|
|---|
| 361 | int_shell3->nfunction(),
|
|---|
| 362 | int_shell4->nfunction(),
|
|---|
| 363 | &redundant_offset,
|
|---|
| 364 | &nonredundant_offset);
|
|---|
| 365 | }
|
|---|
| 366 | }
|
|---|
| 367 |
|
|---|
| 368 | #endif
|
|---|
| 369 |
|
|---|
| 370 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 371 |
|
|---|
| 372 | // Local Variables:
|
|---|
| 373 | // mode: c++
|
|---|
| 374 | // End:
|
|---|
