1 /* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * HF/DFT from Hel 7 * 8 * Copyright © 2015 The Regents of the University of California 9 * All Rights Reserved 10 * 11 * Written by Susi Lehtola, Lawrence Berkeley National Laboratory 12 * 13 * This program is free software; you can redistribute it and/or 14 * modify it under the terms of the GNU General Public License 15 * as published by the Free Software Foundation; either version 2 16 * of the License, or (at your option) any later version. 17 */ 18 19 #include "erifit.h" 20 #include "basis.h" 21 #include "linalg.h" 22 #include "mathf.h" 23 #include "settings.h" 24 #include "eriworker.h" 25 26 extern Settings settings; 27 28 namespace ERIfit { 29 operator <(const bf_pair_t & lhs,const bf_pair_t & rhs)30 bool operator<(const bf_pair_t & lhs, const bf_pair_t & rhs) { 31 return lhs.idx < rhs.idx; 32 } 33 get_basis(BasisSet & basis,const BasisSetLibrary & blib,const ElementBasisSet & orbel)34 void get_basis(BasisSet & basis, const BasisSetLibrary & blib, const ElementBasisSet & orbel) { 35 // Settings needed to form basis set 36 Settings settings0(settings); 37 settings.add_scf_settings(); 38 settings.set_bool("BasisRotate", false); 39 settings.set_string("Decontract", ""); 40 settings.set_bool("UseLM", true); 41 42 // Atoms 43 std::vector<atom_t> atoms(1); 44 atoms[0].el=orbel.get_symbol(); 45 atoms[0].num=0; 46 atoms[0].x=atoms[0].y=atoms[0].z=0.0; 47 atoms[0].Q=0; 48 49 // Form basis set 50 construct_basis(basis,atoms,blib); 51 } 52 orthonormal_ERI_trans(const ElementBasisSet & orbel,double linthr,arma::mat & trans)53 void orthonormal_ERI_trans(const ElementBasisSet & orbel, double linthr, arma::mat & trans) { 54 // Form basis set library 55 BasisSetLibrary blib; 56 blib.add_element(orbel); 57 58 // Form basis set 59 BasisSet basis; 60 get_basis(basis,blib,orbel); 61 62 // Get orthonormal orbitals 63 arma::mat S(basis.overlap()); 64 arma::mat Sinvh(CanonicalOrth(S,linthr)); 65 66 // Sizes 67 size_t Nbf(Sinvh.n_rows); 68 size_t Nmo(Sinvh.n_cols); 69 70 // Fill matrix 71 trans.zeros(Nbf*Nbf,Nmo*Nmo); 72 printf("Size of orthogonal transformation matrix is %i x %i\n",(int) trans.n_rows,(int) trans.n_cols); 73 74 for(size_t iao=0;iao<Nbf;iao++) 75 for(size_t jao=0;jao<Nbf;jao++) 76 for(size_t imo=0;imo<Nmo;imo++) 77 for(size_t jmo=0;jmo<Nmo;jmo++) 78 trans(iao*Nbf+jao,imo*Nmo+jmo)=Sinvh(iao,imo)*Sinvh(jao,jmo); 79 } 80 compute_ERIs(const BasisSet & basis,arma::mat & eris)81 void compute_ERIs(const BasisSet & basis, arma::mat & eris) { 82 // Amount of functions 83 size_t Nbf(basis.get_Nbf()); 84 85 // Get shells in basis set 86 std::vector<GaussianShell> shells(basis.get_shells()); 87 // Get list of shell pairs 88 std::vector<shellpair_t> shpairs(basis.get_unique_shellpairs()); 89 90 // Print basis 91 // basis.print(true); 92 93 // Allocate memory for the integrals 94 eris.zeros(Nbf*Nbf,Nbf*Nbf); 95 printf("Size of integral matrix is %i x %i\n",(int) eris.n_rows,(int) eris.n_cols); 96 97 #ifdef _OPENMP 98 #pragma omp parallel 99 #endif 100 { 101 // Integral worker 102 ERIWorker *eri=new ERIWorker(basis.get_max_am(),basis.get_max_Ncontr()); 103 104 #ifdef _OPENMP 105 #pragma omp for schedule(dynamic,1) 106 #endif 107 for(size_t ip=0;ip<shpairs.size();ip++) 108 for(size_t jp=0;jp<=ip;jp++) { 109 // Shells are 110 size_t is=shpairs[ip].is; 111 size_t js=shpairs[ip].js; 112 size_t ks=shpairs[jp].is; 113 size_t ls=shpairs[jp].js; 114 115 // First functions on shells 116 size_t i0=shells[is].get_first_ind(); 117 size_t j0=shells[js].get_first_ind(); 118 size_t k0=shells[ks].get_first_ind(); 119 size_t l0=shells[ls].get_first_ind(); 120 121 // Amount of functions 122 size_t Ni=shells[is].get_Nbf(); 123 size_t Nj=shells[js].get_Nbf(); 124 size_t Nk=shells[ks].get_Nbf(); 125 size_t Nl=shells[ls].get_Nbf(); 126 127 // Compute integral block 128 eri->compute(&shells[is],&shells[js],&shells[ks],&shells[ls]); 129 // Get array 130 const std::vector<double> *erip=eri->getp(); 131 132 // Store integrals 133 for(size_t ii=0;ii<Ni;ii++) { 134 size_t i=i0+ii; 135 for(size_t jj=0;jj<Nj;jj++) { 136 size_t j=j0+jj; 137 for(size_t kk=0;kk<Nk;kk++) { 138 size_t k=k0+kk; 139 for(size_t ll=0;ll<Nl;ll++) { 140 size_t l=l0+ll; 141 142 // Go through the 8 permutation symmetries 143 double mel=(*erip)[((ii*Nj+jj)*Nk+kk)*Nl+ll]; 144 eris(i*Nbf+j,k*Nbf+l)=mel; 145 if(js!=is) 146 eris(j*Nbf+i,k*Nbf+l)=mel; 147 if(ks!=ls) 148 eris(i*Nbf+j,l*Nbf+k)=mel; 149 if(is!=js && ks!=ls) 150 eris(j*Nbf+i,l*Nbf+k)=mel; 151 152 if(ip!=jp) { 153 eris(k*Nbf+l,i*Nbf+j)=mel; 154 if(js!=is) 155 eris(k*Nbf+l,j*Nbf+i)=mel; 156 if(ks!=ls) 157 eris(l*Nbf+k,i*Nbf+j)=mel; 158 if(is!=js && ks!=ls) 159 eris(l*Nbf+k,j*Nbf+i)=mel; 160 } 161 } 162 } 163 } 164 } 165 } 166 167 // Free memory 168 delete eri; 169 } 170 } 171 compute_ERIs(const ElementBasisSet & orbel,arma::mat & eris)172 void compute_ERIs(const ElementBasisSet & orbel, arma::mat & eris) { 173 // Form basis set library 174 BasisSetLibrary blib; 175 blib.add_element(orbel); 176 177 // Form basis set 178 BasisSet basis; 179 get_basis(basis,blib,orbel); 180 181 // Calculate 182 compute_ERIs(basis,eris); 183 } 184 compute_diag_ERIs(const ElementBasisSet & orbel,arma::mat & eris)185 void compute_diag_ERIs(const ElementBasisSet & orbel, arma::mat & eris) { 186 // Form basis set library 187 BasisSetLibrary blib; 188 blib.add_element(orbel); 189 190 // Form basis set 191 BasisSet basis; 192 get_basis(basis,blib,orbel); 193 194 size_t Nbf(basis.get_Nbf()); 195 196 // Get shells in basis set 197 std::vector<GaussianShell> shells(basis.get_shells()); 198 // Get list of shell pairs 199 std::vector<shellpair_t> shpairs(basis.get_unique_shellpairs()); 200 201 // Print basis 202 // basis.print(true); 203 204 // Allocate memory for the integrals 205 eris.zeros(Nbf,Nbf); 206 printf("Size of integral matrix is %i x %i\n",(int) eris.n_rows,(int) eris.n_cols); 207 208 #ifdef _OPENMP 209 #pragma omp parallel 210 #endif 211 { 212 // Integral worker 213 ERIWorker *eri=new ERIWorker(basis.get_max_am(),basis.get_max_Ncontr()); 214 215 #ifdef _OPENMP 216 #pragma omp for schedule(dynamic,1) 217 #endif 218 for(size_t ip=0;ip<shpairs.size();ip++) { 219 // Shells are 220 size_t is=shpairs[ip].is; 221 size_t js=shpairs[ip].js; 222 223 // First functions on shells 224 size_t i0=shells[is].get_first_ind(); 225 size_t j0=shells[js].get_first_ind(); 226 227 // Amount of functions 228 size_t Ni=shells[is].get_Nbf(); 229 size_t Nj=shells[js].get_Nbf(); 230 231 // Compute integral block 232 eri->compute(&shells[is],&shells[js],&shells[is],&shells[js]); 233 // Get array 234 const std::vector<double> *erip=eri->getp(); 235 236 // Store integrals 237 for(size_t ii=0;ii<Ni;ii++) { 238 size_t i=i0+ii; 239 for(size_t jj=0;jj<Nj;jj++) { 240 size_t j=j0+jj; 241 eris(i,j)=(*erip)[((ii*Nj+jj)*Ni+ii)*Nj+jj]; 242 } 243 } 244 } 245 246 // Free memory 247 delete eri; 248 } 249 } 250 unique_exponent_pairs(const ElementBasisSet & orbel,int am1,int am2,std::vector<std::vector<shellpair_t>> & pairs,std::vector<double> & exps)251 void unique_exponent_pairs(const ElementBasisSet & orbel, int am1, int am2, std::vector< std::vector<shellpair_t> > & pairs, std::vector<double> & exps) { 252 // Form orbital basis set library 253 BasisSetLibrary orblib; 254 orblib.add_element(orbel); 255 256 // Form orbital basis set 257 BasisSet basis; 258 get_basis(basis,orblib,orbel); 259 260 // Get shells 261 std::vector<GaussianShell> shells(basis.get_shells()); 262 // and list of unique shell pairs 263 std::vector<shellpair_t> shpairs(basis.get_unique_shellpairs()); 264 265 // Create the exponent list 266 exps.clear(); 267 for(size_t ip=0;ip<shpairs.size();ip++) { 268 // Shells are 269 size_t is=shpairs[ip].is; 270 size_t js=shpairs[ip].js; 271 272 // Check am 273 if(!( (shells[is].get_am()==am1 && shells[js].get_am()==am2) || (shells[is].get_am()==am2 && shells[js].get_am()==am1))) 274 continue; 275 276 // Check that shells aren't contracted 277 if(shells[is].get_Ncontr()!=1 || shells[js].get_Ncontr()!=1) { 278 ERROR_INFO(); 279 throw std::runtime_error("Must use primitive basis set!\n"); 280 } 281 282 // Exponent value is 283 double zeta=shells[is].get_contr()[0].z + shells[js].get_contr()[0].z; 284 sorted_insertion<double>(exps,zeta); 285 } 286 287 // Create the pair list 288 pairs.resize(exps.size()); 289 for(size_t ip=0;ip<shpairs.size();ip++) { 290 // Shells are 291 size_t is=shpairs[ip].is; 292 size_t js=shpairs[ip].js; 293 294 // Check am 295 if(!( (shells[is].get_am()==am1 && shells[js].get_am()==am2) || (shells[is].get_am()==am2 && shells[js].get_am()==am1))) 296 continue; 297 298 // Pair is 299 double zeta=shells[is].get_contr()[0].z + shells[js].get_contr()[0].z; 300 size_t pos=sorted_insertion<double>(exps,zeta); 301 302 // Insert pair 303 pairs[pos].push_back(shpairs[ip]); 304 } 305 } 306 compute_cholesky_T(const ElementBasisSet & orbel,int am1,int am2,arma::mat & eris,arma::vec & exps_)307 void compute_cholesky_T(const ElementBasisSet & orbel, int am1, int am2, arma::mat & eris, arma::vec & exps_) { 308 // Form basis set library 309 BasisSetLibrary blib; 310 blib.add_element(orbel); 311 // Decontract the basis set 312 blib.decontract(); 313 314 // Form basis set 315 BasisSet basis; 316 get_basis(basis,blib,orbel); 317 318 // Get shells in basis sets 319 std::vector<GaussianShell> shells(basis.get_shells()); 320 321 // Get list of unique exponent pairs 322 std::vector< std::vector<shellpair_t> > upairs; 323 std::vector<double> exps; 324 unique_exponent_pairs(orbel,am1,am2,upairs,exps); 325 326 // Store exponents 327 exps_=arma::conv_to<arma::vec>::from(exps); 328 329 // Allocate memory for the integrals 330 eris.zeros(exps.size(),exps.size()); 331 #ifdef _OPENMP 332 #pragma omp parallel 333 #endif 334 { 335 // Integral worker 336 ERIWorker *eri=new ERIWorker(basis.get_max_am(),basis.get_max_Ncontr()); 337 338 #ifdef _OPENMP 339 #pragma omp for schedule(dynamic,1) 340 #endif 341 // Loop over unique exponent pairs 342 for(size_t iip=0;iip<upairs.size();iip++) 343 for(size_t jjp=0;jjp<=iip;jjp++) { 344 345 // Loop over individual shell pairs in the group 346 for(size_t ip=0;ip<upairs[iip].size();ip++) 347 for(size_t jp=0;jp<upairs[jjp].size();jp++) { 348 // Shells are 349 size_t is=upairs[iip][ip].is; 350 size_t js=upairs[iip][ip].js; 351 size_t ks=upairs[jjp][jp].is; 352 size_t ls=upairs[jjp][jp].js; 353 354 // Amount of functions 355 size_t Ni=shells[is].get_Nbf(); 356 size_t Nj=shells[js].get_Nbf(); 357 size_t Nk=shells[ks].get_Nbf(); 358 size_t Nl=shells[ls].get_Nbf(); 359 360 // Compute integral block 361 eri->compute(&shells[is],&shells[js],&shells[ks],&shells[ls]); 362 // Get array 363 const std::vector<double> *erip=eri->getp(); 364 365 // Store integrals 366 for(size_t ii=0;ii<Ni;ii++) 367 for(size_t jj=0;jj<Nj;jj++) 368 for(size_t kk=0;kk<Nk;kk++) 369 for(size_t ll=0;ll<Nl;ll++) { 370 double mel=std::abs((*erip)[((ii*Nj+jj)*Nk+kk)*Nl+ll]); 371 // mel*=mel; 372 373 eris(iip,jjp)+=mel; 374 if(iip!=jjp) 375 eris(jjp,iip)+=mel; 376 } 377 } 378 } 379 380 // Free memory 381 delete eri; 382 } 383 } 384 compute_fitint(const BasisSetLibrary & fitlib,const ElementBasisSet & orbel,arma::mat & fitint)385 void compute_fitint(const BasisSetLibrary & fitlib, const ElementBasisSet & orbel, arma::mat & fitint) { 386 // Form orbital basis set library 387 BasisSetLibrary orblib; 388 orblib.add_element(orbel); 389 390 // Form orbital basis set 391 BasisSet orbbas; 392 get_basis(orbbas,orblib,orbel); 393 394 // and fitting basis set 395 BasisSet fitbas; 396 get_basis(fitbas,fitlib,orbel); 397 398 // Coulomb normalize the fitting set 399 fitbas.coulomb_normalize(); 400 401 // Get shells in basis sets 402 std::vector<GaussianShell> orbsh(orbbas.get_shells()); 403 std::vector<GaussianShell> fitsh(fitbas.get_shells()); 404 // Get list of shell pairs 405 std::vector<shellpair_t> orbpairs(orbbas.get_unique_shellpairs()); 406 407 // Dummy shell 408 GaussianShell dummy(dummyshell()); 409 410 // Problem sizes 411 const size_t Norb(orbbas.get_Nbf()); 412 const size_t Nfit(fitbas.get_Nbf()); 413 414 // Allocate memory 415 fitint.zeros(Norb*Norb,Nfit); 416 417 #ifdef _OPENMP 418 #pragma omp parallel 419 #endif 420 { 421 // Integral worker 422 int maxam=std::max(orbbas.get_max_am(),fitbas.get_max_am()); 423 ERIWorker *eri=new ERIWorker(maxam,orbbas.get_max_Ncontr()); 424 425 #ifdef _OPENMP 426 #pragma omp for schedule(dynamic) 427 #endif 428 for(size_t ip=0;ip<orbpairs.size();ip++) 429 for(size_t as=0;as<fitsh.size();as++) { 430 // Orbital shells are 431 size_t is=orbpairs[ip].is; 432 size_t js=orbpairs[ip].js; 433 434 // First function is 435 size_t i0=orbsh[is].get_first_ind(); 436 size_t j0=orbsh[js].get_first_ind(); 437 size_t a0=fitsh[as].get_first_ind(); 438 439 // Amount of functions 440 size_t Ni=orbsh[is].get_Nbf(); 441 size_t Nj=orbsh[js].get_Nbf(); 442 size_t Na=fitsh[as].get_Nbf(); 443 444 // Compute integral block 445 eri->compute(&orbsh[is],&orbsh[js],&dummy,&fitsh[as]); 446 // Get array 447 const std::vector<double> *erip=eri->getp(); 448 449 // Store integrals 450 for(size_t ii=0;ii<Ni;ii++) { 451 size_t i=i0+ii; 452 for(size_t jj=0;jj<Nj;jj++) { 453 size_t j=j0+jj; 454 for(size_t aa=0;aa<Na;aa++) { 455 size_t a=a0+aa; 456 457 // Go through the two permutation symmetries 458 double mel=(*erip)[(ii*Nj+jj)*Na+aa]; 459 fitint(i*Norb+j,a)=mel; 460 fitint(j*Norb+i,a)=mel; 461 } 462 } 463 } 464 } 465 466 // Free memory 467 delete eri; 468 } 469 } 470 compute_ERIfit(const BasisSetLibrary & fitlib,const ElementBasisSet & orbel,double linthr,const arma::mat & fitint,arma::mat & fiteri)471 void compute_ERIfit(const BasisSetLibrary & fitlib, const ElementBasisSet & orbel, double linthr, const arma::mat & fitint, arma::mat & fiteri) { 472 // Form orbital basis set library 473 BasisSetLibrary orblib; 474 orblib.add_element(orbel); 475 476 // Form orbital basis set 477 BasisSet orbbas; 478 get_basis(orbbas,orblib,orbel); 479 480 // and fitting basis set 481 BasisSet fitbas; 482 get_basis(fitbas,fitlib,orbel); 483 // Coulomb normalize the fitting set 484 fitbas.coulomb_normalize(); 485 486 if(fitint.n_cols != fitbas.get_Nbf()) 487 throw std::runtime_error("Need to supply fitting integrals for ERIfit!\n"); 488 489 // Get shells in basis sets 490 std::vector<GaussianShell> orbsh(orbbas.get_shells()); 491 std::vector<GaussianShell> fitsh(fitbas.get_shells()); 492 // Get list of shell pairs 493 std::vector<shellpair_t> orbpairs(orbbas.get_unique_shellpairs()); 494 495 // Dummy shell 496 GaussianShell dummy(dummyshell()); 497 498 // Problem size 499 size_t Nfit(fitbas.get_Nbf()); 500 // Overlap matrix 501 arma::mat S(Nfit,Nfit); 502 503 #ifdef _OPENMP 504 #pragma omp parallel 505 #endif 506 { 507 // Integral worker 508 int maxam=std::max(orbbas.get_max_am(),fitbas.get_max_am()); 509 ERIWorker *eri=new ERIWorker(maxam,orbbas.get_max_Ncontr()); 510 511 // Compute the fitting basis overlap 512 #ifdef _OPENMP 513 #pragma omp for 514 #endif 515 for(size_t i=0;i<fitsh.size();i++) 516 for(size_t j=0;j<=i;j++) { 517 // Compute integral block 518 eri->compute(&fitsh[i],&dummy,&fitsh[j],&dummy); 519 // Get array 520 const std::vector<double> *erip=eri->getp(); 521 // Store integrals 522 size_t i0=fitsh[i].get_first_ind(); 523 size_t j0=fitsh[j].get_first_ind(); 524 size_t Ni=fitsh[i].get_Nbf(); 525 size_t Nj=fitsh[j].get_Nbf(); 526 for(size_t ii=0;ii<Ni;ii++) 527 for(size_t jj=0;jj<Nj;jj++) { 528 double mel=(*erip)[ii*Nj+jj]; 529 S(i0+ii,j0+jj)=mel; 530 S(j0+jj,i0+ii)=mel; 531 } 532 } 533 534 // Free memory 535 delete eri; 536 } 537 538 // Do the eigendecomposition 539 arma::vec Sval; 540 arma::mat Svec; 541 eig_sym_ordered(Sval,Svec,S); 542 543 // Count linearly independent vectors 544 size_t Nind=0; 545 for(size_t i=0;i<Sval.n_elem;i++) 546 if(Sval(i)>=linthr) 547 Nind++; 548 // and drop the linearly dependent ones 549 Sval=Sval.subvec(Sval.n_elem-Nind,Sval.n_elem-1); 550 Svec=Svec.cols(Svec.n_cols-Nind,Svec.n_cols-1); 551 552 // Form inverse overlap matrix 553 arma::mat S_inv; 554 S_inv.zeros(Svec.n_rows,Svec.n_rows); 555 for(size_t i=0;i<Sval.n_elem;i++) 556 S_inv+=Svec.col(i)*arma::trans(Svec.col(i))/Sval(i); 557 558 // Fitted ERIs are 559 fiteri=fitint*S_inv*arma::trans(fitint); 560 } 561 compute_diag_ERIfit(const BasisSetLibrary & fitlib,const ElementBasisSet & orbel,double linthr,const arma::mat & fitint,arma::mat & fiteri)562 void compute_diag_ERIfit(const BasisSetLibrary & fitlib, const ElementBasisSet & orbel, double linthr, const arma::mat & fitint, arma::mat & fiteri) { 563 // Form orbital basis set library 564 BasisSetLibrary orblib; 565 orblib.add_element(orbel); 566 567 // Form orbital basis set 568 BasisSet orbbas; 569 get_basis(orbbas,orblib,orbel); 570 571 // and fitting basis set 572 BasisSet fitbas; 573 get_basis(fitbas,fitlib,orbel); 574 // Coulomb normalize the fitting set 575 fitbas.coulomb_normalize(); 576 577 if(fitint.n_rows != orbbas.get_Nbf()*orbbas.get_Nbf()) 578 throw std::runtime_error("Need to supply fitting integrals for ERIfit!\n"); 579 if(fitint.n_cols != fitbas.get_Nbf()) 580 throw std::runtime_error("Need to supply fitting integrals for ERIfit!\n"); 581 582 // Get shells in basis sets 583 std::vector<GaussianShell> orbsh(orbbas.get_shells()); 584 std::vector<GaussianShell> fitsh(fitbas.get_shells()); 585 // Get list of shell pairs 586 std::vector<shellpair_t> orbpairs(orbbas.get_unique_shellpairs()); 587 588 // Dummy shell 589 GaussianShell dummy(dummyshell()); 590 591 // Problem size 592 size_t Nfit(fitbas.get_Nbf()); 593 // Overlap matrix 594 arma::mat S(Nfit,Nfit); 595 596 #ifdef _OPENMP 597 #pragma omp parallel 598 #endif 599 { 600 // Integral worker 601 int maxam=std::max(orbbas.get_max_am(),fitbas.get_max_am()); 602 ERIWorker *eri=new ERIWorker(maxam,orbbas.get_max_Ncontr()); 603 604 // Compute the fitting basis overlap 605 #ifdef _OPENMP 606 #pragma omp for 607 #endif 608 for(size_t i=0;i<fitsh.size();i++) 609 for(size_t j=0;j<=i;j++) { 610 // Compute integral block 611 eri->compute(&fitsh[i],&dummy,&fitsh[j],&dummy); 612 // Get array 613 const std::vector<double> *erip=eri->getp(); 614 // Store integrals 615 size_t i0=fitsh[i].get_first_ind(); 616 size_t j0=fitsh[j].get_first_ind(); 617 size_t Ni=fitsh[i].get_Nbf(); 618 size_t Nj=fitsh[j].get_Nbf(); 619 for(size_t ii=0;ii<Ni;ii++) 620 for(size_t jj=0;jj<Nj;jj++) { 621 double mel=(*erip)[ii*Nj+jj]; 622 S(i0+ii,j0+jj)=mel; 623 S(j0+jj,i0+ii)=mel; 624 } 625 } 626 627 // Free memory 628 delete eri; 629 } 630 631 // Do the eigendecomposition 632 arma::vec Sval; 633 arma::mat Svec; 634 eig_sym_ordered(Sval,Svec,S); 635 636 // Count linearly independent vectors 637 size_t Nind=0; 638 for(size_t i=0;i<Sval.n_elem;i++) 639 if(Sval(i)>=linthr) 640 Nind++; 641 // and drop the linearly dependent ones 642 Sval=Sval.subvec(Sval.n_elem-Nind,Sval.n_elem-1); 643 Svec=Svec.cols(Svec.n_cols-Nind,Svec.n_cols-1); 644 645 // Form inverse overlap matrix 646 arma::mat S_inv; 647 S_inv.zeros(Svec.n_rows,Svec.n_rows); 648 for(size_t i=0;i<Sval.n_elem;i++) 649 S_inv+=Svec.col(i)*arma::trans(Svec.col(i))/Sval(i); 650 651 // Fitted ERIs are 652 size_t Nbf(orbbas.get_Nbf()); 653 fiteri.zeros(Nbf,Nbf); 654 for(size_t i=0;i<Nbf;i++) 655 for(size_t j=0;j<=i;j++) { 656 double el=arma::as_scalar(fitint.row(i*Nbf+j)*S_inv*arma::trans(fitint.row(i*Nbf+j))); 657 fiteri(i,j)=el; 658 fiteri(j,i)=el; 659 } 660 } 661 } 662