1! 2! Copyright (C) 2002-2008 Quantum ESPRESS0 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! 8subroutine dforceb(c0, i, betae, ipol, bec0, ctabin, gqq, gqqm, qmat, dq2, df) 9 10! this subroutine computes the force for electrons 11! in case of Berry,s phase like perturbation 12! it gives the force for the i-th state 13 14! c0 input: Psi^0_i 15! c1 input: Psi^1_i 16! i input: ot computes the force for the i-th state 17! v0 input: the local zeroth order potential 18! v1 input: the local first order potential 19! betae input: the functions beta_iR 20! ipol input:the polarization of nabla_k 21! bec0 input: the factors <beta_iR|Psi^0_v> 22! bec1 input: the factors <beta_iR|Psi^1_v> 23! ctabin input: the inverse-correspondence array g'+(-)1=g 24! gqq input: the factors int dr Beta_Rj*Beta_Ri exp(iGr) 25! gqqm input: the factors int dr Beta_Rj*Beta_Ri exp(iGr) 26! qmat input: 27! dq2 input: factors d^2hxc_ijR 28! df output: force for the i-th state 29 30 31 use kinds, only : DP 32 use gvecw, only: ngw 33 use parameters 34 use electrons_base, only: nx => nbspx, n => nbsp, nspin, f 35 use constants, only : 36 use ions_base, only : nat, nax, nsp, ityp 37 use cell_base, only: at, alat 38 use uspp_param, only: nh, nhm, upf 39 use uspp, only : nkb, nkbus, indv_ijkb0 40 use efield_module, ONLY : ctabin_missing_1,ctabin_missing_2,n_g_missing_m,& 41 & ctabin_missing_rev_1,ctabin_missing_rev_2 42 use mp_global, only: intra_bgrp_comm, nproc_bgrp 43 use mp, only: mp_alltoall 44 use parallel_include 45 46 47 implicit none 48 49 50 complex(DP) c0(ngw, n), betae(ngw,nkb), df(ngw),& 51 & gqq(nhm,nhm,nax,nsp),gqqm(nhm,nhm,nax,nsp),& 52 & qmat(nx,nx) 53 real(DP) bec0(nkb,n), dq2(nat,nhm,nhm,nspin), gmes 54 real(DP), EXTERNAL :: g_mes 55 56 integer i, ipol, ctabin(ngw,2) 57 58! local variables 59 60 integer j,k,ig,iv,jv,ix,jx,is,ia, iss,iss1,mism 61 integer ir,ism,itemp,itempa,jnl,inl 62 complex(DP) ci ,fi, fp, fm 63 real(DP) afr(nkb), dd 64 complex(DP) afrc(nkb) 65 complex(DP), allocatable:: dtemp(:) 66 complex(DP), allocatable :: sndbuf(:,:,:),rcvbuf(:,:,:) 67 integer :: ierr, ip 68 69 allocate( dtemp(ngw)) 70 71 72 ci=(0.d0,1.d0) 73 74 75! now the interaction term 76! first the norm-conserving part 77 78 do ig=1,ngw 79 dtemp(ig)=(0.d0,0.d0) 80 enddo 81 82 do j=1,n 83 do ig=1,ngw 84 if(ctabin(ig,2) .ne. (ngw+1)) then 85 if(ctabin(ig,2).ge.0) then 86 dtemp(ig)=dtemp(ig)+c0(ctabin(ig,2),j)*qmat(j,i) 87 else 88 dtemp(ig)=dtemp(ig)+CONJG(c0(-ctabin(ig,2),j))*qmat(j,i) 89 endif 90 endif 91 enddo 92 do ig=1,ngw 93 if(ctabin(ig,1) .ne. (ngw+1)) then 94 if(ctabin(ig,1).ge.0) then 95 dtemp(ig)=dtemp(ig)-c0(ctabin(ig,1),j)*CONJG(qmat(j,i)) 96 else 97 dtemp(ig)=dtemp(ig)-CONJG(c0(-ctabin(ig,1),j))*conjg(qmat(j,i)) 98 endif 99 endif 100 enddo 101 102#if defined(__MPI) 103 104 if(ipol/=3) then 105!allocate arrays 106 allocate(sndbuf(n_g_missing_m(ipol),2,nproc_bgrp)) 107 sndbuf(:,:,:)=(0.d0,0.d0) 108 allocate(rcvbuf(n_g_missing_m(ipol),2,nproc_bgrp)) 109!copy arrays to snd buf 110 do ip=1,nproc_bgrp 111 do ig=1,n_g_missing_m(ipol) 112 if(ipol==1) then 113 if(ctabin_missing_rev_1(ig,1,ip)>0) then 114 sndbuf(ig,1,ip)=-c0(ctabin_missing_rev_1(ig,1,ip),j)*CONJG(qmat(j,i)) 115 elseif(ctabin_missing_rev_1(ig,1,ip)<0) then 116 sndbuf(ig,1,ip)=-conjg(c0(-ctabin_missing_rev_1(ig,1,ip),j))*CONJG(qmat(j,i)) 117 endif 118 else 119 if(ctabin_missing_rev_2(ig,1,ip)>0) then 120 sndbuf(ig,1,ip)=-c0(ctabin_missing_rev_2(ig,1,ip),j)*CONJG(qmat(j,i)) 121 elseif(ctabin_missing_rev_2(ig,1,ip)<0) then 122 sndbuf(ig,1,ip)=-conjg(c0(-ctabin_missing_rev_2(ig,1,ip),j))*CONJG(qmat(j,i)) 123 endif 124 endif 125 enddo 126 do ig=1,n_g_missing_m(ipol) 127 if(ipol==1) then 128 if(ctabin_missing_rev_1(ig,2,ip)>0) then 129 sndbuf(ig,2,ip)=c0(ctabin_missing_rev_1(ig,2,ip),j)*qmat(j,i) 130 elseif(ctabin_missing_rev_1(ig,2,ip)<0) then 131 sndbuf(ig,2,ip)=conjg(c0(-ctabin_missing_rev_1(ig,2,ip),j))*qmat(j,i) 132 endif 133 else 134 if(ctabin_missing_rev_2(ig,2,ip)>0) then 135 sndbuf(ig,2,ip)=c0(ctabin_missing_rev_2(ig,2,ip),j)*qmat(j,i) 136 elseif(ctabin_missing_rev_2(ig,2,ip)<0) then 137 sndbuf(ig,2,ip)=conjg(c0(-ctabin_missing_rev_2(ig,2,ip),j))*qmat(j,i) 138 endif 139 endif 140 enddo 141 enddo 142 143 144 CALL mp_alltoall( sndbuf, rcvbuf, intra_bgrp_comm ) 145 146!update sca 147 do ip=1,nproc_bgrp 148 do ig=1,n_g_missing_m(ipol) 149 if(ipol==1) then 150 if(ctabin_missing_1(ig,1,ip)/=0) then 151 dtemp(ctabin_missing_1(ig,1,ip))=dtemp(ctabin_missing_1(ig,1,ip))+rcvbuf(ig,1,ip) 152 endif 153 if(ctabin_missing_1(ig,2,ip)/=0) then 154 dtemp(ctabin_missing_1(ig,2,ip))=dtemp(ctabin_missing_1(ig,2,ip))+rcvbuf(ig,2,ip) 155 endif 156 else 157 if(ctabin_missing_2(ig,1,ip)/=0) then 158 dtemp(ctabin_missing_2(ig,1,ip))=dtemp(ctabin_missing_2(ig,1,ip))+rcvbuf(ig,1,ip) 159 endif 160 if(ctabin_missing_2(ig,2,ip)/=0) then 161 dtemp(ctabin_missing_2(ig,2,ip))=dtemp(ctabin_missing_2(ig,2,ip))+rcvbuf(ig,2,ip) 162 endif 163 endif 164 enddo 165 enddo 166 167 168 169 170 171 172 173!deallocate arrays 174 deallocate(rcvbuf,sndbuf) 175 endif 176 177#endif 178 enddo 179 180 gmes = g_mes ( ipol, at, alat ) 181 182 fi=f(i)*ci/(2.d0*gmes) 183 184 do ig=1,ngw 185 df(ig)= fi*dtemp(ig) 186 end do 187 188! now the interacting Vanderbilt term 189! the term (-ie/|G|)(-beta_i'R>gqq(i',j')bec0_jRj'Q^-1_ji+ 190! +beta_i'R>gqqm(i',j')bec0jRj'Q^-1_ij* 191 192 193 194 if(nkb > 0) then 195 do inl=1,nkb 196 afrc(inl)=(0.d0,0.d0) 197 end do 198 199 do ia=1,nat 200 is=ityp(ia) 201 IF(upf(is)%tvanp) THEN 202 do iv=1,nh(is) !loop on projectors 203 do jv=1,nh(is) !loop on projectors 204 inl = indv_ijkb0(ia) + iv 205 jnl = indv_ijkb0(ia) + jv 206 do j=1,n !loop on states 207 afrc(inl)=afrc(inl)+gqq(iv,jv,ia,is)*bec0(jnl,j)*qmat(j,i)& 208 & -CONJG(gqq(jv,iv,ia,is))*bec0(jnl,j)*conjg(qmat(i,j)) 209 210 211 end do 212 end do 213 end do 214 ENDIF 215 end do 216 217 do ig=1,ngw 218 dtemp(ig)=(0.d0,0.d0) 219 end do 220 do inl=1,nkb 221 do ig=1,ngw 222 dtemp(ig)=dtemp(ig)+afrc(inl)*betae(ig,inl) 223 enddo 224 enddo 225! call MXMA 226! & (betae,1,2*ngw,afr,1,nhsax,dtemp,1,2*ngw,2*ngw,nhsa,1) 227 do ig=1,ngw 228 df(ig)=df(ig)+fi*dtemp(ig) 229 end do 230 endif 231 232 deallocate( dtemp) 233 return 234 end subroutine dforceb 235 236 237 238function enberry( detq, ipol ) 239 240 use constants, only : 241 use kinds, only: dp 242 use cell_base, only: alat, at 243 USE electrons_base, ONLY : nspin 244 245 implicit none 246 247 complex(dp), intent (in) :: detq 248 real(dp) :: enberry 249 250 integer ipol 251 real(dp) gmes 252 real(dp), external :: g_mes 253 254 gmes = g_mes ( ipol, at, alat ) 255 enberry = 2.d0/REAL(nspin,DP)*AIMAG(log(detq))/gmes ! take care of sign 256 257 return 258 end function enberry 259 260 261! 262! Copyright (C) 2011 Quantum ESPRESSO group 263! This file is distributed under the terms of the 264! GNU General Public License. See the file `License' 265! in the root directory of the present distribution, 266! or http://www.gnu.org/copyleft/gpl.txt . 267! 268 269FUNCTION g_mes ( ipol, at, alat ) 270 271 USE kinds, ONLY : dp 272 USE constants, ONLY : pi 273 274 IMPLICIT NONE 275 276 INTEGER, INTENT(IN) :: ipol 277 REAL(dp), INTENT(IN) :: at(3,3), alat 278 REAL(dp) :: g_mes 279 280 IF ( ipol < 1 .OR. ipol > 3) CALL errore ( 'gmes','incorrect ipol', 1) 281 g_mes = 2.0_dp*pi/alat/SQRT(at(1,ipol)**2+at(2,ipol)**2+at(3,ipol)**2) 282 283END FUNCTION g_mes 284