alineso.F90 12.4 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
9 10
  SUBROUTINE alineso(eig_id,lapw,&
       mpi,DIMENSION,atoms,sym,kpts,&
11 12
       input,noco,cell,oneD, nk, usdus,rsoc,&
       nsize,nmat, eig_so,zso)
13 14

#include"cpp_double.h"
15
    USE m_types
16 17 18 19 20
    USE m_hsohelp
    USE m_hsoham
    USE m_eig66_io, ONLY : read_eig
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(IN)         :: mpi
21
    TYPE(t_lapw),INTENT(IN)        :: lapw
22 23 24 25 26 27 28
    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
29
    TYPE(t_kpts),INTENT(IN)        :: kpts
30
    TYPE(t_usdus),INTENT(IN)       :: usdus
31
    TYPE(t_rsoc),INTENT(IN)        :: rsoc
32 33 34 35
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: eig_id
    INTEGER, INTENT (IN) :: nk 
Daniel Wortmann's avatar
Daniel Wortmann committed
36
    INTEGER, INTENT (OUT):: nsize,nmat
37 38 39
    !     ..
    !     .. Array Arguments ..
    COMPLEX, INTENT (OUT) :: zso(:,:,:)!(dimension%nbasfcn,2*dimension%neigd,wannierspin)
40
    REAL,    INTENT (OUT) :: eig_so(2*DIMENSION%neigd)
41 42 43 44 45 46
    !-odim
    !+odim
    !     ..
    !     .. Local Scalars ..
    REAL      r2
    INTEGER   i,i1 ,j,jsp,jsp1,k,ne,nn,nn1,nrec,info
47
    INTEGER   idim_c,idim_r,jsp2,nbas,j1,ierr
48
    CHARACTER vectors 
49
    LOGICAL   l_socvec,l_qsgw,l_open,l_real
50
    INTEGER   irec,irecl_qsgw
51
    INTEGER nat_l, extra, nat_start, nat_stop
52 53 54 55
    COMPLEX   cdum
    !     ..
    !     .. Local Arrays ..
    INTEGER :: nsz(2)
56
    REAL    :: eig(DIMENSION%neigd,DIMENSION%jspd),s(3)
57 58
    REAL,   ALLOCATABLE :: rwork(:)
    COMPLEX,ALLOCATABLE :: cwork(:),chelp(:,:,:,:,:)
59
    COMPLEX,ALLOCATABLE :: ahelp(:,:,:,:),bhelp(:,:,:,:)
60 61 62
    COMPLEX,ALLOCATABLE :: zhelp1(:,:),zhelp2(:,:)
    COMPLEX,ALLOCATABLE :: hso(:,:),hsomtx(:,:,:,:)
    COMPLEX,ALLOCATABLE :: sigma_xc_apw(:,:),sigma_xc(:,:)
63
   
64
    TYPE(t_mat)::zmat(dimension%jspd)
65 66 67 68 69 70 71 72 73 74 75 76
    !     ..
    !     .. 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
77
    l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u>0)
78
    zmat%l_real=l_real
79 80
    zMat(1:dimension%jspd)%matsize1=lapw%nv(1:dimension%jspd)+atoms%nlotot
    zmat%matsize2=dimension%neigd
81
   
82 83
    INQUIRE (4649,opened=l_socvec)
    INQUIRE (file='fleur.qsgw',exist=l_qsgw)
84
    if (l_real) THEN
85 86
       ALLOCATE (zmat(1)%data_r(zmat(1)%matsize1,DIMENSION%neigd) )
       zmat(1)%data_r(:,:)= 0.  
87
       if (size(zmat)==2)THEN
88 89
          ALLOCATE(zmat(2)%data_r(zmat(2)%matsize1,DIMENSION%neigd) )
          zmat(2)%data_r=0.0
90
       ENDIF
91
    else
92 93
       ALLOCATE (zmat(1)%data_c(zmat(1)%matsize1,DIMENSION%neigd) )
       zmat(1)%data_c(:,:)= 0.  
94
       if (size(zmat)==2)THEN
95 96
          ALLOCATE(zmat(2)%data_c(zmat(2)%matsize1,DIMENSION%neigd) )
          zmat(2)%data_c=0.0
97
       ENDIF  
98
    endif
99 100 101 102
    zso(:,:,:)= CMPLX(0.,0.)

    DO jsp = 1,input%jspins
       CALL read_eig(&
103
            eig_id,nk,jsp, neig=ne,eig=eig(:,jsp))
104 105 106 107
       IF (judft_was_argument("-debugtime")) THEN
          WRITE(6,*) "Non-SOC ev for nk,jsp:",nk,jsp
          WRITE(6,"(6(f10.6,1x))") eig(:ne,jsp)
       ENDIF
108
       CALL read_eig(&
109 110
               eig_id,nk,jsp,&
               n_start=1,n_end=ne,&
111
               zmat=zmat(jsp))
112 113 114

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

115 116 117 118 119
!!$       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
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139

       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
    !
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
    ! distribution of (abc)cof over atoms
    !
!
! in case of ev-parallelization, now distribute the atoms:
!
      IF (mpi%n_size > 1) THEN
        nat_l = FLOOR(real(atoms%nat)/mpi%n_size)
        extra = atoms%nat - nat_l*mpi%n_size
        nat_start = mpi%n_rank*nat_l + 1 + extra
        nat_stop  = (mpi%n_rank+1)*nat_l + extra
        IF (mpi%n_rank < extra) THEN
          nat_start = nat_start - (extra - mpi%n_rank)
          nat_stop  = nat_stop - (extra - mpi%n_rank - 1)
        ENDIF
      ELSE
        nat_start = 1
        nat_stop  = atoms%nat
      ENDIF
      nat_l = nat_stop - nat_start + 1
    !
160 161
    ! set up A and B coefficients
    !
162 163 164
    ALLOCATE ( ahelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,DIMENSION%neigd,input%jspins) )
    ALLOCATE ( bhelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,DIMENSION%neigd,input%jspins) )
    ALLOCATE ( chelp(-atoms%llod :atoms%llod, DIMENSION%neigd,atoms%nlod,nat_l,input%jspins) )
165
    CALL timestart("alineso SOC: -help") 
166
    write(*,*) nat_start,nat_stop,nat_l
167 168 169
    CALL hsohelp(&
         &             DIMENSION,atoms,sym,&
         &             input,lapw,nsz,&
170
         &             cell,&
171
         &             zmat,usdus,&
172
         &             zso,noco,oneD,&
173
         &             nat_start,nat_stop,nat_l,&
174
         &             ahelp,bhelp,chelp)
175
    CALL timestop("alineso SOC: -help") 
176 177 178
    !
    ! set up hamilton matrix
    !
179
    CALL timestart("alineso SOC: -ham") 
180 181 182 183
#ifdef CPP_MPI
    CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
    ALLOCATE ( hsomtx(DIMENSION%neigd,DIMENSION%neigd,2,2) )
184
    CALL hsoham(atoms,noco,input,nsz,dimension%neigd,chelp,rsoc,ahelp,bhelp,&
185 186
                nat_start,nat_stop,mpi%n_rank,mpi%n_size,mpi%SUB_COMM,&
                hsomtx)
187
    write(*,*) 'after hsoham'
188 189
    DEALLOCATE ( ahelp,bhelp,chelp )
    CALL timestop("alineso SOC: -ham") 
190
    IF (mpi%n_rank==0) THEN
191 192 193 194 195 196 197
    !
    ! add e.v. on diagonal
    !
    !      write(*,*) '!!!!!!!!!!! remove SOC !!!!!!!!!!!!!!'
    !      hsomtx = 0 !!!!!!!!!!!!
    DO jsp = 1,input%jspins
       DO i = 1,nsz(jsp)
198
          hsomtx(i,i,jsp,jsp) = hsomtx(i,i,jsp,jsp) +&
199 200
               &                           CMPLX(eig(i,jsp),0.)
          IF (input%jspins.EQ.1) THEN
201
             hsomtx(i,i,2,2) =  hsomtx(i,i,2,2) +&
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217
                  &                           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)
          !
218
          !write(3333,'(2i3,4e15.8)') jsp,jsp1,hsomtx(jsp,jsp1,8,8),hsomtx(jsp,jsp1,32,109) 
219 220
          DO i = 1,nsz(jsp)
             DO j = 1,nsz(jsp1)
221
                hso(i+nn,j+nn1) = hsomtx(i,j,jsp,jsp1)
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
             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
261 262
             if (l_real) THEN
                sigma_xc(i:i1,j:j1) = &
263 264
                     &        MATMUL (       TRANSPOSE(zmat(1)%data_r(:nbas,:))  ,&
                     &        MATMUL ( sigma_xc_apw,   zmat(1)%data_r(:nbas,:) ) )
265
else
266
             sigma_xc(i:i1,j:j1) = &
267 268
                  &        MATMUL ( CONJG(TRANSPOSE(zmat(1)%data_c(:nbas,:))) ,&
                  &        MATMUL ( sigma_xc_apw,   zmat(1)%data_c(:nbas,:) ) )
269
          endif
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 312 313 314 315 316 317 318 319 320 321 322 323 324 325
             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
326
          IF (input%jspins.EQ.1 .AND..NOT.(input%l_wann.OR.l_socvec)) jsp2=1
327 328 329 330 331 332 333 334 335 336
          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

337
          if (l_real) THEN
338 339
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.d0,0.d0),CMPLX(zmat(jsp)%data_r(:,:)),&
                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(1,1,jsp2),zmat(1)%matsize1)
340
          else
341 342
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.d0,0.d0),zmat(jsp)%data_c(:,:),&
                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(:,:,jsp2),zmat(1)%matsize1)
343
          endif
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363

       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)

364
    DEALLOCATE ( hso )
365
    ENDIF ! (n_rank==0)
366
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
367
    nmat=lapw%nmat
368 369 370 371
    RETURN

  END SUBROUTINE alineso
END MODULE m_alineso