1/* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * HF/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/* Restricted case */ 18#if defined(RESTRICTED) && defined(DFT) 19arma::vec SCF::force_RDFT(rscf_t & sol, const std::vector<double> & occs, const dft_t dft, DFTGrid & grid, DFTGrid & nlgrid, double tol) 20 21#elif defined(RESTRICTED) && defined(HF) 22arma::vec SCF::force_RHF(rscf_t & sol, const std::vector<double> & occs, double tol) 23 24#elif defined(UNRESTRICTED) && defined(DFT) 25arma::vec SCF::force_UDFT(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, const dft_t dft, DFTGrid & grid, DFTGrid & nlgrid, double tol) 26 27#elif defined(UNRESTRICTED) && defined(HF) 28arma::vec SCF::force_UHF(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, double tol) 29 30#elif defined(UNRESTRICTED) && defined(_ROHF) 31arma::vec SCF::force_ROHF(uscf_t & sol, int Nel_alpha, int Nel_beta, double tol) 32#endif 33{ 34 35 arma::mat W; 36#ifdef RESTRICTED 37 W=form_density(sol.E,sol.C,occs); 38#else 39 W=form_density(sol.Ea,sol.Ca,occa)+form_density(sol.Eb,sol.Cb,occb); 40#endif 41 42 // Compute force 43 arma::vec fpul_kin=basisp->kinetic_pulay(sol.P); 44 // interpret_force(fpul_kin).print("Kinetic Pulay"); 45 46 arma::vec fpul_nuc=basisp->nuclear_pulay(sol.P); 47 // interpret_force(fpul_nuc).print("Nuclear Pulay"); 48 49 arma::vec fnuc=basisp->nuclear_der(sol.P); 50 // interpret_force(fnuc).print("Hellman-Feynman"); 51 52 arma::vec forth=basisp->overlap_der(W); 53 // interpret_force(forth).print("Orthonormality"); 54 55 arma::vec frep=basisp->nuclear_force(); 56 // interpret_force(frep).print("Nuclear repulsion"); 57 58#ifdef DFT 59 // Get range separation info 60 double omega, kfull, kshort; 61 range_separation(dft.x_func,omega,kfull,kshort); 62#else 63 double kfull=1.0; 64#endif 65 66 // Exact exchange and Coulomb 67 arma::vec fx_full; 68 fx_full.zeros(fnuc.n_elem); 69 70#ifdef DFT 71 // Short-range exchange 72 arma::vec fx_short; 73 fx_short.zeros(fnuc.n_elem); 74 75 if(kfull==0.0) { 76 77 if(densityfit) { 78 if(kshort!=0.0 || kfull!=0.0) 79 throw std::runtime_error("Forces not implemented for density fitting of exact exchange.\n"); 80 81 fx_full=dfit.forceJ(sol.P); 82 } else { 83 if(!direct) 84 scr.fill(basisp,intthr,verbose); 85 86 fx_full=scr.forceJ(sol.P,tol); 87 } 88 } else { 89#endif 90 91 if(!direct) 92 scr.fill(basisp,intthr,verbose); 93 94#ifdef RESTRICTED 95 fx_full=scr.forceJK(sol.P,tol,kfull); 96#else 97 fx_full=scr.forceJK(sol.Pa,sol.Pb,tol,kfull); 98#endif 99 100#ifdef DFT 101 } 102 103 if(omega != 0.0) { 104 // Set range separation and fill (low-cost compared to forces anyway) 105 scr_rs.set_range_separation(omega,0.0,1.0); 106 scr_rs.fill(basisp,intthr,verbose); 107 108 // Exchange forces 109#ifdef RESTRICTED 110 fx_short=scr_rs.forceK(sol.P,tol,kshort); 111#else 112 fx_short=scr_rs.forceK(sol.Pa,sol.Pb,tol,kshort); 113#endif 114 } 115 116 // Get the DFT contribution 117 arma::vec fxc; 118 fxc.zeros(fnuc.n_elem); 119 if(dft.x_func>0 || dft.c_func>0) { 120#ifdef RESTRICTED 121 fxc=grid.eval_force(dft.x_func,dft.c_func,sol.P); 122#else 123 fxc=grid.eval_force(dft.x_func,dft.c_func,sol.Pa,sol.Pb); 124#endif 125 } 126 127 // Non-local contribution to force 128 if(dft.nl) { 129 arma::vec vv10f(grid.eval_VV10_force(nlgrid,dft.vv10_b,dft.vv10_C,sol.P)); 130 //interpret_force(vv10f).t().print("VV10 force"); 131 fxc+=vv10f; 132 } 133 134#endif // DFT 135 136 arma::vec ftot=fpul_kin+fpul_nuc+fnuc+forth+frep; 137 arma::vec ffull=ftot+fx_full; 138#ifdef DFT 139 // DFT contribution 140 ffull+=fxc+fx_short; 141#endif 142 143 if(lincalc) { 144 // Get rid of numerical noise for linear molecule restriction 145 for(size_t inuc=0;inuc<ffull.n_elem/3;inuc++) { 146 ffull(3*inuc)=0.0; 147 ffull(3*inuc+1)=0.0; 148 } 149 } 150 151 // interpret_force(fx_full).print("Coulomb + exact exchange"); 152 // interpret_force(ftot).print("Total (w.o. 2-electron)"); 153 // interpret_force(ffull).t().print("Total force"); 154 155 return ffull; 156} 157