1 subroutine mp2_lai_fock_uhf_prepar(g_pab_a, g_pab_b, 2 $ g_pij_a, g_pij_b, no_a, no_b, nv_a, nv_b, 3 $ g_p_a, g_p_b, nmo) 4 implicit none 5#include "errquit.fh" 6#include "global.fh" 7#include "mafdecls.fh" 8 integer g_pab_a, g_pab_b, 9 $ g_pij_a, g_pij_b, no_a, no_b, nv_a, nv_b, 10 $ g_p_a, g_p_b, nmo 11c 12c Put the OO and VV blocks of P into the appropriate 13c diagonal blocks of a density matrix which is also 14c allocated here. This is in preparation for building 15c L3+l4. It used to be part of lai_fock_uhf but has 16c been pulled out so that that routine can be used by the RIMP2. 17c 18*ga:1:0 19 if (.not. ga_create(mt_dbl, nmo, nmo, 'mp2_grad: p_a', 20 $ 0, 0, g_p_a)) call errquit('mp2_lai_fock_uhf: GA', nmo*nmo, 21 & GA_ERR) 22*ga:1:0 23 if (.not. ga_create(mt_dbl, nmo, nmo, 'mp2_grad: p_b', 24 $ 0, 0, g_p_b)) call errquit('mp2_lai_fock_uhf: GA', nmo*nmo, 25 & GA_ERR) 26c 27 call ga_zero(g_p_a) 28 call ga_copy_patch('n', 29 $ g_pij_a, 1, no_a, 1, no_a, g_p_a, 1, no_a, 1, no_a) 30 call ga_copy_patch('n', 31 $ g_pab_a, 1, nv_a, 1, nv_a, g_p_a, no_a+1, nmo, no_a+1, nmo) 32c 33 call ga_zero(g_p_b) 34 call ga_copy_patch('n', 35 $ g_pij_b, 1, no_b, 1, no_b, g_p_b, 1, no_b, 1, no_b) 36 call ga_copy_patch('n', 37 $ g_pab_b, 1, nv_b, 1, nv_b, g_p_b, no_b+1, nmo, no_b+1, nmo) 38c 39 end 40 subroutine mp2_lai_fock_uhf_tidy(g_p_a,g_p_b) 41 implicit none 42#include "errquit.fh" 43#include "global.fh" 44 integer g_p_a,g_p_b 45 if (.not. ga_destroy(g_p_a)) callerrquit('mp2_l_f_u_t: GA?',0, 46 & GA_ERR) 47 if (.not. ga_destroy(g_p_b)) callerrquit('mp2_l_f_u_t: GA?',1, 48 & GA_ERR) 49 end 50 subroutine mp2_lai_fock_uhf(geom, basis, 51 $ g_p_a, g_p_b, g_vecs_a, g_vecs_b, 52 $ no_a, no_b, nv_a, nv_b, 53 $ g_lai_a, g_lai_b, rtdb, tol2e) 54* 55* $Id$ 56* 57 implicit none 58#include "errquit.fh" 59#include "global.fh" 60#include "mafdecls.fh" 61#include "geom.fh" 62#include "bas.fh" 63 integer geom, basis ! [input] 64 integer rtdb 65 double precision tol2e 66 integer g_p_a, g_p_b ! [input] GA densities 67 integer g_vecs_a, g_vecs_b ! [input] MO vector handles 68 integer no_a, no_b, nv_a, nv_b ! [input] occ/virt dimensions 69 integer g_lai_a, g_lai_b ! [output] Accumulate results into 70 integer ga_create_atom_blocked 71 external ga_create_atom_blocked 72 73c g_lai += J[Paa] - K[Paa] + J[Pbb] 74c 75 integer g_dens(4), g_fock(4), nfock, nmo, nbf, gtype 76 integer g_tmp, i 77 double precision jfac(4), kfac(4), one, zero, mone 78 logical oskel 79c 80 data nfock /4/ 81 data jfac / 1.0d0, 1.0d0, 1.0d0, 1.0d0/ 82 data kfac /-1.0d0, 0.0d0,-1.0d0, 0.0d0/ 83 data oskel /.false./ 84c 85c 86 one = 1.0d0 87 zero= 0.0d0 88 mone = -1.0d0 89c 90c Allocate space for AO fock and density matrices 91c 92 g_dens(1)=ga_create_atom_blocked(geom,basis,'mp2_grad:da') 93 g_dens(2)=g_dens(1) 94 g_dens(3)=ga_create_atom_blocked(geom,basis,'mp2_grad:db') 95 g_dens(4)=g_dens(3) 96 g_fock(1)=ga_create_atom_blocked(geom,basis,'mp2_grad:fa') 97 g_fock(2)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 98 g_fock(3)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 99 g_fock(4)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 100 do i = 1, 4 101 call ga_zero(g_dens(i)) 102 call ga_zero(g_fock(i)) 103 end do 104 call scf_get_fock_param(rtdb,tol2e) 105 call fock_force_direct(rtdb) ! Force Fock build to be direct 106c 107c Workspace 108c 109 call ga_inquire(g_vecs_a, gtype, nbf, nmo) 110 if (no_a+nv_a.ne.nmo .or. no_b+nv_b.ne.nmo) 111 $ call errquit('mp2_grad: weird nmo?',0, INPUT_ERR) 112c 113*ga:1:0 114 if (.not. ga_create(mt_dbl, nbf, nmo, 'mp2_grad: tmp', 115 $ 0, 0, g_tmp)) call errquit('mp2_lai_fock_uhf: GA', nmo*nmo, 116 & GA_ERR) 117c 118 call ga_dgemm('n', 'n', nbf, nmo, nmo, one, g_vecs_a, 119 $ g_p_a, zero, g_tmp) 120 call ga_dgemm('n', 't', nbf, nbf, nmo, one, g_tmp, 121 $ g_vecs_a, zero, g_dens(1)) 122 call ga_symmetrize(g_dens(1)) 123 call ga_dscal(g_dens(1),mone) 124c 125 call ga_dgemm('n', 'n', nbf, nmo, nmo, one, g_vecs_b, 126 $ g_p_b, zero, g_tmp) 127 call ga_dgemm('n', 't', nbf, nbf, nmo, one, g_tmp, 128 $ g_vecs_b, zero, g_dens(3)) 129 call ga_symmetrize(g_dens(3)) 130 call ga_dscal(g_dens(3),mone) 131c 132 if (.not. ga_destroy(g_tmp)) call errquit('mp2_lfu: ga?',0, 133 & GA_ERR) 134c 135c Make the fock matrices 136c 137 call int_init(rtdb,1,basis) 138 call schwarz_init(geom,basis) 139 140 call fock_2e(geom, basis, nfock, jfac, kfac, 141 $ tol2e, oskel, g_dens, g_fock,.false. ) 142c 143 call int_terminate() 144 call schwarz_tidy() 145 call fock_2e_tidy(rtdb) 146 147 call ga_dadd(one, g_fock(1), one, g_fock(4), g_fock(1)) 148 call ga_dadd(one, g_fock(3), one, g_fock(2), g_fock(3)) 149 call ga_symmetrize(g_fock(1)) 150 call ga_symmetrize(g_fock(3)) 151c 152c Transform back 153c 154 g_tmp = g_dens(1) 155c 156 call ga_zero(g_tmp) 157 call ga_matmul_patch('n','n',one, zero, 158 $ g_fock(1), 1, nbf, 1, nbf, 159 $ g_vecs_a, 1, nbf, no_a+1, nmo, 160 $ g_tmp, 1, nbf, 1, nv_a) 161 call ga_matmul_patch('t','n', one, one, 162 $ g_vecs_a, 1, no_a, 1, nbf, 163 $ g_tmp, 1, nbf, 1, nv_a, 164 $ g_lai_a, 1, no_a, 1, nv_a) 165c 166 call ga_zero(g_tmp) 167 call ga_matmul_patch('n','n',one, zero, 168 $ g_fock(3), 1, nbf, 1, nbf, 169 $ g_vecs_b, 1, nbf, no_b+1, nmo, 170 $ g_tmp, 1, nbf, 1, nv_b) 171 call ga_matmul_patch('t','n', one, one, 172 $ g_vecs_b, 1, no_b, 1, nbf, 173 $ g_tmp, 1, nbf, 1, nv_b, 174 $ g_lai_b, 1, no_b, 1, nv_b) 175c 176 if (.not. ga_destroy(g_dens(1))) call errquit('mp2_lfu: ga?',0, 177 & GA_ERR) 178 if (.not. ga_destroy(g_dens(3))) call errquit('mp2_lfu: ga?',0, 179 & GA_ERR) 180 if (.not. ga_destroy(g_fock(1))) call errquit('mp2_lfu: ga?',0, 181 & GA_ERR) 182 if (.not. ga_destroy(g_fock(2))) call errquit('mp2_lfu: ga?',0, 183 & GA_ERR) 184 if (.not. ga_destroy(g_fock(3))) call errquit('mp2_lfu: ga?',0, 185 & GA_ERR) 186 if (.not. ga_destroy(g_fock(4))) call errquit('mp2_lfu: ga?',0, 187 & GA_ERR) 188c 189 end 190 subroutine mp2_wij_fock_uhf(rtdb, geom, basis, tol2e, 191 $ g_ppq_a, g_ppq_b, 192 $ no_a, no_b, 193 $ g_vecs_a, g_vecs_b, 194 $ g_wij_a, g_wij_b) 195 implicit none 196#include "errquit.fh" 197#include "global.fh" 198#include "mafdecls.fh" 199 integer rtdb, geom, basis ! [input] 200 double precision tol2e ! [input] 201 integer no_a, no_b ! [input] No. of occupied 202 integer g_ppq_a, g_ppq_b ! [input] GA densities 203 integer g_vecs_a, g_vecs_b ! [input] MO vector handles 204 integer g_wij_a, g_wij_b ! [output] Accumulate results into 205 integer ga_create_atom_blocked 206 external ga_create_atom_blocked 207c 208c w_ij += J[Paa] - K[Paa] + J[Pbb] 209c 210 integer g_dens(4), g_fock(4), nfock, nmo, nbf, gtype 211 integer g_tmp1, g_tmp2, i 212 double precision jfac(4), kfac(4), one, zero, mone 213 logical oskel 214c 215 data nfock /4/ 216 data jfac / 1.0d0, 1.0d0, 1.0d0, 1.0d0/ 217 data kfac /-1.0d0, 0.0d0,-1.0d0, 0.0d0/ 218 data oskel /.false./ 219c 220 one = 1.0d0 221 zero= 0.0d0 222 mone = -1.0d0 223c 224c Allocate space for AO fock and density matrices 225c 226 g_dens(1)=ga_create_atom_blocked(geom,basis,'mp2_grad:da') 227 g_dens(2)=g_dens(1) 228 g_dens(3)=ga_create_atom_blocked(geom,basis,'mp2_grad:db') 229 g_dens(4)=g_dens(3) 230 g_fock(1)=ga_create_atom_blocked(geom,basis,'mp2_grad:fa') 231 g_fock(2)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 232 g_fock(3)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 233 g_fock(4)=ga_create_atom_blocked(geom,basis,'mp2_grad:fb') 234 do i = 1, 4 235 call ga_zero(g_dens(i)) 236 call ga_zero(g_fock(i)) 237 end do 238 call scf_get_fock_param(rtdb,tol2e) 239 call fock_force_direct(rtdb) ! Force Fock build to be direct 240c 241c Workspace 242c 243 call ga_inquire(g_vecs_a, gtype, nbf, nmo) 244c 245*ga:1:0 246 if (.not. ga_create(mt_dbl, nbf, nmo, 'mp2_grad: tmp2', 247 $ 0, 0, g_tmp2)) call errquit('mp2_wij_fock_uhf: GA', nmo*nmo, 248 & GA_ERR) 249c 250c Transform the densities and symmetrize 251c 252 call ga_dgemm('n', 'n', nbf, nmo, nmo, one, g_vecs_a, 253 $ g_ppq_a, zero, g_tmp2) 254 call ga_dgemm('n', 't', nbf, nbf, nmo, mone, g_tmp2, 255 $ g_vecs_a, zero, g_dens(1)) 256 call ga_symmetrize(g_dens(1)) 257c 258 call ga_dgemm('n', 'n', nbf, nmo, nmo, one, g_vecs_b, 259 $ g_ppq_b, zero, g_tmp2) 260 call ga_dgemm('n', 't', nbf, nbf, nmo, mone, g_tmp2, 261 $ g_vecs_b, zero, g_dens(3)) 262 call ga_symmetrize(g_dens(3)) 263c 264 if (.not. ga_destroy(g_tmp2)) call errquit('mp2_lfu: ga?',0, 265 & GA_ERR) 266c 267c Make the fock matrices 268c 269 call int_init(rtdb,1,basis) 270 call schwarz_init(geom,basis) 271 272 call fock_2e(geom, basis, nfock, jfac, kfac, 273 $ tol2e, oskel, g_dens, g_fock, .false.) 274c 275 call int_terminate() 276 call schwarz_tidy() 277 call fock_2e_tidy(rtdb) 278 279 call ga_dadd(one, g_fock(1), one, g_fock(4), g_fock(1)) 280 call ga_dadd(one, g_fock(3), one, g_fock(2), g_fock(3)) 281 call ga_symmetrize(g_fock(1)) 282 call ga_symmetrize(g_fock(3)) 283c 284c Transform back ... just want wij 285c 286 g_tmp1 = g_dens(1) 287c 288 call ga_matmul_patch('n','n',one, zero, 289 $ g_fock(1), 1, nbf, 1, nbf, 290 $ g_vecs_a, 1, nbf, 1, no_a, 291 $ g_tmp1, 1, nbf, 1, no_a) 292 call ga_matmul_patch('t','n', one, one, 293 $ g_vecs_a, 1, no_a, 1, nbf, 294 $ g_tmp1, 1, nbf, 1, no_a, 295 $ g_wij_a, 1, no_a, 1, no_a) 296c 297 call ga_matmul_patch('n','n',one, zero, 298 $ g_fock(3), 1, nbf, 1, nbf, 299 $ g_vecs_b, 1, nbf, 1, no_b, 300 $ g_tmp1, 1, nbf, 1, no_b) 301 call ga_matmul_patch('t','n', one, one, 302 $ g_vecs_b, 1, no_b, 1, nbf, 303 $ g_tmp1, 1, nbf, 1, no_b, 304 $ g_wij_b, 1, no_b, 1, no_b) 305c 306 if (.not. ga_destroy(g_dens(1))) call errquit('mp2_lfu: ga?',0, 307 & GA_ERR) 308 if (.not. ga_destroy(g_dens(3))) call errquit('mp2_lfu: ga?',0, 309 & GA_ERR) 310 if (.not. ga_destroy(g_fock(1))) call errquit('mp2_lfu: ga?',0, 311 & GA_ERR) 312 if (.not. ga_destroy(g_fock(2))) call errquit('mp2_lfu: ga?',0, 313 & GA_ERR) 314 if (.not. ga_destroy(g_fock(3))) call errquit('mp2_lfu: ga?',0, 315 & GA_ERR) 316 if (.not. ga_destroy(g_fock(4))) call errquit('mp2_lfu: ga?',0, 317 & GA_ERR) 318c 319 end 320 321