1/* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * DFT from Hel 7 * 8 * Written by Susi Lehtola, 2010-2013 9 * Copyright (c) 2010-2013, Susi Lehtola 10 * 11 * This program is free software; you can redistribute it and/or 12 * modify it under the terms of the GNU General Public License 13 * as published by the Free Software Foundation; either version 2 14 * of the License, or (at your option) any later version. 15 */ 16 17#if defined(RESTRICTED) && defined(DFT) 18void SCF::RDFT(rscf_t & sol, const std::vector<double> & occs, double convthr, const dft_t dft0) 19 20#elif defined(RESTRICTED) && defined(HF) 21void SCF::RHF(rscf_t & sol, const std::vector<double> & occs, double convthr) 22 23#elif defined(UNRESTRICTED) && defined(DFT) 24void SCF::UDFT(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, double convthr, const dft_t dft0) 25 26#elif defined(UNRESTRICTED) && defined(HF) 27void SCF::UHF(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, double convthr) 28 29#elif defined(UNRESTRICTED) && defined(_ROHF) 30void SCF::ROHF(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, double convthr) 31 32#elif defined(DFT) && defined(FULLHOLE) 33size_t XRSSCF::full_hole(uscf_t & sol, double convthr, dft_t dft0, bool xch) 34 35#elif defined(DFT) && defined(HALFHOLE) 36size_t XRSSCF::half_hole(uscf_t & sol, double convthr, dft_t dft0) 37#endif 38{ 39 40 /// Formation of Fock matrix 41#ifdef RESTRICTED 42#if defined(HF) 43#define form_fock(sol) Fock_RHF(sol,occs); nfock++; 44#elif defined(DFT) 45#define form_fock(sol) Fock_RDFT(sol,occs,dft,grid,nlgrid); nfock++; 46#endif 47 48#else 49 // Unrestricted case: first core excited states 50 51#if defined(HALFHOLE) 52#define form_fock(sol) Fock_half_hole(sol,dft,occa,occb,grid,nlgrid); nfock++; 53#elif defined(FULLHOLE) 54#define form_fock(sol) Fock_full_hole(sol,dft,occa,occb,grid,nlgrid,xch); nfock++; 55 56 // Conventional HF and DFT 57#elif defined(HF) 58#define form_fock(sol) Fock_UHF(sol,occa,occb); nfock++; 59#elif defined(_ROHF) 60#define form_fock(sol) Fock_ROHF(sol,occa,occb); nfock++; 61#elif defined(DFT) 62#define form_fock(sol) Fock_UDFT(sol,occa,occb,dft,grid,nlgrid); nfock++; 63#endif 64#endif 65 66 /// Find excited orbital 67#if defined(FULLHOLE) || defined(HALFHOLE) 68#define find_exc(sol) if(spin) {ixc_orb=find_excited_orb(*basisp,coreorb,sol.Cb,noccb);} else {ixc_orb=find_excited_orb(*basisp,coreorb,sol.Ca,nocca);} 69#endif 70 71 /// Occupancy update 72#if defined(FULLHOLE) 73#define upd_occa() if(spin) {occa=norm_occ(nocca);} else { if(xch) {occa=xch_occ(ixc_orb,nocca);} else {occa=fch_occ(ixc_orb,nocca);} } 74#define upd_occb() if(!spin) {occb=norm_occ(noccb);} else { if(xch) {occb=xch_occ(ixc_orb,noccb);} else {occb=fch_occ(ixc_orb,noccb);} } 75#elif defined(HALFHOLE) 76#define upd_occa() if(spin) {occa=norm_occ(nocca);} else { occa=tp_occ(ixc_orb,nocca);} 77#define upd_occb() if(!spin) {occb=norm_occ(noccb);} else { occb=tp_occ(ixc_orb,noccb);} 78#endif 79 80 /// Density formation 81#ifdef RESTRICTED 82#define form_dens(sol) form_density(sol,occs); 83#else 84#define form_dens_ur(sol) form_density(sol,occa,occb); 85 86#if defined(HALFHOLE) || defined(FULLHOLE) 87#define form_dens(sol) find_exc(sol); upd_occa(); upd_occb(); form_dens_ur(sol); 88#else 89#define form_dens(sol) form_dens_ur(sol) 90#endif 91 92#endif // ifdef RESTRICTED 93 94 95#if ( defined(HALFHOLE) || defined(FULLHOLE) ) 96 // Occupation vector of spin up and spin down 97 std::vector<double> occa; 98 std::vector<double> occb; 99#endif 100 101 // Determine number of occupied states 102#if !defined(HALFHOLE) && !defined(FULLHOLE) 103 104#if defined(RESTRICTED) 105 size_t nocc; 106 for(nocc=occs.size()-1;nocc<occs.size();nocc--) 107 if(occs[nocc]>0) 108 break; 109 nocc++; 110#else 111 size_t nocca; 112 for(nocca=occa.size()-1;nocca<occa.size();nocca--) 113 if(occa[nocca]>0) 114 break; 115 nocca++; 116 117 size_t noccb; 118 for(noccb=occb.size()-1;noccb<occb.size();noccb--) 119 if(occb[noccb]>0) 120 break; 121 noccb++; 122#endif 123#endif 124 125 int nfock=0; 126 127 Timer t; 128 Timer ttot; 129 130 // DIIS error 131 double diiserr=DBL_MAX; 132 // Was DIIS succesful? 133 bool diissucc=false; 134 135 // Helper arrays 136 arma::mat orbs; 137 arma::mat Horth; 138 139#ifdef RESTRICTED 140 // DIIS iterator 141 rDIIS diis(S,Sinvh,usediis,diiseps,diisthr,useadiis,verbose,diisorder); 142 // Broyden 143 Broyden broyd(verbose); 144 145 // Coulomb and exchange matrices 146 arma::mat J(Nbf,Nbf), K(Nbf,Nbf); 147 J.zeros(); 148 K.zeros(); 149#else 150 // DIIS iterator 151 uDIIS diis(S,Sinvh,diiscomb,usediis,diiseps,diisthr,useadiis,verbose,diisorder); 152 153 // Broyden 154 Broyden broyd_sum(verbose); 155 Broyden broyd_diff(verbose); 156 157#endif 158 159 // Dipole moment 160 arma::vec dipmom; 161 162 // Change in energy from last iteration 163 double deltaE=0; 164 165 // Maximum and RMS differences of AO density matrix 166 double rmsdiff=0.0, maxdiff=0.0; 167#ifndef RESTRICTED 168 double rmsdiffa=0.0, maxdiffa=0.0; 169 double rmsdiffb=0.0, maxdiffb=0.0; 170#endif 171 172#ifdef RESTRICTED 173 if(sol.P.n_rows!=Nbf) 174#else 175 if(sol.Pa.n_rows!=Nbf || sol.Pb.n_rows!=Nbf) 176#endif 177 { 178 throw std::runtime_error("No starting guess provided for SCF solver!\n"); 179 } 180 181#if defined(FULLHOLE) || defined(HALFHOLE) 182 size_t ixc_orb=0; 183 // X-ray calculations need occupations for exact exchange 184 form_dens(sol); 185#endif 186 187 // Print orbital energies 188 if(verbose) { 189#ifdef RESTRICTED 190 if(sol.E.n_elem) 191 print_E(sol.E,occs,false); 192#else 193 if(sol.Ea.n_elem) { 194 printf("alpha: "); 195 print_E(sol.Ea,occa,false); 196 printf("beta: "); 197 print_E(sol.Eb,occb,false); 198 } 199#endif 200 fflush(stdout); 201 } 202 203#ifdef DFT 204 // Run with non-local correlation? 205 dft_t dft(dft0); 206 207 // Range separation constants 208 double omega, kfull, kshort; 209 range_separation(dft.x_func,omega,kfull,kshort); 210 211 if(verbose) { 212 if(omega!=0.0) { 213 printf("\nUsing range separated exchange with range separation constant omega = % .3f.\n",omega); 214 printf("Using % .3f %% short range and % .3f %% long range exchange.\n",(kfull+kshort)*100,kfull*100); 215 } else if(kfull!=0.0) 216 printf("\nUsing hybrid exchange with % .3f %% of exact exchange.\n",kfull*100); 217 else 218 printf("\nA pure exchange functional used, no exact exchange.\n"); 219 } 220 221 // Evaluators for range separated part 222 if(omega!=0.0) 223 fill_rs(omega); 224 225 DFTGrid grid(basisp,verbose,dft.lobatto); 226 DFTGrid nlgrid(basisp,verbose,dft.lobatto); 227#endif 228 229 230 // Sparse output to stderr 231 if(verbose && maxiter>0) { 232 fprintf(stderr,"Running "); 233#if defined(RESTRICTED) 234 fprintf(stderr,"restricted "); 235#else 236# if defined(_ROHF) 237 fprintf(stderr,"restricted open-shell "); 238#else 239#if defined(UNRESTRICTED) 240 fprintf(stderr,"unrestricted "); 241#endif 242#endif 243#endif 244 245#ifndef DFT 246 fprintf(stderr,"HF "); 247#else 248 if(dft.c_func>0) { 249 // Correlation exists. 250 fprintf(stderr,"%s-%s ",get_keyword(dft.x_func).c_str(),get_keyword(dft.c_func).c_str()); 251 } else 252 fprintf(stderr,"%s ",get_keyword(dft.x_func).c_str()); 253#endif 254 255 fprintf(stderr,"calculation"); 256 if(densityfit) { 257 fprintf(stderr," with density fitting"); 258 } 259 260#ifdef DFT 261#if (defined(FULLHOLE) || defined(HALFHOLE)) 262 if(densityfit) 263 fprintf(stderr," and "); 264 else 265 fprintf(stderr," with "); 266#ifdef FULLHOLE 267 fprintf(stderr,"full core hole"); 268 if(xch) 269 fprintf(stderr," (XCH)"); 270#else // half hole 271 fprintf(stderr,"half core hole"); 272#endif 273#endif 274 275#endif // end DFT clause 276 fprintf(stderr,".\n%4s %16s %10s %9s %9s %9s %10s\n","iter","E","dE","RMS dens","MAX dens","DIIS err","titer (s)"); 277 } 278 fflush(stdout); 279 280#ifdef DFT 281 if(dft.x_func>0 || dft.c_func>0) { 282 if(dft.adaptive) { 283#ifdef RESTRICTED 284 // Form DFT quadrature grid 285 grid.construct(sol.P,dft.gridtol,dft.x_func,dft.c_func); 286#else 287 // Form DFT quadrature grid 288 grid.construct(sol.Pa,sol.Pb,dft.gridtol,dft.x_func,dft.c_func); 289#endif 290 } else { 291 // Fixed size grid 292 grid.construct(dft.nrad,dft.lmax,dft.x_func,dft.c_func,strictint); 293 // Nonlocal grid? 294 if(dft.nl) 295 nlgrid.construct(dft.nlnrad,dft.nllmax,true,false,false,strictint,true); 296 } 297 298 if(verbose) { 299 fflush(stdout); 300 fprintf(stderr,"%-65s %10.3f\n"," DFT grid formation",t.get()); 301 } 302 } 303#endif // DFT 304 305 // Has calculation been converged? 306 bool converged=false; 307 308 // Solution of last iteration 309#ifdef RESTRICTED 310 rscf_t oldsol; 311#else 312 uscf_t oldsol; 313#endif 314 oldsol.en.E=0.0; 315 316 // Linear symmetry? 317 if(lincalc && linfreeze) { 318#ifdef RESTRICTED 319 arma::imat occ(basisp->count_m_occupied(sol.C.cols(0,nocc-1))); 320#else 321 arma::imat occ(basisp->count_m_occupied(sol.Ca.cols(0,nocca-1),sol.Cb.cols(0,noccb-1))); 322#endif 323 occ.save(linoccfname,arma::raw_ascii); 324 if(verbose) 325 printf("Guess occupations saved in %s\n",linoccfname.c_str()); 326 } 327 328 // Pad occupancies with zeros (needed e.g. in Casida routines) 329#if !defined(RESTRICTED) 330 std::vector<double> occar(occa), occbr(occb); 331 while(occar.size()<Sinvh.n_cols) 332 occar.push_back(0.0); 333 while(occbr.size()<Sinvh.n_cols) 334 occbr.push_back(0.0); 335#else 336 std::vector<double> occr(occs); 337 while(occr.size()<Sinvh.n_cols) 338 occr.push_back(0.0); 339#endif 340 341 // Write current matrices to checkpoint. 342 // To use the file effectively, we keep it open for the whole shebang. 343 if(maxiter>0) { 344 chkptp->open(); 345 chkptp->write("tol",intthr); 346 chkptp->write("P",sol.P); 347 chkptp->write(sol.en); 348#if !defined(RESTRICTED) 349 chkptp->write("Ca",sol.Ca); 350 chkptp->write("Cb",sol.Cb); 351 352 chkptp->write("Ea",sol.Ea); 353 chkptp->write("Eb",sol.Eb); 354 355 chkptp->write("Pa",sol.Pa); 356 chkptp->write("Pb",sol.Pb); 357 358 chkptp->write("occa",occar); 359 chkptp->write("occb",occbr); 360 361 // Restricted 362 chkptp->write("Restricted",0); 363#else 364 chkptp->write("C",sol.C); 365 chkptp->write("E",sol.E); 366 chkptp->write("occs",occr); 367 // Unrestricted 368 chkptp->write("Restricted",1); 369#endif 370 chkptp->write("Converged",converged); 371 chkptp->close(); 372 373 } else { 374 Timer tg; 375 // Form the Fock operator. 376 form_fock(sol); 377 // Solve FC=ESC 378 if(verbose) { 379 if(shift==0.0) 380 printf("Solving SCF equations ... "); 381 else 382 printf("Solving SCF equations with level shift %.3f ... ",shift); 383 fflush(stdout); 384 t.set(); 385 } 386 387 // Calculate dipole moment 388 dipmom=dipole_moment(sol.P,*basisp); 389 // Do the diagonalization 390 diagonalize(sol,shift); 391 392 if(verbose) 393 printf("done (%s)\n",t.elapsed().c_str()); 394 } 395 396 // Loop: 397 int iiter=1; 398 399#ifdef DFT 400 // Start out without non-local correlation since it's usually costly 401 // and of minor importance to the density 402 if(dft.nl) 403 dft.nl=false; 404#endif 405 406 while(iiter<=maxiter) { 407 Timer titer; 408 409 if(verbose) 410 printf("\n ******* Iteration %4i ********\n",iiter); 411 412 // Form the Fock operator. 413 form_fock(sol); 414 415 // Compute change of energy 416 deltaE=sol.en.E-oldsol.en.E; 417 418#ifdef RESTRICTED 419 // Update DIIS stack of matrices 420 diis.update(sol.H,sol.P,sol.en.E,diiserr); 421#else 422 // Update DIIS stacks of matrices 423 diis.update(sol.Ha,sol.Hb,sol.Pa,sol.Pb,sol.en.E,diiserr); 424#endif 425 426 if(iiter>1 && usebroyden) { 427#ifdef RESTRICTED 428 // Update Broyden mixer 429 broyd.push_x(MatToVec(oldsol.H)); 430 broyd.push_f(MatToVec(oldsol.H-sol.H)); 431#else 432 // Compute sum and difference 433 arma::mat Hs=sol.Ha+sol.Hb; 434 arma::mat Hd=sol.Ha-sol.Hb; 435 436 arma::mat Hsold=oldsol.Ha+oldsol.Hb; 437 arma::mat Hdold=oldsol.Ha-oldsol.Hb; 438 439 // Update Broyden mixers 440 broyd_sum.push_x(MatToVec(Hsold)); 441 broyd_sum.push_f(MatToVec(Hsold-Hs)); 442 443 broyd_diff.push_x(MatToVec(Hdold)); 444 broyd_diff.push_f(MatToVec(Hdold-Hd)); 445#endif 446 } 447 448 // Solve DIIS 449 try { 450#ifdef RESTRICTED 451 // Solve new matrix 452 diis.solve_F(sol.H); 453#else 454 // Solve new matrices 455 diis.solve_F(sol.Ha,sol.Hb); 456#endif 457 diissucc=true; 458 } catch(std::runtime_error &) { 459 diissucc=false; 460 } 461 462 // Perform Broyden interpolation 463 if(usebroyden && !diissucc && iiter>1) { 464 465 if(verbose) { 466 printf("Performing Broyden interpolation of Fock operator ... "); 467 fflush(stdout); 468 t.set(); 469 } 470 471#ifdef RESTRICTED 472 // Update Hamiltonian 473 sol.H=VecToMat(broyd.update_x(),Nbf,Nbf); 474#else 475 arma::mat Hs=VecToMat(broyd_sum.update_x(),Nbf,Nbf); 476 arma::mat Hd=VecToMat(broyd_diff.update_x(),Nbf,Nbf); 477 478 // Update Hamiltonians 479 sol.Ha=0.5*(Hs+Hd); 480 sol.Hb=0.5*(Hs-Hd); 481#endif 482 483 if(verbose) 484 printf("done (%s)\n",t.elapsed().c_str()); 485 } 486 487 // Save old solution 488 oldsol=sol; 489 490 if(usetrrh) { 491 // Solve FC=ESC 492 if(verbose) { 493 printf("\nSolving TRRH equations.\n"); 494 fflush(stdout); 495 t.set(); 496 } 497 498#ifdef RESTRICTED 499 arma::mat Cnew; 500 arma::vec Enew; 501 TRRH_update(sol.H,sol.C,S,Cnew,Enew,nocc,verbose,trrhmins); 502 503 // Update solution 504 sol.C=Cnew; 505 sol.E=Enew; 506 507 // Check orthonormality of orbitals 508 check_orth(sol.C,S,false); 509#else 510 arma::mat Canew; 511 arma::vec Eanew; 512 TRRH_update(sol.Ha,sol.Ca,S,Canew,Eanew,nocca,verbose,trrhmins); 513 514 arma::mat Cbnew; 515 arma::vec Ebnew; 516 TRRH_update(sol.Hb,sol.Cb,S,Cbnew,Ebnew,noccb,verbose,trrhmins); 517 518 // Update solutions 519 sol.Ca=Canew; 520 sol.Cb=Cbnew; 521 522 sol.Ea=Eanew; 523 sol.Eb=Ebnew; 524 525 // Check orthonormality of orbitals 526 check_orth(sol.Ca,S,false); 527 check_orth(sol.Cb,S,false); 528#endif 529 530 if(verbose) 531 printf("TRRH solved in %s.\n\n",t.elapsed().c_str()); 532 533 } else { 534 // Solve FC=ESC 535 if(verbose) { 536 if(shift==0.0) 537 printf("\nSolving SCF equations ... "); 538 else 539 printf("\nSolving SCF equations with level shift %.3f ... ",shift); 540 fflush(stdout); 541 t.set(); 542 } 543 544 // Do the diagonalization 545 diagonalize(sol,shift); 546 547 if(verbose) 548 printf("done (%s)\n",t.elapsed().c_str()); 549 } 550 551#if !defined(FULLHOLE) && !defined(HALFHOLE) 552 if(lincalc && iiter<=readlinocc) { 553 arma::mat linoccs; 554 linoccs.load(linoccfname,arma::raw_ascii); 555 if(linoccs.n_cols < 3) { 556 throw std::logic_error("Must have at least three columns in occupation data.\n"); 557 } 558 559 arma::vec occnuma, occnumb; 560 // Number of occupied alpha orbitals is first column 561 occnuma=linoccs.col(0); 562 // Number of occupied beta orbitals is second column 563 occnumb=linoccs.col(1); 564 // m value is third column 565 std::vector<arma::uvec> occsym(linoccs.n_rows); 566 for(size_t i=0;i<linoccs.n_rows;i++) { 567 int m=round(linoccs(i,2)); 568 occsym[i]=basisp->m_indices(m); 569 } 570 571 // Check length of specifications 572#ifdef RESTRICTED 573 double numab=0.0; 574 for(size_t i=0;i<occnuma.size();i++) { 575 numab+=occnuma(i); 576 if(occnuma(i) != occnumb(i)) { 577 throw std::logic_error("Alpha occupations don't match beta occupations even though calculation is spin-restricted!\n"); 578 } 579 } 580 if(std::abs(numab-nocc)>10*DBL_EPSILON) { 581 printf("Warning - occupations differ by %f electrons from expected!\n",numab-nocc); 582 } 583#else 584 double numa=0.0, numb=0.0; 585 for(size_t i=0;i<occnuma.size();i++) 586 numa+=occnuma(i); 587 for(size_t i=0;i<occnumb.size();i++) 588 numb+=occnumb(i); 589 if(std::abs(numa-nocca)>10*DBL_EPSILON) { 590 printf("Warning - alpha occupation differs by %f electrons from expected!\n",numa-nocca); 591 } 592 if(std::abs(numb-noccb)>10*DBL_EPSILON) { 593 printf("Warning - beta occupation differs by %f electrons from expected!\n",numb-noccb); 594 } 595#endif 596 597 // Make sure wanted orbitals are occupied 598#ifdef RESTRICTED 599 sol.P=2*enforce_occupations(sol.C,sol.E,S,occnuma,occsym); 600#else 601 sol.Pa=enforce_occupations(sol.Ca,sol.Ea,S,occnuma,occsym); 602 sol.Pb=enforce_occupations(sol.Cb,sol.Eb,S,occnumb,occsym); 603 sol.P=sol.Pa+sol.Pb; 604#endif 605 } else 606#endif 607 { 608 // Form new density matrix 609 form_dens(sol); 610 } 611 612 // Change-of-density matrices 613 arma::mat deltaP=sol.P-oldsol.P; 614#ifndef RESTRICTED 615 arma::mat deltaPa=sol.Pa-oldsol.Pa; 616 arma::mat deltaPb=sol.Pb-oldsol.Pb; 617#endif 618 619 // Compute dipole moment 620 dipmom=dipole_moment(sol.P,*basisp); 621 622 // Compute convergence criteria 623 maxdiff=max_abs(deltaP/2.0); 624 rmsdiff=rms_norm(deltaP/2.0); 625 626 // Convergence checked against 627 double maxdiff_cvd(maxdiff); 628 double rmsdiff_cvd(maxdiff); 629 630#ifndef RESTRICTED 631 maxdiffa=max_abs(deltaPa); 632 maxdiffb=max_abs(deltaPb); 633 634 rmsdiffa=rms_norm(deltaPa); 635 rmsdiffb=rms_norm(deltaPb); 636 637 maxdiff_cvd=std::max(maxdiffa,maxdiffb); 638 rmsdiff_cvd=std::max(rmsdiffa,rmsdiffb); 639#endif 640 641 // Print out status information 642 if(verbose) { 643 printf("\n"); 644 printf("%-30s: % .16e\n","Total energy",sol.en.E); 645 printf("%-30s: % e\n","DIIS error",diiserr); 646 printf("%-30s: % e\n","Energy change",deltaE); 647 printf("%-30s: % e\n","Max total density change",maxdiff); 648 printf("%-30s: % e\n","Max rms density change",rmsdiff); 649#ifndef RESTRICTED 650 printf("%-30s: % e\n","Max total alpha density change",maxdiffa); 651 printf("%-30s: % e\n","Max rms alpha density change",rmsdiffa); 652 printf("%-30s: % e\n","Max total beta density change",maxdiffb); 653 printf("%-30s: % e\n","Max rms beta density change",rmsdiffb); 654#endif 655 printf("Dipole mu = (% 08.8f, % 08.8f, % 08.8f) D\n",dipmom(0)/AUINDEBYE,dipmom(1)/AUINDEBYE,dipmom(2)/AUINDEBYE); 656 657 printf("\nIteration took %s.\n",titer.elapsed().c_str()); 658 fflush(stdout); 659 } 660 661 // Sparse output 662 if(verbose) 663 fprintf(stderr,"%4i % 16.8f % 10.3e %9.3e %9.3e %9.3e %10.3f\n",iiter,sol.en.E,deltaE,rmsdiff_cvd,maxdiff_cvd,diiserr,titer.get()); 664 665#ifdef DFT 666 if(dft0.nl && !dft.nl && diiserr<=1e-3) { 667 if(verbose) { 668 printf("Turning on non-local correlation contributions\n"); 669 fprintf(stderr,"Turning on non-local correlation contributions\n"); 670 } 671 dft.nl=true; 672 diis.clear(); 673 continue; 674 } 675 else 676#endif 677 if(diiserr < convthr) { 678 converged=true; 679 if(verbose) 680 printf("\n ******* Convergence achieved ********\n"); 681 } 682 683 // Write current matrices to checkpoint. 684 // To use the file effectively, we keep it open for the whole shebang. 685 chkptp->open(); 686 chkptp->write("P",sol.P); 687 chkptp->write(sol.en); 688#if !defined(RESTRICTED) 689 chkptp->write("Ca",sol.Ca); 690 chkptp->write("Cb",sol.Cb); 691 692 chkptp->write("Ea",sol.Ea); 693 chkptp->write("Eb",sol.Eb); 694 695 chkptp->write("Ha",sol.Ha); 696 chkptp->write("Hb",sol.Hb); 697 698 chkptp->write("Pa",sol.Pa); 699 chkptp->write("Pb",sol.Pb); 700 701 chkptp->write("occa",occar); 702 chkptp->write("occb",occbr); 703 704 // Restricted 705 chkptp->write("Restricted",0); 706#else 707 chkptp->write("C",sol.C); 708 chkptp->write("E",sol.E); 709 chkptp->write("H",sol.H); 710 chkptp->write("occs",occr); 711 // Unrestricted 712 chkptp->write("Restricted",1); 713#endif 714 chkptp->write("Converged",converged); 715 chkptp->close(); 716 717 // Check convergence 718 if(converged) 719 break; 720 721 iiter++; 722 } // End SCF cycle 723 724 if(verbose) { 725 726 if(converged) { 727 std::string method= 728#ifdef RESTRICTED 729 "R" 730#elif defined(_ROHF) 731 "RO" 732#else 733 "U" 734#endif 735#ifdef DFT 736 "DFT" 737#else 738 "HF" 739#endif 740 ; 741 742 printf("Solution of %s took %s.\n",method.c_str(),ttot.elapsed().c_str()); 743 fprintf(stderr,"%s converged in %s.\n",method.c_str(),ttot.elapsed().c_str()); 744 } 745 746 // Print total energy and its components 747 printf("\n"); 748 printf("%-21s energy: % .16e\n","Kinetic",sol.en.Ekin); 749 printf("%-21s energy: % .16e\n","Nuclear attraction",sol.en.Enuca); 750 printf("%-21s energy: % .16e\n","Total one-electron",sol.en.Eone); 751 printf("%-21s energy: % .16e\n","Nuclear repulsion",sol.en.Enucr); 752 printf("%-21s energy: % .16e\n","Coulomb",sol.en.Ecoul); 753#ifndef DFT 754 printf("%-21s energy: % .16e\n","Exchange",sol.en.Exc); 755#else 756 printf("%-21s energy: % .16e\n","Exchange-correlation",sol.en.Exc); 757 printf("%-21s energy: % .16e\n","Non-local correlation",sol.en.Enl); 758#endif 759 printf("-----------------------------------------------------\n"); 760 printf("%28s: % .16e\n","Total energy",sol.en.E); 761 printf("%28s: % .16e\n","Virial factor",-sol.en.E/sol.en.Ekin); 762 763 printf("\nDipole mu = (% 08.8f, % 08.8f, % 08.8f) D\n",dipmom(0)/AUINDEBYE,dipmom(1)/AUINDEBYE,dipmom(2)/AUINDEBYE); 764 765 printf("\n"); 766 // Print orbital energies 767#ifdef RESTRICTED 768 print_E(sol.E,occs,true); 769#else 770 printf("alpha: "); 771 print_E(sol.Ea,occa,true); 772 printf("beta: "); 773 print_E(sol.Eb,occb,true); 774#endif 775 776#if !defined(FULLHOLE) && !defined(HALFHOLE) 777 if(lincalc) { 778 // Classify occupied orbitals by symmetry 779#ifdef RESTRICTED 780 arma::imat occ(basisp->count_m_occupied(sol.C.cols(0,nocc-1))); 781#else 782 arma::imat occ(basisp->count_m_occupied(sol.Ca.cols(0,nocca-1),sol.Cb.cols(0,noccb-1))); 783#endif 784 785 printf("\nOrbital occupations by symmetry\n"); 786 int msum=0; 787 for(size_t im=0;im<occ.n_rows;im++) { 788 if(occ(im,0)+occ(im,1)>0) 789 printf("m = % i: %2i %2i\n",(int) occ(im,2),(int) occ(im,0),(int) occ(im,1)); 790 msum+=occ(im,2)*(occ(im,0)+occ(im,1)); 791 } 792 printf("Sum of m values is %i ",msum); 793 794 std::string termsymbol; 795 msum=std::abs(msum); 796 if(msum==0) 797 termsymbol="Sigma"; 798 else if(msum==1) 799 termsymbol="Pi"; 800 else if(msum==2) 801 termsymbol="Delta"; 802 else if(msum==3) 803 termsymbol="Phi"; 804 else if(msum==4) 805 termsymbol="Gamma"; 806 else { 807 std::ostringstream oss; 808 oss << "M=" << msum; 809 termsymbol=oss.str(); 810 } 811 printf("i.e. this is a %s state\n",termsymbol.c_str()); 812 } 813#endif 814 } 815 816 if(converged) { 817 if(doforce) { 818 arma::vec f; 819#if defined(RESTRICTED) && defined(HF) 820 f=force_RHF(sol,occs,intthr); 821#endif 822#if defined(UNRESTRICTED) && defined(HF) 823 f=force_UHF(sol,occa,occb,intthr); 824#endif 825#if defined(RESTRICTED) && defined(DFT) 826 f=force_RDFT(sol,occs,dft,grid,nlgrid,intthr); 827#endif 828#if defined(UNRESTRICTED) && defined(DFT) 829 f=force_UDFT(sol,occa,occb,dft,grid,nlgrid,intthr); 830#endif 831#if defined(_ROHF) || defined(FULLHOLE) || defined(HALFHOLE) 832 ERROR_INFO(); 833 throw std::runtime_error("Forces not supported for this method.\n"); 834#endif 835 836 chkptp->write("Force",f); 837 } 838 839 } else if(maxiter>0) { 840 std::ostringstream oss; 841 oss << "Error in function " << __FUNCTION__ << " (file " << __FILE__ << ", near line " << __LINE__ << "): SCF did not converge in "<<maxiter<<" iterations!\n"; 842 throw std::runtime_error(oss.str()); 843 } 844 845#if defined(HALFHOLE) || defined(FULLHOLE) 846 return ixc_orb; 847#endif 848} 849