1! 2! Copyright (C) 2002-2005 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 9subroutine bforceion(fion,tfor,ipol,qmatinv,bec0,becdr,gqq,evalue) 10 11! this subroutine compute the part of force for the ions due to 12! electronic berry phase( see internal notes) 13! it needs becdr 14 15! fion : input, forces on ions 16! tfor : input, if true it computes force 17! at : input, direct lattice vectors, divided by alat 18! ipol : input, electric field polarization 19! qmatinv : input, inverse of Q matrix: Q_i,j=<Psi_i|exp(iG*r)|Psi_j> 20! bec0 : input, factors <beta_iR|Psi_j> 21! becdr : input, factors d<beta_iR|Psi_j>/dR 22! gqq : input, Int_e exp(iG*r)*q_ijR(r) 23! evalue : input, scale of electric field 24 25 use ions_base, only : nax, na, nsp, nat, ityp 26 use uspp_param, only: upf, nh, nhm 27 use uspp, only : nkb, indv_ijkb0 28 use kinds, only : dp 29 use constants, only : 30 use cell_base, only: at, alat 31 use electrons_base, only: nbsp, nbspx, nspin, nbspx_bgrp 32 use mp_global, only: nbgrp 33 34 35 implicit none 36 37 real(dp) evalue 38 complex(dp) qmatinv(nbspx,nbspx),gqq(nhm,nhm,nax,nsp) 39 real(dp) bec0(nkb,nbspx),becdr(nkb,nbspx,3) 40 real(dp) fion(3,*) 41 integer ipol 42 logical tfor 43 44!local variables 45 46 complex(dp) ci, temp, temp1,temp2,temp3 47 real(dp) :: gmes 48 real(dp), external :: g_mes 49 integer iv,jv,ia,is,k,i,j,ilm,jlm,inl,jnl,ism 50 51 if(.not. tfor) return 52 53 if( nbgrp > 1 ) & 54 call errore(' bforceion ', ' parallelization over bands not yet implemented ', 1 ) 55 56 ci = (0.d0,1.d0) 57 gmes = g_mes (ipol, at, alat) 58 59 do ia=1,nat 60 is=ityp(ia) 61 IF(upf(is)%tvanp) THEN 62 do iv= 1,nh(is) 63 do jv=1,nh(is) 64 inl = indv_ijkb0(ia) + iv 65 jnl = indv_ijkb0(ia) + jv 66 temp=(0.d0,0.d0) 67 temp1=(0.d0,0.d0) 68 temp2=(0.d0,0.d0) 69 temp3=(0.d0,0.d0) 70 do i=1,nbsp 71 do j=1,nbsp 72 73 temp = temp + ci*gmes*gqq(iv,jv,ia,is)* &!TAKECARE: sign + due to exp(+iGr) in gqq 74 & bec0(inl,i)*bec0(jnl,j)*qmatinv(j,i) 75 76 temp1 = temp1 + gqq(iv,jv,ia,is)*& 77 & ( becdr(inl,i,1)*bec0(jnl,j)+bec0(inl,i)*becdr(jnl,j,1))*qmatinv(j,i) 78 79 temp2 = temp2 + gqq(iv,jv,ia,is)*& 80 & ( becdr(inl,i,2)*bec0(jnl,j)+bec0(inl,i)*becdr(jnl,j,2))*qmatinv(j,i) 81 82 temp3 = temp3 + gqq(iv,jv,ia,is)*& 83 & ( becdr(inl,i,3)*bec0(jnl,j)+bec0(inl,i)*becdr(jnl,j,3))*qmatinv(j,i) 84 85 86 enddo 87 enddo 88 89 fion(ipol,ia) = fion(ipol,ia) - 2.d0*evalue*AIMAG(temp)/gmes 90 fion(1,ia) = fion(1,ia) - 2.d0*evalue*AIMAG(temp1)/gmes 91 fion(2,ia) = fion(2,ia) - 2.d0*evalue*AIMAG(temp2)/gmes 92 fion(3,ia) = fion(3,ia) - 2.d0*evalue*AIMAG(temp3)/gmes 93 end do 94 end do 95 END IF 96 end do 97 98 return 99end subroutine bforceion 100