TABLE OF CONTENTS
ABINIT/alignph [ Functions ]
NAME
alignph
FUNCTION
Construct linear combinations of the phonon eigendisplacements of degenerate modes in order to align the mode effective charges along the axes of the cartesian frame.
INPUTS
amu(ntypat)=mass of the atoms (atomic mass unit) displ(2,3*natom,3*natom)= the displacements of atoms in cartesian coordinates. The first index means either the real or the imaginary part, The second index runs on the direction and the atoms displaced The third index runs on the modes. d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates mpert =maximum number of ipert natom=number of atoms in unit cell ntypat=number of types of atoms phfrq(3*natom)=phonon frequencies (square root of the dynamical matrix eigenvalues, except if these are negative, and in this case, give minus the square root of the absolute value of the matrix eigenvalues). Hartree units. typat(natom)=integer label of each type of atom (1,2,...)
OUTPUT
displ(2,3*natom,3*natom)= the displacements of atoms in cartesian coordinates. The eigendisplacements of degenerate modes have been aligned along the cartesian axes.
SOURCE
640 subroutine alignph(amu,displ,d2cart,mpert,natom,ntypat,phfrq,typat) 641 642 !Arguments ------------------------------- 643 !scalars 644 integer,intent(in) :: mpert,natom,ntypat 645 !arrays 646 integer,intent(in) :: typat(natom) 647 real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),phfrq(3*natom) 648 real(dp),intent(inout) :: displ(2,3*natom,3*natom) 649 650 !Local variables ------------------------- 651 !scalars 652 integer,parameter :: master=0 653 integer :: i1,idir1,idir2,ii,imode,imodex,imodey,imodez,ipert1 654 real(dp) :: theta 655 !arrays 656 integer,allocatable :: deg(:) 657 real(dp) :: zvec(3,3),zvect(3,3) 658 real(dp),allocatable :: modez(:,:,:),modezabs(:),oscstr(:,:,:),vec(:,:),vect(:,:) 659 660 ! ********************************************************************* 661 662 !Get the oscillator strength and mode effective charge for each mode 663 ABI_MALLOC(oscstr,(2,3,3*natom)) 664 ABI_MALLOC(modez,(2,3,3*natom)) 665 ABI_MALLOC(modezabs,(3*natom)) 666 ABI_MALLOC(vec,(3*natom,3)) 667 ABI_MALLOC(vect,(3*natom,3)) 668 ABI_MALLOC(deg,(3*natom)) 669 670 write(std_out,'(a,a)')ch10,' alignph : before modifying the eigenvectors, mode number and mode effective charges :' 671 do imode=1,3*natom 672 modezabs(imode)=zero 673 do ii=1,2 674 do idir2=1,3 675 oscstr(ii,idir2,imode)=zero 676 modez(ii,idir2,imode)=zero 677 do idir1=1,3 678 do ipert1=1,natom 679 i1=idir1+(ipert1-1)*3 680 oscstr(ii,idir2,imode)=oscstr(ii,idir2,imode)+& 681 & displ(ii,i1,imode)*& 682 & d2cart(1,idir1,ipert1,idir2,natom+2) 683 modez(ii,idir2,imode)=modez(ii,idir2,imode)+& 684 & displ(ii,i1,imode)*& 685 & d2cart(1,idir1,ipert1,idir2,natom+2)*& 686 & sqrt(amu(typat(ipert1))*amu_emass) 687 end do 688 end do 689 if(abs(modez(ii,idir2,imode))>modezabs(imode))modezabs(imode)=abs(modez(ii,idir2,imode)) 690 end do 691 end do 692 write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode) 693 end do 694 695 !Find degenerate modes with non-zero mode effective charge 696 imode = 0 697 do while (imode < 3*natom) 698 imode = imode + 1 699 if (imode == 3*natom) then 700 deg(imode) = 1 701 else if (abs(phfrq(imode) - phfrq(imode+1)) > tol6 .or. modezabs(imode)<tol8 .or. modezabs(imode+1)<tol8) then 702 ! Differ by phonon frequency or zero mode effective charge 703 deg(imode) = 1 704 else 705 deg(imode) = 2 706 if (imode < 3*natom - 1) then 707 if (abs(phfrq(imode) - phfrq(imode+2)) < tol6 .and. modezabs(imode+2)>tol8) then 708 deg(imode) = 3 709 imode = imode + 1 710 end if 711 end if 712 imode = imode + 1 713 end if 714 end do 715 716 717 !In case of a degenerate mode, with non-zero mode effective charge, align the mode effective charge vector along 718 !the axes of the cartesian frame 719 imode = 1 720 do while (imode <= 3*natom) 721 722 write(std_out,'(a,a,i4,a,i2)')ch10,' Mode number ',imode,' has degeneracy ',deg(imode) 723 write(std_out,'(a,3es16.6)') ' Mode effective charge of this mode =',modez(1,:,imode) 724 725 if (deg(imode) == 2) then 726 727 ! Optimize on the x direction 728 write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1) 729 if (abs(modez(1,1,imode)) > tol8) then 730 theta = atan(-modez(1,1,imode+1)/modez(1,1,imode)) 731 vec(:,1) = displ(1,:,imode) 732 vec(:,2) = displ(1,:,imode+1) 733 displ(1,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2) 734 displ(1,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2) 735 end if 736 737 else if (deg(imode) == 3) then 738 739 write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1) 740 write(std_out,'(a,3es16.6)') ' Mode effective charge of next-next mode =',modez(1,:,imode+2) 741 742 ! Before mixing them, select the mode-effective charge vectors as being predominently "x", "y" or "z" type. 743 if(abs(modez(1,1,imode))>abs(modez(1,2,imode))-tol12 .and. & 744 & abs(modez(1,1,imode))>abs(modez(1,3,imode))-tol12) then 745 imodex=imode 746 if(abs(modez(1,2,imode+1))>abs(modez(1,3,imode+1))-tol12)then 747 imodey=imode+1 ; imodez=imode+2 748 else 749 imodez=imode+1 ; imodey=imode+2 750 end if 751 else if(abs(modez(1,2,imode))>abs(modez(1,1,imode))-tol12 .and. & 752 & abs(modez(1,2,imode))>abs(modez(1,3,imode))-tol12) then 753 imodey=imode 754 if(abs(modez(1,1,imode+1))>abs(modez(1,3,imode+1))-tol12)then 755 imodex=imode+1 ; imodez=imode+2 756 else 757 imodez=imode+1 ; imodex=imode+2 758 end if 759 else 760 imodez=imode 761 if(abs(modez(1,1,imode+1))>abs(modez(1,2,imode+1))-tol12)then 762 imodex=imode+1 ; imodey=imode+2 763 else 764 imodey=imode+1 ; imodex=imode+2 765 end if 766 end if 767 vec(:,1)=displ(1,:,imodex) 768 vec(:,2)=displ(1,:,imodey) 769 vec(:,3)=displ(1,:,imodez) 770 zvec(:,1)=modez(1,:,imodex) 771 zvec(:,2)=modez(1,:,imodey) 772 zvec(:,3)=modez(1,:,imodez) 773 774 775 ! Optimize along x : does the first vector has a component along x ? 776 if (abs(zvec(1,1)) > tol8) then 777 ! Optimize on the (1,2) pair of modes along x 778 theta = atan(-zvec(1,2)/zvec(1,1)) 779 zvect(:,:)=zvec(:,:) 780 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,2) 781 zvec(:,2) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,2) 782 vect(:,:)=vec(:,:) 783 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,2) 784 vec(:,2) = sin(theta)*vect(:,1) + cos(theta)*vect(:,2) 785 ! Optimize on the (1,3) pair of modes along x 786 theta = atan(-zvec(1,3)/zvec(1,1)) 787 zvect(:,:)=zvec(:,:) 788 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3) 789 zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3) 790 vect(:,:)=vec(:,:) 791 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3) 792 vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3) 793 if (abs(zvec(2,2)) > tol8) then 794 ! Optimize on the (2,3) pair of modes along y 795 theta = atan(-zvec(2,3)/zvec(2,2)) 796 zvect(:,:)=zvec(:,:) 797 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 798 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 799 vect(:,:)=vec(:,:) 800 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 801 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 802 end if 803 ! Likely, the remaining is not needed ... because the vectors have been ordered in x, y, and z major component ... 804 ! Optimize along x : does the second vector has a component along x ? 805 else if(abs(zvec(1,2)) > tol8) then 806 ! Optimize on the (2,3) pair of modes along x 807 theta = atan(-zvec(1,3)/zvec(1,2)) 808 zvect(:,:)=zvec(:,:) 809 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 810 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 811 vect(:,:)=vec(:,:) 812 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 813 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 814 ! Optimize on the (1,3) pair of modes along y 815 if (abs(zvec(2,1)) > tol8) then 816 theta = atan(-zvec(2,3)/zvec(2,1)) 817 zvect(:,:)=zvec(:,:) 818 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3) 819 zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3) 820 vect(:,:)=vec(:,:) 821 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3) 822 vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3) 823 end if 824 ! We are left with the pair of vectors (2,3) 825 else if (abs(zvec(2,2)) > tol8) then 826 ! Optimize on the (2,3) pair of modes along y 827 theta = atan(-zvec(2,3)/zvec(2,2)) 828 zvect(:,:)=zvec(:,:) 829 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 830 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 831 vect(:,:)=vec(:,:) 832 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 833 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 834 end if 835 836 displ(1,:,imodex)=vec(:,1) 837 displ(1,:,imodey)=vec(:,2) 838 displ(1,:,imodez)=vec(:,3) 839 840 ! Previous coding, from Marek. Apparently, break the orthogonalization of vectors ... 841 ! do ii = 1,3 842 ! coeff(:) = 0._dp 843 ! if (ii == 1) then 844 ! jj = 2 ; kk = 3 845 ! else if (ii == 2) then 846 ! jj = 1 ; kk = 3 847 ! else 848 ! jj = 1 ; kk = 2 849 ! end if 850 ! coeff(ii) = 1._dp 851 ! c1 = modez(1,jj,imode+ii-1) 852 ! c2 = modez(1,jj,imode+jj-1) 853 ! c3 = modez(1,jj,imode+kk-1) 854 ! c4 = modez(1,kk,imode+ii-1) 855 ! c5 = modez(1,kk,imode+jj-1) 856 ! c6 = modez(1,kk,imode+kk-1) 857 ! dtm = c2*c6 - c3*c5 858 ! if (abs(dtm) > tol8) then 859 ! coeff(jj) = (c3*c4 - c1*c6)/dtm 860 ! coeff(kk) = (c1*c5 - c2*c4)/dtm 861 ! end if 862 ! mod_ = sqrt(1._dp + coeff(jj)*coeff(jj) + coeff(kk)*coeff(kk)) 863 ! coeff(:) = coeff(:)/mod_ 864 ! displ(1,:,imode+ii-1) = coeff(1)*vec(1,:) + coeff(2)*vec(2,:) + & 865 !& coeff(3)*vec(3,:) 866 ! end do 867 868 end if ! if deg mode 869 870 imode = imode + deg(imode) 871 872 end do 873 874 write(std_out,'(a,a)')ch10,' alignph : after modifying the eigenvectors, mode number and mode effective charges :' 875 do imode=1,3*natom 876 do ii=1,2 877 do idir2=1,3 878 modez(ii,idir2,imode)=zero 879 do idir1=1,3 880 do ipert1=1,natom 881 i1=idir1+(ipert1-1)*3 882 modez(ii,idir2,imode)=modez(ii,idir2,imode)+& 883 & displ(ii,i1,imode)*& 884 & d2cart(1,idir1,ipert1,idir2,natom+2)*& 885 & sqrt(amu(typat(ipert1))*amu_emass) 886 end do 887 end do 888 end do 889 end do 890 write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode) 891 end do 892 893 ABI_FREE(deg) 894 ABI_FREE(oscstr) 895 ABI_FREE(modez) 896 ABI_FREE(modezabs) 897 ABI_FREE(vec) 898 ABI_FREE(vect) 899 900 end subroutine alignph
ABINIT/ddb_diel [ Functions ]
NAME
ddb_diel
FUNCTION
Get the frequency-dependent dielectric matrix, as well as the oscillator strengths and mode effective charges, and reflectivities (without damping) See the definitions Eq.(53-54) in PRB55, 10355 (1997) [[cite:Gonze1997a]].
INPUTS
amu(ntypat)=mass of the atoms (atomic mass unit) anaddb_dtset= (derived datatype) contains all the input variables matrix (diagonal in the atoms) displ(2,3*natom,3*natom)= the displacements of atoms in cartesian coordinates. The first index means either the real or the imaginary part, The second index runs on the direction and the atoms displaced The third index runs on the modes. d2cart(2,3,mpert,3,mpert)=dynamical matrix, effective charges, dielectric tensor,... all in cartesian coordinates iout=unit number for outputs lst(3*nph2l)=log. of product of frequencies**2, needed to calculate the generalized Lyddane-Sachs-Teller relation at zero frequency mpert =maximum number of ipert natom=number of atoms in unit cell nph2l=input variable from anaddb_dtset, needed to dimension lst phfrq(3*natom)=phonon frequencies (square root of the dynamical matrix eigenvalues, except if these are negative, and in this case, give minus the square root of the absolute value of the matrix eigenvalues). Hartree units. comm=MPI communicator. ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
OUTPUT
fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes, for different direction of the electric field; dielt_rlx(3,3) relaxed ion (zero frequency) dielectric tensor.
NOTES
1. The phonon frequencies phfrq should correspond to the wavevector at Gamma, without any non-analyticities. 2. Should clean for no imaginary part ... This routine should be used only by one processor. 3. frdiel(3,3,nfreq)= frequency-dependent dielectric tensor mode effective charges for the different eigenmodes, for different direction of the electric field
SOURCE
98 subroutine ddb_diel(Crystal,amu,anaddb_dtset,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,& 99 & iout,lst,mpert,natom,nph2l,phfrq,comm,ncid) 100 101 !Arguments ------------------------------- 102 !scalars 103 integer,intent(in) :: iout,mpert,natom,nph2l,comm,ncid 104 type(crystal_t),intent(in) :: Crystal 105 type(anaddb_dataset_type),intent(in) :: anaddb_dtset 106 real(dp),intent(in) :: lst(nph2l+1) 107 108 !arrays 109 real(dp),intent(in) :: amu(Crystal%ntypat),d2cart(2,3,mpert,3,mpert) 110 real(dp),intent(in) :: phfrq(3*natom),epsinf(3,3) 111 real(dp),intent(inout) :: displ(2,3*natom,3*natom) 112 real(dp),intent(out) :: dielt_rlx(3,3),fact_oscstr(2,3,3*natom) 113 114 !Local variables ------------------------- 115 !scalars 116 integer,parameter :: master=0 117 integer :: dieflag,idir1,idir2,ifreq,ii,imode,iphl2,nfreq 118 integer :: nprocs,my_rank,ncerr 119 real(dp) :: afreq,difffr,eps,q2,ucvol 120 character(len=500) :: message 121 !arrays 122 real(dp) :: qphon(3),refl(3) 123 real(dp),allocatable :: frdiel(:,:,:),modez(:,:,:),oscstr(:,:,:,:),dielt_modedecompo(:,:,:) 124 125 ! ********************************************************************* 126 127 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 128 129 dieflag=anaddb_dtset%dieflag 130 nfreq=anaddb_dtset%nfreq 131 132 ucvol = Crystal%ucvol 133 134 135 ! frdiel(3,3,nfreq)= frequency-dependent dielectric tensor 136 ! modez=mode effective charge 137 ABI_MALLOC(frdiel,(3,3,nfreq)) 138 ABI_MALLOC(oscstr,(2,3,3,3*natom)) 139 ABI_MALLOC(modez,(2,3,3*natom)) 140 ABI_MALLOC(dielt_modedecompo,(3,3,3*natom)) 141 142 !In case only the electronic dielectric constant should be printed 143 if (dieflag==2) then 144 call ddb_diel_elec(iout,epsinf) 145 else 146 ! In case the ionic contribution to the dielectric tensor is asked 147 if (dieflag/=2 .and. nph2l==0) then 148 ! Check if the alignement of phonon modes eigenvector is requested from the input flag alphon; 149 ! useful in case of degenerate modes 150 if (anaddb_dtset%alphon > 0) then 151 write(message, '(3a)' )& 152 ' The alphon input variable is non-zero, will mix the degenerate phonon modes',ch10,& 153 ' in order to align the mode effective charges with the cartesian axes.' 154 call wrtout(std_out,message,'COLL') 155 call wrtout(iout,message,'COLL') 156 call alignph(amu,displ,d2cart,mpert,natom,Crystal%ntypat,phfrq,Crystal%typat) 157 end if ! alignment of the phonon eigenvectors 158 159 ! Compute the mode effective charge and oscillator strength 160 call ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank) 161 162 ! Calculation of epsilon_r (Eq.55 PRB 55, 10355) 163 ! Check the acousticity of the three lowest modes, assuming 164 ! that they are ordered correctly 165 if (abs(phfrq(1))>abs(phfrq(4)))then 166 ! This means that there is at least one mode with truly negative frequency 167 write(message, '(12a,4es16.8)' )& 168 'The lowest mode appears to be a "true" negative mode,',ch10,& 169 'and not an acoustic mode. This precludes the computation',ch10,& 170 'of the frequency-dependent dielectric tensor.',ch10,& 171 'Action : likely there is no action to be taken, although you,',ch10,& 172 'could try to raise your convergence parameters (ecut and k-points).',ch10,& 173 'For your information, here are the four lowest frequencies :',ch10,& 174 (phfrq(ii),ii=1,4) 175 ABI_ERROR(message) 176 end if 177 178 ! Compute the relaxed ion dielectric tensor 179 do idir1=1,3 180 do idir2=1,3 181 ! The electronic contribution to epsilon is added 182 dielt_rlx(idir1,idir2)=epsinf(idir1,idir2) 183 ! calculation of the phonon contribution (ionic) to epsilon 184 do imode=4,3*natom 185 ! Note that the acoustic modes are not included : their 186 ! oscillator strength should be exactly zero 187 ! Also, only the real part of oscstr is taken into account: 188 ! the possible imaginary parts of degenerate modes 189 ! will cancel. 190 dielt_rlx(idir1,idir2)=dielt_rlx(idir1,idir2)+& 191 & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol 192 ! Mode decomposition of epsilon 193 if (dieflag==3)then 194 dielt_modedecompo(idir1,idir2,imode)=oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol 195 endif ! mode decompo of epsilon 196 !DEBUG 197 ! if(idir1==1 .and. idir2==2)then 198 ! write(std_out,'(a,i4,a,3es16.6)')'imode=',imode,' dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)=',& 199 !& dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode) 200 ! endif 201 !ENDDEBUG 202 end do ! imode 203 end do ! idir2 204 end do ! idir1 205 206 ! Print the electronic dielectric tensor 207 call ddb_diel_elec(iout,epsinf) 208 209 ! Print the relaxed ion dielectric tensor 210 write(message,'(a,a)') ch10,' Relaxed ion dielectric tensor' 211 call wrtout(std_out,message,'COLL') 212 call wrtout(iout,message,'COLL') 213 214 do idir1=1,3 215 write(message,'(3f16.8)')(dielt_rlx(idir1,idir2),idir2=1,3) 216 call wrtout(std_out,message,'COLL') 217 call wrtout(iout,message,'COLL') 218 end do 219 call wrtout(iout, " ",'COLL') 220 call wrtout(std_out, " ",'COLL') 221 222 ! Mode decompo of epsilon 223 if (dieflag==3) then 224 write(message,'(a,a,a,a)') ch10,' Mode by mode decomposition of the ionic dielectric tensor',& 225 ch10,' (the electronic contribution is not included)' 226 call wrtout(std_out,message,'COLL') 227 call wrtout(iout,message,'COLL') 228 do imode=4,3*natom 229 write(message,'(a,a,i4,a,es14.6,a,a,3f8.3)') ch10,' Mode number ',imode, ' freq = ',phfrq(imode),' Ha', & 230 ' Mode Z* (|x|, |y|, |z|)', (abs(modez(1,idir1,imode)),idir1=1,3) 231 call wrtout(std_out,message,'COLL') 232 call wrtout(iout,message,'COLL') 233 do idir1=1,3 234 ! write(message,'(a,a,i4)') ch10,' Mode number 2',imode 235 write(message,'(3f16.8)')(dielt_modedecompo(idir1,idir2,imode),idir2=1,3) 236 call wrtout(std_out,message,'COLL') 237 call wrtout(iout,message,'COLL') 238 end do 239 end do 240 endif ! mode decompo of epsilon 241 242 ! write the relaxed ion dielectric tensor to the netcdf 243 #ifdef HAVE_NETCDF 244 if (ncid /= nctk_noid) then 245 ncerr = nctk_def_arrays(ncid, [nctkarr_t("emacro_cart_rlx", "dp", & 246 "number_of_cartesian_directions, number_of_cartesian_directions")],defmode=.True.) 247 NCF_CHECK(ncerr) 248 NCF_CHECK(nctk_set_datamode(ncid)) 249 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart_rlx"), dielt_rlx)) 250 end if 251 #endif 252 253 ! Frequency-dependent dielectric tensor case 254 if(dieflag==1) then 255 write(message,'(4a)') ch10,ch10,' Calculate the freq. dep. dielectric constant',ch10 256 call wrtout(std_out,message,'COLL') 257 write(message,'(3a)') ch10,' Frequency dependent dielectric constant:',ch10 258 call wrtout(iout,message,'COLL') 259 260 ! Check the possibility of asking the frequency-dependent 261 ! dielectric tensor (there should be more than one atom in the unit cell) 262 ! EB: this check is not in any automatic test 263 if(natom==1)then 264 write(message, '(6a)' )& 265 ' ddb_diel : WARNING -',ch10,& 266 ' When there is only one atom in the unit cell',ch10,& 267 ' cell, the dielectric tensor is frequency-independent.'!,& 268 ! ' Consequently, dieflag has been reset to 2 . ' 269 call wrtout(std_out,message,'COLL') 270 call wrtout(iout,message,'COLL') 271 end if 272 273 difffr=zero 274 if(nfreq>1)difffr=(anaddb_dtset%frmax-anaddb_dtset%frmin)/(nfreq-1) 275 276 if (nfreq>10 .and. my_rank == master) then 277 write(iout, '(a,a,a,a,a,a,a,a)' )& 278 ' ddb_diel : the number of frequencies is larger',& 279 ' than 10 => I will consider only',ch10,& 280 ' the three principal directions, assume that the tensor',ch10,& 281 ' is diagonalized, and give dielectric constant and ',ch10,& 282 ' reflectivities.' 283 write(iout, '(a,a)' )& 284 ' Frequency(Hartree) Dielectric constant ',& 285 ' Reflectivity ' 286 write(iout, '(a,a)' )& 287 ' x y z',& 288 ' x y z' 289 end if 290 291 ! Loop on frequencies 292 do ifreq=1,nfreq 293 afreq=anaddb_dtset%frmin+difffr*(ifreq-1) 294 do idir1=1,3 295 do idir2=1,3 296 frdiel(idir1,idir2,ifreq)=epsinf(idir1,idir2) 297 do imode=4,3*natom 298 ! Note that the acoustic modes are not included : their 299 ! oscillator strength should be exactly zero 300 ! Also, only the real part of oscstr is taken into account: 301 ! the possible imaginary parts of degenerate modes 302 ! will cancel. 303 frdiel(idir1,idir2,ifreq)=frdiel(idir1,idir2,ifreq)+& 304 & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2-afreq**2)*four_pi/ucvol 305 end do 306 end do 307 end do 308 309 ! Write all this information (actually, there should be a choice of units for the frequencies ... 310 if (nfreq>10) then 311 do idir1=1,3 312 if(frdiel(idir1,idir1,ifreq)<=zero)then 313 refl(idir1)=one 314 else 315 ! See Gervais and Piriou PRB11,3944(1975) [[cite:Gervais1975]]. 316 refl(idir1)=( (sqrt(frdiel(idir1,idir1,ifreq)) -one) /(sqrt(frdiel(idir1,idir1,ifreq)) +one) )**2 317 end if 318 end do 319 if (my_rank == master) then 320 write(iout, '(7es12.4)' )afreq,(frdiel(idir1,idir1,ifreq),idir1=1,3),(refl(idir1),idir1=1,3) 321 end if 322 323 else 324 if (my_rank == master) then 325 write(iout, '(a,es12.4,a)' )' Full dielectric tensor at frequency',afreq,' Hartree' 326 do idir1=1,3 327 write(iout, '(3f16.8)' ) (frdiel(idir1,idir2,ifreq),idir2=1,3) 328 end do 329 write(iout, '(a)' )' ' 330 end if 331 end if 332 333 end do ! End of the loop on frequencies 334 end if ! End the condition on frequency-dependent dielectric tensor (dieflag=1) 335 end if ! dieflag/=2 .and. nph2l==0 336 end if 337 338 !Calculation of the Lyddane-Sachs-Teller value of the dielectric constant at zero frequency 339 if(nph2l/=0 .and. dieflag/=2) then 340 341 ! Prepare the output 342 write(message, '(a,a,a,a)' ) ch10,& 343 ' Generalized Lyddane-Sachs-Teller relation at zero frequency :',ch10,& 344 ' Direction Dielectric constant' 345 call wrtout(std_out,message,'COLL') 346 call wrtout(iout,message,'COLL') 347 348 ! Examine every wavevector in the phonon list 349 do iphl2=1,anaddb_dtset%nph2l 350 351 qphon(1:3)=anaddb_dtset%qph2l(1:3,iphl2) 352 353 if(abs(qphon(1))>DDB_QTOL .or. abs(qphon(2))>DDB_QTOL .or. abs(qphon(3))>DDB_QTOL)then 354 q2=qphon(1)**2+qphon(2)**2+qphon(3)**2 355 eps=qphon(1)**2*epsinf(1,1)+qphon(2)**2*epsinf(2,2)+& 356 & qphon(3)**2*epsinf(3,3)+ 2* ( qphon(1)*qphon(2)*epsinf(1,2)+& 357 & qphon(1)*qphon(3)*epsinf(1,3)+qphon(2)*qphon(3)*epsinf(2,3)) 358 eps=eps/q2*exp(lst(iphl2)-lst(anaddb_dtset%nph2l+1)) 359 if (my_rank == master) then 360 write(iout, '(3f10.5,f16.8)' )qphon,eps 361 write(std_out,'(3f10.5,f16.8)' )qphon,eps 362 end if 363 end if 364 end do 365 end if ! End of the condition of nph2l does not vanish for Lyddane-Sachs-Teller 366 367 368 ABI_FREE(frdiel) 369 ABI_FREE(modez) 370 ABI_FREE(oscstr) 371 ABI_FREE(dielt_modedecompo) 372 373 end subroutine ddb_diel
ABINIT/ddb_diel_elec [ Functions ]
NAME
ddb_diel_elec
FUNCTION
Print the electronic dielectric constant (clamped ions)
INPUTS
iout=unit number for outputs epsinf(3,3)= epsilon^infty = electronic contribution to the dielectric tensor
OUTPUT
SOURCE
395 subroutine ddb_diel_elec(iout,epsinf) 396 397 !Arguments ------------------------------- 398 !scalars 399 integer,intent(in) :: iout 400 !arrays 401 real(dp),intent(in) :: epsinf(3,3) 402 403 !Local variables ------------------------- 404 !scalars 405 integer :: idir1,idir2 406 character(len=500) :: message 407 !arrays 408 409 write(message, '(a,a)' ) ch10,' Electronic dielectric tensor' 410 call wrtout(std_out,message,'COLL') 411 call wrtout(iout,message,'COLL') 412 413 !Compute the electronic contribution to the dielectric tensor 414 !Needs only the perturbations with E-field from the DDB 415 do idir1=1,3 416 ! do idir2=1,3 417 ! epsinf(idir1,idir2)=d2cart(1,idir1,natom+2,idir2,natom+2) 418 ! end do 419 write(message, '(3f16.8)' )(epsinf(idir1,idir2),idir2=1,3) 420 call wrtout(std_out,message,'COLL') 421 call wrtout(iout,message,'COLL') 422 end do 423 call wrtout(iout, " ",'COLL') 424 call wrtout(std_out, " ",'COLL') 425 426 end subroutine ddb_diel_elec
ABINIT/ddb_oscstr [ Functions ]
NAME
ddb_oscstr
FUNCTION
Compute the oscillator strength and the mode effective charge
INPUTS
displ(2,3*natom,3*natom)= the displacements of atoms in cartesian coordinates. The first index means either the real or the imaginary part, The second index runs on the direction and the atoms displaced The third index runs on the modes. d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates iout=unit number for outputs mpert =maximum number of ipert natom=number of atoms in unit cell phfrq(3*natom)=phonon frequencies (square root of the dynamical matrix eigenvalues, except if these are negative, and in this case, give minus the square root of the absolute value of the matrix eigenvalues). Hartree units.
OUTPUT
fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes, for different direction of the electric field. modez(2,3,3*natom)=mode effective charges for the different eigenmodes, for different directions of the electric field, following the definition Eq.(53) in PRB55, 10355 (1997) [[cite:Gonze1997a]] oscstr(2,3,3,3*natom)=oscillator strengths, following the definition Eq.(54) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
SOURCE
467 subroutine ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank) 468 469 !Arguments ------------------------------- 470 !scalars 471 integer,intent(in) :: iout,mpert,natom,ncid,my_rank 472 !arrays 473 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert) 474 real(dp),intent(in) :: phfrq(3*natom) 475 real(dp),intent(inout) :: displ(2,3*natom,3*natom) 476 real(dp),intent(out) :: fact_oscstr(2,3,3*natom),oscstr(2,3,3,3*natom),modez(2,3,3*natom) 477 478 !Local variables ------------------------- 479 !scalars 480 integer,parameter :: master=0 481 integer :: i1,idir1,idir2,imode,ipert1 482 integer :: ncerr 483 real(dp) :: usquare 484 logical :: t_degenerate 485 !arrays 486 character(len=1),allocatable :: metacharacter(:) 487 488 ! ********************************************************************* 489 490 ! Get the factors of the oscillator strength, and the mode effective charge for each mode 491 do imode=1,3*natom 492 usquare=zero 493 do idir1=1,3 494 do ipert1=1,natom 495 i1=idir1+(ipert1-1)*3 496 usquare=usquare+displ(1,i1,imode)*displ(1,i1,imode)+displ(2,i1,imode)*displ(2,i1,imode) 497 end do 498 end do 499 do idir2=1,3 500 fact_oscstr(:,idir2,imode)=zero 501 modez(:,idir2,imode)=zero 502 do idir1=1,3 503 do ipert1=1,natom 504 i1=idir1+(ipert1-1)*3 505 fact_oscstr(:,idir2,imode)=fact_oscstr(:,idir2,imode)+& 506 & displ(:,i1,imode)*d2cart(1,idir1,ipert1,idir2,natom+2) 507 modez(:,idir2,imode)=modez(:,idir2,imode)+& 508 & displ(:,i1,imode)*& 509 & d2cart(1,idir1,ipert1,idir2,natom+2)/sqrt(usquare) 510 end do 511 end do 512 end do 513 end do 514 515 ! Examine the degeneracy of each mode. The portability of the echo of the mode effective 516 ! charges and oscillator strengths is very hard to guarantee. On the contrary, 517 ! the scalar reductions of these quantities are OK. 518 ABI_MALLOC(metacharacter,(3*natom)) 519 do imode=1,3*natom 520 ! The degenerate modes are not portable 521 t_degenerate=.false. 522 if(imode>1)then 523 if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true. 524 end if 525 if(imode<3*natom)then 526 if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true. 527 end if 528 metacharacter(imode)=';' 529 if(t_degenerate)metacharacter(imode)='-' 530 end do 531 532 if (my_rank == master) then 533 ! Write the mode effective charge for each mode 534 write(iout, '(a)' )' ' 535 write(iout, '(a)' )' Mode effective charges ' 536 write(iout, '(a)' )' Mode number. x y z length ' 537 do imode=1,3*natom 538 write(iout, '(a,i6,a,4f13.6)' )metacharacter(imode),imode,' ',(modez(1,idir1,imode),idir1=1,3),& 539 & (sqrt(modez(1,1,imode)**2+modez(1,2,imode)**2+modez(1,3,imode)**2)) 540 end do 541 end if ! master 542 543 ! Get the oscillator strengths 544 do imode=1,3*natom 545 do idir1=1,3 546 do idir2=1,3 547 oscstr(1,idir1,idir2,imode)= & 548 & fact_oscstr(1,idir1,imode)*fact_oscstr(1,idir2,imode) +& 549 & fact_oscstr(2,idir1,imode)*fact_oscstr(2,idir2,imode) 550 if(abs(oscstr(1,idir1,idir2,imode))<tol14)oscstr(1,idir1,idir2,imode)=zero 551 552 !DEBUG 553 ! if(idir1==1 .and. idir2==2)then 554 ! write(std_out,'(a,i4,a,5es16.6)')'imode=',imode,& 555 !& ' oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)=',& 556 !& oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode) 557 ! endif 558 !ENDDEBUG 559 560 oscstr(2,idir1,idir2,imode)= & 561 & fact_oscstr(1,idir1,imode)*fact_oscstr(2,idir2,imode) -& 562 & fact_oscstr(2,idir1,imode)*fact_oscstr(1,idir2,imode) 563 if(abs(oscstr(2,idir1,idir2,imode))<tol14)oscstr(2,idir1,idir2,imode)=zero 564 end do 565 end do 566 end do 567 568 if (my_rank == master) then 569 ! Write the oscillator strength for each mode 570 write(iout, '(a)' )' ' 571 write(iout, '(a)' )' Oscillator strengths (in a.u. ; 1 a.u.=253.2638413 m3/s2). Set to zero if abs()<tol14.' 572 write(iout, '(a)' )' Mode number. xx yy zz xy xz yz trace' 573 do imode=1,3*natom 574 write(iout, '(a,i4,a,7es12.4)' )& 575 & metacharacter(imode),imode,' Real ',(oscstr(1,idir1,idir1,imode),idir1=1,3),& 576 & oscstr(1,1,2,imode), oscstr(1,1,3,imode),oscstr(1,2,3,imode),& 577 & ((oscstr(1,1,1,imode)+oscstr(1,2,2,imode)+oscstr(1,3,3,imode))) 578 write(iout, '(a,a,6es12.4)' )& 579 & metacharacter(imode),' Imag ',(oscstr(2,idir1,idir1,imode),idir1=1,3),& 580 & oscstr(2,1,2,imode),oscstr(2,1,3,imode),oscstr(2,2,3,imode) 581 end do 582 583 ! write the oscillator strength to the netcdf 584 #ifdef HAVE_NETCDF 585 if (ncid /= nctk_noid) then 586 ncerr = nctk_def_arrays(ncid, [nctkarr_t("oscillator_strength", "dp", & 587 "complex, number_of_cartesian_directions, number_of_cartesian_directions, number_of_phonon_modes")], & 588 defmode=.True.) 589 NCF_CHECK(ncerr) 590 NCF_CHECK(nctk_set_datamode(ncid)) 591 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "oscillator_strength"), oscstr)) 592 end if 593 #endif 594 end if 595 596 ABI_FREE(metacharacter) 597 598 end subroutine ddb_oscstr
ABINIT/m_ddb_diel [ Modules ]
NAME
m_ddb_diel
FUNCTION
This module provides routines for the calculation of the dielectric constant (anaddb)
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG,XW, MVeithen, EB) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_ddb_diel 23 24 use defs_basis 25 use m_errors 26 use m_xmpi 27 use m_abicore 28 use m_ddb 29 use m_nctk 30 #ifdef HAVE_NETCDF 31 use netcdf 32 #endif 33 34 use m_anaddb_dataset, only : anaddb_dataset_type 35 use m_crystal, only : crystal_t 36 37 implicit none 38 39 private