1      Subroutine localize_sic(g_movecs, evals)
2c
3C$Id$
4c
5      implicit none
6#include "errquit.fh"
7c
8c     integer noc(2)
9      integer isp, x, me, i_degen, l_degen, n_levels,
10     &        evals(2), ij, ik, n_orbitals, i_level, i_orb, j,
11     &        aux_levels
12      double precision start
13
14      integer g_movecs(2)
15      integer g_uc(4),g_sc,g_over,l_sc,k_sc,l_c,k_c,nocc
16c     integer ncore,nloc,liloc,iiloc,i
17      integer nloc,liloc,iiloc
18      logical geom_num_core
19      external geom_num_core
20c
21#include "mafdecls.fh"
22#include "global.fh"
23#include "cdft.fh"
24#include "stdio.fh"
25#include "util.fh"
26c
27      integer  ga_create_atom_blocked
28      external ga_create_atom_blocked
29c
30      logical oprint_sic
31c
32c******************************************************************************
33c
34c
35      me=ga_nodeid()
36      oprint_sic = util_print('SIC information', print_high)
37c
38      if (.not.ma_push_get
39     &   (MT_Int,nbf,'map of orb to loc',liloc,iiloc))
40     &   call errquit('localize_sic:push_get failed', 911, MA_ERR)
41c
42      do isp=1,ipol
43c
44         nocc=noc(isp)
45         if (.not.MA_Push_Get(MT_Int, nocc, 'levels',
46     &        l_degen, i_degen))
47     &        call errquit('localize_sic: cannot allocate levels',0,
48     &       MA_ERR)
49         if (nocc.ne.0) then
50           aux_levels = 1
51           n_orbitals = 1
52           Int_MB(i_degen + aux_levels - 1) = n_orbitals
53           do j = noc(isp),2,-1
54             start = (dbl_mb(evals(isp) + j - 1) -
55     &                dbl_mb(evals(isp) + j - 2))
56             if (start.le.1.0d-04) then
57               n_orbitals = n_orbitals + 1
58             else
59               n_orbitals = 1
60               aux_levels = aux_levels + 1
61             end if
62             Int_MB(i_degen + aux_levels - 1) = n_orbitals
63           end do
64           n_levels = aux_levels
65         else
66           n_levels = 0
67         end if
68c
69         if (nocc.gt.0) then
70           do x= 1, 4
71             if (.not. ga_create(MT_DBL, nbf, nbf, 'uc',
72     $             nbf, 1, g_uc(x))) call errquit('g_uc_loc_sic',x,
73     &       GA_ERR)
74           end do
75c
76           call int_dip_ga(AO_bas_han, AO_bas_han, g_uc(1), g_uc(2),
77     &                     g_uc(3))
78c
79           if (.not. ga_create(MT_DBL, nbf, nbf, 'sc',
80     $          nbf, 1, g_sc)) call errquit('g_sc_loc_sic',0, GA_ERR)
81           do x = 1, 3
82             call ga_dgemm('n', 'n', nbf, nbf, nbf,
83     $                     1.0d0, g_uc(x), g_movecs(isp), 0.0d0, g_sc)
84             call ga_copy_patch('n',g_sc,1,nbf,1,nbf,g_uc(x),1,
85     &                           nbf,1,nbf)
86           end do
87           g_over  = ga_create_atom_blocked(geom, AO_bas_han,
88     &                                    'g_over_loc_sic')
89           call ga_zero(g_over)
90           call int_1e_ga(AO_bas_han, AO_bas_han, g_over, 'overlap',
91     &                    .false.)
92           call ga_dgemm('n', 'n', nbf, nbf, nbf,
93     $                   1.0d0, g_over, g_movecs(isp), 0.0d0, g_uc(4))
94           if (.not. ma_push_get(mt_dbl, 8*nbf, 'sc', l_sc, k_sc))
95     $          call errquit('ma for sc', 0, MA_ERR)
96           if (.not. ma_push_get(mt_dbl, 8*nbf, 'c', l_c, k_c))
97     $          call errquit('ma for c', 0, MA_ERR)
98c
99c     Localize the occupied orbitals
100c
101c          if (.not. geom_num_core(geom,ncore)) ncore = 0
102c          if (ncore .gt. 0) then
103c             do i = 1, ncore
104c                int_mb(iiloc+i-1) = i
105c             end do
106c             nloc = ncore
107c             call localizeFB(AO_bas_han, dbl_mb(k_c), dbl_mb(k_sc),
108c    $                        nloc, int_mb(iiloc), nbf, nbf,
109c    $                        g_movecs(isp), g_uc)
110c          end if
111c
112c     Localized the occupied
113c
114c          do i = ncore+1, nocc
115c             int_mb(iiloc + i - ncore - 1) = i
116c          end do
117c          nloc = nocc - ncore
118c
119c          ik = noc(isp) -  Int_MB(i_degen) + 1   ! If HOMO is not involved
120           ik = noc(isp) + 1
121           ij = 0
122c          do 100 i_level = 2,n_levels            ! If HOMO is not involved
123           do 100 i_level = 1,n_levels
124             n_orbitals = Int_MB(i_degen + i_level - 1)
125             do 200 i_orb = 1,n_orbitals
126               ik = ik - 1
127               ij = ij + 1
128               int_mb(iiloc + ij - 1) = ik
129  200        continue
130  100      continue
131           nloc = ij
132           if (nloc.gt.0) then
133             if (me.eq.0)
134     $        call util_print_centered(6,
135     $    'Foster-Boys orbital localization for the SIC approximation',
136     &           40, .true.)
137              if (me.eq.0.and.oprint_sic) write(LuOut,*)
138     $                           ' Number of Orbitals to localize',nloc
139              call localizeFB(AO_bas_han, dbl_mb(k_c), dbl_mb(k_sc),
140     $                        nloc, int_mb(iiloc), nbf, nbf,
141     $                        g_movecs(isp), g_uc)
142           end if
143c
144           do x=1,4
145             if (.not. ga_destroy(g_uc(x))) call errquit
146     &        ('localize_sic: could not destroy g_uc', 0, GA_ERR)
147           end do
148           if (.not. ga_destroy(g_sc)) call errquit
149     &        ('localize_sic: could not destroy g_sc', 0, GA_ERR)
150           if (.not. ga_destroy(g_over)) call errquit
151     &        ('localize_sic: could not destroy g_sc', 0, GA_ERR)
152           if (.not.ma_pop_stack(l_c))
153     &        call errquit('localize_sic: cannot pop stack',0, MA_ERR)
154           if (.not.ma_pop_stack(l_sc))
155     &        call errquit('localize_sic: cannot pop stack',0, MA_ERR)
156          end if
157          if (.not.ma_pop_stack(l_degen))
158     &              call errquit('localize_sic: cannot pop stack',0,
159     &       MA_ERR)
160        end do
161        if (.not.ma_pop_stack(liloc))
162     &        call errquit('localize_sic: cannot pop stack',0,
163     &       MA_ERR)
164        return
165      end
166
167