src/symdmatlu.f90

   1
   2 ! Copyright (C) 2007 F. Bultmark, F. Cricchio, L. Nordstrom and J. K. Dewhurst.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 subroutine symdmatlu(dmat)
   7 use modmain
   8 implicit none
   9 ! arguments
  10 complex(8), intent(inout) :: dmat(lmmaxlu,lmmaxlu,nspinor,nspinor,natmtot)
  11 ! local variables
  12 integer isym,lspl,lspn
  13 integer ispn,jspn
  14 integer is,ia,ja,ias,jas
  15 integer lm1,lm2
  16 real(8) det,t1
  17 ! automatic arrays
  18 logical done(natmmax)
  19 ! allocatable arrays
  20 complex(8), allocatable :: zflm(:,:)
  21 complex(8), allocatable :: ulm(:,:,:)
  22 complex(8), allocatable :: su2(:,:,:)
  23 complex(8), allocatable :: dm1(:,:,:,:,:)
  24 complex(8), allocatable :: dm2(:,:)
  25 complex(8), allocatable :: dm3(:,:,:,:)
  26 complex(8), allocatable :: dm4(:,:)
  27 complex(8), allocatable :: dm5(:,:)
  28 if (ldapu.eq.0) return
  29 ! allocate local arrays
  30 allocate(zflm(lmmaxlu,lmmaxlu))
  31 allocate(ulm(lmmaxlu,lmmaxlu,nsymlat))
  32 allocate(su2(nspinor,nspinor,nsymlat))
  33 allocate(dm1(lmmaxlu,lmmaxlu,nspinor,nspinor,natmmax))
  34 allocate(dm2(lmmaxlu,lmmaxlu))
  35 allocate(dm3(lmmaxlu,lmmaxlu,nspinor,nspinor))
  36 allocate(dm4(nspinor,nspinor),dm5(nspinor,nspinor))
  37 ! setup a complex unit matrix for (l,m) components
  38 zflm(:,:)=0.d0
  39 do lm1=1,lmmaxlu
  40   zflm(lm1,lm1)=1.d0
  41 end do
  42 do isym=1,nsymlat
  43 ! construct (l,m) rotation matrix for each lattice symmetry
  44   call rotzflm(symlatc(1,1,isym),lmaxlu,lmmaxlu,lmmaxlu,zflm,ulm(1,1,isym))
  45 ! construct SU(2) matrix for proper rotation of spinor components
  46 ! (note that rotsu2 uses only the proper part of the rotation matrix)
  47   if (spinpol) then
  48     call rotsu2(symlatc(1,1,isym),det,su2(1,1,isym))
  49   end if
  50 end do
  51 t1=1.d0/dble(nsymcrys)
  52 do is=1,nspecies
  53 ! make copy of the density matrices
  54   do ia=1,natoms(is)
  55     ias=idxas(ia,is)
  56     dm1(:,:,:,:,ia)=dmat(:,:,:,:,ias)
  57   end do
  58   done(:)=.false.
  59   do ia=1,natoms(is)
  60     if (.not.done(ia)) then
  61       ias=idxas(ia,is)
  62       dmat(:,:,:,:,ias)=0.d0
  63       do isym=1,nsymcrys
  64         lspl=lsplsymc(isym)
  65         lspn=lspnsymc(isym)
  66 ! equivalent atom index (symmetry rotates atom ja into atom ia)
  67         ja=ieqatom(ia,is,isym)
  68         jas=idxas(ja,is)
  69 ! apply (l,m) symmetry matrix as U*D*conjg(U')
  70         do ispn=1,nspinor
  71           do jspn=1,nspinor
  72             call zgemm('N','N',lmmaxlu,lmmaxlu,lmmaxlu,zone,ulm(1,1,lspl), &
  73              lmmaxlu,dm1(1,1,ispn,jspn,ja),lmmaxlu,zzero,dm2,lmmaxlu)
  74             call zgemm('N','C',lmmaxlu,lmmaxlu,lmmaxlu,zone,dm2,lmmaxlu, &
  75              ulm(1,1,lspl),lmmaxlu,zzero,dm3(1,1,ispn,jspn),lmmaxlu)
  76           end do
  77         end do
  78 ! apply SU(2) symmetry matrix as U*D*conjg(U') and add
  79         if (spinpol) then
  80           do lm1=1,lmmaxlu
  81             do lm2=1,lmmaxlu
  82               dm4(:,:)=dm3(lm1,lm2,:,:)
  83               call z2mm(su2(1,1,lspn),dm4,dm5)
  84               call z2mmct(dm5,su2(1,1,lspn),dm4)
  85               dmat(lm1,lm2,:,:,ias)=dmat(lm1,lm2,:,:,ias)+dm4(:,:)
  86             end do
  87           end do
  88         else
  89           dmat(:,:,1,1,ias)=dmat(:,:,1,1,ias)+dm3(:,:,1,1)
  90         end if
  91 ! end loop over crystal symmetries
  92       end do
  93 ! normalise
  94       dmat(:,:,:,:,ias)=t1*dmat(:,:,:,:,ias)
  95       done(ia)=.true.
  96 ! rotate into equivalent atoms
  97       do isym=1,nsymcrys
  98         ja=ieqatom(ia,is,isym)
  99         if (.not.done(ja)) then
 100           jas=idxas(ja,is)
 101           lspl=lsplsymc(isym)
 102           lspn=lspnsymc(isym)
 103 ! apply (l,m) symmetry matrix as conjg(U')*D*U (rotates atom ia into atom ja)
 104           do ispn=1,nspinor
 105             do jspn=1,nspinor
 106               call zgemm('C','N',lmmaxlu,lmmaxlu,lmmaxlu,zone,ulm(1,1,lspl), &
 107                lmmaxlu,dmat(1,1,ispn,jspn,ias),lmmaxlu,zzero,dm2,lmmaxlu)
 108               call zgemm('N','N',lmmaxlu,lmmaxlu,lmmaxlu,zone,dm2,lmmaxlu, &
 109                ulm(1,1,lspl),lmmaxlu,zzero,dm3(1,1,ispn,jspn),lmmaxlu)
 110             end do
 111           end do
 112 ! apply SU(2) symmetry matrix as conjg(U')*D*U
 113           if (spinpol) then
 114             do lm1=1,lmmaxlu
 115               do lm2=1,lmmaxlu
 116                 dm4(:,:)=dm3(lm1,lm2,:,:)
 117                 call z2mctm(su2(1,1,lspn),dm4,dm5)
 118                 call z2mm(dm5,su2(1,1,lspn),dm4)
 119                 dmat(lm1,lm2,:,:,jas)=dm4(:,:)
 120               end do
 121             end do
 122           else
 123             dmat(:,:,1,1,jas)=dm3(:,:,1,1)
 124           end if
 125           done(ja)=.true.
 126         end if
 127       end do
 128     end if
 129   end do
 130 end do
 131 deallocate(zflm,ulm,su2)
 132 deallocate(dm1,dm2,dm3,dm4,dm5)
 133 return
 134 end subroutine
 135