TABLE OF CONTENTS
- ABINIT/eig2stern
- ABINIT/m_eig2d
- m-eig2d/smeared_delta
- m_eig2d/eig2tot
- m_eig2d/eigr2d_free
- m_eig2d/eigr2d_init
- m_eig2d/eigr2d_ncwrite
- m_eig2d/eigr2d_t
- m_eig2d/elph2_fanddw
- m_eig2d/fan_free
- m_eig2d/fan_init
- m_eig2d/fan_ncwrite
- m_eig2d/fan_t
- m_eig2d/gkk_free
- m_eig2d/gkk_init
- m_eig2d/gkk_ncwrite
- m_eig2d/gkk_t
- m_eig2d/outbsd
ABINIT/eig2stern [ 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.
SOURCE
731 subroutine eig2stern(occ,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0,eigenq,& 732 & eigen1,eig2nkq,elph2_imagden,esmear,gh0c1_pert,gh1c_pert,ieig2rf,istwfk_pert,& 733 & mband,mk1mem,mpert,npert,mpi_enreg,mpw1,nkpt_rbz,npwar1,nspinor,nsppol,smdelta,& 734 & dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 735 736 !Arguments ------------------------------------ 737 !scalars 738 integer,intent(in) :: bdeigrf,dim_eig2nkq,dim_eig2rf,ieig2rf,mband,mk1mem,mpert,mpw1,nkpt_rbz 739 integer,intent(in) :: npert,nspinor,nsppol,smdelta 740 real(dp),intent(in) :: elph2_imagden,esmear 741 type(MPI_type),intent(inout) :: mpi_enreg 742 !arrays 743 integer,intent(in) :: clflg(3,mpert) 744 integer,intent(in) :: istwfk_pert(nkpt_rbz,3,mpert) 745 integer,intent(in) :: npwar1(nkpt_rbz,mpert) 746 real(dp),intent(in) :: cg1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 747 real(dp),intent(in) :: gh0c1_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 748 real(dp),intent(in) :: gh1c_pert(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf,3,mpert) 749 real(dp),intent(inout) :: eigen0(nkpt_rbz*mband*nsppol) 750 real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) 751 real(dp),intent(inout) :: eigenq(nkpt_rbz*mband*nsppol) 752 real(dp),intent(out) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq) 753 real(dp),intent(out),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert) 754 real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:) 755 real(dp), intent(in) :: occ(mband*nkpt_rbz*nsppol) 756 type(dataset_type), intent(in) :: dtset 757 type(hdr_type),intent(in),optional :: hdr_fine,hdr0 758 759 !Local variables------------------------------- 760 !tolerance for non degenerated levels 761 !scalars 762 integer :: band2tot_index,band_index,bandtot_index,iband,icg2,idir1,idir2 763 integer :: ikpt,ipert1,ipert2,isppol,istwf_k,jband,npw1_k,nkpt_sub,ikpt2 764 !integer :: ipw 765 integer :: master,me,spaceworld,ierr 766 integer :: mband_mem 767 !real(dp),parameter :: etol=1.0d-3 768 real(dp),parameter :: etol=1.0d-6 769 !real(dp),parameter :: etol=zero 770 real(dp) :: ar,ai,deltae,den,dot2i,dot2r,dot3i,dot3r,doti,dotr,eig1_i1,eig1_i2 771 real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av 772 real(dp) :: wgt_int 773 real(dp) :: eig2_diar,eigbrd_i,eigbrd_r 774 character(len=500) :: message 775 character(len=500) :: msg 776 !DBSP 777 ! character(len=300000) :: message2 778 !END 779 logical :: test_do_band 780 !arrays 781 integer, allocatable :: nband_rbz(:),icg2_rbz(:,:) 782 integer,pointer :: kpt_fine_sub(:) 783 real(dp) :: tsec(2) 784 real(dp),allocatable :: cwavef(:,:),cwavef2(:,:),center(:),eigen0tmp(:),eigenqtmp(:) 785 real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol) 786 real(dp),allocatable :: gh(:,:),gh1(:,:),ghc(:,:) 787 real(dp),allocatable :: smdfun(:,:) 788 real(dp),pointer :: wgt_sub(:) 789 790 ! ********************************************************************* 791 792 !Init parallelism 793 master =0 794 spaceworld=mpi_enreg%comm_cell 795 me=mpi_enreg%me_kpt 796 797 !Init interpolation method 798 if(present(eigenq_fine))then 799 ABI_MALLOC(center,(3)) 800 end if 801 802 call timab(148,1,tsec) 803 804 if(nsppol==2)then 805 message = 'nsppol=2 is still under development. Be careful when using it ...' 806 ABI_COMMENT(message) 807 end if 808 809 band2tot_index =0 810 bandtot_index=0 811 band_index=0 812 813 !Add scissor shift to eigenenergies 814 if (dtset%dfpt_sciss > tol6 ) then 815 write(msg,'(a,f7.3,2a)')& 816 & ' A scissor operator of ',dtset%dfpt_sciss*Ha_eV,' [eV] has been applied to the eigenenergies',ch10 817 call wrtout(std_out,msg,'COLL') 818 call wrtout(ab_out,msg,'COLL') 819 ABI_MALLOC(eigen0tmp,(nkpt_rbz*mband*nsppol)) 820 ABI_MALLOC(eigenqtmp,(nkpt_rbz*mband*nsppol)) 821 eigen0tmp = eigen0(:) 822 eigenqtmp = eigenq(:) 823 eigen0 = zero 824 eigenq = zero 825 end if 826 827 if(ieig2rf > 0) then 828 eig2nkq(:,:,:,:,:,:,:) = zero 829 end if 830 if(present(eigbrd))then 831 eigbrd(:,:,:,:,:,:,:) = zero 832 end if 833 834 if(xmpi_paral==1) then 835 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol)) 836 ABI_MALLOC(nband_rbz,(nkpt_rbz*nsppol)) 837 if (allocated(mpi_enreg%my_kpttab)) then 838 ABI_FREE(mpi_enreg%my_kpttab) 839 end if 840 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz)) 841 ! Assume the number of bands is the same for all k points. 842 nband_rbz(:)=mband 843 call distrb2(mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg) 844 end if 845 846 icg2=0 847 ipert1=1 ! Suppose that the situation is the same for all perturbations 848 ABI_MALLOC(icg2_rbz,(nkpt_rbz,nsppol)) 849 do isppol=1,nsppol 850 do ikpt=1,nkpt_rbz 851 icg2_rbz(ikpt,isppol)=icg2 852 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) cycle 853 icg2 = icg2 + npwar1(ikpt,ipert1)*nspinor*mband 854 end do 855 end do 856 857 do isppol=1,nsppol 858 do ikpt =1,nkpt_rbz 859 860 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then 861 band2tot_index = band2tot_index + 2*mband**2 862 bandtot_index = bandtot_index + mband 863 cycle 864 end if 865 866 if(present(eigenq_fine))then 867 write(std_out,*) 'Start of the energy denominator interpolation method.' 868 nkpt_sub = 0 869 ! center is the k+q point around which we will average the kpt_fine 870 center = hdr0%kptns(:,ikpt)+ dtset%qptn(:) 871 872 call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub) 873 write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid & 874 & around the k+Q point ',center,' is:',nkpt_sub 875 write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub) 876 end if 877 878 ! Add scissor shift to eigenenergies 879 if (dtset%dfpt_sciss > tol6 ) then 880 do iband=1,mband 881 if (occ(iband+bandtot_index) < tol6) then 882 eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) + dtset%dfpt_sciss 883 eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) + dtset%dfpt_sciss 884 else 885 eigen0(iband+bandtot_index) = eigen0tmp(iband+bandtot_index) 886 eigenq(iband+bandtot_index) = eigenqtmp(iband+bandtot_index) 887 end if 888 end do 889 end if 890 891 892 if(smdelta >0) then !broadening 893 if(.not.allocated(smdfun)) then 894 ABI_MALLOC(smdfun,(mband,mband)) 895 end if 896 smdfun(:,:) = zero 897 do iband=1,mband 898 eigen(iband) = eigen0(iband+bandtot_index) 899 eigen_prime(iband) =eigenq(iband+bandtot_index) 900 end do 901 if(esmear>tol6) then 902 call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun) 903 end if 904 end if 905 icg2=icg2_rbz(ikpt,isppol) 906 907 ipert1=1 ! Suppose all perturbations lead to the same number of planewaves 908 npw1_k = npwar1(ikpt,ipert1) 909 ABI_MALLOC(cwavef,(2,npw1_k*nspinor)) 910 ABI_MALLOC(cwavef2,(2,npw1_k*nspinor)) 911 ABI_MALLOC(gh,(2,npw1_k*nspinor)) 912 ABI_MALLOC(gh1,(2,npw1_k*nspinor)) 913 ABI_MALLOC(ghc,(2,npw1_k*nspinor)) 914 915 do iband=1,bdeigrf 916 917 ! If the k point and band belong to me, compute the contribution 918 test_do_band=.true. 919 if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false. 920 921 if(test_do_band)then 922 923 do ipert1=1,npert 924 925 do idir1=1,3 926 if(clflg(idir1,ipert1)==0)cycle 927 istwf_k = istwfk_pert(ikpt,idir1,ipert1) 928 929 do ipert2=1,npert 930 do idir2=1,3 931 if(clflg(idir2,ipert2)==0)cycle 932 933 eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero 934 935 do jband=1,mband 936 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 937 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 938 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 939 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 940 ! If no interpolation, fallback on to the previous 941 ! implementation 942 if(.not. present(eigenq_fine))then 943 deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index) 944 end if 945 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 946 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 947 948 ! Sum over all active space to retrieve the diagonal gauge 949 if(ieig2rf == 1 .or. ieig2rf ==2 ) then 950 ! if(abs(deltae)>etol) then ! This is commented because 951 ! there is no problem with divergencies with elph2_imag != 0 952 if( present(eigenq_fine))then 953 den_av = zero 954 wgt_int = zero 955 do ikpt2=1,nkpt_sub 956 deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)& 957 & -eigen0(iband+bandtot_index) 958 den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2) 959 wgt_int = wgt_int+wgt_sub(ikpt2) 960 end do 961 den = den_av/wgt_int 962 else 963 if(abs(elph2_imagden) < etol) then 964 if(abs(deltae)>etol) then 965 den=-one/(deltae**2+elph2_imagden**2) 966 else 967 den= zero 968 end if 969 else 970 den=-one/(deltae**2+elph2_imagden**2) 971 end if 972 end if 973 974 ! The following should be the most general implementation of the presence of elph2_imagden 975 ! eig2_diar=eig2_diar+(ar*deltae+ai*elph2_imagden)*den 976 ! eig2_diai=eig2_diai+(ai*deltae-ar*elph2_imagden)*den 977 ! This gives back the implementation without elph2_imagden 978 ! eig2_diar=eig2_diar+ar*deltae*den 979 ! eig2_diai=eig2_diai+ai*deltae*den 980 ! This is what Samuel had implemented 981 ! eig2_diar=eig2_diar+ar*deltae*den 982 ! eig2_diai=eig2_diai+ai*elph2_imagden*den 983 ! Other possibility : throw away the broadening part, that is actually treated separately. 984 if( present(eigenq_fine))then 985 eig2_diar=eig2_diar+ar*den 986 eig2_diai=eig2_diai+ai*den 987 else 988 eig2_diar=eig2_diar+ar*deltae*den 989 eig2_diai=eig2_diai+ai*deltae*den 990 !DBSP 991 ! if (iband+band_index==2 .and. ikpt==1 .and. idir1==1 .and. ipert1==1 .and. idir2==1 .and. ipert2==1) then 992 ! write(message2,*) 'eig2_diar1=',eig2_diar,' ar=',ar,' deltae=',deltae,' den=',den 993 ! call wrtout(std_out,message2,'PERS') 994 ! endif 995 !END 996 997 end if 998 end if ! ieig2rf==1 or 2 999 1000 if(present(eigbrd))then 1001 if(smdelta >0) then !broadening 1002 eigbrd_r = eigbrd_r + ar*smdfun(iband,jband) 1003 eigbrd_i = eigbrd_i + ai*smdfun(iband,jband) 1004 end if 1005 end if 1006 1007 end do !jband 1008 1009 ! Add the contribution of non-active bands, if DFPT calculation (= Sternheimer) 1010 if(ieig2rf == 1 .or. ieig2rf ==3 .or. ieig2rf ==4 .or. ieig2rf==5 ) then 1011 ! if(ieig2rf == 1 ) then 1012 1013 dotr=zero ; doti=zero 1014 dot2r=zero ; dot2i=zero 1015 dot3r=zero ; dot3i=zero 1016 1017 1018 cwavef(:,:) = cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2) 1019 cwavef2(:,:)= cg1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 1020 gh1(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 1021 gh(:,:) = gh1c_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir2,ipert2) 1022 ghc(:,:) = gh0c1_pert(:,1+(iband-1)*npw1_k*nspinor+icg2:iband*npw1_k*nspinor+icg2,idir1,ipert1) 1023 1024 ! 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> 1025 ! They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space. 1026 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,cwavef,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1027 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,gh,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1028 1029 ! This dotprod corresponds to : <Psi(1)nkq|H(0)k+q- E(0)nk|Psi(1)nkq> 1030 ! It is calculated using wavefunctions that are orthogonal to the active space. 1031 ! Should work for metals. (But adiabatic approximation is bad in this case...) 1032 call dotprod_g(dot3r,dot3i,istwf_k,npw1_k*nspinor,2,cwavef,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1033 1034 eig2_diar= eig2_diar + dotr + dot2r + dot3r 1035 eig2_diai= eig2_diai + doti + dot2i + dot3i 1036 1037 end if 1038 1039 ! Store the contribution 1040 if(ieig2rf > 0) then 1041 eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diar 1042 eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eig2_diai 1043 end if 1044 1045 if(present(eigbrd))then 1046 if(smdelta >0) then !broadening 1047 eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r 1048 eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i 1049 end if 1050 end if 1051 1052 end do !idir2 1053 end do !ipert2 1054 end do !idir1 1055 end do !ipert1 1056 1057 end if ! Selection of processor 1058 1059 end do !iband 1060 1061 ABI_FREE(cwavef) 1062 ABI_FREE(cwavef2) 1063 ABI_FREE(gh) 1064 ABI_FREE(gh1) 1065 ABI_FREE(ghc) 1066 band2tot_index = band2tot_index + 2*mband**2 1067 bandtot_index = bandtot_index + mband 1068 1069 if(present(eigenq_fine))then 1070 ABI_FREE(kpt_fine_sub) ! Deallocate the variable 1071 ABI_FREE(wgt_sub) 1072 end if 1073 1074 end do !ikpt 1075 band_index = band_index + mband 1076 end do !isppol 1077 1078 !Accumulate eig2nkq and/or eigbrd 1079 if(xmpi_paral==1) then 1080 if(ieig2rf == 1 .or. ieig2rf == 2) then 1081 call xmpi_sum(eig2nkq,spaceworld,ierr) 1082 if (dtset%dfpt_sciss > tol6 ) then 1083 call xmpi_sum(eigen0,spaceworld,ierr) 1084 call xmpi_sum(eigenq,spaceworld,ierr) 1085 end if 1086 end if 1087 if(present(eigbrd) .and. (ieig2rf == 1 .or. ieig2rf == 2))then 1088 if(smdelta >0) then 1089 call xmpi_sum(eigbrd,spaceworld,ierr) 1090 end if 1091 end if 1092 ABI_FREE(nband_rbz) 1093 ABI_FREE(mpi_enreg%proc_distrb) 1094 ABI_FREE(mpi_enreg%my_kpttab) 1095 end if 1096 1097 if(ieig2rf==1 .or. ieig2rf==2 ) then 1098 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.' 1099 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 1100 do idir1=1,3 1101 do idir2=1,3 1102 ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 1103 ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 1104 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 1105 end do ! idir2 1106 end do ! idir1 1107 end if 1108 1109 if(present(eigbrd))then 1110 if(smdelta >0) then !broadening 1111 write(ab_out,'(a)')' ' 1112 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.' 1113 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 1114 do idir1=1,3 1115 do idir2=1,3 1116 ar=eigbrd(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 1117 ai=eigbrd(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 1118 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 1119 end do 1120 end do !nband 1121 end if 1122 end if 1123 1124 if(allocated(smdfun)) then 1125 ABI_FREE(smdfun) 1126 end if 1127 ABI_FREE(icg2_rbz) 1128 if(present(eigenq_fine))then 1129 ABI_FREE(center) 1130 end if 1131 if (dtset%dfpt_sciss > tol6 ) then 1132 ABI_FREE(eigen0tmp) 1133 ABI_FREE(eigenqtmp) 1134 end if 1135 1136 call timab(148,2,tsec) 1137 1138 end subroutine eig2stern
ABINIT/m_eig2d [ 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-2024 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 .
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 MODULE m_eig2d 25 26 use defs_basis 27 use m_errors 28 use m_abicore 29 use m_nctk 30 #ifdef HAVE_NETCDF 31 use netcdf 32 #endif 33 use m_xmpi 34 use m_ebands 35 use m_cgtools 36 use m_hdr 37 use m_dtset 38 use m_dtfil 39 use m_ddb_hdr 40 use m_ddb 41 42 use defs_datatypes, only : pseudopotential_type, ebands_t 43 use defs_abitypes, only : MPI_type 44 use m_time, only : timab 45 use m_fstrings, only : strcat 46 use m_crystal, only : crystal_init, crystal_t 47 use m_pawtab, only : pawtab_type 48 use m_double_grid,only : kptfine_av 49 use m_mpinfo, only : distrb2, proc_distrb_cycle 50 51 implicit none 52 53 private 54 55 public :: eigr2d_init ! Main creation method of EIG2D.nc files. 56 public :: eigr2d_ncwrite ! Dump the object into NETCDF file. 57 public :: eigr2d_free ! Destruction method. 58 public :: fan_init ! Main creation method of Fan.nc files. 59 public :: fan_ncwrite ! Dump the object into NETCDF file. 60 public :: fan_free ! Destruction method. 61 public :: gkk_init ! Main creation method of GKK.nc files. 62 public :: gkk_ncwrite ! Dump the object into NETCDF file. 63 public :: gkk_free ! Destruction method. 64 65 public :: eig2tot ! This routine calculates the second-order eigenvalues. 66 public :: outbsd ! output bsd file for one perturbation (used for elphon calculations in anaddb) 67 public :: eig2stern 68 public :: elph2_fanddw ! Calculates the zero-point motion corrections
m-eig2d/smeared_delta [ 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})
SOURCE
1846 subroutine smeared_delta(eigen0,eigenq,esmear,mband,smdelta,smdfunc) 1847 1848 !Arguments ------------------------------------ 1849 !scalars 1850 integer,intent(in) :: mband,smdelta 1851 !arrays 1852 real(dp),intent(in) :: eigen0(mband),eigenq(mband),esmear 1853 real(dp),intent(out) :: smdfunc(mband,mband) 1854 1855 !Local variables------------------------------- 1856 !tolerance for non degenerated levels 1857 !scalars 1858 integer :: ii,jj 1859 real(dp) :: aa,dsqrpi,gauss,xx 1860 character(len=500) :: message 1861 1862 ! ********************************************************************* 1863 1864 1865 !--------------------------------------------------------- 1866 !Ordinary (unique) smearing function 1867 !--------------------------------------------------------- 1868 1869 if(smdelta==1)then 1870 1871 ! Fermi-Dirac 1872 do ii=1,mband 1873 do jj= 1,mband 1874 xx= ( eigen0(ii) - eigenq(jj) )/esmear 1875 smdfunc(ii,jj)=0.25_dp/esmear/(cosh(xx/2.0_dp))**2 1876 end do 1877 end do 1878 1879 else if(smdelta==2 .or. smdelta==3)then 1880 1881 ! Cold smearing of Marzari, two values of the "a" parameter being possible 1882 ! first value gives minimization of the bump 1883 if(smdelta==2)aa=-.5634 1884 ! second value gives monotonic occupation function 1885 if(smdelta==3)aa=-.8165 1886 1887 dsqrpi=1.0_dp/sqrt(pi) 1888 do ii=1,mband 1889 do jj=1,mband 1890 xx= ( eigen0(ii) - eigenq(jj) ) / esmear 1891 gauss=dsqrpi*exp(-xx**2)/esmear 1892 smdfunc(ii,jj)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx))) 1893 end do 1894 end do 1895 1896 else if(smdelta==4)then 1897 1898 ! First order Hermite-Gaussian of Paxton and Methfessel 1899 dsqrpi=1.0_dp/sqrt(pi) 1900 do ii=1,mband 1901 do jj=1,mband 1902 xx= ( eigen0(ii) - eigenq (jj) ) / esmear 1903 smdfunc(ii,jj)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)/esmear 1904 end do 1905 end do 1906 1907 else if(smdelta==5)then 1908 1909 ! Gaussian smearing 1910 dsqrpi=1.0_dp/sqrt(pi) 1911 do ii=1,mband 1912 do jj=1,mband 1913 xx= ( eigen0(ii) - eigenq (jj) ) / esmear 1914 smdfunc(ii,jj)=dsqrpi*exp(-xx**2)/esmear 1915 end do 1916 end do 1917 1918 else 1919 write(message, '(a,i0,a)' )' Smdelta= ',smdelta,' is not allowed in smdfunc' 1920 ABI_BUG(message) 1921 end if 1922 1923 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
the contribution of each perturbations pair
to the eigenstate broadening (inverse lifetime)
computed statically (without phonon frequencies).
SOURCE
1205 subroutine eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0,eigenq,eigen1,eig2nkq,& 1206 & elph2_imagden,esmear,ieig2rf,mband,mpert,npert,mpi_enreg,doccde,& 1207 & nkpt_rbz,nsppol,smdelta,rprimd,dtset,occ_rbz,hdr0,eigbrd,eigenq_fine,hdr_fine) 1208 1209 !Arguments ------------------------------------ 1210 !scalars 1211 integer,intent(in) :: bdeigrf,dim_eig2nkq,ieig2rf,mband,mpert,natom,nkpt_rbz 1212 integer,intent(in) :: npert,nsppol,smdelta 1213 real(dp),intent(in) :: elph2_imagden,esmear 1214 type(MPI_type),intent(inout) :: mpi_enreg 1215 type(datafiles_type), intent(in) :: dtfil 1216 type(pseudopotential_type), intent(inout) :: psps 1217 !arrays 1218 type(dataset_type), intent(in) :: dtset 1219 integer,intent(in) :: clflg(3,mpert) 1220 real(dp),intent(in) :: doccde(dtset%mband*dtset%nkpt*dtset%nsppol) 1221 real(dp),intent(in) :: eigen0(nkpt_rbz*mband*nsppol) 1222 real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) 1223 real(dp),intent(in) :: eigenq(nkpt_rbz*mband*nsppol) 1224 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 1225 real(dp),intent(inout) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq) 1226 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 1227 real(dp),intent(inout),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert) 1228 real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:) 1229 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 1230 type(hdr_type),intent(in) :: hdr0 1231 type(hdr_type),intent(in),optional :: hdr_fine 1232 1233 !Local variables------------------------------- 1234 !tolerance for non degenerated levels 1235 !scalars 1236 integer :: band2tot_index,band_index,bantot,bandtot_index,iband,idir1,idir2 1237 integer :: ikpt,ipert1,ipert2,isppol,jband,nkpt_sub,ikpt2,ncid 1238 !integer :: ipw 1239 character(len=fnlen) :: dscrpt,fname 1240 integer :: master,me,spaceworld,ierr 1241 integer :: mband_mem, mpert_ 1242 ! real(dp),parameter :: etol=1.0d-6 1243 real(dp),parameter :: etol=1.0d-7 1244 !real(dp),parameter :: etol=zero 1245 real(dp) :: ar,ai,deltae,den,eig1_i1,eig1_i2,eigen_corr 1246 real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av 1247 real(dp) :: eig2_diar,eigbrd_i,eigbrd_r,wgt_int 1248 !character(len=500) :: message 1249 logical :: remove_inv,test_do_band 1250 type(crystal_t) :: Crystal 1251 type(ebands_t) :: Bands 1252 !type(eigr2d_t) :: eigr2d,eigi2d 1253 type(fan_t) :: fan2d 1254 type(gkk_t) :: gkk2d 1255 type(ddb_hdr_type) :: ddb_hdr 1256 type(ddb_type) :: ddb 1257 !arrays 1258 integer,allocatable :: flg(:,:,:,:) 1259 integer, allocatable :: nband_rbz(:) 1260 integer,pointer :: kpt_fine_sub(:) 1261 real(dp) :: tsec(2) 1262 real(dp),allocatable :: center(:) 1263 real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol) 1264 real(dp),allocatable :: fan(:,:,:,:,:,:,:) 1265 real(dp),allocatable :: gkk(:,:,:,:,:) 1266 real(dp),allocatable :: eig2nkq_tmp(:,:,:,:,:,:,:) 1267 real(dp),allocatable :: smdfun(:,:) 1268 real(dp),pointer :: wgt_sub(:) 1269 1270 ! ********************************************************************* 1271 1272 !Init parallelism 1273 master =0 1274 spaceworld=mpi_enreg%comm_cell 1275 me=mpi_enreg%me_kpt 1276 1277 !Init interpolation method 1278 if(present(eigenq_fine))then 1279 ABI_MALLOC(center,(3)) 1280 end if 1281 1282 call timab(148,1,tsec) 1283 1284 if(nsppol==2)then 1285 ABI_COMMENT('nsppol=2 is still under development. Be careful when using it ...') 1286 end if 1287 1288 band2tot_index =0 1289 bandtot_index=0 1290 band_index=0 1291 1292 if(xmpi_paral==1) then 1293 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol)) 1294 ABI_MALLOC(nband_rbz,(nkpt_rbz*nsppol)) 1295 if (allocated(mpi_enreg%my_kpttab)) then 1296 ABI_FREE(mpi_enreg%my_kpttab) 1297 end if 1298 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz)) 1299 ! Assume the number of bands is the same for all k points. 1300 nband_rbz(:)=mband 1301 call distrb2(mband,mband_mem,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg) 1302 end if 1303 1304 if(ieig2rf == 4 ) then 1305 ABI_MALLOC_OR_DIE(fan,(2*mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq,mband), ierr) 1306 fan(:,:,:,:,:,:,:) = zero 1307 ABI_MALLOC_OR_DIE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr) 1308 eig2nkq_tmp(:,:,:,:,:,:,:) = zero 1309 ! This is not efficient because double the memory. Alternative: use buffer and 1310 ! print part by part. 1311 eig2nkq_tmp = eig2nkq 1312 if(present(eigbrd))then 1313 eigbrd(:,:,:,:,:,:,:)=zero 1314 end if 1315 eigen_corr = 0 1316 end if 1317 1318 if(ieig2rf == 5 ) then 1319 ABI_MALLOC_OR_DIE(gkk,(2*mband*nsppol,dtset%nkpt,3,natom,mband), ierr) 1320 gkk(:,:,:,:,:) = zero 1321 ABI_MALLOC_OR_DIE(eig2nkq_tmp,(2,mband*nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq), ierr) 1322 eig2nkq_tmp(:,:,:,:,:,:,:) = zero 1323 ! This is not efficient because double the memory. Alternative: use buffer and 1324 ! print part by part. 1325 eig2nkq_tmp = eig2nkq 1326 if(present(eigbrd))then 1327 eigbrd(:,:,:,:,:,:,:)=zero 1328 end if 1329 eigen_corr = 0 1330 end if 1331 1332 do isppol=1,nsppol 1333 do ikpt =1,nkpt_rbz 1334 1335 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then 1336 band2tot_index = band2tot_index + 2*mband**2 1337 bandtot_index = bandtot_index + mband 1338 cycle 1339 end if 1340 1341 if(present(eigenq_fine))then 1342 write(std_out,*) 'Start of the energy denominator interpolation method.' 1343 nkpt_sub = 0 1344 ! center is the k+q point around which we will average the kpt_fine 1345 center = hdr0%kptns(:,ikpt)+ dtset%qptn(:) 1346 1347 call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,kpt_fine_sub,nkpt_sub,wgt_sub) 1348 write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid & 1349 & around the k+Q point ',center,' is:',nkpt_sub 1350 write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub) 1351 end if 1352 1353 if(smdelta >0) then !broadening 1354 if(.not.allocated(smdfun)) then 1355 ABI_MALLOC(smdfun,(mband,mband)) 1356 end if 1357 smdfun(:,:) = zero 1358 do iband=1,mband 1359 eigen(iband) = eigen0(iband+bandtot_index) 1360 eigen_prime(iband) =eigenq(iband+bandtot_index) 1361 end do 1362 if(esmear>tol6) then 1363 call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun) 1364 end if 1365 end if 1366 1367 ipert1=1 ! Suppose all perturbations lead to the same number of planewaves 1368 1369 do iband=1,bdeigrf 1370 1371 ! If the k point and band belong to me, compute the contribution 1372 test_do_band=.true. 1373 if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false. 1374 1375 if(test_do_band)then 1376 ! ------------------------------------------------------------------------------------------------------! 1377 ! ------- ieig2rf ==3 : Non dynamic traditional AHC theory with Sternheimer (computed in eig2stern.F90)-! 1378 ! ------------------------------------------------------------------------------------------------------! 1379 ! Note that ieig2rf==4 and ieig2rf==5 also goes into that part only for later printing of the ZPR in the ouput of abinit 1380 ! later in the code 1381 if(ieig2rf==3 .or. ieig2rf==4 .or. ieig2rf==5) then 1382 do ipert1=1,npert 1383 do idir1=1,3 1384 if(clflg(idir1,ipert1)==0) cycle 1385 do ipert2=1,npert 1386 do idir2=1,3 1387 if(clflg(idir2,ipert2)==0)cycle 1388 eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero 1389 do jband=1,mband 1390 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1391 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 1392 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1393 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 1394 ! If no interpolation, fallback on to the previous 1395 ! implementation 1396 if(.not. present(eigenq_fine))then 1397 deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index) 1398 end if 1399 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 1400 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 1401 1402 ! Sum over all active space to retrieve the diagonal gauge 1403 ! if(abs(deltae)>etol) then ! This is commented because 1404 ! there is no problem with divergencies with elph2_imag != 0 1405 if( present(eigenq_fine))then 1406 den_av = zero 1407 wgt_int = zero 1408 do ikpt2=1,nkpt_sub 1409 deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)& 1410 & -eigen0(iband+bandtot_index) 1411 den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2) 1412 wgt_int = wgt_int+wgt_sub(ikpt2) 1413 end do 1414 den = den_av/wgt_int 1415 else 1416 if(abs(elph2_imagden) < etol) then 1417 if(abs(deltae)>etol) then 1418 den=-one/(deltae**2+elph2_imagden**2) 1419 else 1420 den= zero 1421 end if 1422 else 1423 den=-one/(deltae**2+elph2_imagden**2) 1424 end if 1425 end if 1426 1427 if( present(eigenq_fine))then 1428 eig2_diar=eig2_diar+ar*den 1429 eig2_diai=eig2_diai+ai*den 1430 else 1431 eig2_diar=eig2_diar+ar*deltae*den 1432 eig2_diai=eig2_diai+ai*deltae*den 1433 end if 1434 1435 if(present(eigbrd))then 1436 if(smdelta >0) then !broadening 1437 eigbrd_r = eigbrd_r + ar*smdfun(iband,jband) 1438 eigbrd_i = eigbrd_i + ai*smdfun(iband,jband) 1439 end if 1440 end if 1441 end do !jband 1442 1443 ! Store the contribution 1444 eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = & 1445 & eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diar 1446 eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = & 1447 & eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diai 1448 1449 if(present(eigbrd))then 1450 if(smdelta >0) then !broadening 1451 eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r 1452 eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i 1453 end if 1454 end if 1455 1456 end do !idir2 1457 end do !ipert2 1458 end do !idir1 1459 end do !ipert1 1460 end if !ieig2rf 3 1461 1462 ! -------------------------------------------------------------------------------------------! 1463 ! ------- ieig2rf ==4 Dynamic AHC using second quantization and Sternheimer from eig2stern -! 1464 ! -------------------------------------------------------------------------------------------! 1465 if(ieig2rf ==4 ) then 1466 do ipert1=1,npert 1467 do idir1=1,3 1468 if(clflg(idir1,ipert1)==0) cycle 1469 do ipert2=1,npert 1470 do idir2=1,3 1471 if(clflg(idir2,ipert2)==0)cycle 1472 do jband=1,mband 1473 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1474 eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2) 1475 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1476 eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC 1477 ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2 1478 ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2 1479 ! Store the contribution 1480 fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = & 1481 & fan(2*iband-1+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ar 1482 fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) = & 1483 & fan(2*iband+2*band_index,ikpt,idir1,ipert1,idir2,ipert2,jband) + ai 1484 end do !jband 1485 end do !idir2 1486 end do !ipert2 1487 end do !idir1 1488 end do !ipert1 1489 end if !ieig2rf 4 1490 ! --------------------------------------------------------------------------------! 1491 ! ------- ieig2rf ==5 Dynamic AHC with Sternheimer from eig2stern but print GKK -! 1492 ! --------------------------------------------------------------------------------! 1493 if(ieig2rf ==5 ) then 1494 do ipert1=1,npert 1495 do idir1=1,3 1496 if(clflg(idir1,ipert1)==0) cycle 1497 do jband=1,mband 1498 eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1499 eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1) 1500 ! Store the contribution 1501 gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) = & 1502 & gkk(2*iband-1+2*band_index,ikpt,idir1,ipert1,jband) + eig1_r1 1503 gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) = & 1504 & gkk(2*iband+2*band_index,ikpt,idir1,ipert1,jband) + eig1_i1 1505 end do !jband 1506 end do !idir1 1507 end do !ipert1 1508 end if !ieig2rf 5 1509 end if ! Selection of processor 1510 end do !iband 1511 1512 band2tot_index = band2tot_index + 2*mband**2 1513 bandtot_index = bandtot_index + mband 1514 1515 if(present(eigenq_fine))then 1516 ABI_FREE(kpt_fine_sub) ! Deallocate the variable 1517 ABI_FREE(wgt_sub) 1518 end if 1519 1520 end do !ikpt 1521 band_index = band_index + mband 1522 end do !isppol 1523 1524 !Accumulate eig2nkq and/or eigbrd 1525 if(xmpi_paral==1) then 1526 if(ieig2rf == 3) then 1527 call xmpi_sum(eig2nkq,spaceworld,ierr) 1528 end if 1529 if(ieig2rf == 4) then 1530 call xmpi_sum(eig2nkq,spaceworld,ierr) 1531 call xmpi_sum(eig2nkq_tmp,spaceworld,ierr) 1532 call xmpi_sum(fan,spaceworld,ierr) 1533 end if 1534 if(ieig2rf == 5) then 1535 call xmpi_sum(eig2nkq,spaceworld,ierr) 1536 call xmpi_sum(eig2nkq_tmp,spaceworld,ierr) 1537 call xmpi_sum(gkk,spaceworld,ierr) 1538 end if 1539 if(present(eigbrd) .and. (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5))then 1540 if(smdelta >0) then 1541 call xmpi_sum(eigbrd,spaceworld,ierr) 1542 end if 1543 end if 1544 ABI_FREE(nband_rbz) 1545 ABI_FREE(mpi_enreg%proc_distrb) 1546 ABI_FREE(mpi_enreg%my_kpttab) 1547 end if 1548 1549 if(ieig2rf > 2) then 1550 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.' 1551 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 1552 band_index = 0 1553 do isppol=1,dtset%nsppol 1554 do idir1=1,3 1555 do idir2=1,3 1556 ar=eig2nkq(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 1557 ai=eig2nkq(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 1558 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 1559 end do ! idir2 1560 end do ! idir1 1561 band_index = band_index + mband 1562 write(ab_out,'(a)')' ' 1563 end do 1564 end if 1565 1566 if(present(eigbrd))then 1567 if(smdelta >0) then !broadening 1568 write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.' 1569 write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.' 1570 band_index = 0 1571 do isppol=1,dtset%nsppol 1572 do idir1=1,3 1573 do idir2=1,3 1574 ar=eigbrd(1,1+band_index,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero 1575 ai=eigbrd(2,1+band_index,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero 1576 write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 1577 end do 1578 end do 1579 band_index = band_index + mband 1580 write(ab_out,'(a)')' ' 1581 end do 1582 end if 1583 end if 1584 1585 if(allocated(smdfun)) then 1586 ABI_FREE(smdfun) 1587 end if 1588 if(present(eigenq_fine))then 1589 ABI_FREE(center) 1590 end if 1591 1592 if (ieig2rf == 3 .or. ieig2rf == 4 .or. ieig2rf == 5) then 1593 1594 mpert_ = dtset%natom 1595 1596 ! Initialize perturbation flags 1597 ABI_MALLOC(flg,(3,mpert_,3,mpert_)) 1598 flg = one 1599 1600 ! Initialize ddb object 1601 call ddb%init(dtset, 1, mpert_, & 1602 mband=bdeigrf,& 1603 nkpt=nkpt_rbz,& 1604 kpt=dtset%kptns(1:3,1:nkpt_rbz),& 1605 with_d2eig=.true.) 1606 1607 ! Create the ddb header 1608 dscrpt=' Note : temporary (transfer) database ' 1609 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,& 1610 mpert=mpert_,& 1611 xred=xred,occ=occ_rbz,& 1612 mband=bdeigrf / dtset%nsppol,& 1613 nkpt=nkpt_rbz,& 1614 kpt=dtset%kptns(:,1:nkpt_rbz)) 1615 1616 ! Set d2eig data 1617 call ddb%set_qpt(1, dtset%qptn) 1618 call ddb%set_d2eig_reshape(1, eig2nkq, flg) 1619 1620 ! Open the file and write header 1621 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 1622 1623 ! Write d2eig data block 1624 call ddb_hdr%open_write(dtfil%fnameabo_eigr2d, with_psps=1, comm=mpi_enreg%comm_world) 1625 call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world) 1626 1627 ! close and free memory 1628 call ddb_hdr%close() 1629 call ddb_hdr%free() 1630 call ddb%free() 1631 1632 ABI_FREE(flg) 1633 1634 end if 1635 1636 ! print _FAN file for this perturbation. Note that the Fan file will only be produced if 1637 ! abinit is compiled with netcdf. 1638 1639 ! Initialize crystal structure for FAN.nc and GKK.nc files 1640 remove_inv=.false. 1641 if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. 1642 call crystal_init(dtset%amu_orig(:,1),Crystal,dtset%spgroup,dtset%natom,dtset%npsp,psps%ntypat, & 1643 & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,& 1644 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr0%title,& 1645 & dtset%symrel,dtset%tnons,dtset%symafm) 1646 ! Electronic band energies. 1647 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 1648 call ebands_init(bantot,Bands,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 1649 & doccde,eigen0,hdr0%istwfk,hdr0%kptns,& 1650 & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,& 1651 & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,& 1652 & hdr0%cellcharge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, & 1653 & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk) 1654 ! 1655 if(ieig2rf == 4 ) then 1656 ! Output of the Fan.nc file. 1657 #ifdef HAVE_NETCDF 1658 fname = strcat(dtfil%filnam_ds(4),"_FAN.nc") 1659 call fan_init(fan,fan2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 1660 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating FAN file") 1661 NCF_CHECK(crystal%ncwrite(ncid)) 1662 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 1663 call fan_ncwrite(fan2d,dtset%qptn(:),dtset%wtq, ncid) 1664 NCF_CHECK(nf90_close(ncid)) 1665 #else 1666 ABI_ERROR("Dynamical calculation with ieig2rf 4 only work with NETCDF support.") 1667 ABI_UNUSED(ncid) 1668 #endif 1669 ABI_FREE(fan) 1670 ABI_FREE(eig2nkq_tmp) 1671 end if 1672 ! print _GKK.nc file for this perturbation. Note that the GKK file will only be produced if 1673 ! abinit is compiled with netcdf. 1674 if(ieig2rf == 5 ) then 1675 ! Output of the GKK.nc file. 1676 #ifdef HAVE_NETCDF 1677 fname = strcat(dtfil%filnam_ds(4),"_GKK.nc") 1678 call gkk_init(gkk,gkk2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom,3) 1679 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file") 1680 NCF_CHECK(crystal%ncwrite(ncid)) 1681 NCF_CHECK(ebands_ncwrite(Bands, ncid)) 1682 call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid) 1683 NCF_CHECK(nf90_close(ncid)) 1684 #else 1685 ABI_ERROR("Dynamical calculation with ieig2rf 5 only work with NETCDF support.") 1686 ABI_UNUSED(ncid) 1687 #endif 1688 ABI_FREE(gkk) 1689 ABI_FREE(eig2nkq_tmp) 1690 end if 1691 1692 ! print _EIGI2D file for this perturbation 1693 if (ieig2rf /= 5 ) then 1694 if(smdelta>0) then 1695 1696 ABI_MALLOC(flg,(3,mpert_,3,mpert_)) 1697 flg = one 1698 call ddb%init(dtset, 1, mpert_, & 1699 mband=bdeigrf,& 1700 nkpt=nkpt_rbz,& 1701 kpt=dtset%kptns(:,1:nkpt_rbz),& 1702 with_d2eig=.true.) 1703 1704 ! Create the ddb header 1705 dscrpt=' Note : temporary (transfer) database ' 1706 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,& 1707 mpert=mpert_,& 1708 xred=xred,occ=occ_rbz,& 1709 mband=bdeigrf / dtset%nsppol,& 1710 nkpt=nkpt_rbz,& 1711 kpt=dtset%kptns(:,1:nkpt_rbz)) 1712 1713 call ddb%set_qpt(1, dtset%qptn) 1714 call ddb%set_d2eig_reshape(1, eigbrd, flg, blktyp=BLKTYP_d2eig_im) 1715 1716 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 1717 1718 call ddb_hdr%open_write(dtfil%fnameabo_eigi2d, with_psps=1) 1719 call ddb%write_d2eig(ddb_hdr, 1) 1720 1721 call ddb_hdr%close() 1722 call ddb_hdr%free() 1723 call ddb%free() 1724 1725 ABI_FREE(flg) 1726 1727 end if !smdelta 1728 end if 1729 !end if ! master 1730 1731 if (allocated(fan)) then 1732 ABI_FREE(fan) 1733 end if 1734 if (allocated(eig2nkq_tmp)) then 1735 ABI_FREE(eig2nkq_tmp) 1736 end if 1737 if (allocated(gkk)) then 1738 ABI_FREE(gkk) 1739 end if 1740 1741 call crystal%free() 1742 call ebands_free(Bands) 1743 call fan_free(fan2d) 1744 call gkk_free(gkk2d) 1745 1746 call timab(148,2,tsec) 1747 1748 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)
SOURCE
309 subroutine eigr2d_free(eigr2d) 310 311 !Arguments ------------------------------------ 312 !scalars 313 type(eigr2d_t),intent(inout) :: eigr2d 314 ! ************************************************************************* 315 DBG_ENTER("COLL") 316 317 !Deallocate all components of bstruct 318 319 if (allocated(eigr2d%eigr2d)) then 320 ABI_FREE(eigr2d%eigr2d) 321 end if 322 323 DBG_EXIT("COLL") 324 325 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
SOURCE
189 subroutine eigr2d_init(eig2nkq,eigr2d,mband,nsppol,nkpt,natom) 190 191 !Arguments ------------------------------------ 192 !scalars 193 integer,intent(in) ::mband,nsppol,nkpt,natom 194 type(eigr2d_t),intent(out) :: eigr2d 195 !arrays 196 real(dp), intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom) 197 198 ! ************************************************************************* 199 200 eigr2d%mband = mband 201 eigr2d%nsppol = nsppol 202 eigr2d%nkpt = nkpt 203 eigr2d%natom = natom 204 205 ABI_MALLOC(eigr2d%eigr2d ,(2,mband*nsppol,nkpt,3,natom,3,natom)) 206 eigr2d%eigr2d=eig2nkq 207 208 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
SOURCE
228 subroutine eigr2d_ncwrite(eigr2d,iqpt,wtq,ncid) 229 230 !Arguments ------------------------------------ 231 !scalars 232 integer,intent(in) ::ncid 233 real(dp),intent(in) :: iqpt(3),wtq 234 type(eigr2d_t),intent(in) :: eigr2d 235 236 !Local variables------------------------------- 237 #ifdef HAVE_NETCDF 238 integer :: ncerr 239 integer :: cplex,cart_dir,one_dim 240 character(len=200) :: temp 241 242 ! ************************************************************************* 243 244 ! ============================================== 245 ! === Write the dimensions specified by ETSF === 246 ! ============================================== 247 one_dim=1; cplex=2; cart_dir=3 248 249 ncerr = nctk_def_dims(ncid, [& 250 nctkdim_t('max_number_of_states', eigr2d%mband),& 251 nctkdim_t('number_of_spins', eigr2d%nsppol),& 252 nctkdim_t('number_of_kpoints', eigr2d%nkpt),& 253 nctkdim_t('number_of_atoms', eigr2d%natom),& 254 nctkdim_t('number_of_cartesian_directions', cart_dir),& 255 nctkdim_t('current_one_dim', one_dim),& 256 nctkdim_t('cplex', cplex),& 257 nctkdim_t('product_mband_nsppol', eigr2d%mband*eigr2d%nsppol)], defmode=.True.) 258 NCF_CHECK(ncerr) 259 260 temp='cplex,product_mband_nsppol,number_of_kpoints,number_of_cartesian_directions,number_of_atoms,' //& 261 'number_of_cartesian_directions , number_of_atoms' 262 ncerr = nctk_def_arrays(ncid, [& 263 nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'), & 264 nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'), & 265 nctkarr_t('second_derivative_eigenenergies', "dp", temp )]) 266 ! nctkarr_t('second_derivative_eigenenergies', "dp",& 267 ! &'cplex, product_mband_nsppol, number_of_kpoints, number_of_cartesian_directions, number_of_atoms,& 268 ! &number_of_cartesian_directions, number_of_atoms')]) 269 NCF_CHECK(ncerr) 270 271 ! Write data 272 NCF_CHECK(nctk_set_datamode(ncid)) 273 NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt)) 274 NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq)) 275 NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies'), eigr2d%eigr2d)) 276 277 #else 278 ABI_ERROR("ETSF-IO support is not activated. ") 279 #endif 280 281 282 contains 283 integer function vid(vname) 284 character(len=*),intent(in) :: vname 285 vid = nctk_idname(ncid, vname) 286 end function vid 287 288 end subroutine eigr2d_ncwrite
m_eig2d/eigr2d_t [ Types ]
NAME
eig2d_t
FUNCTION
It contains information about the second-order derivative of the eigenenergies wrt atomic displacement
SOURCE
82 type,public :: eigr2d_t 83 84 ! WARNING : if you modify this datatype, please check whether there might be 85 ! creation/destruction/copy routines, 86 ! declared in another part of ABINIT, that might need to take into account your 87 ! modification. 88 89 integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays) 90 integer :: nsppol ! number of spin-polarization 91 integer :: nkpt ! number of k points 92 integer :: natom ! number of atoms 93 94 real(dp),allocatable :: eigr2d(:,:,:,:,:,:,:) 95 ! eigr2d(2,mband*nsppol,nkpt,3,natom,3,natom) 96 ! Second-order derivative of eigenergies (real,im) at each 97 ! spin,band,k-point,dir1,dir2,natom1,natom2 . 98 99 100 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.
SOURCE
1952 subroutine elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_corr,gprimd,mband,natom,nkpt,nsppol,option,phfrq,prtvol) 1953 1954 !Arguments ------------------------------------ 1955 !scalars 1956 integer,intent(in) :: dim_eig2nkq,mband,natom,nkpt,nsppol,option,prtvol 1957 1958 !arrays 1959 real(dp) :: gprimd(3,3) 1960 real(dp),intent(in) :: displ(2*3*natom*3*natom) 1961 real(dp),intent(in) :: eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) 1962 real(dp),intent(in) :: phfrq(3*natom) 1963 real(dp),intent(out) :: eigen_corr(mband*nkpt*nsppol) 1964 1965 !Local variables------------------------------- 1966 !scalars 1967 integer,parameter :: neigs_per_line=6 1968 integer :: iatom1,iatom2,idir1,idir2,iband,ikpt,imode,index,isppol, imin, ii 1969 real(dp) :: d_at1_dir1_re,d_at1_dir1_im 1970 real(dp) :: d_at1_dir2_re,d_at1_dir2_im 1971 real(dp) :: d_at2_dir1_re,d_at2_dir1_im 1972 real(dp) :: d_at2_dir2_re,d_at2_dir2_im 1973 real(dp) :: e2_im,e2_re 1974 real(dp), allocatable :: eigen_corr_mode(:) 1975 character(len=500) :: message 1976 character(len=20) :: eig_format, line_format 1977 !arrays 1978 real(dp) :: displ2cart(2,3,3),displ2red(2,3,3),tmp_displ2(2,3,3) 1979 1980 ! ********************************************************************* 1981 1982 1983 if(option/=1 .and. option/=2)then 1984 write(message,'(a,i0)')' The argument option should be 1 or 2, while it is found that option=',option 1985 ABI_BUG(message) 1986 end if 1987 1988 !printing options 1989 eig_format='f16.8' 1990 write(line_format,'(a,i1,a6,a)') '(',neigs_per_line,eig_format,')' 1991 1992 if (prtvol > 4) then 1993 write(message,'(a,a)')ch10,' ================================================================================' 1994 call wrtout(ab_out_default,message,'COLL') 1995 if (option==1) then 1996 write(message,'(a)') ' ---- Begin Fan contributions to eigenvalues renormalization by mode ----' 1997 call wrtout(ab_out_default,message,'COLL') 1998 else if (option==2) then 1999 write(message,'(a)') ' ---- Begin DDW contributions to eigenvalues renormalization by mode ----' 2000 call wrtout(ab_out_default,message,'COLL') 2001 end if 2002 end if 2003 2004 ABI_MALLOC(eigen_corr_mode,(mband*nkpt*nsppol)) 2005 2006 eigen_corr(:)=zero 2007 do imode=1,3*natom 2008 eigen_corr_mode(:)=zero 2009 2010 if (phfrq(imode)>tol6) then 2011 do iatom1=1,natom 2012 do iatom2=1,natom 2013 2014 do idir1=1,3 2015 do idir2=1,3 2016 ! Compute the mean cartesian displacements 2017 d_at1_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1)))) 2018 d_at1_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1)))) 2019 d_at2_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1)))) 2020 d_at2_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1)))) 2021 2022 if(option==1)then 2023 ! Compute the mean displacement correlation at T=0. 2024 ! Consistent with Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point. 2025 ! but generalized to two different atoms. Note that the complex conjugate is taken on the second direction. 2026 displ2cart(1,idir1,idir2)=(d_at1_dir1_re*d_at2_dir2_re+ & 2027 & d_at1_dir1_im*d_at2_dir2_im )/(two*phfrq(imode)) 2028 displ2cart(2,idir1,idir2)=(d_at1_dir1_im*d_at2_dir2_re- & 2029 & d_at1_dir1_re*d_at2_dir2_im )/(two*phfrq(imode)) 2030 else if(option==2)then 2031 ! Compute the mean square displacement correlation of each atom at T=0, and take mean over iatom1 and iatom2. 2032 ! See Eqs.(7) and (8) of PRB51, 8610 (1995) [[cite:Lee1995]], specialized for the contribution of one q point. 2033 ! Note that the complex conjugate is taken on the second direction. 2034 ! Also, note the overall negative sign, to make it opposite to the Fan term. 2035 d_at1_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1)))) 2036 d_at1_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom1-1 +natom*(imode-1)))) 2037 d_at2_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1)))) 2038 d_at2_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom2-1 +natom*(imode-1)))) 2039 displ2cart(1,idir1,idir2)=-(d_at1_dir1_re*d_at1_dir2_re+ & 2040 & d_at1_dir1_im*d_at1_dir2_im+ & 2041 & d_at2_dir1_re*d_at2_dir2_re+ & 2042 & d_at2_dir1_im*d_at2_dir2_im )/(four*phfrq(imode)) 2043 displ2cart(2,idir1,idir2)=-(d_at1_dir1_im*d_at1_dir2_re- & 2044 & d_at1_dir1_re*d_at1_dir2_im+ & 2045 & d_at2_dir1_im*d_at2_dir2_re- & 2046 & d_at2_dir1_re*d_at2_dir2_im )/(four*phfrq(imode)) 2047 end if 2048 end do 2049 end do 2050 ! Switch to reduced coordinates in two steps 2051 tmp_displ2(:,:,:)=zero 2052 do idir1=1,3 2053 do idir2=1,3 2054 tmp_displ2(:,:,idir1)=tmp_displ2(:,:,idir1)+displ2cart(:,:,idir2)*gprimd(idir2,idir1) 2055 end do 2056 end do 2057 displ2red(:,:,:)=zero 2058 do idir1=1,3 2059 do idir2=1,3 2060 displ2red(:,idir1,:)=displ2red(:,idir1,:)+tmp_displ2(:,idir2,:)*gprimd(idir2,idir1) 2061 end do 2062 end do 2063 ! Compute the T=0 shift due to this q point 2064 do idir1=1,3 2065 do idir2=1,3 2066 do ikpt=1,nkpt 2067 do isppol=1,nsppol 2068 do iband=1,mband 2069 index=iband+mband*(isppol-1 + nsppol*(ikpt-1)) 2070 e2_re=eig2nkq(1,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2) 2071 e2_im=eig2nkq(2,iband+mband*(isppol-1),ikpt,idir1,iatom1,idir2,iatom2) 2072 eigen_corr(index)=eigen_corr(index)+& 2073 & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2) 2074 eigen_corr_mode(index)=eigen_corr_mode(index)+& 2075 & e2_re*displ2red(1,idir1,idir2)-e2_im*displ2red(2,idir1,idir2) 2076 end do ! band 2077 end do ! spin 2078 end do ! kpt 2079 end do ! dir2 2080 end do ! dir1 2081 end do ! atom2 2082 end do ! atom1 2083 end if 2084 2085 if (prtvol > 4) then 2086 ! Print the corrections by mode 2087 write(message,'(a,i1)') ' imode= ',imode 2088 call wrtout(ab_out_default,message,'COLL') 2089 2090 do ikpt=1,nkpt 2091 do isppol=1,nsppol 2092 write(message,'(a,i4,a,i1)')' ikpt= ',ikpt,' ispin= ',isppol 2093 call wrtout(ab_out_default,message,'COLL') 2094 2095 imin = mband * (isppol-1 + nsppol*(ikpt-1)) 2096 do ii=0, (mband-1)/neigs_per_line 2097 write(message, line_format) (eigen_corr_mode(iband+imin), & 2098 & iband = 1 + ii * neigs_per_line, min(mband, (ii+1)*neigs_per_line)) 2099 call wrtout(ab_out_default,message,'COLL') 2100 end do 2101 end do 2102 end do 2103 end if 2104 2105 end do ! mode 2106 2107 if (prtvol > 4) then 2108 if (option==1) then 2109 write(message,'(a)') ' ---- End Fan contribution to eigenvalues renormalization by mode ----' 2110 call wrtout(ab_out_default,message,'COLL') 2111 else if (option==2) then 2112 write(message,'(a)') ' ---- End DDW contribution to eigenvalues renormalization by mode ----' 2113 call wrtout(ab_out_default,message,'COLL') 2114 end if 2115 write(message,'(a,a)')' ================================================================================', ch10 2116 call wrtout(ab_out_default,message,'COLL') 2117 end if 2118 2119 ABI_FREE(eigen_corr_mode) 2120 2121 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)
SOURCE
597 subroutine fan_free(fan2d) 598 599 !Arguments ------------------------------------ 600 !scalars 601 type(fan_t),intent(inout) :: fan2d 602 ! ************************************************************************* 603 DBG_ENTER("COLL") 604 605 !Deallocate all components of bstruct 606 607 if (allocated(fan2d%fan2d)) then 608 ABI_FREE(fan2d%fan2d) 609 end if 610 611 DBG_EXIT("COLL") 612 613 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
SOURCE
350 subroutine fan_init(fan,fan2d,mband,nsppol,nkpt,natom) 351 352 !Arguments ------------------------------------ 353 !scalars 354 integer,intent(in) ::mband,nsppol,nkpt,natom 355 type(fan_t),intent(out) :: fan2d 356 !arrays 357 real(dp), intent(in) :: fan(2*mband*nsppol,nkpt,3,natom,3,natom,mband) 358 ! ************************************************************************* 359 360 fan2d%mband = mband 361 fan2d%nsppol = nsppol 362 fan2d%nkpt = nkpt 363 fan2d%natom = natom 364 365 ABI_MALLOC(fan2d%fan2d,(2*mband*nsppol,nkpt,3,natom,3,natom,mband)) 366 fan2d%fan2d=fan 367 368 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
SOURCE
432 subroutine fan_ncwrite(fan2d,iqpt,wtq,ncid) 433 434 !Arguments ------------------------------------ 435 !scalars 436 integer,intent(in) ::ncid 437 real(dp),intent(in) :: iqpt(3),wtq 438 type(fan_t),intent(in) :: fan2d 439 440 !Local variables------------------------------- 441 #ifdef HAVE_NETCDF 442 integer :: ncerr 443 integer :: cplex,cart_dir,one_dim 444 character(len=200) :: temp 445 446 447 ! ************************************************************************* 448 449 ! ============================================== 450 ! === Write the dimensions specified by ETSF === 451 ! ============================================== 452 one_dim=1; cplex=2; cart_dir=3 453 454 ncerr = nctk_def_dims(ncid, [& 455 nctkdim_t('max_number_of_states',fan2d%mband),& 456 nctkdim_t('number_of_spins',fan2d%nsppol),& 457 nctkdim_t('number_of_kpoints',fan2d%nkpt),& 458 nctkdim_t('number_of_atoms',fan2d%natom),& 459 nctkdim_t('3_number_of_atoms',3*fan2d%natom),& ! TODO: not sure that variables can start with digits 460 nctkdim_t('number_of_cartesian_directions',cart_dir),& 461 nctkdim_t('current_one_dim',one_dim),& 462 nctkdim_t('cplex',cplex),& 463 nctkdim_t('product_mband_nsppol',fan2d%mband*fan2d%nsppol),& 464 nctkdim_t('product_mband_nsppol2',fan2d%mband*fan2d%nsppol*2) & 465 ], defmode=.True.) 466 NCF_CHECK(ncerr) 467 468 temp= 'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions,' //& 469 'number_of_atoms, number_of_cartesian_directions, number_of_atoms, max_number_of_states' 470 ncerr = nctk_def_arrays(ncid, [& 471 nctkarr_t('current_q_point', "dp", 'number_of_cartesian_directions'),& 472 nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'),& 473 nctkarr_t('second_derivative_eigenenergies_actif', "dp", temp )]) 474 ! nctkarr_t('second_derivative_eigenenergies_actif', "dp",& 475 ! &'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions,& 476 ! &number_of_atoms, number_of_cartesian_directions, number_of_atoms, max_number_of_states')]) 477 NCF_CHECK(ncerr) 478 479 ! Write data 480 NCF_CHECK(nctk_set_datamode(ncid)) 481 NCF_CHECK(nf90_put_var(ncid, vid('current_q_point'), iqpt)) 482 NCF_CHECK(nf90_put_var(ncid, vid('current_q_point_weight'), wtq)) 483 NCF_CHECK(nf90_put_var(ncid, vid('second_derivative_eigenenergies_actif'), fan2d%fan2d)) 484 485 #else 486 ABI_ERROR("netcdf support is not activated. ") 487 #endif 488 489 contains 490 integer function vid(vname) 491 character(len=*),intent(in) :: vname 492 vid = nctk_idname(ncid, vname) 493 end function vid 494 495 end subroutine fan_ncwrite
m_eig2d/fan_t [ Types ]
NAME
fan_t
FUNCTION
It contains information about the second-order derivative of the eigenenergies wrt atomic displacement
SOURCE
113 type,public :: fan_t 114 115 ! WARNING : if you modify this datatype, please check whether there might be 116 ! creation/destruction/copy routines, 117 ! declared in another part of ABINIT, that might need to take into account your 118 ! modification. 119 120 integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays) 121 integer :: nsppol ! number of spin-polarization 122 integer :: nkpt ! number of k points 123 integer :: natom ! number of atoms 124 125 real(dp),allocatable :: fan2d(:,:,:,:,:,:,:) 126 ! fan2d(2*mband*nsppol,nkpt,3,natom,3,natom,mband) 127 ! Second-order derivative of the eigenergies (real,im) at each 128 ! ispin,iband(real,im),k-point,dir1,dir2,natom1,natom2,jband 129 130 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)
SOURCE
634 subroutine gkk_free(gkk2d) 635 636 !Arguments ------------------------------------ 637 !scalars 638 type(gkk_t),intent(inout) :: gkk2d 639 ! ************************************************************************* 640 DBG_ENTER("COLL") 641 642 !Deallocate all components of bstruct 643 644 if (allocated(gkk2d%gkk2d)) then 645 ABI_FREE(gkk2d%gkk2d) 646 end if 647 648 DBG_EXIT("COLL") 649 650 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
SOURCE
393 subroutine gkk_init(gkk,gkk2d,mband,nsppol,nkpt,natom,ncart) 394 395 !Arguments ------------------------------------ 396 !scalars 397 integer,intent(in) ::mband,nsppol,nkpt,natom,ncart 398 type(gkk_t),intent(out) :: gkk2d 399 !arrays 400 real(dp), intent(in) :: gkk(2*mband*nsppol,nkpt,ncart,natom,mband) 401 ! ************************************************************************* 402 403 gkk2d%mband = mband 404 gkk2d%nsppol = nsppol 405 gkk2d%nkpt = nkpt 406 gkk2d%natom = natom 407 gkk2d%ncart = ncart 408 409 ABI_MALLOC(gkk2d%gkk2d,(2*mband*nsppol,nkpt,ncart,natom,mband)) 410 gkk2d%gkk2d=gkk 411 412 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
SOURCE
515 subroutine gkk_ncwrite(gkk2d,iqpt,wtq,ncid) 516 517 !Arguments ------------------------------------ 518 !scalars 519 integer,intent(in) ::ncid 520 real(dp), intent(in) :: iqpt(3),wtq 521 type(gkk_t),intent(in) :: gkk2d 522 523 !Local variables------------------------------- 524 #ifdef HAVE_NETCDF 525 integer :: cplex,one_dim,ncerr,vid_ 526 527 ! ************************************************************************* 528 529 ! ============================================== 530 ! === Write the dimensions specified by ETSF === 531 ! ============================================== 532 one_dim=1; cplex=2 533 534 ncerr = nctk_def_dims(ncid, [ & 535 & nctkdim_t('max_number_of_states', gkk2d%mband), & 536 & nctkdim_t('number_of_spins', gkk2d%nsppol), & 537 & nctkdim_t('number_of_kpoints', gkk2d%nkpt), & 538 & nctkdim_t('number_of_atoms_for_gkk', gkk2d%natom), & 539 & nctkdim_t('3_number_of_atoms', 3*gkk2d%natom), & 540 & nctkdim_t('number_of_cartesian_directions_for_gkk', gkk2d%ncart), & 541 & nctkdim_t('current_one_dim', one_dim), & 542 & nctkdim_t('cplex', cplex), & 543 & nctkdim_t('product_mband_nsppol', gkk2d%mband*gkk2d%nsppol), & 544 & nctkdim_t('product_mband_nsppol2', gkk2d%mband*gkk2d%nsppol*2) & 545 & ], defmode=.True.) 546 NCF_CHECK(ncerr) 547 548 !arrays 549 ncerr = nctk_def_arrays(ncid, [& 550 & nctkarr_t('current_q_point', "dp", "number_of_cartesian_directions"), & 551 & nctkarr_t('current_q_point_weight', "dp", 'current_one_dim'), & 552 & nctkarr_t('second_derivative_eigenenergies_actif', "dp", & 553 & 'product_mband_nsppol2, number_of_kpoints, number_of_cartesian_directions_for_gkk,'// & 554 & 'number_of_atoms_for_gkk, max_number_of_states') & 555 & ]) 556 NCF_CHECK(ncerr) 557 558 NCF_CHECK(nctk_set_datamode(ncid)) 559 vid_=vid('current_q_point') 560 NCF_CHECK(nf90_put_var(ncid, vid_, iqpt)) 561 vid_=vid('current_q_point_weight') 562 NCF_CHECK(nf90_put_var(ncid, vid_, wtq)) 563 vid_=vid('second_derivative_eigenenergies_actif') 564 NCF_CHECK(nf90_put_var(ncid, vid_, gkk2d%gkk2d)) 565 566 #else 567 ABI_ERROR("netcdf support is not activated. ") 568 #endif 569 570 contains 571 integer function vid(vname) 572 character(len=*),intent(in) :: vname 573 vid = nctk_idname(ncid, vname) 574 end function vid 575 576 end subroutine gkk_ncwrite
m_eig2d/gkk_t [ Types ]
NAME
gkk_t
FUNCTION
It contains information about the second-order derivative of the eigenenergies wrt atomic displacement
SOURCE
143 type,public :: gkk_t 144 145 ! WARNING : if you modify this datatype, please check whether there might be 146 ! creation/destruction/copy routines, 147 ! declared in another part of ABINIT, that might need to take into account your 148 ! modification. 149 150 integer :: mband ! Max number of bands i.e MAXVAL(nband) (to dimension arrays) 151 integer :: nsppol ! number of spin-polarization 152 integer :: nkpt ! number of k points 153 integer :: natom ! number of atoms 154 integer :: ncart ! number of cartesian directions 155 156 real(dp),allocatable :: gkk2d(:,:,:,:,:) 157 ! gkk2d(2*mband*nsppol,nkpt,ncart,natom,mband) 158 ! Second-order derivative of the eigenergies (real,im) at each 159 ! ispin,iband(real,im),k-point,dir1,natom1,jband 160 161 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 NOTE This function is deprecated. One should write through ddb object instead.
SOURCE
1774 subroutine outbsd(bdeigrf,dtset,eig2nkq,mpert,nkpt_rbz,unitout) 1775 1776 !Arguments ------------------------------------ 1777 !scalars 1778 integer,intent(in) :: bdeigrf,mpert,nkpt_rbz,unitout 1779 type(dataset_type),intent(in) :: dtset 1780 !arrays 1781 real(dp),intent(in) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt_rbz,3,mpert,3,mpert) 1782 1783 !Local variables ------------------------- 1784 !scalars 1785 integer :: bandtot_index,iband,idir1,idir2,ikpt,ipert1,ipert2,isppol 1786 1787 ! ********************************************************************* 1788 1789 1790 !output information in this file 1791 write(unitout,*) 1792 write(unitout,'(a,i8)') ' 2nd eigenvalue derivatives - # elements :', 9*dtset%natom**2 1793 write(unitout,'(a,3es16.8,a)') ' qpt', dtset%qptn(:), ' 1.0' 1794 1795 !output RF eigenvalues 1796 1797 do ikpt=1,nkpt_rbz 1798 ! bandtot_index differs from zero only in the spin-polarized case 1799 bandtot_index=0 1800 write (unitout,'(a,3es16.8)') ' K-point:', dtset%kptns(:,ikpt) 1801 do isppol=1,dtset%nsppol 1802 do iband=1,bdeigrf 1803 write (unitout,'(a,i5)') ' Band:', iband+bandtot_index 1804 ! write (unitout,*) 'ipert1 ','idir1 ','ipert2 ','idir2 ','Real ','Im ' 1805 do ipert2=1,mpert 1806 do idir2=1,3 1807 do ipert1=1,mpert 1808 do idir1=1,3 1809 write (unitout,'(4i4,2d22.14)') idir1,ipert1,idir2,ipert2,& 1810 & eig2nkq(1,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2),& 1811 & eig2nkq(2,iband+bandtot_index,ikpt,idir1,ipert1,idir2,ipert2) 1812 end do !idir2 1813 end do !ipert2 1814 end do !idir1 1815 end do !ipert1 1816 end do !iband 1817 bandtot_index = bandtot_index + dtset%mband 1818 end do !isppol 1819 end do !ikpt 1820 1821 !close bsd file 1822 close (unitout) 1823 1824 end subroutine outbsd