1! 2! Copyright (C) 2002-2007 Quantum ESPRESSO groups 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 9 10!---------------------------------------------------------------------- 11SUBROUTINE vol_clu(rho_real,rho_g,flag) 12!---------------------------------------------------------------------- 13! it computes the volume of the cluster (cluster calculations) starting 14! from the measure of the region of space occupied by the electronic density 15! above a given threshold 16 17 USE kinds, ONLY: dp 18 USE constants, ONLY: pi 19 USE cell_base, ONLY: alat, at, h, omega, tpiba, tpiba2 20 USE electrons_base, ONLY: nspin 21 USE ions_base, ONLY: nsp, amass, nat, ityp 22 USE ions_positions, ONLY: tau0 23 USE gvect, ONLY: g, gg 24 USE cp_main_variables, only: drhor 25 USE control_flags, ONLY: tpre 26 USE fft_base, ONLY: dfftp 27 USE fft_interfaces, ONLY: invfft 28 USE pres_ai_mod, ONLY: rho_thr, n_cntr, cntr, step_rad, fill_vac, & 29 & delta_eps, delta_sigma, axis, & 30 & abisur, dthr, Surf_t, rho_gaus, v_vol, & 31 & posv, xc0, weight, volclu, stress_vol, & 32 & surfclu, n_ele, jellium, R_j, h_j, e_j, & 33 & nelect, P_ext 34 USE mp_world, ONLY: nproc, mpime 35 USE io_global, ONLY: ionode 36 USE mp, ONLY: mp_bcast, mp_sum 37 USE mp_bands, ONLY: intra_bgrp_comm 38 USE fft_rho 39 40 implicit none 41 42 real(kind=8) dx, dxx, xcc(4800) 43 real(kind=8) weight0, wpiu, wmeno, maxr, minr 44 real(kind=8) tau00(3), dist 45 real(kind=8) rho_real(dfftp%nnr,nspin), rhoc 46 real(kind=8) alfa0, sigma, hgt 47 real(kind=8) pos_cry(3), pos_car(3), pos_aux(3) 48 real(kind=8) pos_cry0(3), dpvdh(3,3) 49 real(kind=8) v_d(3) 50 real(kind=8) mtot, rad0, cm(3) 51 real(kind=8) modr, lap 52 real(kind=8) prod, aux1 53 real(kind=8) gxl, xyr, xzr, yzr 54 real(kind=8), allocatable:: vec(:,:,:), aiuto(:,:,:) 55 real(kind=8), allocatable:: drho(:,:), d2rho(:,:) 56 real(kind=8), allocatable:: tauv(:,:) 57 58 complex(kind=8) ci 59 complex(kind=8) sum_sf, aux, auxx, fact, rho_g(dfftp%ngm,nspin) 60 complex(kind=8), allocatable :: rhofill(:), rhotmp(:,:) 61 62 integer ir, ir1, ir2, ir3, is, iss, ia, flag, ierr 63 integer i, j, k, l, ig, cnt, nmin, nmax 64 65#if defined(__MPI) 66 real(kind=8) maxr_p(nproc), minr_p(nproc), maxr_pp, minr_pp 67 integer shift(nproc), incr(nproc), ppp(nproc) 68 integer displs(nproc), ip, me 69#endif 70 if (abisur) allocate(drho(3,dfftp%nnr)) 71 if (abisur) allocate(d2rho(6,dfftp%nnr)) 72 73 call start_clock( 'vol_clu' ) 74 75 ci = (0.d0,1.d0) 76 77#if defined(__MPI) 78 me = mpime + 1 79 do ip=1,nproc 80 ppp(ip) = dfftp%nnp * ( dfftp%nr3p(ip) ) 81 if (ip.eq.1) then 82 shift(ip)=0 83 else 84 shift(ip)=shift(ip-1) + ppp(ip-1) 85 end if 86 end do 87#endif 88 89 sigma = rho_thr/3.d0 !3.d0 90 hgt = 0.0050d0 !5000.d0*rho_thr 91! We smear the step function defining the volume and approximate its derivative 92! with a gaussian. Here we sample the integral of this gaussian. It has to 93! be done once for ever 94! XXX: using an array for xcc() is a big waste. two scalar variables would do. 95 dx = 5.d0*sigma/60.d0 96 if (flag.eq.1) then 97 dxx = dx/40.d0 98 weight(1) = 0.d0 99 xcc(1) = rho_thr - 5.d0*sigma 100 xc0(1) = xcc(1) 101 cnt = 1 102 do i = 2,121 103 weight(i) = weight(i-1) 104 do j = 1,40 105 cnt = cnt + 1 106 xcc(cnt) = xcc(cnt-1) + dxx 107 if (j.eq.40) then 108 xc0(i) = xcc(cnt) 109 end if 110 aux1 = xcc(cnt)-dxx/2.d0-rho_thr 111 weight(i) = weight(i) + 1.d0/(sigma*dsqrt(pi*2.d0)) * & 112 & dxx * dexp(-1.d0*aux1**2/(2.d0*sigma**2)) 113 end do 114 end do 115! This doesn't work yet..... 116 if (jellium) then 117 do ir3 = 1,dfftp%nr3 118 do ir2 = 1,dfftp%nr2 119 do ir1 = 1,dfftp%nr1 120 ir = ir1 + (ir2-1)*dfftp%nr1 + (ir3-1)*dfftp%nr2*dfftp%nr1 121 dist = 0.d0 122 do i = 1,3 123 posv(i,ir) = (DBLE(ir1)-1.0d0)*at(i,1)/DBLE(dfftp%nr1) +& 124 & (DBLE(ir2)-1.0d0)*at(i,2)/DBLE(dfftp%nr2) +& 125 & (DBLE(ir3)-1.0d0)*at(i,3)/DBLE(dfftp%nr3) 126 end do 127 end do 128 end do 129 end do 130 posv(:,:) = posv(:,:)*alat 131 end if 132 end if 133 134 allocate ( tauv(3,nat) ) 135 do ia = 1,nat 136 is=ityp(ia) 137 ! alfa(is) = step_rad(is)/2.d0 (not used) 138 do k = 1,3 139 tauv(k,ia) = tau0(k,ia) 140 end do 141 end do 142 143 stress_vol = 0.d0 144 dpvdh = 0.d0 145 146! Now we compute the volume and other quantities 147 148 volclu = 0.d0 149 n_ele = 0.d0 150 surfclu = 0.d0 151 152! Let's add rhops to fill possible holes in the valence charge density on top 153! of the ions 154 155 allocate(rhotmp(dfftp%ngm,nspin)) 156 rhotmp = (0.d0,0.d0) 157 158 if (nspin.eq.1) then 159 do ig = 1,dfftp%ngm 160 rhotmp(ig,1)=rho_g(ig,1) 161 end do 162 else 163 do ig = 1,dfftp%ngm 164 do iss = 1,2 165 rhotmp(ig,iss) = rho_g(ig,iss) 166 end do 167 end do 168 end if 169 170! To fill the vacuum inside hollow structures 171 172 if (fill_vac) then 173 allocate(rhofill(dfftp%ngm)) 174 rhofill = 0.d0 175 do k = 1,3 176 cm(k) = 0.d0 177 mtot = 0.d0 178 do ia = 1,nat 179 cm(k) = cm(k) + tauv(k,ia)*amass(ityp(ia)) 180 mtot = mtot + amass(ityp(ia)) 181 end do 182 cm(k) = cm(k)/mtot 183 end do 184 end if 185 186 if (fill_vac) then 187 do i = 1,n_cntr 188 do ia = 1,nat 189 is = ityp(ia) 190 if (cntr(is)) then 191 rad0 = step_rad(is) + DBLE(i)*delta_sigma 192 alfa0 = rad0/2.d0 193 do k = 1,3 194 if (k.ne.axis) then 195 tau00(k) = (tauv(k,ia)-cm(k))* & 196 & (1.d0-delta_eps*DBLE(i))+cm(k) 197 else 198 tau00(k) = tauv(k,ia) 199 end if 200 end do 201 do ig = 1,dfftp%ngm 202 prod = 0.d0 203 do k = 1,3 204 prod = prod + g(k,ig)*tau00(k) 205 end do 206 prod = prod*tpiba 207 fact = CMPLX(cos(prod),-1.d0*sin(prod),kind=DP) 208 aux = alfa0*hgt*EXP(-(0.50d0*alfa0**2*gg(ig)*tpiba2)) 209 rhofill(ig) = rhofill(ig) + aux*fact 210 end do 211 end if 212 end do 213 end do 214 if (nspin.eq.1) then 215 do ig=1,dfftp%ngm 216 rhotmp(ig,1) = rhotmp(ig,1) + rhofill(ig) 217 end do 218 else 219 do ig = 1,dfftp%ngm 220 do iss = 1,2 221 rhotmp(ig,iss) = rhotmp(ig,iss) + 0.5d0*rhofill(ig) 222 end do 223 end do 224 end if 225 end if 226 227 if (fill_vac) then 228 deallocate(rhofill) 229 end if 230 231 IF (abisur) THEN 232 DO iss = 1, nspin 233 CALL fft_gradient_g2r( dfftp, rhotmp(1,iss), g, drho(1,iss) ) 234 CALL fft_hessian_g2r ( dfftp, rhotmp, g, d2rho(1,iss) ) 235 END DO 236 END IF 237 238 CALL rho_g2r( dfftp, rhotmp, rho_gaus ) 239 deallocate(rhotmp) 240 241 e_j = 0.d0 242 243 do ir = 1,dfftp%nnr 244 245 v_vol(ir) = 0.d0 246 247 if (jellium) then 248#if defined(__MPI) 249 do j = 1,3 250 pos_aux(j) = posv(j,ir+shift(me)) 251 end do 252#else 253 do j = 1,3 254 pos_aux(j) = posv(j,ir) 255 end do 256#endif 257 dist = 0.d0 258 do j = 1,3 259 dist = dist + (pos_aux(j) - 0.5d0*(at(j,1)+at(j,2)+at(j,3)))**2 260 end do 261 dist = dsqrt(dist)*alat 262 if (dist.ge.R_j) then 263 v_vol(ir) = - nelect/dist 264 v_vol(ir) = 0.d0 265 else 266! The last term in the internal potential is for its continuity 267 v_vol(ir) = + 0.5d0*nelect*dist**2/R_j**3 & 268 - 1.5d0*nelect/R_j 269 v_vol(ir) = - h_j 270 end if 271 if (nspin.eq.1) then 272 e_j = e_j + v_vol(ir) * rho_real(ir,1) * omega / & 273 & DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 274 else 275 e_j = e_j + v_vol(ir) * & 276 ( rho_real(ir,1) + rho_real(ir,2) ) * omega / & 277 & DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 278 end if 279 end if 280 281 rhoc = rho_gaus(ir) 282! Volume and surface 283 if (rhoc.gt.rho_thr+5.d0*sigma) then 284 weight0 = 1.d0 285 wpiu = 1.d0 286 i = int((rhoc-rho_thr-dthr+5.d0*sigma)/dx) + 1 287 if (i.gt.120) then 288 wmeno = 1.d0 289 else 290 wmeno = weight(i) + (weight(i+1)-weight(i)) * & 291 & (rhoc-rho_thr-dthr-DBLE(i-1)*dx+5.d0*sigma)/dx 292 end if 293 go to 79 294 end if 295! Volume and surface 296 k = int((rhoc-rho_thr+5.d0*sigma)/dx) + 1 297 weight0 = weight(k) + (weight(k+1)-weight(k)) * & 298 (rhoc-rho_thr+5.d0*sigma-DBLE(k-1)*dx)/dx 299 if (abisur) then 300 if (rhoc-rho_thr+dthr.gt.5.d0*sigma) then 301 wpiu = weight0 302 i = int((rhoc-rho_thr-dthr+5.d0*sigma)/dx) + 1 303 wmeno = weight(i)+(weight(i+1)-weight(i))* & 304 & (rhoc-rho_thr-dthr+5.d0*sigma-DBLE(i-1)*dx)/dx 305 else if (rho_thr+dthr-rhoc.gt.5.d0*sigma) then 306 wmeno = 0.d0 307 i = int((rhoc-rho_thr+dthr+5.d0*sigma)/dx) + 1 308 wpiu = weight0 309 else 310 i = int((rhoc-rho_thr+dthr+5.d0*sigma)/dx) + 1 311 wpiu = weight0 312 i = int((rhoc-rho_thr-dthr+5.d0*sigma)/dx) + 1 313 wmeno = weight(i)+(weight(i+1)-weight(i))* & 314 & (rhoc-rho_thr-dthr+5.d0*sigma-DBLE(i-1)*dx)/dx 315 end if 316 end if 317 79 continue 318 if (nspin.eq.1) then 319 n_ele = n_ele + weight0 * rho_real(ir,1) 320 else 321 n_ele = n_ele + weight0 * (rho_real(ir,1) + rho_real(ir,2)) 322 end if 323 volclu = volclu + weight0 324 v_vol(ir) = v_vol(ir) + P_ext /(sigma*dsqrt(pi*2.d0)) * & 325 & dexp(-1.d0*(rhoc-rho_thr)**2/(2.d0*sigma**2)) 326 if (tpre) then 327 do k = 1,3 328 do j = 1,3 329 do is = 1,nspin 330 dpvdh(k,j) = dpvdh(k,j) + & 331 & v_vol(ir)*drhor(ir,is,k,j)*omega/ & 332 & DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 333 end do 334 end do 335 end do 336 end if 337 338 if (abisur) then 339 ! for d2rho: 1=xx, 2=xy, 3=yy, 4=xz, 5=yz, 6=zz 340 modr = drho(1,ir)**2 + drho(2,ir)**2 + drho(3,ir)**2 341 lap = d2rho(1,ir) + d2rho(3,ir) + d2rho(6,ir) 342 gxl = drho(1,ir)**2*d2rho(1,ir) + & 343 drho(2,ir)**2*d2rho(3,ir) + & 344 drho(3,ir)**2*d2rho(6,ir) 345 xyr = 2.d0*d2rho(2,ir)*drho(1,ir)*drho(2,ir) 346 xzr = 2.d0*d2rho(4,ir)*drho(1,ir)*drho(3,ir) 347 yzr = 2.d0*d2rho(5,ir)*drho(2,ir)*drho(3,ir) 348 modr = dsqrt(modr) 349 surfclu = surfclu + (wpiu-wmeno)*modr 350 v_vol(ir) = v_vol(ir) -1.d0*Surf_t/dthr * (wpiu-wmeno) * & 351 & (lap/modr - (gxl + xyr + xzr + yzr)/modr**3) 352 end if 353 354 end do 355 356 call mp_sum(volclu,intra_bgrp_comm) 357 call mp_sum(n_ele,intra_bgrp_comm) 358 if (jellium) call mp_sum(e_j,intra_bgrp_comm) 359 call mp_sum(surfclu,intra_bgrp_comm) 360 call mp_sum(dpvdh,intra_bgrp_comm) 361 362 volclu = volclu * omega / DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 363 n_ele = n_ele * omega / DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 364 surfclu = surfclu * omega / DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) / dthr 365 do i = 1,3 366 do j = 1,3 367 stress_vol(i,j) = dpvdh(i,1)*h(j,1) + dpvdh(i,2)*h(j,2) + & 368 & dpvdh(i,3)*h(j,3) 369 end do 370 end do 371 372 deallocate( tauv ) 373 if ( abisur ) deallocate( drho ) 374 if ( abisur ) deallocate( d2rho ) 375 376 call stop_clock( 'vol_clu' ) 377 378END SUBROUTINE vol_clu 379