1! 2! Copyright (C) 2017 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!---------------------------------------------------------------------------- 10SUBROUTINE stres_mgga( sigmaxc ) 11 !---------------------------------------------------------------------------- 12 ! 13 ! Analytic stress tensor contribution from metagga is added to sigmaxc 14 ! 15 USE kinds, ONLY : DP 16 USE control_flags, ONLY : gamma_only 17 USE noncollin_module, ONLY : noncolin 18 USE cell_base, ONLY : alat, at, bg, omega, tpiba 19 USE gvect, ONLY : g 20 USE scf, ONLY : rho, v 21 USE wavefunctions, ONLY : evc 22 USE funct, ONLY : dft_is_meta 23 USE klist, ONLY : nks, xk, ngk 24 USE buffers, ONLY : get_buffer 25 USE io_files, ONLY : iunwfc, nwordwfc 26 USE wvfct, ONLY : nbnd, npwx, wg 27 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 28 USE fft_interfaces, ONLY : fwfft, invfft 29 USE fft_base, ONLY : dfftp, dffts 30 USE mp, ONLY : mp_sum 31 USE mp_pools, ONLY : inter_pool_comm 32 USE mp_bands, ONLY : intra_bgrp_comm 33 ! 34 IMPLICIT NONE 35 ! 36 REAL(DP), INTENT(INOUT) :: sigmaxc(3,3) 37 ! 38 ! Internal variables 39 ! 40 INTEGER :: ix, iy, ir, iss, ipol, incr, ibnd, ik, npw 41 INTEGER :: ipol2xy(3,3) 42 !! ipol2xy(i,j) = ipol2x(j,i) is a collapsed symmetric index 43 DATA ipol2xy / 1, 2, 3, 2, 4, 5, 3, 5, 6/ 44 REAL(DP), PARAMETER :: epsr = 1.0d-6, epsg = 1.0d-10, e2 = 2.d0 45 COMPLEX(DP), ALLOCATABLE :: gradwfc (:,:), crosstaus(:,:,:) 46 REAL(DP) :: w1, w2, delta, sigma_mgga(3,3) 47 ! 48 if ( .not. dft_is_meta() ) return 49 ! 50 current_spin=1 51 ! 52 ! 53 ! Stop if something is not yet implemented 54 ! 55 if (noncolin) call errore('stres_mgga', & 56 'noncollinear stress + meta-GGA not implemented',1) 57 ! 58 ! Initialization of a set of variables 59 ! 60 allocate (gradwfc( dffts%nnr, 3)) 61 allocate (crosstaus( dffts%nnr,6,nspin)) 62 ! 63 ! For gamma_only efficiency 64 ! 65 incr=1 66 IF ( gamma_only ) incr=2 67 ! 68 crosstaus(:,:,:) = 0.d0 69 gradwfc(:,:) = 0.d0 70 ! 71 ! Loop over the k points 72 ! 73 k_loop: DO ik = 1, nks 74 ! 75 ! 76 IF ( lsda ) current_spin = isk(ik) 77 ! 78 npw = ngk(ik) 79 ! 80 ! Read the wavefunctions 81 ! 82 IF ( nks > 1 ) THEN 83 ! 84 CALL get_buffer ( evc, nwordwfc, iunwfc, ik ) 85 ! 86 END IF 87 ! 88 do ibnd = 1, nbnd, incr 89 ! 90 ! w1, w2: weights for each k point and band 91 ! 92 w1 = wg(ibnd,ik) / omega 93 ! 94 IF ( (ibnd < nbnd) .and. (gamma_only) ) THEN 95 ! 96 ! ... two ffts at the same time 97 ! 98 w2 = wg(ibnd+1,ik) / omega 99 ! 100 ELSE 101 ! 102 w2 = w1 103 ! 104 END IF 105 ! 106 ! Gradient of the wavefunction in real space 107 ! 108 CALL wfc_gradient( ibnd, ik, npw, gradwfc ) 109 ! 110 ! Cross terms of kinetic energy density 111 ! 112 do ix=1,3 113 ! 114 do iy=1,ix 115 ! 116 ipol = ipol2xy(iy,ix) 117 ! 118 do ir=1,dffts%nnr 119 ! 120 crosstaus(ir,ipol,current_spin) = crosstaus(ir,ipol,current_spin) +& 121 2.0_DP*w1*DBLE(gradwfc(ir,ix))*DBLE(gradwfc(ir,iy)) +& 122 2.0_DP*w2*AIMAG(gradwfc(ir,ix))*AIMAG(gradwfc(ir,iy)) 123 ! 124 end do 125 ! 126 end do 127 ! 128 end do 129 ! 130 end do 131 ! 132 END DO k_loop 133 ! 134 call mp_sum( crosstaus, inter_pool_comm ) 135 ! 136 ! gradwfc not used anymore 137 ! 138 deallocate (gradwfc) 139 ! 140 sigma_mgga(:,:) = 0.D0 141 ! 142 ! metagga contribution to the stress tensor 143 ! 144 do iss=1,nspin 145 ! 146 do ix=1,3 147 ! 148 do iy=1,3 149 ! 150 delta=0. 151 if (ix==iy) delta=1. 152 ! 153 do ir=1,dffts%nnr 154 ! 155 sigma_mgga(ix,iy) = sigma_mgga(ix,iy) + v%kin_r(ir,iss) & 156 * ( rho%kin_r(ir,iss) * delta & 157 + crosstaus(ir,ipol2xy(ix,iy),iss) ) 158 ! 159 end do 160 ! 161 ! 162 end do 163 ! 164 end do 165 ! 166 end do 167 deallocate( crosstaus ) 168 ! 169 call mp_sum( sigma_mgga, intra_bgrp_comm ) 170 sigmaxc(:,:) = sigmaxc(:,:) + sigma_mgga(:,:) / & 171 (dffts%nr1 * dffts%nr2 * dffts%nr3) 172 ! 173 return 174 ! 175END SUBROUTINE stres_mgga 176 177SUBROUTINE wfc_gradient ( ibnd, ik, npw, gradpsi ) 178 ! 179 ! Returns the gradient of the wavefunction in real space 180 ! Slightly adapted from sum_bands.f90 181 ! 182 USE kinds, ONLY : DP 183 USE control_flags, ONLY : gamma_only 184 USE wavefunctions, ONLY : psic, evc 185 USE wvfct, ONLY : npwx, nbnd 186 USE cell_base, ONLY : omega, tpiba 187 USE klist, ONLY : xk, igk_k 188 USE gvect, ONLY : g 189 USE fft_base, ONLY : dffts 190 USE fft_interfaces, ONLY : invfft 191 ! 192 IMPLICIT NONE 193 ! 194 INTEGER :: ibnd, ik, npw 195 COMPLEX(DP) :: gradpsi(dffts%nnr,3) 196 ! 197 ! Internal variables 198 ! 199 REAL(DP) :: kplusg(npwx) 200 INTEGER :: ipol 201 ! 202 ! Compute the gradient of the wavefunction in reciprocal space 203 ! 204 IF ( gamma_only ) THEN 205 ! 206 DO ipol=1,3 207 ! 208 psic(:) = ( 0.D0, 0.D0 ) 209 ! 210 kplusg (1:npw) = (xk(ipol,ik)+g(ipol,igk_k(1:npw,ik))) * tpiba 211 ! 212 IF ( ibnd < nbnd ) THEN 213 ! 214 ! ... two ffts at the same time 215 ! 216 psic(dffts%nl(1:npw)) = CMPLX(0d0, kplusg(1:npw),kind=DP)* & 217 ( evc(1:npw,ibnd) + & 218 ( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1) ) 219 ! 220 psic(dffts%nlm(1:npw)) = CMPLX(0d0,-kplusg(1:npw),kind=DP) * & 221 CONJG( evc(1:npw,ibnd) - & 222 ( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1) ) 223 ! 224 ELSE 225 ! 226 psic(dffts%nl(1:npw)) = CMPLX(0d0, kplusg(1:npw),kind=DP)* & 227 evc(1:npw,ibnd) 228 ! 229 psic(dffts%nlm(1:npw)) = CMPLX(0d0,-kplusg(1:npw),kind=DP) * & 230 CONJG( evc(1:npw,ibnd) ) 231 ! 232 END IF 233 ! 234 ! Gradient of the wavefunction in real space 235 ! 236 CALL invfft ('Wave', psic, dffts) 237 ! 238 gradpsi(:,ipol) = psic 239 ! 240 END DO 241 ! 242 ELSE 243 ! 244 DO ipol=1,3 245 ! 246 psic(:) = ( 0.D0, 0.D0 ) 247 ! 248 kplusg (1:npw) = (xk(ipol,ik)+g(ipol,igk_k(1:npw,ik))) * tpiba 249 psic(dffts%nl(igk_k(1:npw,ik))) = CMPLX(0d0, kplusg(1:npw),kind=DP)* & 250 evc(1:npw,ibnd) 251 ! 252 ! Gradient of the wavefunction in real space 253 ! 254 CALL invfft ('Wave', psic, dffts) 255 ! 256 gradpsi(:,ipol) = psic 257 ! 258 END DO 259 ! 260 END IF 261 ! 262END SUBROUTINE wfc_gradient 263