TABLE OF CONTENTS


ABINIT/alignph [ Functions ]

[ Top ] [ 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.

PARENTS

      ddb_diel

CHILDREN

SOURCE

566 subroutine alignph(amu,displ,d2cart,mpert,natom,ntypat,phfrq,typat)
567 
568 
569 !This section has been created automatically by the script Abilint (TD).
570 !Do not modify the following lines by hand.
571 #undef ABI_FUNC
572 #define ABI_FUNC 'alignph'
573 !End of the abilint section
574 
575  implicit none
576 
577 !Arguments -------------------------------
578 !scalars
579  integer,intent(in) :: mpert,natom,ntypat
580 !arrays
581  integer,intent(in) :: typat(natom)
582  real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),phfrq(3*natom)
583  real(dp),intent(inout) :: displ(2,3*natom,3*natom)
584 
585 !Local variables -------------------------
586 !scalars
587  integer :: i1,idir1,idir2,ii,imode,imodex,imodey,imodez,ipert1
588  real(dp) :: theta
589 !arrays
590  integer,allocatable :: deg(:)
591  real(dp) :: zvec(3,3),zvect(3,3)
592  real(dp),allocatable :: modez(:,:,:),modezabs(:),oscstr(:,:,:),vec(:,:),vect(:,:)
593 
594 ! *********************************************************************
595 
596 !Get the oscillator strength and mode effective charge for each mode
597  ABI_ALLOCATE(oscstr,(2,3,3*natom))
598  ABI_ALLOCATE(modez,(2,3,3*natom))
599  ABI_ALLOCATE(modezabs,(3*natom))
600  ABI_ALLOCATE(vec,(3*natom,3))
601  ABI_ALLOCATE(vect,(3*natom,3))
602  ABI_ALLOCATE(deg,(3*natom))
603 
604  write(std_out,'(a,a)')ch10,' alignph : before modifying the eigenvectors, mode number and mode effective charges :'
605  do imode=1,3*natom
606    modezabs(imode)=zero
607    do ii=1,2
608      do idir2=1,3
609        oscstr(ii,idir2,imode)=zero
610        modez(ii,idir2,imode)=zero
611        do idir1=1,3
612          do ipert1=1,natom
613            i1=idir1+(ipert1-1)*3
614            oscstr(ii,idir2,imode)=oscstr(ii,idir2,imode)+&
615 &           displ(ii,i1,imode)*&
616 &           d2cart(1,idir1,ipert1,idir2,natom+2)
617            modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
618 &           displ(ii,i1,imode)*&
619 &           d2cart(1,idir1,ipert1,idir2,natom+2)*&
620 &           sqrt(amu(typat(ipert1))*amu_emass)
621          end do
622        end do
623        if(abs(modez(ii,idir2,imode))>modezabs(imode))modezabs(imode)=abs(modez(ii,idir2,imode))
624      end do
625    end do
626    write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
627  end do
628 
629 !Find degenerate modes with non-zero mode effective charge
630  imode = 0
631  do while (imode < 3*natom)
632    imode = imode + 1
633    if (imode == 3*natom) then
634      deg(imode) = 1
635    else if (abs(phfrq(imode) - phfrq(imode+1)) > tol6 .or. modezabs(imode)<tol8 .or. modezabs(imode+1)<tol8) then
636 !    Differ by phonon frequency or zero mode effective charge
637      deg(imode) = 1
638    else
639      deg(imode) = 2
640      if (imode < 3*natom - 1) then
641        if (abs(phfrq(imode) - phfrq(imode+2)) < tol6 .and. modezabs(imode+2)>tol8) then
642          deg(imode) = 3
643          imode = imode + 1
644        end if
645      end if
646      imode = imode + 1
647    end if
648  end do
649 
650 
651 !In case of a degenerate mode, with non-zero mode effective charge, align the mode effective charge vector along
652 !the axes of the cartesian frame
653  imode = 1
654  do while (imode <= 3*natom)
655 
656    write(std_out,'(a,a,i4,a,i2)')ch10,' Mode number ',imode,' has degeneracy ',deg(imode)
657    write(std_out,'(a,3es16.6)') ' Mode effective charge of this mode =',modez(1,:,imode)
658 
659    if (deg(imode) == 2) then
660 
661 !    Optimize on the x direction
662      write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
663      if (abs(modez(1,1,imode)) > tol8) then
664        theta = atan(-modez(1,1,imode+1)/modez(1,1,imode))
665        vec(:,1) = displ(1,:,imode)
666        vec(:,2) = displ(1,:,imode+1)
667        displ(1,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2)
668        displ(1,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2)
669      end if
670 
671    else if (deg(imode) == 3) then
672 
673      write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
674      write(std_out,'(a,3es16.6)') ' Mode effective charge of next-next mode =',modez(1,:,imode+2)
675 
676 !    Before mixing them, select the mode-effective charge vectors as being predominently "x", "y" or "z" type.
677      if(abs(modez(1,1,imode))>abs(modez(1,2,imode))-tol12 .and. &
678 &     abs(modez(1,1,imode))>abs(modez(1,3,imode))-tol12) then
679        imodex=imode
680        if(abs(modez(1,2,imode+1))>abs(modez(1,3,imode+1))-tol12)then
681          imodey=imode+1 ; imodez=imode+2
682        else
683          imodez=imode+1 ; imodey=imode+2
684        end if
685      else if(abs(modez(1,2,imode))>abs(modez(1,1,imode))-tol12 .and. &
686 &       abs(modez(1,2,imode))>abs(modez(1,3,imode))-tol12) then
687        imodey=imode
688        if(abs(modez(1,1,imode+1))>abs(modez(1,3,imode+1))-tol12)then
689          imodex=imode+1 ; imodez=imode+2
690        else
691          imodez=imode+1 ; imodex=imode+2
692        end if
693      else
694        imodez=imode
695        if(abs(modez(1,1,imode+1))>abs(modez(1,2,imode+1))-tol12)then
696          imodex=imode+1 ; imodey=imode+2
697        else
698          imodey=imode+1 ; imodex=imode+2
699        end if
700      end if
701      vec(:,1)=displ(1,:,imodex)
702      vec(:,2)=displ(1,:,imodey)
703      vec(:,3)=displ(1,:,imodez)
704      zvec(:,1)=modez(1,:,imodex)
705      zvec(:,2)=modez(1,:,imodey)
706      zvec(:,3)=modez(1,:,imodez)
707 
708 
709 !    Optimize along x : does the first vector has a component along x ?
710      if (abs(zvec(1,1)) > tol8) then
711 !      Optimize on the (1,2) pair of modes along x
712        theta = atan(-zvec(1,2)/zvec(1,1))
713        zvect(:,:)=zvec(:,:)
714        zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,2)
715        zvec(:,2) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,2)
716        vect(:,:)=vec(:,:)
717        vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,2)
718        vec(:,2) = sin(theta)*vect(:,1) + cos(theta)*vect(:,2)
719 !      Optimize on the (1,3) pair of modes along x
720        theta = atan(-zvec(1,3)/zvec(1,1))
721        zvect(:,:)=zvec(:,:)
722        zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
723        zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
724        vect(:,:)=vec(:,:)
725        vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
726        vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
727        if (abs(zvec(2,2)) > tol8) then
728 !        Optimize on the (2,3) pair of modes along y
729          theta = atan(-zvec(2,3)/zvec(2,2))
730          zvect(:,:)=zvec(:,:)
731          zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
732          zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
733          vect(:,:)=vec(:,:)
734          vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
735          vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
736        end if
737 !    Likely, the remaining is not needed ... because the vectors have been ordered in x, y, and z major component ...
738 !    Optimize along x : does the second vector has a component along x ?
739      else if(abs(zvec(1,2)) > tol8) then
740 !      Optimize on the (2,3) pair of modes along x
741        theta = atan(-zvec(1,3)/zvec(1,2))
742        zvect(:,:)=zvec(:,:)
743        zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
744        zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
745        vect(:,:)=vec(:,:)
746        vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
747        vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
748 !      Optimize on the (1,3) pair of modes along y
749        if (abs(zvec(2,1)) > tol8) then
750          theta = atan(-zvec(2,3)/zvec(2,1))
751          zvect(:,:)=zvec(:,:)
752          zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
753          zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
754          vect(:,:)=vec(:,:)
755          vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
756          vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
757        end if
758 !    We are left with the pair of vectors (2,3)
759      else if (abs(zvec(2,2)) > tol8) then
760 !      Optimize on the (2,3) pair of modes along y
761        theta = atan(-zvec(2,3)/zvec(2,2))
762        zvect(:,:)=zvec(:,:)
763        zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
764        zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
765        vect(:,:)=vec(:,:)
766        vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
767        vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
768      end if
769 
770      displ(1,:,imodex)=vec(:,1)
771      displ(1,:,imodey)=vec(:,2)
772      displ(1,:,imodez)=vec(:,3)
773 
774 !    Previous coding, from Marek. Apparently, break the orthogonalization of vectors ...
775 !    do ii = 1,3
776 !      coeff(:) = 0._dp
777 !      if (ii == 1) then
778 !        jj = 2 ; kk = 3
779 !      else if (ii == 2) then
780 !        jj = 1 ; kk = 3
781 !      else
782 !        jj = 1 ; kk = 2
783 !      end if
784 !      coeff(ii) = 1._dp
785 !      c1 = modez(1,jj,imode+ii-1)
786 !      c2 = modez(1,jj,imode+jj-1)
787 !      c3 = modez(1,jj,imode+kk-1)
788 !      c4 = modez(1,kk,imode+ii-1)
789 !      c5 = modez(1,kk,imode+jj-1)
790 !      c6 = modez(1,kk,imode+kk-1)
791 !      dtm = c2*c6 - c3*c5
792 !      if (abs(dtm) > tol8) then
793 !        coeff(jj) = (c3*c4 - c1*c6)/dtm
794 !        coeff(kk) = (c1*c5 - c2*c4)/dtm
795 !      end if
796 !      mod_ = sqrt(1._dp + coeff(jj)*coeff(jj) + coeff(kk)*coeff(kk))
797 !      coeff(:) = coeff(:)/mod_
798 !      displ(1,:,imode+ii-1) = coeff(1)*vec(1,:) + coeff(2)*vec(2,:) + &
799 !&       coeff(3)*vec(3,:)
800 !    end do
801 
802    end if ! if deg mode
803 
804    imode = imode + deg(imode)
805 
806  end do
807 
808  write(std_out,'(a,a)')ch10,' alignph : after modifying the eigenvectors, mode number and mode effective charges :'
809  do imode=1,3*natom
810    do ii=1,2
811      do idir2=1,3
812        modez(ii,idir2,imode)=zero
813        do idir1=1,3
814          do ipert1=1,natom
815            i1=idir1+(ipert1-1)*3
816            modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
817 &           displ(ii,i1,imode)*&
818 &           d2cart(1,idir1,ipert1,idir2,natom+2)*&
819 &           sqrt(amu(typat(ipert1))*amu_emass)
820          end do
821        end do
822      end do
823    end do
824    write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
825  end do
826 
827  ABI_DEALLOCATE(deg)
828  ABI_DEALLOCATE(oscstr)
829  ABI_DEALLOCATE(modez)
830  ABI_DEALLOCATE(modezabs)
831  ABI_DEALLOCATE(vec)
832  ABI_DEALLOCATE(vect)
833 
834 end subroutine alignph

ABINIT/ddb_diel [ Functions ]

[ Top ] [ 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

PARENTS

      anaddb

CHILDREN

      alignph,wrtout

SOURCE

108 subroutine ddb_diel(Crystal,amu,anaddb_dtset,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
109 & iout,lst,mpert,natom,nph2l,phfrq,comm,ncid)
110 
111 
112 !This section has been created automatically by the script Abilint (TD).
113 !Do not modify the following lines by hand.
114 #undef ABI_FUNC
115 #define ABI_FUNC 'ddb_diel'
116 !End of the abilint section
117 
118  implicit none
119 
120 !Arguments -------------------------------
121 !scalars
122  integer,intent(in) :: iout,mpert,natom,nph2l,comm,ncid
123  type(crystal_t),intent(in) :: Crystal
124  type(anaddb_dataset_type),intent(in) :: anaddb_dtset
125 !arrays
126  real(dp),intent(in) :: amu(Crystal%ntypat),d2cart(2,3,mpert,3,mpert),lst(nph2l)
127  real(dp),intent(in) :: phfrq(3*natom)
128  real(dp),intent(inout) :: displ(2,3*natom,3*natom)
129  real(dp),intent(out) :: dielt_rlx(3,3),epsinf(3,3),fact_oscstr(2,3,3*natom)
130 
131 !Local variables -------------------------
132 !scalars
133  integer,parameter :: master=0
134  integer :: dieflag,i1,idir1,idir2,ifreq,ii,imode,ipert1,iphl2,nfreq
135  integer :: nprocs,my_rank,ncerr
136  real(dp) :: afreq,difffr,eps,lst0,q2,usquare,ucvol
137  character(len=500) :: message
138  logical :: t_degenerate
139 !arrays
140  real(dp) :: qphon(3),refl(3)
141  real(dp),allocatable :: frdiel(:,:,:),modez(:,:,:),oscstr(:,:,:,:)
142  character(len=1),allocatable :: metacharacter(:)
143 
144 ! *********************************************************************
145 
146  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
147 
148  dieflag=anaddb_dtset%dieflag
149  nfreq=anaddb_dtset%nfreq
150 
151  ucvol = Crystal%ucvol
152 
153 !Check the possibility of asking the frequency-dependent
154 !dielectric tensor (there should be more than one atom in the unit cell)
155  if(natom==1)then
156    write(message, '(6a)' )&
157 &   ' ddb_diel : WARNING -',ch10,&
158 &   '  When there is only one atom in the unit cell',ch10,&
159 &   '  cell, the dielectric tensor is frequency-independent.',&
160 &   '  Consequently, dieflag has been reset to 2 . '
161    call wrtout(std_out,message,'COLL')
162    call wrtout(iout,message,'COLL')
163    dieflag=2
164  end if
165 
166 !frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
167 !modez(2,3,3*natom)=mode effective charges for the different eigenmodes,
168 !for different directions of the electric field, following
169 !the definition Eq.(53) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
170 !fact_oscstr(2,3,3*natom)=factors of the oscillator strengths
171 !for the different eigenmodes,
172 !for different direction of the electric field
173  ABI_ALLOCATE(frdiel,(3,3,nfreq))
174  ABI_ALLOCATE(modez,(2,3,3*natom))
175 !oscstr(2,3,3,3*natom)=oscillator strengths, following
176 !the definition Eq.(54) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
177  ABI_ALLOCATE(oscstr,(2,3,3,3*natom))
178 
179 ! write(std_out,'(a)')' Enter ddb_diel : displ ='
180 ! do imode=1,3*natom
181 !   write(std_out,'(a,i4,a,12es16.6)')'imode=',imode,' displ(:,:,imode)',displ(:,:,imode)
182 ! end do
183 
184 !In case the frequency-dependent dielectric tensor is asked
185  if(dieflag==1 .or. dieflag==3 .or. dieflag==4)then
186 
187    if (anaddb_dtset%alphon > 0) then
188      write(message, '(3a)' )&
189 &     ' The alphon input variable is non-zero, will mix the degenerate phonon modes',ch10,&
190 &     ' in order to align the mode effective charges with the cartesian axes.'
191      call wrtout(std_out,message,'COLL')
192      call wrtout(iout,message,'COLL')
193      call alignph(amu,displ,d2cart,mpert,natom,Crystal%ntypat,phfrq,Crystal%typat)
194    end if
195 
196 !  Get the factors of the oscillator strength, and the mode effective charge for each mode
197    do imode=1,3*natom
198      usquare=zero
199      do idir1=1,3
200        do ipert1=1,natom
201          i1=idir1+(ipert1-1)*3
202          usquare=usquare+displ(1,i1,imode)*displ(1,i1,imode)+displ(2,i1,imode)*displ(2,i1,imode)
203        end do
204      end do
205      do idir2=1,3
206        fact_oscstr(:,idir2,imode)=zero
207        modez(:,idir2,imode)=zero
208        do idir1=1,3
209          do ipert1=1,natom
210            i1=idir1+(ipert1-1)*3
211            fact_oscstr(:,idir2,imode)=fact_oscstr(:,idir2,imode)+&
212 &           displ(:,i1,imode)*d2cart(1,idir1,ipert1,idir2,natom+2)
213            modez(:,idir2,imode)=modez(:,idir2,imode)+&
214 &           displ(:,i1,imode)*&
215 &           d2cart(1,idir1,ipert1,idir2,natom+2)/sqrt(usquare)
216          end do
217        end do
218      end do
219 
220 !DEBUG
221 !           write(std_out,'(a,i4,a,12es16.6)')'imode=',imode,&
222 !&           ' displ(:,:,imode)',displ(:,:,imode)
223 !           write(std_out,'(a,i4,a,6es16.6)')'imode=',imode,&
224 !&           ' fact_oscstr(:,:,imode)=',fact_oscstr(:,:,imode)
225 !ENDDEBUG
226 
227    end do
228 
229 !  Examine the degeneracy of each mode. The portability of the echo of the mode effective
230 !  charges and oscillator strengths is very hard to guarantee. On the contrary,
231 !  the scalar reductions of these quantities are OK.
232    ABI_ALLOCATE(metacharacter,(3*natom))
233    do imode=1,3*natom
234 !    The degenerate modes are not portable
235      t_degenerate=.false.
236      if(imode>1)then
237        if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
238      end if
239      if(imode<3*natom)then
240        if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
241      end if
242      metacharacter(imode)=';'
243      if(t_degenerate)metacharacter(imode)='-'
244    end do
245 
246    if (my_rank == master) then
247      !  Write the mode effective charge for each mode
248      write(iout, '(a)' )'  '
249      write(iout, '(a)' )' Mode effective charges '
250      write(iout, '(a)' )' Mode number.     x               y               z '
251      do imode=1,3*natom
252        write(iout, '(a,i4,3f16.3)' )metacharacter(imode),imode,(modez(1,idir1,imode),idir1=1,3)
253      end do
254 
255      !  Write the mode effective charge length for each mode
256      write(iout, '(a)' )'  '
257      write(iout, '(a)' )' Length of mode effective charge for each phonon mode :'
258      do imode=1,3*natom,5
259        if (3*natom-imode<5) then
260          write(iout, '(1x,5es14.6)') (sqrt(modez(1,1,ii)**2+modez(1,2,ii)**2+modez(1,3,ii)**2),ii=imode,3*natom)
261        else
262          write(iout, '(1x,5es14.6)') (sqrt(modez(1,1,ii)**2+modez(1,2,ii)**2+modez(1,3,ii)**2),ii=imode,imode+4)
263        end if
264      end do
265    end if ! master
266 
267 !  Get the oscillator strengths
268    do imode=1,3*natom
269      do idir1=1,3
270        do idir2=1,3
271          oscstr(1,idir1,idir2,imode)= &
272 &         fact_oscstr(1,idir1,imode)*fact_oscstr(1,idir2,imode) +&
273 &         fact_oscstr(2,idir1,imode)*fact_oscstr(2,idir2,imode)
274          if(abs(oscstr(1,idir1,idir2,imode))<tol14)oscstr(1,idir1,idir2,imode)=zero
275 
276 !DEBUG
277 !         if(idir1==1 .and. idir2==2)then
278 !           write(std_out,'(a,i4,a,5es16.6)')'imode=',imode,&
279 !&           ' oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)=',&
280 !&            oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)
281 !         endif
282 !ENDDEBUG
283 
284          oscstr(2,idir1,idir2,imode)= &
285 &         fact_oscstr(1,idir1,imode)*fact_oscstr(2,idir2,imode) -&
286 &         fact_oscstr(2,idir1,imode)*fact_oscstr(1,idir2,imode)
287          if(abs(oscstr(2,idir1,idir2,imode))<tol14)oscstr(2,idir1,idir2,imode)=zero
288        end do
289      end do
290    end do
291 
292    if (my_rank == master) then
293      !  Write the oscillator strength for each mode
294      write(iout, '(a)' )'  '
295      write(iout, '(a)' )' Oscillator strengths (in a.u. ; 1 a.u.=253.2638413 m3/s2). Set to zero if abs()<tol14.'
296      write(iout, '(a)' )' Mode number.       xx          yy          zz          xy          xz          yz '
297      do imode=1,3*natom
298        write(iout, '(a,i4,a,6es12.4)' )&
299 &       metacharacter(imode),imode,'     Real  ',(oscstr(1,idir1,idir1,imode),idir1=1,3),&
300 &       oscstr(1,1,2,imode), oscstr(1,1,3,imode),oscstr(1,2,3,imode)
301        write(iout, '(a,a,6es12.4)' )&
302 &       metacharacter(imode),'         Imag  ',(oscstr(2,idir1,idir1,imode),idir1=1,3),&
303 &       oscstr(2,1,2,imode),oscstr(2,1,3,imode),oscstr(2,2,3,imode)
304      end do
305 
306      ! write the oscillator strength to the netcdf
307 #ifdef HAVE_NETCDF
308      if (ncid /= nctk_noid) then
309        ncerr = nctk_def_arrays(ncid, [nctkarr_t("oscillator_strength", "dp", &
310        "complex, number_of_cartesian_directions, number_of_cartesian_directions, number_of_phonon_modes")], &
311        defmode=.True.)
312        NCF_CHECK(ncerr)
313        NCF_CHECK(nctk_set_datamode(ncid))
314        NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "oscillator_strength"), oscstr))
315      end if
316 #endif
317 
318      !  Write the trace of oscillator strength (real part only) for each mode
319      write(iout, '(a)' )'  '
320      write(iout, '(a)' )' Trace of oscillator strength, for each phonon mode :'
321      do imode=1,3*natom,5
322        if (3*natom-imode<5) then
323          write(iout, '(1x,5es14.6)') ((oscstr(1,1,1,ii)+oscstr(1,2,2,ii)+oscstr(1,3,3,ii)),ii=imode,3*natom)
324        else
325          write(iout, '(1x,5es14.6)') ((oscstr(1,1,1,ii)+oscstr(1,2,2,ii)+oscstr(1,3,3,ii)),ii=imode,imode+4)
326        end if
327      end do
328    end if
329 
330    ABI_DEALLOCATE(metacharacter)
331 
332  end if !  end the condition on frequency-dependent dielectric tensor
333 
334 !In case the electronic dielectric tensor is needed
335  if(dieflag==1.or.dieflag==2.or.dieflag==3 .or. dieflag==4)then
336 
337    write(message, '(a,a)' ) ch10,' Electronic dielectric tensor'
338    call wrtout(std_out,message,'COLL')
339    call wrtout(iout,message,'COLL')
340 
341    do idir1=1,3
342      do idir2=1,3
343        epsinf(idir1,idir2)=d2cart(1,idir1,natom+2,idir2,natom+2)
344      end do
345      write(message, '(3f16.8)' )(epsinf(idir1,idir2),idir2=1,3)
346      call wrtout(std_out,message,'COLL')
347      call wrtout(iout,message,'COLL')
348    end do
349    call wrtout(iout, " ",'COLL')
350    call wrtout(std_out, " ",'COLL')
351 
352  end if
353 
354 !Only in case the frequency-dependent dielectric tensor is needed
355  if(dieflag==1 .or. dieflag==3 .or. dieflag==4) then
356 !  Check the acousticity of the three lowest modes, assuming
357 !  that they are ordered correctly
358    if (abs(phfrq(1))>abs(phfrq(4)))then
359 !    This means that there is at least one mode with truly negative frequency
360      write(message, '(12a,4es16.8)' )&
361 &     'The lowest mode appears to be a "true" negative mode,',ch10,&
362 &     'and not an acoustic mode. This precludes the computation',ch10,&
363 &     'of the frequency-dependent dielectric tensor.',ch10,&
364 &     'Action : likely there is no action to be taken, although you,',ch10,&
365 &     'could try to raise your convergence parameters (ecut and k-points).',ch10,&
366 &     'For your information, here are the four lowest frequencies :',ch10,&
367 &     (phfrq(ii),ii=1,4)
368      MSG_ERROR(message)
369    end if
370 
371 !  Extract the relaxed ion dielectric tensor
372    do idir1=1,3
373      do idir2=1,3
374        dielt_rlx(idir1,idir2)=epsinf(idir1,idir2)
375        do imode=4,3*natom
376 !        Note that the acoustic modes are not included : their
377 !        oscillator strength should be exactly zero
378 !        Also, only the real part of oscstr is taken into account:
379 !        the possible imaginary parts of degenerate modes
380 !        will cancel.
381          dielt_rlx(idir1,idir2)=dielt_rlx(idir1,idir2)+&
382 &         oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol
383 !DEBUG
384 !         if(idir1==1 .and. idir2==2)then
385 !           write(std_out,'(a,i4,a,3es16.6)')'imode=',imode,' dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)=',&
386 !&            dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)
387 !         endif
388 !ENDDEBUG
389        end do
390      end do
391    end do
392 
393    write(message,'(a,a)') ch10,' Relaxed ion dielectric tensor'
394    call wrtout(std_out,message,'COLL')
395    call wrtout(iout,message,'COLL')
396 
397    do idir1=1,3
398      write(message,'(3f16.8)')(dielt_rlx(idir1,idir2),idir2=1,3)
399      call wrtout(std_out,message,'COLL')
400      call wrtout(iout,message,'COLL')
401    end do
402    call wrtout(iout, " ",'COLL')
403    call wrtout(std_out, " ",'COLL')
404 
405    ! write the relaxed ion dielectric tensor to the netcdf
406 #ifdef HAVE_NETCDF
407    if (ncid /= nctk_noid) then
408      ncerr = nctk_def_arrays(ncid, [nctkarr_t("emacro_cart_rlx", "dp", &
409      "number_of_cartesian_directions, number_of_cartesian_directions")],defmode=.True.)
410      NCF_CHECK(ncerr)
411      NCF_CHECK(nctk_set_datamode(ncid))
412      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart_rlx"), dielt_rlx))
413    end if
414 #endif
415 
416  end if
417 
418 !Only in case the frequency-dependent dielectric tensor is needed
419  if(dieflag==1) then
420 
421    difffr=zero
422    if(nfreq>1)difffr=(anaddb_dtset%frmax-anaddb_dtset%frmin)/(nfreq-1)
423 
424    if (nfreq>10 .and. my_rank == master) then
425      write(iout, '(a,a,a,a,a,a,a,a)' )&
426 &     ' ddb_diel : the number of frequencies is larger',&
427 &     ' than 10 => I will consider only',ch10,&
428 &     ' the three principal directions, assume that the tensor',ch10,&
429 &     ' is diagonalized, and give dielectric constant and ',ch10,&
430 &     ' reflectivities.'
431      write(iout, '(a,a)' )&
432 &     ' Frequency(Hartree)    Dielectric constant   ',&
433 &     '             Reflectivity    '
434      write(iout, '(a,a)' )&
435 &     '                     x           y          z',&
436 &     '          x        y        z'
437    end if
438 
439 !  Loop on frequencies
440    do ifreq=1,nfreq
441      afreq=anaddb_dtset%frmin+difffr*(ifreq-1)
442      do idir1=1,3
443        do idir2=1,3
444          frdiel(idir1,idir2,ifreq)=epsinf(idir1,idir2)
445          do imode=4,3*natom
446 !          Note that the acoustic modes are not included : their
447 !          oscillator strength should be exactly zero
448 !          Also, only the real part of oscstr is taken into account:
449 !          the possible imaginary parts of degenerate modes
450 !          will cancel.
451            frdiel(idir1,idir2,ifreq)=frdiel(idir1,idir2,ifreq)+&
452 &           oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2-afreq**2)*four_pi/ucvol
453          end do
454        end do
455      end do
456 
457      ! Write all this information (actually, there should be a choice of units for the frequencies ...
458      if (nfreq>10) then
459        do idir1=1,3
460          if(frdiel(idir1,idir1,ifreq)<=zero)then
461            refl(idir1)=one
462          else
463 !          See Gervais and Piriou PRB11,3944(1975) [[cite:Gervais1975]].
464            refl(idir1)=( (sqrt(frdiel(idir1,idir1,ifreq)) -one) /(sqrt(frdiel(idir1,idir1,ifreq)) +one) )**2
465          end if
466        end do
467        if (my_rank == master) then
468          write(iout, '(7es12.4)' )afreq,(frdiel(idir1,idir1,ifreq),idir1=1,3),(refl(idir1),idir1=1,3)
469        end if
470 
471      else
472        if (my_rank == master) then
473          write(iout, '(a,es12.4,a)' )' Full dielectric tensor at frequency',afreq,' Hartree'
474          do idir1=1,3
475            write(iout, '(3es16.8)' ) (frdiel(idir1,idir2,ifreq),idir2=1,3)
476          end do
477          write(iout, '(a)' )' '
478        end if
479      end if
480 
481    end do ! End of the loop on frequencies
482  end if ! End the condition on frequency-dependent dielectric tensor
483 
484 !Calculation of the Lyddane-Sachs-Teller value of the dielectric constant at zero frequency
485  if(anaddb_dtset%nph2l/=0 .and.dieflag==1)then
486 
487 !  Get the log of product of the square of the frequencies without non-analyticities.
488    lst0=zero
489    do imode=4,3*natom
490      lst0=lst0+2*log(phfrq(imode))
491    end do
492 
493 !  Prepare the output
494    write(message, '(a,a,a,a)' ) ch10,&
495 &   ' Generalized Lyddane-Sachs-Teller relation at zero frequency :',ch10,&
496 &   ' Direction                     Dielectric constant'
497    call wrtout(std_out,message,'COLL')
498    call wrtout(iout,message,'COLL')
499 
500 !  Examine every wavevector in the phonon list
501    do iphl2=1,anaddb_dtset%nph2l
502      qphon(1:3)=anaddb_dtset%qph2l(1:3,iphl2)
503      if( abs(qphon(1))>DDB_QTOL .or. abs(qphon(2))>DDB_QTOL .or. abs(qphon(3))>DDB_QTOL)then
504        q2=qphon(1)**2+qphon(2)**2+qphon(3)**2
505        eps=qphon(1)**2*epsinf(1,1)+qphon(2)**2*epsinf(2,2)+&
506 &       qphon(3)**2*epsinf(3,3)+ 2* ( qphon(1)*qphon(2)*epsinf(1,2)+&
507 &       qphon(1)*qphon(3)*epsinf(1,3)+qphon(2)*qphon(3)*epsinf(2,3))
508        eps=eps/q2*exp(lst(iphl2)-lst0)
509        if (my_rank == master) then
510          write(iout, '(3f10.5,es18.8)' )qphon,eps
511          write(std_out,'(3f10.5,es18.8)' )qphon,eps
512        end if
513      end if
514    end do
515  end if ! End of the condition of nph2l does not vanish
516 
517  ABI_DEALLOCATE(frdiel)
518  ABI_DEALLOCATE(modez)
519  ABI_DEALLOCATE(oscstr)
520 
521 end subroutine ddb_diel

ABINIT/m_ddb_diel [ Modules ]

[ Top ] [ Modules ]

NAME

 m_ddb_diel

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_ddb_diel
27 
28  use defs_basis
29  use m_errors
30  use m_xmpi
31  use m_abicore
32  use m_ddb
33  use m_nctk
34 #ifdef HAVE_NETCDF
35  use netcdf
36 #endif
37 
38  use m_anaddb_dataset, only : anaddb_dataset_type
39  use m_crystal,        only : crystal_t
40 
41  implicit none
42 
43  private