1! 2! Copyright (C) 2002-2007 Quantum ESPRESSO 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 9 10!=----------------------------------------------------------------------------=! 11 subroutine core_charge_ftr( tpre ) 12!=----------------------------------------------------------------------------=! 13 ! 14 ! Compute the fourier transform of the core charge, from the radial 15 ! mesh to the reciprocal space 16 ! 17 use kinds, ONLY : DP 18 use ions_base, ONLY : nsp 19 use atom, ONLY : rgrid 20 use uspp, ONLY : nlcc_any 21 use uspp_param, ONLY : upf 22 use smallbox_gvec, ONLY : ngb, gb 23 use small_box, ONLY : omegab, tpibab 24 use pseudo_base, ONLY : compute_rhocg 25 use pseudopotential, ONLY : tpstab, rhoc1_sp, rhocp_sp 26 use cell_base, ONLY : omega, tpiba2, tpiba 27 USE splines, ONLY : spline 28 use gvect, ONLY : gg, gstart 29 USE core, ONLY : rhocb, rhocg, drhocg 30 USE fft_base, ONLY: dfftp 31 ! 32 IMPLICIT NONE 33 ! 34 LOGICAL, INTENT(IN) :: tpre 35 ! 36 INTEGER :: is, ig 37 REAL(DP) :: xg, cost1 38 ! 39 ! 40 IF( .NOT. nlcc_any ) RETURN 41 ! 42 IF( .NOT. ALLOCATED( rgrid ) ) & 43 CALL errore( ' core_charge_ftr ', ' rgrid not allocated ', 1 ) 44 IF( .NOT. ALLOCATED( upf ) ) & 45 CALL errore( ' core_charge_ftr ', ' upf not allocated ', 1 ) 46 ! 47 do is = 1, nsp 48 ! 49 if( upf(is)%nlcc ) then 50 ! 51 CALL compute_rhocg( rhocb(:,is), rhocb(:,is), rgrid(is)%r, & 52 rgrid(is)%rab, upf(is)%rho_atc(:), gb, omegab, tpibab**2, & 53 rgrid(is)%mesh, ngb, 0 ) 54 ! 55 IF( tpre ) THEN 56 ! 57 IF( tpstab ) THEN 58 ! 59 cost1 = 1.0d0/omega 60 ! 61 IF( gstart == 2 ) THEN 62 rhocg (1,is) = rhoc1_sp(is)%y( 1 ) * cost1 63 drhocg(1,is) = 0.0d0 64 END IF 65 DO ig = gstart, SIZE( rhocg, 1 ) 66 xg = SQRT( gg(ig) ) * tpiba 67 rhocg (ig,is) = spline( rhoc1_sp(is), xg ) * cost1 68 drhocg(ig,is) = spline( rhocp_sp(is), xg ) * cost1 69 END DO 70 ! 71 ELSE 72 73 CALL compute_rhocg( rhocg(:,is), drhocg(:,is), rgrid(is)%r, & 74 rgrid(is)%rab, upf(is)%rho_atc(:), gg, & 75 omega, tpiba2, rgrid(is)%mesh, dfftp%ngm, 1 ) 76 77 END IF 78 ! 79 END IF 80 ! 81 endif 82 ! 83 end do 84 85 return 86 end subroutine core_charge_ftr 87 88 89 90!----------------------------------------------------------------------- 91 subroutine add_cc( rhoc, rhog, rhor ) 92!----------------------------------------------------------------------- 93! 94! add core correction to the charge density for exch-corr calculation 95! 96 USE kinds, ONLY: DP 97 use electrons_base, only: nspin 98 use control_flags, only: iverbosity 99 use io_global, only: stdout 100 use mp_global, only: intra_bgrp_comm 101 use cell_base, only: omega 102 USE mp, ONLY: mp_sum 103 104 ! this isn't really needed, but if I remove it, ifc 7.1 105 ! gives an "internal compiler error" 106 use gvect, only: gstart 107 USE fft_interfaces, ONLY: fwfft 108 USE fft_base, ONLY: dfftp 109 USE fft_helper_subroutines, ONLY: fftx_add_threed2oned_gamma 110! 111 implicit none 112 ! 113 REAL(DP), INTENT(IN) :: rhoc( dfftp%nnr ) 114 REAL(DP), INTENT(INOUT):: rhor( dfftp%nnr, nspin ) 115 COMPLEX(DP), INTENT(INOUT):: rhog( dfftp%ngm, nspin ) 116 ! 117 COMPLEX(DP), ALLOCATABLE :: wrk1( : ) 118! 119 integer :: ig, ir, iss, isup, isdw 120 REAL(DP) :: rsum 121 ! 122 IF( iverbosity > 1 ) THEN 123 rsum = SUM( rhoc ) * omega / DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 124 CALL mp_sum( rsum, intra_bgrp_comm ) 125 WRITE( stdout, 10 ) rsum 12610 FORMAT( 3X, 'Core Charge = ', D14.6 ) 127 END IF 128 ! 129 ! In r-space: 130 ! 131 if ( nspin .eq. 1 ) then 132 iss=1 133 call daxpy(dfftp%nnr,1.d0,rhoc,1,rhor(1,iss),1) 134 else 135 isup=1 136 isdw=2 137 call daxpy(dfftp%nnr,0.5d0,rhoc,1,rhor(1,isup),1) 138 call daxpy(dfftp%nnr,0.5d0,rhoc,1,rhor(1,isdw),1) 139 end if 140 ! 141 ! rhoc(r) -> rhoc(g) (wrk1 is used as work space) 142 ! 143 allocate( wrk1( dfftp%nnr ) ) 144 145 wrk1(:) = rhoc(:) 146 147 call fwfft('Rho',wrk1, dfftp ) 148 ! 149 ! In g-space: 150 ! 151 if (nspin.eq.1) then 152 CALL fftx_add_threed2oned_gamma( dfftp, wrk1, rhog(:,iss) ) 153 else 154 wrk1 = wrk1 * 0.5d0 155 CALL fftx_add_threed2oned_gamma( dfftp, wrk1, rhog(:,isup) ) 156 CALL fftx_add_threed2oned_gamma( dfftp, wrk1, rhog(:,isdw) ) 157 end if 158 159 deallocate( wrk1 ) 160! 161 return 162 end subroutine add_cc 163 164 165! 166!----------------------------------------------------------------------- 167 subroutine force_cc(irb,eigrb,vxc,fion1) 168!----------------------------------------------------------------------- 169! 170! core correction force: f = \int V_xc(r) (d rhoc(r)/d R_i) dr 171! same logic as in newd - uses box grid. For parallel execution: 172! the sum over node contributions is done in the calling routine 173! 174 USE kinds, ONLY: DP 175 use electrons_base, only: nspin 176 use smallbox_gvec, only: gxb, ngb 177 use smallbox_subs, only: fft_oned2box, boxdotgrid 178 use cell_base, only: omega 179 use ions_base, only: nsp, na, nat, ityp 180 use small_box, only: tpibab 181 use uspp_param, only: upf 182 use core, only: rhocb 183 use fft_interfaces, only: invfft 184 use fft_base, only: dfftb, dfftp 185 use gvect, only: gstart 186 187 implicit none 188 189! input 190 integer, intent(in) :: irb(3,nat) 191 complex(dp), intent(in):: eigrb(ngb,nat) 192 real(dp), intent(in) :: vxc(dfftp%nnr,nspin) 193! output 194 real(dp), intent(inout):: fion1(3,nat) 195! local 196 integer :: iss, ix, ig, is, ia, nfft, isa 197 real(dp) :: fac, res 198 complex(dp) ci, facg 199 complex(dp), allocatable :: qv(:), fg1(:), fg2(:) 200 real(dp), allocatable :: fcc(:,:) 201 202#if defined(_OPENMP) 203 INTEGER :: itid, mytid, ntids, omp_get_thread_num, omp_get_num_threads 204 EXTERNAL :: omp_get_thread_num, omp_get_num_threads 205#endif 206! 207 call start_clock( 'forcecc' ) 208 ci = (0.d0,1.d0) 209 210 fac = omega/DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3*nspin) 211 212!$omp parallel default(none) & 213!$omp shared(nsp, na, ngb, eigrb, dfftb, irb, ci, rhocb, & 214!$omp gxb, nat, fac, upf, vxc, nspin, tpibab, fion1, ityp ) & 215!$omp private(mytid, ntids, is, ia, nfft, ig, qv, fg1, fg2, itid, res, ix, fcc, facg, iss ) 216 217 218 allocate( fcc( 3, nat ) ) 219 allocate( qv( dfftb%nnr ) ) 220 allocate( fg1( ngb ) ) 221 allocate( fg2( ngb ) ) 222 223 fcc(:,:) = 0.d0 224 225#if defined(_OPENMP) 226 mytid = omp_get_thread_num() ! take the thread ID 227 ntids = omp_get_num_threads() ! take the number of threads 228 itid = 0 229#endif 230 231 do ia = 1, nat 232 233 is = ityp(ia) 234 235 if( .not. upf(is)%nlcc ) then 236 cycle 237 end if 238 239 nfft = 1 240 241#if defined(__MPI) 242 if ( ( dfftb%np3( ia ) <= 0 ) .OR. ( dfftb%np2( ia ) <= 0 ) ) & 243 CYCLE 244#endif 245 246#if defined(_OPENMP) 247 IF ( mytid /= itid ) THEN 248 itid = MOD( itid + 1, ntids ) 249 CYCLE 250 ELSE 251 itid = MOD( itid + 1, ntids ) 252 END IF 253#endif 254 255 do ix=1,3 256 if (nfft.eq.2) then 257 do ig=1,ngb 258 facg = tpibab*CMPLX(0.d0,gxb(ix,ig),kind=DP)*rhocb(ig,is) 259 fg1(ig) = eigrb(ig,ia )*facg 260 fg2(ig) = eigrb(ig,ia+1)*facg 261 end do 262 CALL fft_oned2box( qv, fg1, fg2 ) 263 else 264 do ig=1,ngb 265 facg = tpibab*CMPLX(0.d0,gxb(ix,ig),kind=DP)*rhocb(ig,is) 266 fg1(ig) = eigrb(ig,ia)*facg 267 end do 268 CALL fft_oned2box( qv, fg1 ) 269 end if 270! 271 call invfft( qv, dfftb, ia ) 272 ! 273 ! note that a factor 1/2 is hidden in fac if nspin=2 274 ! 275 do iss=1,nspin 276 res = boxdotgrid(irb(:,ia),1,qv,vxc(:,iss)) 277 fcc(ix,ia) = fcc(ix,ia) + fac * res 278 if (nfft.eq.2) then 279 res = boxdotgrid(irb(:,ia+1),2,qv,vxc(:,iss)) 280 fcc(ix,ia+1) = fcc(ix,ia+1) + fac * res 281 end if 282 end do 283 end do 284 end do 285 286! 287!$omp critical 288 do ia = 1, nat 289 fion1(:,ia) = fion1(:,ia) + fcc(:,ia) 290 end do 291!$omp end critical 292 293 deallocate( qv ) 294 deallocate( fg1 ) 295 deallocate( fg2 ) 296 deallocate( fcc ) 297 298!$omp end parallel 299! 300 call stop_clock( 'forcecc' ) 301 302 return 303 end subroutine force_cc 304 305 306! 307!----------------------------------------------------------------------- 308 subroutine set_cc( irb, eigrb, rhoc ) 309!----------------------------------------------------------------------- 310! 311! Calculate core charge contribution in real space, rhoc(r) 312! Same logic as for rhov: use box grid for core charges 313! 314 use kinds, only: dp 315 use ions_base, only: nsp, na, nat, ityp 316 use uspp_param, only: upf 317 use smallbox_gvec, only: ngb 318 use smallbox_subs, only: fft_oned2box, box2grid 319 use control_flags, only: iprint 320 use core, only: rhocb 321 use fft_interfaces, only: invfft 322 use fft_base, only: dfftb, dfftp 323 324 implicit none 325! input 326 integer, intent(in) :: irb(3,nat) 327 complex(dp), intent(in):: eigrb(ngb,nat) 328! output 329 real(dp), intent(out) :: rhoc(dfftp%nnr) 330! local 331 integer nfft, ig, is, ia, isa 332 complex(dp) ci 333 complex(dp), allocatable :: wrk1(:) 334 complex(dp), allocatable :: qv(:), fg1(:), fg2(:) 335 336#if defined(_OPENMP) 337 INTEGER :: itid, mytid, ntids, omp_get_thread_num, omp_get_num_threads 338 EXTERNAL :: omp_get_thread_num, omp_get_num_threads 339#endif 340! 341 call start_clock( 'set_cc' ) 342 ci=(0.d0,1.d0) 343 344 allocate( wrk1 ( dfftp%nnr ) ) 345 wrk1 (:) = (0.d0, 0.d0) 346! 347!$omp parallel default(none) & 348!$omp shared(nsp, na, ngb, eigrb, dfftb, irb, ci, rhocb, & 349!$omp nat, upf, wrk1, ityp ) & 350!$omp private(mytid, ntids, is, ia, nfft, ig, isa, qv, fg1, fg2, itid ) 351 352 allocate( qv ( dfftb%nnr ) ) 353 allocate( fg1 ( ngb ) ) 354 allocate( fg2 ( ngb ) ) 355! 356 isa = 0 357 358#if defined(_OPENMP) 359 mytid = omp_get_thread_num() ! take the thread ID 360 ntids = omp_get_num_threads() ! take the number of threads 361 itid = 0 362#endif 363 364 do ia = 1, nat 365 ! 366 is = ityp(ia) 367 368 if (.not.upf(is)%nlcc) then 369 cycle 370 end if 371 ! 372#if defined(__MPI) 373 nfft=1 374 if ( ( dfftb%np3( ia ) <= 0 ) .OR. ( dfftb%np2 ( ia ) <= 0 ) ) cycle 375#endif 376 377#if defined(_OPENMP) 378 IF ( mytid /= itid ) THEN 379 itid = MOD( itid + 1, ntids ) 380 CYCLE 381 ELSE 382 itid = MOD( itid + 1, ntids ) 383 END IF 384#endif 385 386 if(nfft.eq.2)then 387 fg1 = eigrb(1:ngb,ia )*rhocb(1:ngb,is) 388 fg2 = eigrb(1:ngb,ia+1)*rhocb(1:ngb,is) 389 CALL fft_oned2box( qv, fg1, fg2 ) 390 else 391 fg1 = eigrb(1:ngb,ia )*rhocb(1:ngb,is) 392 CALL fft_oned2box( qv, fg1 ) 393 endif 394! 395 call invfft( qv, dfftb, ia ) 396! 397 call box2grid(irb(:,ia),1,qv,wrk1) 398 if (nfft.eq.2) call box2grid(irb(:,ia+1),2,qv,wrk1) 399! 400 end do 401! 402 deallocate( qv ) 403 deallocate( fg1 ) 404 deallocate( fg2 ) 405 406!$omp end parallel 407 408 call dcopy( dfftp%nnr, wrk1, 2, rhoc, 1 ) 409 410 deallocate( wrk1 ) 411! 412 call stop_clock( 'set_cc' ) 413! 414 return 415 end subroutine set_cc 416 417