alineso.F90 12.1 KB
Newer Older
1 2 3 4 5 6 7 8
MODULE m_alineso
  USE m_juDFT
  !---------------------------------------------------------------------- 
  ! Set up SO-hamiltonian for 2nd variation (hsohelp and hsoham) and 
  ! diagonalize by lapack routines.
  ! Eigenvalues and vectors (eig_so and zso) are returned 
  !----------------------------------------------------------------------
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
9 10
  SUBROUTINE alineso(eig_id,lapw,&
       mpi,DIMENSION,atoms,sym,kpts,&
11 12 13 14
       input,noco,cell,oneD,&
       rsopp,rsoppd,rsopdp,rsopdpd,nk,&
       rsoplop,rsoplopd,rsopdplo,rsopplo,rsoploplop,&
       usdus,soangl,&
15
       nsize,nmat,&
16 17 18 19 20 21 22 23 24
       eig_so,zso)

#include"cpp_double.h"
    USE m_hsohelp
    USE m_hsoham
    USE m_eig66_io, ONLY : read_eig
    USE m_types
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(IN)         :: mpi
Daniel Wortmann's avatar
Daniel Wortmann committed
25
    TYPE(t_lapw),INTENT(IN)        :: lapw
26 27 28 29 30 31 32
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_oneD),INTENT(IN)        :: oneD
    TYPE(t_input),INTENT(IN)       :: input
    TYPE(t_noco),INTENT(IN)        :: noco
    TYPE(t_sym),INTENT(IN)         :: sym
    TYPE(t_cell),INTENT(IN)        :: cell
    TYPE(t_atoms),INTENT(IN)       :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
33
    TYPE(t_kpts),INTENT(IN)        :: kpts
34 35 36 37 38
    TYPE(t_usdus),INTENT(IN)       :: usdus
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: eig_id
    INTEGER, INTENT (IN) :: nk 
Daniel Wortmann's avatar
Daniel Wortmann committed
39
    INTEGER, INTENT (OUT):: nsize,nmat
40 41
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
42 43 44 45 46 47 48 49 50
    REAL,    INTENT (IN) :: rsopp  (atoms%ntype,atoms%lmaxd,2,2)
    REAL,    INTENT (IN) :: rsoppd (atoms%ntype,atoms%lmaxd,2,2)
    REAL,    INTENT (IN) :: rsopdp (atoms%ntype,atoms%lmaxd,2,2)
    REAL,    INTENT (IN) :: rsopdpd(atoms%ntype,atoms%lmaxd,2,2)
    REAL,    INTENT (IN) :: rsoplop (atoms%ntype,atoms%nlod,2,2)
    REAL,    INTENT (IN) :: rsoplopd(atoms%ntype,atoms%nlod,2,2)
    REAL,    INTENT (IN) :: rsopdplo(atoms%ntype,atoms%nlod,2,2)
    REAL,    INTENT (IN) :: rsopplo (atoms%ntype,atoms%nlod,2,2)
    REAL,    INTENT (IN) :: rsoploplop(atoms%ntype,atoms%nlod,atoms%nlod,2,2)
51 52
    COMPLEX, INTENT (IN) :: soangl(atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2)
    COMPLEX, INTENT (OUT) :: zso(:,:,:)!(dimension%nbasfcn,2*dimension%neigd,wannierspin)
53
    REAL,    INTENT (OUT) :: eig_so(2*DIMENSION%neigd)
54 55 56 57 58 59 60 61
    !-odim
    !+odim
    !     ..
    !     .. Local Scalars ..
    REAL      r2
    INTEGER   i,i1 ,j,jsp,jsp1,k,ne,nn,nn1,nrec,info
    INTEGER   idim_c,idim_r,jsp2,nbas,j1
    CHARACTER vectors 
62
    LOGICAL   l_socvec,l_qsgw,l_open,l_real
63 64 65 66 67
    INTEGER   irec,irecl_qsgw
    COMPLEX   cdum
    !     ..
    !     .. Local Arrays ..
    INTEGER :: nsz(2)
68
    REAL    :: eig(DIMENSION%neigd,DIMENSION%jspd),s(3)
69 70 71 72 73 74
    REAL,   ALLOCATABLE :: rwork(:)
    COMPLEX,ALLOCATABLE :: cwork(:),chelp(:,:,:,:,:)
    COMPLEX,ALLOCATABLE :: ahelp(:,:,:,:,:),bhelp(:,:,:,:,:)
    COMPLEX,ALLOCATABLE :: zhelp1(:,:),zhelp2(:,:)
    COMPLEX,ALLOCATABLE :: hso(:,:),hsomtx(:,:,:,:)
    COMPLEX,ALLOCATABLE :: sigma_xc_apw(:,:),sigma_xc(:,:)
75 76
   
    TYPE(t_zMAT)::zmat(dimension%jspd)
77 78 79 80 81 82 83 84 85 86 87 88
    !     ..
    !     .. External Subroutines ..
    EXTERNAL CPP_LAPACK_cheev
    !     ..
    !     .. External Functions ..
    COMPLEX  CPP_BLAS_cdotu,CPP_BLAS_cdotc
    EXTERNAL CPP_BLAS_cdotu,CPP_BLAS_cdotc
    !     ..

    !     read from eigenvalue and -vector file
    !

Daniel Wortmann's avatar
Daniel Wortmann committed
89
    l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u>0)
90
    zmat%l_real=l_real
91
    zMat(1:dimension%jspd)%nbasfcn=lapw%nv(1:dimension%jspd)+atoms%nlotot
92 93
    zmat%nbands=dimension%neigd
   
94 95
    INQUIRE (4649,opened=l_socvec)
    INQUIRE (file='fleur.qsgw',exist=l_qsgw)
96
    if (l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
97
       ALLOCATE (zmat(1)%z_r(zmat(1)%nbasfcn,DIMENSION%neigd) )
98 99
       zmat(1)%z_r(:,:)= 0.  
       if (size(zmat)==2)THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
100
          ALLOCATE(zmat(2)%z_r(zmat(2)%nbasfcn,DIMENSION%neigd) )
101 102
          zmat(2)%z_r=0.0
       ENDIF
103
    else
Daniel Wortmann's avatar
Daniel Wortmann committed
104
       ALLOCATE (zmat(1)%z_c(zmat(1)%nbasfcn,DIMENSION%neigd) )
105 106
       zmat(1)%z_c(:,:)= 0.  
       if (size(zmat)==2)THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
107
          ALLOCATE(zmat(2)%z_c(zmat(2)%nbasfcn,DIMENSION%neigd) )
108 109
          zmat(2)%z_c=0.0
       ENDIF  
110
    endif
111 112 113 114
    zso(:,:,:)= CMPLX(0.,0.)

    DO jsp = 1,input%jspins
       CALL read_eig(&
115
            eig_id,nk,jsp, neig=ne,eig=eig(:,jsp))
116
       CALL read_eig(&
117 118
               eig_id,nk,jsp,&
               n_start=1,n_end=ne,&
119
               zmat=zmat(jsp))
120 121
       write(6,*) "jspin=",jsp,",nk=",nk
       write(6,"(5f12.4)") eig(:ne,jsp)	       
122 123 124

       ! write(*,*) 'process',irank,' reads ',nk

Daniel Wortmann's avatar
Daniel Wortmann committed
125 126 127 128 129
!!$       DO i = 1, lapw%nv(1)
!!$          s = lapw%bkpt +lapw%gvec(:,i,1)
!!$          r2 = DOT_PRODUCT(s,MATMUL(s,cell%bbmat))
!!$          lapw%rk(i,1) = SQRT(r2)
!!$       ENDDO
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151

       IF (ne.GT.DIMENSION%neigd) THEN
          WRITE (6,'(a13,i4,a8,i4)') 'alineso: ne=',ne,' > dimension%neigd=',DIMENSION%neigd
          CALL juDFT_error("alineso: ne > neigd",calledby="alineso")
       ENDIF
       nsz(jsp) = ne
    ENDDO
    !
    ! set up size of e.v. problem in second variation: nsize
    !
    nsize = 0
    DO jsp = 1,input%jspins
       IF (input%jspins.EQ.1) THEN
          nsize = 2*nsz(jsp)
          nsz(2) = nsz(1)
       ELSE
          nsize = nsize + nsz(jsp)
       ENDIF
    ENDDO
    !
    ! set up A and B coefficients
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
152 153 154
    ALLOCATE ( ahelp(-atoms%lmaxd:atoms%lmaxd,atoms%lmaxd,atoms%nat,DIMENSION%neigd,DIMENSION%jspd) )
    ALLOCATE ( bhelp(-atoms%lmaxd:atoms%lmaxd,atoms%lmaxd,atoms%nat,DIMENSION%neigd,DIMENSION%jspd) )
    ALLOCATE ( chelp(-atoms%llod :atoms%llod, DIMENSION%neigd,atoms%nlod,atoms%nat ,DIMENSION%jspd) )
155 156 157 158
    CALL timestart("alineso SOC: -help") 
    CALL hsohelp(&
         &             DIMENSION,atoms,sym,&
         &             input,lapw,nsz,&
Daniel Wortmann's avatar
Daniel Wortmann committed
159
         &             cell,&
160
         &             zmat,usdus,&
161 162 163 164 165 166 167
         &             zso,noco,oneD,&
         &             ahelp,bhelp,chelp)
    CALL timestop("alineso SOC: -help") 
    !
    ! set up hamilton matrix
    !

168
    CALL timestart("alineso SOC: -ham") 
169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
    ALLOCATE ( hsomtx(2,2,DIMENSION%neigd,DIMENSION%neigd) )
    CALL hsoham(&
         &            atoms,noco,input,nsz,chelp,&
         &            rsoplop,rsoplopd,rsopdplo,rsopplo,rsoploplop,&
         &            ahelp,bhelp,rsopp,rsoppd,rsopdp,rsopdpd,soangl,&
         &            hsomtx)
    DEALLOCATE ( ahelp,bhelp,chelp )
    CALL timestop("alineso SOC: -ham") 
    !
    ! add e.v. on diagonal
    !
    !      write(*,*) '!!!!!!!!!!! remove SOC !!!!!!!!!!!!!!'
    !      hsomtx = 0 !!!!!!!!!!!!
    DO jsp = 1,input%jspins
       DO i = 1,nsz(jsp)
          hsomtx(jsp,jsp,i,i) = hsomtx(jsp,jsp,i,i) +&
               &                           CMPLX(eig(i,jsp),0.)
          IF (input%jspins.EQ.1) THEN
             hsomtx(2,2,i,i) =  hsomtx(2,2,i,i) +&
                  &                           CMPLX(eig(i,jsp),0.)
          ENDIF
       ENDDO
    ENDDO

    !
    !  resort H-matrix 
    !
    ALLOCATE ( hso(2*DIMENSION%neigd,2*DIMENSION%neigd) )
    DO jsp = 1,2
       DO jsp1 = 1,2
          IF (jsp.EQ.1) nn = 0 
          IF (jsp1.EQ.1) nn1 = 0
          IF (jsp.EQ.2) nn = nsz(1)
          IF (jsp1.EQ.2) nn1 = nsz(1)
          !
          DO i = 1,nsz(jsp)
             DO j = 1,nsz(jsp1)
                hso(i+nn,j+nn1) = hsomtx(jsp,jsp1,i,j)
             ENDDO
          ENDDO
          !
       ENDDO
    ENDDO
    DEALLOCATE ( hsomtx )

    !
    !  add Sigma-vxc (QSGW)
    !
    IF( l_qsgw ) THEN
       nbas = lapw%nv(1) + atoms%nlotot
       WRITE(*,'(A,I3,A,I5,A)') 'Read fleur.qsgw  (',nk,',',nbas,')'
       IF( DIMENSION%jspd .EQ. 2 ) STOP 'alineso: GW+noco not implemented.'
       ALLOCATE ( sigma_xc(2*nsz(1),2*nsz(1)) )        
       ALLOCATE ( sigma_xc_apw(nbas,nbas) )
       INQUIRE(667,opened=l_open)
       IF( .NOT.l_open ) THEN
          IF( nk.NE.1 ) STOP 'unit 667 not opened but not at 1st k'
          OPEN(667,file='fleur.qsgw',form='unformatted')
       ELSE IF( nk.EQ.1) THEN
          REWIND(667)
       ENDIF
       DO jsp1 = 1,2
          DO jsp2 = 1,jsp1
             IF(jsp1.EQ.jsp2) THEN
                READ(667) ((sigma_xc_apw(i,j),j=1,i),i=1,nbas)
                DO i = 1,nbas
                   DO j = 1,i-1
                      sigma_xc_apw(j,i) = CONJG(sigma_xc_apw(i,j))
                   ENDDO
                ENDDO
             ELSE
                READ(667) sigma_xc_apw
             ENDIF
             !            write(*,*) 'lo part set to zero!'
             !            sigma_xc_apw(nv+1:,nv+1:) = 0
             i  = nsz(1) * (jsp1-1) + 1 ; i1 = nsz(1) * jsp1
             j  = nsz(1) * (jsp2-1) + 1 ; j1 = nsz(1) * jsp2
246 247
             if (l_real) THEN
                sigma_xc(i:i1,j:j1) = &
248 249
                     &        MATMUL (       TRANSPOSE(zmat(1)%z_r(:nbas,:))  ,&
                     &        MATMUL ( sigma_xc_apw,   zmat(1)%z_r(:nbas,:) ) )
250
else
251
             sigma_xc(i:i1,j:j1) = &
252 253
                  &        MATMUL ( CONJG(TRANSPOSE(zmat(1)%z_c(:nbas,:))) ,&
                  &        MATMUL ( sigma_xc_apw,   zmat(1)%z_c(:nbas,:) ) )
254
          endif
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
             hso(i:i1,j:j1) = hso(i:i1,j:j1) + CONJG(sigma_xc(i:i1,j:j1))
             IF(jsp1.NE.jsp2) THEN
                sigma_xc(j:j1,i:i1) = TRANSPOSE(CONJG(sigma_xc(i:i1,j:j1)))
                hso(j:j1,i:i1) = hso(j:j1,i:i1)+CONJG(sigma_xc(j:j1,i:i1))
             ENDIF
          ENDDO
       ENDDO
       DEALLOCATE ( sigma_xc_apw )
    ENDIF

    !
    ! diagonalize the hamiltonian using library-routines
    !
    idim_c = 4*DIMENSION%neigd
    idim_r = 6*DIMENSION%neigd

    CALL timestart("alineso SOC: -diag") 

    ALLOCATE ( cwork(idim_c),rwork(idim_r) )

    IF (input%eonly) THEN
       vectors= 'N'
    ELSE
       vectors= 'V'
    ENDIF
    CALL CPP_LAPACK_cheev(vectors,'U',nsize,&
         &                      hso,2*DIMENSION%neigd,&
         &                      eig_so,&
         &                      cwork, idim_c, rwork, &
         &                      info)

    IF (info.NE.0) WRITE (6,FMT=8000) info
8000 FORMAT (' AFTER CPP_LAPACK_cheev: info=',i4)
    CALL timestop("alineso SOC: -diag") 

    DEALLOCATE ( cwork,rwork )

    IF (input%eonly) THEN
       IF(l_socvec)  CALL juDFT_error&
            &        ("EONLY set. Vectors not calculated.",calledby ="alineso")
    ELSE
       ALLOCATE ( zhelp2(DIMENSION%neigd,2*DIMENSION%neigd) )
       !
       ! proj. back to G - space: old eigenvector 'z' to new one 'Z'
       !                                 +
       !  s      ---    s                | z(G,j,s) ...   z(ig,i,jsp)
       ! Z (G) = >     z  (G) * C (i,j)  | Z(G,j,s) ... zso(ig,j,jsp)
       !  j      ---    i                | C(i,j)   ... hso(i ,j)
       !          i                      +
       ! reorder new e.w.  in 2x2 spin space : zhelp(,1),zhelp(,2)
       !

       DO i1 = 1,2

          jsp = i1
          jsp2= i1
          IF (input%jspins.EQ.1) jsp = 1
312
          IF (input%jspins.EQ.1 .AND..NOT.(input%l_wann.OR.l_socvec)) jsp2=1
313 314 315 316 317 318 319 320 321 322
          IF (i1.EQ.1) nn = 0
          IF (i1.EQ.2) nn = nsz(1)

          zhelp2(:,:) = 0.d0
          DO j = 1,nsize
             DO i = 1,nsz(jsp)
                zhelp2(i,j) =  CONJG(hso(i+nn,j))
             ENDDO
          ENDDO  ! j

323
          if (l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
324 325
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%nbasfcn,2*dimension%neigd,dimension%neigd,CMPLX(1.d0,0.d0),CMPLX(zmat(jsp)%z_r(:,:)),&
                  zmat(1)%nbasfcn, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(1,1,jsp2),zmat(1)%nbasfcn)
326
          else
Daniel Wortmann's avatar
Daniel Wortmann committed
327 328
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%nbasfcn,2*dimension%neigd,dimension%neigd, CMPLX(1.d0,0.d0),zmat(jsp)%z_c(:,:),&
                  zmat(1)%nbasfcn, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(:,:,jsp2),zmat(1)%nbasfcn)
329
          endif
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349

       ENDDO    !isp

       IF(l_socvec) THEN
          !RS: write SOC vectors to SOCVEC
          WRITE(4649) lapw%nmat,nsize,input%jspins,nsz,2*DIMENSION%neigd,CONJG(hso)
          !CF: write qsgw
          IF(l_qsgw) THEN
             nn = 2*nsz(1)
             sigma_xc = MATMUL ( TRANSPOSE(hso(:nn,:nn)) ,&
                  &                 MATMUL ( sigma_xc , CONJG(hso(:nn,:nn)) ) )
             WRITE(1014) nn
             WRITE(1014) ((sigma_xc(i,j),i=1,j),j=1,nn)
             DEALLOCATE ( sigma_xc )
          ENDIF
       ENDIF

       DEALLOCATE ( zhelp2 )
    ENDIF ! (.NOT.input%eonly)

350
    DEALLOCATE ( hso )
351
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
352
    nmat=lapw%nmat
353 354 355 356
    RETURN

  END SUBROUTINE alineso
END MODULE m_alineso