alineso.F90 11.6 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,input,noco,&
                     cell,oneD, nk, usdus,rsoc,nsize,nmat, eig_so,zso)
11 12

#include"cpp_double.h"
13
    USE m_types
14 15 16 17 18
    USE m_hsohelp
    USE m_hsoham
    USE m_eig66_io, ONLY : read_eig
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(IN)         :: mpi
19
    TYPE(t_lapw),INTENT(IN)        :: lapw
20 21 22 23 24 25 26
    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
27
    TYPE(t_kpts),INTENT(IN)        :: kpts
28
    TYPE(t_usdus),INTENT(IN)       :: usdus
29
    TYPE(t_rsoc),INTENT(IN)        :: rsoc
30 31 32 33
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: eig_id
    INTEGER, INTENT (IN) :: nk 
Daniel Wortmann's avatar
Daniel Wortmann committed
34
    INTEGER, INTENT (OUT):: nsize,nmat
35 36 37
    !     ..
    !     .. Array Arguments ..
    COMPLEX, INTENT (OUT) :: zso(:,:,:)!(dimension%nbasfcn,2*dimension%neigd,wannierspin)
38
    REAL,    INTENT (OUT) :: eig_so(2*DIMENSION%neigd)
39 40 41 42 43 44
    !-odim
    !+odim
    !     ..
    !     .. Local Scalars ..
    REAL      r2
    INTEGER   i,i1 ,j,jsp,jsp1,k,ne,nn,nn1,nrec,info
45
    INTEGER   idim_c,idim_r,jsp2,nbas,j1,ierr
46
    CHARACTER vectors 
47
    LOGICAL   l_socvec,l_qsgw,l_open,l_real
48
    INTEGER   irec,irecl_qsgw
49
    INTEGER nat_l, extra, nat_start, nat_stop
50 51 52 53
    COMPLEX   cdum
    !     ..
    !     .. Local Arrays ..
    INTEGER :: nsz(2)
54
    REAL    :: eig(DIMENSION%neigd,input%jspins),s(3)
55 56
    REAL,   ALLOCATABLE :: rwork(:)
    COMPLEX,ALLOCATABLE :: cwork(:),chelp(:,:,:,:,:)
57
    COMPLEX,ALLOCATABLE :: ahelp(:,:,:,:),bhelp(:,:,:,:)
58 59 60
    COMPLEX,ALLOCATABLE :: zhelp1(:,:),zhelp2(:,:)
    COMPLEX,ALLOCATABLE :: hso(:,:),hsomtx(:,:,:,:)
    COMPLEX,ALLOCATABLE :: sigma_xc_apw(:,:),sigma_xc(:,:)
61
   
62
    TYPE(t_mat)::zmat(input%jspins)
63 64 65 66 67 68 69
    !     ..
    !     .. External Subroutines ..
    EXTERNAL CPP_LAPACK_cheev

    !     read from eigenvalue and -vector file
    !

Daniel Wortmann's avatar
Daniel Wortmann committed
70
    l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u>0)
71
    zmat%l_real=l_real
72
    zMat(1:input%jspins)%matsize1=lapw%nv(1:input%jspins)+atoms%nlotot
73
    zmat%matsize2=dimension%neigd
74
   
75 76
    INQUIRE (4649,opened=l_socvec)
    INQUIRE (file='fleur.qsgw',exist=l_qsgw)
77
    if (l_real) THEN
78 79
       ALLOCATE (zmat(1)%data_r(zmat(1)%matsize1,DIMENSION%neigd) )
       zmat(1)%data_r(:,:)= 0.  
80
       if (size(zmat)==2)THEN
81 82
          ALLOCATE(zmat(2)%data_r(zmat(2)%matsize1,DIMENSION%neigd) )
          zmat(2)%data_r=0.0
83
       ENDIF
84
    else
85 86
       ALLOCATE (zmat(1)%data_c(zmat(1)%matsize1,DIMENSION%neigd) )
       zmat(1)%data_c(:,:)= 0.  
87
       if (size(zmat)==2)THEN
88 89
          ALLOCATE(zmat(2)%data_c(zmat(2)%matsize1,DIMENSION%neigd) )
          zmat(2)%data_c=0.0
90
       ENDIF  
91
    endif
92 93 94
    zso(:,:,:)= CMPLX(0.,0.)

    DO jsp = 1,input%jspins
95
       CALL read_eig(eig_id,nk,jsp, neig=ne,eig=eig(:,jsp))
96 97 98 99
       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
100
       CALL read_eig(eig_id,nk,jsp,list=[(i,i=1,ne)],zmat=zmat(jsp))
101 102 103

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

104 105 106 107 108
!!$       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
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128

       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
    !
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
    ! 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
148

149
    ! set up A and B coefficients
150 151 152
    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))
153
    CALL timestart("alineso SOC: -help") 
154 155
    CALL hsohelp(DIMENSION,atoms,sym,input,lapw,nsz,cell,zmat,usdus,&
                 zso,noco,oneD,nat_start,nat_stop,nat_l,ahelp,bhelp,chelp)
156
    CALL timestop("alineso SOC: -help") 
157

158
    ! set up hamilton matrix
159
    CALL timestart("alineso SOC: -ham") 
160 161 162
#ifdef CPP_MPI
    CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
163
    ALLOCATE (hsomtx(DIMENSION%neigd,DIMENSION%neigd,2,2))
164
    CALL hsoham(atoms,noco,input,nsz,dimension%neigd,chelp,rsoc,ahelp,bhelp,&
165 166
                nat_start,nat_stop,mpi%n_rank,mpi%n_size,mpi%SUB_COMM,hsomtx)
    DEALLOCATE (ahelp,bhelp,chelp)
167
    CALL timestop("alineso SOC: -ham") 
168
    IF (mpi%n_rank==0) THEN
169

170 171 172 173 174
    ! add e.v. on diagonal
    !      write(*,*) '!!!!!!!!!!! remove SOC !!!!!!!!!!!!!!'
    !      hsomtx = 0 !!!!!!!!!!!!
    DO jsp = 1,input%jspins
       DO i = 1,nsz(jsp)
175
          hsomtx(i,i,jsp,jsp) = hsomtx(i,i,jsp,jsp) + CMPLX(eig(i,jsp),0.)
176
          IF (input%jspins.EQ.1) THEN
177
             hsomtx(i,i,2,2) =  hsomtx(i,i,2,2) + CMPLX(eig(i,jsp),0.)
178 179 180 181 182 183 184
          ENDIF
       ENDDO
    ENDDO

    !
    !  resort H-matrix 
    !
185
    ALLOCATE (hso(2*DIMENSION%neigd,2*DIMENSION%neigd))
186 187 188 189 190 191
    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)
192

193
          !write(3333,'(2i3,4e15.8)') jsp,jsp1,hsomtx(jsp,jsp1,8,8),hsomtx(jsp,jsp1,32,109) 
194 195
          DO i = 1,nsz(jsp)
             DO j = 1,nsz(jsp1)
196
                hso(i+nn,j+nn1) = hsomtx(i,j,jsp,jsp1)
197 198 199 200
             ENDDO
          ENDDO
       ENDDO
    ENDDO
201
    DEALLOCATE (hsomtx)
202 203 204 205

    !
    !  add Sigma-vxc (QSGW)
    !
206
    IF(l_qsgw) THEN
207 208
       nbas = lapw%nv(1) + atoms%nlotot
       WRITE(*,'(A,I3,A,I5,A)') 'Read fleur.qsgw  (',nk,',',nbas,')'
209
       IF( input%jspins .EQ. 2 ) STOP 'alineso: GW+noco not implemented.'
210 211
       ALLOCATE (sigma_xc(2*nsz(1),2*nsz(1)))
       ALLOCATE (sigma_xc_apw(nbas,nbas))
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
       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
235 236
             if (l_real) THEN
                sigma_xc(i:i1,j:j1) = &
237 238
                              MATMUL (       TRANSPOSE(zmat(1)%data_r(:nbas,:))  ,&
                              MATMUL ( sigma_xc_apw,   zmat(1)%data_r(:nbas,:) ) )
239
else
240
             sigma_xc(i:i1,j:j1) = &
241 242
                           MATMUL ( CONJG(TRANSPOSE(zmat(1)%data_c(:nbas,:))) ,&
                           MATMUL ( sigma_xc_apw,   zmat(1)%data_c(:nbas,:) ) )
243
          endif
244 245 246 247 248 249 250
             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
251
       DEALLOCATE (sigma_xc_apw)
252 253 254 255 256 257 258 259 260 261
    ENDIF

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

    CALL timestart("alineso SOC: -diag") 

262
    ALLOCATE (cwork(idim_c),rwork(idim_r))
263 264 265 266 267 268

    IF (input%eonly) THEN
       vectors= 'N'
    ELSE
       vectors= 'V'
    ENDIF
269 270
    CALL CPP_LAPACK_cheev(vectors,'U',nsize,hso,2*DIMENSION%neigd,eig_so,&
                          cwork, idim_c, rwork, info)
271 272 273 274
    IF (info.NE.0) WRITE (6,FMT=8000) info
8000 FORMAT (' AFTER CPP_LAPACK_cheev: info=',i4)
    CALL timestop("alineso SOC: -diag") 

275
    DEALLOCATE (cwork,rwork)
276 277

    IF (input%eonly) THEN
278
       IF(l_socvec) CALL juDFT_error("EONLY set. Vectors not calculated.",calledby ="alineso")
279
    ELSE
280
       ALLOCATE (zhelp2(DIMENSION%neigd,2*DIMENSION%neigd))
281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
       !
       ! 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
296
          IF (input%jspins.EQ.1 .AND..NOT.(input%l_wann.OR.l_socvec)) jsp2=1
297 298 299
          IF (i1.EQ.1) nn = 0
          IF (i1.EQ.2) nn = nsz(1)

Matthias Redies's avatar
Matthias Redies committed
300
          zhelp2(:,:) = 0.0
301 302 303 304 305 306
          DO j = 1,nsize
             DO i = 1,nsz(jsp)
                zhelp2(i,j) =  CONJG(hso(i+nn,j))
             ENDDO
          ENDDO  ! j

307
          if (l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
308 309
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),&
                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
310
          else
Matthias Redies's avatar
Matthias Redies committed
311 312
             CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),&
                  zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(:,:,jsp2),zmat(1)%matsize1)
313
          endif
314 315 316 317 318 319 320 321 322 323 324 325 326

       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)
327
             DEALLOCATE (sigma_xc)
328 329 330
          ENDIF
       ENDIF

331
       DEALLOCATE (zhelp2)
332 333
    ENDIF ! (.NOT.input%eonly)

334
    DEALLOCATE ( hso )
335
    ENDIF ! (n_rank==0)
336
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
337
    nmat=lapw%nmat
338 339 340 341
    RETURN

  END SUBROUTINE alineso
END MODULE m_alineso