LCOV - code coverage report
Current view: top level - eigen_soc - alineso.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 113 155 72.9 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_alineso
       2             :   USE m_juDFT
       3             :   !---------------------------------------------------------------------- 
       4             :   ! Set up SO-hamiltonian for 2nd variation (hsohelp and hsoham) and 
       5             :   ! diagonalize by lapack routines.
       6             :   ! Eigenvalues and vectors (eig_so and zso) are returned 
       7             :   !----------------------------------------------------------------------
       8             : CONTAINS
       9          40 :   SUBROUTINE alineso(eig_id,lapw,mpi,DIMENSION,atoms,sym,kpts,input,noco,&
      10          40 :                      cell,oneD, nk, usdus,rsoc,nsize,nmat, eig_so,zso)
      11             : 
      12             : #include"cpp_double.h"
      13             :     USE m_types
      14             :     USE m_hsohelp
      15             :     USE m_hsoham
      16             :     USE m_eig66_io, ONLY : read_eig
      17             :     IMPLICIT NONE
      18             :     TYPE(t_mpi),INTENT(IN)         :: mpi
      19             :     TYPE(t_lapw),INTENT(IN)        :: lapw
      20             :     TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      21             :     TYPE(t_oneD),INTENT(IN)        :: oneD
      22             :     TYPE(t_input),INTENT(IN)       :: input
      23             :     TYPE(t_noco),INTENT(IN)        :: noco
      24             :     TYPE(t_sym),INTENT(IN)         :: sym
      25             :     TYPE(t_cell),INTENT(IN)        :: cell
      26             :     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             :     !     .. Scalar Arguments ..
      32             :     INTEGER, INTENT (IN) :: eig_id
      33             :     INTEGER, INTENT (IN) :: nk 
      34             :     INTEGER, INTENT (OUT):: nsize,nmat
      35             :     !     ..
      36             :     !     .. Array Arguments ..
      37             :     COMPLEX, INTENT (OUT) :: zso(:,:,:)!(dimension%nbasfcn,2*dimension%neigd,wannierspin)
      38             :     REAL,    INTENT (OUT) :: eig_so(2*DIMENSION%neigd)
      39             :     !-odim
      40             :     !+odim
      41             :     !     ..
      42             :     !     .. Local Scalars ..
      43             :     REAL      r2
      44             :     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             :     COMPLEX   cdum
      51             :     !     ..
      52             :     !     .. Local Arrays ..
      53             :     INTEGER :: nsz(2)
      54          80 :     REAL    :: eig(DIMENSION%neigd,input%jspins),s(3)
      55          40 :     REAL,   ALLOCATABLE :: rwork(:)
      56          40 :     COMPLEX,ALLOCATABLE :: cwork(:),chelp(:,:,:,:,:)
      57          40 :     COMPLEX,ALLOCATABLE :: ahelp(:,:,:,:),bhelp(:,:,:,:)
      58          40 :     COMPLEX,ALLOCATABLE :: zhelp1(:,:),zhelp2(:,:)
      59          40 :     COMPLEX,ALLOCATABLE :: hso(:,:),hsomtx(:,:,:,:)
      60          40 :     COMPLEX,ALLOCATABLE :: sigma_xc_apw(:,:),sigma_xc(:,:)
      61             :    
      62         460 :     TYPE(t_mat)::zmat(input%jspins)
      63             :     !     ..
      64             :     !     .. External Subroutines ..
      65             :     EXTERNAL CPP_LAPACK_cheev
      66             : 
      67             :     !     read from eigenvalue and -vector file
      68             :     !
      69             : 
      70          40 :     l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u>0)
      71         116 :     zmat%l_real=l_real
      72          40 :     zMat(1:input%jspins)%matsize1=lapw%nv(1:input%jspins)+atoms%nlotot
      73          40 :     zmat%matsize2=dimension%neigd
      74             :    
      75          40 :     INQUIRE (4649,opened=l_socvec)
      76          40 :     INQUIRE (file='fleur.qsgw',exist=l_qsgw)
      77          40 :     if (l_real) THEN
      78           4 :        ALLOCATE (zmat(1)%data_r(zmat(1)%matsize1,DIMENSION%neigd) )
      79          68 :        zmat(1)%data_r(:,:)= 0.  
      80           4 :        if (size(zmat)==2)THEN
      81           4 :           ALLOCATE(zmat(2)%data_r(zmat(2)%matsize1,DIMENSION%neigd) )
      82          68 :           zmat(2)%data_r=0.0
      83             :        ENDIF
      84             :     else
      85          36 :        ALLOCATE (zmat(1)%data_c(zmat(1)%matsize1,DIMENSION%neigd) )
      86        1448 :        zmat(1)%data_c(:,:)= 0.  
      87          36 :        if (size(zmat)==2)THEN
      88          32 :           ALLOCATE(zmat(2)%data_c(zmat(2)%matsize1,DIMENSION%neigd) )
      89        1180 :           zmat(2)%data_c=0.0
      90             :        ENDIF  
      91             :     endif
      92         116 :     zso(:,:,:)= CMPLX(0.,0.)
      93             : 
      94         116 :     DO jsp = 1,input%jspins
      95          76 :        CALL read_eig(eig_id,nk,jsp, neig=ne,eig=eig(:,jsp))
      96          76 :        IF (judft_was_argument("-debugtime")) THEN
      97           0 :           WRITE(6,*) "Non-SOC ev for nk,jsp:",nk,jsp
      98           0 :           WRITE(6,"(6(f10.6,1x))") eig(:ne,jsp)
      99             :        ENDIF
     100          76 :        CALL read_eig(eig_id,nk,jsp,list=[(i,i=1,ne)],zmat=zmat(jsp))
     101             : 
     102             :        ! write(*,*) 'process',irank,' reads ',nk
     103             : 
     104             : !!$       DO i = 1, lapw%nv(1)
     105             : !!$          s = lapw%bkpt +lapw%gvec(:,i,1)
     106             : !!$          r2 = DOT_PRODUCT(s,MATMUL(s,cell%bbmat))
     107             : !!$          lapw%rk(i,1) = SQRT(r2)
     108             : !!$       ENDDO
     109             : 
     110          76 :        IF (ne.GT.DIMENSION%neigd) THEN
     111           0 :           WRITE (6,'(a13,i4,a8,i4)') 'alineso: ne=',ne,' > dimension%neigd=',DIMENSION%neigd
     112           0 :           CALL juDFT_error("alineso: ne > neigd",calledby="alineso")
     113             :        ENDIF
     114         116 :        nsz(jsp) = ne
     115             :     ENDDO
     116             :     !
     117             :     ! set up size of e.v. problem in second variation: nsize
     118             :     !
     119          40 :     nsize = 0
     120         116 :     DO jsp = 1,input%jspins
     121         116 :        IF (input%jspins.EQ.1) THEN
     122           4 :           nsize = 2*nsz(jsp)
     123           4 :           nsz(2) = nsz(1)
     124             :        ELSE
     125          72 :           nsize = nsize + nsz(jsp)
     126             :        ENDIF
     127             :     ENDDO
     128             :     !
     129             :     ! distribution of (abc)cof over atoms
     130             :     !
     131             : !
     132             : ! in case of ev-parallelization, now distribute the atoms:
     133             : !
     134          40 :       IF (mpi%n_size > 1) THEN
     135           4 :         nat_l = FLOOR(real(atoms%nat)/mpi%n_size)
     136           4 :         extra = atoms%nat - nat_l*mpi%n_size
     137           4 :         nat_start = mpi%n_rank*nat_l + 1 + extra
     138           4 :         nat_stop  = (mpi%n_rank+1)*nat_l + extra
     139           4 :         IF (mpi%n_rank < extra) THEN
     140           2 :           nat_start = nat_start - (extra - mpi%n_rank)
     141           2 :           nat_stop  = nat_stop - (extra - mpi%n_rank - 1)
     142             :         ENDIF
     143             :       ELSE
     144          36 :         nat_start = 1
     145          36 :         nat_stop  = atoms%nat
     146             :       ENDIF
     147          40 :       nat_l = nat_stop - nat_start + 1
     148             : 
     149             :     ! set up A and B coefficients
     150          40 :     ALLOCATE (ahelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,DIMENSION%neigd,input%jspins))
     151          40 :     ALLOCATE (bhelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,DIMENSION%neigd,input%jspins))
     152          40 :     ALLOCATE (chelp(-atoms%llod :atoms%llod, DIMENSION%neigd,atoms%nlod,nat_l,input%jspins))
     153          40 :     CALL timestart("alineso SOC: -help") 
     154             :     CALL hsohelp(DIMENSION,atoms,sym,input,lapw,nsz,cell,zmat,usdus,&
     155          40 :                  zso,noco,oneD,nat_start,nat_stop,nat_l,ahelp,bhelp,chelp)
     156          40 :     CALL timestop("alineso SOC: -help") 
     157             : 
     158             :     ! set up hamilton matrix
     159          40 :     CALL timestart("alineso SOC: -ham") 
     160             : #ifdef CPP_MPI
     161          40 :     CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
     162             : #endif
     163          40 :     ALLOCATE (hsomtx(DIMENSION%neigd,DIMENSION%neigd,2,2))
     164             :     CALL hsoham(atoms,noco,input,nsz,dimension%neigd,chelp,rsoc,ahelp,bhelp,&
     165          40 :                 nat_start,nat_stop,mpi%n_rank,mpi%n_size,mpi%SUB_COMM,hsomtx)
     166          40 :     DEALLOCATE (ahelp,bhelp,chelp)
     167          40 :     CALL timestop("alineso SOC: -ham") 
     168          40 :     IF (mpi%n_rank==0) THEN
     169             : 
     170             :     ! add e.v. on diagonal
     171             :     !      write(*,*) '!!!!!!!!!!! remove SOC !!!!!!!!!!!!!!'
     172             :     !      hsomtx = 0 !!!!!!!!!!!!
     173         111 :     DO jsp = 1,input%jspins
     174        2685 :        DO i = 1,nsz(jsp)
     175        2574 :           hsomtx(i,i,jsp,jsp) = hsomtx(i,i,jsp,jsp) + CMPLX(eig(i,jsp),0.)
     176        2647 :           IF (input%jspins.EQ.1) THEN
     177         198 :              hsomtx(i,i,2,2) =  hsomtx(i,i,2,2) + CMPLX(eig(i,jsp),0.)
     178             :           ENDIF
     179             :        ENDDO
     180             :     ENDDO
     181             : 
     182             :     !
     183             :     !  resort H-matrix 
     184             :     !
     185          38 :     ALLOCATE (hso(2*DIMENSION%neigd,2*DIMENSION%neigd))
     186         114 :     DO jsp = 1,2
     187         418 :        DO jsp1 = 1,2
     188         152 :           IF (jsp.EQ.1) nn = 0 
     189         152 :           IF (jsp1.EQ.1) nn1 = 0
     190         152 :           IF (jsp.EQ.2) nn = nsz(1)
     191         152 :           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        5772 :           DO i = 1,nsz(jsp)
     195      236368 :              DO j = 1,nsz(jsp1)
     196      236216 :                 hso(i+nn,j+nn1) = hsomtx(i,j,jsp,jsp1)
     197             :              ENDDO
     198             :           ENDDO
     199             :        ENDDO
     200             :     ENDDO
     201          38 :     DEALLOCATE (hsomtx)
     202             : 
     203             :     !
     204             :     !  add Sigma-vxc (QSGW)
     205             :     !
     206          38 :     IF(l_qsgw) THEN
     207           0 :        nbas = lapw%nv(1) + atoms%nlotot
     208           0 :        WRITE(*,'(A,I3,A,I5,A)') 'Read fleur.qsgw  (',nk,',',nbas,')'
     209           0 :        IF( input%jspins .EQ. 2 ) STOP 'alineso: GW+noco not implemented.'
     210           0 :        ALLOCATE (sigma_xc(2*nsz(1),2*nsz(1)))
     211           0 :        ALLOCATE (sigma_xc_apw(nbas,nbas))
     212           0 :        INQUIRE(667,opened=l_open)
     213           0 :        IF( .NOT.l_open ) THEN
     214           0 :           IF( nk.NE.1 ) STOP 'unit 667 not opened but not at 1st k'
     215           0 :           OPEN(667,file='fleur.qsgw',form='unformatted')
     216           0 :        ELSE IF( nk.EQ.1) THEN
     217           0 :           REWIND(667)
     218             :        ENDIF
     219           0 :        DO jsp1 = 1,2
     220           0 :           DO jsp2 = 1,jsp1
     221           0 :              IF(jsp1.EQ.jsp2) THEN
     222           0 :                 READ(667) ((sigma_xc_apw(i,j),j=1,i),i=1,nbas)
     223           0 :                 DO i = 1,nbas
     224           0 :                    DO j = 1,i-1
     225           0 :                       sigma_xc_apw(j,i) = CONJG(sigma_xc_apw(i,j))
     226             :                    ENDDO
     227             :                 ENDDO
     228             :              ELSE
     229           0 :                 READ(667) sigma_xc_apw
     230             :              ENDIF
     231             :              !            write(*,*) 'lo part set to zero!'
     232             :              !            sigma_xc_apw(nv+1:,nv+1:) = 0
     233           0 :              i  = nsz(1) * (jsp1-1) + 1 ; i1 = nsz(1) * jsp1
     234           0 :              j  = nsz(1) * (jsp2-1) + 1 ; j1 = nsz(1) * jsp2
     235           0 :              if (l_real) THEN
     236             :                 sigma_xc(i:i1,j:j1) = &
     237             :                               MATMUL (       TRANSPOSE(zmat(1)%data_r(:nbas,:))  ,&
     238           0 :                               MATMUL ( sigma_xc_apw,   zmat(1)%data_r(:nbas,:) ) )
     239             : else
     240             :              sigma_xc(i:i1,j:j1) = &
     241             :                            MATMUL ( CONJG(TRANSPOSE(zmat(1)%data_c(:nbas,:))) ,&
     242           0 :                            MATMUL ( sigma_xc_apw,   zmat(1)%data_c(:nbas,:) ) )
     243             :           endif
     244           0 :              hso(i:i1,j:j1) = hso(i:i1,j:j1) + CONJG(sigma_xc(i:i1,j:j1))
     245           0 :              IF(jsp1.NE.jsp2) THEN
     246           0 :                 sigma_xc(j:j1,i:i1) = TRANSPOSE(CONJG(sigma_xc(i:i1,j:j1)))
     247           0 :                 hso(j:j1,i:i1) = hso(j:j1,i:i1)+CONJG(sigma_xc(j:j1,i:i1))
     248             :              ENDIF
     249             :           ENDDO
     250             :        ENDDO
     251           0 :        DEALLOCATE (sigma_xc_apw)
     252             :     ENDIF
     253             : 
     254             :     !
     255             :     ! diagonalize the hamiltonian using library-routines
     256             :     !
     257          38 :     idim_c = 4*DIMENSION%neigd
     258          38 :     idim_r = 6*DIMENSION%neigd
     259             : 
     260          38 :     CALL timestart("alineso SOC: -diag") 
     261             : 
     262          38 :     ALLOCATE (cwork(idim_c),rwork(idim_r))
     263             : 
     264          38 :     IF (input%eonly) THEN
     265           0 :        vectors= 'N'
     266             :     ELSE
     267          38 :        vectors= 'V'
     268             :     ENDIF
     269             :     CALL CPP_LAPACK_cheev(vectors,'U',nsize,hso,2*DIMENSION%neigd,eig_so,&
     270          38 :                           cwork, idim_c, rwork, info)
     271          38 :     IF (info.NE.0) WRITE (6,FMT=8000) info
     272             : 8000 FORMAT (' AFTER CPP_LAPACK_cheev: info=',i4)
     273          38 :     CALL timestop("alineso SOC: -diag") 
     274             : 
     275          38 :     DEALLOCATE (cwork,rwork)
     276             : 
     277          38 :     IF (input%eonly) THEN
     278           0 :        IF(l_socvec) CALL juDFT_error("EONLY set. Vectors not calculated.",calledby ="alineso")
     279             :     ELSE
     280          38 :        ALLOCATE (zhelp2(DIMENSION%neigd,2*DIMENSION%neigd))
     281             :        !
     282             :        ! proj. back to G - space: old eigenvector 'z' to new one 'Z'
     283             :        !                                 +
     284             :        !  s      ---    s                | z(G,j,s) ...   z(ig,i,jsp)
     285             :        ! Z (G) = >     z  (G) * C (i,j)  | Z(G,j,s) ... zso(ig,j,jsp)
     286             :        !  j      ---    i                | C(i,j)   ... hso(i ,j)
     287             :        !          i                      +
     288             :        ! reorder new e.w.  in 2x2 spin space : zhelp(,1),zhelp(,2)
     289             :        !
     290             : 
     291         114 :        DO i1 = 1,2
     292             : 
     293          76 :           jsp = i1
     294          76 :           jsp2= i1
     295          76 :           IF (input%jspins.EQ.1) jsp = 1
     296          76 :           IF (input%jspins.EQ.1 .AND..NOT.(input%l_wann.OR.l_socvec)) jsp2=1
     297          76 :           IF (i1.EQ.1) nn = 0
     298          76 :           IF (i1.EQ.2) nn = nsz(1)
     299             : 
     300        5620 :           zhelp2(:,:) = 0.0
     301        5620 :           DO j = 1,nsize
     302      236292 :              DO i = 1,nsz(jsp)
     303      236216 :                 zhelp2(i,j) =  CONJG(hso(i+nn,j))
     304             :              ENDDO
     305             :           ENDDO  ! j
     306             : 
     307         114 :           if (l_real) THEN
     308             :              CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),&
     309           8 :                   zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
     310             :           else
     311             :              CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),&
     312          68 :                   zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(:,:,jsp2),zmat(1)%matsize1)
     313             :           endif
     314             : 
     315             :        ENDDO    !isp
     316             : 
     317          38 :        IF(l_socvec) THEN
     318             :           !RS: write SOC vectors to SOCVEC
     319           0 :           WRITE(4649) lapw%nmat,nsize,input%jspins,nsz,2*DIMENSION%neigd,CONJG(hso)
     320             :           !CF: write qsgw
     321           0 :           IF(l_qsgw) THEN
     322           0 :              nn = 2*nsz(1)
     323             :              sigma_xc = MATMUL ( TRANSPOSE(hso(:nn,:nn)) ,&
     324           0 :                   &                 MATMUL ( sigma_xc , CONJG(hso(:nn,:nn)) ) )
     325           0 :              WRITE(1014) nn
     326           0 :              WRITE(1014) ((sigma_xc(i,j),i=1,j),j=1,nn)
     327           0 :              DEALLOCATE (sigma_xc)
     328             :           ENDIF
     329             :        ENDIF
     330             : 
     331          38 :        DEALLOCATE (zhelp2)
     332             :     ENDIF ! (.NOT.input%eonly)
     333             : 
     334          38 :     DEALLOCATE ( hso )
     335             :     ENDIF ! (n_rank==0)
     336             :     !
     337          40 :     nmat=lapw%nmat
     338          40 :     RETURN
     339             : 
     340         268 :   END SUBROUTINE alineso
     341             : END MODULE m_alineso

Generated by: LCOV version 1.13