1c 2c 3c Form 4c 5c F = h core Hamiltonian 6c ij ij 7c 8c rs rs 9c F = h + sum ( 2J - K ) inactive Fock element 10c rs rs i ii ii 11c 12c 13 subroutine mcscf_fcore( basis, nbf, nclosed, nact, 14 $ g_movecs, g_coul, g_exch, g_fcore ) 15* 16* $Id$ 17* 18 implicit none 19#include "errquit.fh" 20#include "mafdecls.fh" 21#include "global.fh" 22#include "util.fh" 23 integer basis 24 integer nbf 25 integer nclosed 26 integer nact 27 integer g_movecs 28 integer g_coul 29 integer g_exch 30 integer g_fcore 31 integer nn, t, u, tu, jlo, jhi, k_j, k_k, ld1 32 logical ga_check_JKblocked 33 external ga_check_JKblocked 34c 35c 1e Hamiltononian 36c 37 if (util_print('fcore', print_never)) then 38 write(6,*) ' MOVECS in fcore ' 39 call ga_print(g_movecs) 40 endif 41 call moints_1e( nbf, basis, g_movecs, g_fcore) 42 if (util_print('fcore', print_never)) call ga_print(g_fcore) 43c 44c Inactive-active Coulomb interation 45c 46 nn = nbf*nbf 47 if (.not.ga_check_JKblocked(g_coul,nact,nbf,jlo,jhi)) 48 $ call errquit('mcscf_fcore: wrong distrib. for Coulomb',0, 49 & GA_ERR) 50 do t=1,nact 51 do u=1,t 52 tu = ((t-1)*t)/2 + u 53 if ((tu.ge.jlo).and.(tu.le.jhi)) then 54 call ga_access(g_coul,1,nn,tu,tu,k_j,ld1) 55 call mcscf_fcore01( nbf, nclosed, nact, t, u, 2.d0, 56 $ dbl_mb(k_j), g_fcore) 57 call ga_release(g_coul,1,nn,tu,tu) 58 endif 59 enddo 60 enddo 61c 62c Inactive-active Exchange interaction 63c 64 if (.not.ga_check_JKblocked(g_exch,nact,nbf,jlo,jhi)) 65 $ call errquit('mcscf_fcore: wrong distrib. for exchange',0, 66 & GA_ERR) 67 do t=1,nact 68 do u=1,t 69 tu = ((t-1)*t)/2 + u 70 if ((tu.ge.jlo).and.(tu.le.jhi)) then 71 call ga_access(g_exch,1,nn,tu,tu,k_k,ld1) 72 call mcscf_fcore01( nbf, nclosed, nact, t, u, -1.d0, 73 $ dbl_mb(k_k), g_fcore) 74 call ga_release(g_exch,1,nn,tu,tu) 75 endif 76 enddo 77 enddo 78c 79c Properly synchronize 80c 81 call ga_sync() 82c 83c Done 84c 85 return 86 end 87 88 89 90 91 92 93 subroutine mcscf_fcore01( nbf, nclosed, nact, t, u, fact, 94 $ xmo, g_fcore ) 95 implicit none 96#include "global.fh" 97 integer nbf 98 integer nclosed 99 integer nact 100 integer t, u 101 double precision fact 102 double precision xmo(nbf,nbf) 103 integer g_fcore 104 double precision xx 105 integer i, tt, uu 106 107 xx = 0.d0 108 do i=1,nclosed 109 xx = xx + xmo(i,i) 110 enddo 111 tt = nclosed + t 112 uu = nclosed + u 113 call ga_acc(g_fcore, tt, tt, uu, uu, xx, 1, fact ) 114 if (t.ne.u) call ga_acc(g_fcore, uu, uu, tt, tt, xx, 1, fact ) 115 return 116 end 117 118 119 120 121 122 123 124 125 126 127 128 129 130c 131c Evaluate inactive Fock piece separately 132c (see mcscf_fock) 133c 134c 135 subroutine mcscf_ifock( geom, basis, nbf, nclosed, nact, 136 $ oskel, tol2e, g_movecs, eone, etwo, 137 $ ecore, g_ifock ) 138 implicit none 139#include "errquit.fh" 140#include "mafdecls.fh" 141#include "global.fh" 142#include "msgids.fh" 143#include "rtdb.fh" 144#include "bas.fh" 145#include "util.fh" 146#include "geom.fh" 147 integer geom, basis ! [input] Geometry and basis handles 148 integer nbf ! [input] Number of basis functions 149 integer nclosed ! [input] Number of closed shells 150 integer nact ! [input] Number of open shells 151 logical oskel ! [input] Symmetry toggle 152 double precision tol2e ! [input] Integral tolerance 153 integer g_movecs ! [input] MO coefficients 154 double precision eone, etwo, ecore ! [output] Energy components 155 integer g_ifock ! [output] Inactive Fock matrix 156c 157 integer nset 158 parameter(nset=1) 159c 160 integer g_cdens, g_tmp 161 double precision e2ii 162 integer i 163 double precision xx 164 integer iv_dens(2), iv_fock(2) 165 double precision jfac(2), kfac(2) 166 integer ga_create_atom_blocked 167 external ga_create_atom_blocked 168 data jfac/1.0d0, 1.d0/, kfac/-0.5d0, -0.5d0/ 169c 170c 171c 172 g_tmp = ga_create_atom_blocked(geom, basis, 'temp1') 173 g_cdens = ga_create_atom_blocked(geom, basis, 'closed dens') 174 call ga_dgemm( 'n', 't', nbf, nbf, nclosed, 2.d0, 175 $ g_movecs, g_movecs, 0.d0, g_cdens ) 176c 177c One-electron component 178c 179 call ga_zero(g_ifock) 180 call int_1e_ga( basis, basis, g_ifock, 'kinetic', oskel) 181 call int_1e_ga( basis, basis, g_ifock, 'potential', oskel) 182 if (oskel) call sym_symmetrize(geom, basis, .false., g_ifock) 183 eone = ga_ddot(g_ifock,g_cdens) 184c 185c Two-electron component of the AO Fock matrices 186c 187 call ga_zero(g_tmp) 188 iv_dens(1) = g_cdens 189 iv_fock(1) = g_tmp 190 call fock_2e( geom, basis, nset, jfac, kfac, tol2e, oskel, 191 $ iv_dens, iv_fock, .false. ) 192c 193c Symmetrize Fock AO components 194c 195 if (oskel) 196 $ call sym_symmetrize(geom, basis, .false., g_tmp ) 197 e2ii = ga_ddot(g_tmp,g_cdens) 198c 199c MO integral contribution 200c 201 etwo = e2ii*0.5d0 202 if ((ga_nodeid().eq.0).and. 203 $ (util_print('energy trace',print_debug))) then 204 write(6,911) e2ii*0.5d0 205 911 format( 'energy components eii: ',f12.6) 206 endif 207c 208c Inactive core energy 209c 210 call ga_dadd( 1.d0, g_tmp, 1.d0, g_ifock, g_cdens ) 211 call two_index_transf( g_cdens, g_movecs, g_movecs, 212 $ g_tmp, g_ifock ) 213 ecore = 0.d0 214 do i=ga_nodeid()+1,nclosed,ga_nnodes() 215 call ga_get(g_ifock,i,i,i,i,xx,1) 216 ecore = ecore + 2.d0*xx 217 enddo 218 call ga_sync() 219 call ga_dgop(msg_mcscf_ifocktrace,ecore,1,'+') 220 ecore = ecore - 0.5d0*e2ii 221c 222c Clean up 223c 224 if (.not.ga_destroy(g_tmp)) 225 $ call errquit('mcscf_etr: cannot destroy',0, GA_ERR) 226 if (.not.ga_destroy(g_cdens)) 227 $ call errquit('mcscf_etr: cannot destroy',0, GA_ERR) 228 return 229 end 230