1! 2! Copyright (C) 2005-2010 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 subroutine dforce_meta (c,ca,df,da, psi,iss1,iss2,fi,fip) 9!----------------------------------------------------------------------- 10!computes: the generalized force df=cmplx(dfr,dfi) acting on the i-th 11! electron state at the gamma point of the brillouin zone 12! represented by the vector c=cmplx(cr,ci) 13! 14! contribution from metaGGA 15 use kinds, only: dp 16 use gvect, only : g 17 use gvecw, only : ngw 18 use cell_base, only : tpiba2 19 USE metagga_cp, ONLY : kedtaus 20 USE fft_interfaces, ONLY : fwfft, invfft 21 USE fft_base, ONLY : dffts 22 USE fft_helper_subroutines, ONLY : c2psi_gamma, fftx_psi2c_gamma 23! 24 implicit none 25! 26 complex(dp) c(ngw), ca(ngw), df(ngw), da(ngw),psi(dffts%nnr) 27 complex(dp), allocatable :: dc(:), dca(:) 28 integer iss1, iss2 29 real(dp) fi, fip 30! local variables 31 integer ir,ig, ipol !metagga 32 complex(dp) fp,fm,ci 33! 34! 35 ci=(0.0d0,1.0d0) 36 allocate( dc( ngw ) ) 37 allocate( dca( ngw ) ) 38! 39 do ipol = 1, 3 40 41 dc(:) = ci*g(ipol,1:ngw)*c(:) 42 dca(:) = ci*g(ipol,1:ngw)*ca(:) 43 CALL c2psi_gamma( dffts, psi, dc, dca ) 44 CALL invfft( 'Wave', psi, dffts ) 45 46! on smooth grids--> grids for charge density 47 48 do ir=1, dffts%nnr 49 psi(ir) = CMPLX (kedtaus(ir,iss1)*DBLE(psi(ir)), & 50 kedtaus(ir,iss2)*AIMAG(psi(ir)),kind=DP) 51 end do 52 call fwfft('Wave',psi, dffts ) 53 CALL fftx_psi2c_gamma( dffts, psi, dc, dca ) 54 do ig=1,ngw 55 df(ig)= df(ig) - ci*fi *tpiba2*g(ipol,ig) * dc(ig) 56 da(ig)= da(ig) - ci*fip*tpiba2*g(ipol,ig) * dca(ig) 57 end do 58 end do 59 60 deallocate( dc ) 61 deallocate( dca ) 62! 63 return 64 end subroutine dforce_meta 65!----------------------------------------------------------------------- 66! 67!----------------------------------------------------------------------- 68 subroutine kedtauofr_meta (c) 69!----------------------------------------------------------------------- 70! 71 use kinds, only: dp 72 use control_flags, only: tpre 73 use gvecw, only: ngw 74 use gvect, only: g 75 use cell_base, only : omega, tpiba, ainv 76 use electrons_base, only: nx => nbspx, n => nbsp, f, ispin, nspin 77 use constants, only: pi, fpi 78! 79 use dener 80 use metagga_cp, ONLY : kedtaur, kedtaus, kedtaug, crosstaus, gradwfc, & 81 dkedtaus 82 USE fft_interfaces, ONLY: fwfft, invfft 83 USE fft_base, ONLY: dffts, dfftp 84 USE fft_helper_subroutines, ONLY : c2psi_gamma 85 USE fft_rho 86 87 implicit none 88 89 complex(dp) :: c(ngw,nx) 90 complex(dp), allocatable :: psis( : ) 91 complex(dp), allocatable :: dc( : ), dca( : ) 92 93! local variables 94 integer iss, isup, isdw, iss1, iss2, ios, i, ir, ig 95 integer ipol, ix,iy, ipol2xy(3,3) 96 real(dp) sa1, sa2 97 complex(dp) ci,fp,fm 98! 99 ALLOCATE( psis( dffts%nnr ) ) 100 ALLOCATE( dc( ngw ) ) 101 ALLOCATE( dca( ngw ) ) 102! 103 ci=(0.0d0,1.0d0) 104 kedtaur(:,:)=0.d0 105 kedtaus(:,:)=0.d0 106 kedtaug(:,:)=(0.d0,0.d0) 107 if(tpre) crosstaus(:,:,:)=0.d0 108 109! 110! 111! warning! trhor and thdyn are not compatible yet! 112! 113! important: if n is odd then nx must be .ge.n+1 and c(*,n+1)=0. 114! 115 if (mod(n,2).ne.0) then 116 c(1:ngw,n+1)=(0.d0,0.d0) 117 endif 118 ! 119 do i=1,n,2 120 iss1=ispin(i) 121 sa1=f(i)/omega 122 if (i.ne.n) then 123 iss2=ispin(i+1) 124 sa2=f(i+1)/omega 125 else 126 iss2=iss1 127 sa2=0.0d0 128 end if 129 130 do ipol = 1, 3 131 psis( : ) = (0.d0,0.d0) 132 ! gradient of wfc in real space 133 do ig=1,ngw 134 dc( ig ) = ci * tpiba * g(ipol,ig) * c(ig,i) 135 dca( ig ) = ci * tpiba * g(ipol,ig) * c(ig,i+1) 136 end do 137 CALL c2psi_gamma( dffts, psis, dc, dca ) 138 call invfft('Wave',psis, dffts ) 139 ! on smooth grids--> grids for charge density 140 do ir=1, dffts%nnr 141 kedtaus(ir,iss1)=kedtaus(ir,iss1)+0.5d0*sa1*DBLE(psis(ir))**2 142 kedtaus(ir,iss2)=kedtaus(ir,iss2)+0.5d0*sa2*AIMAG(psis(ir))**2 143 end do 144 if(tpre) then 145 do ir=1, dffts%nnr 146 gradwfc(ir,ipol)=psis(ir) 147 end do 148 end if 149 end do 150 if(tpre) then 151 ipol=1 152 do ix=1,3 153 do iy=1,ix 154 ipol2xy(ix,iy)=ipol 155 ipol2xy(iy,ix)=ipol 156 do ir=1,dffts%nnr 157 crosstaus(ir,ipol,iss1) = crosstaus(ir,ipol,iss1) +& 158 sa1*DBLE(gradwfc(ir,ix))*DBLE(gradwfc(ir,iy)) 159 crosstaus(ir,ipol,iss2) = crosstaus(ir,ipol,iss2) +& 160 sa2*AIMAG(gradwfc(ir,ix))*AIMAG(gradwfc(ir,iy)) 161 end do 162 ipol=ipol+1 163 end do 164 end do 165 end if 166 167 ! d kedtaug / d h 168 if(tpre) then 169 do iss=1,nspin 170 do ix=1,3 171 do iy=1,3 172 do ir=1,dffts%nnr 173 dkedtaus(ir,ix,iy,iss)=-kedtaus(ir,iss)*ainv(iy,ix)& 174 -crosstaus(ir,ipol2xy(1,ix),iss)*ainv(iy,1)& 175 -crosstaus(ir,ipol2xy(2,ix),iss)*ainv(iy,2)& 176 -crosstaus(ir,ipol2xy(3,ix),iss)*ainv(iy,3) 177 end do 178 end do 179 end do 180 end do 181 end if !end metagga 182 ! 183 end do 184 185 DEALLOCATE( dca ) 186 DEALLOCATE( dc ) 187 DEALLOCATE( psis ) 188 189! kinetic energy density (kedtau) in g-space (kedtaug) 190 191 CALL rho_r2g( dffts, kedtaus, kedtaug ) 192! 193 kedtaug(dffts%ngm+1:,:) = 0.0d0 194 195 CALL rho_g2r( dfftp, kedtaug, kedtaur ) 196! 197 return 198 end subroutine kedtauofr_meta 199! 200! 201!----------------------------------------------------------------------- 202 subroutine vofrho_meta ( ) 203!----------------------------------------------------------------------- 204! computes: the one-particle potential v in real space, 205! the total energy etot, 206! the forces fion acting on the ions, 207! the derivative of total energy to cell parameters h 208! rhor input : electronic charge on dense real space grid 209! (plus core charge if present) 210! rhog input : electronic charge in g space (up to density cutoff) 211! rhos input : electronic charge on smooth real space grid 212! rhor output: total potential on dense real space grid 213! rhos output: total potential on smooth real space grid 214! 215 use kinds, only: dp 216 use control_flags, only: thdyn, tpre, tfor, tprnfor 217 use io_global, only: stdout 218 use ions_base, only: nsp, na, nat 219 use cell_base, only: omega 220 use electrons_base, only: nspin 221 use constants, only: pi, fpi 222 use energies, only: etot, eself, enl, ekin, epseu, esr, eht, exc 223 use local_pseudo, only: vps, rhops 224 use core 225 use smallbox_gvec 226 use dener 227 use mp, ONLY : mp_sum 228 use mp_global, ONLY : intra_bgrp_comm 229 use metagga_cp, ONLY : kedtaur, kedtaug, kedtaus, dkedtaus 230 USE fft_interfaces, ONLY: fwfft, invfft 231 USE fft_base, ONLY: dffts, dfftp 232 USE fft_rho 233! 234 implicit none 235! 236 integer iss, isup, isdw, ig, ir,i,j,k,is, ia 237 real(dp) dkedxc(3,3) !metagga 238 complex(dp) fp, fm, ci 239! 240 ci=(0.d0,1.d0) 241! 242! =================================================================== 243! calculation exchange and correlation energy and potential 244! ------------------------------------------------------------------- 245! if (nlcc.gt.0) call add_cc(rhoc,rhog,rhor) 246! 247#ifdef VARIABLECELL 248! call exch_corr_h(nspin,rhog,rhor,exc,dxc) 249#else 250! call exch_corr(nspin,rhog,rhor,exc) 251#endif 252! 253! rhor contains the xc potential in r-space 254! 255! =================================================================== 256! fourier transform of xc potential to g-space (dense grid) 257! ------------------------------------------------------------------- 258! 259 CALL rho_r2g( dfftp, kedtaur, kedtaug ) 260! 261 CALL rho_g2r ( dffts, kedtaug, kedtaus ) 262 263 !calculate dkedxc in real space on smooth grids !metagga 264 if(tpre) then 265 do iss=1,nspin 266 do j=1,3 267 do i=1,3 268 dkedxc(i,j)=0.d0 269 do ir=1,dffts%nnr 270 !2.d0 : because kedtau = 0.5d0 d_Exc/d_kedtau 271 dkedxc(i,j)= dkedxc(i,j)+kedtaus(ir,iss)*2.d0*& 272 dkedtaus(ir,i,j,iss) 273 end do 274 end do 275 end do 276 end do 277 ! 278 call mp_sum( dkedxc, intra_bgrp_comm ) 279 ! 280 do j=1,3 281 do i=1,3 282 dxc(i,j) = dxc(i,j) + omega/(dffts%nr1*dffts%nr2*dffts%nr3)*dkedxc(i,j) 283 end do 284 end do 285 end if 286 return 287 end subroutine vofrho_meta 288!----------------------------------------------------------------------- 289