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