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