1! 2! Copyright (C) 2002-2005 FPMD-CPV 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! 8MODULE efield_module 9 10 USE kinds, ONLY : DP 11 12 IMPLICIT NONE 13 SAVE 14 15 logical :: tefield = .FALSE. 16 logical :: tefield2 = .FALSE. 17 integer :: epol = 3 !direction electric field 18 real(kind=DP) :: efield = 0.d0 !intensity electric field 19 real(kind=DP) :: efield2 =0.d0 20 real(kind=DP) evalue!strenght of electric field 21 real(kind=DP) evalue2 22 integer epol2,ipolp2 23 integer ipolp !direction of electric field 24 25 real(kind=DP) :: pberryel = 0.0d0, pberryion = 0.0d0 26 real(kind=DP) :: pberryel2 = 0.0d0, pberryion2 = 0.0d0 27 28!*** 29!*** Berry phase 30!*** 31 integer, allocatable:: ctable(:,:,:)!correspondence tables for diff. polarization 32 integer, allocatable:: ctabin(:,:,:)!inverse correspondence table 33 complex(DP), allocatable:: qmat(:,:)!inverse of matrix Q, for Barry's phase 34 complex(DP), allocatable:: gqq(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr 35 complex(DP), allocatable:: gqqm(:,:,:,:)! the same with exp(-iGr) 36 complex(DP), allocatable:: gqq0(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr, at Gamma 37 complex(DP), allocatable:: gqqm0(:,:,:,:)! the same with exp(-iGr), at Gamma 38 complex(DP), allocatable:: df(:) 39 integer, allocatable:: ctable2(:,:,:)!correspondence tables for diff. polarization 40 integer, allocatable:: ctabin2(:,:,:)!inverse correspondence table 41 complex(DP), allocatable:: qmat2(:,:)!inverse of matrix Q, for Barry's phase 42 complex(DP), allocatable:: gqq2(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr 43 complex(DP), allocatable:: gqqm2(:,:,:,:)! the same with exp(-iGr) 44 complex(DP), allocatable:: gqq02(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr, at Gamma 45 complex(DP), allocatable:: gqqm02(:,:,:,:)! the same with exp(-iGr), at Gamma 46 complex(DP) detq 47 complex(DP) detq2 48 real(DP) cdzp(3),cdzm(3), cdz0(3)!centers of ionic charges 49 50!for parallelization for direcions 1 and 2 51 integer :: n_g_missing_p(2)!number of g vector with correspondence G-->G+1 is missing 52 integer :: n_g_missing_m(2)!number of g vector with correspondence G-->G-1 is missing 53 integer, allocatable :: whose_is_g(:) !correspondence G(plane waves, global) ---> processor 54 integer, allocatable :: ctable_missing_1(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall 55 !n_g_missing*nproc 56 integer, allocatable :: ctable_missing_rev_1(:,:,:)!missing_g --> G (plane waves local) 57 integer, allocatable :: ctable_missing_2(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall 58 !n_g_missing*nproc 59 integer, allocatable :: ctable_missing_rev_2(:,:,:)!missing_g --> G (plane waves local) 60 integer, allocatable :: ctabin_missing_1(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall 61 !n_g_missing*nproc 62 integer, allocatable :: ctabin_missing_rev_1(:,:,:)!missing_g --> G (plane waves local) 63 integer, allocatable :: ctabin_missing_2(:,:,:)!correspondence G(plane waves local)--> array for mpi_alltoall 64 !n_g_missing*nproc 65 integer, allocatable :: ctabin_missing_rev_2(:,:,:)!missing_g --> G (plane waves local) 66 67CONTAINS 68 69 70 SUBROUTINE efield_init( epol_ , efield_ ) 71 USE kinds, ONLY: DP 72 REAL(DP), INTENT(IN) :: efield_ 73 INTEGER, INTENT(IN) :: epol_ 74 epol = epol_ 75 efield = efield_ 76 RETURN 77 END SUBROUTINE efield_init 78 79 SUBROUTINE efield_info( ) 80 USE io_global, ONLY: ionode,stdout 81 if(ionode) write (stdout,401) epol, efield 82 83401 format (/4x,'=====================================' & 84 & /4x,'| BERRY PHASE ELECTRIC FIELD 1 ' & 85 & /4x,'=====================================' & 86 & /4x,'| direction =',i10,' ' & 87 & /4x,'| intensity =',f10.5,' a.u. ' & 88 & /4x,'=====================================') 89 90 RETURN 91 END SUBROUTINE efield_info 92 93 94 SUBROUTINE efield_berry_setup( eigr, tau0 ) 95 USE io_global, ONLY: ionode,stdout 96 IMPLICIT NONE 97 COMPLEX(DP), INTENT(IN) :: eigr(:,:) 98 REAL(DP), INTENT(IN) :: tau0(:,:) 99 if(ionode) write(stdout,'("Initialize Berry phase electric field")') 100 ipolp = epol 101 evalue = efield 102!set up for parallel calculations 103 104#if defined(__MPI) 105 call find_whose_is_g 106 call gtable_missing 107 call gtable_missing_inv 108#endif 109 110 call gtable(ipolp,ctable(1,1,ipolp)) 111 call gtablein(ipolp,ctabin(1,1,ipolp)) 112 call qqberry2(gqq0,gqqm0,ipolp)!for Vanderbilt pps 113 call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp) 114 !the following line was to keep the center of charge fixed 115 !when performing molecular dynamics in the presence of an electric 116 !field 117 !call cofcharge(tau0,cdz0) 118 119 RETURN 120 END SUBROUTINE efield_berry_setup 121 122 123 SUBROUTINE efield_update( eigr ) 124 IMPLICIT NONE 125 COMPLEX(DP), INTENT(IN) :: eigr(:,:) 126 call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp) 127 RETURN 128 END SUBROUTINE efield_update 129 130 131 SUBROUTINE allocate_efield( ngw, ngw_g, nx, nhx, nax, nsp ) 132 IMPLICIT NONE 133 INTEGER, INTENT(IN) :: ngw, ngw_g, nx, nhx, nax, nsp 134 allocate( ctable(ngw,2,3)) 135 allocate( ctabin(ngw,2,3)) 136 allocate( qmat(nx,nx)) 137 allocate( gqq(nhx,nhx,nax,nsp)) 138 allocate( gqqm(nhx,nhx,nax,nsp)) 139 allocate( df(ngw)) 140 allocate( gqq0(nhx,nhx,nax,nsp)) 141 allocate( gqqm0(nhx,nhx,nax,nsp)) 142 allocate( whose_is_g(ngw_g)) 143 144 RETURN 145 END SUBROUTINE allocate_efield 146 147 148 SUBROUTINE deallocate_efield( ) 149 IMPLICIT NONE 150 IF( allocated( ctable ) ) deallocate( ctable ) 151 IF( allocated( ctabin ) ) deallocate( ctabin ) 152 IF( allocated( qmat ) ) deallocate( qmat ) 153 IF( allocated( gqq ) ) deallocate( gqq ) 154 IF( allocated( gqqm ) ) deallocate( gqqm ) 155 IF( allocated( df ) ) deallocate( df ) 156 IF( allocated( gqq0 ) ) deallocate( gqq0 ) 157 IF( allocated( gqqm0 ) ) deallocate( gqqm0 ) 158 IF( allocated( whose_is_g) ) deallocate(whose_is_g) 159 IF( allocated( ctable_missing_1) ) deallocate( ctable_missing_1) 160 IF( allocated( ctable_missing_2) ) deallocate( ctable_missing_2) 161 IF( allocated( ctable_missing_rev_1) ) deallocate( ctable_missing_rev_1) 162 IF( allocated( ctable_missing_rev_1) ) deallocate( ctable_missing_rev_2) 163 IF( allocated( ctabin_missing_1) ) deallocate( ctabin_missing_1) 164 IF( allocated( ctabin_missing_2) ) deallocate( ctabin_missing_2) 165 IF( allocated( ctabin_missing_rev_1) ) deallocate( ctabin_missing_rev_1) 166 IF( allocated( ctabin_missing_rev_1) ) deallocate( ctabin_missing_rev_2) 167 RETURN 168 END SUBROUTINE deallocate_efield 169 170 171 SUBROUTINE berry_energy( enb, enbi, bec, cm, fion ) 172 USE ions_positions, ONLY: tau0 173 USE control_flags, ONLY: tfor, tprnfor 174 IMPLICIT NONE 175 real(DP), intent(out) :: enb, enbi 176 real(DP) :: bec(:,:) 177 real(DP) :: fion(:,:) 178 complex(DP) :: cm(:,:) 179 real(dp), external :: enberry 180 181 call qmatrixd(cm,bec,ctable(1,1,ipolp),gqq,qmat,detq,ipolp) 182 enb = enberry( detq, ipolp ) 183 call berryion(tau0,fion,tfor.or.tprnfor,ipolp,evalue,enbi) 184 pberryel=enb 185 pberryion=enbi 186 enb=enb*evalue 187 enbi=enbi*evalue 188 END SUBROUTINE berry_energy 189 190 191 SUBROUTINE dforce_efield (bec,i,cm,c2,c3,rhos) 192 USE uspp, ONLY: betae => vkb, deeq 193 USE gvecw, ONLY: ngw 194 IMPLICIT NONE 195 complex(DP), intent(out) :: c2(:), c3(:) 196 complex(DP), intent(in) :: cm(:,:) 197 REAL(DP) :: rhos(:,:) 198 real(DP) :: bec(:,:) 199 integer :: i 200 integer :: ig 201 call dforceb (cm, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) 202 do ig=1,ngw 203 c2(ig)=c2(ig)+evalue*df(ig) 204 enddo 205 call dforceb (cm, i+1, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) 206 do ig=1,ngw 207 c3(ig)=c3(ig)+evalue*df(ig) 208 enddo 209 END SUBROUTINE dforce_efield 210 211 SUBROUTINE efield_init2( epol_ , efield_ ) 212 USE kinds, ONLY: DP 213 REAL(DP), INTENT(IN) :: efield_ 214 INTEGER, INTENT(IN) :: epol_ 215 epol2 = epol_ 216 efield2 = efield_ 217 RETURN 218 END SUBROUTINE efield_init2 219 220 SUBROUTINE efield_info2( ) 221 USE io_global, ONLY: ionode,stdout 222 if(ionode) write (stdout,402) epol2, efield2 223 224402 format (/4x,'=====================================' & 225 & /4x,'| BERRY PHASE ELECTRIC FIELD 2 ' & 226 & /4x,'=====================================' & 227 & /4x,'| direction =',i10,' ' & 228 & /4x,'| intensity =',f10.5,' a.u. ' & 229 & /4x,'=====================================') 230 231 RETURN 232 END SUBROUTINE efield_info2 233 234 235 SUBROUTINE efield_berry_setup2( eigr, tau0 ) 236 USE io_global, ONLY: ionode,stdout 237 IMPLICIT NONE 238 COMPLEX(DP), INTENT(IN) :: eigr(:,:) 239 REAL(DP), INTENT(IN) :: tau0(:,:) 240 if(ionode) write(stdout,'("Initialize Berry phase electric field")') 241 ipolp2 = epol2 242 evalue2 = efield2 243 call gtable(ipolp2,ctable2(1,1,ipolp2)) 244 call gtablein(ipolp2,ctabin2(1,1,ipolp2)) 245 call qqberry2(gqq02,gqqm02,ipolp2)!for Vanderbilt pps 246 call qqupdate(eigr,gqqm02,gqq2,gqqm2,ipolp2) 247 !the following line was to keep the center of charge fixed 248 !when performing molecular dynamics in the presence of an electric 249 !field 250 !call cofcharge(tau0,cdz0) 251 RETURN 252 END SUBROUTINE efield_berry_setup2 253 254 255 SUBROUTINE efield_update2( eigr ) 256 IMPLICIT NONE 257 COMPLEX(DP), INTENT(IN) :: eigr(:,:) 258 call qqupdate(eigr,gqqm02,gqq2,gqqm2,ipolp2) 259 RETURN 260 END SUBROUTINE efield_update2 261 262 263 SUBROUTINE allocate_efield2( ngw, nx, nhx, nax, nsp ) 264 IMPLICIT NONE 265 INTEGER, INTENT(IN) :: ngw, nx, nhx, nax, nsp 266 allocate( ctable2(ngw,2,3)) 267 allocate( ctabin2(ngw,2,3)) 268 allocate( qmat2(nx,nx)) 269 allocate( gqq2(nhx,nhx,nax,nsp)) 270 allocate( gqqm2(nhx,nhx,nax,nsp)) 271 allocate( gqq02(nhx,nhx,nax,nsp)) 272 allocate( gqqm02(nhx,nhx,nax,nsp)) 273 RETURN 274 END SUBROUTINE allocate_efield2 275 276 277 SUBROUTINE deallocate_efield2( ) 278 IMPLICIT NONE 279 IF( allocated( ctable2 ) ) deallocate( ctable2 ) 280 IF( allocated( ctabin2 ) ) deallocate( ctabin2 ) 281 IF( allocated( qmat2 ) ) deallocate( qmat2 ) 282 IF( allocated( gqq2 ) ) deallocate( gqq2 ) 283 IF( allocated( gqqm2 ) ) deallocate( gqqm2 ) 284 IF( allocated( gqq02 ) ) deallocate( gqq02 ) 285 IF( allocated( gqqm02 ) ) deallocate( gqqm02 ) 286 RETURN 287 END SUBROUTINE deallocate_efield2 288 289 290 SUBROUTINE berry_energy2( enb, enbi, bec, cm, fion ) 291 USE ions_positions, ONLY: tau0 292 USE control_flags, ONLY: tfor, tprnfor 293 IMPLICIT NONE 294 real(DP), intent(out) :: enb, enbi 295 real(DP) :: bec(:,:) 296 real(DP) :: fion(:,:) 297 complex(DP) :: cm(:,:) 298 real(dp), external :: enberry 299 300 call qmatrixd(cm,bec,ctable2(1,1,ipolp2),gqq2,qmat2,detq2,ipolp2) 301 enb = enberry( detq2, ipolp2 ) 302 call berryion(tau0,fion,tfor.or.tprnfor,ipolp2,evalue2,enbi) 303 pberryel2=enb 304 pberryion2=enbi 305 enb=enb*evalue2 306 enbi=enbi*evalue2 307 END SUBROUTINE berry_energy2 308 309 310 SUBROUTINE dforce_efield2 (bec,i,cm,c2,c3,rhos) 311 USE uspp, ONLY: betae => vkb, deeq 312 USE gvecw, ONLY: ngw 313 IMPLICIT NONE 314 complex(DP), intent(out) :: c2(:), c3(:) 315 complex(DP), intent(in) :: cm(:,:) 316 REAL(DP) :: rhos(:,:) 317 real(DP) :: bec(:,:) 318 integer :: i 319 integer :: ig 320 call dforceb (cm, i, betae, ipolp2, bec ,ctabin2(1,1,ipolp2), gqq2, gqqm2, qmat2, deeq, df) 321 do ig=1,ngw 322 c2(ig)=c2(ig)+evalue2*df(ig) 323 enddo 324 call dforceb (cm, i+1, betae, ipolp2, bec ,ctabin2(1,1,ipolp2), gqq2, gqqm2, qmat2, deeq, df) 325 do ig=1,ngw 326 c3(ig)=c3(ig)+evalue2*df(ig) 327 enddo 328 END SUBROUTINE dforce_efield2 329 330END MODULE efield_module 331