./72_response/corrmetalwf1.F90

ABINIT/corrmetalwf1

ABINIT/corrmetalwf1 [ Functions ]

NAME

 corrmetalwf1

FUNCTION

 Response function calculation only:
 Correct 1st-order wave-function, taking into account "metallic" occupations.
 1st-order WF orthogonal to C_n,k+q, restore the "active space" content of the first-order WF.

 Copyright (C) 2009-2018 ABINIT group (MT)
 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,mcgq)=planewave coefficients of wavefunctions at k+q
  cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
  cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
  cwaveprj(natom,nspinor)= 1st-order wave-function before correction
                           projected on NL projectors (PAW)
  eig1(2*nband**2)=first-order eigenvalues (hartree)
  fermie1=derivative of fermi energy wrt (strain) perturbation
  ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
                      (C1 before correction)
  iband=index of current band
  ibgq=shift to be applied on the location of data in the array cprjq
  icgq=shift to be applied on the location of data in the array cgq
  istwf_k=option parameter that describes the storage of wfs
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nband=number of bands
  npw1=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  occ(nband)=occupation number for each band for each k.
  rocceig(nband,nband)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
    if this ratio has been attributed to the band n (second argument), zero otherwise
  timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
  usepaw=flag for PAW

OUTPUT

  cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
  cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
                            projected on NL projectors (PAW)
  edocc(nband)=correction to 2nd-order total energy coming from changes of occupations
  wf_corrected=flag put to 1 if input cwave1 is effectively different from output cwavef

NOTES

  Was part of dfpt_vtowfk before.

PARENTS

      dfpt_vtowfk

CHILDREN

      cg_zcopy,dotprod_g,pawcprj_copy,pawcprj_zaxpby,timab

SOURCE

 64 #if defined HAVE_CONFIG_H
 65 #include "config.h"
 66 #endif
 67 
 68 #include "abi_common.h"
 69 
 70 subroutine corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,edocc,eig1,fermie1,ghc,iband, &
 71 &          ibgq,icgq,istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,nspinor,occ,rocceig,timcount,&
 72 &          usepaw,wf_corrected)
 73 
 74  use defs_basis
 75  use defs_abitypes
 76  use m_errors
 77  use m_profiling_abi
 78  use m_cgtools
 79 
 80  use m_pawcprj, only : pawcprj_type, pawcprj_copy, pawcprj_zaxpby
 81 
 82 !This section has been created automatically by the script Abilint (TD).
 83 !Do not modify the following lines by hand.
 84 #undef ABI_FUNC
 85 #define ABI_FUNC 'corrmetalwf1'
 86  use interfaces_18_timing
 87 !End of the abilint section
 88 
 89  implicit none
 90 
 91 !Arguments ------------------------------------
 92 !scalars
 93  integer,intent(in) :: iband,ibgq,icgq,istwf_k,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw
 94  integer,intent(out) :: wf_corrected
 95  real(dp),intent(in) :: fermie1
 96  type(MPI_type),intent(in) :: mpi_enreg
 97 !arrays
 98  real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor)
 99  real(dp),intent(in) :: eig1(2*nband**2),ghc(2,npw1*nspinor),occ(nband),rocceig(nband,nband)
100  real(dp),intent(out) :: cwave1(2,npw1*nspinor),edocc(nband)
101  type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw)
102  type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i
103 
104 !Local variables-------------------------------
105 !scalars
106  integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii
107  real(dp) :: facti,factr,invocc
108 !arrays
109  real(dp) :: tsec(2)
110  real(dp),allocatable :: cwcorr(:,:)
111 
112 ! *********************************************************************
113 
114  DBG_ENTER("COLL")
115 
116  call timab(214+timcount,1,tsec)
117 
118 !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf).
119 !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 .
120 
121 !First copy input WF into output WF
122  wf_corrected=0
123  call cg_zcopy(npw1*nspinor,cwavef,cwave1)
124 
125  if (usepaw==1) then
126    call pawcprj_copy(cwaveprj,cwaveprj1)
127  end if
128 
129 !Correct WF only for occupied states
130  if (abs(occ(iband)) > tol8) then
131    invocc=one/occ(iband)
132 
133    edocc(iband)=zero
134 
135 !  Loop over WF at k+q subspace
136    do ibandkq=1,nband
137 
138 !    Select bands with variable occupation
139      if (abs(rocceig(ibandkq,iband))>tol8) then
140 
141        wf_corrected=1
142 
143        index_eig1=2*ibandkq-1+(iband-1)*2*nband
144        index_cgq=npw1*nspinor*(ibandkq-1)+icgq
145 
146        if(ibandkq==iband) then
147          factr=rocceig(ibandkq,iband)*invocc*(eig1(index_eig1)-fermie1)
148        else
149          factr=rocceig(ibandkq,iband)*invocc*eig1(index_eig1)
150        end if
151        facti= rocceig(ibandkq,iband)*invocc*eig1(index_eig1+1)
152 
153 !      Apply correction to 1st-order WF
154 !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwave1,facti,factr,index_cgq,npw1,nspinor)
155        do ii=1,npw1*nspinor
156          cwave1(1,ii)=cwave1(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
157          cwave1(2,ii)=cwave1(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq))
158        end do
159 
160 !      In the PAW case, also apply correction to projected WF
161        if (usepaw==1) then
162          index_cprjq=nspinor*(ibandkq-1)+ibgq
163          call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1)
164        end if
165 
166 !      The factor of two is needed because we compute the 2DTE, and not E(2)
167        edocc(iband)=edocc(iband)-two*(factr*eig1(index_eig1)+facti*eig1(index_eig1+1))
168 
169      end if ! Variable occupations
170    end do ! Loop over k+q subspace
171  end if ! occupied states
172 
173 !In the PAW case, compute <Psi^(1)_ortho|H-Eig0_k.S|Psi^(1)_parallel> contribution to 2DTE
174  if (usepaw==1.and.wf_corrected==1) then
175    ABI_ALLOCATE(cwcorr,(2,npw1*nspinor))
176 !$OMP WORKSHARE
177    cwcorr(:,:)=cwave1(:,:)-cwavef(:,:)
178 !$OMP END WORKSHARE
179    call dotprod_g(factr,facti,istwf_k,npw1*nspinor,1,cwcorr,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
180    edocc(iband)=edocc(iband)+four*factr
181    ABI_DEALLOCATE(cwcorr)
182  end if
183 
184  call timab(214+timcount,2,tsec)
185 
186  DBG_EXIT("COLL")
187 
188 end subroutine corrmetalwf1

TABLE OF CONTENTS

ABINIT/corrmetalwf1 [ Functions ]