TABLE OF CONTENTS


ABINIT/eig2stern [ Functions ]

[ Top ] [ Functions ]

NAME

 eig2stern

FUNCTION

 This routine calculates the second-order eigenvalues.
 The output eig2nkq is this quantity for the input k points.

INPUTS

  bdeigrf = number of bands for which to calculate the second-order eigenvalues.
  clflg(3,mpert)= array on calculated perturbations for eig2rf.
  dim_eig2nkq = 1 if eig2nkq is to be computed.
  cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = first-order wf in G
            space for each perturbation. The wavefunction is orthogonal to the
            active space.
  gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert) = matrix containing the
            vector:  <G|H(0)|psi(1)>, for each perturbation.
  gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol,3,mpert)) = matrix containing the
            vector:  <G|H(1)|n,k>, for each perturbation. The wavefunction is
            orthogonal to the active space.
  eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the
            electronic eigenvalues (optional).
  eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points:
            <k,n'|H(0)|k,n'> (hartree).
  eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
            <k+Q,n'|H(0)|k+Q,n'> (hartree).
  eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order:
            <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
  eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
            the electronic eigenvalues.
  elph2_imagden = imaginary part of the denominator of the sum-over-state expression
            for the electronic eigenenergy shift due to second-order electron-phonon
            interation.
  ieig2rf = integer for calculation type.
  indsym(4,nsym,natom) = indirect indexing array for atom labels
            (not used yet, but will be used with symmetries).
  istwfk_pert(nkpt_rbz,3,mpert) = integer for choice of storage of wavefunction at
            each k point for each perturbation.
  mband = maximum number of bands.
  mk1mem = maximum number of k points which can fit in memory (RF data);
            0 if use disk.
  mpert = maximum number of perturbations.
  natom = number of atoms in the unit cell.
  npert = number of phonon perturbations, without taking into account directions:
            natom.
  nsym = number of symmetries (not used yet).
  mpi_enreg = information about MPI parallelization.
  mpw1 = maximum number of planewaves used to represent first-order wavefunctions.
  nkpt_rbz = number of k-points for each perturbation.
  npwar1(nkpt_rbz,mpert) = number of planewaves at k-point for first-order.
  nspinor = number of spinorial components of the wavefunctions.
  nsppol = 1 for unpolarized, 2 for spin-polarized.
  occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
  smdelta = integer controling the calculation of electron lifetimes.
  symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From littlegroup_q (not used yet).
  symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
            (not used yet).
  symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
  timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise
            (not in use yet).
  dtset = OPTIONAL, dataset structure containing the input variable of the
            calculation. This is required to use the k-interpolation routine.
  eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
            at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
            fine grid. This information is read from the WF dense k-grid file.
  hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
            variable is required for the k-interpolation routine.
  hdr0     = OPTIONAL, header of the GS WF file of the corse k-point grid. This
            variable is required for the k-interpolation routine.

OUTPUT

  eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the
            second-order eigenvalues: E^{(2),diag}_{k,q,j}.
  eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
            electron lifetimes.

PARENTS

      dfpt_looppert

CHILDREN

      distrb2,dotprod_g,kptfine_av,smeared_delta,timab,wrtout,xmpi_sum

SOURCE

 899 subroutine eig2stern(occ,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0,eigenq,&
 900 &  eigen1,eig2nkq,elph2_imagden,esmear,gh0c1_pert,gh1c_pert,ieig2rf,istwfk_pert,&
 901 &  mband,mk1mem,mpert,npert,mpi_enreg,mpw1,nkpt_rbz,npwar1,nspinor,nsppol,smdelta,&
 902 &  dtset,eigbrd,eigenq_fine,hdr_fine,hdr0)
 903 
 904 
 905 !This section has been created automatically by the script Abilint (TD).
 906 !Do not modify the following lines by hand.
 907 #undef ABI_FUNC
 908 #define ABI_FUNC 'eig2stern'
 909 !End of the abilint section
 910 
 911  implicit none
 912 
 913 !Arguments ------------------------------------
 914 !scalars
 915  integer,intent(in) :: bdeigrf,dim_eig2nkq,dim_eig2rf,ieig2rf,mband,mk1mem,mpert,mpw1,nkpt_rbz
 916  integer,intent(in) :: npert,nspinor,nsppol,smdelta
 917  real(dp),intent(in) :: elph2_imagden,esmear
 918  type(MPI_type),intent(inout) :: mpi_enreg
 919 !arrays
 920  integer,intent(in) :: clflg(3,mpert)
 921  integer,intent(in) :: istwfk_pert(nkpt_rbz,3,mpert)
 922  integer,intent(in) :: npwar1(nkpt_rbz,mpert)
 923  real(dp),intent(in) :: cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
 924  real(dp),intent(in) :: gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
 925  real(dp),intent(in) :: gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert)
 926  real(dp),intent(inout) :: eigen0(nkpt_rbz*mband*nsppol)
 927  real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert)
 928  real(dp),intent(inout) :: eigenq(nkpt_rbz*mband*nsppol)
 929  real(dp),intent(out) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq)
 930  real(dp),intent(out),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)
 931  real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:)
 932  real(dp), intent(in) :: occ(mband*nkpt_rbz*nsppol)
 933  type(dataset_type), intent(in) :: dtset
 934  type(hdr_type),intent(in),optional :: hdr_fine,hdr0
 935 
 936 !Local variables-------------------------------
 937 !tolerance for non degenerated levels
 938 !scalars
 939  integer :: band2tot_index,band_index,bandtot_index,iband,icg2,idir1,idir2
 940  integer :: ikpt,ipert1,ipert2,isppol,istwf_k,jband,npw1_k,nkpt_sub,ikpt2
 941 !integer :: ipw
 942  integer :: master,me,spaceworld,ierr
 943 !real(dp),parameter :: etol=1.0d-3
 944  real(dp),parameter :: etol=1.0d-6
 945 !real(dp),parameter :: etol=zero
 946  real(dp) :: ar,ai,deltae,den,dot2i,dot2r,dot3i,dot3r,doti,dotr,eig1_i1,eig1_i2
 947  real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av
 948  real(dp) :: wgt_int
 949  real(dp) :: eig2_diar,eigbrd_i,eigbrd_r
 950  character(len=500) :: message
 951  character(len=500) :: msg
 952 !DBSP
 953 ! character(len=300000) :: message2
 954 !END
 955  logical :: test_do_band
 956 !arrays
 957  integer, allocatable :: nband_rbz(:),icg2_rbz(:,:)
 958  integer,pointer      :: kpt_fine_sub(:)
 959  real(dp)             :: tsec(2)
 960  real(dp),allocatable :: cwavef(:,:),cwavef2(:,:),center(:),eigen0tmp(:),eigenqtmp(:)
 961  real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol)
 962  real(dp),allocatable :: gh(:,:),gh1(:,:),ghc(:,:)
 963  real(dp),allocatable :: smdfun(:,:)
 964  real(dp),pointer     :: wgt_sub(:)
 965 
 966 ! *********************************************************************
 967 
 968 !Init parallelism
 969  master =0
 970  spaceworld=mpi_enreg%comm_cell
 971  me=mpi_enreg%me_kpt
 972 !DEBUG
 973 !write(std_out,*)' eig2stern : enter '
 974 !write(std_out,*)' mpw1=',mpw1
 975 !write(std_out,*)' mband=',mband
 976 !write(std_out,*)' nsppol=',nsppol
 977 !write(std_out,*)' nkpt_rbz=',nkpt_rbz
 978 !write(std_out,*)' npert=',npert
 979 !ENDDEBUG
 980 
 981 !Init interpolation method
 982  if(present(eigenq_fine))then
 983    ABI_ALLOCATE(center,(3))
 984  end if
 985 
 986  call timab(148,1,tsec)
 987 
 988  if(nsppol==2)then
 989    message = 'nsppol=2 is still under development. Be careful when using it ...'
 990    MSG_COMMENT(message)
 991  end if
 992 
 993  band2tot_index =0
 994  bandtot_index=0
 995  band_index=0
 996 
 997 !Add scissor shift to eigenenergies
 998  if (dtset%dfpt_sciss > tol6 ) then
 999    write(msg,'(a,f7.3,2a)')&
1000 &   ' A scissor operator of ',dtset%dfpt_sciss*Ha_eV,' [eV] has been applied to the eigenenergies',ch10
1001    call wrtout(std_out,msg,'COLL')
1002    call wrtout(ab_out,msg,'COLL')
1003    ABI_ALLOCATE(eigen0tmp,(nkpt_rbz*mband*nsppol))
1004    ABI_ALLOCATE(eigenqtmp,(nkpt_rbz*mband*nsppol))
1005    eigen0tmp =   eigen0(:)
1006    eigenqtmp =   eigenq(:)
1007    eigen0 = zero
1008    eigenq = zero
1009  end if
1010 
1011  if(ieig2rf > 0) then
1012    eig2nkq(:,:,:,:,:,:,:) = zero
1013  end if
1014  if(present(eigbrd))then
1015    eigbrd(:,:,:,:,:,:,:) = zero
1016  end if
1017 
1018  if(xmpi_paral==1) then
1019    ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol))
1020    ABI_ALLOCATE(nband_rbz,(nkpt_rbz*nsppol))
1021    if (allocated(mpi_enreg%my_kpttab)) then
1022      ABI_DEALLOCATE(mpi_enreg%my_kpttab)
1023    end if
1024    ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz))
1025 !  Assume the number of bands is the same for all k points.
1026    nband_rbz(:)=mband
1027    call distrb2(mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg)
1028  end if
1029 
1030  icg2=0
1031  ipert1=1 ! Suppose that the situation is the same for all perturbations
1032  ABI_ALLOCATE(icg2_rbz,(nkpt_rbz,nsppol))
1033  do isppol=1,nsppol
1034    do ikpt=1,nkpt_rbz
1035      icg2_rbz(ikpt,isppol)=icg2
1036      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) cycle
1037      icg2 = icg2 + npwar1(ikpt,ipert1)*nspinor*mband
1038    end do
1039  end do
1040 
1041  do isppol=1,nsppol
1042    do ikpt =1,nkpt_rbz
1043 
1044      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then
1045        band2tot_index = band2tot_index + 2*mband**2
1046        bandtot_index = bandtot_index + mband
1047        cycle
1048      end if
1049 
1050      if(present(eigenq_fine))then
1051        write(std_out,*) 'Start of the energy denominator interpolation method.'
1052        nkpt_sub = 0
1053 !      center is the k+q point around which we will average the kpt_fine
1054        center = hdr0%kptns(:,ikpt)+ dtset%qptn(:)
1055 
1056        call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub)
1057        write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid &
1058 &       around the k+Q point ',center,' is:',nkpt_sub
1059        write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub)
1060      end if
1061 
1062 !    Add scissor shift to eigenenergies
1063      if (dtset%dfpt_sciss > tol6 ) then
1064        do iband=1,mband
1065          if (occ(iband+bandtot_index) < tol6) then
1066            eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) + dtset%dfpt_sciss
1067            eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) + dtset%dfpt_sciss
1068          else
1069            eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index)
1070            eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index)
1071          end if
1072        end do
1073      end if
1074 
1075 
1076      if(smdelta >0) then   !broadening
1077        if(.not.allocated(smdfun))  then
1078          ABI_ALLOCATE(smdfun,(mband,mband))
1079        end if
1080        smdfun(:,:) = zero
1081        do iband=1,mband
1082          eigen(iband) = eigen0(iband+bandtot_index)
1083          eigen_prime(iband) =eigenq(iband+bandtot_index)
1084        end do
1085        if(esmear>tol6) then
1086          call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun)
1087        end if
1088      end if
1089      icg2=icg2_rbz(ikpt,isppol)
1090 
1091      ipert1=1 ! Suppose all perturbations lead to the same number of planewaves
1092      npw1_k = npwar1(ikpt,ipert1)
1093      ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor))
1094      ABI_ALLOCATE(cwavef2,(2,npw1_k*nspinor))
1095      ABI_ALLOCATE(gh,(2,npw1_k*nspinor))
1096      ABI_ALLOCATE(gh1,(2,npw1_k*nspinor))
1097      ABI_ALLOCATE(ghc,(2,npw1_k*nspinor))
1098 
1099      do iband=1,bdeigrf
1100 
1101 !      If the k point and band belong to me, compute the contribution
1102        test_do_band=.true.
1103        if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false.
1104 
1105        if(test_do_band)then
1106 
1107          do ipert1=1,npert
1108 
1109            do idir1=1,3
1110              if(clflg(idir1,ipert1)==0)cycle
1111              istwf_k = istwfk_pert(ikpt,idir1,ipert1)
1112 
1113              do ipert2=1,npert
1114                do idir2=1,3
1115                  if(clflg(idir2,ipert2)==0)cycle
1116 
1117                  eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero
1118 
1119                  do jband=1,mband
1120                    eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1121                    eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
1122                    eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1123                    eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
1124 !                  If no interpolation, fallback on to the previous
1125 !                  implementation
1126                    if(.not. present(eigenq_fine))then
1127                      deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index)
1128                    end if
1129                    ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
1130                    ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
1131 
1132 !                  Sum over all active space to retrieve the diagonal gauge
1133                    if(ieig2rf == 1 .or. ieig2rf ==2 ) then
1134 !                    if(abs(deltae)>etol) then ! This is commented because
1135 !                    there is no problem with divergencies with elph2_imag != 0
1136                      if( present(eigenq_fine))then
1137                        den_av = zero
1138                        wgt_int = zero
1139                        do ikpt2=1,nkpt_sub
1140                          deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)&
1141 &                         -eigen0(iband+bandtot_index)
1142                          den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2)
1143                          wgt_int = wgt_int+wgt_sub(ikpt2)
1144                        end do
1145                        den = den_av/wgt_int
1146                      else
1147                        if(abs(elph2_imagden) < etol) then
1148                          if(abs(deltae)>etol) then
1149                            den=-one/(deltae**2+elph2_imagden**2)
1150                          else
1151                            den= zero
1152                          end if
1153                        else
1154                          den=-one/(deltae**2+elph2_imagden**2)
1155                        end if
1156                      end if
1157 
1158 !                    The following should be the most general implementation of the presence of elph2_imagden
1159 !                    eig2_diar=eig2_diar+(ar*deltae+ai*elph2_imagden)*den
1160 !                    eig2_diai=eig2_diai+(ai*deltae-ar*elph2_imagden)*den
1161 !                    This gives back the implementation without elph2_imagden
1162 !                    eig2_diar=eig2_diar+ar*deltae*den
1163 !                    eig2_diai=eig2_diai+ai*deltae*den
1164 !                    This is what Samuel had implemented
1165 !                    eig2_diar=eig2_diar+ar*deltae*den
1166 !                    eig2_diai=eig2_diai+ai*elph2_imagden*den
1167 !                    Other possibility : throw away the broadening part, that is actually treated separately.
1168                      if( present(eigenq_fine))then
1169                        eig2_diar=eig2_diar+ar*den
1170                        eig2_diai=eig2_diai+ai*den
1171                      else
1172                        eig2_diar=eig2_diar+ar*deltae*den
1173                        eig2_diai=eig2_diai+ai*deltae*den
1174 !DBSP
1175 !                       if (iband+band_index==2 .and. ikpt==1 .and. idir1==1 .and. ipert1==1 .and. idir2==1 .and. ipert2==1) then
1176 !                         write(message2,*) 'eig2_diar1=',eig2_diar,' ar=',ar,' deltae=',deltae,' den=',den
1177 !                         call wrtout(std_out,message2,'PERS')
1178 !                       endif
1179 !END
1180 
1181                      end if
1182                    end if ! ieig2rf==1 or 2
1183 
1184                    if(present(eigbrd))then
1185                      if(smdelta >0) then   !broadening
1186                        eigbrd_r = eigbrd_r + ar*smdfun(iband,jband)
1187                        eigbrd_i = eigbrd_i + ai*smdfun(iband,jband)
1188                      end if
1189                    end if
1190 
1191                  end do !jband
1192 
1193 !                Add the contribution of non-active bands, if DFPT calculation (= Sternheimer)
1194                  if(ieig2rf == 1 .or. ieig2rf ==3 .or. ieig2rf ==4 .or. ieig2rf==5 ) then
1195 !                  if(ieig2rf == 1   ) then
1196 
1197                    dotr=zero ; doti=zero
1198                    dot2r=zero ; dot2i=zero
1199                    dot3r=zero ; dot3i=zero
1200 
1201 
1202                    cwavef(:,:) = cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2)
1203                    cwavef2(:,:)= cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
1204                    gh1(:,:)    = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
1205                    gh(:,:)     = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2)
1206                    ghc(:,:)    = gh0c1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1)
1207 
1208 !                  The first two dotprod corresponds to:  <Psi(1)nkq|H(1)k+q,k|Psi(0)nk> and <Psi(0)nk|H(1)k,k+q|Psi(1)nkq>
1209 !                  They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space.
1210                    call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,cwavef,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1211                    call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,gh,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1212 
1213 !                  This dotprod corresponds to : <Psi(1)nkq|H(0)k+q- E(0)nk|Psi(1)nkq>
1214 !                  It is calculated using wavefunctions that are orthogonal to the active space.
1215 !                  Should work for metals. (But adiabatic approximation is bad in this case...)
1216                    call dotprod_g(dot3r,dot3i,istwf_k,npw1_k*nspinor,2,cwavef,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1217 
1218                    eig2_diar= eig2_diar + dotr + dot2r + dot3r
1219                    eig2_diai= eig2_diai + doti + dot2i + dot3i
1220 
1221                  end if
1222 
1223 !                Store the contribution
1224                  if(ieig2rf > 0) then
1225                    eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diar
1226                    eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diai
1227                  end if
1228 
1229                  if(present(eigbrd))then
1230                    if(smdelta >0) then   !broadening
1231                      eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r
1232                      eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i
1233                    end if
1234                  end if
1235 
1236                end do !idir2
1237              end do !ipert2
1238            end do  !idir1
1239          end do   !ipert1
1240 
1241        end if ! Selection of processor
1242 
1243      end do !iband
1244 
1245      ABI_DEALLOCATE(cwavef)
1246      ABI_DEALLOCATE(cwavef2)
1247      ABI_DEALLOCATE(gh)
1248      ABI_DEALLOCATE(gh1)
1249      ABI_DEALLOCATE(ghc)
1250      band2tot_index = band2tot_index + 2*mband**2
1251      bandtot_index = bandtot_index + mband
1252 
1253      if(present(eigenq_fine))then
1254        ABI_DEALLOCATE(kpt_fine_sub) ! Deallocate the variable
1255        ABI_DEALLOCATE(wgt_sub)
1256      end if
1257 
1258    end do    !ikpt
1259    band_index = band_index + mband
1260  end do !isppol
1261 
1262 !Accumulate eig2nkq and/or eigbrd
1263  if(xmpi_paral==1) then
1264    if(ieig2rf == 1 .or. ieig2rf == 2) then
1265      call xmpi_sum(eig2nkq,spaceworld,ierr)
1266      if (dtset%dfpt_sciss > tol6 ) then
1267        call xmpi_sum(eigen0,spaceworld,ierr)
1268        call xmpi_sum(eigenq,spaceworld,ierr)
1269      end if
1270    end if
1271    if(present(eigbrd) .and. (ieig2rf == 1 .or. ieig2rf == 2))then
1272      if(smdelta >0) then
1273        call xmpi_sum(eigbrd,spaceworld,ierr)
1274      end if
1275    end if
1276    ABI_DEALLOCATE(nband_rbz)
1277    ABI_DEALLOCATE(mpi_enreg%proc_distrb)
1278    ABI_DEALLOCATE(mpi_enreg%my_kpttab)
1279  end if
1280 
1281  if(ieig2rf==1 .or. ieig2rf==2 ) then
1282    write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.'
1283    write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1284    do idir1=1,3
1285      do idir2=1,3
1286        ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1287        ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1288        write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1289      end do ! idir2
1290    end do ! idir1
1291  end if
1292 
1293  if(present(eigbrd))then
1294    if(smdelta >0) then   !broadening
1295      write(ab_out,'(a)')' '
1296      write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.'
1297      write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1298      do idir1=1,3
1299        do idir2=1,3
1300          ar=eigbrd(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1301          ai=eigbrd(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1302          write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1303        end do
1304      end do !nband
1305    end if
1306  end if
1307 
1308  if(allocated(smdfun))  then
1309    ABI_DEALLOCATE(smdfun)
1310  end if
1311  ABI_DEALLOCATE(icg2_rbz)
1312  if(present(eigenq_fine))then
1313    ABI_DEALLOCATE(center)
1314  end if
1315  if (dtset%dfpt_sciss > tol6 ) then
1316    ABI_DEALLOCATE(eigen0tmp)
1317    ABI_DEALLOCATE(eigenqtmp)
1318  end if
1319 
1320  call timab(148,2,tsec)
1321 
1322 !DEBUG
1323 !write(std_out,*)' eig2stern: exit'
1324 !ENDDEBUG
1325 
1326 end subroutine eig2stern

ABINIT/m_eig2d [ Modules ]

[ Top ] [ Modules ]

NAME

  m_eig2d

FUNCTION

  This module contains utilities to analyze and retrieve information
  from the second order derivative of the eigen-energies wrt
  displacements.

COPYRIGHT

 Copyright (C) 2014-2018 ABINIT group (SP, PB, XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

  dfpt_looppert

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 MODULE m_eig2d
31 
32  use defs_basis
33  use defs_datatypes
34  use defs_abitypes
35  use m_errors
36  use m_abicore
37  use m_nctk
38 #ifdef HAVE_NETCDF
39  use netcdf
40 #endif
41  use m_xmpi
42  use m_ebands
43  use m_cgtools
44 
45  use m_time,       only : timab
46  use m_fstrings,   only : strcat
47  use m_crystal,    only : crystal_init, crystal_free, crystal_t
48  use m_crystal_io, only : crystal_ncwrite
49  use m_pawtab,     only : pawtab_type
50  use m_ddb,        only : DDB_VERSION
51  use m_ddb_hdr,    only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
52  use m_double_grid,only : kptfine_av
53  use m_mpinfo,     only : distrb2, proc_distrb_cycle
54 
55  implicit none
56 
57  private
58 
59  public :: eigr2d_init              ! Main creation method of EIG2D.nc files.
60  public :: eigr2d_ncwrite           ! Dump the object into NETCDF file.
61  public :: eigr2d_free              ! Destruction method.
62  public :: fan_init                 ! Main creation method of Fan.nc files.
63  public :: fan_ncwrite              ! Dump the object into NETCDF file.
64  public :: fan_free                 ! Destruction method.
65  public :: gkk_init                 ! Main creation method of GKK.nc files.
66  public :: gkk_ncwrite              ! Dump the object into NETCDF file.
67  public :: gkk_free                 ! Destruction method.
68 
69  public :: eig2tot                  ! This routine calculates the second-order eigenvalues.
70  public :: outbsd                   ! output bsd file for one perturbation (used for elphon calculations in anaddb)
71  public :: eig2stern
72  public :: elph2_fanddw             ! Calculates the zero-point motion corrections

m-eig2d/smeared_delta [ Functions ]

[ Top ] [ Functions ]

NAME

 smeared_delta

FUNCTION

 This subroutine calculates the smeared delta that weights matrix elements:
 \delta (\epsilon_{kn}-\epsilon_{k+Q,n'})

INPUTS

 eigen0(mband*nsppol) : Eigenvalues at point K
 eigenq(mband*nsppol)  : Eigenvalues at point K+Q
 mband : maximum number of bands
 smdelta : Variable controlling the smearinf scheme

OUTPUT

 smdfunc(mband,mband) : Smeared delta function weight corresponding to \delta(\epsilon_{n,k} - \epsilon_{n',k+Q})

PARENTS

      eig2stern,eig2tot

CHILDREN

SOURCE

2072 subroutine smeared_delta(eigen0,eigenq,esmear,mband,smdelta,smdfunc)
2073 
2074 
2075 !This section has been created automatically by the script Abilint (TD).
2076 !Do not modify the following lines by hand.
2077 #undef ABI_FUNC
2078 #define ABI_FUNC 'smeared_delta'
2079 !End of the abilint section
2080 
2081  implicit none
2082 
2083 !Arguments ------------------------------------
2084 !scalars
2085  integer,intent(in) :: mband,smdelta
2086 !arrays
2087  real(dp),intent(in) :: eigen0(mband),eigenq(mband),esmear
2088  real(dp),intent(out) :: smdfunc(mband,mband)
2089 
2090 !Local variables-------------------------------
2091 !tolerance for non degenerated levels
2092 !scalars
2093  integer :: ii,jj
2094  real(dp) :: aa,dsqrpi,gauss,xx
2095  character(len=500) :: message
2096 
2097 ! *********************************************************************
2098 
2099 
2100 !---------------------------------------------------------
2101 !Ordinary (unique) smearing function
2102 !---------------------------------------------------------
2103 
2104  if(smdelta==1)then
2105 
2106 !  Fermi-Dirac
2107    do ii=1,mband
2108      do jj= 1,mband
2109        xx= ( eigen0(ii) - eigenq(jj) )/esmear
2110        smdfunc(ii,jj)=0.25_dp/esmear/(cosh(xx/2.0_dp))**2
2111      end do
2112    end do
2113 
2114  else if(smdelta==2 .or. smdelta==3)then
2115 
2116 !  Cold smearing of Marzari, two values of the "a" parameter being possible
2117 !  first value gives minimization of the bump
2118    if(smdelta==2)aa=-.5634
2119 !  second value gives monotonic occupation function
2120    if(smdelta==3)aa=-.8165
2121 
2122    dsqrpi=1.0_dp/sqrt(pi)
2123    do ii=1,mband
2124      do jj=1,mband
2125        xx= ( eigen0(ii) - eigenq(jj) ) / esmear
2126        gauss=dsqrpi*exp(-xx**2)/esmear
2127        smdfunc(ii,jj)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx)))
2128      end do
2129    end do
2130 
2131  else if(smdelta==4)then
2132 
2133 !  First order Hermite-Gaussian of Paxton and Methfessel
2134    dsqrpi=1.0_dp/sqrt(pi)
2135    do ii=1,mband
2136      do jj=1,mband
2137        xx= ( eigen0(ii) - eigenq (jj) ) / esmear
2138        smdfunc(ii,jj)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)/esmear
2139      end do
2140    end do
2141 
2142  else if(smdelta==5)then
2143 
2144 !  Gaussian smearing
2145    dsqrpi=1.0_dp/sqrt(pi)
2146    do ii=1,mband
2147      do jj=1,mband
2148        xx= ( eigen0(ii) - eigenq (jj) ) / esmear
2149        smdfunc(ii,jj)=dsqrpi*exp(-xx**2)/esmear
2150      end do
2151    end do
2152 
2153  else
2154    write(message, '(a,i0,a)' )'  Smdelta= ',smdelta,' is not allowed in smdfunc'
2155    MSG_BUG(message)
2156  end if
2157 
2158 end subroutine smeared_delta

m_eig2d/eig2tot [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 eig2tot

FUNCTION

 This routine calculates the second-order eigenvalues.
 The output eig2nkq is this quantity for the input k points.

INPUTS

  bdeigrf = number of bands for which to calculate the second-order eigenvalues.
  clflg(3,mpert)= array on calculated perturbations for eig2rf.
  dim_eig2nkq = 1 if eig2nkq is to be computed.
  eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the
            electronic eigenvalues (optional).
  eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points:
            <k,n'|H(0)|k,n'> (hartree).
  eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
            <k+Q,n'|H(0)|k+Q,n'> (hartree).
  eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order:
            <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
  eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
            the electronic eigenvalues.
  elph2_imagden = imaginary part of the denominator of the sum-over-state expression
            for the electronic eigenenergy shift due to second-order electron-phonon
            interation.
  ieig2rf = integer for calculation type.
  indsym(4,nsym,natom) = indirect indexing array for atom labels
            (not used yet, but will be used with symmetries).
  mband = maximum number of bands.
  mpert = maximum number of perturbations.
  natom = number of atoms in the unit cell.
  npert = number of phonon perturbations, without taking into account directions:
            natom.
  nsym = number of symmetries (not used yet).
  mpi_enreg = information about MPI parallelization.
  nkpt_rbz = number of k-points for each perturbation.
  nsppol = 1 for unpolarized, 2 for spin-polarized.
  smdelta = integer controling the calculation of electron lifetimes.
  symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From littlegroup_q (not used yet).
  symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
            (not used yet).
  symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
  timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise
            (not in use yet).
  dtset = OPTIONAL, dataset structure containing the input variable of the
            calculation. This is required to use the k-interpolation routine.
  eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
            at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
            fine grid. This information is read from the WF dense k-grid file.
  hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
            variable is required for the k-interpolation routine.
  hdr0     = header of the GS WF file of the corse k-point grid.

OUTPUT

  eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the
            second-order eigenvalues: E^{(2),diag}_{k,q,j}.
  eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
            electron lifetimes.

PARENTS

      respfn

CHILDREN

      crystal_free,crystal_init,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write
      distrb2,ebands_free,ebands_init,eigr2d_free,eigr2d_init,eigr2d_ncwrite
      fan_free,fan_init,fan_ncwrite,gkk_free,gkk_init,gkk_ncwrite,kptfine_av
      outbsd,smeared_delta,timab,xmpi_sum

SOURCE

1400 subroutine eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0,eigenq,eigen1,eig2nkq,&
1401 &  elph2_imagden,esmear,ieig2rf,mband,mpert,npert,mpi_enreg,doccde,&
1402 &  nkpt_rbz,nsppol,smdelta,rprimd,dtset,occ_rbz,hdr0,eigbrd,eigenq_fine,hdr_fine)
1403 
1404 
1405 !This section has been created automatically by the script Abilint (TD).
1406 !Do not modify the following lines by hand.
1407 #undef ABI_FUNC
1408 #define ABI_FUNC 'eig2tot'
1409 !End of the abilint section
1410 
1411  implicit none
1412 
1413 !Arguments ------------------------------------
1414 !scalars
1415  integer,intent(in) :: bdeigrf,dim_eig2nkq,ieig2rf,mband,mpert,natom,nkpt_rbz
1416  integer,intent(in) :: npert,nsppol,smdelta
1417  real(dp),intent(in) :: elph2_imagden,esmear
1418  type(MPI_type),intent(inout) :: mpi_enreg
1419  type(datafiles_type), intent(in) :: dtfil
1420  type(pseudopotential_type), intent(inout) :: psps
1421 !arrays
1422  type(dataset_type), intent(in) :: dtset
1423  integer,intent(in) :: clflg(3,mpert)
1424  real(dp),intent(in) :: doccde(dtset%mband*dtset%nkpt*dtset%nsppol)
1425  real(dp),intent(in) :: eigen0(nkpt_rbz*mband*nsppol)
1426  real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert)
1427  real(dp),intent(in) :: eigenq(nkpt_rbz*mband*nsppol)
1428  real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
1429  real(dp),intent(inout) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq)
1430  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1431  real(dp),intent(inout),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)
1432  real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:)
1433  type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
1434  type(hdr_type),intent(in) :: hdr0
1435  type(hdr_type),intent(in),optional :: hdr_fine
1436 
1437 !Local variables-------------------------------
1438 !tolerance for non degenerated levels
1439 !scalars
1440  integer :: band2tot_index,band_index,bantot,bandtot_index,iband,idir1,idir2
1441  integer :: ikpt,ipert1,ipert2,isppol,jband,nkpt_sub,ikpt2,unitout,ncid
1442 !integer :: ipw
1443  character(len=fnlen) :: dscrpt,fname
1444  integer :: master,me,spaceworld,ierr
1445 ! real(dp),parameter :: etol=1.0d-6
1446  real(dp),parameter :: etol=1.0d-7
1447 !real(dp),parameter :: etol=zero
1448  real(dp) :: ar,ai,deltae,den,eig1_i1,eig1_i2,eigen_corr
1449  real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av
1450  real(dp) :: eig2_diar,eigbrd_i,eigbrd_r,wgt_int
1451  character(len=500) :: message
1452  logical :: remove_inv,test_do_band
1453  type(crystal_t) :: Crystal
1454  type(ebands_t)  :: Bands
1455  type(eigr2d_t)  :: eigr2d,eigi2d
1456  type(fan_t)     :: fan2d
1457  type(gkk_t)     :: gkk2d
1458  type(ddb_hdr_type) :: ddb_hdr
1459 !arrays
1460  integer, allocatable :: nband_rbz(:)
1461  integer,pointer      :: kpt_fine_sub(:)
1462  real(dp)             :: tsec(2)
1463  real(dp),allocatable :: center(:)
1464  real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol)
1465  real(dp),allocatable :: fan(:,:,:,:,:,:,:)
1466  real(dp),allocatable :: gkk(:,:,:,:,:)
1467  real(dp),allocatable :: eig2nkq_tmp(:,:,:,:,:,:,:)
1468  real(dp),allocatable :: smdfun(:,:)
1469  real(dp),pointer     :: wgt_sub(:)
1470 
1471 ! *********************************************************************
1472 
1473 !Init parallelism
1474  master =0
1475  spaceworld=mpi_enreg%comm_cell
1476  me=mpi_enreg%me_kpt
1477 !DEBUG
1478 !write(std_out,*)' eig2tot : enter '
1479 !write(std_out,*)' mband=',mband
1480 !write(std_out,*)' nsppol=',nsppol
1481 !write(std_out,*)' nkpt_rbz=',nkpt_rbz
1482 !write(std_out,*)' npert=',npert
1483 !ENDDEBUG
1484 
1485 !Init interpolation method
1486  if(present(eigenq_fine))then
1487    ABI_ALLOCATE(center,(3))
1488  end if
1489 
1490  call timab(148,1,tsec)
1491 
1492  if(nsppol==2)then
1493    message = 'nsppol=2 is still under development. Be careful when using it ...'
1494    MSG_COMMENT(message)
1495  end if
1496 
1497  band2tot_index =0
1498  bandtot_index=0
1499  band_index=0
1500 
1501  if(xmpi_paral==1) then
1502    ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol))
1503    ABI_ALLOCATE(nband_rbz,(nkpt_rbz*nsppol))
1504    if (allocated(mpi_enreg%my_kpttab)) then
1505      ABI_DEALLOCATE(mpi_enreg%my_kpttab)
1506    end if
1507    ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz))
1508 !  Assume the number of bands is the same for all k points.
1509    nband_rbz(:)=mband
1510    call distrb2(mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg)
1511  end if
1512 
1513  if(ieig2rf == 4 ) then
1514    ABI_STAT_ALLOCATE(fan,(2*mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq,mband), ierr)
1515    ABI_CHECK(ierr==0, "out-of-memory in fan")
1516    fan(:,:,:,:,:,:,:) = zero
1517    ABI_STAT_ALLOCATE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr)
1518    ABI_CHECK(ierr==0, "out-of-memory in eig2nkq_tmp")
1519    eig2nkq_tmp(:,:,:,:,:,:,:) = zero
1520 !  This is not efficient because double the memory. Alternative: use buffer and
1521 !  print part by part.
1522    eig2nkq_tmp = eig2nkq
1523    if(present(eigbrd))then
1524      eigbrd(:,:,:,:,:,:,:)=zero
1525    end if
1526    eigen_corr = 0
1527  end if
1528 
1529  if(ieig2rf == 5 ) then
1530    ABI_STAT_ALLOCATE(gkk,(2*mband*nsppol,dtset%nkpt,3,natom,mband), ierr)
1531    ABI_CHECK(ierr==0, "out-of-memory in gkk")
1532    gkk(:,:,:,:,:) = zero
1533    ABI_STAT_ALLOCATE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr)
1534    ABI_CHECK(ierr==0, "out-of-memory in eig2nkq_tmp")
1535    eig2nkq_tmp(:,:,:,:,:,:,:) = zero
1536 !  This is not efficient because double the memory. Alternative: use buffer and
1537 !  print part by part.
1538    eig2nkq_tmp = eig2nkq
1539    if(present(eigbrd))then
1540      eigbrd(:,:,:,:,:,:,:)=zero
1541    end if
1542    eigen_corr = 0
1543  end if
1544 
1545  do isppol=1,nsppol
1546    do ikpt =1,nkpt_rbz
1547 
1548      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then
1549        band2tot_index = band2tot_index + 2*mband**2
1550        bandtot_index = bandtot_index + mband
1551        cycle
1552      end if
1553 
1554      if(present(eigenq_fine))then
1555        write(std_out,*) 'Start of the energy denominator interpolation method.'
1556        nkpt_sub = 0
1557 !      center is the k+q point around which we will average the kpt_fine
1558        center = hdr0%kptns(:,ikpt)+ dtset%qptn(:)
1559 
1560        call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub)
1561        write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid &
1562 &       around the k+Q point ',center,' is:',nkpt_sub
1563        write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub)
1564      end if
1565 
1566      if(smdelta >0) then   !broadening
1567        if(.not.allocated(smdfun))  then
1568          ABI_ALLOCATE(smdfun,(mband,mband))
1569        end if
1570        smdfun(:,:) = zero
1571        do iband=1,mband
1572          eigen(iband) = eigen0(iband+bandtot_index)
1573          eigen_prime(iband) =eigenq(iband+bandtot_index)
1574        end do
1575        if(esmear>tol6) then
1576          call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun)
1577        end if
1578      end if
1579 
1580      ipert1=1 ! Suppose all perturbations lead to the same number of planewaves
1581 
1582      do iband=1,bdeigrf
1583 
1584 !      If the k point and band belong to me, compute the contribution
1585        test_do_band=.true.
1586        if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false.
1587 
1588        if(test_do_band)then
1589 !        ------------------------------------------------------------------------------------------------------!
1590 !        ------- ieig2rf ==3 : Non dynamic traditional AHC theory with Sternheimer (computed in eig2stern.F90)-!
1591 !        ------------------------------------------------------------------------------------------------------!
1592 !        Note that ieig2rf==4 and ieig2rf==5 also goes into that part only for later printing of the ZPR in the ouput of abinit
1593 !        later in the code
1594          if(ieig2rf==3 .or. ieig2rf==4 .or. ieig2rf==5) then
1595            do ipert1=1,npert
1596              do idir1=1,3
1597                if(clflg(idir1,ipert1)==0) cycle
1598                do ipert2=1,npert
1599                  do idir2=1,3
1600                    if(clflg(idir2,ipert2)==0)cycle
1601                    eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero
1602                    do jband=1,mband
1603                      eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1604                      eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
1605                      eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1606                      eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
1607 !                    If no interpolation, fallback on to the previous
1608 !                    implementation
1609                      if(.not. present(eigenq_fine))then
1610                        deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index)
1611                      end if
1612                      ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
1613                      ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
1614 
1615 !                    Sum over all active space to retrieve the diagonal gauge
1616 !                    if(abs(deltae)>etol) then ! This is commented because
1617 !                    there is no problem with divergencies with elph2_imag != 0
1618                      if( present(eigenq_fine))then
1619                        den_av = zero
1620                        wgt_int = zero
1621                        do ikpt2=1,nkpt_sub
1622                          deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)&
1623 &                         -eigen0(iband+bandtot_index)
1624                          den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2)
1625                          wgt_int = wgt_int+wgt_sub(ikpt2)
1626                        end do
1627                        den = den_av/wgt_int
1628                      else
1629                        if(abs(elph2_imagden) < etol) then
1630                          if(abs(deltae)>etol) then
1631                            den=-one/(deltae**2+elph2_imagden**2)
1632                          else
1633                            den= zero
1634                          end if
1635                        else
1636                          den=-one/(deltae**2+elph2_imagden**2)
1637                        end if
1638                      end if
1639 
1640                      if( present(eigenq_fine))then
1641                        eig2_diar=eig2_diar+ar*den
1642                        eig2_diai=eig2_diai+ai*den
1643                      else
1644                        eig2_diar=eig2_diar+ar*deltae*den
1645                        eig2_diai=eig2_diai+ai*deltae*den
1646                      end if
1647 
1648                      if(present(eigbrd))then
1649                        if(smdelta >0) then   !broadening
1650                          eigbrd_r = eigbrd_r + ar*smdfun(iband,jband)
1651                          eigbrd_i = eigbrd_i + ai*smdfun(iband,jband)
1652                        end if
1653                      end if
1654                    end do !jband
1655 
1656 !                  Store the contribution
1657                    eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
1658 &                   eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diar
1659                    eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
1660 &                   eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diai
1661 
1662                    if(present(eigbrd))then
1663                      if(smdelta >0) then   !broadening
1664                        eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r
1665                        eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i
1666                      end if
1667                    end if
1668 
1669                  end do !idir2
1670                end do !ipert2
1671              end do  !idir1
1672            end do   !ipert1
1673          end if !ieig2rf 3
1674 
1675 !        -------------------------------------------------------------------------------------------!
1676 !        ------- ieig2rf ==4  Dynamic AHC using second quantization and Sternheimer from eig2stern -!
1677 !        -------------------------------------------------------------------------------------------!
1678          if(ieig2rf ==4 ) then
1679            do ipert1=1,npert
1680              do idir1=1,3
1681                if(clflg(idir1,ipert1)==0) cycle
1682                do ipert2=1,npert
1683                  do idir2=1,3
1684                    if(clflg(idir2,ipert2)==0)cycle
1685                    do jband=1,mband
1686                      eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1687                      eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
1688                      eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1689                      eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
1690                      ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
1691                      ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
1692 !                  Store the contribution
1693                      fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = &
1694 &                     fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ar
1695                      fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = &
1696 &                     fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ai
1697                    end do !jband
1698                  end do !idir2
1699                end do !ipert2
1700              end do  !idir1
1701            end do   !ipert1
1702          end if !ieig2rf 4
1703 !        --------------------------------------------------------------------------------!
1704 !        ------- ieig2rf ==5  Dynamic AHC with Sternheimer from eig2stern but print GKK -!
1705 !        --------------------------------------------------------------------------------!
1706          if(ieig2rf ==5 ) then
1707            do ipert1=1,npert
1708              do idir1=1,3
1709                if(clflg(idir1,ipert1)==0) cycle
1710                do jband=1,mband
1711                  eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1712                  eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
1713 !              Store the contribution
1714                  gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) = &
1715 &                 gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) + eig1_r1
1716                  gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) = &
1717 &                 gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) + eig1_i1
1718                end do !jband
1719              end do  !idir1
1720            end do   !ipert1
1721          end if !ieig2rf 5
1722        end if ! Selection of processor
1723      end do !iband
1724 
1725      band2tot_index = band2tot_index + 2*mband**2
1726      bandtot_index = bandtot_index + mband
1727 
1728      if(present(eigenq_fine))then
1729        ABI_DEALLOCATE(kpt_fine_sub) ! Deallocate the variable
1730        ABI_DEALLOCATE(wgt_sub)
1731      end if
1732 
1733    end do    !ikpt
1734    band_index = band_index + mband
1735  end do !isppol
1736 
1737 !Accumulate eig2nkq and/or eigbrd
1738  if(xmpi_paral==1) then
1739    if(ieig2rf == 3) then
1740      call xmpi_sum(eig2nkq,spaceworld,ierr)
1741    end if
1742    if(ieig2rf == 4) then
1743      call xmpi_sum(eig2nkq,spaceworld,ierr)
1744      call xmpi_sum(eig2nkq_tmp,spaceworld,ierr)
1745      call xmpi_sum(fan,spaceworld,ierr)
1746    end if
1747    if(ieig2rf == 5) then
1748      call xmpi_sum(eig2nkq,spaceworld,ierr)
1749      call xmpi_sum(eig2nkq_tmp,spaceworld,ierr)
1750      call xmpi_sum(gkk,spaceworld,ierr)
1751    end if
1752    if(present(eigbrd) .and. (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5))then
1753      if(smdelta >0) then
1754        call xmpi_sum(eigbrd,spaceworld,ierr)
1755      end if
1756    end if
1757    ABI_DEALLOCATE(nband_rbz)
1758    ABI_DEALLOCATE(mpi_enreg%proc_distrb)
1759    ABI_DEALLOCATE(mpi_enreg%my_kpttab)
1760  end if
1761 
1762  if(ieig2rf > 2) then
1763    write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.'
1764    write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1765    band_index = 0
1766    do isppol=1,dtset%nsppol
1767      do idir1=1,3
1768        do idir2=1,3
1769          ar=eig2nkq(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1770          ai=eig2nkq(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1771          write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1772        end do ! idir2
1773      end do ! idir1
1774      band_index = band_index + mband
1775      write(ab_out,'(a)')' '
1776    end do
1777  end if
1778 
1779  if(present(eigbrd))then
1780    if(smdelta >0) then   !broadening
1781      write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.'
1782      write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
1783      band_index = 0
1784      do isppol=1,dtset%nsppol
1785        do idir1=1,3
1786          do idir2=1,3
1787            ar=eigbrd(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
1788            ai=eigbrd(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
1789            write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
1790          end do
1791        end do
1792        band_index = band_index + mband
1793        write(ab_out,'(a)')' '
1794      end do
1795    end if
1796  end if
1797 
1798  if(allocated(smdfun))  then
1799    ABI_DEALLOCATE(smdfun)
1800  end if
1801  if(present(eigenq_fine))then
1802    ABI_DEALLOCATE(center)
1803  end if
1804 
1805  master=0
1806  if (me==master) then
1807 !  print _EIGR2D file for this perturbation in the case of ieig2rf 3 or 4 or 5
1808    if (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5) then
1809 
1810      dscrpt=' Note : temporary (transfer) database '
1811      unitout = dtfil%unddb
1812 
1813      call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
1814 &     1,xred=xred,occ=occ_rbz)
1815 
1816      call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigr2d, unitout)
1817 
1818      call ddb_hdr_free(ddb_hdr)
1819 
1820    end if
1821    if(ieig2rf == 3 ) then
1822      call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout)
1823    end if
1824    if(ieig2rf == 4 .or. ieig2rf == 5 ) then
1825      call outbsd(bdeigrf,dtset,eig2nkq_tmp,dtset%natom,nkpt_rbz,unitout)
1826    end if
1827 !  Output of the EIGR2D.nc file.
1828    fname = strcat(dtfil%filnam_ds(4),"_EIGR2D.nc")
1829 !  Crystalline structure.
1830    remove_inv=.false.
1831    if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true.
1832    call crystal_init(dtset%amu_orig(:,1),Crystal,dtset%spgroup,dtset%natom,dtset%npsp,psps%ntypat, &
1833 &   dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1834 &   dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr0%title,&
1835 &   dtset%symrel,dtset%tnons,dtset%symafm)
1836 !  Electronic band energies.
1837    bantot= dtset%mband*dtset%nkpt*dtset%nsppol
1838    call ebands_init(bantot,Bands,dtset%nelect,doccde,eigen0,hdr0%istwfk,hdr0%kptns,&
1839 &   hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
1840 &   hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,&
1841 &   hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, &
1842 &   hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk)
1843 
1844 !  Second order derivative EIGR2D (real and Im)
1845    if(ieig2rf == 3 ) then
1846      call eigr2d_init(eig2nkq,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
1847    end if
1848    if(ieig2rf == 4 .or. ieig2rf == 5 ) then
1849      call eigr2d_init(eig2nkq_tmp,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
1850    end if
1851 #ifdef HAVE_NETCDF
1852    NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGR2D file")
1853    NCF_CHECK(crystal_ncwrite(Crystal,ncid))
1854    NCF_CHECK(ebands_ncwrite(Bands, ncid))
1855    call eigr2d_ncwrite(eigr2d,dtset%qptn(:),dtset%wtq,ncid)
1856    NCF_CHECK(nf90_close(ncid))
1857 #else
1858    ABI_UNUSED(ncid)
1859 #endif
1860 
1861 !  print _FAN file for this perturbation. Note that the Fan file will only be produced if
1862 !  abinit is compiled with netcdf.
1863    if(ieig2rf == 4 ) then
1864 !    Output of the Fan.nc file.
1865 #ifdef HAVE_NETCDF
1866      fname = strcat(dtfil%filnam_ds(4),"_FAN.nc")
1867      call fan_init(fan,fan2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
1868      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating FAN file")
1869      NCF_CHECK(crystal_ncwrite(Crystal, ncid))
1870      NCF_CHECK(ebands_ncwrite(Bands, ncid))
1871      call fan_ncwrite(fan2d,dtset%qptn(:),dtset%wtq, ncid)
1872      NCF_CHECK(nf90_close(ncid))
1873 #else
1874      MSG_ERROR("Dynamical calculation with ieig2rf 4 only work with NETCDF support.")
1875      ABI_UNUSED(ncid)
1876 #endif
1877      ABI_DEALLOCATE(fan)
1878      ABI_DEALLOCATE(eig2nkq_tmp)
1879    end if
1880 !  print _GKK.nc file for this perturbation. Note that the GKK file will only be produced if
1881 !  abinit is compiled with netcdf.
1882    if(ieig2rf == 5 ) then
1883 !    Output of the GKK.nc file.
1884 #ifdef HAVE_NETCDF
1885      fname = strcat(dtfil%filnam_ds(4),"_GKK.nc")
1886      call gkk_init(gkk,gkk2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom,3)
1887      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
1888      NCF_CHECK(crystal_ncwrite(Crystal, ncid))
1889      NCF_CHECK(ebands_ncwrite(Bands, ncid))
1890      call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid)
1891      NCF_CHECK(nf90_close(ncid))
1892 #else
1893      MSG_ERROR("Dynamical calculation with ieig2rf 5 only work with NETCDF support.")
1894      ABI_UNUSED(ncid)
1895 #endif
1896      ABI_DEALLOCATE(gkk)
1897      ABI_DEALLOCATE(eig2nkq_tmp)
1898    end if
1899 !  print _EIGI2D file for this perturbation
1900    if (ieig2rf /= 5 ) then
1901      if(smdelta>0) then
1902        unitout = dtfil%unddb
1903        dscrpt=' Note : temporary (transfer) database '
1904 
1905        call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
1906 &       1,xred=xred,occ=occ_rbz)
1907 
1908        call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigi2d, unitout)
1909 
1910        call ddb_hdr_free(ddb_hdr)
1911 
1912        call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout)
1913 
1914 !      Output of the EIGI2D.nc file.
1915        fname = strcat(dtfil%filnam_ds(4),"_EIGI2D.nc")
1916 !      Broadening EIGI2D (real and Im)
1917        call eigr2d_init(eigbrd,eigi2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
1918 #ifdef HAVE_NETCDF
1919        NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGI2D file")
1920        NCF_CHECK(crystal_ncwrite(Crystal, ncid))
1921        NCF_CHECK(ebands_ncwrite(Bands, ncid))
1922        call eigr2d_ncwrite(eigi2d,dtset%qptn(:),dtset%wtq,ncid)
1923        NCF_CHECK(nf90_close(ncid))
1924 #else
1925        ABI_UNUSED(ncid)
1926 #endif
1927      end if !smdelta
1928    end if
1929  end if
1930 
1931  if (allocated(fan)) then
1932    ABI_DEALLOCATE(fan)
1933  end if
1934  if (allocated(eig2nkq_tmp)) then
1935    ABI_DEALLOCATE(eig2nkq_tmp)
1936  end if
1937  if (allocated(gkk)) then
1938    ABI_DEALLOCATE(gkk)
1939  end if
1940 
1941  call crystal_free(Crystal)
1942  call ebands_free(Bands)
1943  call eigr2d_free(eigr2d)
1944  call eigr2d_free(eigi2d)
1945  call fan_free(fan2d)
1946  call gkk_free(gkk2d)
1947 
1948 
1949  call timab(148,2,tsec)
1950 !DEBUG
1951 !write(std_out,*)' eig2tot: exit'
1952 !ENDDEBUG
1953 
1954 end subroutine eig2tot

m_eig2d/eigr2d_free [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 eigr2d_free

FUNCTION

 Deallocates the components of the eigr2d_t structured datatype

INPUTS

  eigr2d<eigr2d_t>=The data type to be deallocated.

OUTPUT

  Deallocate the dynamic arrays in the ebands_t type.
  (only deallocate)

PARENTS

      dfpt_looppert,eig2tot

CHILDREN

SOURCE

356 subroutine eigr2d_free(eigr2d)
357 
358 
359 !This section has been created automatically by the script Abilint (TD).
360 !Do not modify the following lines by hand.
361 #undef ABI_FUNC
362 #define ABI_FUNC 'eigr2d_free'
363 !End of the abilint section
364 
365  implicit none
366 
367 !Arguments ------------------------------------
368 !scalars
369  type(eigr2d_t),intent(inout) :: eigr2d
370 ! *************************************************************************
371 DBG_ENTER("COLL")
372 
373 !Deallocate all components of bstruct
374 
375  if (allocated(eigr2d%eigr2d)) then
376    ABI_FREE(eigr2d%eigr2d)
377  end if
378 
379  DBG_EXIT("COLL")
380 
381 end subroutine eigr2d_free

m_eig2d/eigr2d_init [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 eigr2d_init

FUNCTION

 This subroutine initializes the eigr2d_t structured datatype

INPUTS

 mbands=maximum number of bands
 nkpt=number of k points
 nsppol=1 for unpolarized, 2 for spin-polarized
 natom=number of atoms
 eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom)=second-order derivative of the
     eigenenergies wrt phononic displacements

OUTPUT

 eigr2d<eigr2d_t>=the eigr2d_t datatype

PARENTS

      dfpt_looppert,eig2tot

CHILDREN

SOURCE

198 subroutine eigr2d_init(eig2nkq,eigr2d,mband,nsppol,nkpt,natom)
199 
200 
201 !This section has been created automatically by the script Abilint (TD).
202 !Do not modify the following lines by hand.
203 #undef ABI_FUNC
204 #define ABI_FUNC 'eigr2d_init'
205 !End of the abilint section
206 
207  implicit none
208 
209 !This section has been created automatically by the script Abilint (TD).
210 !Do not modify the following lines by hand.
211 #undef ABI_FUNC
212 #define ABI_FUNC 'eigr2d_init'
213 !End of the abilint section
214 
215 !Arguments ------------------------------------
216 !scalars
217  integer,intent(in) ::mband,nsppol,nkpt,natom
218  type(eigr2d_t),intent(out) :: eigr2d
219 !arrays
220  real(dp), intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom)
221 
222 ! *************************************************************************
223 
224  eigr2d%mband = mband
225  eigr2d%nsppol = nsppol
226  eigr2d%nkpt = nkpt
227  eigr2d%natom = natom
228 
229  ABI_MALLOC(eigr2d%eigr2d  ,(2,mband*nsppol,nkpt,3,natom,3,natom))
230  eigr2d%eigr2d=eig2nkq
231 
232 end subroutine eigr2d_init

m_eig2d/eigr2d_ncwrite [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 eigr2d_ncwrite

FUNCTION

  Writes the content of a eigr2d_t object to a NETCDF file
  according to the ETSF-IO specifications.

INPUTS

  ncid =NC file handle

OUTPUT

PARENTS

      dfpt_looppert,eig2tot

CHILDREN

SOURCE

257 subroutine eigr2d_ncwrite(eigr2d,iqpt,wtq,ncid)
258 
259 
260 !This section has been created automatically by the script Abilint (TD).
261 !Do not modify the following lines by hand.
262 #undef ABI_FUNC
263 #define ABI_FUNC 'eigr2d_ncwrite'
264 !End of the abilint section
265 
266  implicit none
267 
268 !Arguments ------------------------------------
269 !scalars
270  integer,intent(in) ::ncid
271  real(dp),intent(in) :: iqpt(3),wtq
272  type(eigr2d_t),intent(in) :: eigr2d
273 
274 !Local variables-------------------------------
275 #ifdef HAVE_NETCDF
276  integer :: ncerr
277  integer :: cplex,cart_dir,one_dim
278 
279 ! *************************************************************************
280 
281  ! ==============================================
282  ! === Write the dimensions specified by ETSF ===
283  ! ==============================================
284  one_dim=1; cplex=2; cart_dir=3
285 
286  ncerr = nctk_def_dims(ncid, [&
287    nctkdim_t('max_number_of_states', eigr2d%mband),&
288    nctkdim_t('number_of_spins', eigr2d%nsppol),&
289    nctkdim_t('number_of_kpoints', eigr2d%nkpt),&
290    nctkdim_t('number_of_atoms', eigr2d%natom),&
291    nctkdim_t('number_of_cartesian_directions', cart_dir),&
292    nctkdim_t('current_one_dim', one_dim),&
293    nctkdim_t('cplex', cplex),&
294    nctkdim_t('product_mband_nsppol', eigr2d%mband*eigr2d%nsppol)], defmode=.True.)
295  NCF_CHECK(ncerr)
296 
297  ncerr = nctk_def_arrays(ncid, [&
298 &  nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'),&
299 &  nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'),&
300 &  nctkarr_t('second_derivative_eigenenergies', "dp", &
301 &    'cplex, product_mband_nsppol, number_of_kpoints, number_of_cartesian_directions, number_of_atoms,&
302 & number_of_cartesian_directions, number_of_atoms')])
303  NCF_CHECK(ncerr)
304 
305 ! Write data
306  NCF_CHECK(nctk_set_datamode(ncid))
307  NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt))
308  NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq))
309  NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies'), eigr2d%eigr2d))
310 
311 #else
312  MSG_ERROR("ETSF-IO support is not activated. ")
313 #endif
314 
315 
316 contains
317  integer function vid(vname)
318 
319 
320 !This section has been created automatically by the script Abilint (TD).
321 !Do not modify the following lines by hand.
322 #undef ABI_FUNC
323 #define ABI_FUNC 'vid'
324 !End of the abilint section
325 
326    character(len=*),intent(in) :: vname
327    vid = nctk_idname(ncid, vname)
328  end function vid
329 
330 end subroutine eigr2d_ncwrite

m_eig2d/eigr2d_t [ Types ]

[ Top ] [ m_eig2d ] [ Types ]

NAME

 eig2d_t

FUNCTION

 It contains information about the second-order derivative of the
 eigenenergies wrt atomic displacement

SOURCE

 86  type,public :: eigr2d_t
 87 
 88 ! WARNING : if you modify this datatype, please check whether there might be
 89 ! creation/destruction/copy routines,
 90 ! declared in another part of ABINIT, that might need to take into account your
 91 ! modification.
 92 
 93   integer :: mband                 ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
 94   integer :: nsppol                ! number of spin-polarization
 95   integer :: nkpt                  ! number of k points
 96   integer :: natom                 ! number of atoms
 97 
 98   real(dp),allocatable :: eigr2d(:,:,:,:,:,:,:)
 99   ! eigr2d(2,mband*nsppol,nkpt,3,natom,3,natom)
100   ! Second-order derivative of eigenergies (real,im) at each
101   ! spin,band,k-point,dir1,dir2,natom1,natom2 .
102 
103 
104  end type eigr2d_t

m_eig2d/elph2_fanddw [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 elph2_fanddw

FUNCTION

 This routine calculates the zero-point motion corrections
 due to the Fan term or to the DDW term..

INPUTS

  dim_eig2nkq=1 if eig2nkq is to be computed
  displ(2*3*natom*3*natom)=the displacements of atoms in cartesian coordinates.
  eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=one half second derivatives of the electronic eigenvalues
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  mband= maximum number of bands
  natom= number of atoms in the unit cell
  nkpt= number of k-points
  nsppol= 1 for unpolarized, 2 for spin-polarized
  option 1 for Fan term, 2 for DDW term
  phfrq(3*natom)=phonon frequencies
  (prtvol > 4) if the mode decomposition is to be printed

OUTPUT

  eigen_corr(mband*nkpt*nsppol)= T=0 correction to the electronic eigenvalues, due to the Fan term.

PARENTS

      respfn

CHILDREN

      wrtout

SOURCE

2193 subroutine elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_corr,gprimd,mband,natom,nkpt,nsppol,option,phfrq,prtvol)
2194 
2195 
2196 !This section has been created automatically by the script Abilint (TD).
2197 !Do not modify the following lines by hand.
2198 #undef ABI_FUNC
2199 #define ABI_FUNC 'elph2_fanddw'
2200 !End of the abilint section
2201 
2202  implicit none
2203 
2204 !Arguments ------------------------------------
2205 !scalars
2206  integer,intent(in) :: dim_eig2nkq,mband,natom,nkpt,nsppol,option,prtvol
2207 
2208 !arrays
2209  real(dp) :: gprimd(3,3)
2210  real(dp),intent(in) :: displ(2*3*natom*3*natom)
2211  real(dp),intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)
2212  real(dp),intent(in) :: phfrq(3*natom)
2213  real(dp),intent(out) :: eigen_corr(mband*nkpt*nsppol)
2214 
2215 !Local variables-------------------------------
2216 !scalars
2217  integer,parameter :: neigs_per_line=6
2218  integer :: iatom1,iatom2,idir1,idir2,iband,ikpt,imode,index,isppol, imin, ii
2219  real(dp) :: d_at1_dir1_re,d_at1_dir1_im
2220  real(dp) :: d_at1_dir2_re,d_at1_dir2_im
2221  real(dp) :: d_at2_dir1_re,d_at2_dir1_im
2222  real(dp) :: d_at2_dir2_re,d_at2_dir2_im
2223  real(dp) :: e2_im,e2_re
2224  real(dp), allocatable :: eigen_corr_mode(:)
2225  character(len=500) :: message
2226  character(len=20) :: eig_format, line_format
2227 !arrays
2228  real(dp) :: displ2cart(2,3,3),displ2red(2,3,3),tmp_displ2(2,3,3)
2229 
2230 ! *********************************************************************
2231 
2232 !DEBUG
2233 !write(std_out,*)' elph2_fanddw : enter '
2234 !write(std_out,*)' option=',option
2235 !ENDDEBUG
2236 
2237  if(option/=1 .and. option/=2)then
2238    write(message,'(a,i0)')' The argument option should be 1 or 2, while it is found that option=',option
2239    MSG_BUG(message)
2240  end if
2241 
2242  !printing options
2243  eig_format='f16.8'
2244  write(line_format,'(a,i1,a6,a)') '(',neigs_per_line,eig_format,')'
2245 
2246  if (prtvol > 4) then
2247    write(message,'(a,a)')ch10,' ================================================================================'
2248    call wrtout(ab_out_default,message,'COLL')
2249    if (option==1) then
2250      write(message,'(a)') ' ---- Begin Fan contributions to eigenvalues renormalization by mode ----'
2251      call wrtout(ab_out_default,message,'COLL')
2252    else if (option==2) then
2253      write(message,'(a)') ' ---- Begin DDW contributions to eigenvalues renormalization by mode ----'
2254      call wrtout(ab_out_default,message,'COLL')
2255    end if
2256  end if
2257 
2258  ABI_ALLOCATE(eigen_corr_mode,(mband*nkpt*nsppol))
2259 
2260  eigen_corr(:)=zero
2261  do imode=1,3*natom
2262    eigen_corr_mode(:)=zero
2263 
2264 !  DEBUG
2265 !  write(std_out,*)' Contribution of mode ',imode,' with frequency=',phfrq(imode),' and displacements :'
2266 !  write(std_out,'(2f14.7)' ) displ(1+2*3*natom*(imode-1):2*3*natom*imode)
2267 !  ENDDEBUG
2268 
2269    if (phfrq(imode)>tol6) then
2270      do iatom1=1,natom
2271        do iatom2=1,natom
2272 !        DEBUG
2273 !        write(std_out,*)' iatom1,iatom2=',iatom1,iatom2
2274 !        ENDDEBUG
2275 
2276          do idir1=1,3
2277            do idir2=1,3
2278 !            Compute the mean cartesian displacements
2279              d_at1_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
2280              d_at1_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
2281              d_at2_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
2282              d_at2_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
2283 
2284 !            DEBUG
2285 !            write(std_out,*)' idir1,idir2=',iatom1,iatom2,idir1,idir2
2286 !            write(std_out,'(a,4f12.5)' )' d_at1_dir1 re,d_at2_dir2 re=',d_at1_dir1_re,d_at2_dir2_re
2287 !            ENDDEBUG
2288 
2289              if(option==1)then
2290 !              Compute the mean displacement correlation at T=0.
2291 !              Consistent with Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point.
2292 !              but generalized to two different atoms. Note that the complex conjugate is taken on the second direction.
2293                displ2cart(1,idir1,idir2)=(d_at1_dir1_re*d_at2_dir2_re+ &
2294 &               d_at1_dir1_im*d_at2_dir2_im )/(two*phfrq(imode))
2295                displ2cart(2,idir1,idir2)=(d_at1_dir1_im*d_at2_dir2_re- &
2296 &               d_at1_dir1_re*d_at2_dir2_im )/(two*phfrq(imode))
2297              else if(option==2)then
2298 !              Compute the mean square displacement correlation of each atom at T=0, and take mean over iatom1 and iatom2.
2299 !              See Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point.
2300 !              Note that the complex conjugate is taken on the second direction.
2301 !              Also, note the overall negative sign, to make it opposite to the Fan term.
2302                d_at1_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1))))
2303                d_at1_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1))))
2304                d_at2_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1))))
2305                d_at2_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1))))
2306                displ2cart(1,idir1,idir2)=-(d_at1_dir1_re*d_at1_dir2_re+ &
2307 &               d_at1_dir1_im*d_at1_dir2_im+ &
2308 &               d_at2_dir1_re*d_at2_dir2_re+ &
2309 &               d_at2_dir1_im*d_at2_dir2_im )/(four*phfrq(imode))
2310                displ2cart(2,idir1,idir2)=-(d_at1_dir1_im*d_at1_dir2_re- &
2311 &               d_at1_dir1_re*d_at1_dir2_im+ &
2312 &               d_at2_dir1_im*d_at2_dir2_re- &
2313 &               d_at2_dir1_re*d_at2_dir2_im )/(four*phfrq(imode))
2314              end if
2315            end do
2316          end do
2317 !        Switch to reduced coordinates in two steps
2318          tmp_displ2(:,:,:)=zero
2319          do idir1=1,3
2320            do idir2=1,3
2321              tmp_displ2(:,:,idir1)=tmp_displ2(:,:,idir1)+displ2cart(:,:,idir2)*gprimd(idir2,idir1)
2322            end do
2323          end do
2324          displ2red(:,:,:)=zero
2325          do idir1=1,3
2326            do idir2=1,3
2327              displ2red(:,idir1,:)=displ2red(:,idir1,:)+tmp_displ2(:,idir2,:)*gprimd(idir2,idir1)
2328            end do
2329          end do
2330 !        Compute the T=0 shift due to this q point
2331          do idir1=1,3
2332            do idir2=1,3
2333              do ikpt=1,nkpt
2334                do isppol=1,nsppol
2335                  do iband=1,mband
2336                    index=iband+mband*(isppol-1 + nsppol*(ikpt-1))
2337                    e2_re=eig2nkq(1,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2)
2338                    e2_im=eig2nkq(2,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2)
2339                    eigen_corr(index)=eigen_corr(index)+&
2340 &                   e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2)
2341                    eigen_corr_mode(index)=eigen_corr_mode(index)+&
2342 &                   e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2)
2343                  end do  ! band
2344                end do  ! spin
2345              end do  ! kpt
2346            end do  ! dir2
2347          end do  ! dir1
2348        end do  ! atom2
2349      end do  ! atom1
2350    end if
2351 
2352    if (prtvol > 4) then
2353      ! Print the corrections by mode
2354      write(message,'(a,i1)') ' imode= ',imode
2355      call wrtout(ab_out_default,message,'COLL')
2356 
2357      do ikpt=1,nkpt
2358        do isppol=1,nsppol
2359          write(message,'(a,i4,a,i1)')' ikpt= ',ikpt,' ispin= ',isppol
2360          call wrtout(ab_out_default,message,'COLL')
2361 
2362          imin = mband * (isppol-1 + nsppol*(ikpt-1))
2363          do ii=0, (mband-1)/neigs_per_line
2364            write(message, line_format) (eigen_corr_mode(iband+imin), &
2365 &           iband = 1 + ii * neigs_per_line, min(mband, (ii+1)*neigs_per_line))
2366            call wrtout(ab_out_default,message,'COLL')
2367          end do
2368        end do
2369      end do
2370    end if
2371 
2372  end do  ! mode
2373 
2374  if (prtvol > 4) then
2375    if (option==1) then
2376      write(message,'(a)') ' ---- End Fan contribution to eigenvalues renormalization by mode ----'
2377      call wrtout(ab_out_default,message,'COLL')
2378    else if (option==2) then
2379      write(message,'(a)') ' ---- End DDW contribution to eigenvalues renormalization by mode ----'
2380      call wrtout(ab_out_default,message,'COLL')
2381    end if
2382    write(message,'(a,a)')' ================================================================================', ch10
2383    call wrtout(ab_out_default,message,'COLL')
2384  end if
2385 
2386  ABI_DEALLOCATE(eigen_corr_mode)
2387 
2388 !DEBUG
2389 !write(std_out,*)' elph2_fanddw : exit'
2390 !write(std_out,*)' eigen_corr(1)=',eigen_corr(1)
2391 !ENDDEBUG
2392 
2393 end subroutine elph2_fanddw

m_eig2d/fan_free [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 fan_free

FUNCTION

 Deallocates the components of the fan_t structured datatype

INPUTS

  fan2d<fan_t>=The data type to be deallocated.

OUTPUT

  Deallocate the dynamic arrays in the fan_t type.
  (only deallocate)

PARENTS

      eig2tot

CHILDREN

SOURCE

736 subroutine fan_free(fan2d)
737 
738 
739 !This section has been created automatically by the script Abilint (TD).
740 !Do not modify the following lines by hand.
741 #undef ABI_FUNC
742 #define ABI_FUNC 'fan_free'
743 !End of the abilint section
744 
745  implicit none
746 
747 !Arguments ------------------------------------
748 !scalars
749  type(fan_t),intent(inout) :: fan2d
750 ! *************************************************************************
751 DBG_ENTER("COLL")
752 
753 !Deallocate all components of bstruct
754 
755  if (allocated(fan2d%fan2d)) then
756    ABI_FREE(fan2d%fan2d)
757  end if
758 
759  DBG_EXIT("COLL")
760 
761 end subroutine fan_free

m_eig2d/fan_init [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 fan_init

FUNCTION

 This subroutine initializes the fan_t structured datatype

INPUTS

 mbands=maximum number of bands
 nkpt=number of k points
 nsppol=1 for unpolarized, 2 for spin-polarized
 natom=number of atoms
 fan2d(2*mband*nsppol,nkpt,3,natom,3,natom,mband*nsppol)=second-order derivative of the
     eigenenergies wrt phononic displacements

OUTPUT

 fan2d<fan_t>=the fan_t datatype

SIDE EFFECTS

PARENTS

      eig2tot

CHILDREN

SOURCE

411 subroutine fan_init(fan,fan2d,mband,nsppol,nkpt,natom)
412 
413 
414 !This section has been created automatically by the script Abilint (TD).
415 !Do not modify the following lines by hand.
416 #undef ABI_FUNC
417 #define ABI_FUNC 'fan_init'
418 !End of the abilint section
419 
420  implicit none
421 
422 !This section has been created automatically by the script Abilint (TD).
423 !Do not modify the following lines by hand.
424 #undef ABI_FUNC
425 #define ABI_FUNC 'fan_init'
426 !End of the abilint section
427 
428 !Arguments ------------------------------------
429 !scalars
430  integer,intent(in) ::mband,nsppol,nkpt,natom
431  type(fan_t),intent(out) :: fan2d
432 !arrays
433  real(dp), intent(in) :: fan(2*mband*nsppol,nkpt,3,natom,3,natom,mband)
434 ! *************************************************************************
435 
436  fan2d%mband = mband
437  fan2d%nsppol = nsppol
438  fan2d%nkpt = nkpt
439  fan2d%natom = natom
440 
441  ABI_MALLOC(fan2d%fan2d,(2*mband*nsppol,nkpt,3,natom,3,natom,mband))
442  fan2d%fan2d=fan
443 
444 end subroutine fan_init

m_eig2d/fan_ncwrite [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 fan_ncwrite

FUNCTION

  Writes the content of a fan_t object to a NETCDF file
  according to the ETSF-IO specifications.

INPUTS

  ncid =NC file handle

OUTPUT

PARENTS

      eig2tot

CHILDREN

SOURCE

533 subroutine fan_ncwrite(fan2d,iqpt,wtq,ncid)
534 
535 
536 !This section has been created automatically by the script Abilint (TD).
537 !Do not modify the following lines by hand.
538 #undef ABI_FUNC
539 #define ABI_FUNC 'fan_ncwrite'
540 !End of the abilint section
541 
542  implicit none
543 
544 !Arguments ------------------------------------
545 !scalars
546  integer,intent(in) ::ncid
547  real(dp),intent(in) :: iqpt(3),wtq
548  type(fan_t),intent(in) :: fan2d
549 
550 !Local variables-------------------------------
551 #ifdef HAVE_NETCDF
552  integer :: ncerr
553  integer :: cplex,cart_dir,one_dim
554 
555 ! *************************************************************************
556 
557  ! ==============================================
558  ! === Write the dimensions specified by ETSF ===
559  ! ==============================================
560  one_dim=1; cplex=2; cart_dir=3
561 
562  ncerr = nctk_def_dims(ncid, [&
563    nctkdim_t('max_number_of_states',fan2d%mband),&
564    nctkdim_t('number_of_spins',fan2d%nsppol),&
565    nctkdim_t('number_of_kpoints',fan2d%nkpt),&
566    nctkdim_t('number_of_atoms',fan2d%natom),&
567    nctkdim_t('3_number_of_atoms',3*fan2d%natom),&     ! TODO: not sure that variables can start with digits
568    nctkdim_t('number_of_cartesian_directions',cart_dir),&
569    nctkdim_t('current_one_dim',one_dim),&
570    nctkdim_t('cplex',cplex),&
571    nctkdim_t('product_mband_nsppol',fan2d%mband*fan2d%nsppol),&
572    nctkdim_t('product_mband_nsppol2',fan2d%mband*fan2d%nsppol*2) &
573  ], defmode=.True.)
574  NCF_CHECK(ncerr)
575 
576  ncerr = nctk_def_arrays(ncid, [&
577    nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'),&
578    nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'),&
579    nctkarr_t('second_derivative_eigenenergies_actif', "dp", &
580 &'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions, &
581 &number_of_atoms, number_of_cartesian_directions, number_of_atoms, max_number_of_states')])
582  NCF_CHECK(ncerr)
583 
584 ! Write data
585  NCF_CHECK(nctk_set_datamode(ncid))
586  NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt))
587  NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq))
588  NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies_actif'), fan2d%fan2d))
589 
590 #else
591  MSG_ERROR("netcdf support is not activated. ")
592 #endif
593 
594 contains
595  integer function vid(vname)
596 
597 !This section has been created automatically by the script Abilint (TD).
598 !Do not modify the following lines by hand.
599 #undef ABI_FUNC
600 #define ABI_FUNC 'vid'
601 !End of the abilint section
602 
603    character(len=*),intent(in) :: vname
604    vid = nctk_idname(ncid, vname)
605  end function vid
606 
607 end subroutine fan_ncwrite

m_eig2d/fan_t [ Types ]

[ Top ] [ m_eig2d ] [ Types ]

NAME

 fan_t

FUNCTION

 It contains information about the second-order derivative of the
 eigenenergies wrt atomic displacement

SOURCE

117  type,public :: fan_t
118 
119 ! WARNING : if you modify this datatype, please check whether there might be
120 ! creation/destruction/copy routines,
121 ! declared in another part of ABINIT, that might need to take into account your
122 ! modification.
123 
124   integer :: mband                 ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
125   integer :: nsppol                ! number of spin-polarization
126   integer :: nkpt                  ! number of k points
127   integer :: natom                 ! number of atoms
128 
129   real(dp),allocatable :: fan2d(:,:,:,:,:,:,:)
130   ! fan2d(2*mband*nsppol,nkpt,3,natom,3,natom,mband)
131   ! Second-order derivative of the eigenergies (real,im) at each
132   ! ispin,iband(real,im),k-point,dir1,dir2,natom1,natom2,jband
133 
134  end type fan_t

m_eig2d/gkk_free [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 gkk_free

FUNCTION

 Deallocates the components of the gkk_t structured datatype

INPUTS

  gkk2d<gkk_t>=The data type to be deallocated.

OUTPUT

  Deallocate the dynamic arrays in the gkk_t type.
  (only deallocate)

PARENTS

      dfpt_looppert,eig2tot,m_gkk

CHILDREN

SOURCE

787 subroutine gkk_free(gkk2d)
788 
789 
790 !This section has been created automatically by the script Abilint (TD).
791 !Do not modify the following lines by hand.
792 #undef ABI_FUNC
793 #define ABI_FUNC 'gkk_free'
794 !End of the abilint section
795 
796  implicit none
797 
798 !Arguments ------------------------------------
799 !scalars
800  type(gkk_t),intent(inout) :: gkk2d
801 ! *************************************************************************
802 DBG_ENTER("COLL")
803 
804 !Deallocate all components of bstruct
805 
806  if (allocated(gkk2d%gkk2d)) then
807    ABI_FREE(gkk2d%gkk2d)
808  end if
809 
810  DBG_EXIT("COLL")
811 
812 end subroutine gkk_free

m_eig2d/gkk_init [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 gkk_init

FUNCTION

 This subroutine initializes the gkk_t structured datatype

INPUTS

 mbands=maximum number of bands
 nkpt=number of k points
 nsppol=1 for unpolarized, 2 for spin-polarized
 natom=number of atoms
 gkk2d(2*mband*nsppol,nkpt,3,natom,mband*nsppol)=second-order derivative of the
     eigenenergies wrt phononic displacements

OUTPUT

 gkk2d<gkk_t>=the gkk_t datatype

SIDE EFFECTS

PARENTS

      dfpt_looppert,eig2tot,m_gkk

CHILDREN

SOURCE

474 subroutine gkk_init(gkk,gkk2d,mband,nsppol,nkpt,natom,ncart)
475 
476 
477 !This section has been created automatically by the script Abilint (TD).
478 !Do not modify the following lines by hand.
479 #undef ABI_FUNC
480 #define ABI_FUNC 'gkk_init'
481 !End of the abilint section
482 
483  implicit none
484 
485 !This section has been created automatically by the script Abilint (TD).
486 !Do not modify the following lines by hand.
487 #undef ABI_FUNC
488 #define ABI_FUNC 'gkk_init'
489 !End of the abilint section
490 
491 !Arguments ------------------------------------
492 !scalars
493  integer,intent(in) ::mband,nsppol,nkpt,natom,ncart
494  type(gkk_t),intent(out) :: gkk2d
495 !arrays
496  real(dp), intent(in) :: gkk(2*mband*nsppol,nkpt,ncart,natom,mband)
497 ! *************************************************************************
498 
499  gkk2d%mband = mband
500  gkk2d%nsppol = nsppol
501  gkk2d%nkpt = nkpt
502  gkk2d%natom = natom
503  gkk2d%ncart = ncart
504 
505  ABI_MALLOC(gkk2d%gkk2d,(2*mband*nsppol,nkpt,ncart,natom,mband))
506  gkk2d%gkk2d=gkk
507 
508 end subroutine gkk_init

m_eig2d/gkk_ncwrite [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 gkk_ncwrite

FUNCTION

  Writes the content of a gkk_t object to a NETCDF file
  according to the ETSF-IO specifications.

INPUTS

  ncid =NC file handle

OUTPUT

PARENTS

      dfpt_looppert,eig2tot,m_gkk

CHILDREN

SOURCE

632 subroutine gkk_ncwrite(gkk2d,iqpt,wtq,ncid)
633 
634 
635 !This section has been created automatically by the script Abilint (TD).
636 !Do not modify the following lines by hand.
637 #undef ABI_FUNC
638 #define ABI_FUNC 'gkk_ncwrite'
639 !End of the abilint section
640 
641  implicit none
642 
643 !Arguments ------------------------------------
644 !scalars
645  integer,intent(in) ::ncid
646  real(dp), intent(in) :: iqpt(3),wtq
647  type(gkk_t),intent(in) :: gkk2d
648 
649 !Local variables-------------------------------
650 #ifdef HAVE_NETCDF
651  integer :: cplex,one_dim,ncerr,vid_
652 
653 ! *************************************************************************
654 
655  ! ==============================================
656  ! === Write the dimensions specified by ETSF ===
657  ! ==============================================
658  one_dim=1; cplex=2
659 
660  ncerr = nctk_def_dims(ncid, [ &
661 &   nctkdim_t('max_number_of_states', gkk2d%mband), &
662 &   nctkdim_t('number_of_spins', gkk2d%nsppol), &
663 &   nctkdim_t('number_of_kpoints', gkk2d%nkpt), &
664 &   nctkdim_t('number_of_atoms_for_gkk', gkk2d%natom), &
665 &   nctkdim_t('3_number_of_atoms', 3*gkk2d%natom), &
666 &   nctkdim_t('number_of_cartesian_directions_for_gkk', gkk2d%ncart), &
667 &   nctkdim_t('current_one_dim', one_dim), &
668 &   nctkdim_t('cplex', cplex), &
669 &   nctkdim_t('product_mband_nsppol', gkk2d%mband*gkk2d%nsppol), &
670 &   nctkdim_t('product_mband_nsppol2', gkk2d%mband*gkk2d%nsppol*2) &
671 & ], defmode=.True.)
672  NCF_CHECK(ncerr)
673 
674 !arrays
675  ncerr = nctk_def_arrays(ncid, [&
676 &   nctkarr_t('current_q_point', "dp", "number_of_cartesian_directions"), &
677 &   nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'), &
678 &   nctkarr_t('second_derivative_eigenenergies_actif', "dp", &
679 &     'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions_for_gkk,'// &
680 &     'number_of_atoms_for_gkk, max_number_of_states') &
681 & ])
682  NCF_CHECK(ncerr)
683 
684  NCF_CHECK(nctk_set_datamode(ncid))
685  vid_=vid('current_q_point')
686  NCF_CHECK(nf90_put_var(ncid, vid_, iqpt))
687  vid_=vid('current_q_point_weight')
688  NCF_CHECK(nf90_put_var(ncid, vid_, wtq))
689  vid_=vid('second_derivative_eigenenergies_actif')
690  NCF_CHECK(nf90_put_var(ncid, vid_, gkk2d%gkk2d))
691 
692 #else
693  MSG_ERROR("netcdf support is not activated. ")
694 #endif
695 
696 contains
697  integer function vid(vname)
698 
699 
700 !This section has been created automatically by the script Abilint (TD).
701 !Do not modify the following lines by hand.
702 #undef ABI_FUNC
703 #define ABI_FUNC 'vid'
704 !End of the abilint section
705 
706    character(len=*),intent(in) :: vname
707    vid = nctk_idname(ncid, vname)
708  end function vid
709 
710 end subroutine gkk_ncwrite

m_eig2d/gkk_t [ Types ]

[ Top ] [ m_eig2d ] [ Types ]

NAME

 gkk_t

FUNCTION

 It contains information about the second-order derivative of the
 eigenenergies wrt atomic displacement

SOURCE

147  type,public :: gkk_t
148 
149 ! WARNING : if you modify this datatype, please check whether there might be
150 ! creation/destruction/copy routines,
151 ! declared in another part of ABINIT, that might need to take into account your
152 ! modification.
153 
154   integer :: mband                 ! Max number of bands i.e MAXVAL(nband) (to dimension arrays)
155   integer :: nsppol                ! number of spin-polarization
156   integer :: nkpt                  ! number of k points
157   integer :: natom                 ! number of atoms
158   integer :: ncart                 ! number of cartesian directions
159 
160   real(dp),allocatable :: gkk2d(:,:,:,:,:)
161   ! gkk2d(2*mband*nsppol,nkpt,ncart,natom,mband)
162   ! Second-order derivative of the eigenergies (real,im) at each
163   ! ispin,iband(real,im),k-point,dir1,natom1,jband
164 
165  end type gkk_t

m_eig2d/outbsd [ Functions ]

[ Top ] [ m_eig2d ] [ Functions ]

NAME

 outbsd

FUNCTION

 output bsd file for one perturbation (used for elphon calculations in anaddb)

INPUTS

  bdeigrf=number of bands for which the derivatives of the eigenvalues have been computed
  dtset = dataset variable for run flags
  eig2nkq= second ordre eigenvalue (or electron lifetime) that must be printed out
  mpert= maximum number of perturbations
  nkpt_rbz= number of k-points for perturbation
  unitout= writting unit of file

 OUTPUTS
  to file

PARENTS

      dfpt_looppert,eig2tot

CHILDREN

SOURCE

1982 subroutine outbsd(bdeigrf,dtset,eig2nkq,mpert,nkpt_rbz,unitout)
1983 
1984 
1985 !This section has been created automatically by the script Abilint (TD).
1986 !Do not modify the following lines by hand.
1987 #undef ABI_FUNC
1988 #define ABI_FUNC 'outbsd'
1989 !End of the abilint section
1990 
1991  implicit none
1992 
1993 !Arguments ------------------------------------
1994 !scalars
1995  integer,intent(in) :: bdeigrf,mpert,nkpt_rbz,unitout
1996  type(dataset_type),intent(in) :: dtset
1997 !arrays
1998  real(dp),intent(in) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt_rbz,3,mpert,3,mpert)
1999 
2000 !Local variables -------------------------
2001 !scalars
2002  integer :: bandtot_index,iband,idir1,idir2,ikpt,ipert1,ipert2,isppol
2003 
2004 ! *********************************************************************
2005 
2006 !DEBUG
2007 !write(std_out,*)' outbsd : enter'
2008 !write(std_out,*)' eig2nkq(1,1,1,1,1,1)=',eig2nkq(1,1,1,1,1,1,1)
2009 !ENDDEBUG
2010 
2011 !output information in this file
2012  write(unitout,*)
2013  write(unitout,'(a,i8)') ' 2nd eigenvalue derivatives   - # elements :', 9*dtset%natom**2
2014  write(unitout,'(a,3es16.8,a)') ' qpt', dtset%qptn(:), ' 1.0'
2015 
2016 !output RF eigenvalues
2017 
2018  do ikpt=1,nkpt_rbz
2019 !  bandtot_index differs from zero only in the spin-polarized case
2020    bandtot_index=0
2021    write (unitout,'(a,3es16.8)') ' K-point:', dtset%kptns(:,ikpt)
2022    do isppol=1,dtset%nsppol
2023      do iband=1,bdeigrf
2024        write (unitout,'(a,i5)') ' Band:', iband+bandtot_index
2025 !      write (unitout,*) 'ipert1     ','idir1     ','ipert2     ','idir2    ','Real    ','Im    '
2026        do ipert2=1,mpert
2027          do idir2=1,3
2028            do ipert1=1,mpert
2029              do idir1=1,3
2030                write (unitout,'(4i4,2d22.14)') idir1,ipert1,idir2,ipert2,&
2031 &               eig2nkq(1,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2),&
2032 &               eig2nkq(2,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2)
2033              end do !idir2
2034            end do !ipert2
2035          end do !idir1
2036        end do !ipert1
2037      end do !iband
2038      bandtot_index = bandtot_index + dtset%mband
2039    end do !isppol
2040  end do !ikpt
2041 
2042 !close bsd file
2043  close (unitout)
2044 
2045 end subroutine outbsd