## 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


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
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
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
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
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 (C) 2014-2018 ABINIT group (SP, PB, XG)
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


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
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
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
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
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