1*
2* $Id$
3*
4      subroutine lcao_init_dn(ispin,ne,n2ft3d,dn,phi)
5      implicit none
6#include "errquit.fh"
7      integer ispin,ne(2)
8      integer n2ft3d
9      real*8 dn(n2ft3d,2)
10      complex*16 phi(*)
11
12#include "bafdecls.fh"
13#include "global.fh"
14
15*     **** local variables ****
16      logical value
17      integer i,k,nbasis,ms
18      integer nx,ny,nz
19      real*8  sum,scal,scal1,scal2,dv
20      real*8 dnscal(2)
21      integer tmp(2)
22
23
24*     ***** external functions ****
25      integer  aorbs_nbasis
26      real*8   aorbs_weight,lattice_omega,util_random
27      external aorbs_nbasis
28      external aorbs_weight,lattice_omega,util_random
29
30      dnscal(1) = dble(ne(1))/dble(ne(1)+ne(2))
31      dnscal(2) = dble(ne(2))/dble(ne(1)+ne(2))
32
33      call D3dB_nx(1,nx)
34      call D3dB_ny(1,ny)
35      call D3dB_nz(1,nz)
36      scal1 = 1.0d0/dble(nx*ny*nz)
37      scal2 = 1.0d0/lattice_omega()
38      dv    = scal1*lattice_omega()
39
40      value = BA_push_get(mt_dbl,n2ft3d,'tmp',tmp(2),tmp(1))
41      if (.not. value)
42     >  call errquit('lcao_init_dn:out of stack memory',0, MA_ERR)
43
44
45      call dcopy(ispin*n2ft3d,0.0d0,0,dn,1)
46      nbasis = aorbs_nbasis()
47      do i=1,nbasis
48
49*        **** get phi1 ****
50         call aorbs_aorb(i,phi)
51         call Pack_c_Copy(1,phi,dbl_mb(tmp(1)))
52         call Pack_c_unpack(1,dbl_mb(tmp(1)))
53         call D3dB_cr_pfft3b(1,1,dbl_mb(tmp(1)))
54
55        do ms=1,ispin
56           do k=1,n2ft3d
57              scal = aorbs_weight(i)*scal2*dnscal(ms)
58              dn(i,ms) = dn(i,ms) + scal*(dbl_mb(tmp(1)+k-1)**2)
59           end do
60        end do
61
62      end do
63
64
65*     **** randomize dn(*,2) if ispin=2 and ne(1)==ne(2) ****
66      if ((ispin.eq.2).and.(ne(1).eq.ne(2))) then
67           scal = util_random(9)
68           do k=1,n2ft3d
69              scal = (0.5d0-util_random(0))/dsqrt(dble(n2ft3d))
70              dn(i,2) = dn(i,2) + scal
71           end do
72      end if
73
74
75*     **** normalize densities ***
76      do ms=1,ispin
77        call D3dB_r_Zero_Ends(1,dn(1,ms))
78        call D3dB_r_dsum(1,dn(1,ms),sum)
79        sum = sum*dv
80        sum = dble(ne(ms))/sum
81c        call D3dB_r_SMul(1,sum,dn(1,ms),dn(1,ms))
82        call D3dB_r_SMul1(1,sum,dn(1,ms))
83      end do
84
85      value = BA_pop_stack(tmp(2))
86      if (.not. value)
87     >  call errquit('lcao_init_dn:popping stack memory',0, MA_ERR)
88      return
89      end
90
91
92