1!---------------------------------------------------------------------
2subroutine dfx_new(dchi0, vx)
3!---------------------------------------------------------------------
4! OEP calculation of the exchange potential
5! --------------------------------------------------------------------
6#undef DEBUG
7
8   use constants, only: pi
9   use kinds,     only: DP
10   use ld1_parameters, only : nwfx
11   use ld1inc,    only: nwf, nspin, oc, rho, psi, isw, grid
12   use radial_grids, only: ndmx, hartree
13   implicit none
14   !
15   ! I/O variables
16   !
17   real(DP) :: dchi0(ndmx,nwfx)
18   real(DP),intent(out) ::  vx(ndmx,2) ! the oep exchange only potential
19                                       ! for the two spin configurations
20
21   !
22   ! local variables
23   !
24   integer, parameter :: niterx = 12 ! 6, 12, 24
25   real(DP) :: drho0(ndmx,2),appchim1(ndmx,2), vslater(ndmx,2)
26   real(DP) :: int_0_inf_dr
27   real(DP) :: vvx(ndmx,2,niterx),drhox(ndmx,2,niterx), dvh(ndmx,2,niterx), &
28                    dvh0(ndmx,2), aux(ndmx), drho1(ndmx,2), dvh1(ndmx,2)
29   real(DP) :: a(niterx,niterx), &
30               inva(niterx,niterx), & ! inverse of a, i guess
31               b1(niterx), b2(niterx), c, c1, work(niterx), x(niterx), uno
32   integer :: iwork(niterx), info, iterx
33   integer :: i, j, jter, k, nu, is
34   real(DP) :: third, fac, capel
35   logical :: first = .true.
36   save first
37
38!-set the left hand side of the equation
39   call drho0ofvx(drho0,dchi0)
40
41   do is=1,nspin
42      call hartree(0,2,grid%mesh,grid,drho0(1,is),dvh0(1,is))
43   end do
44
45   aux(1:grid%mesh) = drho0(1:grid%mesh,1)*dvh0(1:grid%mesh,1)
46   if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
47                               drho0(1:grid%mesh,2)*dvh0(1:grid%mesh,2) )
48   c = int_0_inf_dr(aux,grid,grid%mesh,2)
49
50   if (first) then
51!      first = .false.
52!-slater exchange potential
53      do is=1,nspin
54         vslater(1:grid%mesh,is) = 0.d0
55         do nu=1,nwf
56            if (isw(nu) == is) vslater(1:grid%mesh,is) = vslater(1:grid%mesh,is) + &
57                                    oc(nu)*psi(1:grid%mesh,1,nu)*dchi0(1:grid%mesh,nu)
58         end do
59         vslater(1:grid%mesh,is) = vslater(1:grid%mesh,is) / rho(1:grid%mesh,is)
60!         write (*,*) vslater(1:grid%mesh,is)
61       end do
62   else
63      vslater(1:grid%mesh,1:nspin) = vx(1:grid%mesh,1:nspin)
64   end if
65!- is a reasonable starting guess
66   call drhoofv(drho1,vslater)
67
68   drho1(1:grid%mesh,1:nspin) = drho1(1:grid%mesh,1:nspin) - drho0(1:grid%mesh,1:nspin)
69   do is=1,nspin
70      call hartree(0,2,grid%mesh,grid,drho1(1,is),dvh1(1,is))
71   end do
72
73   aux(1:grid%mesh) = drho1(1:grid%mesh,1) * dvh1(1:grid%mesh,1)
74   if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
75                               drho1(1:grid%mesh,2) * dvh1(1:grid%mesh,2) )
76   c1 = int_0_inf_dr(aux,grid,grid%mesh,2)
77
78#ifdef DEBUG
79   write (*,'(a,2f16.10)') "C, C1 ", c, c1
80#endif
81
82!- simple Thomas-Fermi approximation to \chi^-1
83   third = 1.0d0/3.0d0
84   fac   = -(0.75d0*pi)**(2.0d0/3.0d0)
85
86   do is =1, nspin
87      appchim1(1:grid%mesh,is) = fac/(grid%r(1:grid%mesh)* &
88                            (nspin*rho(1:grid%mesh,is)*grid%r(1:grid%mesh))**third)
89   end do
90
91   drhox(1:grid%mesh,1:nspin,1) = drho1(1:grid%mesh,1:nspin)
92
93   if (c1 < 1.d-12 ) then
94      vx(1:grid%mesh,1:nspin) = vslater(1:grid%mesh,1:nspin)
95      return
96   end if
97!- ITERATE !
98   do iterx =1,niterx
99!- set a new normalized correction vector vvx = chim1*drho/norm
100
101      vvx(1:grid%mesh,1:nspin,iterx) = appchim1(1:grid%mesh,1:nspin) * &
102                                             drhox(1:grid%mesh,1:nspin,iterx)
103      do is=1,nspin
104         call hartree(0,2,grid%mesh,grid,vvx(1,is,iterx),dvh(1,is,iterx))
105      end do
106
107      aux(1:grid%mesh) =vvx(1:grid%mesh,1,iterx) * dvh(1:grid%mesh,1,iterx)
108      if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) +  &
109                                  vvx(1:grid%mesh,2,iterx) * dvh(1:grid%mesh,2,iterx) )
110      capel = int_0_inf_dr(aux,grid,grid%mesh,2)
111#ifdef DEBUG
112      write (*,*) "norm ", capel
113#endif
114      if (capel > 0) then
115         capel = 1.d0/sqrt(capel)
116         vvx(1:grid%mesh,1:nspin,iterx) = vvx(1:grid%mesh,1:nspin,iterx) * capel
117      end if
118!- compute the corresponding drho
119      call drhoofv(drhox(1,1,iterx),vvx(1,1,iterx) )
120      do is =1,nspin
121         call hartree(0,2,grid%mesh,grid,drhox(1,is,iterx),dvh(1,is,iterx))
122      end do
123
124      aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx) * dvh1(1:grid%mesh,1)
125      if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
126                                  drhox(1:grid%mesh,2,iterx) * dvh1(1:grid%mesh,2) )
127      b1(iterx) = int_0_inf_dr(aux,grid,grid%mesh,2)
128
129      aux(1:grid%mesh) = drho1(1:grid%mesh,1) * dvh(1:grid%mesh,1,iterx)
130      if (nspin==2) aux(1:grid%mesh) =  2.d0 * ( aux(1:grid%mesh) + &
131                                   drho1(1:grid%mesh,2) * dvh(1:grid%mesh,2,iterx) )
132      b2(iterx) = int_0_inf_dr(aux,grid,grid%mesh,2)
133
134      do jter =1,iterx
135
136         aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx) * dvh(1:grid%mesh,1,jter)
137         if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
138                                     drhox(1:grid%mesh,2,iterx)*dvh(1:grid%mesh,2,jter) )
139
140         a(iterx,jter) = int_0_inf_dr(aux,grid,grid%mesh,2)
141
142         aux(1:grid%mesh) = drhox(1:grid%mesh,1,jter) * dvh(1:grid%mesh,1,iterx)
143         if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
144                                     drhox(1:grid%mesh,2,jter)*dvh(1:grid%mesh,2,iterx) )
145         a(jter,iterx) = int_0_inf_dr(aux,grid,grid%mesh,2)
146
147      end do
148
149      capel = 0.d0
150
151      do i=1,iterx
152         do j=1,iterx
153            capel = capel + abs(a(i,j)-a(j,i))
154            a(i,j) = 0.5d0*(a(i,j)+a(j,i))
155            a(j,i) = a(i,j)
156         end do
157      end do
158
159
160#ifdef DEBUG
161      write (*,'(a,12f16.10)') "CAPEL a  ", capel
162#endif
163
164      inva = a
165
166
167
168      CALL DSYTRF('U', iterx, inva, niterx, iwork, work, niterx,info)
169      if (info.ne.0) stop 'factorization'
170      CALL DSYTRI( 'U', iterx, inva, niterx, iwork, work, info )
171      if (info.ne.0) stop 'DSYTRI'
172      forall ( i=1:iterx, j=1:iterx, j>i ) inva(j,i) = inva(i,j)
173
174!      write (*,*) "INVA "
175!      write (*,'(12f16.10)') inva
176      capel = 0.d0
177      do i=1,iterx
178         do j=1,iterx
179            uno = 0.d0
180            do k=1,iterx
181               uno = uno + a(i,k)*inva(k,j)
182            end do
183            if (i.eq.j) uno = uno - 1.d0
184            capel = capel + abs(uno)
185         end do
186      end do
187#ifdef DEBUG
188      write (*,'(a,12f16.10)') "CAPEL uno", capel
189#endif
190
191      x = 0.d0
192      capel = c1
193      do i=1,iterx
194         do j=1,iterx
195            x(i) = x(i) - inva(i,j) * 0.5d0*(b1(j)+b2(j))
196         end do
197         capel = capel + x(i) * (b1(i)+b2(i))
198         do j =1,i
199            capel = capel + x(i)*a(i,j)*x(j)
200         end do
201         do j =1,i-1
202            capel = capel + x(i)*a(j,i)*x(j)
203         end do
204      end do
205!      write (*,'(a,12f16.10)') "X ", x
206#ifdef DEBUG
207      write (*,*) "capel       ", capel
208#endif
209
210
211      vx(1:grid%mesh,1:nspin) = vslater(1:grid%mesh,1:nspin)
212      do j=1,iterx
213         vx(1:grid%mesh,1:nspin) = vx(1:grid%mesh,1:nspin) + x(j) * vvx(1:grid%mesh,1:nspin,j)
214      end do
215
216      if (iterx.eq.niterx) return
217
218      call drhoofv(drhox(1,1,iterx+1), vx )
219      drhox(1:grid%mesh,1:nspin,iterx+1) = drhox(1:grid%mesh,1:nspin,iterx+1) - &
220                                      drho0(1:grid%mesh,1:nspin)
221      do is=1,nspin
222         call hartree(0,2,grid%mesh,grid,drhox(1,1,iterx+1),dvh(1,1,iterx+1))
223      end do
224
225      aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx+1)*dvh(1:grid%mesh,1,iterx+1)
226      if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + &
227                                  drhox(1:grid%mesh,1,iterx+1)*dvh(1:grid%mesh,1,iterx+1))
228      capel = int_0_inf_dr(aux,grid,grid%mesh,2)
229#ifdef DEBUG
230      write (*,*) "capel-check ", capel
231#endif
232!
233   end do
234!
235   return
236end subroutine dfx_new
237