1c
2c     $Id$
3c
4c****f* uhf/uhf_fock_setup
5c
6c     NAME
7c
8c       uhf_fock_setup -- set up pre-factors for Coulomb and exchange
9c       terms, density matrices and Fock matrices for the UHF
10c       shell_fock_build call
11c
12c     FUNCTION
13c
14c       This routine takes the multi-dimensional density and Fock
15c       matrices corresponding to alpha and beta spin channel
16c       quantities and expands them to alpha Coulomb and exchange
17c       matrices and beta Coulomb and exchange matrices.
18c
19c       After the call to shell_fock_build the sister routine
20c       uhf_fock_finish will combine the resulting matrices to
21c       construct alpha Fock matrices and beta Fock matrices.
22c
23c       The mathematics used in uhf_fock_finish imposes an interesting
24c       constraint on how we copy the input Fock matrices here in
25c       uhf_fock_setup. Looking for example at the alpha-Fock matrices,
26c       these are constructed as
27c
28c         alpha-Fock = alpha-Coulomb + beta-Coulomb - alpha-Exchange
29c
30c       This means that alpha-Coulomb and alpha-exchange should be
31c       initialized with the beta input Fock matrices (beta-Fock0 for
32c       short), likewise beta-Coulomb and beta-exchange should be
33c       initialized with alpha-Fock0. To verify this simply substitute
34c       those quantities in the above expression to get:
35c
36c         alpha-Fock = beta-Fock0 + alpha-Fock0 - beta-Fock0
37c                    = alpha-Fock0
38c
39c       to see that this way the input alpha-Fock matrices contribute
40c       to the output alpha-Fock matrices. There is no problem with
41c       respect to the Coulomb and exchange contributions that
42c       shell_fock_build has added, because they appear in the right
43c       places.
44c
45c       Schematically the operation of this routine can be presented
46c       as (ndens and nfock are the corresponding values on input):
47c
48c          Input                        Output
49c          =====                        ======
50c
51c          -------                      -------
52c            ^                            ^
53c            |                          ndens/2  alpha-density matrices
54c          ndens/2  alpha-density         v
55c            |      matrices            -------
56c            v                            ^
57c          -------                ===>  ndens/2  alpha-density matrices
58c            ^                            v
59c            |                          -------
60c          ndens/2  beta-density          ^
61c            |      matrices            ndens/2  beta-density matrices
62c            v                            v
63c          -------                      -------
64c                                         ^
65c                                       ndens/2  beta-density matrices
66c                                         v
67c                                       -------
68c
69c
70c                                       -------
71c                                         ^
72c          -------                      nfock/2  alpha-Coulomb matrices
73c            ^                            v
74c            |                          -------
75c          nfock/2  alpha-Fock            ^
76c            |      matrices         >  nfock/2  alpha-exchange matrices
77c            v                    \ /     v
78c          -------                 X    -------
79c            ^                    / \     ^
80c            |                       >  nfock/2  beta-Coulomb matrices
81c          nfock/2  beta-Fock             v
82c            |      matrices            -------
83c            v                            ^
84c          -------                      nfock/2  beta-exchange matrices
85c                                         v
86c                                       -------
87c
88c       The pre-factors for the Coulomb terms (jfac) and the exchange
89c       terms (kfac) undergo a corresponding copying operation as well.
90c
91c     SYNOPSIS
92c
93      subroutine uhf_fock_setup(g_dens, g_fock, jfac, kfac,
94     &                          ndens, nfock, nbf, oprint)
95c
96c     SOURCE
97c
98      implicit none
99c
100#include "errquit.fh"
101#include "global.fh"
102#include "mafdecls.fh"
103#include "stdio.fh"
104c
105c  This routine takes the multidimensional density and the fock matrices
106c  that are set up for alpha and beta and expand them to be alpha and
107c  beta/coulomb and exchange (2 -> 4) so that the correct summation can
108c  be performed after the shell_fock_build call.
109c
110c     INPUTS
111c
112      integer g_dens       ! [input/output] density
113      integer g_fock       ! [input/output] Fock matrices
114      integer ndens        ! [input/output] number of density matrices
115      integer nfock        ! [input/output] number of Fock matrices
116      integer nbf          ! [input] number of basis functions
117      double precision jfac(ndens) ! [input/output] Coulomb scale factors
118      double precision kfac(ndens) ! [input/output] Exchange scale factors
119      logical oprint       ! [input] Print output
120c
121c     SOURCE
122c
123      integer dims(3), chunk(3)
124      integer alo(3), ahi(3), blo(3), bhi(3)
125      integer g_dens2, g_fock2
126      integer i, noffset, noffset2, noffset3
127      double precision zero
128      parameter (zero = 0.0D0)
129c
130c     First set up jfac and kfac correctly.  This assumes that there is
131c     already enough space in jfac and kfac to expand from 2 to 4.
132c
133      noffset = ndens/2
134      noffset2 = noffset*2
135      noffset3 = noffset*3
136      do i = 1, noffset
137        jfac(noffset3+i) = zero
138        jfac(noffset2+i) = jfac(noffset+i)
139        jfac(noffset+i)  = zero
140c       jfac(i)          = already set correctly
141        kfac(noffset3+i) = kfac(noffset+i)
142        kfac(noffset2+i) = zero
143        kfac(noffset+i)  = kfac(i)
144        kfac(i)          = zero
145      enddo
146c
147c     Next expand the density.  We create a new ga and substitute it
148c     for the original.
149c
150      dims(1) = ndens * 2
151      dims(2) = nbf
152      dims(3) = nbf
153      chunk(1) = dims(1)
154      chunk(2) = -1
155      chunk(3) = -1
156      if (.not. nga_create (MT_DBL, 3, dims, 'Density',chunk,
157     &     g_dens2)) call errquit
158     &     ('uhf_fock_setup: could not allocate g_dens2',555, GA_ERR)
159c
160c     Copy "backwards" so that we don't loose any data
161c     Copy g_dens{a(1:ndens/2),b(1:ndens/2)}
162c     --> g_dens2{a(1:ndens/2),a(1:ndens/2),b(1:ndens/2),b(1:ndens/2)}
163c
164      alo(1) = noffset + 1 ! Beta
165      ahi(1) = ndens
166      alo(2) = 1
167      ahi(2) = nbf
168      alo(3) = 1
169      ahi(3) = nbf
170      blo(1) = noffset*3 + 1
171      bhi(1) = noffset*4
172      blo(2) = 1
173      bhi(2) = nbf
174      blo(3) = 1
175      bhi(3) = nbf
176      call nga_copy_patch('N',g_dens,alo,ahi,g_dens2,blo,bhi)
177      blo(1) = noffset*2 + 1
178      bhi(1) = noffset*3
179      call nga_copy_patch('N',g_dens,alo,ahi,g_dens2,blo,bhi)
180      alo(1) = 1       ! Alpha
181      ahi(1) = noffset
182      blo(1) = noffset + 1
183      bhi(1) = ndens
184      call nga_copy_patch('N',g_dens,alo,ahi,g_dens2,blo,bhi)
185      blo(1) = 1
186      bhi(1) = noffset
187      call nga_copy_patch('N',g_dens,alo,ahi,g_dens2,blo,bhi)
188c
189c     Now get rid of original density
190c
191      if (.not. ga_destroy(g_dens)) call errquit
192     *   ('uhf_fock_setup: failed to free g_dens', ndens, GA_ERR)
193c
194c     Assign g_dens to be the new matrix and update ndens
195c
196      g_dens = g_dens2
197      ndens = ndens * 2
198c
199c     Now expand the fock matrices doing what we did above for
200c     the densities.
201c     Copy g_fock{a(1:nfock/2),b(1:nfock/2)}
202c     --> g_fock2{b(1:nfock/2),b(1:nfock/2),a(1:nfock/2),a(1:nfock/2)}
203c
204      dims(1) = nfock * 2
205      chunk(1) = dims(1)
206      if (.not. nga_create (MT_DBL, 3, dims, 'Fock matrices',
207     &    chunk, g_fock2)) call errquit
208     &     ('uhf_fock_setup: could not allocate g_fock2',555, GA_ERR)
209c
210      noffset = nfock/2
211      alo(1) = 1             ! Alpha input Fock matrices
212      ahi(1) = noffset
213      blo(1) = noffset*3 + 1 ! Beta exchange matrices
214      bhi(1) = noffset*4
215      call nga_copy_patch('N',g_fock,alo,ahi,g_fock2,blo,bhi)
216      blo(1) = noffset*2 + 1 ! Beta Coulomb matrices
217      bhi(1) = noffset*3
218      call nga_copy_patch('N',g_fock,alo,ahi,g_fock2,blo,bhi)
219      alo(1) = noffset + 1   ! Beta input Fock matrices
220      ahi(1) = nfock
221      blo(1) = noffset + 1   ! Alpha exchange matrices
222      bhi(1) = nfock
223      call nga_copy_patch('N',g_fock,alo,ahi,g_fock2,blo,bhi)
224      blo(1) = 1             ! Alpha Coulomb matrices
225      bhi(1) = noffset
226      call nga_copy_patch('N',g_fock,alo,ahi,g_fock2,blo,bhi)
227c
228      if (.not. ga_destroy(g_fock)) call errquit
229     *   ('uhf_fock_setup: failed to free g_fock', nfock, GA_ERR)
230c
231      g_fock = g_fock2
232      nfock = nfock * 2
233c
234      if (oprint) then
235        if (ga_nodeid().eq.0) then
236          write(LuOut,*)'Density and Fock matrices after uhf_fock_setup'
237          call util_flush(LuOut)
238        endif
239        do i=1,ndens
240           ahi(1)=i
241           alo(1)=i
242           call nga_print_patch(g_dens,alo,ahi,0)
243        enddo
244        call ga_print(g_fock)
245      endif
246c
247      return
248      end
249c
250c******
251