TABLE OF CONTENTS


ABINIT/m_matlu [ Modules ]

[ Top ] [ Modules ]

NAME

  m_matlu

FUNCTION

COPYRIGHT

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

INPUTS

OUTPUT

NOTES

  subroutines in this module must never call
   a subroutine of m_oper, m_green, m_self
   in order to avoid circular dependancies

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 
29 #include "abi_common.h"
30 
31 MODULE m_matlu
32 
33  use defs_basis
34  use m_errors
35  use m_abicore
36 
37  use m_hide_lapack,  only : xginv
38  use m_geometry, only : symredcart
39 
40  implicit none
41 
42  private
43 
44  public :: init_matlu
45  public :: inverse_matlu
46  public :: destroy_matlu
47  public :: diff_matlu
48  public :: add_matlu
49  public :: print_matlu
50  public :: sym_matlu
51  public :: copy_matlu
52  public :: gather_matlu
53  public :: zero_matlu
54  public :: trace_matlu
55  public :: diag_matlu
56  public :: rotate_matlu
57  public :: shift_matlu
58  public :: checkdiag_matlu
59  public :: checkreal_matlu
60  public :: prod_matlu
61  public :: conjg_matlu
62  public :: ln_matlu
63  public :: slm2ylm_matlu
64  public :: fac_matlu
65  public :: printplot_matlu
66  public :: identity_matlu
67  public :: magmomforb_matlu
68  public :: magmomfspin_matlu
69  public :: magmomfzeeman_matlu
70  public :: chi_matlu

m_matlu/add_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 add_matlu

FUNCTION

INPUTS

  maltu1 <type(matlu_type)>= density matrix matlu1 in the local orbital basis and related variables
  maltu2 <type(matlu_type)>= density matrix matlu2 in the local orbital basis and related variables
  natom = number of atoms
  sign_matlu2= 1 add matlu1 and matlu2
              -1 substract matlu2 to matlu1

OUTPUT

  maltu3 <type(matlu_type)>= density matrix matlu3, sum/substract matlu1 and matlu2

SOURCE

 982 subroutine add_matlu(matlu1,matlu2,matlu3,natom,sign_matlu2)
 983 
 984  use defs_basis
 985  use m_paw_dmft, only : paw_dmft_type
 986  use m_crystal, only : crystal_t
 987  implicit none
 988 
 989 !Arguments ------------------------------------
 990 !type
 991  integer,intent(in) :: natom,sign_matlu2
 992  type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
 993  type(matlu_type), intent(inout) :: matlu3(natom) !vz_i
 994 
 995 !local variables-------------------------------
 996  integer :: iatom
 997 ! *********************************************************************
 998 
 999  do iatom = 1 , natom
1000    matlu3(iatom)%mat=matlu1(iatom)%mat+float(sign_matlu2)*matlu2(iatom)%mat
1001  enddo ! natom
1002 
1003 end subroutine add_matlu

m_matlu/checkdiag_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 checkdiag_matlu

FUNCTION

 Check that matlu is diagonal in the orbital index with given precision

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.
  tol : precision

SIDE EFFECTS

NOTES

SOURCE

2254  subroutine checkdiag_matlu(matlu,natom,tol,nondiag)
2255  use defs_basis
2256  use defs_wvltypes
2257  use m_crystal, only : crystal_t
2258  implicit none
2259 
2260 !Arguments ------------------------------------
2261 !scalars
2262  real(dp),intent(in) :: tol
2263  integer, intent(in) :: natom
2264  logical, intent(out) :: nondiag
2265 !arrays
2266  type(matlu_type),intent(inout) :: matlu(natom)
2267 !Local variables-------------------------------
2268 !scalars
2269  integer :: iatom,im,im1,ispinor,ispinor1,isppol
2270  integer :: lpawu
2271 !arrays
2272 !************************************************************************
2273  nondiag=.false.
2274  do iatom=1,natom
2275    lpawu=matlu(iatom)%lpawu
2276    if(lpawu.ne.-1) then
2277      do im=1,2*lpawu+1
2278        do im1=1,2*lpawu+1
2279          do isppol=1,matlu(1)%nsppol
2280            do ispinor=1,matlu(1)%nspinor
2281 !              if(im/=im1) write(std_out,*) "im,im1",im,im1,matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor)
2282  !            if(present(nondiag).eqv..false.) then
2283  !              if(im/=im1.and.(abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor))>tol))  then
2284  !                write(message,'(5i5)') im,im1,isppol,ispinor,ispinor
2285  !                call wrtout(std_out,message,'COLL')
2286  !                write(message,'(a,3e16.5)')" checkdiag_matlu: Warning ",matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor),tol
2287  !                call wrtout(std_out,message,'COLL')
2288  !                if(.not.present(opt)) ABI_ERROR("not present(opt)")
2289  !                if(matlu(1)%nspinor==1) ABI_ERROR("matlu%nspinor==1")
2290  !              endif
2291 !             endif
2292              do ispinor1=1,matlu(1)%nspinor
2293  !              if(present(nondiag)) then
2294                  if((im/=im1.or.ispinor/=ispinor1)&
2295 &                         .and.abs(real(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)))>tol) then
2296                    nondiag=.true.
2297                   ! write(6,*) "NONDIAG", matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)
2298                  endif
2299                !if(ispinor/=ispinor1.and.(abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))>tol))  then
2300                !  write(message,'(a,3e16.5)')" checkdiag_matlu :i Warning ",matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1),tol
2301                !  call wrtout(std_out,message,'COLL')
2302                !  write(message,'(5i5)') im,im1,isppol,ispinor,ispinor
2303                !  call wrtout(std_out,message,'COLL')
2304                !  if(matlu(1)%nspinor==1) ABI_ERROR("matlu%nspinor==1")
2305                !endif
2306              enddo
2307            enddo ! ispinor
2308          enddo ! isppol
2309        enddo ! im1
2310      enddo ! im
2311    endif ! lpawu
2312  enddo ! iatom
2313 
2314  end subroutine checkdiag_matlu

m_matlu/checkreal_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 checkreal_matlu

FUNCTION

 Check that matlu is real in the orbital index with given precision

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.
  tol : precision

SIDE EFFECTS

NOTES

SOURCE

2153  subroutine checkreal_matlu(matlu,natom,tol)
2154  use defs_basis
2155  use defs_wvltypes
2156  use m_crystal, only : crystal_t
2157  implicit none
2158 
2159 !Arguments ------------------------------------
2160 !scalars
2161  real(dp),intent(in) :: tol
2162  integer, intent(in) :: natom
2163 !arrays
2164  type(matlu_type),intent(inout) :: matlu(natom)
2165 !Local variables-------------------------------
2166 !scalars
2167  integer :: iatom,im,im1,ispinor,ispinor1,isppol
2168  integer :: lpawu
2169  character(len=500) :: message
2170  real(dp) :: maximag,maxoffdiag,maximagdiag
2171 !arrays
2172 !************************************************************************
2173  maximag=zero
2174  maximagdiag=zero
2175  maxoffdiag=zero
2176  do iatom=1,natom
2177    lpawu=matlu(iatom)%lpawu
2178    if(lpawu.ne.-1) then
2179      do im=1,2*lpawu+1
2180        do im1=1,2*lpawu+1
2181          do isppol=1,matlu(1)%nsppol
2182            do ispinor=1,matlu(1)%nspinor
2183              do ispinor1=1,matlu(1)%nspinor
2184                if(abs(aimag(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)))>maximag) then
2185                  maximag=abs(aimag(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)))
2186                  if(im==im1.and.ispinor==ispinor1) maximagdiag=abs(aimag(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)))
2187                endif
2188                if((im/=im1.or.ispinor/=ispinor1).and.abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))>maxoffdiag) then
2189                  maxoffdiag=abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
2190                endif
2191              enddo
2192            enddo ! ispinor
2193          enddo ! isppol
2194        enddo ! im1
2195      enddo ! im
2196    endif ! lpawu
2197  enddo ! iatom
2198  if (maximagdiag>tol) then
2199    write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2200 &   ' Diagonal part of the occupation matrix is complex: the imaginary part ',&
2201 &     maximagdiag,' is larger than',tol,ch10  &
2202 &    , "The calculation cannot handle it : check that your calculation is meaningfull"
2203    ABI_ERROR(message)
2204  endif
2205  if (maximag>tol) then
2206    write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2207 &   ' Off diag occupation matrix is complex: the imaginary part ',maximag,' is larger than',tol,ch10&
2208     , "Check that your calculation is meaningfull"
2209    ABI_WARNING(message)
2210  else
2211    write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2212 &   ' Occupation matrix is real: the imaginary part ',maximag,' is lower than',tol
2213    ABI_COMMENT(message)
2214  endif
2215  if (maxoffdiag>tol) then
2216    write(message,'(3x,2a,e12.4,a,e12.4,6a)') ch10,&
2217 &   ' Occupation matrix is non diagonal : the maximum off-diag part ',maxoffdiag,' is larger than',tol,ch10&
2218 &    , "The corresponding non diagonal elements will be neglected in the Weiss/Hybridization functions",ch10&
2219 &    , "(Except if dmft_solv=8,9 where these elements are taken into accounts)",ch10&
2220 &    , "This is an approximation"
2221    ABI_WARNING(message)
2222  else
2223    write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2224 &   ' Occupation matrix is diagonal : the off-diag part ',maxoffdiag,' is lower than',tol
2225    ABI_COMMENT(message)
2226  endif
2227 
2228  end subroutine checkreal_matlu

m_matlu/chg_repr_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 chg_repr_matlu

FUNCTION

 Change representation of density matrix (useful for nspinor=2)

COPYRIGHT

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

INPUTS

  glocspsp(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: density matrix in the spin spin representation
  glocnm(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: density matrix in the magnetization representation
  natom=number of atoms in cell.
  option= 1 glocspsp is input, glocnm is computed
  option= -1 glocspsp is computed, glocnm is input
  prtopt= option to define level of printing

OUTPUT

  glocspsp(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: density matrix in the spin spin representation
  glocnm(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: density matrix in the magnetization representation

SIDE EFFECTS

NOTES

SOURCE

1036  subroutine chg_repr_matlu(glocspsp,glocnm,natom,option,prtopt)
1037  use defs_basis
1038  use defs_wvltypes
1039  use m_crystal, only : crystal_t
1040  implicit none
1041 
1042 !Arguments ------------------------------------
1043 !scalars
1044  integer,intent(in) :: natom,option,prtopt
1045 !arrays
1046  type(matlu_type),intent(inout) :: glocspsp(natom)
1047  type(matlu_type),intent(inout) :: glocnm(natom)
1048 !Local variables-------------------------------
1049 !scalars
1050  integer :: iatom,isppol,lpawu,m1,m2,ndim,nsppol,mu
1051  complex(dpc) :: ci
1052  character(len=500) :: message
1053 
1054 ! DBG_ENTER("COLL")
1055 
1056  ci=j_dpc
1057 
1058 !==  Compute density matrix in density magnetization representation
1059  if (option==1) then
1060   nsppol=glocspsp(1)%nsppol
1061   do isppol=1,nsppol
1062    do iatom=1,natom
1063     if(glocspsp(iatom)%lpawu/=-1) then
1064      ndim=2*glocspsp(iatom)%lpawu+1
1065      do m1=1,ndim
1066       do m2=1,ndim
1067        glocnm(iatom)%mat(m1,m2,isppol,1,1)=glocspsp(iatom)%mat(m1,m2,isppol,1,1) &
1068 &                                         +glocspsp(iatom)%mat(m1,m2,isppol,2,2)
1069        glocnm(iatom)%mat(m1,m2,isppol,4,1)=glocspsp(iatom)%mat(m1,m2,isppol,1,1) &
1070 &                                         -glocspsp(iatom)%mat(m1,m2,isppol,2,2)
1071        glocnm(iatom)%mat(m1,m2,isppol,2,1)=glocspsp(iatom)%mat(m1,m2,isppol,1,2) &
1072 &                                         +glocspsp(iatom)%mat(m1,m2,isppol,2,1)
1073        glocnm(iatom)%mat(m1,m2,isppol,3,1)= &
1074 &                               cmplx((aimag(glocspsp(iatom)%mat(m1,m2,isppol,2,1))   &
1075 &                                     -aimag(glocspsp(iatom)%mat(m1,m2,isppol,1,2))), &
1076 &                                    (-real(glocspsp(iatom)%mat(m1,m2,isppol,2,1))+  &
1077 &                                      real(glocspsp(iatom)%mat(m1,m2,isppol,1,2))),kind=dp)
1078       enddo  ! m2
1079      enddo ! m1
1080      if(abs(prtopt)>=3) then
1081       write(message,'(a)') "        -- in n, m repr "
1082       call wrtout(std_out,  message,'COLL')
1083       do mu=1,4
1084        do m1=1,ndim
1085         write(message,'(8x,(14(2f9.5,2x)))')(glocnm(iatom)%mat(m1,m2,isppol,mu,1),m2=1,ndim)
1086         call wrtout(std_out,  message,'COLL')
1087        enddo ! m1
1088        write(message,'(a)') ch10
1089        call wrtout(std_out,  message,'COLL')
1090       enddo ! mu
1091      endif ! prtopt >3
1092     endif ! lpawu/=-1
1093    enddo
1094   enddo
1095 
1096 !==  Compute back density matrix in upup dndn updn dnup representation
1097  else  if (option==-1) then
1098   isppol=1
1099   do iatom=1,natom
1100    if(glocnm(iatom)%lpawu/=-1) then
1101     lpawu=glocnm(iatom)%lpawu
1102     ndim=2*glocnm(iatom)%lpawu+1
1103     do m1=1, 2*lpawu+1
1104      do m2=1, 2*lpawu+1
1105       glocspsp(iatom)%mat(m1,m2,isppol,1,1)=half*(glocnm(iatom)%mat(m1,m2,isppol,1,1)+glocnm(iatom)%mat(m1,m2,isppol,4,1))
1106       glocspsp(iatom)%mat(m1,m2,isppol,2,2)=half*(glocnm(iatom)%mat(m1,m2,isppol,1,1)-glocnm(iatom)%mat(m1,m2,isppol,4,1))
1107       glocspsp(iatom)%mat(m1,m2,isppol,1,2)=half*(glocnm(iatom)%mat(m1,m2,isppol,2,1)-ci*glocnm(iatom)%mat(m1,m2,isppol,3,1))
1108       glocspsp(iatom)%mat(m1,m2,isppol,2,1)=half*(glocnm(iatom)%mat(m1,m2,isppol,2,1)+ci*glocnm(iatom)%mat(m1,m2,isppol,3,1))
1109      end do  ! m2
1110     end do ! m1
1111     if(abs(prtopt)>6) then
1112      write(message,'(a)') "        -- in spin spin repr "
1113      call wrtout(std_out,  message,'COLL')
1114      do mu=1,4
1115       do m1=1,ndim
1116        write(message,'(8x,14(2f9.5,2x))')(glocspsp(iatom)%mat(m1,m2,isppol,mu,1),m2=1,ndim)
1117        call wrtout(std_out,  message,'COLL')
1118       enddo
1119       write(message,'(a)') ch10
1120       call wrtout(std_out,  message,'COLL')
1121      enddo
1122     endif !prtopt>3
1123    endif ! lpawu/=-1
1124   end do ! iatom
1125  else
1126   message = "stop in chg_repr_matlu"
1127   ABI_ERROR(message)
1128  endif
1129 
1130 
1131 ! DBG_EXIT("COLL")
1132 
1133  end subroutine chg_repr_matlu

m_matlu/chi_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 chi_matlu

FUNCTION

 return the matrix of dimension [(2*ll+1)]**4 in the Ylm basis
 with the matrix elements of the orbital (option=1), spin (option=2) and
 total (option=3) angular momentum for z direction. It is used by the 
 QMC for the correlation function of the magnetic moment

COPYRIGHT

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

INPUTS

 matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: 
 natom :: number of atoms
 option = 1 :: Orbital angular momentum along z axis
        = 2 :: 2*Spin angluar momentum alonf z axis
        = 3 :: total angular momentum along z axis
 optptr > 2 :: print angular matrix elements 

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: quantity in Ylm basis

SIDE EFFECTS

NOTES

SOURCE

3677  subroutine chi_matlu(matlu,natom,option,optprt)
3678  use defs_basis
3679  use defs_wvltypes
3680  implicit none
3681 
3682 !Arguments ------------------------------------
3683 !scalars
3684  integer, intent(in) :: natom,option,optprt
3685 !arrays
3686  type(matlu_type), intent(inout) :: matlu(natom)
3687 !Local variables-------------------------------
3688 !scalars
3689  integer :: iatom,im
3690  integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,tndim
3691  character(len=500) :: message
3692  real(dp) :: xj
3693 !arrays
3694  integer, allocatable :: ind_msml(:,:)
3695 ! type(coeff2c_type), allocatable :: gathermatlu(:)
3696  type(coeff2c_type), allocatable :: muchi(:)
3697 !************************************************************************
3698 
3699  !=====================================
3700  ! Allocate matrices 
3701  !=====================================
3702     
3703  ABI_MALLOC(muchi,(natom))
3704  do iatom=1,natom
3705    if(matlu(iatom)%lpawu.ne.-1) then
3706      tndim=2*(2*matlu(iatom)%lpawu+1)
3707      ABI_MALLOC(muchi(iatom)%value,(tndim,tndim))
3708      muchi(iatom)%value=czero
3709    end if
3710  end do 
3711 
3712  do iatom=1,natom
3713    lpawu=matlu(iatom)%lpawu
3714    if(lpawu.ne.-1) then
3715      ll=lpawu
3716      lcor=lpawu
3717      ABI_MALLOC(ind_msml,(2,-ll:ll))
3718      ind_msml=czero
3719      !=====================================
3720      !Orbital angular momentum matrix along z axis
3721      !=====================================
3722      if(option==1) then
3723        jc1=0
3724        do ms1=1,2
3725          do ml1=-ll,ll
3726           jc1=jc1+1
3727           ind_msml(ms1,ml1)=jc1
3728          end do
3729        end do
3730 
3731        jc1=0
3732        do ms1=-1,1
3733          do ml1=-ll,ll
3734             jc1=jc1+1
3735             if(jc1 <tndim+1) then
3736                 muchi(iatom)%value(jc1,jc1) = ml1
3737             endif
3738           end do
3739         end do
3740 
3741      !=====================================
3742      !Spin angular momentum matrix along z axis
3743      !up spin is first
3744      !=====================================
3745      else if(option==2) then
3746         jc1=0
3747        do ms1=1,2
3748          do ml1=-ll,ll
3749           jc1=jc1+1
3750           ind_msml(ms1,ml1)=jc1
3751          end do
3752        end do
3753 
3754        jc1=0
3755        do ms1=-1,1
3756          xj=float(ms1)+half 
3757          do ml1=-ll,ll
3758             jc1=jc1+1
3759             if(jc1<tndim+1) then
3760                 if(xj < 0.0 ) then
3761                         muchi(iatom)%value(jc1,jc1) = -2*xj
3762                 else if(xj > 0.0) then
3763                         muchi(iatom)%value(jc1,jc1) = -2*xj
3764 
3765                 end if
3766              endif
3767          end do
3768        end do
3769 
3770      !=====================================
3771      !Total angular momentum matrix along z axis
3772      !up spin is first
3773      !=====================================
3774      else if(option==3) then
3775        jc1=0
3776        do ms1=-1,1
3777          xj=float(ms1)+half 
3778          do ml1=-ll,ll
3779             jc1=jc1+1
3780             if(jc1<tndim+1) then
3781                 if(xj < 0.0 ) then
3782                         muchi(iatom)%value(jc1,jc1) = ml1-2*xj
3783                 else if(xj > 0.0) then
3784                         muchi(iatom)%value(jc1,jc1) = ml1-2*xj
3785                 end if
3786             endif
3787          end do
3788        end do
3789      end if
3790 
3791      if(optprt>2) then
3792         if(option==1) then
3793           write(message,'(a)') "Orbital angular momentum matrix elements in |m_l,m_s> basis"
3794         else if(option==2) then
3795           write(message,'(a)') "Spin angular momentum matrix elements in |m_l,m_s> basis"
3796         else if(option==3) then
3797           write(message,'(a)') "Zeeman angular momentum matrix elements in |m_l,m_s> basis"
3798         end if
3799         call wrtout(std_out,message,"COLL")
3800         do im=1,2*(ll*2+1)
3801           write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muchi(iatom)%value(im,jm),jm=1,2*(ll*2+1))
3802           call wrtout(std_out,message,"COLL")
3803         end do
3804      end if
3805 
3806    ABI_FREE(ind_msml) 
3807 
3808      !=====================================
3809      ! Reshape matrix into matlu format 
3810      !=====================================
3811 
3812     call gather_matlu(matlu,muchi(iatom),natom=1,option=-1,prtopt=1)
3813 
3814 
3815    end if !lpawu
3816  end do !atom
3817 
3818      !=====================================
3819      ! Deallocate gathermatlu 
3820      !=====================================
3821 
3822  do iatom=1,natom
3823    if(matlu(iatom)%lpawu.ne.-1) then
3824      ABI_FREE(muchi(iatom)%value)
3825    end if
3826  end do
3827  ABI_FREE(muchi)
3828  
3829  end subroutine chi_matlu

m_matlu/conjg_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 conjg_matlu

FUNCTION

 conjugate of input matlu

COPYRIGHT

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

INPUTS

  matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity

SIDE EFFECTS

NOTES

SOURCE

2409  subroutine conjg_matlu(matlu1,natom)
2410  use defs_basis
2411  use defs_wvltypes
2412  implicit none
2413 
2414 !Arguments ------------------------------------
2415 !scalars
2416  integer, intent(in) :: natom
2417 !arrays
2418  type(matlu_type), intent(inout) :: matlu1(natom)
2419 !Local variables-------------------------------
2420 !scalars
2421  integer :: iatom,im1,im2,ispinor2,ispinor1,isppol
2422  integer :: lpawu
2423 !arrays
2424 !************************************************************************
2425  do iatom=1,natom
2426    lpawu=matlu1(iatom)%lpawu
2427    if(lpawu.ne.-1) then
2428      do isppol=1,matlu1(1)%nsppol
2429        do ispinor1=1,matlu1(1)%nspinor
2430          do ispinor2=1,matlu1(1)%nspinor
2431            do im1=1,2*lpawu+1
2432              do im2=1,2*lpawu+1
2433                matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)= &
2434 &               conjg(matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2))
2435              enddo ! im2
2436            enddo ! im1
2437          enddo ! ispinor2
2438        enddo ! ispinor1
2439      enddo ! isppol
2440    endif ! lpawu
2441  enddo ! iatom
2442 
2443  end subroutine conjg_matlu

m_matlu/copy_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 copy_matlu

FUNCTION

  Copy matlu1 into matlu2

INPUTS

  maltu1 <type(matlu_type)>= density matrix matlu1 in the local orbital basis and related variables
  natom = number of atoms

OUTPUT

  maltu2 <type(matlu_type)>= density matrix matlu2 in the local orbital basis and related variables

SOURCE

294 subroutine copy_matlu(nmat1,nmat2,natom,opt_diag,opt_non_diag,opt_re)
295 
296  use defs_basis
297  implicit none
298 
299 !Arguments ------------------------------------
300 !type
301  integer, intent(in) :: natom
302  type(matlu_type),intent(in) :: nmat1(natom)
303  type(matlu_type),intent(inout) :: nmat2(natom) !vz_i
304  integer, optional, intent(in) :: opt_diag,opt_non_diag,opt_re
305 
306 !Local variables-------------------------------
307  integer :: iatom,isppol,im1,im2,ispinor,ispinor1,tndim
308 ! *********************************************************************
309 
310  do isppol=1,nmat1(1)%nsppol
311    do iatom=1,natom
312      if(nmat1(iatom)%lpawu.ne.-1) then
313        tndim=(2*nmat1(iatom)%lpawu+1)
314        do im1=1,tndim
315          do im2=1,tndim
316            do ispinor=1,nmat1(1)%nspinor
317              do ispinor1=1,nmat1(1)%nspinor
318                if(present(opt_diag)) then
319                  nmat2(iatom)%mat(im1,im1,isppol,ispinor,ispinor)=nmat1(iatom)%mat(im1,im1,isppol,ispinor,ispinor)
320                else if(present(opt_non_diag)) then
321                  if(im1/=im2.or.ispinor/=ispinor1) then
322                    nmat2(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)=nmat1(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)
323                  endif
324                else
325                  nmat2(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)=nmat1(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)
326                  if(present(opt_re)) nmat2(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)=&
327 &                         cmplx(real(nmat1(iatom)%mat(im1,im2,isppol,ispinor,ispinor1)),0.d0,kind=dp)
328                endif
329              enddo
330            enddo
331          enddo
332        enddo
333      endif ! lpawu=-1
334    enddo ! iatom
335  enddo ! isppol
336 !   do iatom=1,natom
337 !    lpawu=nmat1(iatom)%lpawu
338 !    if(lpawu.ne.-1) then
339 !     nmat2(iatom)%mat=nmat1(iatom)%mat
340 !    endif
341 !   enddo
342 
343 
344 end subroutine copy_matlu

m_matlu/destroy_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 destroy_matlu

FUNCTION

  deallocate matlu

INPUTS

  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables
  natom = number of atoms

OUTPUT

SOURCE

253 subroutine destroy_matlu(matlu,natom)
254 
255  use defs_basis
256  use m_crystal, only : crystal_t
257  implicit none
258 
259 !Arguments ------------------------------------
260 !scalars
261  integer, intent(in) :: natom
262  type(matlu_type),intent(inout) :: matlu(natom)
263 
264 !Local variables-------------------------------
265  integer :: iatom
266 
267 ! *********************************************************************
268 
269  do iatom=1,natom
270   if ( allocated(matlu(iatom)%mat) )   then
271     ABI_FREE(matlu(iatom)%mat)
272   end if
273  enddo
274 
275 end subroutine destroy_matlu

m_matlu/diag_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 diag_matlu

FUNCTION

 Diagonalize matlu matrix

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to diagonalize
  natom=number of atoms
  prtopt= option to define level of printing
  nsppol_imp= if 1, one can diagonalize with the same matrix the Up
   and Dn matlu matrix. It is convenient because one can thus have the
 same interaction matrix for up and dn spins.

OUTPUT

  matlu_diag(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: diagonalized density matrix
  eigvectmatlu(natom) <type(coeff2c_type)> = Eigenvectors corresponding to the diagonalization

SIDE EFFECTS

NOTES

SOURCE

1390  subroutine diag_matlu(matlu,matlu_diag,natom,prtopt,eigvectmatlu,nsppol_imp,checkstop,optreal,test)
1391  use defs_basis
1392  use defs_wvltypes
1393  use m_crystal, only : crystal_t
1394  use m_matrix,         only : blockdiago_fordsyev,blockdiago_forzheev
1395  implicit none
1396 
1397 !Arguments ------------------------------------
1398 !scalars
1399  integer, intent(in) :: natom
1400  integer, intent(in) :: prtopt
1401 !arrays
1402  type(matlu_type),intent(inout) :: matlu(natom)
1403  type(matlu_type),intent(inout) :: matlu_diag(natom) !vz_i
1404  type(coeff2c_type),optional,intent(inout) :: eigvectmatlu(natom,matlu(1)%nsppol) !vz_i
1405  integer,optional,intent(in) :: nsppol_imp
1406  logical,optional,intent(in) :: checkstop
1407  integer,optional,intent(in) :: optreal
1408  integer,optional,intent(in) :: test
1409 !Local variables-------------------------------
1410 !scalars
1411  integer :: iatom,im1,im2,im3,imc,imc1,info,ispinor,ispinor1,isppol,lwork,tndim,lworkr
1412  integer :: nsppol,nsppolimp,nspinor
1413  logical :: checkstop_in,blockdiag
1414  character(len=500) :: message
1415 !arrays
1416  type(coeff2c_type),allocatable :: gathermatlu(:)
1417  real(dp),allocatable :: eig(:),rwork(:),work(:),valuer(:,:)!,valuer2(:,:)
1418  !real(dp),allocatable :: valuer3(:,:),valuer4(:,:)
1419 ! real(dp),allocatable :: eigvec(:,:)
1420  complex(dpc),allocatable :: zwork(:)
1421  logical :: donotdiag,print_temp_mat2
1422  complex(dpc),allocatable :: temp_mat(:,:)
1423  complex(dpc),allocatable :: temp_mat2(:,:)
1424 !debug complex(dpc),allocatable :: temp_mat3(:,:)
1425 !************************************************************************
1426 
1427  call zero_matlu(matlu_diag,natom)
1428  nsppol=matlu(1)%nsppol
1429  if(present(nsppol_imp)) then
1430    nsppolimp=nsppol_imp
1431  else
1432    nsppolimp=nsppol
1433  endif
1434  if(present(checkstop)) then
1435   checkstop_in=checkstop
1436  else
1437   checkstop_in=.true.
1438  endif
1439  if(present(test)) then
1440   blockdiag=(test==8)
1441  else
1442   blockdiag=.false.
1443  endif
1444  nspinor=matlu(1)%nspinor
1445 
1446  donotdiag=.true.
1447  donotdiag=.false.
1448 ! ===========================
1449 ! Check is diagonalization is necessary and how
1450 ! ===========================
1451  do isppol=1,matlu(1)%nsppol
1452    do iatom=1,natom
1453       if(matlu(iatom)%lpawu.ne.-1) then
1454        tndim=(2*matlu(iatom)%lpawu+1)
1455         do im1=1,tndim
1456           do im2=1,tndim
1457             do ispinor=1,nspinor
1458               do ispinor1=1,nspinor
1459                 if(abs(matlu(iatom)%mat(im1,im2,isppol,ispinor,ispinor1))>tol8.and.&
1460 &                 (im1/=im2.or.ispinor/=ispinor1)) then
1461 !                 if matrix is diagonal: do not diagonalize
1462                   donotdiag=.false.
1463                   exit
1464                 endif
1465               enddo
1466             enddo
1467           enddo
1468         enddo
1469       endif
1470    enddo
1471  enddo
1472 
1473  if(donotdiag) then
1474    do isppol=1,matlu(1)%nsppol
1475      do iatom=1,natom
1476         if(matlu(iatom)%lpawu.ne.-1) then
1477           tndim=(2*matlu(iatom)%lpawu+1)
1478           eigvectmatlu(iatom,isppol)%value(:,:)=czero
1479           do im1=1,tndim
1480             eigvectmatlu(iatom,isppol)%value(im1,im1)=cone
1481           enddo
1482         endif
1483      enddo
1484    enddo
1485    call copy_matlu(matlu,matlu_diag,natom)
1486    write(message,'(a)')  "   Diagonalisation of matlu will not be performed"
1487    call wrtout(std_out,message,'COLL')
1488    return
1489  endif
1490 
1491 
1492 ! Below, gathermatlu is used to store eigenvectors (after diagonalization)
1493 ! For nsppol=2, and if nsppolimp=1, these eigenvectors computed for isppol=1, and applied to
1494 ! rotate matrix for isppol=2. It is the reason why the sum below is only
1495 ! from 1 to nsppolimp !
1496 
1497 
1498 
1499  do isppol=1,nsppolimp !
1500 ! ===========================
1501 ! Define gathermatlu
1502 ! ===========================
1503    ABI_MALLOC(gathermatlu,(natom))
1504    do iatom=1,natom
1505      if(matlu(iatom)%lpawu.ne.-1) then
1506        tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1507        ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
1508        gathermatlu(iatom)%value=czero
1509      endif
1510    enddo
1511    if(nsppol==1.and.nspinor==2) then
1512      call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
1513    else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
1514      do iatom=1,natom
1515        if(matlu(iatom)%lpawu.ne.-1) then
1516          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1517          do im1=1,tndim
1518            do im2=1,tndim
1519              gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,isppol,1,1)
1520            enddo
1521          enddo
1522        endif
1523      enddo
1524    endif
1525 ! ===========================
1526 ! Diagonalize
1527 ! ===========================
1528    do iatom=1,natom
1529      if(matlu(iatom)%lpawu.ne.-1) then
1530        if(isppol<=nsppolimp) then
1531          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1532 !debug       allocate(temp_mat2(tndim,tndim))
1533 !debug       temp_mat2=zero
1534          lwork=2*tndim-1
1535          ABI_MALLOC(rwork,(3*tndim-2))
1536          rwork = zero
1537 
1538          lworkr=tndim*(tndim+2)*2
1539          ABI_MALLOC(work,(lworkr))
1540          work = zero
1541          ABI_MALLOC(valuer,(tndim,tndim))
1542 !         ABI_MALLOC(valuer2,(tndim,tndim))
1543 !         ABI_MALLOC(valuer3,(tndim,tndim))
1544 !         ABI_MALLOC(valuer4,(tndim,tndim))
1545          ABI_MALLOC(zwork,(lwork))
1546 !         valuer2=zero
1547 !         valuer3=zero
1548 !         valuer4=zero
1549          zwork = czero
1550          ABI_MALLOC(eig,(tndim))
1551          eig = zero
1552          info = 0
1553          if(prtopt>=4) then
1554            write(message,'(a,i4,a,i4)')  "       BEFORE DIAGONALIZATION for atom",iatom,"  and isppol",isppol
1555            call wrtout(std_out,message,'COLL')
1556            do im1=1,tndim
1557              write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1558 &             (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
1559              call wrtout(std_out,message,'COLL')
1560            end do
1561          endif
1562 !debug       temp_mat2(:,:)=gathermatlu(iatom)%value(:,:)
1563 !           write(std_out,*)"diag"
1564          if(present(optreal).and.maxval(abs(aimag(gathermatlu(iatom)%value(:,:))))<tol6) then
1565            write(message,'(a,2x,a,e9.3,a)') ch10,"Imaginary part of Local Hamiltonian is lower than ",&
1566 &                   tol8, ": the real matrix is used"
1567            call wrtout(std_out,message,'COLL')
1568            valuer=real(gathermatlu(iatom)%value,kind=dp)
1569 !           write(message,'(a)') ch10
1570 !           call wrtout(std_out,message,'COLL')
1571 !           write(message,'(a,i4,a,i4)')  "BEFORE valuer for atom",iatom,"  and isppol",isppol
1572 !           call wrtout(std_out,message,'COLL')
1573 !           do im1=1,tndim
1574 !             write(message,'(2(1x,18(1x,"(",f20.15,",",f20.15,")")))')&
1575 !&             (valuer(im1,im2),im2=1,tndim)
1576 !             call wrtout(std_out,message,'COLL')
1577 !           end do
1578 !           do im1=1,tndim
1579 !             valuer(im1,im1)=real(im1,kind=dp)*0.00000000001_dp+valuer(im1,im1)
1580 !           enddo
1581 !           write(message,'(a)') ch10
1582 !           call wrtout(std_out,message,'COLL')
1583 !           write(message,'(a,i4,a,i4)')  "BEFORE valuer for atom",iatom,"  and isppol",isppol
1584 !           call wrtout(std_out,message,'COLL')
1585 !           do im1=1,tndim
1586 !             write(message,'(2(1x,18(1x,f20.15,f20.15)))')&
1587 !&             (valuer(im1,im2),im2=1,tndim)
1588 !             call wrtout(std_out,message,'COLL')
1589 !           end do
1590            !call dsyev('v','u',tndim,valuer,tndim,eig,work,lworkr,info)
1591            if (blockdiag) then
1592              call blockdiago_fordsyev(valuer,tndim,eig)
1593            else
1594              call dsyev('v','u',tndim,valuer,tndim,eig,work,lworkr,info)
1595            endif
1596 !!       For reproductibility
1597 !           ! valuer2: eigenvector for the perturb matrix
1598 !           valuer2=real(gathermatlu(iatom)%value,kind=dp)
1599 !           do im1=1,tndim
1600 !             valuer2(im1,im1)=float(im1)*0.00000000001+valuer2(im1,im1)
1601 !           enddo
1602 !           call dsyev('v','u',tndim,valuer2,tndim,eig,work,lworkr,info)
1603 !           write(message,'(a)') ch10
1604 !           call wrtout(std_out,message,'COLL')
1605 !           write(message,'(a,i4,a,i4)')  "       valuer2 for atom",iatom,"  and isppol",isppol
1606 !           call wrtout(std_out,message,'COLL')
1607 !           do im1=1,tndim
1608 !             write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1609 !&             (valuer2(im1,im2),im2=1,tndim)
1610 !             call wrtout(std_out,message,'COLL')
1611 !           end do
1612 !           call dgemm('n','n',tndim,tndim,tndim,cone,valuer,tndim,&
1613 !&            valuer2,tndim,czero,valuer3                ,tndim)
1614 !           call dgemm('c','n',tndim,tndim,tndim,cone,valuer2,tndim,&
1615 !&            valuer3                   ,tndim,czero,valuer4,tndim)
1616 !           ! valuer4: compute unpert matrix in the basis of the
1617 !           ! perturb basis
1618 !           write(message,'(a)') ch10
1619 !           call wrtout(std_out,message,'COLL')
1620 !           write(message,'(a,i4,a,i4)')  "BEFORE valuer4 for atom",iatom,"  and isppol",isppol
1621 !           call wrtout(std_out,message,'COLL')
1622 !           do im1=1,tndim
1623 !             write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1624 !&             (valuer4(im1,im2),im2=1,tndim)
1625 !             call wrtout(std_out,message,'COLL')
1626 !           end do
1627 !           call dsyev('v','u',tndim,valuer4,tndim,eig,work,lworkr,info)
1628 !           ! valuer4: Diago valuer4 (nearly diag)
1629 !           write(message,'(a)') ch10
1630 !           call wrtout(std_out,message,'COLL')
1631 !           write(message,'(a,i4,a,i4)')  "AFTER  valuer4 for atom",iatom,"  and isppol",isppol
1632 !           call wrtout(std_out,message,'COLL')
1633 !           do im1=1,tndim
1634 !             write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1635 !&             (valuer4(im1,im2),im2=1,tndim)
1636 !             call wrtout(std_out,message,'COLL')
1637 !           end do
1638 !           call dgemm('n','n',tndim,tndim,tndim,cone,valuer2,tndim,&
1639 !&            valuer4,tndim,czero,valuer                ,tndim)
1640            !write(6,*) "INFO",info
1641            gathermatlu(iatom)%value=cmplx(valuer,0.d0,kind=dp)
1642 !           write(message,'(a,i4,a,i4)')  "AFTER valuer for atom",iatom,"  and isppol",isppol
1643 !           call wrtout(std_out,message,'COLL')
1644 !           do im1=1,tndim
1645 !             write(message,'(2(1x,18(1x,"(",f20.15,",",f20.15,")")))')&
1646 !&             (valuer(im1,im2),im2=1,tndim)
1647 !             call wrtout(std_out,message,'COLL')
1648 !           end do
1649          else
1650            if(present(optreal).and.maxval(abs(aimag(gathermatlu(iatom)%value(:,:))))>tol8) then
1651              write(message,'(a)') " Local hamiltonian in correlated basis is complex"
1652              ABI_COMMENT(message)
1653            endif
1654            call zheev('v','u',tndim,gathermatlu(iatom)%value,tndim,eig,zwork,lwork,rwork,info)
1655            !call blockdiago_forzheev(gathermatlu(iatom)%value,tndim,eig)
1656          endif
1657          if(prtopt>=3) then
1658            write(message,'(a)') ch10
1659            call wrtout(std_out,message,'COLL')
1660            write(message,'(a,i4,a,i4)')  "       EIGENVECTORS for atom",iatom,"  and isppol",isppol
1661            call wrtout(std_out,message,'COLL')
1662            do im1=1,tndim
1663              write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1664 &             (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
1665              call wrtout(std_out,message,'COLL')
1666            end do
1667           ! do im1=1,tndim
1668           !   xcheck=czero
1669           !   do im3=1,tndim
1670           !     do im2=1,tndim
1671           !       xcheck=xcheck+gathermatlu(iatom)%value(im1,im2)*conjg(gathermatlu(iatom)%value(im2,im3))
1672           !     end do
1673           !   end do
1674           !   write(6,*) "check",im3,im1,xcheck
1675           ! end do
1676          endif
1677 !       write(std_out,*) "eig",eig
1678 ! ===========================
1679 ! Put eigenvalue in matlu_diag
1680 ! ===========================
1681          imc=0
1682          do ispinor=1,nspinor
1683            do im1=1,2*matlu(iatom)%lpawu+1
1684              imc=imc+1
1685              matlu_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)=eig(imc)
1686 !debug             temp_mat2(imc,imc)=eig(imc)
1687            enddo
1688          enddo
1689          if(prtopt>=2) then
1690 !            write(message,'(a,12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1691 !&            ch10,(eig(im1),im1=1,tndim)
1692 !             call wrtout(std_out,message,'COLL')
1693            !call wrtout(std_out,message,'COLL')
1694             !write(std_out,*) "EIG", eig
1695          endif
1696          ABI_FREE(zwork)
1697          ABI_FREE(rwork)
1698          ABI_FREE(work)
1699          ABI_FREE(valuer)
1700 !         ABI_FREE(valuer2)
1701 !         ABI_FREE(valuer3)
1702 !         ABI_FREE(valuer4)
1703          ABI_FREE(eig)
1704 !     endif
1705 !   enddo
1706 ! ===========================
1707 ! Keep eigenvectors gathermatlu
1708 ! ===========================
1709          if (present(eigvectmatlu)) then
1710            tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1711            eigvectmatlu(iatom,isppol)%value(:,:)=gathermatlu(iatom)%value(:,:)
1712 !           write(std_out,*) "eigvect in diag_matlu"
1713 !           do im1=1,tndim
1714 !             write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1715 !&             (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
1716 !             call wrtout(std_out,message,'COLL')
1717 !           end do
1718          endif
1719 
1720 
1721 ! ==================================================================
1722        endif
1723        if(nsppolimp==1.and.matlu(1)%nsppol==2) then
1724 ! If necessary rotate levels for this other spin, assuming the same
1725 ! rotation matrix: it have to be checked afterwards that the matrix is
1726 ! diagonal
1727 ! ===================================================================
1728          ABI_MALLOC(temp_mat,(tndim,tndim))
1729          ABI_MALLOC(temp_mat2,(tndim,tndim))
1730          temp_mat(:,:)=czero
1731 !        input matrix: gathermatlu
1732 !        rotation matrix: eigvectmatlu
1733 !        intermediate matrix: temp_mat
1734 !        result matrix: temp_mat2
1735          do im1=1,tndim
1736            do im2=1,tndim
1737              gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,2,1,1)
1738            enddo
1739          enddo
1740          if(prtopt>=3) then
1741            write(message,'(a,i4,a,i4)')  "       GATHERMATLU for atom",iatom," inside if nsppolimp==1"
1742            call wrtout(std_out,message,'COLL')
1743            do im1=1,tndim
1744              write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1745 &             (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
1746              call wrtout(std_out,message,'COLL')
1747            end do
1748          endif
1749          call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,&
1750 &          eigvectmatlu(iatom,1)%value,tndim,czero,temp_mat                ,tndim)
1751          call zgemm('c','n',tndim,tndim,tndim,cone,eigvectmatlu(iatom,1)%value,tndim,&
1752 &          temp_mat                   ,tndim,czero,temp_mat2,tndim)
1753          eigvectmatlu(iatom,2)%value=eigvectmatlu(iatom,1)%value
1754          imc=0
1755          print_temp_mat2=.false.
1756          do ispinor=1,nspinor
1757            do im1=1,2*matlu(iatom)%lpawu+1
1758              imc=imc+1
1759              imc1=0
1760              do ispinor1=1,nspinor
1761                do im2=1,2*matlu(iatom)%lpawu+1
1762                  imc1=imc1+1
1763                  matlu_diag(iatom)%mat(im1,im2,2,ispinor,ispinor1)=temp_mat2(imc,imc1)
1764                  if (imc/=imc1.and.(abs(temp_mat2(imc,imc1))>tol5)) then
1765                    write(message,'(3a,i4,2f16.4)') ch10,'diag_matlu= Matrix for spin number 2 obtained with', &
1766 &                   ' eigenvectors from diagonalization for spin nb 1 is non diagonal for atom:',iatom,&
1767 &                    abs(temp_mat2(imc,imc1)),tol5
1768                    call wrtout(std_out,message,'COLL')
1769                    if(.not.checkstop_in.and.(abs(temp_mat2(imc,imc1))>0.1)) print_temp_mat2=.true.
1770                    if(checkstop_in) print_temp_mat2=.true.
1771                  endif
1772 !                 write(std_out,*) temp_mat2(imc,imc1)
1773 !debug             temp_mat2(imc,imc)=eig(imc)
1774                enddo
1775              enddo
1776            enddo
1777          enddo
1778          if(print_temp_mat2.and.prtopt>=3) then
1779            write(message,'(a)')  "       temp_mat2"
1780            call wrtout(std_out,message,'COLL')
1781            do im1=1,tndim
1782              write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1783 &             (temp_mat2(im1,im3),im3=1,tndim)
1784              call wrtout(std_out,message,'COLL')
1785            end do
1786            if(iatom==2) then
1787              ABI_ERROR("iatom==2")
1788            end if
1789          endif
1790          ABI_FREE(temp_mat)
1791          ABI_FREE(temp_mat2)
1792        endif
1793 
1794 
1795      endif ! end lpawu/=-1
1796 
1797 !!  for check only
1798 !debug     if(matlu(iatom)%lpawu.ne.-1) then
1799 !debug       allocate(temp_mat(tndim,tndim))
1800 !debug       allocate(temp_mat3(tndim,tndim))
1801 !debug           do im1=1,tndim
1802 !debug             do im2=1,tndim
1803 !debug!               rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom,isppol)%value(im1,im2)
1804 !debug               temp_mat3(im1,im2)=conjg(gathermatlu(iatom)%value(im2,im1))
1805 !debug             enddo
1806 !debug           enddo
1807 !debug       temp_mat(:,:)=czero
1808 !debug!      input matrix: temp_mat2
1809 !debug!      rotation matrix: gathermatlu
1810 !debug!      intermediate matrix: temp_mat
1811 !debug!      result matrix: temp_mat2
1812 !debug       call zgemm('n','c',tndim,tndim,tndim,cone,temp_mat2   ,tndim,&
1813 !debug&        temp_mat3,tndim,czero,temp_mat                ,tndim)
1814 !debug       call zgemm('n','n',tndim,tndim,tndim,cone,temp_mat3,tndim,&
1815 !debug&        temp_mat                   ,tndim,czero,temp_mat2,tndim)
1816 !debug!       call zgemm('n','c',tndim,tndim,tndim,cone,temp_mat2   ,tndim,&
1817 !debug!&        gathermatlu(iatom)%value,tndim,czero,temp_mat                ,tndim)
1818 !debug!       call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,&
1819 !debug!&        temp_mat                   ,tndim,czero,temp_mat2,tndim)
1820 !debug         write(std_out,*) "result"
1821 !debug         do im1=1,tndim
1822 !debug           write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1823 !debug&           (temp_mat2(im1,im2),im2=1,tndim)
1824 !debug           call wrtout(std_out,message,'COLL')
1825 !debug         end do
1826 !debug       deallocate(temp_mat)
1827 !debug       deallocate(temp_mat3)
1828 !debug     endif ! lpawu
1829 
1830    enddo  ! iatom
1831 ! End loop over atoms
1832 ! ===========================
1833    do iatom=1,natom
1834      if(matlu(iatom)%lpawu.ne.-1) then
1835 !debug       deallocate(temp_mat2)
1836        ABI_FREE(gathermatlu(iatom)%value)
1837      endif
1838    enddo
1839    ABI_FREE(gathermatlu)
1840  enddo ! isppol
1841 
1842  end subroutine diag_matlu

m_matlu/diff_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 diff_matlu

FUNCTION

INPUTS

  char1 = character describing matlu1
  char2 = character describing matlu2
  matlu1(natom) <type(matlu_type)>= density matrix 1 in the local orbital basis and related variables
  matlu2(natom) <type(matlu_type)>= density matrix 2 in the local orbital basis and related variables
  natom = number of atoms
  option =1      if diff> toldiff , stop
          0      print diff and toldiff
          else   do not test and do not print
  toldiff = maximum value for the difference between matlu1 and matlu2

OUTPUT

SOURCE

874 subroutine diff_matlu(char1,char2,matlu1,matlu2,natom,option,toldiff,ierr,zero_or_one)
875 
876  use defs_basis
877  use m_paw_dmft, only : paw_dmft_type
878  use m_crystal, only : crystal_t
879  use m_io_tools,           only : flush_unit
880  implicit none
881 
882 !Arguments ------------------------------------
883 !type
884  integer,intent(in) :: natom,option
885  type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
886  character(len=*), intent(in) :: char1,char2
887  real(dp),intent(in) :: toldiff
888  integer,intent(out), optional :: ierr
889  integer,intent(in), optional :: zero_or_one
890 
891 !local variables-------------------------------
892  integer :: iatom,idiff,ispinor,ispinor1,isppol,m1,m,lpawu,nspinor,nsppol
893  real(dp) :: matludiff
894  character(len=500) :: message
895 ! *********************************************************************
896  nsppol=matlu1(1)%nsppol
897  nspinor=matlu1(1)%nspinor
898 
899  matludiff=zero
900  idiff=0
901  do iatom = 1 , natom
902   lpawu=matlu1(iatom)%lpawu
903   if(lpawu/=-1) then
904    do isppol = 1 , nsppol
905     do ispinor = 1 , nspinor
906      do ispinor1 = 1, nspinor
907       do m1 = 1 , 2*lpawu+1
908        do m = 1 ,  2*lpawu+1
909         idiff=idiff+1
910         matludiff=matludiff+ &
911 &        sqrt( real(matlu1(iatom)%mat(m1,m,isppol,ispinor,ispinor1)      &
912 &            -      matlu2(iatom)%mat(m1,m,isppol,ispinor,ispinor1))**2  &
913 &            +aimag(matlu1(iatom)%mat(m1,m,isppol,ispinor,ispinor1)      &
914 &            -      matlu2(iatom)%mat(m1,m,isppol,ispinor,ispinor1))**2  )
915 !       write(std_out,*) m,m1,matlu1(iatom)%mat(m1,m,isppol,ispinor,ispinor1),matlu2(iatom)%mat(m1,m,isppol,ispinor,ispinor1),matludiff
916        enddo
917       enddo
918      end do ! ispinor1
919     end do ! ispinor
920    enddo ! isppol
921   endif ! lpawu/=1
922  enddo ! natom
923  if(.not.present(zero_or_one)) matludiff=matludiff/float(idiff)
924 
925  if(option==1.or.option==0) then
926   if( matludiff < toldiff ) then
927    write(message,'(5a,6x,3a,4x,e12.4,a,e12.4)') ch10,&
928 &   '   ** Differences between ',trim(char1),' and ',ch10,trim(char2),' are small enough:',&
929 &   ch10,matludiff,' is lower than',toldiff
930    call wrtout(std_out,message,'COLL')
931    if(present(ierr)) ierr=0
932   else 
933    write(message,'(5a,3x,3a,3x,e12.4,a,e12.4)') ch10,&
934 &   'Differences between ',trim(char1),' and ',ch10,trim(char2),' is too large:',&
935 &   ch10,matludiff,' is larger than',toldiff
936    ABI_WARNING(message)
937 !   write(message,'(8a,4x,e12.4,a,e12.4)') ch10,"  Matrix for ",trim(char1)
938    write(message,'(a,3x,a)') ch10,trim(char1)
939    call wrtout(std_out,message,'COLL')
940    call print_matlu(matlu1,natom,prtopt=1,opt_diag=-1)
941    write(message,'(a,3x,a)') ch10,trim(char2)
942    call wrtout(std_out,message,'COLL')
943    call print_matlu(matlu2,natom,prtopt=1,opt_diag=-1)
944    if (present(zero_or_one).and.(mod(matludiff,1.d0)< toldiff)) then
945      write(message,'(a,3x,a)') ch10," The norm is not identity for this k-point but&
946     & is compatible with a high symmetry point"
947      call wrtout(std_out,message,'COLL')
948    else if(present(zero_or_one)) then
949      write(message,'(a,3x,a)') ch10," The norm is not identity for this k-point but&
950     & might be compatible with a high symmetry point: it should be checked"
951      call wrtout(std_out,message,'COLL')
952    else 
953      if(option==1) then
954        call flush_unit(std_out)
955        ABI_ERROR("option==1, aborting now!")
956      end if
957    end if
958    if(present(ierr)) ierr=-1
959   endif
960  endif
961 
962 end subroutine diff_matlu

m_matlu/fac_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 fac_matlu

FUNCTION

 shift matlu matrix

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.
  shift= shift of the diagonal part.

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: shifted matrix

SIDE EFFECTS

NOTES

SOURCE

2670  subroutine fac_matlu(matlu,natom,fac)
2671  use defs_basis
2672  implicit none
2673 
2674 !Arguments ------------------------------------
2675 !scalars
2676  integer, intent(in) :: natom
2677 !arrays
2678  type(matlu_type),intent(inout) :: matlu(natom)
2679  complex(dpc),intent(in)      :: fac
2680 !Local variables-------------------------------
2681 !scalars
2682  integer :: iatom,im,im1,ispinor,ispinor1,isppol
2683  integer :: lpawu
2684 ! character(len=500) :: message
2685 !arrays
2686 !************************************************************************
2687  do iatom=1,natom
2688    lpawu=matlu(iatom)%lpawu
2689    if(lpawu.ne.-1) then
2690      do im=1,2*lpawu+1
2691        do im1=1,2*lpawu+1
2692          do isppol=1,matlu(1)%nsppol
2693            do ispinor=1,matlu(1)%nspinor
2694              do ispinor1=1,matlu(1)%nspinor
2695                matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=&
2696 &                matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)*fac
2697              enddo
2698            enddo
2699          enddo
2700        enddo
2701      enddo
2702    endif
2703  enddo
2704 
2705  end subroutine fac_matlu

m_matlu/gather_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 gather_matlu

FUNCTION

 Create new array from matlu

COPYRIGHT

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

INPUTS

  gloc(natom) <type(matlu_type)>        = density matrix in the spin spin representation
  gatherloc(natom) <type(coeff2c_type)> = density matrix where spin and angular momentum are gathered in the same index
  natom=number of atoms in cell.
  option= 1 go from gloc to gathergloc
  option= -1 go from gathergloc to gloc
  prtopt= option to define level of printing

OUTPUT

  gloc(natom) <type(matlu_type)>        = density matrix in the spin spin representation
  gatherloc(natom) <type(coeff2c_type)> = density matrix where spin and angular momentum are gathered in the same index

SIDE EFFECTS

SOURCE

1275  subroutine gather_matlu(gloc,gathergloc,natom,option,prtopt)
1276  use defs_basis
1277  use defs_wvltypes
1278  use m_crystal, only : crystal_t
1279  implicit none
1280 
1281 ! type  matlus_type
1282 !  SEQUENCE
1283 !  complex(dpc), pointer :: mat(:,:)
1284 ! end type matlus_type
1285 
1286 !Arguments ------------------------------------
1287 !scalars
1288  integer,intent(in) :: natom,option,prtopt
1289  type(coeff2c_type), intent(inout) :: gathergloc(natom)
1290  type(matlu_type),intent(inout) :: gloc(natom)
1291 !Local variables-------------------------------
1292 !scalars
1293  integer :: iatom,im1,im2,ispinor,ispinor1,isppol,isppol1
1294  integer :: jc1,jc2,ml1,ml2,ndim,nspinor,nsppol,tndim
1295  character(len=500) :: message
1296 
1297 ! DBG_ENTER("COLL")
1298  nsppol=gloc(1)%nsppol
1299  nspinor=gloc(1)%nspinor
1300 
1301  do iatom=1,natom
1302    if(gloc(iatom)%lpawu.ne.-1) then
1303 !==-------------------------------------
1304 
1305      ndim=2*gloc(iatom)%lpawu+1
1306      tndim=nsppol*nspinor*ndim
1307 
1308 !== Put norm into array "gathergloc"
1309      jc1=0
1310      do isppol=1,nsppol
1311        do ispinor=1,nspinor
1312          do ml1=1,ndim
1313            jc1=jc1+1
1314            jc2=0
1315            do isppol1=1,nsppol
1316              do ispinor1=1,nspinor
1317                do ml2=1,ndim
1318                  jc2=jc2+1
1319                  if(option==1) then
1320                    if(isppol==isppol1) then
1321                      gathergloc(iatom)%value(jc1,jc2)=gloc(iatom)%mat(ml1,ml2,isppol,ispinor,ispinor1)
1322                    endif
1323                  else if(option==-1) then
1324                    if(isppol==isppol1) then
1325                      gloc(iatom)%mat(ml1,ml2,isppol,ispinor,ispinor1)=gathergloc(iatom)%value(jc1,jc2)
1326                    endif
1327                  endif
1328                enddo
1329              enddo ! ispinor1
1330            enddo ! isppol1
1331          enddo
1332        enddo !ispinor
1333      enddo ! isppol
1334    endif
1335  enddo ! iatom
1336  if(option==1.and.prtopt==3) then
1337    do iatom=1,natom
1338      if(gloc(iatom)%lpawu.ne.-1) then
1339        tndim=nsppol*nspinor*(2*gloc(iatom)%lpawu+1)
1340        write(message,'(2a,i5)') ch10,' (gathermatlu:) For atom', iatom
1341        call wrtout(std_out,message,'COLL')
1342        do im1=1,tndim
1343          write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1344 &         (gathergloc(iatom)%value(im1,im2),im2=1,tndim)
1345          call wrtout(std_out,message,'COLL')
1346        end do
1347      endif
1348    enddo ! iatom
1349  else if(option==-1.and.prtopt==3) then
1350    call print_matlu(gloc,natom,prtopt)
1351  endif
1352 
1353 
1354 
1355 ! DBG_EXIT("COLL")
1356 
1357  end subroutine gather_matlu

m_matlu/identity_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 identity_matlu

FUNCTION

 Make the matlu the identity

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2806  subroutine identity_matlu(matlu,natom)
2807  use defs_basis
2808  implicit none
2809 
2810 !Arguments ------------------------------------
2811 !scalars
2812  integer, intent(in) :: natom
2813 !arrays
2814  type(matlu_type),intent(inout) :: matlu(natom)
2815 !Local variables-------------------------------
2816 !scalars
2817  integer :: iatom,im,ispinor,isppol
2818  integer :: ndim
2819 ! character(len=500) :: message
2820 !arrays
2821 !************************************************************************
2822  ! not yet tested and used
2823  do iatom=1,natom
2824    if(matlu(iatom)%lpawu.ne.-1) then
2825      ndim=2*matlu(iatom)%lpawu+1
2826      do isppol=1,matlu(iatom)%nsppol
2827        do im=1,ndim
2828          do ispinor=1,matlu(iatom)%nspinor
2829            matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)= cone
2830          enddo ! ispinor
2831        enddo ! im
2832      enddo
2833    endif ! lpawu
2834  enddo ! iatom
2835  end subroutine identity_matlu

m_matlu/init_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 init_matlu

FUNCTION

  Allocate variables used in type matlu_type.

INPUTS

  natom = number of atoms
  nspinor = number of spinorial components
  nsppol = number of polarisation components
  lpawu_natom(natom) = value of lpawu for all atoms
  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables

 OUTPUTS
  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables

SOURCE

140 subroutine init_matlu(natom,nspinor,nsppol,lpawu_natom,matlu)
141 
142  use defs_basis
143  use m_crystal, only : crystal_t
144  implicit none
145 
146 !Arguments ------------------------------------
147 !scalars
148  integer :: natom,nspinor,nsppol
149 !type
150  integer, intent(in) :: lpawu_natom(natom)
151  type(matlu_type), intent(inout) :: matlu(natom)
152 !Local variables ------------------------------------
153  integer :: iatom,lpawu
154 
155 !************************************************************************
156 
157 ! matlu%mband       = mband
158 ! matlu%dmftbandf   = dmftbandf
159 ! matlu%dmftbandi   = dmftbandi
160 ! matlu%nkpt        = nkpt
161 ! matlu%mbandc  = 0
162  matlu%nsppol      = nsppol
163  matlu%nspinor     = nspinor
164  do iatom=1,natom
165   lpawu=lpawu_natom(iatom)
166   matlu(iatom)%lpawu=lpawu
167   if(lpawu.ne.-1) then
168    ABI_MALLOC(matlu(iatom)%mat,(2*lpawu+1,2*lpawu+1,nsppol,nspinor,nspinor))
169    matlu(iatom)%mat=czero
170   else
171    ABI_MALLOC(matlu(iatom)%mat,(0,0,nsppol,nspinor,nspinor))
172   endif
173  enddo
174 
175 end subroutine init_matlu

m_matlu/inverse_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 inverse_matlu

FUNCTION

 Inverse local quantity.

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to inverse
  natom=number of atoms in cell.
  prtopt= option to define level of printing

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: inverse of input matrix

SIDE EFFECTS

NOTES

SOURCE

802  subroutine inverse_matlu(matlu,natom,prtopt)
803  use defs_basis
804  use defs_wvltypes
805  use m_crystal, only : crystal_t
806  implicit none
807 
808 !Arguments ------------------------------------
809 !scalars
810  integer, intent(in) :: natom
811  integer, intent(in) :: prtopt
812 !arrays
813  type(matlu_type),intent(inout) :: matlu(natom)
814 !Local variables-------------------------------
815  integer :: iatom,tndim
816  integer :: nsppol,nspinor
817 !scalars
818  type(coeff2c_type),allocatable :: gathermatlu(:)
819  !************************************************************************
820 
821 
822  nspinor=matlu(1)%nspinor
823  nsppol=matlu(1)%nsppol
824  if(prtopt>0) then
825  endif
826  ABI_MALLOC(gathermatlu,(natom))
827  do iatom=1,natom
828    if(matlu(iatom)%lpawu.ne.-1) then
829      tndim=nsppol*nspinor*(2*matlu(iatom)%lpawu+1)
830      ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
831      gathermatlu(iatom)%value=czero
832    endif
833  enddo
834 
835  call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
836  do iatom=1,natom
837    if(matlu(iatom)%lpawu.ne.-1) then
838      tndim=nsppol*nspinor*(2*matlu(iatom)%lpawu+1)
839      !call matcginv_dpc(gathermatlu(iatom)%value,tndim,tndim)
840      call xginv(gathermatlu(iatom)%value,tndim)
841    endif
842  enddo
843  call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
844 
845  do iatom=1,natom
846    if(matlu(iatom)%lpawu.ne.-1) then
847      ABI_FREE(gathermatlu(iatom)%value)
848    endif
849  enddo
850  ABI_FREE(gathermatlu)
851  end subroutine inverse_matlu

m_matlu/ln_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 ln_matlu

FUNCTION

 Compute the logarithm of matlu (only if diagonal for the moment)

COPYRIGHT

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

INPUTS

  matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity

SIDE EFFECTS

NOTES

SOURCE

2467  subroutine ln_matlu(matlu1,natom)
2468  use defs_basis
2469  use defs_wvltypes
2470  implicit none
2471 
2472 !Arguments ------------------------------------
2473 !scalars
2474  integer, intent(in) :: natom
2475 !arrays
2476  type(matlu_type), intent(inout) :: matlu1(natom)
2477 !Local variables-------------------------------
2478 !scalars
2479  integer :: iatom,im,ispinor,isppol
2480  integer :: lpawu
2481  character(len=500) :: message
2482 !arrays
2483 !************************************************************************
2484  !call checkdiag_matlu(matlu1,natom,tol8)
2485  do iatom=1,natom
2486    lpawu=matlu1(iatom)%lpawu
2487    if(lpawu.ne.-1) then
2488      do isppol=1,matlu1(1)%nsppol
2489        do ispinor=1,matlu1(1)%nspinor
2490          do im=1,2*lpawu+1
2491            if( real(matlu1(iatom)%mat(im,im,isppol,ispinor,ispinor))<zero) then
2492              write(message,'(2a,2es13.5,a)') ch10," ln_matlu: PROBLEM " &
2493 &             , matlu1(iatom)%mat(im,im,isppol,ispinor,ispinor)
2494              ABI_ERROR(message)
2495            endif
2496            matlu1(iatom)%mat(im,im,isppol,ispinor,ispinor)= &
2497 &           log(matlu1(iatom)%mat(im,im,isppol,ispinor,ispinor))
2498          enddo ! im
2499        enddo ! ispinor
2500      enddo ! isppol
2501    endif ! lpawu
2502  enddo ! iatom
2503 
2504  end subroutine ln_matlu

m_matlu/magmomforb_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 magmomforb_matlu

FUNCTION

 return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
 with the matrix of orbital angular momentum element for the x,y and z direction. 
 Option gives the direction of the magnetic moment x==1, y==2 and z==3. 

COPYRIGHT

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

INPUTS

 matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
 natom :: number of atoms
 option = 1 :: x axis
        = 2 :: y axis
        = 3 :: z axis
 optptr > 2 :: print orbital angular matrix elements and resulting product

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product 

SIDE EFFECTS

NOTES

SOURCE

2871  subroutine magmomforb_matlu(matlu,mu,natom,option,optprt)
2872  use defs_basis
2873  use defs_wvltypes
2874  implicit none
2875 
2876 !Arguments ------------------------------------
2877 !scalars
2878  integer, intent(in) :: natom,option,optprt
2879  complex(dpc), allocatable, intent(out) :: mu
2880 !arrays
2881  type(matlu_type), intent(inout) :: matlu(natom)
2882 !Local variables-------------------------------
2883 !scalars
2884  integer :: iatom,im,ispinor,isppol,ispinor2
2885  integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim
2886  character(len=500) :: message
2887  real(dp) :: xj
2888 !arrays
2889  complex(dpc),allocatable :: mat_out_c(:,:)
2890 ! integer, allocatable :: ind_msml(:,:)
2891  complex(dpc), allocatable :: temp_mat(:,:)
2892  type(coeff2c_type), allocatable :: gathermatlu(:)
2893  type(coeff2c_type), allocatable :: muorb(:)
2894 !************************************************************************
2895 
2896  !=====================================
2897  ! Allocate matrices 
2898  !=====================================
2899     
2900  ABI_MALLOC(gathermatlu,(natom))
2901  ABI_MALLOC(muorb,(natom))
2902  do iatom=1,natom
2903    if(matlu(iatom)%lpawu.ne.-1) then
2904      tndim=2*(2*matlu(iatom)%lpawu+1)
2905      ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
2906      gathermatlu(iatom)%value=czero
2907      ABI_MALLOC(muorb(iatom)%value,(tndim,tndim))
2908      muorb(iatom)%value=czero
2909    end if
2910  end do 
2911 
2912  do iatom=1,natom
2913    lpawu=matlu(iatom)%lpawu
2914    if(lpawu.ne.-1) then
2915      ll=lpawu
2916      lcor=lpawu
2917      !=====================================
2918      !build orbital angular momentum matrix along x axis
2919      !=====================================
2920      if(option==1) then
2921 
2922        jc1=0
2923        do ms1 =-1,1
2924          xj=float(ms1)+half
2925          do ml1 = -ll,ll
2926             jc1=jc1+1
2927             if(jc1 == 1) then
2928                muorb(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
2929             endif
2930             if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
2931                muorb(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
2932                muorb(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
2933             else if(jc1== 2*(2*ll+1)) then
2934                muorb(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
2935             end if
2936           end do
2937         end do
2938 
2939      !=====================================
2940      !build orbital angular momentum matrix along y axis
2941      !=====================================
2942      else if(option==2) then
2943 
2944        jc1=0
2945        do ms1 =-1,1
2946          xj=float(ms1)+half
2947          do ml1 = -ll,ll
2948             jc1=jc1+1
2949             if(jc1 == 1) then
2950                muorb(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
2951             endif
2952             if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
2953                muorb(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
2954                muorb(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
2955             else if(jc1 == 2*(2*ll+1)) then
2956                muorb(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
2957             end if
2958           end do
2959         end do
2960 
2961      !=====================================
2962      !build orbital angular momentum matrix along z axis
2963      !=====================================
2964      else if(option==3) then
2965        jc1=0
2966        do ms1=-1,1
2967          do ml1=-ll,ll
2968             jc1=jc1+1
2969             if(jc1 < tndim+1) then
2970               muorb(iatom)%value(jc1,jc1) = ml1
2971             endif
2972           end do
2973         end do
2974      end if
2975 
2976      if(optprt>2) then
2977         write(message,'(a,i4)') "Orbital angular momentum matrix elements in |m_l,m_s> basis for axis=", option
2978         call wrtout(std_out,message,"COLL")
2979         do im=1,2*(ll*2+1)
2980           write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muorb(iatom)%value(im,jm),jm=1,2*(ll*2+1))
2981           call wrtout(std_out,message,"COLL")
2982         end do
2983      end if
2984 
2985    end if !lpawu
2986  end do !atom
2987 
2988      !=====================================
2989      ! Reshape input Ylm matlu in one 14x14 matrix 
2990      !=====================================
2991      
2992  call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
2993 
2994 !!printing for debug
2995 ! write(std_out,*) "gathermatlu in magmomforb"
2996 ! do im1=1,tndim
2997 !   write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2998 !        (gathermatlu(1)%value(im1,im2),im2=1,tndim)
2999 !   call wrtout(std_out,message,'COLL')
3000 ! end do
3001 
3002      !=====================================
3003      ! Matrix product of Occ and muorb 
3004      !=====================================
3005  do iatom=1,natom
3006    if(matlu(iatom)%lpawu.ne.-1) then
3007      tndim=2*(2*matlu(iatom)%lpawu+1)
3008      ABI_MALLOC(temp_mat,(tndim,tndim))
3009 
3010      call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,muorb(iatom)%value,tndim,czero,temp_mat,tndim)
3011      
3012      gathermatlu(iatom)%value=temp_mat
3013      ABI_FREE(temp_mat)
3014 
3015      !=====================================
3016      ! Trace of matrix product
3017      !=====================================
3018 
3019    mu=czero
3020    do im1=1,tndim
3021      do im2=1,tndim
3022        if(im1==im2) then
3023          mu = mu + gathermatlu(iatom)%value(im1,im2)
3024        end if
3025      end do
3026    end do
3027 
3028      !=====================================
3029      ! Reshape product matrix into matlu format 
3030      !=====================================
3031 
3032      call gather_matlu(matlu(iatom),gathermatlu(iatom),iatom,option=-1,prtopt=1)
3033 
3034      !=====================================
3035      ! Print matlu
3036      !=====================================
3037 
3038      if(optprt>2) then
3039        ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1)) 
3040        do isppol=1,matlu(1)%nsppol
3041          do ispinor=1,matlu(1)%nspinor
3042            do ispinor2=1,matlu(1)%nspinor
3043              mat_out_c(:,:) = matlu(iatom)%mat(:,:,isppol,ispinor,ispinor2)
3044 
3045              write(message,'(2a, i2, a, i2, a, i2)') ch10,"Orbital angular momentum matrix, isppol=", isppol, ", ispinor=",&
3046 &            ispinor,", ispinor2=", ispinor2
3047              call wrtout(std_out,message,'COLL')
3048              do im1=1,ll*2+1
3049                write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3050       &         (mat_out_c(im1,im2),im2=1,ll*2+1)
3051                call wrtout(std_out,message,'COLL')
3052              end do
3053 
3054            end do ! ispinor2
3055          end do ! ispinor
3056        end do ! isppol
3057        ABI_FREE(mat_out_c)
3058      endif
3059 
3060    end if !lpawu
3061  end do !atom 
3062 
3063      !=====================================
3064      ! Deallocate gathermatlu 
3065      !=====================================
3066 
3067  do iatom=1,natom
3068    if(matlu(iatom)%lpawu.ne.-1) then
3069      ABI_FREE(gathermatlu(iatom)%value)
3070      ABI_FREE(muorb(iatom)%value)
3071    end if
3072  end do
3073  ABI_FREE(gathermatlu)
3074  ABI_FREE(muorb)
3075 
3076  end subroutine magmomforb_matlu

m_matlu/magmomfspin_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 magmomfspin_matlu

FUNCTION

 return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
 with the matrix of spin angular momentum element for the x,y and z direction. 
 Option gives the direction of the magnetic moment x==1, y==2 and z==3. 

COPYRIGHT

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

INPUTS

 matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
 natom :: number of atoms
 option = 1 :: x axis
        = 2 :: y axis
        = 3 :: z axis
 optptr > 2 :: print spin angular matrix elements and resulting product

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product

SIDE EFFECTS

NOTES

SOURCE

3114  subroutine magmomfspin_matlu(matlu,mu,natom,option,optprt)
3115  use defs_basis
3116  use defs_wvltypes
3117  implicit none
3118 
3119 !Arguments ------------------------------------
3120 !scalars
3121  integer, intent(in) :: natom,option,optprt
3122  complex(dpc), allocatable, intent(out) :: mu
3123 !arrays
3124  type(matlu_type), intent(inout) :: matlu(natom)
3125 !Local variables-------------------------------
3126 !scalars
3127  integer :: iatom,im,ispinor,isppol,ispinor2
3128  integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim
3129  character(len=500) :: message
3130  real(dp) :: xj
3131 !arrays
3132  complex(dpc),allocatable :: mat_out_c(:,:)
3133  integer, allocatable :: ind_msml(:,:)
3134  complex(dpc), allocatable :: temp_mat(:,:)
3135  type(coeff2c_type), allocatable :: gathermatlu(:)
3136  type(coeff2c_type), allocatable :: muspin(:)
3137 !************************************************************************
3138 
3139  !=====================================
3140  ! Allocate matrices 
3141  !=====================================
3142     
3143  ABI_MALLOC(gathermatlu,(natom))
3144  ABI_MALLOC(muspin,(natom))
3145  do iatom=1,natom
3146    if(matlu(iatom)%lpawu.ne.-1) then
3147      tndim=2*(2*matlu(iatom)%lpawu+1)
3148      ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
3149      gathermatlu(iatom)%value=czero
3150      ABI_MALLOC(muspin(iatom)%value,(tndim,tndim))
3151      muspin(iatom)%value=czero
3152    end if
3153  end do 
3154 
3155  do iatom=1,natom
3156    lpawu=matlu(iatom)%lpawu
3157    if(lpawu.ne.-1) then
3158      ll=lpawu
3159      lcor=lpawu
3160      ABI_MALLOC(ind_msml,(2,-ll:ll))
3161      ind_msml=czero
3162      !=====================================
3163      !build spin angular momentum matrix along x axis
3164      !=====================================
3165      if(option==1) then
3166        jc1=0
3167        do ms1=1,2
3168          do ml1=-ll,ll
3169           jc1=jc1+1
3170           ind_msml(ms1,ml1)=jc1
3171          end do
3172        end do
3173 
3174        jc1=0
3175        do ms1 =-1,1
3176          xj=float(ms1)+half 
3177          do ml1 = -ll,ll
3178             jc1=jc1+1
3179             if(xj < 0.0 ) then
3180                muspin(iatom)%value(ind_msml(2,ml1),jc1) = 0.5
3181             else if(xj > 0.0) then
3182                muspin(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 0.5
3183             end if
3184           end do
3185         end do
3186 
3187      !=====================================
3188      !build spin angular momentum matrix along y axis
3189      !up spin is first
3190      !=====================================
3191      else if(option==2) then
3192         jc1=0
3193        do ms1=1,2
3194          do ml1=-ll,ll
3195           jc1=jc1+1
3196           ind_msml(ms1,ml1)=jc1
3197          end do
3198        end do
3199 
3200        jc1=0
3201        do ms1 =-1,1
3202          xj=float(ms1)+half 
3203          do ml1 = -ll,ll
3204             jc1=jc1+1
3205             if(xj < 0.0 ) then
3206                muspin(iatom)%value(ind_msml(2,ml1),jc1) = cmplx(zero,0.5,kind=dp)
3207             else if(xj > 0.0) then
3208                muspin(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = cmplx(zero,-0.5,kind=dp)
3209             end if
3210           end do
3211         end do
3212 
3213      !=====================================
3214      !build spin angular momentum matrix along z axis
3215      !up spin is first
3216      !=====================================
3217      else if(option==3) then
3218        jc1=0
3219        do ms1=-1,1
3220          xj=float(ms1)+half 
3221          do ml1=-ll,ll
3222             jc1=jc1+1
3223             if(jc1 < tndim+1) then
3224               if(xj < 0.0 ) then
3225                  muspin(iatom)%value(jc1,jc1) = -xj
3226               else if(xj > 0.0) then
3227                  muspin(iatom)%value(jc1,jc1) = -xj
3228               end if
3229             endif
3230          end do
3231        end do
3232      end if
3233 
3234      if(optprt>2) then
3235         write(message,'(a,i4)') "Spin angular momentum matrix elements in |m_l,m_s> basis for axis", option
3236         call wrtout(std_out,message,"COLL")
3237         do im=1,2*(ll*2+1)
3238           write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muspin(iatom)%value(im,jm),jm=1,2*(ll*2+1))
3239           call wrtout(std_out,message,"COLL")
3240         end do
3241      end if
3242 
3243    ABI_FREE(ind_msml) 
3244    end if !lpawu
3245  end do !atom
3246 
3247      !=====================================
3248      ! Reshape input Ylm matlu in one 14x14 matrix 
3249      !=====================================
3250      
3251  call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
3252 
3253 !!printing for debug
3254 !! write(std_out,*) "gathermatlu in magmomfspin"
3255 !! do im1=1,tndim
3256 !!   write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3257 !!        (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3258 !!   call wrtout(std_out,message,'coll')
3259 !! end do
3260 
3261      !=====================================
3262      ! Matrix product of Occ and muspin 
3263      !=====================================
3264  do iatom=1,natom
3265    if(matlu(iatom)%lpawu.ne.-1) then
3266      tndim=2*(2*matlu(iatom)%lpawu+1)
3267      ABI_MALLOC(temp_mat,(tndim,tndim))
3268 
3269      call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,muspin(iatom)%value,tndim,czero,temp_mat,tndim)
3270      
3271      gathermatlu(iatom)%value=temp_mat
3272      ABI_FREE(temp_mat)
3273 
3274      !!printing for debug
3275      !!write(std_out,*) "gathermatlu in magmomfspin after product"
3276      !!do im1=1,tndim
3277      !!   write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3278      !!        (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3279      !!   call wrtout(std_out,message,'coll')
3280      !!end do
3281 
3282      !=====================================
3283      ! Trace of matrix product
3284      !=====================================
3285 
3286    mu=czero
3287    do im1=1,tndim
3288      do im2=1,tndim
3289        if(im1==im2) then
3290          mu = mu + gathermatlu(iatom)%value(im1,im2)
3291        end if
3292      end do
3293    end do
3294 
3295      !=====================================
3296      ! Reshape product matrix into matlu format 
3297      !=====================================
3298 
3299     call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
3300 
3301 
3302      !=====================================
3303      ! Print matlu
3304      !=====================================
3305      if(optprt>2) then
3306        ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
3307        do isppol=1,matlu(1)%nsppol
3308          do ispinor=1,matlu(1)%nspinor
3309            do ispinor2=1,matlu(1)%nspinor
3310              mat_out_c(:,:) = matlu(iatom)%mat(:,:,isppol,ispinor,ispinor2)
3311 
3312              write(message,'(2a, i2, a, i2, a, i2)') ch10,"Spin angular momentum matrix, isppol=", isppol, ", ispinor=", ispinor,&
3313 &             ", ispinor2=", ispinor2
3314              call wrtout(std_out,message,'COLL')
3315              do im1=1,ll*2+1
3316                write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3317       &         (mat_out_c(im1,im2),im2=1,ll*2+1)
3318                call wrtout(std_out,message,'COLL')
3319              end do
3320 
3321            end do ! im
3322          end do ! ispinor
3323        end do ! isppol
3324        ABI_FREE(mat_out_c)
3325      endif
3326 
3327    end if !lpawu
3328  end do !atom
3329 
3330      !=====================================
3331      ! Deallocate gathermatlu 
3332      !=====================================
3333 
3334  do iatom=1,natom
3335    if(matlu(iatom)%lpawu.ne.-1) then
3336      ABI_FREE(gathermatlu(iatom)%value)
3337      ABI_FREE(muspin(iatom)%value)
3338    end if
3339  end do
3340  ABI_FREE(gathermatlu)
3341  ABI_FREE(muspin)
3342  
3343  end subroutine magmomfspin_matlu

m_matlu/magmomfzeeman_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 magmomfspin_matlu

FUNCTION

 return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
 with the matrix of Zeeman angular momentum element (L_u + 2*S_u) for the x,y and z direction. 
 Option gives the direction of the magnetic moment x==1, y==2 and z==3. 

COPYRIGHT

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

INPUTS

 matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
 natom :: number of atoms
 option = 1 :: x axis 
        = 2 :: y axis 
        = 3 :: z axis 
 optptr > 2 :: print Zeeman angular matrix elements and resulting product

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product

SIDE EFFECTS

NOTES

SOURCE

3381  subroutine magmomfzeeman_matlu(matlu,mu,natom,option,optprt)
3382  use defs_basis
3383  use defs_wvltypes
3384  implicit none
3385 
3386 !Arguments ------------------------------------
3387 !scalars
3388  integer, intent(in) :: natom,option,optprt
3389  complex(dpc), allocatable, intent(out) :: mu
3390 !arrays
3391  type(matlu_type), intent(inout) :: matlu(natom)
3392 !Local variables-------------------------------
3393 !scalars
3394  integer :: iatom,im,ispinor,isppol,ispinor2
3395  integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim
3396  character(len=500) :: message
3397  real(dp) :: xj
3398 !arrays
3399  complex(dpc),allocatable :: mat_out_c(:,:)
3400  integer, allocatable :: ind_msml(:,:)
3401  complex(dpc), allocatable :: temp_mat(:,:)
3402  type(coeff2c_type), allocatable :: gathermatlu(:)
3403  type(coeff2c_type), allocatable :: muzeeman(:)
3404 !************************************************************************
3405 
3406  !=====================================
3407  ! Allocate matrices 
3408  !=====================================
3409     
3410  ABI_MALLOC(gathermatlu,(natom))
3411  ABI_MALLOC(muzeeman,(natom))
3412  do iatom=1,natom
3413    if(matlu(iatom)%lpawu.ne.-1) then
3414      tndim=2*(2*matlu(iatom)%lpawu+1)
3415      ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
3416      gathermatlu(iatom)%value=czero
3417      ABI_MALLOC(muzeeman(iatom)%value,(tndim,tndim))
3418      muzeeman(iatom)%value=czero
3419    end if
3420  end do 
3421 
3422  do iatom=1,natom
3423    lpawu=matlu(iatom)%lpawu
3424    if(lpawu.ne.-1) then
3425      ll=lpawu
3426      lcor=lpawu
3427      ABI_MALLOC(ind_msml,(2,-ll:ll))
3428      ind_msml=czero
3429      !=====================================
3430      !build Zeeman angular momentum matrix along x axis
3431      !=====================================
3432      if(option==1) then
3433        jc1=0
3434        do ms1=1,2
3435          do ml1=-ll,ll
3436           jc1=jc1+1
3437           ind_msml(ms1,ml1)=jc1
3438          end do
3439        end do
3440 
3441        jc1=0
3442        do ms1 =-1,1
3443          xj=float(ms1)+half 
3444          do ml1 = -ll,ll
3445             jc1=jc1+1
3446             if(xj < 0.0 ) then
3447                muzeeman(iatom)%value(ind_msml(2,ml1),jc1) = 2*0.5
3448             else if(xj > 0.0) then
3449                muzeeman(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 2*0.5
3450             end if
3451             if(jc1 == 1) then
3452                muzeeman(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3453             endif
3454             if(jc1 > 1 .and.jc1 < 2*(2*ll+1)) then
3455                muzeeman(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3456                muzeeman(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3457             else if(jc1 == 2*(2*ll+1)) then
3458                muzeeman(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3459             end if
3460           end do
3461         end do
3462 
3463      !=====================================
3464      !build Zeeman angular momentum matrix along y axis
3465      !up spin is first
3466      !=====================================
3467      else if(option==2) then
3468         jc1=0
3469        do ms1=1,2
3470          do ml1=-ll,ll
3471           jc1=jc1+1
3472           ind_msml(ms1,ml1)=jc1
3473          end do
3474        end do
3475 
3476        jc1=0
3477        do ms1 =-1,1
3478          xj=float(ms1)+half 
3479          do ml1 = -ll,ll
3480             jc1=jc1+1
3481             if(xj < 0.0 ) then
3482                muzeeman(iatom)%value(ind_msml(2,ml1),jc1) = 2*cmplx(zero,0.5,kind=dp)
3483             else if(xj > 0.0) then
3484                muzeeman(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 2*cmplx(zero,-0.5,kind=dp)
3485             end if
3486             if(jc1 == 1) then
3487                muzeeman(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
3488             endif
3489             if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
3490                muzeeman(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
3491                muzeeman(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
3492             else if(jc1== 2*(2*ll+1)) then
3493                muzeeman(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp) 
3494             end if
3495           end do
3496         end do
3497 
3498 
3499      !=====================================
3500      !build Zeeman angular momentum matrix along z axis
3501      !up spin is first
3502      !=====================================
3503      else if(option==3) then
3504        jc1=0
3505        do ms1=-1,1
3506          xj=float(ms1)+half 
3507          do ml1=-ll,ll
3508             jc1=jc1+1
3509             if(jc1 < tndim+1) then
3510               if(xj < 0.0 ) then
3511                  muzeeman(iatom)%value(jc1,jc1) = ml1-2*xj
3512               else if(xj > 0.0) then
3513                  muzeeman(iatom)%value(jc1,jc1) = ml1-2*xj
3514               end if
3515             endif
3516          end do
3517        end do
3518      end if
3519 
3520 
3521      if(optprt>2) then
3522         write(message,'(a,i4)') "Zeeman angular momentum matrix elements in |m_l,m_s> basis for axis", option
3523         call wrtout(std_out,message,"COLL")
3524         do im=1,2*(ll*2+1)
3525           write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muzeeman(iatom)%value(im,jm),jm=1,2*(ll*2+1))
3526           call wrtout(std_out,message,"COLL")
3527         end do
3528      end if
3529 
3530    ABI_FREE(ind_msml) 
3531    end if !lpawu
3532  end do !atom
3533 
3534      !=====================================
3535      ! Reshape input Ylm matlu in one 14x14 matrix 
3536      !=====================================
3537     
3538  !ABI_MALLOC(gathermatlu,(natom))
3539  !do iatom=1,natom
3540  !  if(matlu(iatom)%lpawu.ne.-1) then
3541  !    tndim=2*(2*matlu(iatom)%lpawu+1)
3542  !    ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
3543  !    gathermatlu(iatom)%value=czero
3544  !  end if
3545  !end do 
3546      
3547  call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
3548 
3549 !!printing for debug
3550 !! write(std_out,*) "gathermatlu in magmomfspin"
3551 !! do im1=1,tndim
3552 !!   write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3553 !!        (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3554 !!   call wrtout(std_out,message,'coll')
3555 !! end do
3556 
3557      !=====================================
3558      ! Matrix product of Occ and muzeeman 
3559      !=====================================
3560  do iatom=1,natom
3561    if(matlu(iatom)%lpawu.ne.-1) then
3562      tndim=2*(2*matlu(iatom)%lpawu+1)
3563      ABI_MALLOC(temp_mat,(tndim,tndim))
3564 
3565      call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,muzeeman(iatom)%value,tndim,czero,temp_mat,tndim)
3566      
3567      gathermatlu(iatom)%value=temp_mat
3568      ABI_FREE(temp_mat)
3569 
3570      !!printing for debug
3571      !!write(std_out,*) "gathermatlu in magmomfspin after product"
3572      !!do im1=1,tndim
3573      !!   write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3574      !!        (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3575      !!   call wrtout(std_out,message,'coll')
3576      !!end do
3577 
3578      !=====================================
3579      ! Trace of matrix product
3580      !=====================================
3581 
3582    mu=czero
3583    do im1=1,tndim
3584      do im2=1,tndim
3585        if(im1==im2) then
3586          mu = mu + gathermatlu(iatom)%value(im1,im2)
3587        end if
3588      end do
3589    end do
3590 
3591      !=====================================
3592      ! Reshape product matrix into matlu format 
3593      !=====================================
3594 
3595     call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
3596 
3597 
3598      !=====================================
3599      ! Print matlu
3600      !=====================================
3601      if(optprt>2) then
3602        ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
3603        do isppol=1,matlu(1)%nsppol
3604          do ispinor=1,matlu(1)%nspinor
3605            do ispinor2=1,matlu(1)%nspinor
3606              mat_out_c(:,:) = matlu(iatom)%mat(:,:,isppol,ispinor,ispinor2)
3607 
3608              write(message,'(2a, i2, a, i2, a, i2)') ch10,"Zeeman angular momentum matrix, isppol=", isppol, ", ispinor=",&
3609 &            ispinor,", ispinor2=", ispinor2
3610              call wrtout(std_out,message,'COLL')
3611              do im1=1,ll*2+1
3612                write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3613       &         (mat_out_c(im1,im2),im2=1,ll*2+1)
3614                call wrtout(std_out,message,'COLL')
3615              end do
3616 
3617            end do ! im
3618          end do ! ispinor
3619        end do ! isppol
3620        ABI_FREE(mat_out_c)
3621      endif ! optprt
3622 
3623    end if !lpawu
3624  end do !atom
3625 
3626      !=====================================
3627      ! Deallocate gathermatlu 
3628      !=====================================
3629 
3630  do iatom=1,natom
3631    if(matlu(iatom)%lpawu.ne.-1) then
3632      ABI_FREE(gathermatlu(iatom)%value)
3633      ABI_FREE(muzeeman(iatom)%value)
3634    end if
3635  end do
3636  ABI_FREE(gathermatlu)
3637  ABI_FREE(muzeeman)
3638  
3639  end subroutine magmomfzeeman_matlu

m_matlu/matlu_type [ Types ]

[ Top ] [ m_matlu ] [ Types ]

NAME

  matlu_type

FUNCTION

  This structured datatype contains an matrix for the correlated subspace

SOURCE

 83  type, public :: matlu_type ! for each atom
 84 
 85   integer :: lpawu
 86   ! value of the angular momentum for correlated electrons
 87 
 88 !  integer :: natom
 89    ! number of atoms (given for each atom, not useful..could be changed)
 90 !
 91 !  integer :: mband
 92 !  ! Number of bands
 93 !
 94 !  integer :: mbandc
 95 !  ! Total number of bands in the Kohn-Sham Basis for PAW+DMFT
 96 !
 97 !  integer :: natpawu         ! Number of correlated atoms
 98 !
 99 !  integer :: nkpt
100 !  ! Number of k-point in the IBZ.
101   character(len=12) :: whichmatlu
102   ! describe the type of local matrix computed (greenDFT, etc..)
103 !
104   integer :: nspinor
105   ! Number of spinorial component
106 !
107   integer :: nsppol
108   ! Number of polarization
109 
110   complex(dpc), allocatable :: mat(:,:,:,:,:)
111   ! local quantity
112 
113  end type matlu_type
114 
115 !----------------------------------------------------------------------
116 
117 
118 CONTAINS  !========================================================================================

m_matlu/print_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 print_matlu

FUNCTION

INPUTS

  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables
  natom= number of atoms
  prtopt= option for printing
  opt_diag=   0   print non diagonal matrix (real or complex according to nspinor)
             -1   print non diagonal complex matrix
            >=1   print diagonal matrix (real or complex according to nspinor)
  opt_ab_out=  0  print matrix on std_out
             /=0  print matrix on ab_out

OUTPUT

SOURCE

367 subroutine print_matlu(matlu,natom,prtopt,opt_diag,opt_ab_out,opt_exp,argout,compl)
368 
369  use defs_basis
370  use m_crystal, only : crystal_t
371  implicit none
372 
373 !Arguments ------------------------------------
374 !type
375  integer, intent(in):: natom
376  type(matlu_type),intent(in) :: matlu(natom)
377  integer, intent(in) :: prtopt
378  integer, optional, intent(in) :: opt_diag,opt_ab_out,opt_exp,argout,compl
379 !Local variables-------------------------------
380  integer :: iatom,ispinor,ispinor1,isppol,m1,m,lpawu,nspinor,nsppol,optdiag,optab_out,arg_out
381  character(len=500) :: message
382  character(len=4) :: mode_paral
383  type(matlu_type), allocatable :: matlu_tmp(:)
384  character(len=9),parameter :: dspinm(2,2)= RESHAPE((/"n        ","mx       ","my       ","mz       "/),(/2,2/))
385  logical :: testcmplx
386  real(dp) :: noccspin
387 ! *********************************************************************
388  mode_paral='COLL'
389  if(present(opt_diag)) then
390    optdiag=opt_diag
391  else
392    optdiag=0
393  endif
394  if(present(opt_ab_out)) then
395    optab_out=opt_ab_out
396  else
397    optab_out=0
398  endif
399  if(optab_out==0) then
400    arg_out=std_out
401  else
402    arg_out=ab_out
403  endif
404  if(present(argout)) then
405   arg_out=argout
406   mode_paral='PERS'
407  endif
408  nspinor=matlu(1)%nspinor
409  nsppol=matlu(1)%nsppol
410  testcmplx=(nspinor==2)
411  if(present(compl)) testcmplx=(nspinor==2).or.(compl==1)
412 
413  do iatom = 1 , natom
414    lpawu=matlu(iatom)%lpawu
415    if(lpawu/=-1) then
416      write(message,'(2a,i4)')  ch10,'   -------> For Correlated Atom', iatom
417      call wrtout(arg_out,  message,mode_paral)
418      do isppol = 1 , nsppol
419        if(present(opt_ab_out).and.nsppol==2) then
420          noccspin=zero
421          do m1=1, 2*lpawu +1
422            noccspin=noccspin+REAL(matlu(iatom)%mat(m1,m1,isppol,1,1))
423          enddo
424          !write(message,fmt='(7x,a,i3,a,f10.5)') ". Occ. for lpawu and for spin",isppol," =",noccspin
425          !call wrtout(arg_out, message,mode_paral)
426        endif
427      enddo
428 
429      do isppol = 1 , nsppol
430        if(nspinor == 1) then
431          write(message,'(a,10x,a,i3,i3)')  ch10,'-- polarization spin component',isppol
432          call wrtout(arg_out,  message,mode_paral)
433        endif
434        do ispinor = 1 , nspinor
435          do ispinor1 = 1, nspinor
436            if(nspinor == 2) then
437              write(message,'(a,10x,a,i3,i3)')  ch10,'-- spin components',ispinor,ispinor1
438              call wrtout(arg_out,  message,mode_paral)
439            endif
440            if(optdiag<=0) then
441              do m1=1, 2*lpawu+1
442                if(optdiag==0) then
443                  if(nspinor==1.and.abs(prtopt)>0)  then
444                    if(present(opt_exp)) then
445                      write(message,'(5x,20e24.14)') (REAL(matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
446 !                     call wrtout(arg_out,  message,mode_paral)
447 !                     write(message,'(5x,20e20.14)') (REAL(sqrt(matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1))),m=1,2*lpawu+1)
448 !                     call wrtout(arg_out,  message,mode_paral)
449 !                     write(message,'(5x,20e20.14)') (REAL(1.d0/sqrt(matlu(iatom)%mat(m,m,isppol,ispinor,ispinor1))),m=1,2*lpawu+1)
450                    else
451                      write(message,'(5x,20f10.5)') (REAL(matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
452                    endif
453                  endif
454                  if(testcmplx.and.abs(prtopt)>0) then
455                    if(present(opt_exp)) then
456                      if(opt_exp==2) then
457                        write(message,'(5x,14(2e18.10,1x))') ((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
458                      else
459                        write(message,'(5x,14(2e14.4,2x))') ((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
460                      endif
461                    else
462                      write(message,'(5x,14(2f9.5,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
463                    endif
464 !&                  write(message,'(5x,14(2f15.11,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
465                  endif
466                else if(optdiag==-1) then
467                  write(message,'(5x,14(2f10.5,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
468                endif
469                call wrtout(arg_out,  message,mode_paral)
470              enddo
471            elseif (optdiag>=1) then
472              if(nspinor==1.and.abs(prtopt)>0) &
473 &             write(message,'(5x,20f10.5)') (REAL(matlu(iatom)%mat(m,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
474              if(testcmplx.and.abs(prtopt)>0) &
475 &             write(message,'(5x,14(2f9.5,2x))')((matlu(iatom)%mat(m,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
476 !            write(std_out,'(5x,14(2f9.5,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
477              call wrtout(arg_out,  message,mode_paral)
478            endif
479          end do ! ispinor1
480        end do ! ispinor
481        if(nspinor==2.and.prtopt>=5) then
482          ABI_MALLOC(matlu_tmp,(0:natom))
483          call init_matlu(natom,nspinor,nsppol,matlu(1:natom)%lpawu,matlu_tmp)
484          matlu_tmp(iatom)%mat(m1,m,isppol,1,1)= matlu(iatom)%mat(m1,m,isppol,1,1)+matlu(iatom)%mat(m1,m,isppol,2,2)
485          matlu_tmp(iatom)%mat(m1,m,isppol,2,2)= matlu(iatom)%mat(m1,m,isppol,1,1)+matlu(iatom)%mat(m1,m,isppol,2,2)
486          matlu_tmp(iatom)%mat(m1,m,isppol,1,2)= matlu(iatom)%mat(m1,m,isppol,1,2)+matlu(iatom)%mat(m1,m,isppol,2,1)
487          matlu_tmp(iatom)%mat(m1,m,isppol,2,1)= &
488 &           (matlu(iatom)%mat(m1,m,isppol,1,2)+matlu(iatom)%mat(m1,m,isppol,2,1))*cmplx(0_dp,1_dp,kind=dp)
489          do ispinor = 1 , nspinor
490            do ispinor1 = 1, nspinor
491              write(message,'(a,10x,a,a)')  ch10,'-- spin components',dspinm(ispinor,ispinor1)
492              call wrtout(arg_out,  message,mode_paral)
493              write(message,'(5x,14(2f9.5,2x))')((matlu_tmp(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
494              call wrtout(arg_out,  message,mode_paral)
495            end do ! ispinor1
496          end do ! ispinor
497          call destroy_matlu(matlu_tmp,natom)
498          ABI_FREE(matlu_tmp)
499        endif
500      enddo ! isppol
501 !     if(nsppol==1.and.nspinor==1) then
502 !       write(message,'(a,10x,a,i3,a)')  ch10,'-- polarization spin component',isppol+1,' is identical'
503 !       call wrtout(arg_out,  message,mode_paral)
504 !     endif
505    endif ! lpawu/=1
506  enddo ! natom
507 
508 
509 end subroutine print_matlu

m_matlu/printplot_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 printplot_matlu

FUNCTION

 shift matlu matrix

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.
  shift= shift of the diagonal part.

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: shifted matrix

SIDE EFFECTS

NOTES

SOURCE

2734  subroutine printplot_matlu(matlu,natom,freq,char1,units,imre)
2735  use defs_basis
2736  use m_fstrings,       only : int2char4
2737  implicit none
2738 
2739 !Arguments ------------------------------------
2740 !scalars
2741  integer, intent(in) :: natom,units
2742  integer, optional, intent(in) :: imre
2743  real(dp), intent(in) :: freq
2744 !arrays
2745  type(matlu_type),intent(inout) :: matlu(natom)
2746  character(len=*), intent(in) :: char1
2747 !Local variables-------------------------------
2748 !scalars
2749  integer :: iatom,im,im1,ispinor,ispinor1,isppol
2750  integer :: lpawu,unitnb
2751  character(len=4) :: tag_at
2752  character(len=fnlen) :: tmpfil,tmpfilre,tmpfilim
2753 ! character(len=500) :: message
2754 !arrays
2755 !************************************************************************
2756  ! not yet tested and used
2757  do iatom=1,natom
2758    lpawu=matlu(iatom)%lpawu
2759    if(lpawu.ne.-1) then
2760      unitnb=units+iatom
2761      call int2char4(iatom,tag_at)
2762      if(present(imre)) then
2763        tmpfilre = trim(char1)//tag_at//"re"
2764        tmpfilim = trim(char1)//tag_at//"im"
2765        open (unit=unitnb+10,file=trim(tmpfilre),status='unknown',form='formatted')
2766        open (unit=unitnb+20,file=trim(tmpfilim),status='unknown',form='formatted')
2767        write(unitnb+10,'(400e26.16)') freq,(((((real(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)),&
2768        im=1,2*lpawu+1),ispinor=1,matlu(1)%nspinor),im1=1,2*lpawu+1),ispinor1=1,matlu(1)%nspinor),isppol=1,matlu(1)%nsppol)
2769        write(unitnb+20,'(400e26.16)') freq,(((((aimag(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)),&
2770        im=1,2*lpawu+1),ispinor=1,matlu(1)%nspinor),im1=1,2*lpawu+1),ispinor1=1,matlu(1)%nspinor),isppol=1,matlu(1)%nsppol)
2771      else
2772        tmpfil = trim(char1)//tag_at
2773        open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
2774        write(unitnb,'(400e26.16)') freq,(((((matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1),&
2775        im=1,2*lpawu+1),ispinor=1,matlu(1)%nspinor),im1=1,2*lpawu+1),ispinor1=1,matlu(1)%nspinor),isppol=1,matlu(1)%nsppol)
2776      endif
2777    endif
2778  enddo
2779  end subroutine printplot_matlu

m_matlu/prod_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 prod_matlu

FUNCTION

 Do the matrix product of two matlus

COPYRIGHT

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

INPUTS

  matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity
  matlu2(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity

OUTPUT

  matlu3(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: output quantity

SIDE EFFECTS

NOTES

SOURCE

2342  subroutine prod_matlu(matlu1,matlu2,matlu3,natom)
2343  use defs_basis
2344  use defs_wvltypes
2345  use m_crystal, only : crystal_t
2346  implicit none
2347 
2348 !Arguments ------------------------------------
2349 !scalars
2350  integer, intent(in) :: natom
2351 !arrays
2352  type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
2353  type(matlu_type), intent(inout) :: matlu3(natom)
2354 !Local variables-------------------------------
2355 !scalars
2356  integer :: iatom,im1,im2,im3,ispinor1,ispinor2,ispinor3,isppol
2357  integer :: lpawu
2358 !arrays
2359 !************************************************************************
2360  call zero_matlu(matlu3,natom)
2361  do iatom=1,natom
2362    lpawu=matlu1(iatom)%lpawu
2363    if(lpawu.ne.-1) then
2364      do isppol=1,matlu1(1)%nsppol
2365        do ispinor1=1,matlu1(1)%nspinor
2366          do ispinor2=1,matlu1(1)%nspinor
2367            do ispinor3=1,matlu1(1)%nspinor
2368              do im1=1,2*lpawu+1
2369                do im2=1,2*lpawu+1
2370                  do im3=1,2*lpawu+1
2371                    matlu3(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)= &
2372 &                    matlu3(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)+ &
2373 &                    matlu1(iatom)%mat(im1,im3,isppol,ispinor1,ispinor3)*&
2374 &                    matlu2(iatom)%mat(im3,im2,isppol,ispinor3,ispinor2)
2375                  enddo ! im3
2376                enddo ! im2
2377              enddo ! im1
2378            enddo ! ispinor3
2379          enddo ! ispinor2
2380        enddo ! ispinor1
2381      enddo ! isppol
2382    endif ! lpawu
2383  enddo ! iatom
2384 
2385  end subroutine prod_matlu

m_matlu/rotate_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 rotate_matlu

FUNCTION

 Rotate matlu matrix

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  rot_mat(natom) <type(coeff2c_type)> = Rotation matrix (usually from diag_matlu)
  natom=number of atoms in cell.
  prtopt= option to define level of printing

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: Rotated matrix

SIDE EFFECTS

NOTES

SOURCE

1872  subroutine rotate_matlu(matlu,rot_mat,natom,prtopt,inverse)
1873  use defs_basis
1874  use defs_wvltypes
1875  use m_crystal, only : crystal_t
1876  implicit none
1877 
1878 !Arguments ------------------------------------
1879 !scalars
1880  integer, intent(in) :: natom
1881  integer, intent(in) :: prtopt
1882  integer, intent(in) :: inverse
1883 !arrays
1884  type(matlu_type),intent(inout) :: matlu(natom)
1885  type(coeff2c_type),optional,intent(inout) :: rot_mat(natom,matlu(1)%nsppol)
1886 !Local variables-------------------------------
1887 !scalars
1888  integer :: iatom,im1,im2,isppol
1889  integer :: nsppol,nspinor,tndim
1890 !arrays
1891  type(coeff2c_type),allocatable :: gathermatlu(:)
1892 ! type(coeff2c_type),allocatable :: rot_mat_orig(:,:)
1893  type(coeff2c_type),allocatable :: rot_mat_orig(:)
1894  complex(dpc),allocatable :: temp_mat(:,:)
1895 !************************************************************************
1896  if(prtopt==1) then
1897  endif
1898  nsppol=matlu(1)%nsppol
1899  nspinor=matlu(1)%nspinor
1900 ! ABI_MALLOC(rot_mat_orig,(natom,matlu(1)%nsppol))
1901 
1902  do isppol=1,nsppol
1903 
1904 ! ===========================
1905 ! Define gathermatlu and rot_mat_orig and allocate
1906 ! ===========================
1907    ABI_MALLOC(rot_mat_orig,(natom))
1908    ABI_MALLOC(gathermatlu,(natom))
1909    do iatom=1,natom
1910      if(matlu(iatom)%lpawu.ne.-1) then
1911        tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1912        ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
1913        gathermatlu(iatom)%value=czero
1914 !       ABI_MALLOC(rot_mat_orig(iatom,isppol)%value,(tndim,tndim))
1915 !       rot_mat_orig(iatom,isppol)%value(:,:)=rot_mat(iatom,isppol)%value(:,:)
1916        ABI_MALLOC(rot_mat_orig(iatom)%value,(tndim,tndim))
1917        rot_mat_orig(iatom)%value(:,:)=rot_mat(iatom,isppol)%value(:,:)
1918      endif
1919    enddo
1920    if(nsppol==1.and.nspinor==2) then
1921      call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
1922    else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
1923      do iatom=1,natom
1924        if(matlu(iatom)%lpawu.ne.-1) then
1925          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1926          do im1=1,tndim
1927            do im2=1,tndim
1928              gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,isppol,1,1)
1929            enddo
1930          enddo
1931        endif
1932      enddo
1933    endif
1934         ! write(std_out,*) "gathermatlu in rotate matlu"
1935         ! do im1=1,tndim
1936         !   write(message,'(12(1x,18(1x,"(",e17.10,",",e17.10,")")))')&
1937         !    (gathermatlu(1)%value(im1,im2),im2=1,tndim)
1938         !   call wrtout(std_out,message,'COLL')
1939         ! end do
1940 
1941 ! ===========================
1942 ! If necessary, invert rot_mat
1943 ! ===========================
1944    if(inverse==1) then
1945      do iatom=1,natom
1946        if(matlu(iatom)%lpawu.ne.-1) then
1947          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1948            do im1=1,tndim
1949              do im2=1,tndim
1950 !               rot_mat(iatom,isppol)%value(im1,im2)=conjg(rot_mat_orig(iatom,isppol)%value(im2,im1))
1951                rot_mat(iatom,isppol)%value(im1,im2)=conjg(rot_mat_orig(iatom)%value(im2,im1))
1952              enddo
1953            enddo
1954        endif ! lpawu
1955      enddo ! iatom
1956    endif
1957         ! write(std_out,*) "rot_mat_orig "
1958         ! do im1=1,tndim
1959         !   write(message,'(12(1x,18(1x,"(",e18.10,",",e18.10,")")))')&
1960         ! &   (rot_mat_orig(1)%value(im1,im2),im2=1,tndim)
1961         !   call wrtout(std_out,message,'COLL')
1962         ! end do
1963         ! write(std_out,*) "rot_mat "
1964         ! do im1=1,tndim
1965         !   write(message,'(12(1x,18(1x,"(",e18.10,",",e18.10,")")))')&
1966         ! &   (rot_mat(1,1)%value(im1,im2),im2=1,tndim)
1967         !   call wrtout(std_out,message,'COLL')
1968         ! end do
1969 
1970 ! ===========================
1971 ! Rotate
1972 ! ===========================
1973    ABI_MALLOC(temp_mat,(tndim,tndim))
1974    do iatom=1,natom
1975      if(matlu(iatom)%lpawu.ne.-1) then
1976        tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1977        temp_mat(:,:)=czero
1978 !      input matrix: gathermatlu
1979 !      rotation matrix: rot_mat
1980 !      intermediate matrix: temp_mat
1981 !      result matrix: gathermatlu
1982        ! temp_mat = gathermatlu * conjg(rot_mat)
1983        call zgemm('n','c',tndim,tndim,tndim,cone,gathermatlu(iatom)%value   ,tndim,&
1984 &        rot_mat(iatom,isppol)%value,tndim,czero,temp_mat                ,tndim)
1985        ! gathermatlu = rot_mat * temp_mat = rot_mat * gathermatlu * conjg(rot_mat)
1986        call zgemm('n','n',tndim,tndim,tndim,cone,rot_mat(iatom,isppol)%value,tndim,&
1987 &        temp_mat                   ,tndim,czero,gathermatlu(iatom)%value,tndim)
1988      endif ! lpawu
1989    enddo ! iatom
1990    !do iatom=1,natom
1991    !  if(matlu(iatom)%lpawu.ne.-1) then
1992    !    write(std_out,*) "temp_mat in rotate_matlu 2"
1993    !    do im1=1,tndim
1994    !      write(message,'(12(1x,18(1x,"(",f17.10,",",f17.10,")")))')&
1995    !&       (temp_mat(im1,im2),im2=1,tndim)
1996    !      call wrtout(std_out,message,'COLL')
1997    !    end do
1998    !  endif
1999    !enddo
2000    !do iatom=1,natom
2001    !  if(matlu(iatom)%lpawu.ne.-1) then
2002    !    write(std_out,*) "gathermatlu in rotate_matlu 2"
2003    !    do im1=1,tndim
2004    !      write(message,'(12(1x,18(1x,"(",f17.10,",",f17.10,")")))')&
2005    !&       (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
2006    !      call wrtout(std_out,message,'COLL')
2007    !    end do
2008    !  endif
2009    !enddo
2010    ABI_FREE(temp_mat)
2011      !ABI_ERROR("Aborting now")
2012 
2013 ! Choose inverse rotation: reconstruct correct rot_mat from rot_mat_orig
2014 ! ========================================================================
2015    if(inverse==1) then
2016      do iatom=1,natom
2017        if(matlu(iatom)%lpawu.ne.-1) then
2018          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2019            do im1=1,tndim
2020              do im2=1,tndim
2021 !               rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom,isppol)%value(im1,im2)
2022                rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom)%value(im1,im2)
2023              enddo
2024            enddo
2025        endif ! lpawu
2026      enddo ! iatom
2027    endif
2028 
2029 ! ===========================
2030 ! Put data into matlu(iatom)
2031 ! ===========================
2032    if(nsppol==1.and.nspinor==2) then
2033      call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
2034    else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
2035      do iatom=1,natom
2036        if(matlu(iatom)%lpawu.ne.-1) then
2037          tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2038          do im1=1,tndim
2039            do im2=1,tndim
2040              matlu(iatom)%mat(im1,im2,isppol,1,1)= gathermatlu(iatom)%value(im1,im2)
2041            enddo
2042          enddo
2043        endif
2044      enddo
2045    endif ! test nsppol/nspinor
2046 ! ===========================
2047 ! Deallocations
2048 ! ===========================
2049    do iatom=1,natom
2050      if(matlu(iatom)%lpawu.ne.-1) then
2051        ABI_FREE(gathermatlu(iatom)%value)
2052 !       ABI_FREE(rot_mat_orig(iatom,isppol)%value)
2053        ABI_FREE(rot_mat_orig(iatom)%value)
2054      endif
2055    enddo
2056    ABI_FREE(gathermatlu)
2057    ABI_FREE(rot_mat_orig)
2058  enddo ! isppol
2059 
2060  end subroutine rotate_matlu

m_matlu/shift_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 shift_matlu

FUNCTION

 shift matlu matrix

COPYRIGHT

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

INPUTS

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity to rotate
  natom=number of atoms in cell.
  shift= shift of the diagonal part.

OUTPUT

  matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: shifted matrix

SIDE EFFECTS

NOTES

SOURCE

2089  subroutine shift_matlu(matlu,natom,shift,signe)
2090  use defs_basis
2091  use defs_wvltypes
2092  use m_crystal, only : crystal_t
2093  implicit none
2094 
2095 !Arguments ------------------------------------
2096 !scalars
2097  integer, intent(in) :: natom
2098 !arrays
2099  type(matlu_type),intent(inout) :: matlu(natom)
2100  complex(dpc),intent(in) :: shift(natom)
2101  integer, optional,intent(in) :: signe
2102 !Local variables-------------------------------
2103 !scalars
2104  integer :: iatom,im,ispinor,isppol
2105  integer :: lpawu,signe_used
2106 ! character(len=500) :: message
2107 !arrays
2108 !************************************************************************
2109  signe_used=1
2110  if(present(signe)) then
2111    if(signe==-1) signe_used=-1
2112  endif
2113  do iatom=1,natom
2114    lpawu=matlu(iatom)%lpawu
2115    if(lpawu.ne.-1) then
2116      do im=1,2*lpawu+1
2117        do isppol=1,matlu(1)%nsppol
2118          do ispinor=1,matlu(1)%nspinor
2119            matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=&
2120 &           matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)-signe_used*shift(iatom)
2121          enddo
2122        enddo
2123      enddo
2124    endif
2125  enddo
2126 
2127  end subroutine shift_matlu

m_matlu/slm2ylm_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 slm2ylm_matlu

FUNCTION

 Transform mat from Slm to Ylm basis or vice versa

COPYRIGHT

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

INPUTS

  matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity
  natom :: number of atoms
  option=1 go from Slm to Ylm basis
  option=2 go from Ylm to Slm basis

SIDE EFFECTS

NOTES

SOURCE

2530  subroutine slm2ylm_matlu(matlu,natom,option,optprt)
2531  use defs_basis
2532  use defs_wvltypes
2533  implicit none
2534 
2535 !Arguments ------------------------------------
2536 !scalars
2537  integer, intent(in) :: natom,option,optprt
2538 !arrays
2539  type(matlu_type), intent(inout) :: matlu(natom)
2540 !Local variables-------------------------------
2541 !scalars
2542  integer :: iatom,im,ispinor,isppol,ispinor2
2543  integer :: lpawu,ll,mm,jm,ii,jj,im1,im2
2544  character(len=500) :: message
2545  real(dp) :: onem
2546  complex(dpc),allocatable :: slm2ylm(:,:)
2547  complex(dpc),allocatable :: mat_inp_c(:,:)
2548  complex(dpc),allocatable :: mat_out_c(:,:)
2549  complex(dpc) :: tmp2
2550  real(dp),parameter :: invsqrt2=one/sqrt2
2551 !arrays
2552 !************************************************************************
2553 
2554  do iatom=1,natom
2555    lpawu=matlu(iatom)%lpawu
2556    if(lpawu.ne.-1) then
2557      ll=lpawu
2558      ABI_MALLOC(slm2ylm,(2*ll+1,2*ll+1))
2559      slm2ylm=czero
2560      do im=1,2*ll+1
2561        mm=im-ll-1;jm=-mm+ll+1
2562        onem=dble((-1)**mm)
2563        if (mm> 0) then
2564          slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
2565          slm2ylm(jm,im)= cmplx(invsqrt2,     zero,kind=dp)
2566        end if
2567        if (mm==0) then
2568          slm2ylm(im,im)=cone
2569        end if
2570        if (mm< 0) then
2571          slm2ylm(im,im)= cmplx(zero,     invsqrt2,kind=dp)
2572          slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
2573        end if
2574      end do
2575      if(optprt>2) then
2576        write(message,'(2a)') ch10,"SLM2YLM matrix"
2577        call wrtout(std_out,message,'COLL')
2578        do im1=1,ll*2+1
2579          write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2580 &         (slm2ylm(im1,im2),im2=1,ll*2+1)
2581          call wrtout(std_out,message,'COLL')
2582        end do
2583      endif
2584      do isppol=1,matlu(1)%nsppol
2585        do ispinor=1,matlu(1)%nspinor
2586          do ispinor2=1,matlu(1)%nspinor
2587            ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
2588            ABI_MALLOC(mat_inp_c,(2*ll+1,2*ll+1))
2589            mat_inp_c(:,:) = matlu(iatom)%mat(:,:,isppol,ispinor,ispinor2)
2590            mat_out_c=czero
2591 
2592            if(optprt>2) then
2593              write(message,'(2a, i2, a, i2, a, i2)') ch10,"SLM input matrix, isppol=", isppol, ", ispinor=", ispinor,& 
2594 &             ", ispinor2=", ispinor2
2595              call wrtout(std_out,message,'COLL')
2596              do im1=1,ll*2+1
2597                write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2598 &               (mat_inp_c(im1,im2),im2=1,ll*2+1)
2599                call wrtout(std_out,message,'COLL')
2600              end do
2601            endif
2602 
2603            do jm=1,2*ll+1
2604              do im=1,2*ll+1
2605                tmp2=czero
2606                do ii=1,2*ll+1
2607                  do jj=1,2*ll+1
2608                    if(option==1) then
2609                      tmp2=tmp2+mat_inp_c(ii,jj)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))
2610                    else if(option==2) then
2611                      tmp2=tmp2+mat_inp_c(ii,jj)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))
2612                    end if
2613                  end do
2614                end do
2615                mat_out_c(im,jm)=tmp2
2616              end do
2617            end do
2618 
2619            if(optprt>2) then
2620              write(message,'(2a, i2, a, i2, a, i2)') ch10,"YLM output matrix, isppol=", isppol, ", ispinor=", ispinor,&
2621 &             ", ispinor2=", ispinor2
2622              call wrtout(std_out,message,'COLL')
2623              do im1=1,ll*2+1
2624                write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2625       &         (mat_out_c(im1,im2),im2=1,ll*2+1)
2626                call wrtout(std_out,message,'COLL')
2627              end do
2628            endif
2629 
2630            matlu(iatom)%mat(:,:,isppol,ispinor,ispinor2)=mat_out_c(:,:)
2631            ABI_FREE(mat_out_c)
2632            ABI_FREE(mat_inp_c)
2633          enddo ! im
2634        enddo ! ispinor
2635      enddo ! isppol
2636      ABI_FREE(slm2ylm)
2637    endif ! lpawu
2638  enddo ! iatom
2639 
2640 
2641  end subroutine slm2ylm_matlu

m_matlu/sym_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 sym_matlu

FUNCTION

 Symetrise local quantity.

COPYRIGHT

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

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  gloc(natom) <type(matlu_type)>= density matrix in the local orbital basis and related variables
  pawang <type(pawang)>=paw angular mesh and related data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

  gloc(natom) <type(matlu_type)>= density matrix symetrized in the local orbital basis and related variables

SIDE EFFECTS

NOTES

SOURCE

539  subroutine sym_matlu(cryst_struc,gloc,pawang,paw_dmft)
540 
541  use defs_basis
542 ! use defs_wvltypes
543  use m_pawang, only : pawang_type
544  use m_crystal, only : crystal_t
545  use m_paw_dmft, only: paw_dmft_type
546 
547  implicit none
548 
549 !Arguments ------------------------------------
550 !scalars
551  type(crystal_t),intent(in) :: cryst_struc
552  type(pawang_type),intent(in) :: pawang
553  type(paw_dmft_type), intent(in) :: paw_dmft
554 !arrays
555  type(matlu_type),intent(inout) :: gloc(cryst_struc%natom)
556 !scalars
557 !Local variables-------------------------------
558 !scalars
559  integer :: at_indx,iatom,irot,ispinor,ispinor1,isppol,lpawu,m1,m2,m3,m4,mu
560  integer :: natom,ndim,nsppol,nspinor,nu,t2g,m1s,m2s,m3s,m4s,lpawu_zarot,x2my2d
561  complex(dpc) :: sumrho,summag(3),rotmag(3),ci
562  real(dp) :: zarot2
563 !arrays
564 ! complex(dpc),allocatable :: glocnm(:,:,:,:,:)
565  type(matlu_type),allocatable :: glocnm(:)
566 ! complex(dpc),allocatable :: glocnms(:,:,:,:,:)
567  type(matlu_type),allocatable :: glocnms(:)
568  type(matlu_type),allocatable :: glocsym(:)
569  real(dp),allocatable :: symrec_cart(:,:,:)
570  integer :: mt2g(3),mx2my2d
571  mt2g(1)=1
572  mt2g(2)=2
573  mt2g(3)=4
574  mx2my2d=5
575  t2g=paw_dmft%dmftqmc_t2g
576  x2my2d=paw_dmft%dmftqmc_x2my2d
577 
578 ! DBG_ENTER("COLL")
579 
580  ci=cone
581  nspinor=gloc(1)%nspinor
582  nsppol=gloc(1)%nsppol
583  natom=cryst_struc%natom
584 
585  ABI_MALLOC(glocnm,(natom))
586  ABI_MALLOC(glocnms,(natom))
587  ABI_MALLOC(glocsym,(natom))
588  call init_matlu(natom,nspinor,nsppol,gloc(1:natom)%lpawu,glocnm)
589  call init_matlu(natom,nspinor,nsppol,gloc(1:natom)%lpawu,glocnms)
590  call init_matlu(natom,nspinor,nsppol,gloc(1:natom)%lpawu,glocsym)
591 
592 
593 !=========  Case nspinor ==1 ========================
594 
595  if (nspinor==1) then
596   ispinor=1
597   ispinor1=1
598   do iatom=1,cryst_struc%natom
599    do isppol=1,nsppol
600     if(gloc(iatom)%lpawu/=-1) then
601      lpawu=gloc(iatom)%lpawu
602      do m1=1, 2*lpawu+1
603       do m2=1, 2*lpawu+1
604        do irot=1,cryst_struc%nsym
605         at_indx=cryst_struc%indsym(4,irot,iatom)
606         do m3=1, 2*lpawu+1
607          do m4=1, 2*lpawu+1
608           if(t2g==1) then
609            m1s=mt2g(m1)
610            m2s=mt2g(m2)
611            m3s=mt2g(m3)
612            m4s=mt2g(m4)
613            lpawu_zarot=2
614           else if (x2my2d==1) then
615            m1s=mx2my2d
616            m2s=mx2my2d
617            m3s=mx2my2d
618            m4s=mx2my2d
619            lpawu_zarot=2
620           else
621            m1s=m1
622            m2s=m2
623            m3s=m3
624            m4s=m4
625            lpawu_zarot=lpawu
626           endif
627           zarot2=pawang%zarot(m3s,m1s,lpawu_zarot+1,irot)*pawang%zarot(m4s,m2s,lpawu_zarot+1,irot)
628           glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)=&
629 &          glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)&
630 &          +gloc(at_indx)%mat(m3,m4,isppol,ispinor,ispinor1)*zarot2
631          end do  ! m3
632         end do  ! m4
633        end do  ! irot
634        glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)=&
635 &       glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)/real(cryst_struc%nsym,kind=dp)
636       end do ! m2
637      end do ! m1
638     endif ! lpawu/=-1
639    end do ! isppol
640   end do ! iatom
641 !==  Put glocsym into gloc
642   do iatom=1,cryst_struc%natom
643     if(gloc(iatom)%lpawu/=-1) then
644       gloc(iatom)%mat=glocsym(iatom)%mat
645 !      gloc(iatom)%mat(:,:,1,:,:)=(glocsym(iatom)%mat(:,:,1,:,:) &
646 !&      + glocsym(iatom)%mat(:,:,2,:,:))/two
647 !      gloc(iatom)%mat(:,:,2,:,:)= gloc(iatom)%mat(:,:,1,:,:)
648 !      write(std_out,*) "WARNING: SYM non mag"
649 !      write(ab_out,*) "WARNING: SYM non mag"
650     endif
651   end do ! iatom
652 
653 !=========  Case nspinor ==2 ========================
654 
655  else if (nspinor==2) then
656 
657 !== Allocate temporary arrays
658   do iatom=1,cryst_struc%natom
659    if(gloc(iatom)%lpawu/=-1) then
660     ndim=2*gloc(iatom)%lpawu+1
661     ABI_FREE(glocnm(iatom)%mat)
662     ABI_FREE(glocnms(iatom)%mat)
663     ABI_FREE(glocsym(iatom)%mat)
664     ABI_MALLOC(glocnm(iatom)%mat,(ndim,ndim,nsppol,4,1))
665     ABI_MALLOC(glocnms(iatom)%mat,(ndim,ndim,nsppol,4,1))
666     ABI_MALLOC(glocsym(iatom)%mat,(ndim,ndim,nsppol,2,2))
667    endif
668   enddo
669   ABI_MALLOC(symrec_cart,(3,3,cryst_struc%nsym))
670 
671 !==  Compute symrec_cart
672   do irot=1,cryst_struc%nsym
673    call symredcart(cryst_struc%gprimd,cryst_struc%rprimd,symrec_cart(:,:,irot),cryst_struc%symrec(:,:,irot))
674   end do
675 
676 !==  Compute density matrix in density and magnetization representation
677   call chg_repr_matlu(gloc,glocnm,cryst_struc%natom,option=1,prtopt=1)
678 
679 !==  Do the sum over symetrized density matrix (in n,m repr)
680   isppol=1
681   do iatom=1,cryst_struc%natom
682    if(gloc(iatom)%lpawu/=-1) then
683     lpawu=gloc(iatom)%lpawu
684     ndim=2*gloc(iatom)%lpawu+1
685     do m1=1, 2*lpawu+1
686      do m2=1, 2*lpawu+1
687       sumrho=czero
688       rotmag=czero
689       do irot=1,cryst_struc%nsym
690        summag=czero
691        at_indx=cryst_struc%indsym(4,irot,iatom)
692        do m3=1, 2*lpawu+1
693         do m4=1, 2*lpawu+1
694           if(t2g==1) then
695            m1s=mt2g(m1)
696            m2s=mt2g(m2)
697            m3s=mt2g(m3)
698            m4s=mt2g(m4)
699            lpawu_zarot=2
700           else if (x2my2d==1) then
701            m1s=mx2my2d
702            m2s=mx2my2d
703            m3s=mx2my2d
704            m4s=mx2my2d
705            lpawu_zarot=2
706           else
707            m1s=m1
708            m2s=m2
709            m3s=m3
710            m4s=m4
711            lpawu_zarot=lpawu
712           endif
713          zarot2=pawang%zarot(m3s,m2s,lpawu_zarot+1,irot)*pawang%zarot(m4s,m1s,lpawu_zarot+1,irot)
714          sumrho=sumrho +  glocnm(at_indx)%mat(m4,m3,isppol,1,1)  * zarot2
715          do mu=1,3
716           summag(mu)=summag(mu) + glocnm(at_indx)%mat(m4,m3,isppol,mu+1,1) * zarot2
717          enddo
718         end do ! m3
719        end do !m4
720 
721 !       ==  special case of magnetization
722        do nu=1,3
723         do mu=1,3
724          rotmag(mu)=rotmag(mu)+symrec_cart(mu,nu,irot)*summag(nu)
725         end do
726        end do
727 !      write(std_out,'(a,3i4,2x,3(2f10.5,2x))') "rotmag",irot,m1,m2,(rotmag(mu),mu=1,3)
728       end do ! irot
729 
730 !       ==  Normalizes sum
731       sumrho=sumrho/real(cryst_struc%nsym,kind=dp)
732 !        sumrho=glocnm(isppol,1,iatom,m1,m2) ! test without sym
733       glocnms(iatom)%mat(m1,m2,isppol,1,1)=sumrho
734       do mu=1,3
735        rotmag(mu)=rotmag(mu)/real(cryst_struc%nsym,kind=dp)
736 !          rotmag(mu)=glocnm(isppol,mu+1,iatom,m1,m2) ! test without sym
737        glocnms(iatom)%mat(m1,m2,isppol,mu+1,1)=rotmag(mu)
738       enddo
739      end do  ! m2
740     end do ! m1
741    endif ! lpawu/=-1
742   end do ! iatom
743 
744 !==  Compute back density matrix in upup dndn updn dnup representation
745   call chg_repr_matlu(glocsym,glocnms,cryst_struc%natom,option=-1,prtopt=1)
746 
747 !==  Put glocsym into gloc
748   do iatom=1,cryst_struc%natom
749     if(gloc(iatom)%lpawu/=-1) then
750       gloc(iatom)%mat=glocsym(iatom)%mat
751 !      gloc(iatom)%mat(:,:,1,:,:)=(glocsym(iatom)%mat(:,:,1,:,:) &
752 !&      + glocsym(iatom)%mat(:,:,2,:,:))/two
753 !      gloc(iatom)%mat(:,:,2,:,:)= gloc(iatom)%mat(:,:,1,:,:)
754 !      write(std_out,*) "WARNING: SYM non mag"
755 !      write(ab_out,*) "WARNING: SYM non mag"
756     endif
757   end do ! iatom
758 
759   ABI_FREE(symrec_cart)
760  endif
761 
762  call destroy_matlu(glocnm,cryst_struc%natom)
763  call destroy_matlu(glocnms,cryst_struc%natom)
764  call destroy_matlu(glocsym,cryst_struc%natom)
765  ABI_FREE(glocnm)
766  ABI_FREE(glocnms)
767  ABI_FREE(glocsym)
768 !==============end of nspinor==2 case ===========
769 
770 
771 ! DBG_EXIT("COLL")
772 
773  end subroutine sym_matlu

m_matlu/trace_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 trace_matlu

FUNCTION

  Compute the trace of the matlu matrix

INPUTS

  maltu(natom) <type(matlu_type)>= density matrix in the
               local orbital basis and related variables
  natom = number of atoms

OUTPUT

  trace_loc(natom,nsppol+1)= trace for each atoms and each polarization,
                             trace_loc(iatom,nsppol+1) is
                             the full trace over polarization also.

SOURCE

1154  subroutine trace_matlu(matlu,natom,trace_loc,itau)
1155 
1156  use defs_basis
1157  implicit none
1158 
1159 !Arguments ------------------------------------
1160 !type
1161  integer, intent(in) :: natom
1162  type(matlu_type), intent(in) :: matlu(natom)
1163  real(dp),intent(inout),target,optional :: trace_loc(natom,matlu(1)%nsppol+1)
1164  integer, intent(in),optional :: itau
1165 
1166 !local variables-------------------------------
1167  integer :: iatom,isppol,ispinor,m,lpawu
1168  integer :: nsppol,nspinor
1169  real(dp), pointer :: traceloc(:,:)=>null()
1170  character(len=500) :: message
1171 ! *********************************************************************
1172  nsppol=matlu(1)%nsppol
1173  nspinor=matlu(1)%nspinor
1174  if(present(trace_loc)) then
1175    traceloc=>trace_loc
1176  else
1177    ABI_MALLOC(traceloc,(natom,matlu(1)%nsppol+1))
1178  endif
1179 
1180  traceloc=zero
1181  do iatom = 1 , natom
1182    lpawu=matlu(iatom)%lpawu
1183    if(lpawu/=-1) then
1184       write(message,'(2a,i4)')  ch10,'   -------> For Correlated Atom', iatom
1185      if(.not.present(itau)) then
1186        call wrtout(std_out,  message,'COLL')
1187      end if
1188      if(present(itau)) then
1189        if (itau>0) then
1190          call wrtout(std_out,  message,'COLL')
1191        end if
1192      endif
1193      do isppol = 1 , nsppol
1194        do ispinor = 1 , nspinor
1195          do m = 1 ,  2*lpawu+1
1196            traceloc(iatom,isppol)=traceloc(iatom,isppol)+&
1197 &           matlu(iatom)%mat(m,m,isppol,ispinor,ispinor)
1198          enddo
1199        enddo
1200       traceloc(iatom,nsppol+1)=traceloc(iatom,nsppol+1)+traceloc(iatom,isppol)
1201      enddo
1202      if(nsppol==1.and.nspinor==1)  traceloc(iatom,nsppol+1)=traceloc(iatom,nsppol+1)*two
1203      if(.not.present(itau)) then
1204        write(message,'(8x,a,f12.6)')   'Nb of Corr. elec. from G(w) is:'&
1205 &       ,traceloc(iatom,nsppol+1)
1206        call wrtout(std_out,  message,'COLL')
1207      endif
1208      if(present(itau)) then
1209        if(itau==1) then
1210          write(message,'(8x,a,f12.6)')   'Nb of Corr. elec. from G(tau) is:'&
1211 &         ,traceloc(iatom,nsppol+1)
1212          call wrtout(std_out,  message,'COLL')
1213        else if(itau==-1) then
1214          write(message,'(8x,a,f12.6)')   'Nb: Sum of the values of G0(tau=0-) is:'&
1215 &       ,traceloc(iatom,nsppol+1)
1216          call wrtout(std_out,  message,'COLL')
1217        else if(itau==4) then
1218          write(message,'(8x,a,f12.6)')   'Trace of matlu matrix is:'&
1219 &         ,traceloc(iatom,nsppol+1)
1220          call wrtout(std_out,  message,'COLL')
1221        endif
1222      endif
1223    endif
1224  enddo
1225  do iatom = 1 , natom
1226    lpawu=matlu(iatom)%lpawu
1227    if(lpawu/=-1) then
1228      !! MAG
1229      if(nsppol>1) then
1230     ! if(nsppol>1.and.present(itau)) then
1231     !   if(itau==1) then
1232          write(message,'(8x,a,f12.6)')   'DMFT Cor. El. Mag: ',traceloc(iatom,2)-traceloc(iatom,1)
1233          call wrtout(std_out,  message,'COLL')
1234     !   endif
1235      endif
1236 
1237    endif
1238  enddo
1239  if(.not.present(trace_loc)) then
1240   ABI_FREE(traceloc)
1241   traceloc => null()
1242  endif
1243 
1244  end subroutine trace_matlu

m_matlu/zero_matlu [ Functions ]

[ Top ] [ m_matlu ] [ Functions ]

NAME

 zero_matlu

FUNCTION

  zero_matlu

INPUTS

  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables
  natom = number of atoms

OUTPUT

  maltu <type(matlu_type)>= density matrix in the local orbital basis and related variables

SOURCE

194 subroutine zero_matlu(matlu,natom,onlynondiag)
195 
196  use defs_basis
197  implicit none
198 
199 !Arguments ------------------------------------
200 !scalars
201  integer, intent(in) :: natom
202  type(matlu_type),intent(inout) :: matlu(natom)
203  integer, optional, intent(in) :: onlynondiag
204 !Local variables-------------------------------
205  integer :: iatom,im,im1,ispinor,ispinor1,isppol,ndim
206 
207 !*********************************************************************
208 
209  if(present(onlynondiag)) then
210    do iatom=1,natom
211      if(matlu(iatom)%lpawu.ne.-1) then
212        do ispinor=1,matlu(iatom)%nspinor
213          ndim=(2*matlu(iatom)%lpawu+1)
214          do im=1,ndim
215            do im1=1,ndim
216              do ispinor1=1,matlu(iatom)%nspinor
217                if(im/=im1.or.ispinor/=ispinor1) then
218                  do isppol=1,matlu(iatom)%nsppol
219                    matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=czero
220                  enddo
221                end if
222              end do
223            end do
224          end do
225        end do
226      endif
227    enddo
228  else
229    do iatom=1,natom
230     matlu(iatom)%mat=czero
231    enddo
232  endif
233 
234 
235 end subroutine zero_matlu