1c $Id$
2
3*     *********************************
4*     *                               *
5*     *         nwpw_hartree_init     *
6*     *                               *
7*     *********************************
8      subroutine nwpw_hartree_init(nkatm,
9     >                        nprj,nbasis,psp_type,lmax,
10     >                        nprj_max0,l_prj,m_prj,b_prj,
11     <                        hartree_tag)
12      implicit none
13      integer nkatm
14      integer nprj(*),nbasis(*),psp_type(*),lmax(*)
15      integer nprj_max0
16      integer l_prj(nprj_max0,*),m_prj(nprj_max0,*),b_prj(nprj_max0,*)
17      integer hartree_tag
18
19
20#include "bafdecls.fh"
21#include "errquit.fh"
22#include "nwpw_hartree.fh"
23
24c     **** local variables ****
25      logical ok
26      integer ii,ia,nsize
27      integer ic
28      integer l,m,li,mi,lj,mj,bi,bj,iprj,jprj,lm
29      integer li1,mi1,lj1,mj1,bi1,bj1,iprj1,jprj1
30      integer i_p,i_t
31      integer i_tlm,i_plm,i_rlm
32      integer matr_ptr
33      integer nb
34      real*8  tmp_theta,cs_theta,tmp_phi,tmp_gaunt
35
36c     **** external functions ****
37      integer  psi_data_get_chnk
38      external psi_data_get_chnk
39      real*8   rtheta_lm,drtheta_lm,nwpw_gaunt
40      external rtheta_lm,drtheta_lm,nwpw_gaunt
41
42      call nwpw_timing_start(4)
43      ok = .true.
44
45*      ***** set up index arrays for nwpw_hartree_solve2 *****
46      ok = BA_alloc_get(mt_int,nkatm,"nindx_hartree",
47     >                  nindx_hartree(2),nindx_hartree(1))
48      ok = ok.and.
49     >     BA_alloc_get(mt_int,nkatm,"shift_hartree",
50     >                  shift_hartree(2),shift_hartree(1))
51      if (.not.ok)
52     >call errquit("nwpw_hartree_init:error allocating work arrays",0,0)
53
54      nsize = 0
55      do ia=1,nkatm
56         nb = nbasis(ia)
57         matr_ptr = psi_data_get_chnk(hartree_tag(ia))
58         int_mb(shift_hartree(1)+ia-1) = nsize
59         do jprj = 1,nprj(ia)
60            lj = l_prj(jprj,ia)
61            mj = m_prj(jprj,ia)
62            bj = b_prj(jprj,ia)
63            do iprj = 1,nprj(ia)
64               li = l_prj(iprj,ia)
65               mi = m_prj(iprj,ia)
66               bi = b_prj(iprj,ia)
67               do jprj1 = 1,nprj(ia)
68                  lj1 = l_prj(jprj1,ia)
69                  mj1 = m_prj(jprj1,ia)
70                  bj1 = b_prj(jprj1,ia)
71                  do iprj1 = 1,nprj(ia)
72                     li1 = l_prj(iprj1,ia)
73                     mi1 = m_prj(iprj1,ia)
74                     bi1 = b_prj(iprj1,ia)
75                     do l=0,2*lmax(ia)
76                        do m=-l,l
77                           tmp_gaunt =
78     >                        nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
79     >                       *nwpw_gaunt(.false.,l,m,li1,mi1,lj1,mj1)
80     >                       *dbl_mb(matr_ptr
81     >                              + (bi-1)
82     >                              + (bj-1)*nb
83     >                              + (bi1-1)*nb*nb
84     >                              + (bj1-1)*nb*nb*nb
85     >                              + l*nb*nb*nb*nb)
86                           if (dabs(tmp_gaunt).gt.1.0d-11) then
87                              nsize = nsize + 1
88                           end if
89                        end do
90                     end do
91
92                  end do
93               end do
94
95            end do
96         end do
97         int_mb(nindx_hartree(1)+ia-1)= nsize
98     >                                - int_mb(shift_hartree(1)+ia-1)
99      end do
100
101      ok = BA_alloc_get(mt_int,nsize,"iprj_hartree",
102     >                 iprj_hartree(2),iprj_hartree(1))
103      ok = ok.and.
104     >     BA_alloc_get(mt_int,nsize,"jprj_hartree",
105     >                 jprj_hartree(2),jprj_hartree(1))
106      ok = ok.and.
107     >     BA_alloc_get(mt_int,nsize,"iprj1_hartree",
108     >                 iprj1_hartree(2),iprj1_hartree(1))
109      ok = ok.and.
110     >     BA_alloc_get(mt_int,nsize,"jprj1_hartree",
111     >                 jprj1_hartree(2),jprj1_hartree(1))
112      ok = ok.and.
113     >     BA_alloc_get(mt_dbl,nsize,"coeff_hartree",
114     >                 coeff_hartree(2),coeff_hartree(1))
115      if (.not.ok)
116     >call errquit("nwpw_hartree_init:error allocating work arrays",0,0)
117
118
119      nsize = 0
120      do ia=1,nkatm
121         nb = nbasis(ia)
122         matr_ptr = psi_data_get_chnk(hartree_tag(ia))
123         do jprj = 1,nprj(ia)
124            lj = l_prj(jprj,ia)
125            mj = m_prj(jprj,ia)
126            bj = b_prj(jprj,ia)
127            do iprj = 1,nprj(ia)
128               li = l_prj(iprj,ia)
129               mi = m_prj(iprj,ia)
130               bi = b_prj(iprj,ia)
131               do jprj1 = 1,nprj(ia)
132                  lj1 = l_prj(jprj1,ia)
133                  mj1 = m_prj(jprj1,ia)
134                  bj1 = b_prj(jprj1,ia)
135                  do iprj1 = 1,nprj(ia)
136                     li1 = l_prj(iprj1,ia)
137                     mi1 = m_prj(iprj1,ia)
138                     bi1 = b_prj(iprj1,ia)
139                     do l=0,2*lmax(ia)
140                        do m=-l,l
141                           tmp_gaunt =
142     >                        nwpw_gaunt(.false.,l,m,li,mi,lj,mj)
143     >                       *nwpw_gaunt(.false.,l,m,li1,mi1,lj1,mj1)
144     >                       *dbl_mb(matr_ptr
145     >                              + (bi-1)
146     >                              + (bj-1)*nb
147     >                              + (bi1-1)*nb*nb
148     >                              + (bj1-1)*nb*nb*nb
149     >                              + l*nb*nb*nb*nb)
150                           if (dabs(tmp_gaunt).gt.1.0d-11) then
151                              int_mb(iprj_hartree(1)+nsize)  = iprj
152                              int_mb(jprj_hartree(1)+nsize)  = jprj
153                              int_mb(iprj1_hartree(1)+nsize) = iprj1
154                              int_mb(jprj1_hartree(1)+nsize) = jprj1
155                              nsize = nsize + 1
156                           end if
157                        end do
158                     end do
159
160                  end do
161               end do
162
163            end do
164         end do
165      end do
166
167      call nwpw_timing_end(4)
168      return
169      end
170
171*     ********************************************
172*     *                                          *
173*     *             nwpw_hartree_solve           *
174*     *                                          *
175*     ********************************************
176      subroutine nwpw_hartree_solve(ia,
177     >                              ispin,ne,nprj,sw1,sw2)
178      implicit none
179      integer ii,ia
180      integer ispin,ne(2),nprj
181      real*8  sw1(ne(1)+ne(2),nprj)
182      real*8  sw2(ne(1)+ne(2),nprj)
183
184#include "bafdecls.fh"
185#include "errquit.fh"
186#include "nwpw_hartree.fh"
187
188      integer shift
189
190      call nwpw_timing_start(4)
191      call nwpw_timing_start(22)
192
193*     **** hartree potential  non-local matrix elements ****
194      shift = int_mb(shift_hartree(1)+ia-1)
195      call nwpw_hartree_sw1tosw2(ispin,ne,nprj,sw1,sw2,
196     >                           int_mb(nindx_hartree(1)+ia-1),
197     >                           int_mb(iprj_hartree(1)+shift),
198     >                           int_mb(jprj_hartree(1)+shift),
199     >                           int_mb(iprj1_hartree(1)+shift),
200     >                           int_mb(jprj1_hartree(1)+shift),
201     >                           dbl_mb(coeff_hartree(1)+shift))
202
203      call nwpw_timing_end(4)
204      call nwpw_timing_end(22)
205
206      return
207      end
208
209
210*     *****************************************
211*     *                                       *
212*     *             nwpw_hartree_end          *
213*     *                                       *
214*     *****************************************
215      subroutine nwpw_xc_end()
216      implicit none
217
218#include "bafdecls.fh"
219#include "errquit.fh"
220#include "nwpw_hartree.fh"
221
222      logical ok
223
224      call nwpw_timing_start(4)
225      ok = .true.
226      ok = ok.and.BA_free_heap(coeff_hartree(2))
227      ok = ok.and.BA_free_heap(jprj_hartree(2))
228      ok = ok.and.BA_free_heap(iprj_hartree(2))
229      ok = ok.and.BA_free_heap(jprj1_hartree(2))
230      ok = ok.and.BA_free_heap(iprj1_hartree(2))
231      ok = ok.and.BA_free_heap(shift_hartree(2))
232      ok = ok.and.BA_free_heap(nindx_hartree(2))
233
234      if (.not.ok)
235     > call errquit("nwpw_hartree_end: error freeing heap",0,MA_ERR)
236
237      call nwpw_timing_end(4)
238      return
239      end
240
241c     *********************************************
242c     *                                           *
243c     *           nwpw_hartree_sw1tosw2           *
244c     *                                           *
245c     *********************************************
246      subroutine nwpw_hartree_sw1tosw2(ispin,ne,nprj,sw1,sw2,
247     >                                 nindx,
248     >                                 iprj_hartree, jprj_hartree,
249     >                                 iprj1_hartree,jprj1_hartree,
250     >                                 coeff_hartree)
251      implicit none
252      integer ispin,ne(2),nprj
253      real*8  sw1(ne(1)+ne(2),nprj)
254      real*8  sw2(ne(1)+ne(2),nprj)
255      integer n1dgrid,nbasis,lmax2
256      integer nindx,iprj_hartree(*),jprj_hartree(*)
257      integer iprj1_hartree(*),jprj1_hartree(*)
258      real*8  coeff_hartree(*)
259
260      integer n,i,iprj,jprj,iprj1,jprj1
261      real*8  coeff,w,scal
262
263*     **** external functions ****
264      real*8   lattice_omega
265      external lattice_omega
266
267      call nwpw_timing_start(21)
268      scal = 1.0d0/dsqrt(lattice_omega())
269      scal = scal/lattice_omega()
270
271*     ***init to zero***
272      do i=1,nindx
273         iprj  = iprj_hartree(i)
274         jprj  = jprj_hartree(i)
275         iprj1 = iprj1_hartree(i)
276         jprj1 = jprj1_hartree(i)
277         coeff = coeff_hartree(i)*scal
278         w = 0.0d0
279         do n=1,ne(1)+ne(2)
280            w = w + sw1(n,iprj1)*sw1(n,jprj1)
281         end do
282         do n=1,ne(1)+ne(2)
283            sw2(n,iprj) = sw2(n,iprj) + coeff*w*sw1(n,jprj)
284         end do
285      end do
286      call nwpw_timing_end(21)
287      return
288      end
289