1C> \ingroup wfn1_nxt 2C> @{ 3C> 4C> \brief Compute the energy and gradient of the current state vector 5C> 6C> Computes the energy and gradient of the state vector given the 7C> various sets of integrals. This routine can generate various 8C> energy expressions including: 9C> 1. Hartree-Fock 10C> 2. Density Functional Theory 11C> 3. Density Matrix Functional Theory 12C> 4. Entropy model 13C> 14 subroutine wfn1_nxt_energy_gradient(rtdb,geom,nbf,nea,neb,h1,eri, 15 & erix,ov,ov12,state,grad,Etot) 16 implicit none 17c 18#include "mafdecls.fh" 19#include "errquit.fh" 20#include "wfn1_nxt.fh" 21#include "geom.fh" 22c 23 integer rtdb !< [Input] The RTDB handle (needed for the DFT terms) 24 integer geom !< [Input] The geometry handle (needed for the 25 !< nuclear repulsion energy) 26 integer nbf !< [Input] The number of basis functions 27 integer nea !< [Input] The number of alpha electrons 28 integer neb !< [Input] The number of beta electrons 29c 30 double precision h1(nbf,nbf) !< [Input] The 1-electron 31 !< integrals 32 double precision eri(nbf,nbf,nbf,nbf) !< [Input] The Coulomb 33 !< integrals 34 double precision erix(nbf,nbf,nbf,nbf) !< [Input] The exchange 35 !< integrals 36 double precision ov(nbf,nbf) !< [Input] The overlap 37 !< integrals 38 double precision ov12(nbf,nbf) !< [Input] The 1/sqrt(overlap) 39c 40 double precision state(8*nbf*nbf) !< [Input] The state vector 41 double precision grad(8*nbf*nbf) !< [Output] The gradient 42cDEBUG 43c double precision gral(8*nbf*nbf) !< [Scratch] The gradient 44cDEBUG 45c 46 double precision Etot !< [Output] The total energy 47c 48 integer l_occa, k_occa !< Memory for alpha electron occupation 49 !< numbers 50 integer l_occb, k_occb !< Memory for beta electron occupation 51 !< numbers 52c 53 integer l_da, k_da !< Memory for the alpha electron density 54 !< matrix 55 integer l_db, k_db !< Memory for the beta electron density 56 !< matrix 57c 58 integer l_fa, k_fa !< Memory for the alpha electron Fock 59 !< matrix 60 integer l_fb, k_fb !< Memory for the beta electron Fock 61 !< matrix 62c 63 integer l_ta, k_ta !< Scratch matrix 64 integer l_tb, k_tb !< Scratch matrix 65c 66 integer k_tv !< Scratch vector 67c 68 double precision Enuc !< The nuclear repulsion energy 69 double precision E1el !< The 1-electron energy 70 double precision E2elC !< The 2-electron Coulomb energy 71 double precision E2elX !< The 2-electron eXchange energy 72 double precision ELo !< The occupation function Lagrangian 73 double precision ELn !< The natural orbital Lagrangian 74c 75 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Da',l_da,k_da)) 76 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 77 & //"Da",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 78 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Db',l_db,k_db)) 79 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 80 & //"Db",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 81c 82 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Fa',l_fa,k_fa)) 83 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 84 & //"Fa",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 85 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Fb',l_fb,k_fb)) 86 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 87 & //"Fb",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 88c 89 if (.not.ma_push_get(MT_DBL,nbf,'occ a',l_occa,k_occa)) 90 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 91 & //"occa",ma_sizeof(MT_DBL,nbf,MT_BYTE),MA_ERR) 92 if (.not.ma_push_get(MT_DBL,nbf,'occ b',l_occb,k_occb)) 93 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 94 & //"occb",ma_sizeof(MT_DBL,nbf,MT_BYTE),MA_ERR) 95 k_tv = k_occa 96c 97 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Ta',l_ta,k_ta)) 98 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 99 & //"Ta",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 100 if (.not.ma_push_get(MT_DBL,nbf*nbf,'Tb',l_tb,k_tb)) 101 & call errquit("wfn1_nxt_energy_gradient: could not allocate " 102 & //"Tb",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 103c 104c Zero the gradient 105c 106 call dfill(wfn1_nxt_size1(),0.0d0,grad,1) 107c 108c Get the nuclear repulsion energy 109c 110 if (.not.geom_nuc_rep_energy(geom,Enuc)) 111 & call errquit("wfn1_energy: Enuc failed",0,GEOM_ERR) 112c 113c Compute the 1-electron energy and gradient 114c 115cDEBUG 116c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 117cDEBUG 118 call wfn1_nxt_1el_eg(nbf,nea,neb,h1, 119 & dbl_mb(k_occa),dbl_mb(k_occb),dbl_mb(k_da),dbl_mb(k_db), 120 & state(wfn1_nxt_aocc()),state(wfn1_nxt_anat()), 121 & state(wfn1_nxt_bocc()),state(wfn1_nxt_bnat()), 122 & grad(wfn1_nxt_aocc()),grad(wfn1_nxt_anat()), 123 & grad(wfn1_nxt_bocc()),grad(wfn1_nxt_bnat()),E1el) 124cDEBUG 125c write(*,*)'*** 1-electron energy' 126c call wfn1_nxt_print_state(gral,nbf) 127c call daxpy(wfn1_nxt_size1(),1.0d0,gral,1,grad,1) 128c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 129cDEBUG 130 call wfn1_nxt_2el_C_eg(nbf,nea,neb,eri, 131 & dbl_mb(k_occa),dbl_mb(k_occb),dbl_mb(k_da),dbl_mb(k_db), 132 & state(wfn1_nxt_aocc()),state(wfn1_nxt_anat()), 133 & state(wfn1_nxt_bocc()),state(wfn1_nxt_bnat()), 134 & grad(wfn1_nxt_aocc()),grad(wfn1_nxt_anat()), 135 & grad(wfn1_nxt_bocc()),grad(wfn1_nxt_bnat()),E2elC, 136 & dbl_mb(k_fa),dbl_mb(k_fb)) 137cDEBUG 138c write(*,*)'*** 2-electron Coulomb energy' 139c call wfn1_nxt_print_state(gral,nbf) 140c call daxpy(wfn1_nxt_size1(),1.0d0,gral,1,grad,1) 141c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 142cDEBUG 143 call wfn1_nxt_2el_X_eg(nbf,nea,neb,wfn1_nxt_x_pow,erix, 144 & dbl_mb(k_occa),dbl_mb(k_occb),dbl_mb(k_da),dbl_mb(k_db), 145 & state(wfn1_nxt_aocc()),state(wfn1_nxt_anat()), 146 & state(wfn1_nxt_bocc()),state(wfn1_nxt_bnat()), 147 & grad(wfn1_nxt_aocc()),grad(wfn1_nxt_anat()), 148 & grad(wfn1_nxt_bocc()),grad(wfn1_nxt_bnat()),E2elX, 149 & dbl_mb(k_fa),dbl_mb(k_fb)) 150cDEBUG 151c write(*,*)'*** 2-electron eXchange energy' 152c call wfn1_nxt_print_state(gral,nbf) 153c call daxpy(wfn1_nxt_size1(),1.0d0,gral,1,grad,1) 154c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 155cDEBUG 156c 157c call wfn1_nxt_Lo_eg(nbf, 158c & state(wfn1_nxt_aocc()),state(wfn1_nxt_bocc()), 159c & state(wfn1_nxt_aoccl()),state(wfn1_nxt_boccl()), 160c & grad(wfn1_nxt_aocc()),grad(wfn1_nxt_bocc()), 161c & grad(wfn1_nxt_aoccl()),grad(wfn1_nxt_boccl()),ELo) 162c call wfn1_nxt_Ln_eg(nbf,ov, 163c & state(wfn1_nxt_anat()),state(wfn1_nxt_bnat()), 164c & state(wfn1_nxt_anatl()),state(wfn1_nxt_bnatl()), 165c & grad(wfn1_nxt_anat()),grad(wfn1_nxt_bnat()), 166c & grad(wfn1_nxt_anatl()),grad(wfn1_nxt_bnatl()),ELn, 167c & dbl_mb(k_tv)) 168c 169 call wfn1_nxt_Po_eg(nbf, 170 & state(wfn1_nxt_aocc()),state(wfn1_nxt_bocc()), 171 & state(wfn1_nxt_aoccl()),state(wfn1_nxt_boccl()), 172 & grad(wfn1_nxt_aocc()),grad(wfn1_nxt_bocc()), 173 & grad(wfn1_nxt_aoccl()),grad(wfn1_nxt_boccl()),ELo, 174 & dbl_mb(k_ta),dbl_mb(k_tb)) 175cDEBUG 176c write(*,*)'*** occupation function penalty' 177c call wfn1_nxt_print_state(gral,nbf) 178c call daxpy(wfn1_nxt_size1(),1.0d0,gral,1,grad,1) 179c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 180cDEBUG 181 call wfn1_nxt_Pn_eg(nbf,ov, 182 & state(wfn1_nxt_anat()),state(wfn1_nxt_bnat()), 183 & state(wfn1_nxt_anatl()),state(wfn1_nxt_bnatl()), 184 & grad(wfn1_nxt_anat()),grad(wfn1_nxt_bnat()), 185 & grad(wfn1_nxt_anatl()),grad(wfn1_nxt_bnatl()),ELn, 186 & dbl_mb(k_tv), 187 & dbl_mb(k_ta),dbl_mb(k_tb)) 188cDEBUG 189c write(*,*)'*** natural orbital penalty' 190c call wfn1_nxt_print_state(gral,nbf) 191c call daxpy(wfn1_nxt_size1(),1.0d0,gral,1,grad,1) 192c call dfill(wfn1_nxt_size1(),0.0d0,gral,1) 193cDEBUG 194c 195 if (.not.ma_chop_stack(l_da)) 196 & call errquit("wfn1_nxt_energy_gradient: could not deallocate " 197 & //"Da",ma_sizeof(MT_DBL,nbf*nbf,MT_BYTE),MA_ERR) 198c 199 Etot = Enuc + E1el + E2elC + E2elX + ELo + ELn 200cDEBUG 201 write(*,*)'Enuc = ',Enuc 202 write(*,*)'E1el = ',E1el 203 write(*,*)'E2elC = ',E2elC 204 write(*,*)'E2elX = ',E2elX 205 write(*,*)'EPo = ',ELo 206 write(*,*)'EPn = ',ELn 207 write(*,*)'Etot = ',Etot 208cDEBUG 209c 210 end 211C> 212C> @} 213