Skip to content

Commit

Permalink
V2.3.0 - The user can now opt for considering or not considering the …
Browse files Browse the repository at this point in the history 
…atomic positions when cheking the symmetries of the pc and SC. BandUP can now suggest prim. cells for the SC and reference unit cell provided by the user. A new colormap (flame) has been added to the plotting tool (still testing), thanks to Manfred Matena. Some changes in the symmetry-related routines. Major changes in the pre-processing tool to find the SC-Kpts. All changes are backwards compatible.
  • Loading branch information
Paulo V C Medeiros committed Jun 11, 2014
1 parent bb78cb5 commit 9c76479
Show file tree
Hide file tree
Showing 72 changed files with 2,408 additions and 2,307 deletions.
9 changes: 8 additions & 1 deletion README
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ System requirements:
===================================================================================
Please mind that:
===================================================================================
* If you use BandUP (or any part of it) in your work, please be so kind as to
state that you've used the BandUP code and don't forget to cite our paper,
Phys. Rev. B 89, 041407(R) (2014).
* Although I've tried to make the compilation as simple as possible - and it has
indeed worked fine so far with many combinations of compilers, computers, and
operating systems, I cannot guarantee that it will always work smoothly.
Expand All @@ -61,6 +64,10 @@ Please mind that:
using other versions of ifort/icc (and *maybe* even other non-intel compilers),
but I cannot, unfortunatelly, guarantee that it will work in 100% of the cases.
* Of course, I'll try my best to help you out in case you do have problems.
* Last, but not least: Always check the results with a critical eye, specially
if they don't look the way you think they are supposed to. Please notify me if
weird stuff happens and you think it's the code's fault!


Let us know if you have problems!
Thanks! Let us know if you have problems!

10 changes: 7 additions & 3 deletions src/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,15 @@ FCFLAGS += -diag-disable remark
current_dir = $(shell pwd)
BandUP_src_folder = $(current_dir)
spglib_folder = $(BandUP_src_folder)/external/spglib-1.5.2
vasp_related_modules_folder = $(BandUP_src_folder)/vasp_related_modules

# libraries needed for linking
LIBS = $(spglib_folder)/src/.libs/libsymspg.a

# Folders containing needed sources
vpath %.f90 $(spglib_folder)/example
vpath %.f90 $(BandUP_src_folder)
vpath %.f90 $(vasp_related_modules_folder)

# List of executables to be built within the package
SRC = main_BandUP.f90
Expand All @@ -43,9 +45,11 @@ all: $(EXE)
math_mod.o: spglib_f08.o
band_unfolding_routines_mod.o: math_mod.o
read_wavecar_mod.o: math_mod.o general_io_mod.o
io_routines_mod.o: math_mod.o strings_mod.o general_io_mod.o read_wavecar_mod.o
$(OBJ): spglib_f08.o math_mod.o io_routines_mod.o read_wavecar_mod.o band_unfolding_routines_mod.o
$(EXE): spglib_f08.o math_mod.o band_unfolding_routines_mod.o general_io_mod.o read_wavecar_mod.o strings_mod.o io_routines_mod.o $(OBJ)
read_vasp_files_mod.o: math_mod.o general_io_mod.o strings_mod.o
write_vasp_files_mod.o: math_mod.o general_io_mod.o strings_mod.o
io_routines_mod.o: math_mod.o strings_mod.o general_io_mod.o read_vasp_files_mod.o write_vasp_files_mod.o read_wavecar_mod.o
$(OBJ): spglib_f08.o math_mod.o io_routines_mod.o read_wavecar_mod.o read_vasp_files_mod.o write_vasp_files_mod.o band_unfolding_routines_mod.o
$(EXE): spglib_f08.o math_mod.o band_unfolding_routines_mod.o general_io_mod.o read_wavecar_mod.o read_vasp_files_mod.o write_vasp_files_mod.o strings_mod.o io_routines_mod.o $(OBJ)

# ======================================================================
# General rules
Expand Down
111 changes: 60 additions & 51 deletions src/band_unfolding_routines_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,19 +33,20 @@ module band_unfolding

CONTAINS

subroutine get_geom_unfolding_relations(GUR,list_of_SCKPTS,pckpts_to_be_checked,B_matrix_SC,verbose)
subroutine get_geom_unfolding_relations(GUR,list_of_SCKPTS,pckpts_to_be_checked,crystal_SC,verbose)
!! Copyright (C) 2013, 2014 Paulo V. C. Medeiros
implicit none
type(geom_unfolding_relations_for_each_SCKPT), intent(out) :: GUR !! Geometric Unfolding Relations
type(selected_pcbz_directions), intent(in) :: pckpts_to_be_checked
type(vec3d), dimension(:), intent(in) :: list_of_SCKPTS
real(kind=dp), dimension(1:3,1:3), intent(in) :: B_matrix_SC
type(crystal_3D), intent(in) :: crystal_SC
logical, intent(in), optional :: verbose
integer :: nkpts, n_selec_pcbz_dirs, i_SCKPT,ipc_kpt,ieqv_SCKPT,nsym,isym, &
integer :: nkpts, n_selec_pcbz_dirs, i_SCKPT,ipc_kpt,ieqv_SCKPT,isym, &
i_selec_pcbz_dir,i_needed_dirs,alloc_stat
integer, dimension(:), allocatable :: n_dirs_for_EBS_along_pcbz_dir,n_pckpts_dirs
logical, dimension(:,:,:), allocatable :: pc_kpt_already_folded
real(kind=dp), dimension(1:3) :: pc_kpt, current_SCKPT,SCKPT_eqv_to_current_SCKPT,trial_folding_G
real(kind=dp), dimension(1:3,1:3) :: B_matrix_SC
type(star), dimension(:), allocatable :: SKPTS_eqv_to_SKPT
type(symmetry_operation), dimension(:), allocatable :: symops
logical :: print_stuff
Expand All @@ -58,9 +59,11 @@ subroutine get_geom_unfolding_relations(GUR,list_of_SCKPTS,pckpts_to_be_checked,
if(print_stuff)then
write(*,"(A)",advance="no")"Verifying geometric unfolding relations between pcbz and SCBZ wave-vectors... "
endif
call get_symm(nsym=nsym, symops=symops, symprec=1D-6, lattice=B_matrix_SC)
call get_star(star_of_pt=SKPTS_eqv_to_SKPT,points=list_of_SCKPTS,lattice=B_matrix_SC, &
tol_for_vec_equality=tol_for_vec_equality,reduce_to_bz=.TRUE.)
B_matrix_SC = crystal_SC%rec_latt_vecs
call get_star(star_of_pt=SKPTS_eqv_to_SKPT, points=list_of_SCKPTS, crystal=crystal_SC, symm_ops=symops, &
tol_for_vec_equality=default_tol_for_vec_equality, &
symprec=default_symprec, reduce_to_bz=.TRUE., &
try_to_reduce_to_a_pc=.FALSE.) ! fails if try_to_reduce_to_a_pc=.TRUE.
!! Allocating and initializing table
nkpts = size(list_of_SCKPTS)
n_selec_pcbz_dirs = size(pckpts_to_be_checked%selec_pcbz_dir(:))
Expand All @@ -74,20 +77,27 @@ subroutine get_geom_unfolding_relations(GUR,list_of_SCKPTS,pckpts_to_be_checked,
enddo
GUR%n_folding_pckpts = 0
GUR%n_pckpts = 0
deallocate(GUR%SCKPT, stat=alloc_stat)
allocate(GUR%SCKPT(1:nkpts))
deallocate(GUR%SCKPT_used_for_unfolding, stat=alloc_stat)
allocate(GUR%SCKPT_used_for_unfolding(1:nkpts))
do i_SCKPT=1,nkpts
deallocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir, stat=alloc_stat)
allocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(1:n_selec_pcbz_dirs))
GUR%SCKPT_used_for_unfolding(i_SCKPT) = .FALSE.
do i_selec_pcbz_dir=1,n_selec_pcbz_dirs
deallocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir, stat=alloc_stat)
allocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(1:n_dirs_for_EBS_along_pcbz_dir(i_selec_pcbz_dir)))
do i_needed_dirs=1,n_dirs_for_EBS_along_pcbz_dir(i_selec_pcbz_dir)
deallocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt, stat=alloc_stat)
allocate(GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(1:n_pckpts_dirs(i_selec_pcbz_dir)))
do ipc_kpt=1,n_pckpts_dirs(i_selec_pcbz_dir)
GUR%n_pckpts = GUR%n_pckpts + 1
pc_kpt(:) = pckpts_to_be_checked%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%coords(:)
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%coords(:) = pc_kpt(:)
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Scoords(:) = pc_kpt(:)
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%folds = .FALSE.
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Sfolding_vec(:) = 0.0d0
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Sfolding_vec(:) = 0.0_dp
enddo
enddo
enddo
Expand All @@ -109,22 +119,30 @@ subroutine get_geom_unfolding_relations(GUR,list_of_SCKPTS,pckpts_to_be_checked,
do ieqv_SCKPT=1,SKPTS_eqv_to_SKPT(i_SCKPT) % neqv
SCKPT_eqv_to_current_SCKPT(:) = SKPTS_eqv_to_SKPT(i_SCKPT) % eqv_pt(ieqv_SCKPT) % coord(:)
trial_folding_G(:) = pc_kpt(:) - SCKPT_eqv_to_current_SCKPT(:)
if(vec_in_latt(vec=trial_folding_G, latt=B_matrix_SC,tolerance=tol_for_vec_equality))then
if(vec_in_latt(vec=trial_folding_G, latt=B_matrix_SC,tolerance=default_tol_for_vec_equality))then
GUR%SCKPT_used_for_unfolding(i_SCKPT) = .TRUE.
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%folds = .TRUE.
GUR%n_folding_pckpts = GUR%n_folding_pckpts + 1
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%coords_actual_unfolding_K = &
SCKPT_eqv_to_current_SCKPT(:)
isym = SKPTS_eqv_to_SKPT(i_SCKPT) % eqv_pt(ieqv_SCKPT) % symop
pc_kpt = pt_eqv_by_point_group_symop(point=pc_kpt,symops=symops,isym=isym, fractional_coords=.FALSE.,invert_symop=.TRUE.)
! Message from Paulo:
! The prefix "S" means "symmetrized". This is a little trick I came up with
! that allows me to use the coefficients of a SC wavefunction psi(K',n) to
! calculate the spectral weights associated with a SC wavefunction psi(K,n),
! where K' = SK and S is a symmetry operation of the crystal's point group.
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Scoords(:) = pc_kpt(:)
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Sfolding_vec(:) = pc_kpt(:) - current_SCKPT(:)
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%Sfolding_vec(:) = &
pc_kpt(:) - current_SCKPT(:)
pc_kpt_already_folded(i_selec_pcbz_dir,i_needed_dirs,ipc_kpt) = .TRUE.
exit ! It has already folded. No need to keep looking for this pckpt.
else
GUR%SCKPT(i_SCKPT)%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(i_needed_dirs)%pckpt(ipc_kpt)%folds = .FALSE.
endif
enddo ! Loop over all the SCKPTS eqv. to the one on the WF file by symm. ops. of the SC
enddo ! Loop over the pckpts along the dirs. of the pcbz that are eqv. to the selec. ones by symm. ops. of the pcbz and NOT by symm. ops. of the SCBZ
enddo ! Loop over the pckpts along the dirs. of the pcbz that are eqv. to the selec. ones by
! symm. ops. of the pcbz and NOT by symm. ops. of the SCBZ
enddo ! Loop over all dirs. of the pcbz that are eqv. to the selec. ones by symm. ops. of the pcbz and NOT by symm. ops. of the SCBZ
enddo ! Loop over the selected pcbz directions
enddo ! Loop over the SCKPTS found on the wavefunction file
Expand Down Expand Up @@ -187,21 +205,34 @@ subroutine select_coeffs_to_calc_spectral_weights(selected_coeff_indices, &
real(kind=dp), intent(in) :: tol_for_vec_equality
integer :: alloc_stat, ig
real(kind=dp), dimension(1:3) :: SC_G, trial_pc_g,g
real(kind=dp) :: tol

tol = 0.1_dp * tol_for_vec_equality
deallocate(selected_coeff_indices, stat=alloc_stat)
!$omp parallel do &
!$ schedule(guided) default(none) &
!$ shared(nplane,iall_G,B_matrix_SC,folding_G,b_matrix_pc,tol_for_vec_equality,selected_coeff_indices) &
!$ private(ig,g,SC_G,trial_pc_g)
do ig=1, nplane
SC_G = iall_G(1,ig)*B_matrix_SC(1,:) + iall_G(2,ig)*B_matrix_SC(2,:) + iall_G(3,ig)*B_matrix_SC(3,:)
trial_pc_g = SC_G(:) - folding_G(:)
if(vec_in_latt(vec=trial_pc_g, latt=b_matrix_pc, tolerance=tol_for_vec_equality))then
!$omp critical
call append(item=ig,list=selected_coeff_indices) ! pc_g + folding_G has to belong to the
!$omp end critical
endif
do while((.not. allocated(selected_coeff_indices)) .and. &
(tol <= max_tol_for_vec_equality))
tol = 10_dp * tol
!$omp parallel do &
!$ schedule(guided) default(none) &
!$ shared(nplane,iall_G,B_matrix_SC,folding_G,b_matrix_pc,tol,selected_coeff_indices) &
!$ private(ig,g,SC_G,trial_pc_g)
do ig=1, nplane
SC_G = iall_G(1,ig)*B_matrix_SC(1,:) + iall_G(2,ig)*B_matrix_SC(2,:) + iall_G(3,ig)*B_matrix_SC(3,:)
trial_pc_g = SC_G(:) - folding_G(:)
if(vec_in_latt(vec=trial_pc_g, latt=b_matrix_pc, tolerance=tol))then
!$omp critical
call append(item=ig,list=selected_coeff_indices) ! pc_g + folding_G has to belong to the
!$omp end critical
endif
enddo
enddo
if(abs(tol - tol_for_vec_equality) > epsilon(1.0_dp) .and. allocated(selected_coeff_indices))then
write(*,'(A)') &
' WARNING (select_coeffs_to_calc_spectral_weights): Problems selecting the coeffs to calculate the spec. weights for the current pc-kpt.'
write(*,'(2(A,f0.6),A)') &
' The tolerace for testing vector equality had to be increased from ',tol_for_vec_equality,' to ',tol,'.'
write(*,'(A)')' The results might not be entirely correct.'
endif

end subroutine select_coeffs_to_calc_spectral_weights

Expand All @@ -221,15 +252,15 @@ subroutine calc_spectral_weights(spectral_weight,coeff,selected_coeff_indices,ad
n_bands_SC_calculation = size(coeff,dim=2)
deallocate(spectral_weight, stat=alloc_stat)
allocate(spectral_weight(1:n_bands_SC_calculation))
spectral_weight(:) = 0.0d0
spectral_weight(:) = 0.0_dp
!$omp parallel do &
!$ schedule(guided) default(none) &
!$ private(iband,i,icoeff) &
!$ shared(n_bands_SC_calculation,selected_coeff_indices,spectral_weight,coeff)
do iband = 1, n_bands_SC_calculation
do i=1, size(selected_coeff_indices)
icoeff = selected_coeff_indices(i)
spectral_weight(iband) = spectral_weight(iband) + abs(coeff(icoeff,iband))**2.0d0
spectral_weight(iband) = spectral_weight(iband) + abs(coeff(icoeff,iband))**2.0_dp
enddo
enddo

Expand Down Expand Up @@ -257,7 +288,7 @@ subroutine calc_spectral_function(SF_at_pckpt,energies,SC_calc_ener,spectral_wei
n_SC_bands = size(spectral_weight)
deallocate(SF_at_pckpt, stat=alloc_stat)
allocate(SF_at_pckpt(1:size(energies)))
SF_at_pckpt(:) = 0.0d0
SF_at_pckpt(:) = 0.0_dp
!$omp parallel do &
!$ schedule(guided) default(none) &
!$ private(iener,E,i_SC_band,E_SC_band) &
Expand Down Expand Up @@ -294,7 +325,7 @@ subroutine calc_delta_N_pckpt(delta_N_pckpt,energies,SC_calc_ener,spectral_weigh
n_SC_bands = size(spectral_weight)
deallocate(delta_N_pckpt, stat=alloc_stat)
allocate(delta_N_pckpt(1:size(energies)))
delta_N_pckpt(:) = 0.0d0
delta_N_pckpt(:) = 0.0_dp
!$omp parallel do &
!$ schedule(guided) default(none) &
!$ private(iener,E,i_SC_band,E_SC_band) &
Expand All @@ -305,8 +336,8 @@ subroutine calc_delta_N_pckpt(delta_N_pckpt,energies,SC_calc_ener,spectral_weigh
E_SC_band = SC_calc_ener(i_SC_band)
delta_N_pckpt(iener) = delta_N_pckpt(iener) + &
spectral_weight(i_SC_band)*integral_delta_x_minus_x0(x0=E_SC_band, &
lower_lim=E-0.5d0*delta_e, &
upper_lim=E+0.5d0*delta_e, &
lower_lim=E-0.5_dp*delta_e, &
upper_lim=E+0.5_dp*delta_e, &
std_dev=sigma)
enddo
enddo
Expand Down Expand Up @@ -359,7 +390,7 @@ subroutine get_delta_Ns_for_output(delta_N_only_selected_dirs,delta_N_symm_avrgd
do i_selec_pcbz_dir=1,size(delta_N%selec_pcbz_dir(:))
delta_N_only_selected_dirs%pcbz_dir(i_selec_pcbz_dir) = delta_N%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(1)
do ipc_kpt=1, size(delta_N%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(1)%pckpt(:))
avrgd_dNs(:) = 0.0d0
avrgd_dNs(:) = 0.0_dp
allocate(delta_N_symm_avrgd_for_EBS%pcbz_dir(i_selec_pcbz_dir)%pckpt(ipc_kpt)%dN(1:nener))
do i_needed_dirs=1,size(delta_N%selec_pcbz_dir(i_selec_pcbz_dir)%needed_dir(:))
weight = dirs_req_for_symmavgd_EBS_along_pcbz_dir(i_selec_pcbz_dir)%irr_dir(i_needed_dirs)%weight
Expand All @@ -372,28 +403,6 @@ subroutine get_delta_Ns_for_output(delta_N_only_selected_dirs,delta_N_symm_avrgd
end subroutine get_delta_Ns_for_output


subroutine append(item, list)
implicit none
integer, intent(in) :: item
integer, dimension(:), allocatable, intent(inout) :: list
integer, dimension(:), allocatable :: new_list

if(.not.allocated(list))then
allocate(list(1:1))
list(1) = item
else
allocate(new_list(1:size(list)+1))
new_list(1:size(list)) = list
new_list(size(list)+1) = item
deallocate(list)
allocate(list(1:size(new_list)))
list = new_list
deallocate(new_list)
endif

end subroutine append


end module band_unfolding


12 changes: 11 additions & 1 deletion src/general_io_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,19 @@
module general_io
implicit none
PRIVATE
PUBLIC :: available_io_unit
PUBLIC :: available_io_unit, package_version, file_header_BandUP, file_header_BandUP_short, &
file_for_pc_reduced_to_prim_cell, file_for_SC_reduced_to_prim_cell

character(len=30), parameter :: package_version="2.3.0, 2014-06-11"
character(len=127), parameter :: file_header_BandUP="# File created by BandUP - Band Unfolding code for Plane-wave based calculations, &
V"//trim(adjustl(package_version)), &
file_header_BandUP_short="# File created by BandUP, V"//trim(adjustl(package_version)), &
file_for_pc_reduced_to_prim_cell="BandUP_suggestion_of_pc_for_your_reference_unit_cell.POSCAR", &
file_for_SC_reduced_to_prim_cell="BandUP_suggestion_of_smaller_SC_based_on_your_input_SC.POSCAR"
!! Functions and subroutines
CONTAINS


function available_io_unit(min_unit,max_unit) result(unit_num)
! Returns a number unit_num which can be safely used in statements like
! open(file=unit_num)
Expand Down
Loading

0 comments on commit 9c76479

Please sign in to comment.