Logo Search packages:      
Sourcecode: abinit version File versions  Download package

berryphase.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/berryphase
!! NAME
!! berryphase
!!
!!
!! FUNCTION
!! This routine is called in scfcv.f to compute the electronic Berry Phase
!! polarization and the ionic contribution to the polarization
!! Work for nsppol=1 or 2 ,but only accept nspinor=1, and mkmem=nkpt
!! or 0, kptopt = 2 or 3
!!
!! COPYRIGHT
!! Copyright (C) 2000-2007 ABINIT  group (NSAI,XG,MKV)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! INPUTS
!! atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f)
!! bdberry(4)=band limits for Berry phase contributions,
!!  spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
!! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
!! gprimd(3,3)=reciprocal space dimensional primitive translations
!! istwfk(nkpt_)=input option parameter that describes the storage of wfs
!! kberry(3,nberry)= different delta k for Berry phases, in unit of kptrlatt
!!  only kberry(1:3,1:nberry) is relevant
!! kg(3,mpw*mkmem)=reduced planewave coordinates
!! kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
!!               kpt_ sampels half the BZ if time-reversal symetrie is used
!! kptopt=2 when time-reversal symetrie is used
!!       =3 when time-reversal symetrie is not used
!! kptrlatt(3,3)=k-point lattice specification
!! mband=maximum number of bands
!! mgfft=maximum size of 1D FFTs
!! mkmem=number of k points which can fit in memory; set to 0 if use disk
!! mpw=maximum dimensioned size of npw
!! natom=number of atoms in cell
!! nattyp(ntypat)= # atoms of each type.
!! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
!! nberry=number of Berry phases to be computed
!! nkpt=number of k points
!! npwarr(nkpt)=number of planewaves in basis at this k point
!! nspinor=number of spinorial components of the wavefunctions
!! nsppol=1 for unpolarized, 2 for spin-polarized
!! ntypat=number of types of atoms in unit cell
!! nkpt_=number of k points generated by ABINIT, (see kpt_)
!! rprimd(3,3)=dimensional real space primitive translations (bohr)
!! typat(natom)=type integer for each atom in cell
!! ucvol=unit cell volume in bohr**3.
!! unkg=unit number for (k+G) basis sphere data
!! wffnow=struct info for wf disk file
!! xred(3,natom)=reduced atomic coordinates
!! zion(ntypat)=valence charge of each type of atom
!!
!! OUTPUT
!!  (the polarization is printed)
!!
!! SIDE EFFECTS
!!
!!
!! TODO
!!  Cleaning, checking for rules.
!!  Should allow for time-reversal symmetry (istwfk)
!!  Should use randac to scan rapidly the wf file
!!
!! NOTES
!! Local Variables:
!!  cmatrix(:,:,:)= overlap matrix of size maxband*maxband
!!  cg_index(:,:,:)= unpacked cg index array for specific band,
!!   k point and polarization.
!!  det(2,2)= intermediate output of Lapack routine zgedi.f
!!  determinant(:,:)= determinant of cmatrix
!!  det_average(2)=  averaged det_string over all strings
!!  det_string(:,:)= determinant product of cmatrices along each string
!!  dk(3)= step taken to the next k mesh point along the kberry direction
!!  dkptnext(3)= step between the next and current k point
!!  dphase= phase angle computed from rel_string(2)
!!  gpard(3)= dimensionalreciprocal lattice vector G along which the
!!          polarization is computed
!!  firstk(mpw)/secondk(mpw)= planewaves indexing array of a first
!!          k point and a second k point
!!  kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
!!          planewave and k point
!!  kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
!!  kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT code
!!                 kpt_flag(ikpt) gives the indices of the k-point related
!!                   to ikpt by time reversal symetrie
!!  kpt_mark(nkpt)= 0, if k point is unmarked; 1, if k point has been marked
!!  maxband/minband= control the minimum and maximum band calculated in the
!!           overlap matrix
!!  nkstr= number of k points per string
!!  npw_k= npwarr(ikpt), number of planewaves in basis at this k point
!!  nstr= number of k point strings
!!  nkpt=number of k points in the whole BZ
!!  phase0=  phase angle computed from det_average(2)
!!  polberry(:)= berry phase of each string (2/nsppol)*(phase0+dphase)/two_pi
!!  polb(isppol) = total berry phase polarization for each spin
!!  polbtot= total berry phase polarization
!!  polion=  ionic polarization for each ion
!!  politot= total ionic polarization
!!  poltot=  total polarization =  polbtot + politot
!!  rel_string(2)= det_string(2)/det_average(2)
!!  shift_g(nkpt)= .true. if the k point should be shifted by a G vector;
!!          .false. if not
!!  tr(2)=variable that changes k to -k
!!                              G to -G
!!                              $c_g$ to $c_g^*$
!!          when time-reversal symetrie is used
!!  xcart(3,natom)= cartesian coordinates of atoms (bohr)
!!  xcart_reindex(:,:,:)= unpack xcart for each atomic species and number
!!           of atoms for each species
!!
!! WARNING
!! This routine is not yet memory optimized
!! It might be also rather time-consuming, since there is a
!! double loop on the number of plane waves.
!! xred is not modified by xredxcart evenf if it is declared inout (TD).
!!
!! PARENTS
!!      elpolariz
!!
!! CHILDREN
!!      cgedi,cgefa,hdr_skip,leave_new,matr3inv,rwwf,wrtout,xredxcart,zgedi
!!      zgefa
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine berryphase(atindx1,bdberry,cg,gprimd,istwfk,kberry,kg,kpt_,&
&  kptopt,kptrlatt,mband,&
&  mgfft,mkmem,mpw,natom,nattyp,nband,nberry,npwarr,nspinor,nsppol,ntypat,&
&  nkpt_,rprimd,typat,ucvol,unkg,wffnow,xred,zion)

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_11util
 use interfaces_12geometry
 use interfaces_13io_mpi
#else
 use defs_interfaces
#endif
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: kptopt,mband,mgfft,mkmem,mpw,natom,nberry,nkpt_,nspinor
 integer,intent(in) :: nsppol,ntypat,unkg
 real(dp),intent(in) :: ucvol
 type(wffile_type),intent(inout) :: wffnow
!arrays
 integer,intent(in) :: atindx1(natom),bdberry(4),istwfk(nkpt_),kberry(3,nberry)
 integer,intent(in) :: kg(3,mpw*mkmem),kptrlatt(3,3),nattyp(ntypat)
 integer,intent(in) :: nband(nkpt_*nsppol),npwarr(nkpt_),typat(natom)
 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gprimd(1:3,1:3)
 real(dp),intent(in) :: kpt_(3,nkpt_),rprimd(3,3),zion(ntypat)
 real(dp),intent(inout) :: xred(3,natom)

!Local variables -------------------------
!scalars
 integer :: band_in,cg_index_iband,cg_index_jband,flag,flag1,formeig,iatom
 integer :: iattyp,iband,iberry,icg,idir,ierr,ii,ikpt,ikpt2,index,index1,info
 integer :: ipw,isp,isppol,istr,istwf_k,itypat,iunmark,jband,jj,jkpt,jkstr
 integer :: jkstr_ori,jpw,jsppol,lkstr,lkstr_ori,lkstr_ori_,maxband,mcg_disk
 integer :: minband,muig,nband_k,ndisk,nkpt,nkstr,npw_k,nstr,option,read_k
 integer :: tim_rwwf
 real(dp) :: det_mod,dphase,ecut,fac,gmod,phase0,pol,polbtot,polion,politot
 real(dp) :: poltot
 character(len=500) :: message
!arrays
 integer :: dg(3),firstk(mpw),kpt_flag(2*nkpt_),kpt_mark(2*nkpt_),secondk(mpw)
 integer,allocatable :: cg_index(:,:,:),ikpt_dk(:),ikstr(:,:),ipvt(:)
 integer,allocatable :: kg_dum(:,:),kg_jl(:,:,:),kg_kpt(:,:,:)
 real(dp) :: cg_temp(mpw,2),det(2,2),det_average(2),diffk(3),dk(3),gpard(3)
 real(dp) :: klattice(3,3),kptrlattr(3,3),polb(nsppol),rel_string(2),tr(2)
 real(dp) :: xcart(3,natom)
 real(dp),allocatable :: cg_disk(:,:,:),cmatrix(:,:,:),det_string(:,:)
 real(dp),allocatable :: det_tmp(:,:),determinant(:,:),eig_dum(:),kpt(:,:)
 real(dp),allocatable :: occ_dum(:),polberry(:),xcart_reindex(:,:,:)
 real(dp),allocatable :: zgwork(:,:)
 logical,allocatable :: shift_g(:)

! ***********************************************************************
!DEBUG
!write(6,*)' berryphase : enter '
!ENDDEBUG

 if(nspinor==2)then
  write(6,*)' berryphase : does not yet work for nspinor=2'
  stop
 end if

 if(maxval(istwfk(:))/=1)then
  write(message, '(a,a,a,a,a,a)' )ch10,&
&  ' berryphase : BUG -',ch10,&
&  '  Sorry, this routine does not work yet with istwfk/=1.',ch10,&
&  '  This should have been tested previously ...'
  call wrtout(6,message,'COLL')
  call leave_new('COLL')
 end if

!change8: set up the whole k point grid in the case where kptopt = 2
 if (kptopt==3) then
  nkpt = nkpt_
  allocate(kpt(3,nkpt))
  kpt(:,:)=kpt_(:,:)
 else if (kptopt==2) then
  nkpt = nkpt_*2
  allocate(kpt(3,nkpt))
  do ikpt = 1,nkpt/2
   kpt_flag(ikpt) = 0
   kpt(:,ikpt)=kpt_(:,ikpt)
  end do
  index = 0
  do ikpt = (nkpt/2+1),nkpt
   flag1 = 0
   do jkpt = 1, nkpt/2
    if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
&       (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
&      .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
&          (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
&      .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
&          (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
     flag1 = 1
     index = index + 1
     exit
    end if
   end do
   if (flag1==0) then
    kpt_flag(ikpt-index)=ikpt-nkpt/2
    kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
   end if
  end do
  nkpt = nkpt - index
 end if

!change8

 allocate(shift_g(nkpt))
 allocate(kg_dum(3,0))

!Compute primitive vectors of the k point lattice
!Copy to real(dp)
 kptrlattr(:,:)=kptrlatt(:,:)
!Go to reciprocal space (in reduced coordinates)
 call matr3inv(kptrlattr,klattice)

 do iberry=1,nberry

! Calculate dimensional recip lattice vector along which P is calculated
! dk =  step to the nearest k point along that direction
! in reduced coordinates
  dk(:)=kberry(1,iberry)*klattice(:,1)+&
&       kberry(2,iberry)*klattice(:,2)+&
&       kberry(3,iberry)*klattice(:,3)
  gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
  gmod=sqrt(dot_product(gpard,gpard))

! *****************************************************************************
! Select the k grid  points along the kberry direction
! dk =  step to the nearest k point along that direction

! For each k point, find k_prim such that k_prim= k + dk mod(G)
! where G is a vector of the reciprocal lattice
  allocate(ikpt_dk(nkpt))
  shift_g(:)= .false.
  do ikpt=1,nkpt
   do ikpt2=1,nkpt
    diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
    if(sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
     ikpt_dk(ikpt)=ikpt2
     if(sum(diffk(:))>=3*tol8)shift_g(ikpt2) = .true.
     exit
    end if
   end do
  end do

!DEBUG
!do ikpt = 1,nkpt
!write(100,*)'ikpt_dk = ',ikpt_dk(ikpt)
!if (shift_g(ikpt))then
! write(100,*)'true'
!else
! write(100,*)'false'
!end if
!write(100,*)''
!end do
!ENDDEBUG

! Find the string length, starting from k point 1
! (all strings must have the same number of points)
  nkstr=1
  ikpt2=1
  do ikpt=1,nkpt
   ikpt2=ikpt_dk(ikpt2)
   if(ikpt2==1)exit
   nkstr=nkstr+1
  end do

! Check that the string length is a divisor of nkpt
  if(mod(nkpt,nkstr)/=0)then
   write(message,'(a,a,a,a,i5,a,i7)')ch10,&
&   ' berryphase: BUG -',ch10,&
&   '  The string length=',nkstr,', is not a divisor of nkpt=',nkpt
   call wrtout(6,message,'COLL')
  end if
  nstr=nkpt/nkstr

  write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
&  ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,&
&  dk(:),' (in reduced coordinates)',ch10,&
&  gpard(1:3),' (in cartesian coordinates - atomic units)'
  call wrtout(ab_out,message,'COLL')
  call wrtout(6,message,'COLL')

  write(message,'(a,i5,a,a,i5)')&
&  ' Number of strings: ',nstr,ch10,&
&  ' Number of k points in string:', nkstr
  call wrtout(6,message,'COLL')

  if(nsppol==1)then
   write(message, '(a,i5,a,i5)')&
&   ' From band number',bdberry(1),'  to band number',bdberry(2)
  else
   write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
&   ' From band number',bdberry(1),'  to band number',bdberry(2),' for spin up,',&
&   ch10,&
&   ' from band number',bdberry(3),'  to band number',bdberry(4),' for spin down.'
  end if
  call wrtout(ab_out,message,'COLL')
  call wrtout(6,message,'COLL')

!DEBUG
! write(6,*)' berryphase : find nkpt,nkstr,nstr=',nkpt,nkstr,nstr
! stop
!ENDDEBUG

! Build the different strings
  allocate(ikstr(nkstr,nstr))

  iunmark=1
  kpt_mark(:)=0
  do istr = 1, nstr
   do while(kpt_mark(iunmark)/=0)
    iunmark = iunmark + 1
   end do
   ikstr(1, istr) = iunmark
   kpt_mark(iunmark)=1
   do jkstr = 2, nkstr
    ikstr(jkstr,istr)=ikpt_dk(ikstr(jkstr-1,istr))
    kpt_mark(ikstr(jkstr,istr))=1
   end do
  end do ! istr

!DEBUG
!do istr = 1,nstr
!   do jkstr = 1,nkstr
!   if (shift_g(ikstr(jkstr,istr))) then
!      write(99,*) ikstr(jkstr,istr),'true'
!   else
!      write(99,*) ikstr(jkstr,istr),'false'
!   end if
!   end do
!end do
!ENDDEBUG

  deallocate(ikpt_dk)
! DEBUG!
! write(100,*) 'list all the k points strings:'
! do istr=1,nstr
!  write(100,*) (ikstr(jkstr,istr),jkstr=1,nkstr)
! end do
! ENDDEBUG!

! *****************************************************************************
! Find the location of each wavefunction
  allocate(cg_index(mband,nkpt,nsppol))

  icg = 0
  do isppol=1,nsppol
   do ikpt=1,nkpt_
    nband_k=nband(ikpt+(isppol-1)*nkpt_)
    npw_k=npwarr(ikpt)
    do iband=1,nband_k
     cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
    end do
    icg=icg+npw_k*nspinor*nband(ikpt)
   end do
  end do

!change5
 if (mkmem/=0) then
! Find the planewave vectors and their indexes for each k point
  allocate(kg_kpt(3,mpw*nspinor,nkpt_))
  kg_kpt(:,:,:) = 0
  index1 = 0
  do ikpt=1,nkpt_
   npw_k=npwarr(ikpt)
   do ipw=1,npw_k*nspinor
    kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
   end do
   index1=index1+npw_k*nspinor
  end do
 end if            !change5
!*****************************************************************************
  allocate(det_string(2, nstr), det_tmp(2, nstr) )
  allocate(polberry(nstr))

!change1
  if (mkmem==0) then
   call hdr_skip(wffnow,ierr)

!  Should Define offsets, in case of MPI I/O : also next calls !!!
!  For the time being, does not work in parallel
!  formeig=0
!  call xdefineOff(formeig,wffnow,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt)

   mcg_disk=mpw*nspinor*mband
   allocate(cg_disk(2,mcg_disk,2))
   rewind unkg
   allocate(kg_jl(3,mpw,2))
  end if

! Initialize berry phase polarization for each spin and the total one
  polb(1:nsppol) = 0.0_dp
  polbtot=0.0_dp

! Loop over spins
  do isppol=1,nsppol

   minband=bdberry(2*isppol-1)
   maxband=bdberry(2*isppol)

   if(minband<1)then
    write(message,'(a,a,a,a,i5,a)')ch10,&
&    ' berryphase : BUG - ',ch10,&
&    '  The band limit minband=',minband,', is lower than 0.'
    call wrtout(6,message,'COLL')
    call leave_new('COLL')
   end if
   if(maxband<1)then
    write(message,'(a,a,a,a,i5,a)')ch10,&
&    ' berryphase : BUG - ',ch10,&
&    '  The band limit maxband=',maxband,', is lower than 0.'
    call wrtout(6,message,'COLL')
    call leave_new('COLL')
   end if
   if(maxband<minband)then
    write(message,'(4a,i5,a,i5)')ch10,&
&    ' berryphase : BUG - ',ch10,&
&    '  maxband=',maxband,', is lower than minband=',minband
    call wrtout(6,message,'COLL')
    call leave_new('COLL')
   end if

!  Initialize det_string and det_average
   det_string(1, 1:nstr) = 1.0_dp; det_string(2, 1:nstr) = 0.0_dp
   det_average(1:2)=0.0_dp; det_average(2)=0.0_dp

!  Loop over strings
   do istr = 1, nstr

!change7
    read_k = 0

!   DEBUG!
!   write(100,'(a,i4)') 'This is in string', istr
!   ENDDEBUG!

!   Loop over k points per string
    allocate(determinant(2, nkstr))

    do jkstr=1,nkstr

     allocate(cmatrix(2,maxband,maxband))
     if(jkstr < nkstr) then
      lkstr=jkstr+1
     else
      lkstr= jkstr+1-nkstr
     end if
     jkstr_ori=ikstr(jkstr,istr)
     lkstr_ori=ikstr(lkstr,istr)

!change9
     lkstr_ori_=lkstr_ori
     tr(1) = 1.0_dp
     tr(2) = 1.0_dp
     if (kptopt==2) then
      if (read_k == 0) then
       if (kpt_flag(jkstr_ori)/=0) then
        tr(1) = -1.0_dp
        jkstr_ori = kpt_flag(jkstr_ori)
       end if
       if (kpt_flag(lkstr_ori)/=0) then
        tr(2) = -1.0_dp
        lkstr_ori = kpt_flag(lkstr_ori)
       end if
      else           !read_k
       if (kpt_flag(jkstr_ori)/=0) then
        tr(-1*read_k+3) = -1.0_dp
        jkstr_ori = kpt_flag(jkstr_ori)
       end if
       if (kpt_flag(lkstr_ori)/=0) then
        tr(read_k) = -1.0_dp
        lkstr_ori = kpt_flag(lkstr_ori)
       end if
      end if       !read_k
     end if           !kptopt
!change9

     nband_k=nband(jkstr_ori+(isppol-1)*nkpt_)
     if(nband_k<maxband)then
      write(message,'(4a,i5,a,i5)')ch10,&
&      ' berryphase : BUG - ',ch10,&
&      '  maxband=',maxband,', is larger than nband(j,isppol)=',nband_k
      call wrtout(6,message,'COLL')
      call leave_new('COLL')
     end if
     nband_k=nband(lkstr_ori+(isppol-1)*nkpt_)
     if(nband_k<maxband)then
      write(message,'(4a,i5,a,i5)')ch10,&
&      ' berryphase : BUG - ',ch10,&
&      '  maxband=',maxband,', is larger than nband(l,isppol)=',nband_k
      call wrtout(6,message,'COLL')
      call leave_new('COLL')
     end if

!change2
     if (mkmem==0) then

!if read_k = 0,read first k point of string
!******************************************
      if (read_k==0) then
       read_k = 1
       npw_k = npwarr(jkstr_ori)
       rewind unkg
       index = 1
       do while(index < jkstr_ori)
        read(unkg)
        read(unkg)
        read(unkg)
        index = index + 1
       end do
       read(unkg) ndisk
       read(unkg)
       read(unkg) ((kg_jl(ii,muig,read_k),ii=1,3),muig=1,npw_k)
       tim_rwwf = 0
       allocate(eig_dum(mband),occ_dum(mband))

       call hdr_skip(wffnow,ierr)
       do isp=1,nsppol
        do ikpt=1,nkpt_
         if(isp==isppol .and. ikpt==jkstr_ori)then
          option=-2  ! will read cg only
         else
          option=-1  ! will skip the records
         end if
!        This should be changed : one should not call a segment of the cg_disk array....
         call rwwf(cg_disk(:,:,read_k),eig_dum,0,0,0,jkstr_ori,isppol,kg_dum,mband,mcg_disk,nband_k, &
&         nband_k,npw_k,nspinor,occ_dum,option,0,tim_rwwf,wffnow)
         if(option==-2)exit
        end do
        if(option==-2)exit
       end do
       deallocate(eig_dum,occ_dum)
       read_k = 2
      end if           !read_k

!     read the next k point
!     *********************
      if (read_k /= 0) then
       npw_k = npwarr(lkstr_ori)
       rewind unkg
       index = 1
       do while(index < lkstr_ori)
        read(unkg)
        read(unkg)
        read(unkg)
        index = index + 1
       end do
       read(unkg) ndisk
       read(unkg)
       read(unkg) ((kg_jl(ii,muig,read_k),ii=1,3),muig=1,npw_k)
       tim_rwwf = 0
       allocate(eig_dum(mband),occ_dum(mband))

       call hdr_skip(wffnow,ierr)
       do isp=1,nsppol
        do ikpt=1,nkpt_
         if(isp==isppol .and. ikpt==lkstr_ori)then
          option=-2  ! will read cg only
         else
          option=-1  ! will skip the records
         end if
!        This should be changed : one should not call a segment of the cg_disk array....
         call rwwf(cg_disk(:,:,read_k),eig_dum,0,0,0,lkstr_ori,isppol,kg_dum,mband,mcg_disk,nband_k, &
&         nband_k,npw_k,nspinor,occ_dum,option,0,tim_rwwf,wffnow)
         if(option==-2)exit
        end do
        if(option==-2)exit
       end do
       deallocate(eig_dum,occ_dum)
      end if           !read_k
     end if         !change2 (mkmem)

     if (jkstr==1) read_k = 2
!    Compute the overlap matrix <u_k|u_k+b>
     cmatrix(1:2,1:maxband,1:maxband)=0.0_dp
     jj = read_k
     ii = -1*read_k+3
     if(shift_g(lkstr_ori_) .eqv. .false.) then
!     Change3
      if (mkmem == 0) then
       do ipw=1,npwarr(jkstr_ori)
        do jpw=1,npwarr(lkstr_ori)
!        Check if  Fourier components of jkstr and jkstr+1 matches

         if((tr(ii)*kg_jl(1,ipw,ii)==tr(jj)*kg_jl(1,jpw,jj))&
&         .and.(tr(ii)*kg_jl(2,ipw,ii) == tr(jj)*kg_jl(2,jpw,jj))&
&         .and.(tr(ii)*kg_jl(3,ipw,ii) == tr(jj)*kg_jl(3,jpw,jj)))&
&          then
          do iband=minband,maxband
           cg_index_iband= (iband-1)*npwarr(jkstr_ori)
           do jband=minband,maxband
            cg_index_jband= (jband-1)*npwarr(lkstr_ori)

            cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
&            cg_disk(1,ipw+cg_index_iband,ii)*cg_disk(1,jpw+cg_index_jband,jj)+&
&            tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*tr(jj)*cg_disk(2,jpw+cg_index_jband,jj)
            cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
&            cg_disk(1,ipw+cg_index_iband,ii)*tr(jj)*cg_disk(2,jpw+cg_index_jband,jj)-&
&            tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*cg_disk(1,jpw+cg_index_jband,jj)

           end do !jband
          end do !iband
          exit  !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
         end if

        end do ! jpw
       end do ! ipw
      else                 !change3
       do ipw=1,npwarr(jkstr_ori)
        do jpw=1,npwarr(lkstr_ori)

!        Check if  Fourier components of jkstr and jkstr+1 matches

         if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori))&
&         .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori))&
&         .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)))&
&          then

          do iband=minband,maxband
           cg_index_iband=cg_index(iband,jkstr_ori,isppol)
           do jband=minband,maxband
            cg_index_jband=cg_index(jband,lkstr_ori,isppol)

            cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
&            cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+&
&            tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)
            cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
&            cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-&
&            tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)

           end do !jband
          end do !iband
          exit  !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
         end if

        end do ! jpw
       end do ! ipw

      end if             !change3

!     But there is a special pair of k points which involves the shift of a
!     G vector

     else

      dg(:) = -1*nint(tr(jj)*kpt(:,lkstr_ori)-tr(ii)*kpt(:,jkstr_ori)-dk(:))

!DEBUG
!write(100,*)dg
!write(100,*)kberry(:,iberry)
!write(100,*)''
!ENDDEBUG

!change4
      if (mkmem==0) then

       do ipw=1,npwarr(jkstr_ori)
        do jpw=1,npwarr(lkstr_ori)

!        Check if  Fourier components of jkstr and jkstr+1
!        matches by comparing the G vectors

         if((tr(ii)*kg_jl(1,ipw,ii)==tr(jj)*kg_jl(1,jpw,jj)-dg(1))&
&         .and.(tr(ii)*kg_jl(2,ipw,ii) == tr(jj)*kg_jl(2,jpw,jj)-dg(2))&
&         .and.(tr(ii)*kg_jl(3,ipw,ii) == tr(jj)*kg_jl(3,jpw,jj)-dg(3)))&
&          then

          do iband=minband,maxband
           cg_index_iband=(iband-1)*npwarr(jkstr_ori)

           do jband=minband,maxband
            cg_index_jband=(jband-1)*npwarr(lkstr_ori)

            cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
&            cg_disk(1,ipw+cg_index_iband,ii)*cg_disk(1,jpw+cg_index_jband,jj)+&
&            tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*tr(jj)*cg_disk(2,jpw+cg_index_jband,jj)
            cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
&            cg_disk(1,ipw+cg_index_iband,ii)*tr(jj)*cg_disk(2,jpw+cg_index_jband,jj)-&
&            tr(ii)*cg_disk(2,ipw+cg_index_iband,ii)*cg_disk(1,jpw+cg_index_jband,jj)

           end do ! jband
          end do ! iband
         exit  !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
         end if
        end do ! jpw
       end do ! ipw
      else              !change4
       do ipw=1,npwarr(jkstr_ori)
        do jpw=1,npwarr(lkstr_ori)

!        Check if  Fourier components of jkstr and jkstr+1
!        matches by comparing the G vectors

         if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori)-dg(1))&
&         .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori)-dg(2))&
&         .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)-dg(3)))&
&          then

          do iband=minband,maxband
           cg_index_iband=cg_index(iband,jkstr_ori,isppol)

           do jband=minband,maxband
            cg_index_jband=cg_index(jband,lkstr_ori,isppol)

            cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
&            cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+&
&            tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)
            cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
&            cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-&
&            tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)

           end do ! jband
          end do ! iband
         exit  !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
         end if
        end do ! jpw
       end do ! ipw
      end if                       !change4
     end if

!    Compute the determinant of cmatrix(1:2,minband:maxband, minband:maxband)

     band_in = maxband - minband + 1

     allocate(ipvt(maxband))
     allocate(zgwork(2,1:maxband))

!    Last argument of zgedi means calculate determinant only.
#if defined T3E
     call cgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
     call cgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,10)
#else
     call zgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
     call zgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,10)
#endif

     deallocate(zgwork,ipvt)

     fac=exp(log(10._dp)*det(1,2))
     determinant(1, jkstr) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - &
&      det(2,1)*sin(log(10._dp)*det(2,2)))
     determinant(2, jkstr) = fac*(det(1,1)*sin(log(10._dp)*det(2,2)) + &
&      det(2,1)*cos(log(10._dp)*det(2,2)))
!    DEBUG!
!     write(100,*) 'det',jkstr,lkstr,'=', determinant(1:2,jkstr)
!    ENDDEBUG!

     det_tmp(1,istr) = det_string(1,istr)*determinant(1,jkstr) - &
&      det_string(2,istr)*determinant(2,jkstr)
     det_tmp(2,istr) = det_string(1,istr)*determinant(2,jkstr) + &
&      det_string(2,istr)*determinant(1,jkstr)
     det_string(1:2,istr) = det_tmp(1:2,istr)

     deallocate(cmatrix)

!    Close loop over k points along string
     read_k = -1*read_k + 3             ! read_k=2 <-> read_k=1
    end do

!   DEBUG!
!    write(100,*) 'det_string =',  det_string(1:2,istr)
!    write(100,*)
!   ENDDEBUG!

    det_average(1) = det_average(1) + det_string(1,istr)/nstr
    det_average(2) = det_average(2) + det_string(2,istr)/nstr

    deallocate(determinant)

!  Close loop over strings
   end do


!*****************************************************************************
!  Calculate the electronic contribution to the polarization

   write(message,'(a,a)')ch10,&
&   ' Compute the electronic contribution to polarization'
   call wrtout(6,message,'COLL')

!  First berry phase that corresponds to det_average
   phase0 = atan2(det_average(2),det_average(1))
   det_mod = det_average(1)**2+det_average(2)**2

!  Then berry phase that corresponds to each string relative to the average
   do istr = 1, nstr
    rel_string(1) = (det_string(1,istr)*det_average(1) + &
       det_string(2,istr)*det_average(2))/det_mod
    rel_string(2) = (det_string(2,istr)*det_average(1) - &
       det_string(1,istr)*det_average(2))/det_mod
    dphase = atan2(rel_string(2),rel_string(1))
    polberry(istr) = (2.0_dp/nsppol)*(phase0+dphase)/two_pi
    polb(isppol) = polb(isppol) + polberry(istr)/nstr
   end do

!  Output berry phase polarization
   write(message,'(a,10x,a,10x,a)')ch10,&
&   'istr','polberry(istr)'
   call wrtout(6,message,'COLL')
   do istr=1,nstr
    write(message,'(10x,i4,7x,e15.9)')istr,polberry(istr)
    call wrtout(6,message,'COLL')
   end do

   write(message,'(9x,a,7x,e15.9,1x,a,i4,a,a)')&
&   'total',polb(isppol),'(isppol=',isppol,')',ch10
   call wrtout(6,message,'COLL')

   polbtot=polbtot+polb(isppol)

  end do ! isppol

  deallocate(polberry)
  deallocate(det_tmp, det_string)
  deallocate(ikstr, cg_index)
!change6
  if (mkmem /=0) deallocate(kg_kpt)
  if (mkmem == 0) deallocate(cg_disk,kg_jl)
! *****************************************************************************
! Reindex xcart according to atom and type
  call xredxcart(natom,1,rprimd,xcart,xred)
  allocate(xcart_reindex(3,natom,ntypat))
  index=1
  do itypat=1,ntypat
   do iattyp=1,nattyp(itypat)
    iatom=atindx1(index)
    xcart_reindex(1:3,iattyp,itypat) = xcart(1:3,iatom)
    index = index+1
   end do
  end do

! Compute the ionic contribution to the polarization
  politot = 0.0_dp
  write(message,'(a)')' Compute the ionic contributions'
  call wrtout(6,message,'COLL')

  write(message,'(a,2x,a,2x,a,15x,a)')ch10,&
&  'itypat', 'iattyp', 'polion'
  call wrtout(6,message,'COLL')

  do itypat=1,ntypat
   do iattyp=1,nattyp(itypat)
    polion=zion(itypat)*nkstr*&
&    dot_product(xcart_reindex(1:3,iattyp,itypat),gpard(1:3))
!   Fold into interval (-1,1)
    polion=polion-2._dp*nint(polion/2.0_dp)
    politot=politot+polion
    write(message,'(2x,i2,5x,i2,10x,e15.9)') itypat,iattyp,polion
    call wrtout(6,message,'COLL')
   end do
  end do

! Fold into interval (-1,1) again
  politot=politot-2.0_dp*nint(politot/2.0_dp)

  write(message,'(9x,a,7x,es19.9)') 'total',politot
  call wrtout(6,message,'COLL')

  deallocate(xcart_reindex)

! Compute the total polarizations

  poltot=politot+polbtot

  write(message,'(a,a)')ch10,&
&  ' Summary of the results'
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

  write(message,'(a,es19.9)')&
&  ' Electronic Berry phase ' ,polbtot
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

  write(message,'(a,es19.9)') &
&  '            Ionic phase ', politot
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

  write(message,'(a,es19.9)') &
&  '            Total phase ', poltot
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

  poltot=poltot-2.0_dp*nint(poltot/2._dp)
  write(message,'(a,es19.9)') &
&  '    Remapping in [-1,1] ', poltot
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

! Transform the phase into a polarization
  fac = 1._dp/(gmod*nkstr)
  fac = fac/ucvol
  pol = fac*poltot

  write(message,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,&
&  '           Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,&
&  '           Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,&
&        ' C/m^2',ch10
  call wrtout(6,message,'COLL')
  call wrtout(ab_out,message,'COLL')

 end do ! iberry

 deallocate(shift_g)
 deallocate(kpt,kg_dum)
end subroutine berryphase
!!***

Generated by  Doxygen 1.6.0   Back to index