1!
2! Copyright (C) 2004 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8#undef DEBUG
9!---------------------------------------------------------------
10subroutine c6_dft (mesh, zed, grid)
11   !--------------------------------------------------------------------
12   !
13   use kinds,      only : DP
14   use constants,  only : e2, pi, fpi, BOHR_RADIUS_ANGS
15   use ld1inc,     only : lsd, nwf, oc, nn, ll, isw, psi, enl, vpot,vxt,vh, &
16                          enne, latt, rho
17   use radial_grids, only: radial_grid_type, ndmx
18   !
19   implicit none
20   !
21   ! I/O variables
22   !
23   type(radial_grid_type), intent(in) :: grid
24   integer mesh , mesh_save
25   real(DP) :: zed
26   !
27   ! local variables
28   !
29   logical :: csi, l_add_tf_term
30   real(DP) :: vnew(ndmx,2), rhoc1(ndmx), ze2, fac, vme(ndmx)
31   real(DP) :: rho_save(ndmx,2)
32   real(DP) :: error, error2, e, charge, beta, u, alpha, dalpha, c6, du1, &
33               du2, factor, thresh
34   real(DP), allocatable :: y(:), yy(:), sqr(:)
35   real(DP), allocatable :: dvpot(:), dvscf(:), drho(:), dvhx(:), dvxc(:), pp(:)
36   complex(DP), allocatable :: dy(:), drho_old(:)
37   integer i, is, n, l, iu, Nu, Nc, counter, nstop, nerr
38
39   allocate ( y(mesh),yy(mesh),sqr(mesh) )
40   allocate ( dvpot(mesh),dvscf(mesh),drho(mesh),dvhx(mesh),dvxc(mesh),pp(mesh) )
41   allocate ( dy(mesh), drho_old(mesh) )
42   !
43   write(6,'(/,/,/,5x,20(''-''),'' Compute C6 from polarizability.'',10(''-''),/)')
44   !
45   if (mesh.ne.grid%mesh) call errore('c6_dft',' mesh dimension is not as expected',1)
46   counter = 1
47   do i = 1, mesh
48      if (rho(i,1) .gt. 1.0d-30) counter = counter + 1
49   enddo
50   mesh_save = mesh
51   mesh = counter
52
53   if (lsd .ne. 0) call errore ('c6_dft', 'implemented only for non-magnetic ions', lsd)
54   csi = .true.
55   do i = 1, nwf
56      csi = csi .and. ( ((ll(i).eq.0) .and. (oc(i).eq.2 )) .or. &
57                        ((ll(i).eq.1) .and. (oc(i).eq.6 )) .or. &
58                        ((ll(i).eq.2) .and. (oc(i).eq.10)) .or. &
59                        ((ll(i).eq.3) .and. (oc(i).eq.14)) )
60   enddo
61   if (.not. csi) call errore ('c6_dft', 'implemented only for closed-shell ions', 1)
62!
63
64   n = 1
65   l = 0
66   e = -1.d-7
67   charge=0.d0
68   ze2 = - zed * e2
69   thresh = 1.d-10
70!
71   rho_save =  rho
72   rho=0.0_dp
73   do n=1,nwf
74      do i=1,mesh
75         rho(i,isw(n))=rho(i,isw(n))+oc(n)*(psi(i,1,n)**2+psi(i,2,n)**2)
76      enddo
77   enddo
78
79   error = 0.d0
80   do i=1,mesh
81      error = error + abs( rho(i,1)-rho_save(i,1) ) * grid%r2(i) * grid%dx
82      error = error + abs( rho(i,2)-rho_save(i,2) ) * grid%r2(i) * grid%dx
83   end do
84
85   if (error > 1.d-8) then
86      write (*,*) error
87      call errore('c6_dft','charge density rho from last vnew is inaccurate',1)
88   end if
89
90   rhoc1=0.d0
91   call new_potential(ndmx,mesh,grid,zed,vxt,lsd,.false.,latt,enne,rhoc1,rho,vh,vnew,0)
92   error = 0.d0
93   do i=1,mesh
94      error = error + abs( vpot(i,1)-vnew(i,1) ) * grid%r2(i) * grid%dx
95      error = error + abs( vpot(i,2)-vnew(i,2) ) * grid%r2(i) * grid%dx
96   end do
97   write (*,*) "Vpot-Vnew", error
98
99   nerr = 0
100   do n=1,nwf
101      if (oc(n) >= 0.0_dp) then
102         is=isw(n)
103         call ascheq (nn(n),ll(n),enl(n),mesh,grid,vnew(1,is),ze2,&
104                      thresh,psi(1,1,n),nstop)
105         nerr = nerr + nstop
106         write (*,'(4i3,2f10.5,i5)') n, nn(n),ll(n),isw(n),oc(n),enl(n),nstop
107      else
108         enl(n)=0.0_dp
109         psi(:,:,n)=0.0_dp
110      end if
111   end do
112
113!  from now on rho is the REAL rho w/o the volume element
114   do i = 1, mesh
115      rho(i,1) = rho(i,1) / (fpi*grid%r(i)**2)
116   end do
117!
118! initialize external perturbation (electric field)
119!
120   call init_dpot(grid%r,mesh,dvpot)
121!
122! derivative of xc-potential
123!
124   call dvxc_dn(mesh, rho, dvxc)
125!
126!write(*,'(1PE20.12)')sum(abs(dvxc))
127!stop
128   write(6,'(5x,''Frequency dependent polarizability is written into freq-pol.dat'',/)')
129
130   c6    = 0.0d0
131   alpha = 0.0d0
132
133   open(1, file = 'freq-pol-dft.dat')
134   write (1,'(15x,"    u          alpha(angstrom)       alpha(a.u.)  ",/)')
135   !
136   Nu  = 230
137   Nc  = 50
138   du1 = 0.1d0
139   du2 = 0.25d0
140   u   = -du1
141   !
142   do iu=0,Nu
143      !
144      if (iu .le. 50) then
145         u = u + du1
146      else
147         u = u + du2
148      endif
149      !
150      if (iu.eq.0) then
151         do i=1,mesh
152            dvscf(i) = dvpot(i)
153            drho_old(i) = 0.d0
154         end do
155      end if
156      beta = 0.05
157      dalpha = 1.0d+99
158      alpha = 0.d0
159      counter = 0
160      do while (dalpha > 1.d-9)
161         counter =  counter + 1
162         !
163         ! solve Sternheimer equation for the auxiliary wavefunction
164         !
165         drho = 0.d0
166         do n=1,nwf
167            do l = 1 + ll(n), max( 1 - ll(n), 0 ), - 2
168!               write (*,*) l, ll(n)
169               y(1:mesh) = psi(1:mesh,1,n)/grid%sqr(1:mesh)
170               vme(:) = vnew(:,isw(n)) - enl(n)
171               call sternheimer(u,l,ll(n),mesh,grid%dx,grid%r,grid%sqr,grid%r2,vme,zed,y,dvscf,dy)
172               fac = 2.0d0 * (2.d0 * ll(n) + 1.d0 )
173               if (ll(n)==1 .and. l==2) fac = fac * 2.d0/3.d0
174               if (ll(n)==1 .and. l==0) fac = fac * 1.d0/3.d0
175               if (ll(n)==2 .and. l==3) fac = fac * 3.d0/5.d0
176               if (ll(n)==2 .and. l==1) fac = fac * 2.d0/5.d0
177               call inc_drho_of_r(mesh, grid%dx, grid%r, grid%r2, y, dy, fac, drho)
178#ifdef DEBUG
179         write (*,*) "========================", n, l
180         write (*,*) "y(1:3)"
181         write (*,*) y(1:3)
182         write (*,*) "dy(1:3)"
183         write (*,*) dy(1:3)
184         write (*,*) "drho(1:3)"
185         write (*,*) drho(1:3)
186#endif
187
188            end do
189         end do
190#ifdef DEBUG
191         write (*,*) "========================"
192         write (*,*) "drho(1:3)"
193         write (*,*) drho(1:3)
194         write (*,*) "drho(20:22)"
195         write (*,*) drho(20:22)
196         write (*,*) "drho(40:42)"
197         write (*,*) drho(40:42)
198         write (*,*) "drho(mesh-2:mesh)"
199         write (*,*) drho(mesh-2:mesh)
200#endif
201         !
202         ! compute dv of drho (w/o the TF term)
203         !
204         l_add_tf_term = .false.
205         call dv_of_drho(mesh, grid%dx, grid%r,grid%r2,rho,drho,dvhx,dvxc,pp,l_add_tf_term)
206
207#ifdef DEBUG
208         write (*,*) "========================"
209         write (*,*) "dvhx(1:3)"
210         write (*,*) dvhx(1:3)
211         write (*,*) "dvhx(20:22)"
212         write (*,*) dvhx(20:22)
213         write (*,*) "dvhx(40:42)"
214         write (*,*) dvhx(40:42)
215         write (*,*) "dvhx(mesh-2:mesh)"
216         write (*,*) dvhx(mesh-2:mesh)
217         write (*,*) "========================"
218         write (*,*) "pp(1:3)"
219         write (*,*) pp(1:3)
220         write (*,*) "pp(20:22)"
221         write (*,*) pp(20:22)
222         write (*,*) "pp(40:42)"
223         write (*,*) pp(40:42)
224         write (*,*) "pp(mesh-2:mesh)"
225         write (*,*) pp(mesh-2:mesh)
226#endif
227         !
228         ! mix
229         !
230         error = 0.d0
231         error2 = 0.d0
232         do i=1,mesh
233            dvscf(i) = dvscf(i) + beta * (dvpot(i)+dvhx(i) -dvscf(i))
234            error = error + abs (drho(i) -drho_old(i))
235            error2 = error2 + abs (drho(i) -drho_old(i))* grid%r(i) * grid%dx
236            drho_old(i) = drho(i)
237         end do
238         dalpha = abs(alpha + pp(mesh))
239         alpha = -pp(mesh)
240!         write (*,'(4e16.6)') alpha, dalpha, error, error2
241
242      end do
243
244      write (1,'(17x, f8.4, 3x, 1PE14.6, 9x, 1PE14.6)') u, pp(mesh)*BOHR_RADIUS_ANGS**3, pp(mesh)
245      if (iu .eq. 0) &
246      write (6,'(5x, "Static polarizability: ", f10.5, " (in angstrom^3)   --->", f10.5,&
247                & "  (in e^2a0^3)")') pp(mesh)*BOHR_RADIUS_ANGS**3, pp(mesh)
248
249      if (iu .eq. 0)                  factor = 0.5d0 * du1
250      if (iu .gt. 0 .and. iu .lt. Nc) factor = du1
251      if (iu .eq. Nc)                 factor = 0.5d0 * ( du1 + du2)
252      if (iu .gt. Nc .and. iu .lt. Nu) factor = du2
253      if (iu .eq. Nu)                 factor = 0.5d0 * du2
254      c6 = c6 + factor*alpha*alpha
255
256   end do
257
258   c6 = c6 * 3.d0 / pi
259
260   write (*,'(/, 5x, a, f12.6)') "C6 coefficient in units [e2*a0**5]", c6/e2
261   !
262   write(6,'(/,5x,20(''-''),'' End of C6 calculation '',20(''-''),/)')
263
264   deallocate ( dy )
265   deallocate ( y, yy, sqr )
266   deallocate ( dvpot, dvscf, drho, dvhx, pp )
267
268   return
269end subroutine c6_dft
270
271!--------------------------------------------------------------------
272subroutine inc_drho_of_r(mesh, dx, r, r2, y, dy, fac, drho)
273   !--------------------------------------------------------------------
274   ! compute the first order variation of the density from
275   ! the zeroth and first order auxiliary wavefunctions y and dy
276   !
277   use constants, only : e2, pi, fpi
278   implicit none
279   !
280   ! I/O vaiables
281   !
282   integer mesh
283   real (kind=8) :: dx, fac, r(mesh), r2(mesh), y(mesh), drho(mesh)
284   complex (kind=8) :: dy(mesh)
285   ! local variables
286   integer i
287
288   do i=1,mesh
289      drho(i) = drho(i) + fac * 2.d0 * y(i) * real(dy(i)) * r(i) / (fpi*r2(i))
290   end do
291
292   return
293end subroutine
294
295