! this subroutine returns the berry phase energy
! = L/2*Pi*Im(log Sum_R exp(i*(2pi/L)*R_i*rho_i))
! of the ions and the constant force on the ions
! now only for orthorombic primitive cell

subroutine berryion( tau0,fion, tfor,ipol,evalue,enbi)

!  tau0    : input, positions of ions
!  fion    : input,output, forces on ions
!  tfor    : input, flag for force calculation
!  ipol    : input, electric field polarization
!  evalue  : input, scale for electric field
!  enbi    : output, berry phase energy of the ions

     

  use ions
  use parm
  use elct
  use cnst

  implicit none

  real(kind=8) tau0(3,nax,nsx)
  real(kind=8) fion(3,nax,nsx)
  real(kind=8) enbi, evalue
  integer ipol
  logical tfor

!local variables
  real(kind=8) gmes
  real(kind=8) pola
  integer is, ia
  complex(kind=8) temp, ci

  temp = (0.,0.)
  ci = (0.,1.)

  if(ipol.eq.1) then
     gmes=a1(1)**2+a1(2)**2+a1(3)**2
     gmes=2*pi/SQRT(gmes)
  endif
  if(ipol.eq.2) then
     gmes=a2(1)**2+a2(2)**2+a2(3)**2
     gmes=2*pi/SQRT(gmes)
  endif
  if(ipol.eq.3) then
     gmes=a3(1)**2+a3(2)**2+a3(3)**2
     gmes=2*pi/SQRT(gmes)
  endif
  pola=0.d0
  do is=1,nsp
     do ia=1,na(is)

!this force term is along ipol-direction
        if( tfor) then
           fion(ipol,ia,is)=fion(ipol,ia,is)+evalue*zv(is)
        endif
            
        temp = temp - ci*gmes*tau0(ipol,ia,is)*zv(is)
        pola=pola+evalue*zv(is)*tau0(ipol,ia,is)!this is just the center of ionic charge
     enddo
  enddo

  enbi=aimag(log(exp(temp)))/gmes!this sounds stupid it's just a Riemann plane
!  write(6,*) 'Pola  :', pola!ATTENZIONE
  return
end subroutine berryion

            
!-------------------------------------------------------------------------
      subroutine cofcharge(tau,cdz)
!-----------------------------------------------------------------------
!this subroutine gives the center of the ionic charge


      use ions
!
      implicit none
      real(kind=8) tau(3,nax,nsp), cdz(3)
! local variables
      real(kind=8) zmas
      integer is,i,ia
!
      zmas=0.0
      do is=1,nsp
         zmas=zmas+na(is)*zv(is)
      end do
!
      do i=1,3
         cdz(i)=0.0
         do is=1,nsp
            do ia=1,na(is)
               cdz(i)=cdz(i)+tau(i,ia,is)*zv(is)
            end do
         end do
         cdz(i)=cdz(i)/zmas
      end do
      write(6,*) 'Center of charge', cdz(3)!ATTENZIONE
!
      return
      end
!


!-----------------------------------------------
      subroutine keptwfc( c0, t1, t2, t3)
!-----------------------------------------------
! this subroutine translate the electronic wfcs
!   of (t1,t2,t3) in order to keep the center steady
! c0, input/output: electronic wfc
! t1,t2,t3: traslation vector


       
        use ions
        use elct, only:ngw,n
        use parm
        use cnst
        use gvec
        use control
        implicit none

        complex(kind=8) c0(ngw,n), t1, t2, t3

       
        complex(kind=8) ei1(-nr1:nr1), ei2(-nr2:nr2),      &
     &                ei3(-nr3:nr3)
        integer i,j,k, ig
        real(kind=8) taup(3), sum, ar1,ar2,ar3, tau(3)
        complex(kind=8) ctep1,ctep2,ctep3,ctem1,ctem2,ctem3
        
           
        tau(1)=t1
        tau(2)=t2
        tau(3)=t3
        
        do i=1,3
           sum=0.d0
           do j=1,3
              sum=sum+ainv(i,j)*tau(j)
           end do
           taup(i)=sum
!
! tau0=x1*a1+x2*a2+x3*a3 => taup(1)=x1=tau0*b1 and so on
!
        end do
!
        ar1=2.d0*pi*taup(1)
        ctep1=cmplx(cos(ar1),-sin(ar1))
        ctem1=conjg(ctep1)
        ei1( 0)=cmplx(1.d0,0.d0)
        ei1( 1)=ctep1
        ei1(-1)=ctem1
        do i=2,nr1-1
           ei1( i)=ei1( i-1)*ctep1
           ei1(-i)=ei1(-i+1)*ctem1
        end do
        !
        ar2=2.d0*pi*taup(2)
        ctep2=cmplx(cos(ar2),-sin(ar2))
        ctem2=conjg(ctep2)
        ei2( 0)=cmplx(1.d0,0.d0)
        ei2( 1)=ctep2
        ei2(-1)=ctem2
        do j=2,nr2-1
           ei2( j)=ei2( j-1)*ctep2
           ei2(-j)=ei2(-j+1)*ctem2
        end do
            !
        ar3=2.d0*pi*taup(3)
        ctep3=cmplx(cos(ar3),-sin(ar3))
        ctem3=conjg(ctep3)
        ei3( 0)=cmplx(1.d0,0.d0)
        ei3( 1)=ctep3
        ei3(-1)=ctem3
        do k=2,nr3-1
           ei3( k)=ei3( k-1)*ctep3
           ei3(-k)=ei3(-k+1)*ctem3
        end do

        do i=1,n
           do ig=1,ngw
              c0(ig,i)=c0(ig,i)*ei1(in1p(ig))*ei2(in2p(ig))*ei3(in3p(ig))
           end do
        end do

        return
        end

!----------------------------------------------------
        subroutine noforce(fion, ipol)
!----------------------------------------------------

! this subroutine adds a electric force, in order
! to keep steady the center of mass along the electric
! field direction

          use ions

          implicit none

          real(kind=8) fion(3,nax,nsx)
          integer ipol!el. field polarization


          integer i,ia,is
          real(kind=8) fcm!force appplied on center of mass
          real(kind=8) tch!total charge

          fcm=0.d0
          tch=0.d0
          do is=1,nsp
             do ia=1,na(is)
                fcm=fcm+fion(ipol,ia,is)     
                tch=tch+zv(is)
             enddo             
          enddo
          write(6,*) 'Forza su ioni in ipol:', fcm!ATTENZIONE 
          fcm=fcm/tch
          do is=1,nsp
             do ia=1,na(is)
                fion(ipol,ia,is)=fion(ipol,ia,is)-fcm*zv(is)
             enddo
          enddo
   
          return
        end subroutine noforce