TABLE OF CONTENTS


ABINIT/qmatrix [ Functions ]

[ Top ] [ Functions ]

NAME

 qmatrix

FUNCTION

 calculation of the inverse of the overlap matrix

COPYRIGHT

 Copyright (C) 2004-2018 ABINIT group (XW).
 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

 cg(2,mpw*nspinor*mband*mkmem*nsppol) = planewave coefficients of wavefunctions
 RF wavefunctions at k,q.
 dtefield = variables related to response Berry-phase calculation
 ikpt = the index of the current k point
 mband =  maximum number of bands
 mkmem = maximum number of k-points in core memory
 mpw = maximum number of plane waves
 mpw1 = maximum number of plane waves for response wavefunctions
 nkpt = number of k points
 npwarr(nkpt) = number of planewaves in basis and boundary at this k point
 npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
 nspinor = 1 for scalar wfs, 2 for spinor wfs
 nsppol = 1 for unpolarized, 2 for spin-polarized
 pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices

OUTPUT

 qmat(2,dtefield%nband_occ,dtefield%nband_occ,nkpt,2,3) =
 inverse of the overlap matrix

PARENTS

      dfpt_scfcv

CHILDREN

      dzgedi,dzgefa,overlap_g

SOURCE

 44 #if defined HAVE_CONFIG_H
 45 #include "config.h"
 46 #endif
 47 
 48 #include "abi_common.h"
 49 
 50 subroutine qmatrix(cg,dtefield,qmat,mpw,mpw1,mkmem,mband,npwarr,nkpt,nspinor,nsppol,pwindall)
 51 
 52  use defs_basis
 53  use m_profiling_abi
 54  use m_efield
 55 
 56  use m_cgtools,   only : overlap_g
 57  use m_abilasi,   only : dzgedi, dzgefa
 58 
 59 !This section has been created automatically by the script Abilint (TD).
 60 !Do not modify the following lines by hand.
 61 #undef ABI_FUNC
 62 #define ABI_FUNC 'qmatrix'
 63 !End of the abilint section
 64 
 65  implicit none
 66 
 67 !Arguments ----------------------------------------
 68 !scalars
 69  integer,intent(in) :: mband,mkmem,mpw,mpw1,nkpt,nspinor,nsppol
 70  type(efield_type),intent(in) :: dtefield
 71 !arrays
 72  integer,intent(in) :: npwarr(nkpt),pwindall(max(mpw,mpw1)*mkmem,8,3)
 73  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 74  real(dp),intent(out) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
 75 
 76 !Local variables -------------------------
 77 !scalars
 78  integer :: iband,icg,icg1,idir,ifor,ikpt,ikpt2,info,jband,job
 79  integer :: npw_k1,npw_k2,pwmax,pwmin
 80  integer :: isppol
 81  real(dp) :: doti,dotr
 82 !arrays
 83  integer,allocatable :: ipvt(:),pwind_k(:)
 84  real(dp) :: det(2,2)
 85  real(dp),allocatable :: sinv(:,:,:),smat_k(:,:,:),vect1(:,:),vect2(:,:)
 86  real(dp),allocatable :: zgwork(:,:)
 87 
 88 ! *************************************************************************
 89 
 90  ABI_ALLOCATE(ipvt,(dtefield%mband_occ))
 91  ABI_ALLOCATE(sinv,(2,dtefield%mband_occ,dtefield%mband_occ))
 92  ABI_ALLOCATE(zgwork,(2,dtefield%mband_occ))
 93  ABI_ALLOCATE(vect1,(2,0:mpw))
 94  ABI_ALLOCATE(vect2,(2,0:mpw))
 95  ABI_ALLOCATE(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ))
 96  ABI_ALLOCATE(pwind_k,(max(mpw,mpw1)))
 97  vect1(:,0) = zero ; vect2(:,0) = zero
 98 
 99  job = 11
100 
101 !**************
102 !loop over k points
103  do isppol = 1, nsppol
104    do ikpt = 1, nkpt
105      npw_k1 = npwarr(ikpt)
106      icg  = dtefield%cgindex(ikpt,1)
107 
108      do idir = 1, 3
109 
110        do ifor = 1, 2
111 
112          ikpt2 = dtefield%ikpt_dk(ikpt,ifor,idir)
113          npw_k2 = npwarr(ikpt2)
114          icg1 = dtefield%cgindex(ikpt2,1)
115          pwind_k(1:npw_k1) = pwindall((ikpt-1)*max(mpw,mpw1)+1:(ikpt-1)*max(mpw,mpw1)+npw_k1,ifor,idir)
116 
117          do jband = 1, dtefield%nband_occ(isppol)
118            vect2(:,1:npw_k2) = &
119 &           cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
120            if (npw_k2 < mpw) vect2(:,npw_k2+1:mpw) = zero
121 
122            do iband = 1, dtefield%nband_occ(isppol)
123 
124              pwmin = (iband-1)*npw_k1*nspinor
125              pwmax = pwmin + npw_k1*nspinor
126              vect1(:,1:npw_k1) = &
127 &             cg(:,icg + 1 + pwmin:icg + pwmax)
128              if (npw_k1 < mpw) vect1(:,npw_k1+1:mpw) = zero
129              call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,&
130 &             vect1,vect2)
131              smat_k(1,iband,jband) = dotr
132              smat_k(2,iband,jband) = doti
133 
134            end do    ! iband
135 
136          end do    !jband
137 
138          sinv(:,:,:) = smat_k(:,:,:)
139 
140          call dzgefa(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,info)
141          call dzgedi(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,det,zgwork,job)
142 
143          qmat(:,:,:,ikpt,ifor,idir) = sinv(:,:,:)
144 
145 
146        end do
147 
148      end do
149 
150    end do  !end loop over k
151  end do
152 
153  ABI_DEALLOCATE(ipvt)
154  ABI_DEALLOCATE(sinv)
155  ABI_DEALLOCATE(zgwork)
156  ABI_DEALLOCATE(vect1)
157  ABI_DEALLOCATE(vect2)
158  ABI_DEALLOCATE(smat_k)
159  ABI_DEALLOCATE(pwind_k)
160 
161 end subroutine qmatrix