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