LCOV - code coverage report
Current view: top level - eigen - eigen.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 69 102 67.6 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : !<@brief
       7             : !<The eigen routine sets up and solves the generalized eigenvalue problem
       8             : !More description at end of file
       9             : MODULE m_eigen
      10             :    USE m_juDFT
      11             :    IMPLICIT NONE
      12             : CONTAINS
      13             :    !>The eigenvalue problem is constructed and solved in this routine. The following steps are performed:
      14             :    !> 1. Preparation: generate energy parameters, open eig-file
      15             :    !> 2. CALL to mt_setup() : this constructs the local Hamiltonian (i.e. the Hamiltonian in the \f$ u,\dot u, u_{lo} \f$ basis) LDA+U is also added here
      16             :    !> 3. within the (collinear)spin and k-point loop: CALL to eigen_hssetup() to generate the matrices, CALL to eigen_diag() to perform diagonalization
      17             :    !> 4. writing (saving) of eigenvectors
      18             :    !>
      19             :    !> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat.
      20             :    !>@author D. Wortmann
      21         340 :    SUBROUTINE eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,&
      22             :                     cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inden,v,vx)
      23             : 
      24             : #include"cpp_double.h"
      25             :       USE m_types
      26             :       USE m_eigen_hssetup
      27             :       USE m_pot_io
      28             :       USE m_eigen_diag
      29             :       USE m_add_vnonlocal
      30             :       USE m_subvxc
      31             :       !USE m_hsefunctional
      32             :       USE m_mt_setup
      33             :       USE m_util
      34             :       USE m_io_hybrid
      35             :       !USE m_icorrkeys
      36             :       USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig
      37             :       USE m_xmlOutput
      38             : #ifdef CPP_MPI
      39             :       USE m_mpi_bc_potden
      40             : #endif
      41             :       USE m_symmetrize_matrix
      42             :       USE m_unfold_band_kpts !used for unfolding bands
      43             :       USE m_types_mpimat
      44             : 
      45             :       IMPLICIT NONE
      46             :       TYPE(t_results),INTENT(INOUT):: results
      47             :       CLASS(t_xcpot),INTENT(IN)    :: xcpot
      48             :       TYPE(t_mpi),INTENT(IN)       :: mpi
      49             :       TYPE(t_dimension),INTENT(IN) :: DIMENSION
      50             :       TYPE(t_oneD),INTENT(IN)      :: oneD
      51             :       TYPE(t_hybrid),INTENT(INOUT) :: hybrid
      52             :       TYPE(t_enpara),INTENT(INOUT) :: enpara
      53             :       TYPE(t_input),INTENT(IN)     :: input
      54             :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
      55             :       TYPE(t_noco),INTENT(IN)      :: noco
      56             :       TYPE(t_banddos),INTENT(IN)   :: banddos
      57             :       TYPE(t_sym),INTENT(IN)       :: sym
      58             :       TYPE(t_stars),INTENT(IN)     :: stars
      59             :       TYPE(t_cell),INTENT(IN)      :: cell
      60             :       TYPE(t_kpts),INTENT(INOUT)   :: kpts
      61             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
      62             :       TYPE(t_atoms),INTENT(IN)     :: atoms
      63             :       TYPE(t_potden),INTENT(IN)    :: inden,vx
      64             :       TYPE(t_potden),INTENT(INOUT) :: v    !u_setup will modify the potential matrix
      65             : 
      66             : #ifdef CPP_MPI
      67             :       INCLUDE 'mpif.h'
      68             : #endif
      69             : 
      70             : !    EXTERNAL MPI_BCAST    !only used by band_unfolding to broadcast the gvec
      71             : 
      72             :       ! Scalar Arguments
      73             :       INTEGER,INTENT(IN)    :: iter
      74             :       INTEGER,INTENT(IN)    :: eig_id
      75             : 
      76             :       ! Local Scalars
      77             :       INTEGER jsp,nk,nred,ne_all,ne_found
      78             :       INTEGER ne, nk_i
      79             :       INTEGER isp,i,j,err
      80             :       LOGICAL l_wu,l_file,l_real,l_zref
      81             :       INTEGER :: solver=0
      82             :       ! Local Arrays
      83             :       INTEGER              :: ierr
      84         680 :       INTEGER              :: neigBuffer(kpts%nkpt,input%jspins)
      85             : 
      86         680 :       COMPLEX              :: unfoldingBuffer(SIZE(results%unfolding_weights,1),kpts%nkpt,input%jspins) ! needed for unfolding bandstructure mpi case
      87             : 
      88         340 :       REAL,    ALLOCATABLE :: bkpt(:)
      89         340 :       REAL,    ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
      90             : 
      91             :       INTEGER                   :: jsp_m, i_kpt_m, i_m
      92             : 
      93         340 :       TYPE(t_tlmplm)            :: td
      94         340 :       TYPE(t_usdus)             :: ud
      95         340 :       TYPE(t_lapw)              :: lapw
      96         680 :       CLASS(t_mat), ALLOCATABLE :: zMat
      97        1360 :       CLASS(t_mat), ALLOCATABLE :: hmat,smat
      98         340 :       CLASS(t_mat), ALLOCATABLE :: smat_unfold !used for unfolding bandstructure
      99             : 
     100             :       ! Variables for HF or hybrid functional calculation
     101             :       INTEGER                   :: comm(kpts%nkpt),irank2(kpts%nkpt),isize2(kpts%nkpt), dealloc_stat
     102             :       character(len=300)        :: errmsg
     103             : 
     104         340 :       call ud%init(atoms,input%jspins)
     105         340 :       ALLOCATE(eig(DIMENSION%neigd))
     106         340 :       ALLOCATE(bkpt(3))
     107         340 :       ALLOCATE(eigBuffer(DIMENSION%neigd,kpts%nkpt,input%jspins))
     108             : 
     109         340 :       l_real=sym%invs.AND..NOT.noco%l_noco
     110             : 
     111             :       ! check if z-reflection trick can be used
     112         340 :       l_zref=(sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco)
     113         340 :       IF (mpi%n_size > 1) l_zref = .FALSE.
     114             : 
     115             : #ifdef CPP_MPI
     116         340 :       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,v)
     117             : #endif
     118             : 
     119             :       !IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),(/iter,v%iter/),&
     120             :       !                                                RESHAPE((/19,13,5,5/),(/2,2/)))
     121             : 
     122             :       ! Set up and solve the eigenvalue problem
     123             :       !   loop over spins
     124             :       !     set up k-point independent t(l'm',lm) matrices
     125         340 :       CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,v,mpi,results,DIMENSION,td,ud)
     126             : 
     127         948 :       neigBuffer = 0
     128         340 :       results%neig = 0
     129         340 :       results%eig = 1.0e300
     130         948 :       eigBuffer = 1.0e300
     131         948 :       unfoldingBuffer = CMPLX(0.0,0.0)
     132             : 
     133         700 :       DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
     134        1888 :          k_loop:DO nk_i = 1,size(mpi%k_list)
     135        1188 :             nk=mpi%k_list(nk_i)
     136             :             ! Set up lapw list
     137        1188 :             CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
     138        1188 :             call timestart("Setup of H&S matrices")
     139             :             CALL eigen_hssetup(jsp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,sym,&
     140        1188 :                                stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
     141        1188 :             CALL timestop("Setup of H&S matrices")
     142             : 
     143        1188 :             IF(hybrid%l_hybrid.OR.input%l_rdmft) THEN
     144             : 
     145           0 :                DO i = 1, hmat%matsize1
     146           0 :                   DO j = 1, i
     147           0 :                      IF (hmat%l_real) THEN
     148           0 :                         IF ((i.LE.5).AND.(j.LE.5)) THEN
     149           0 :                            WRITE(1233,'(2i7,2f15.8)') i, j, hmat%data_r(i,j), hmat%data_r(j,i)
     150             :                         END IF
     151             :                      ELSE
     152           0 :                         IF ((i.LE.5).AND.(j.LE.5)) THEN
     153           0 :                            WRITE(1233,'(2i7,4f15.8)') i, j, hmat%data_c(i,j), hmat%data_c(j,i)
     154             :                         END IF
     155             :                      ENDIF
     156             :                   END DO
     157             :                END DO
     158             : 
     159             :                ! Write overlap matrix smat to direct access file olap
     160           0 :                print *,"Wrong overlap matrix used, fix this later"
     161           0 :                CALL write_olap(smat,(jsp-1)*kpts%nkpt+nk) ! Note: At this moment this only works without MPI parallelization
     162             :             END IF ! hybrid%l_hybrid.OR.input%l_rdmft
     163             : 
     164        1188 :             IF(hybrid%l_hybrid) THEN
     165           0 :                PRINT *,"TODO"
     166             : !             STOP "TODO"
     167           0 :                PRINT *,"BASIS:", lapw%nv(jsp), atoms%nlotot
     168           0 :                IF (hybrid%l_addhf) CALL add_Vnonlocal(nk,lapw,atoms,hybrid,dimension,kpts,jsp,results,xcpot,noco,hmat)
     169             : 
     170           0 :                IF(hybrid%l_subvxc) THEN
     171             :                   CALL subvxc(lapw,kpts%bk(:,nk),DIMENSION,input,jsp,v%mt(:,0,:,:),atoms,ud,hybrid,enpara%el0,enpara%ello0,&
     172           0 :                               sym,cell,sphhar,stars,xcpot,mpi,oneD,hmat,vx)
     173             :                END IF
     174             :             END IF ! hybrid%l_hybrid
     175             : 
     176        1188 :             l_wu=.FALSE.
     177        1188 :             ne_all=DIMENSION%neigd
     178        1188 :             IF(ne_all < 0) ne_all = lapw%nmat
     179        1188 :             IF(ne_all > lapw%nmat) ne_all = lapw%nmat
     180             : 
     181             :             !Try to symmetrize matrix
     182        1188 :             CALL symmetrize_matrix(mpi,noco,kpts,nk,hmat,smat)
     183             : 
     184        1188 :             IF (banddos%unfoldband) THEN
     185             :                select type(smat)
     186             :                type is (t_mat)
     187           0 :                   allocate(t_mat::smat_unfold)
     188           0 :                   select type(smat_unfold)
     189             :                   type is (t_mat)
     190           0 :                      smat_unfold=smat
     191             :                   end select
     192             :                type is (t_mpimat)
     193           0 :                   allocate(t_mpimat::smat_unfold)
     194           0 :                   select type(smat_unfold)
     195             :                   type is (t_mpimat)
     196           0 :                      smat_unfold=smat
     197             :                   end select
     198             :                end select
     199             :             END IF
     200             : 
     201             :             ! Solve generalized eigenvalue problem.
     202             :             !     ne_all ... number of eigenpairs searched (and found) on this node
     203             :             !                on input, overall number of eigenpairs searched,
     204             :             !                on output, local number of eigenpairs found
     205             :             !     eig ...... all eigenvalues, output
     206             :             !     zMat ..... local eigenvectors, output
     207        1188 :             CALL eigen_diag(solver,hmat,smat,ne_all,eig,zMat,nk,jsp,iter)
     208             :               
     209        1188 :             CALL smat%free()
     210        1188 :             CALL hmat%free()
     211        1188 :             DEALLOCATE(hmat,smat, stat=dealloc_stat, errmsg=errmsg)
     212        1188 :             if(dealloc_stat /= 0) call juDFT_error("deallocate failed for hmat or smat",&
     213           0 :                                                    hint=errmsg, calledby="eigen.F90")
     214             : 
     215             :             ! Output results
     216        1188 :             CALL timestart("EV output")
     217             : #if defined(CPP_MPI)
     218             :             ! Collect number of all eigenvalues
     219        1188 :             ne_found=ne_all
     220        1188 :             CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, mpi%sub_comm,ierr)
     221        1188 :             ne_all=MIN(DIMENSION%neigd,ne_all)
     222             : #else
     223             :             ne_found=ne_all
     224             : #endif
     225        1188 :             IF (.NOT.zMat%l_real) THEN
     226         824 :                zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
     227             :             END IF
     228        1188 :             IF (mpi%n_rank == 0) THEN
     229             :                 ! Only process 0 writes out the value of ne_all and the
     230             :                 ! eigenvalues. 
     231             :                 CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
     232         746 :                            eig(:ne_all),n_start=mpi%n_size,n_end=mpi%n_rank,zMat=zMat)
     233         746 :                 eigBuffer(:ne_all,nk,jsp) = eig(:ne_all)
     234             :             ELSE
     235             :                 CALL write_eig(eig_id, nk,jsp,ne_found,&
     236         442 :                            n_start=mpi%n_size,n_end=mpi%n_rank,zMat=zMat)
     237             :             ENDIF
     238        1188 :             neigBuffer(nk,jsp) = ne_found
     239             : #if defined(CPP_MPI)
     240             :             ! RMA synchronization
     241        1188 :             CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
     242             : #endif
     243        1188 :             CALL timestop("EV output")
     244             : 
     245        1188 :             IF (banddos%unfoldband) THEN
     246           0 :                IF(modulo (kpts%nkpt,mpi%n_size).NE.0) call juDFT_error("number kpts needs to be multiple of number mpi threads",& 
     247           0 :                    hint=errmsg, calledby="eigen.F90")
     248           0 :                CALL calculate_plot_w_n(banddos,cell,kpts,smat_unfold,zMat,lapw,nk,jsp,eig,results,input,atoms,unfoldingBuffer,mpi)
     249           0 :                CALL smat_unfold%free()
     250           0 :                DEALLOCATE(smat_unfold, stat=dealloc_stat, errmsg=errmsg)
     251           0 :                if(dealloc_stat /= 0) call juDFT_error("deallocate failed for smat_unfold",&
     252           0 :                                                       hint=errmsg, calledby="eigen.F90")
     253             :             END IF
     254             : 
     255        1188 :             call zMat%free()
     256        1548 :             deallocate(zMat)
     257             :          END DO  k_loop
     258             :       END DO ! spin loop ends
     259             : 
     260             : #ifdef CPP_MPI
     261         340 :       IF (banddos%unfoldband) THEN
     262           0 :          results%unfolding_weights = CMPLX(0.0,0.0)
     263           0 :        CALL MPI_ALLREDUCE(unfoldingBuffer,results%unfolding_weights,SIZE(results%unfolding_weights,1)*SIZE(results%unfolding_weights,2)*SIZE(results%unfolding_weights,3),CPP_MPI_COMPLEX,MPI_SUM,mpi%mpi_comm,ierr)
     264             :       END IF
     265         340 :       CALL MPI_ALLREDUCE(neigBuffer,results%neig,kpts%nkpt*input%jspins,MPI_INTEGER,MPI_SUM,mpi%mpi_comm,ierr)
     266         340 :       CALL MPI_ALLREDUCE(eigBuffer(:dimension%neigd,:,:),results%eig(:dimension%neigd,:,:),dimension%neigd*kpts%nkpt*input%jspins,MPI_DOUBLE_PRECISION,MPI_MIN,mpi%mpi_comm,ierr)
     267         340 :       CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
     268             : 
     269             : #else
     270             :       results%neig(:,:) = neigBuffer(:,:)
     271             :       results%eig(:dimension%neigd,:,:) = eigBuffer(:dimension%neigd,:,:)
     272             :       results%unfolding_weights(:,:,:) = unfoldingBuffer(:,:,:)
     273             : #endif
     274             : 
     275             :       !IF (hybrid%l_hybrid.OR.hybrid%l_calhf) CALL close_eig(eig_id)
     276             : 
     277         340 :       IF( input%jspins .EQ. 1 .AND. hybrid%l_hybrid ) THEN
     278           0 :          results%te_hfex%valence = 2*results%te_hfex%valence
     279           0 :          IF(hybrid%l_calhf) results%te_hfex%core = 2*results%te_hfex%core
     280             :       END IF
     281         340 :       enpara%epara_min = MINVAL(enpara%el0)
     282         340 :       enpara%epara_min = MIN(MINVAL(enpara%ello0),enpara%epara_min)
     283        1700 :    END SUBROUTINE eigen
     284           0 : END MODULE m_eigen
     285             : 

Generated by: LCOV version 1.13