search

band-unfolding / bandup

BandUP: Band Unfolding code for Plane-wave based calculations
http://www.ifm.liu.se/theomod/compphys/band-unfolding
GNU General Public License v3.0
98 stars 50 forks source link

spin polarized cased #4

Open mikeatm opened 6 years ago

mikeatm commented 6 years ago

Looking at your code, here and here it seems spin polarized cases are not implemented, (else this would segfault: bandup/src/external/espresso-5.1_modules_for_BandUP/Modules/qexml.f90 line 4026,ik_eff would lead to an overflow, going beyond number of kpoints, meaning this was considered but not implemented), what is needed to introduce proper support for spin polarized cases?

Thc02 commented 5 years ago

Yeah, It seems to be that BandUp does not work with QE spin polarized calculations. I have got an error during unfolding process.

dt1125 commented 5 years ago

It seems that BANDUP code is not implemented for spin polarized calculations.

paulovcmedeiros commented 5 years ago

Hello!

First, sorry for not replying very quickly to this. I am working on other projects now and it has been extremely difficult to maintain the code without help (this is actually why I made the code free and open from the beginning, so people in the community could not only use it, but also collaborate with coding).

Back to your question: I use QE's own routines to read QE wavefunctions (I didn't write those myself). In general, however, you select the spin channel to be parsed by bandup by using the command line option "-spin_channel", e.g., "./bandup -spin_channel 2" (you can run the code with the "-h" option to get a list of all supported command line arguments).

Does that work? If not, then I can assist anyone interested in helping fixing this.

Thc02 commented 5 years ago

Indeed, I run band up with the spin channel option but i got: ########################################################################################################################

FROM IOTK LIBRARY, VERSION 1.2.0

UNRECOVERABLE ERROR (ierr=1)

ERROR IN: iotk_scan_end (iotk_scan.f90:241)

CVS Revision: 1.23

foundl

########################################################################################################################

I recompiled the iotk module in QE but this error persist. I will be glad to collaborate solving this issue and extend this code to other ab-initio packages.

i31june commented 4 years ago

Replacing the subroutine "read_qe_evec_file" in the read_qe_wavefunctions_mod.f90 by the following lines would fix the issue involving the spin-polarized case. Best, Pilkwang Kim, SNU, Korea


subroutine read_qe_evc_file(wf, args, ikpt, read_coeffs, iostat)
implicit none
! Mandatory input and output variables
type(pw_wavefunction), intent(inout) :: wf
type(comm_line_args), intent(in) :: args
integer, intent(in) :: ikpt
! Optional input and output variables
logical, optional, intent(in) :: read_coeffs
integer, optional, intent(out) :: iostat
! Local variables
integer, dimension(:,:), allocatable :: G_frac
integer :: nkpts, n_pw, n_bands, n_bands_up, n_bands_down, n_spin, n_spinor, &
           input_file_unit, ios, alloc_stat, i, j, ipw, i_spinor
real(kind=dp) :: e_fermi, ef_up, ef_dw, alat, encut
real(kind=dp), dimension(1:3,1:3) :: A_matrix, B_matrix
real(kind=dp), dimension(:,:), allocatable :: cart_coords_all_kpts, eigenvales, occupations
logical :: is_spinor, two_efs, read_coefficients
character(len=256) :: prefix, outdir, xml_data_file_folder, xml_data_file, &
                      length_units_kpts, length_units, e_units, encut_units

integer :: jkpt

    if(present(iostat)) iostat = -1
    read_coefficients = .TRUE.
    if(present(read_coeffs)) read_coefficients = read_coeffs

    input_file_unit = available_io_unit()
    outdir = trim(adjustl(args%qe_outdir))
    prefix = trim(adjustl(args%qe_prefix))
    xml_data_file_folder = trim(outdir) // '/' // trim(prefix) // '.save'
    xml_data_file = trim(xml_data_file_folder) // "/data-file.xml"

    call qexml_init(input_file_unit, dir=xml_data_file_folder) 
    call qexml_openfile(xml_data_file, "read", ierr=ios)
    if(ios/=0)then
        write(*,'(3A)')'ERROR (read_qe_evc_file): Problems opening the XML data-file "', &
                        trim(adjustl(xml_data_file)),'"!'
        return
    endif

    call qexml_read_bands_info(nbnd=n_bands, nbnd_up=n_bands_up, nbnd_down=n_bands_down, &
                             nspin=n_spin, noncolin=is_spinor, &
                             ef=e_fermi, two_fermi_energies=two_efs, &
                             ef_up=ef_up, ef_dw=ef_dw, &
                             energy_units=e_units, ierr=ios)

    call qexml_read_bz(num_k_points=nkpts, ierr=ios)
    if(n_spin==2) nkpts=nkpts*2
    if(wf%i_spin==2) then
    jkpt=ikpt+nkpts/2
    else
    jkpt=ikpt
    endif

    if(ios/=0)then
        write(*,'(A)')'ERROR (read_qe_evc_file): Could not determine the number of kpts!'
        return
    endif

    call qexml_read_cell(alat=alat, a_units=length_units, &
                         a1=A_matrix(1,:), a2=A_matrix(2,:), a3=A_matrix(3,:), &
                         ierr=ios)
    if(ios/=0)then
        write(*,'(A)')"ERROR (read_qe_evc_file): Could not read lattice info!"
        return
    endif

    ! Converting length units. BandUP works with Angstroms.
    alat = to_angstrom(alat, length_units)
    do j=1,3
        do i=1,3
            A_matrix(i,j) = to_angstrom(A_matrix(i,j), length_units)
        enddo
    enddo
    length_units = 'Angstrom'
    call get_rec_latt(latt=A_matrix, rec_latt=B_matrix)

    ! Getting list of Kpts in cartesian coordinates
    deallocate(cart_coords_all_kpts, stat=alloc_stat)
    allocate(cart_coords_all_kpts(1:3, 1:nkpts))
    call qexml_read_bz(xk=cart_coords_all_kpts, k_units=length_units_kpts, ierr=ios)
    if(ios/=0) then
        write(*,'(A)')'ERROR (read_qe_evc_file): Could not read the list of kpts!'
        return
    endif

    if(n_spin==2) then 
    do i=1,nkpts/2
    cart_coords_all_kpts(:,i+nkpts/2)=cart_coords_all_kpts(:,i)
    enddo
    endif
    cart_coords_all_kpts(:,:) = (twopi/alat) * cart_coords_all_kpts(:,:)

    if(ios/=0)then
        write(*,'(A)')'ERROR (read_qe_evc_file): Could not read info about bands!'
        return
    endif

    deallocate(eigenvales, stat=alloc_stat)
    deallocate(occupations, stat=alloc_stat)
    allocate(eigenvales(1:n_bands, 1:nkpts), occupations(1:n_bands, 1:nkpts))
    if(n_spin==2) then
    call qexml_read_bands_pw(num_k_points=nkpts/2, nkstot=nkpts, lsda=(n_spin==2), nbnd=n_bands, &
                             lkpoint_dir=.TRUE., filename='default', &
                             et=eigenvales, wg=occupations , ierr=ios)
    else 
    call qexml_read_bands_pw(num_k_points=nkpts, nkstot=nkpts, lsda=(n_spin==2), nbnd=n_bands, &
                             lkpoint_dir=.TRUE., filename='default', &
                             et=eigenvales, wg=occupations , ierr=ios)
    endif
    if(ios/=0)then
        write(*,'(A)')'ERROR (read_qe_evc_file): Could not read eigenvalues and occupations!'
        return
    endif

    call qexml_read_planewaves(ecutwfc=encut, cutoff_units=encut_units, ierr=ios)
    if(ios/=0)then
        write(*,'(A)')"ERROR (read_qe_evc_file): Could not read general info regarding wavefunctions!"
        return
    endif

    n_spinor = 1
    if(is_spinor) n_spinor = 2

    if(read_coefficients)then
        ! Getting number of pw for current kpt
        call qexml_read_gk(ik=ikpt, npwk=n_pw, ierr=ios)
        if(ios/=0)then
            write(*,'(A)')"ERROR (read_qe_evc_file): Could not determine the number of plane-waves for the current Kpt."
            return
        endif
        wf%n_pw = n_pw

        ! Reading the G vectors for the current kpt
        deallocate(G_frac, stat=alloc_stat)
        allocate(G_frac(1:3,1:n_pw))
        call qexml_read_gk(ik=ikpt, igk=G_frac, ierr=ios)
        if(ios/=0)then
            write(*,'(A)')"ERROR (read_qe_evc_file): Could not read G-vectors for the current Kpt."
            return
        endif

        ! Reading pw coeffs
        deallocate(wf%pw_coeffs, stat=alloc_stat)
        allocate(wf%pw_coeffs(1:n_spinor, 1:n_pw, 1:n_bands))
        do i_spinor=1, n_spinor
            call get_pw_coeffs_from_evc_file(ikpt, xml_data_file_folder, &
                                             wf, ios, spinor_comp=i_spinor)
            if(ios/=0)then
                write(*,'(A)')"ERROR (read_qe_evc_file): Could not read the plane-wave coefficients for the current KPT."
                return
            endif
        enddo
    endif

    ! Done. Closing file now.
    call qexml_closefile("read", ierr=ios)
    if(ios/=0)then
        write(*,'(A)')"WARNING (read_qe_evc_file): Error closing xml data-file."
    endif

    ! Transferring data to the wf object
    wf%is_spinor = is_spinor
    wf%n_spinor = n_spinor
    wf%n_spin = n_spin
    wf%n_bands = n_bands
    wf%n_bands_up = n_bands_up
    wf%n_bands_down = n_bands_down

    wf%encut = to_ev(encut, encut_units)
    wf%e_fermi = to_ev(e_fermi, e_units)
    wf%e_fermi_up = to_ev(ef_up, e_units)
    wf%e_fermi_down = to_ev(ef_dw, e_units)
    wf%two_efs = two_efs
    if(.not. two_efs)then
        wf%e_fermi_up = wf%e_fermi
        wf%e_fermi_down = wf%e_fermi
    endif

    deallocate(wf%band_occupations, wf%band_energies, stat=alloc_stat)
    allocate(wf%band_occupations(1:n_bands), wf%band_energies(1:n_bands))
    wf%band_occupations(:) = occupations(:,jkpt)
    do j=1, size(eigenvales, dim=1)
        wf%band_energies(j) = to_ev(eigenvales(j,jkpt), e_units)
    enddo

    wf%A_matrix = A_matrix
    wf%B_matrix = B_matrix
    wf%Vcell = abs(triple_product(A_matrix(1,:), A_matrix(2,:), A_matrix(3,:)))
    wf%kpt_cart_coords = cart_coords_all_kpts(:,jkpt)
    wf%kpt_frac_coords = coords_cart_vec_in_new_basis(wf%kpt_cart_coords, new_basis=B_matrix)

    if(read_coefficients)then
        ! To avoid allocating too much memory, 
        ! the plane=wave coeffs are read directly into the wf

        ! Transferring G-vecs
        deallocate(wf%G_frac, wf%G_cart, stat=alloc_stat)
        allocate(wf%G_frac(1:n_pw), wf%G_cart(1:n_pw))
        do ipw=1, n_pw
            wf%G_frac(ipw)%coord(:) = G_frac(:,ipw)
            wf%G_cart(ipw)%coord(:) = 0.0_dp
            do j=1,3
                wf%G_cart(ipw)%coord(:) = wf%G_cart(ipw)%coord(:) + &
                                          real(G_frac(j,ipw), kind=dp)*B_matrix(j,:)
            enddo
        enddo
        deallocate(G_frac, stat=alloc_stat)
    endif

    if(present(iostat)) iostat = 0

end subroutine read_qe_evc_file
ahzeeshan commented 3 years ago

Related to this issue: is bandup implemented for non-collinear spin case?

stepan-tsirkin commented 3 years ago

Dear @ahzeeshan, yes, bandup is implemented for the non-collinear spin case, when the two spin channels are mixed by spin-orbit coupling (wavefunctions are spinors) That was the snece of the paper http://dx.doi.org/10.1103/PhysRevB.91.041116 . When the spin channels are separated (spin-polarised) the unfolding procedure is the same as for scalar (spinless) wavefunctions. The only problem could be in reading wavefunctions and consistency with the particular version of the abinitio code.

stepan-tsirkin commented 2 years ago

Spin pilarized calculations with QE can be performed with banduppy https://github.com/band-unfolding/banduppy using the parameter spin_channel of BandStructure class. If needed, such option may be easily added for VASP.