!-*- mode: F90 -*-! !------------------------------------------------------------! ! Copyright (C) 2026 Wannier Developer Group ! ! ! ! This library is free software; you can redistribute it ! ! and/or modify it under the terms of the GNU Lesser General ! ! Public License as published by the Free Software ! ! Foundation; either version 2.1 of the License, or (at your ! ! option) any later version. ! ! ! ! This library is distributed in the hope that it will be ! ! useful,but WITHOUT ANY WARRANTY; without even the implied ! ! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ! ! PURPOSE. See the GNU Lesser General Public License for ! ! more details. ! ! ! ! You should have received a copy of the GNU Lesser General ! ! Public License along with this library; if not, see ! ! . ! ! ! ! The webpage of the Wannier90 code is ! ! . ! ! ! ! The Wannier90 code is hosted on GitHub ! ! ! !------------------------------------------------------------! ! ! ! w90_postw90_common: routines used throughout postw90 ! ! ! !------------------------------------------------------------! module w90_postw90_common !! This contains the common variables and procedures needed to set up a Wannier !! interpolatation calculation for any physical property use w90_constants, only: dp use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal, & set_error_input, set_error_fatal, set_error_file implicit none private public :: pw90common_fourier_R_to_k public :: pw90common_fourier_R_to_k_new public :: pw90common_fourier_R_to_k_new_second_d public :: pw90common_fourier_R_to_k_new_second_d_TB_conv public :: pw90common_fourier_R_to_k_vec public :: pw90common_fourier_R_to_k_vec_dadb public :: pw90common_fourier_R_to_k_vec_dadb_TB_conv public :: pw90common_get_occ public :: pw90common_kmesh_spacing public :: pw90common_wanint_data_dist public :: pw90common_wanint_get_kpoint_file public :: pw90common_wanint_setup public :: pw90common_wanint_w90_wannier90_readwrite_dist interface pw90common_kmesh_spacing module procedure kmesh_spacing_singleinteger module procedure kmesh_spacing_mesh end interface pw90common_kmesh_spacing contains !================================================! ! PUBLIC PROCEDURES ! Public procedures have names starting with wanint_ !================================================! subroutine pw90common_wanint_setup(num_wann, print_output, real_lattice, mp_grid, ws_region, & ws_distance, effective_model, wigner_seitz, wannier_data, & stdout, seedname, timer, error, comm) !================================================! ! !! Setup data ready for interpolation ! !================================================! use w90_constants, only: dp use w90_io, only: io_stopwatch_start, io_stopwatch_stop use w90_types, only: ws_region_type, ws_distance_type, wannier_data_type, & print_output_type, timer_list_type use w90_comms, only: mpirank, w90_comm_type, comms_bcast use w90_postw90_types, only: wigner_seitz_type use w90_ws_distance, only: ws_translate_dist type(print_output_type), intent(in) :: print_output type(wigner_seitz_type), intent(inout) :: wigner_seitz type(ws_distance_type), intent(inout) :: ws_distance type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(timer_list_type), intent(inout) :: timer type(w90_comm_type), intent(in) :: comm type(w90_error_type), allocatable, intent(out) :: error real(kind=dp), intent(in) :: real_lattice(3, 3) integer, intent(in) :: num_wann integer, intent(in) :: stdout integer, intent(in) :: mp_grid(3) logical, intent(in) :: effective_model character(len=50), intent(in) :: seedname integer :: ierr, ir, file_unit, num_wann_loc logical :: on_root = .false. if (mpirank(comm) == 0) on_root = .true. ! Find nrpts, the number of points in the Wigner-Seitz cell if (effective_model) then if (on_root) then ! nrpts is read from file, together with num_wann open (newunit=file_unit, file=trim(seedname)//'_HH_R.dat', form='formatted', & status='old', err=101) read (file_unit, *) !header read (file_unit, *) num_wann_loc if (num_wann_loc /= num_wann) then call set_error_fatal(error, 'Inconsistent values of num_wann in '//trim(seedname) & //'_HH_R.dat and '//trim(seedname)//'.win', comm) return end if read (file_unit, *) wigner_seitz%nrpts close (file_unit) end if call comms_bcast(wigner_seitz%nrpts, 1, error, comm) if (allocated(error)) return else call wignerseitz(print_output, real_lattice, mp_grid, ws_region, wigner_seitz, stdout, & .true., timer, error, comm) if (allocated(error)) return end if ! Now can allocate several arrays allocate (wigner_seitz%irvec(3, wigner_seitz%nrpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec in pw90common_wanint_setup', comm) return end if wigner_seitz%irvec = 0 allocate (wigner_seitz%crvec(3, wigner_seitz%nrpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating crvec in pw90common_wanint_setup', comm) return end if wigner_seitz%crvec = 0.0_dp allocate (wigner_seitz%ndegen(wigner_seitz%nrpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating ndegen in pw90common_wanint_setup', comm) return end if wigner_seitz%ndegen = 0 ! Also rpt_origin, so that when effective_model=.true it is not ! passed to get_HH_R without being initialized. wigner_seitz%rpt_origin = 0 ! If effective_model, this is done in get_HH_R if (.not. effective_model) then ! Set up the lattice vectors on the Wigner-Seitz supercell ! where the Wannier functions live call wignerseitz(print_output, real_lattice, mp_grid, ws_region, wigner_seitz, stdout, & .false., timer, error, comm) if (allocated(error)) return ! Convert from reduced to Cartesian coordinates do ir = 1, wigner_seitz%nrpts ! Note that 'real_lattice' stores the lattice vectors as *rows* wigner_seitz%crvec(:, ir) = matmul(transpose(real_lattice), wigner_seitz%irvec(:, ir)) end do if (ws_region%use_ws_distance) then call ws_translate_dist(ws_distance, ws_region, num_wann, wannier_data%centres, real_lattice, & mp_grid, wigner_seitz%nrpts, wigner_seitz%irvec, error, comm) if (allocated(error)) return end if call wigner_seitz_opt_setup(ws_distance, ws_region, wigner_seitz, num_wann, real_lattice, & wigner_seitz%nrpts, error, comm) end if return 101 call set_error_file(error, 'Error in pw90common_wanint_setup: problem opening file '// & trim(seedname)//'_HH_R.dat', comm) return end subroutine pw90common_wanint_setup !================================================! subroutine pw90common_wanint_get_kpoint_file(kpoint_dist, error, comm) !================================================! ! !! read kpoints from kpoint.dat and distribute ! !================================================! use w90_constants, only: dp use w90_io, only: io_date, io_time use w90_comms, only: mpirank, mpisize, w90_comm_type, comms_bcast use w90_postw90_types, only: kpoint_dist_type ! arguments type(kpoint_dist_type), intent(inout) :: kpoint_dist type(w90_comm_type), intent(in) :: comm type(w90_error_type), allocatable, intent(out) :: error ! local variables integer :: i, ierr, my_node_id, num_nodes, k_unit !real(kind=dp) :: sum logical :: on_root = .false. real(kind=dp), allocatable :: kt(:, :), wt(:) !temp for read/dist integer :: ik, nk, ir, nkloc, off my_node_id = mpirank(comm) num_nodes = mpisize(comm) if (my_node_id == 0) on_root = .true. if (on_root) then open (newunit=k_unit, file='kpoint.dat', status='old', form='formatted', err=106) read (k_unit, *) kpoint_dist%num_int_kpts end if call comms_bcast(kpoint_dist%num_int_kpts, 1, error, comm) if (allocated(error)) return allocate (kpoint_dist%num_int_kpts_on_node(0:num_nodes - 1)) kpoint_dist%num_int_kpts_on_node(:) = kpoint_dist%num_int_kpts/num_nodes kpoint_dist%max_int_kpts_on_node = kpoint_dist%num_int_kpts & - (num_nodes - 1)*(kpoint_dist%num_int_kpts/num_nodes) kpoint_dist%num_int_kpts_on_node(0) = kpoint_dist%max_int_kpts_on_node ! if(my_node_id < num_int_kpts- num_int_kpts_on_node*num_nodes) num_int_kpts_on_node= num_int_kpts_on_node+1 allocate (kpoint_dist%int_kpts(3, kpoint_dist%max_int_kpts_on_node), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating max_int_kpts_on_node in w90_wannier90_readwrite_read_um', comm) return end if kpoint_dist%int_kpts = 0.0_dp allocate (kpoint_dist%weight(kpoint_dist%max_int_kpts_on_node), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating weight in w90_wannier90_readwrite_read_um', comm) return end if kpoint_dist%weight = 0.0_dp !sum = 0.0_dp nk = kpoint_dist%num_int_kpts allocate (kt(3, nk)) allocate (wt(nk)) if (on_root) then do ik = 1, nk read (k_unit, *) (kt(i, ik), i=1, 3), wt(ik) end do end if call comms_bcast(kt(1, 1), 3*nk, error, comm) if (allocated(error)) return call comms_bcast(wt(1), nk, error, comm) if (allocated(error)) return nkloc = kpoint_dist%num_int_kpts_on_node(my_node_id) off = 1 do ir = 0, my_node_id - 1 off = off + kpoint_dist%num_int_kpts_on_node(ir) end do call dcopy(3*nkloc, kt(1, off), 1, kpoint_dist%int_kpts, 1) call dcopy(nkloc, wt(off), 1, kpoint_dist%weight, 1) deallocate (kt) deallocate (wt) ! if (on_root) then ! do loop_nodes = 1, num_nodes - 1 ! do loop_kpt = 1, kpoint_dist%num_int_kpts_on_node(loop_nodes) ! read (k_unit, *) (kpoint_dist%int_kpts(i, loop_kpt), i=1, 3), kpoint_dist%weight(loop_kpt) ! sum = sum + kpoint_dist%weight(loop_kpt) ! end do ! ! call comms_no_sync_send(kpoint_dist%int_kpts(1, 1), 3*kpoint_dist%num_int_kpts_on_node(loop_nodes), loop_nodes, error, comm) ! if (allocated(error)) return ! call comms_no_sync_send(kpoint_dist%weight(1), kpoint_dist%num_int_kpts_on_node(loop_nodes), loop_nodes, error, comm) ! if (allocated(error)) return ! end do ! do loop_kpt = 1, kpoint_dist%num_int_kpts_on_node(0) ! read (k_unit, *) (kpoint_dist%int_kpts(i, loop_kpt), i=1, 3), kpoint_dist%weight(loop_kpt) ! sum = sum + kpoint_dist%weight(loop_kpt) ! end do !! print*,'rsum',sum ! end if ! ! if (.not. on_root) then ! call comms_no_sync_recv(kpoint_dist%int_kpts(1, 1), 3*kpoint_dist%num_int_kpts_on_node(my_node_id), 0, error, comm) ! if (allocated(error)) return ! call comms_no_sync_recv(kpoint_dist%weight(1), kpoint_dist%num_int_kpts_on_node(my_node_id), 0, error, comm) ! if (allocated(error)) return ! end if return 106 call set_error_file(error, 'Error: Problem opening file kpoint.dat in pw90common_wanint_get_kpoint_file', comm) return end subroutine pw90common_wanint_get_kpoint_file !================================================! subroutine pw90common_wanint_w90_wannier90_readwrite_dist(print_output, ws_region, kmesh_info, kpt_latt, num_kpts, & dis_manifold, w90_system, fermi_energy_list, num_bands, & num_wann, eigval, mp_grid, real_lattice, & pw90_calculation, scissors_shift, effective_model, & pw90_spin, pw90_band_deriv_degen, pw90_kpath, & pw90_kslice, pw90_dos, pw90_berry, pw90_spin_hall, & pw90_gyrotropic, pw90_geninterp, pw90_boltzwann, & eig_found, error, comm) !================================================! ! !! distribute the parameters across processors !! NOTE: we only send the ones postw90 uses, not all in w90 ! !================================================! use w90_constants, only: dp use w90_io, only: io_date, io_time use w90_comms, only: mpirank, w90_comm_type, comms_bcast use w90_types use w90_postw90_types, only: pw90_calculation_type, pw90_spin_mod_type, & pw90_band_deriv_degen_type, pw90_kpath_mod_type, pw90_kslice_mod_type, pw90_dos_mod_type, & pw90_berry_mod_type, pw90_spin_hall_type, pw90_gyrotropic_type, pw90_geninterp_mod_type, & pw90_boltzwann_type type(print_output_type), intent(inout) :: print_output type(ws_region_type), intent(inout) :: ws_region type(w90_system_type), intent(inout) :: w90_system type(kmesh_info_type), intent(inout) :: kmesh_info type(dis_manifold_type), intent(inout) :: dis_manifold type(pw90_calculation_type), intent(inout) :: pw90_calculation type(pw90_spin_mod_type), intent(inout) :: pw90_spin type(pw90_band_deriv_degen_type), intent(inout) :: pw90_band_deriv_degen type(pw90_kpath_mod_type), intent(inout) :: pw90_kpath type(pw90_kslice_mod_type), intent(inout) :: pw90_kslice type(pw90_dos_mod_type), intent(inout) :: pw90_dos type(pw90_berry_mod_type), intent(inout) :: pw90_berry type(pw90_spin_hall_type), intent(inout) :: pw90_spin_hall type(pw90_gyrotropic_type), intent(inout) :: pw90_gyrotropic type(pw90_geninterp_mod_type), intent(inout) :: pw90_geninterp type(pw90_boltzwann_type), intent(inout) :: pw90_boltzwann type(w90_comm_type), intent(in) :: comm type(w90_error_type), allocatable, intent(out) :: error real(kind=dp), allocatable, intent(inout) :: kpt_latt(:, :) real(kind=dp), allocatable, intent(inout) :: fermi_energy_list(:) real(kind=dp), pointer, intent(inout) :: eigval(:, :) real(kind=dp), intent(inout) :: real_lattice(3, 3) real(kind=dp), intent(inout) :: scissors_shift integer, intent(inout) :: num_kpts integer, intent(inout) :: num_bands integer, intent(inout) :: num_wann integer, intent(inout) :: mp_grid(3) logical, intent(inout) :: eig_found logical, intent(inout) :: effective_model integer :: ierr integer :: iprintroot integer :: fermi_n, kdotp_nbands logical :: on_root = .false. if (mpirank(comm) == 0) on_root = .true. call comms_bcast(effective_model, 1, error, comm) if (allocated(error)) return call comms_bcast(eig_found, 1, error, comm) if (allocated(error)) return if (.not. effective_model) then call comms_bcast(mp_grid(1), 3, error, comm) if (allocated(error)) return call comms_bcast(num_kpts, 1, error, comm) if (allocated(error)) return call comms_bcast(num_bands, 1, error, comm) if (allocated(error)) return end if call comms_bcast(num_wann, 1, error, comm) if (allocated(error)) return call comms_bcast(print_output%timing_level, 1, error, comm) if (allocated(error)) return !______________________________________ ! not so pretty solution to setting iprint to zero on non-root processes iprintroot = print_output%iprint print_output%iprint = 0 call comms_bcast(print_output%iprint, 1, error, comm) if (allocated(error)) return if (on_root) print_output%iprint = iprintroot !______________________________________ call comms_bcast(ws_region%ws_distance_tol, 1, error, comm) if (allocated(error)) return call comms_bcast(ws_region%ws_search_size(1), 3, error, comm) if (allocated(error)) return ! call comms_bcast(num_atoms,1) ! Ivo: not used in postw90, right? ! call comms_bcast(num_species,1) ! Ivo: not used in postw90, right? call comms_bcast(real_lattice(1, 1), 9, error, comm) if (allocated(error)) return !call comms_bcast(recip_lattice(1, 1), 9, error, comm) !call comms_bcast(real_metric(1, 1), 9) !call comms_bcast(recip_metric(1, 1), 9) !call comms_bcast(cell_volume, 1, error, comm) call comms_bcast(pw90_dos%energy_step, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%smearing%use_adaptive, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%smearing%type_index, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%kmesh%spacing, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%kmesh%mesh(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%smearing%adaptive_max_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%smearing%fixed_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%smearing%adaptive_prefactor, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%num_project, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_calculation%berry, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%task, len(pw90_berry%task), error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kmesh%spacing, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kmesh%mesh(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%curv_adpt_kmesh, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%curv_adpt_kmesh_thresh, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%curv_unit, len(pw90_berry%curv_unit), error, comm) if (allocated(error)) return ! Tsirkin call comms_bcast(pw90_calculation%gyrotropic, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%task, len(pw90_gyrotropic%task), error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%kmesh%spacing, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%kmesh%mesh(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%eigval_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%nfreq, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%degen_thresh, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%num_bands, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%box(1, 1), 9, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%box_corner(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%smearing%use_adaptive, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%smearing%fixed_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%smearing%type_index, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%smearing%max_arg, 1, error, comm) if (allocated(error)) return call comms_bcast(w90_system%spinors, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%freq_scan, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%alpha, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%beta, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%gamma, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%bandshift, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%bandshift_firstband, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%bandshift_energyshift, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin_hall%method, len(pw90_spin_hall%method), error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_smearing%use_adaptive, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_smearing%adaptive_prefactor, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_smearing%adaptive_max_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_smearing%fixed_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_smearing%type_index, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_eigval_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_nfreq, 1, error, comm) if (allocated(error)) return fermi_n = 0 if (on_root) then if (allocated(fermi_energy_list)) fermi_n = size(fermi_energy_list) end if call comms_bcast(fermi_n, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%energy_min, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%energy_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin%kmesh%spacing, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin%kmesh%mesh(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%wanint_kpoint_file, 1, error, comm) if (allocated(error)) return call comms_bcast(dis_manifold%win_min, 1, error, comm) if (allocated(error)) return call comms_bcast(dis_manifold%win_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%sc_eta, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%sc_w_thr, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%sc_phase_conv, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%sc_use_eta_corr, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%tetrahedron_method, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%tetrahedron_higher_correction, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%tetrahedron_cutoff, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%tetrahedron_avoid_degeneracy, 1, error, comm) if (allocated(error)) return ! ---------------------------------------------- ! ! New input variables in development ! !call comms_bcast(print_output%devel_flag, len(print_output%devel_flag), error, comm) call comms_bcast(pw90_calculation%spin_moment, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin%axis_polar, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_spin%axis_azimuth, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_calculation%spin_decomp, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_band_deriv_degen%use_degen_pert, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_band_deriv_degen%degen_thr, 1, error, comm) if (allocated(error)) return call comms_bcast(w90_system%num_valence_bands, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_calculation%dos, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%task, len(pw90_dos%task), error, comm) if (allocated(error)) return call comms_bcast(pw90_calculation%kpath, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_kpath%task, len(pw90_kpath%task), error, comm) if (allocated(error)) return call comms_bcast(pw90_kpath%bands_colour, len(pw90_kpath%bands_colour), error, comm) if (allocated(error)) return call comms_bcast(pw90_calculation%kslice, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%task, len(pw90_kslice%task), error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%corner(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%b1(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%b2(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%kmesh2d(1), 2, error, comm) if (allocated(error)) return call comms_bcast(pw90_kslice%fermi_lines_colour, len(pw90_kslice%fermi_lines_colour), error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%transl_inv, 1, error, comm) if (allocated(error)) return call comms_bcast(w90_system%num_elec_per_state, 1, error, comm) if (allocated(error)) return call comms_bcast(scissors_shift, 1, error, comm) if (allocated(error)) return ! Do these have to be broadcasted? (Plots done on root node only) ! ! call comms_bcast(bands_num_points,1) ! call comms_bcast(bands_num_spec_points,1) ! if(allocated(bands_spec_points)) & ! call comms_bcast(bands_spec_points(1,1),3*bands_num_spec_points) ! if(allocated(bands_label)) & ! call comms_bcast(bands_label(:),len(bands_label(1))*bands_num_spec_points) ! ---------------------------------------------- call comms_bcast(pw90_calculation%geninterp, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_geninterp%alsofirstder, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_geninterp%single_file, 1, error, comm) if (allocated(error)) return ! BoltzWann variables call comms_bcast(pw90_calculation%boltzwann, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%calc_also_dos, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dir_num_2d, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_energy_step, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_energy_min, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_energy_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_smearing%use_adaptive, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_smearing%fixed_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_smearing%adaptive_prefactor, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_smearing%adaptive_max_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%mu_min, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%mu_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%mu_step, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%temp_min, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%temp_max, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%temp_step, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%kmesh%spacing, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%kmesh%mesh(1), 3, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%tdf_energy_step, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%relax_time, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%tdf_smearing%use_adaptive, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%tdf_smearing%fixed_width, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%tdf_smearing%type_index, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%dos_smearing%type_index, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%bandshift, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%bandshift_firstband, 1, error, comm) if (allocated(error)) return call comms_bcast(pw90_boltzwann%bandshift_energyshift, 1, error, comm) if (allocated(error)) return call comms_bcast(ws_region%use_ws_distance, 1, error, comm) if (allocated(error)) return ! These variables are different from the ones above in that they are ! allocatable, and in w90_wannier90_readwrite_read they were allocated on the root node only if (.not. on_root) then allocate (fermi_energy_list(fermi_n), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating fermi_energy_read in postw90_w90_wannier90_readwrite_dist', comm) return end if allocate (pw90_berry%kubo_freq_list(pw90_berry%kubo_nfreq), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating kubo_freq_list in postw90_w90_wannier90_readwrite_dist', comm) return end if allocate (pw90_gyrotropic%band_list(pw90_gyrotropic%num_bands), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating gyrotropic_band_list in postw90_w90_wannier90_readwrite_dist', comm) return end if allocate (pw90_gyrotropic%freq_list(pw90_gyrotropic%nfreq), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating gyrotropic_freq_list in postw90_w90_wannier90_readwrite_dist', comm) return end if allocate (pw90_dos%project(pw90_dos%num_project), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating dos_project in postw90_w90_wannier90_readwrite_dist', comm) return end if if (.not. effective_model) then if (eig_found) then allocate (eigval(num_bands, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating eigval in postw90_w90_wannier90_readwrite_dist', comm) return end if end if allocate (kpt_latt(3, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating kpt_latt in postw90_w90_wannier90_readwrite_dist', comm) return end if end if end if kdotp_nbands = 0 if (on_root) then if (allocated(pw90_berry%kdotp_bands)) kdotp_nbands = size(pw90_berry%kdotp_bands) end if call comms_bcast(kdotp_nbands, 1, error, comm) if (kdotp_nbands > 0) then if (.not. on_root) then allocate (pw90_berry%kdotp_bands(kdotp_nbands), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating kdotp_bands in postw90_param_dist', comm) return end if end if call comms_bcast(pw90_berry%kdotp_bands(1), kdotp_nbands, error, comm) end if if (fermi_n > 0) then call comms_bcast(fermi_energy_list(1), fermi_n, error, comm) if (allocated(error)) return end if call comms_bcast(pw90_gyrotropic%freq_list(1), pw90_gyrotropic%nfreq, error, comm) if (allocated(error)) return call comms_bcast(pw90_gyrotropic%band_list(1), pw90_gyrotropic%num_bands, error, comm) if (allocated(error)) return call comms_bcast(pw90_berry%kubo_freq_list(1), pw90_berry%kubo_nfreq, error, comm) if (allocated(error)) return call comms_bcast(pw90_dos%project(1), pw90_dos%num_project, error, comm) if (allocated(error)) return if (.not. effective_model) then if (eig_found) then call comms_bcast(eigval(1, 1), num_bands*num_kpts, error, comm) if (allocated(error)) return end if call comms_bcast(kpt_latt(1, 1), 3*num_kpts, error, comm) if (allocated(error)) return end if ! kmesh: only nntot,wb, and bk are needed to evaluate the WF matrix ! elements of the position operator in reciprocal space. For the ! extra matrix elements entering the orbital magnetization, also ! need nnlist. In principle could only broadcast those four variables if (.not. effective_model) then call comms_bcast(kmesh_info%nnh, 1, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%nntot, 1, error, comm) if (allocated(error)) return if (.not. on_root) then allocate (kmesh_info%nnlist(num_kpts, kmesh_info%nntot), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating nnlist in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if allocate (kmesh_info%neigh(num_kpts, kmesh_info%nntot/2), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating neigh in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if allocate (kmesh_info%nncell(3, num_kpts, kmesh_info%nntot), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating nncell in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if allocate (kmesh_info%wb(kmesh_info%nntot), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating wb in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if allocate (kmesh_info%bka(3, kmesh_info%nntot/2), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating bka in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if allocate (kmesh_info%bk(3, kmesh_info%nntot, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating bk in pw90common_wanint_w90_wannier90_readwrite_dist', comm) return end if end if call comms_bcast(kmesh_info%nnlist(1, 1), num_kpts*kmesh_info%nntot, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%neigh(1, 1), num_kpts*kmesh_info%nntot/2, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%nncell(1, 1, 1), 3*num_kpts*kmesh_info%nntot, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%wb(1), kmesh_info%nntot, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%bka(1, 1), 3*kmesh_info%nntot/2, error, comm) if (allocated(error)) return call comms_bcast(kmesh_info%bk(1, 1, 1), 3*kmesh_info%nntot*num_kpts, error, comm) if (allocated(error)) return end if end subroutine pw90common_wanint_w90_wannier90_readwrite_dist !================================================! subroutine pw90common_wanint_data_dist(num_wann, num_kpts, num_bands, u_matrix_opt, u_matrix, & dis_manifold, wannier_data, scissors_shift, v_matrix, & num_valence_bands, have_disentangled, error, comm) !================================================! ! !! Distribute the um and chk files ! !================================================! use w90_constants, only: dp, cmplx_0 use w90_io, only: io_date, io_time use w90_types, only: dis_manifold_type, wannier_data_type use w90_comms, only: w90_comm_type, mpirank, comms_bcast implicit none type(dis_manifold_type), intent(inout) :: dis_manifold type(wannier_data_type), intent(inout) :: wannier_data type(w90_comm_type), intent(in) :: comm type(w90_error_type), allocatable, intent(out) :: error integer, intent(in) :: num_valence_bands integer, intent(in) :: num_wann, num_kpts, num_bands real(kind=dp), intent(in) :: scissors_shift complex(kind=dp), allocatable :: v_matrix(:, :, :) complex(kind=dp), allocatable, intent(inout) :: u_matrix_opt(:, :, :), u_matrix(:, :, :) logical, intent(inout) :: have_disentangled integer :: ierr, loop_kpt, m, i, j logical :: on_root = .false. if (mpirank(comm) == 0) on_root = .true. ! if (.not. on_root) then ! ! wannier_centres is allocated in w90_wannier90_readwrite_read, so only on root node ! ! It is then read in w90_wannier90_readwrite_read_chpkt ! ! Therefore, now we need to allocate it on all nodes, and then broadcast it ! allocate (wannier_data%centres(3, num_wann), stat=ierr) ! if (ierr /= 0) then ! call set_error_alloc(error, 'Error allocating wannier_centres in pw90common_wanint_data_dist', comm) ! return ! endif ! end if ! call comms_bcast(wannier_data%centres(1, 1), 3*num_wann, error, comm) ! if (allocated(error)) return ! ------------------- ! Ivo: added 8april11 ! ------------------- ! ! Calculate the matrix that describes the combined effect of ! disentanglement and maximal localization. This is the combination ! that is most often needed for interpolation purposes ! ! Allocate on all nodes allocate (v_matrix(num_bands, num_wann, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating v_matrix in pw90common_wanint_data_dist', comm) return end if ! u_matrix and u_matrix_opt are stored on root only if (on_root) then if (.not. have_disentangled) then v_matrix(1:num_wann, :, :) = u_matrix(1:num_wann, :, :) else v_matrix = cmplx_0 do loop_kpt = 1, num_kpts do j = 1, num_wann do m = 1, dis_manifold%ndimwin(loop_kpt) do i = 1, num_wann v_matrix(m, j, loop_kpt) = v_matrix(m, j, loop_kpt) & + u_matrix_opt(m, i, loop_kpt)*u_matrix(i, j, loop_kpt) end do end do end do end do end if if (allocated(u_matrix_opt)) deallocate (u_matrix_opt) if (.not. (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp)) then if (allocated(u_matrix)) deallocate (u_matrix) end if end if call comms_bcast(v_matrix(1, 1, 1), num_bands*num_wann*num_kpts, error, comm) if (allocated(error)) return if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then if (.not. on_root .and. .not. allocated(u_matrix)) then allocate (u_matrix(num_wann, num_wann, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating u_matrix in pw90common_wanint_data_dist', comm) return end if end if call comms_bcast(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts, error, comm) end if call comms_bcast(have_disentangled, 1, error, comm) if (allocated(error)) return if (have_disentangled) then if (.not. on_root) then if (.not. allocated(dis_manifold%lwindow)) then allocate (dis_manifold%lwindow(num_bands, num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating lwindow in pw90common_wanint_data_dist', comm) return end if end if if (.not. allocated(dis_manifold%ndimwin)) then allocate (dis_manifold%ndimwin(num_kpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating ndimwin in pw90common_wanint_data_dist', comm) return end if end if end if ! call comms_bcast(u_matrix_opt(1,1,1),num_bands*num_wann*num_kpts) call comms_bcast(dis_manifold%lwindow(1, 1), num_bands*num_kpts, error, comm) if (allocated(error)) return call comms_bcast(dis_manifold%ndimwin(1), num_kpts, error, comm) if (allocated(error)) return end if end subroutine pw90common_wanint_data_dist !================================================ subroutine pw90common_get_occ(ef, eig, occ, num_wann) !================================================! ! !! Compute the electronic occupancy ! !================================================! use w90_constants, only: dp ! arguments integer, intent(in) :: num_wann real(kind=dp), intent(in) :: eig(num_wann) !! Eigenvalues real(kind=dp), intent(in) :: ef !! Fermi level real(kind=dp), intent(out) :: occ(num_wann) !! Occupancy of states ! local variables integer :: i ! State occupancies ! if(smear_temp < eps7) then ! ! Use a step function occupancy (T=0) ! occ(:) = 0.0_dp do i = 1, num_wann if (eig(i) < ef) occ(i) = 1.0_dp end do ! else ! ! Use a Fermi-Dirac occupancy (T=smear_temp, in Kelvin) ! ! k_B.T in electron-volts ! ! kt=k_B_SI*smear_temp/elem_charge_SI ! do i=1,num_wann ! occ(i)=1.0_dp/(exp((eig(i)-ef)/kt)+1.0_dp) ! end do ! end if end subroutine pw90common_get_occ !================================================ function kmesh_spacing_singleinteger(num_points, recip_lattice) !================================================ !! Set up the value of the interpolation mesh spacing, needed for !! adaptive smearing [see Eqs. (34-35) YWVS07]. Choose it as the largest of !! the three Delta_k's for each of the primitive translations b1, b2, and b3 !================================================ integer, intent(in) :: num_points real(kind=dp), intent(in) :: recip_lattice(3, 3) real(kind=dp) :: kmesh_spacing_singleinteger integer :: i real(kind=dp) :: Delta_k_i(3) ! NOTE: The vectors b_i are stored as *rows* in recip_lattice (argh!). ! Hence I believe Jonathan's original code confused rows with columns ! when computing Delta_k, which he called 'rspace' ! (See my e-mail of 20Sept07) ! do i = 1, 3 Delta_k_i(i) = sqrt(dot_product(recip_lattice(i, :), recip_lattice(i, :)))/num_points end do kmesh_spacing_singleinteger = maxval(Delta_k_i) end function kmesh_spacing_singleinteger function kmesh_spacing_mesh(mesh, recip_lattice) !! Same as kmesh_spacing_singleinteger, but for a kmesh with three !! different mesh samplings along the three directions integer, dimension(3), intent(in) :: mesh real(kind=dp), intent(in) :: recip_lattice(3, 3) real(kind=dp) :: kmesh_spacing_mesh integer :: i real(kind=dp) :: Delta_k_i(3) do i = 1, 3 Delta_k_i(i) = sqrt(dot_product(recip_lattice(i, :), recip_lattice(i, :)))/mesh(i) end do kmesh_spacing_mesh = maxval(Delta_k_i) end function kmesh_spacing_mesh !================================================! subroutine pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, OO, & OO_R, kpt, real_lattice, mp_grid, alpha, num_wann, & error, comm) !================================================! ! !! For alpha=0: !! O_ij(R) --> O_ij(k) = sum_R e^{+ik.R}*O_ij(R) !! !! For alpha=1,2,3: !! sum_R [cmplx_i*R_alpha*e^{+ik.R}*O_ij(R)] !! where R_alpha is a Cartesian component of R !! ***REMOVE EVENTUALLY*** (replace with pw90common_fourier_R_to_k_new) ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: wannier_data_type, ws_region_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(wigner_seitz_type), intent(in) :: wigner_seitz type(ws_distance_type), intent(inout) :: ws_distance type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: num_wann integer, intent(in) :: mp_grid(3) integer, intent(in) :: alpha real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :) complex(kind=dp), intent(out) :: OO(:, :) ! local variables integer :: ir, i, j, ideg real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac OO(:, :) = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: rdotk = twopi*dot_product(kpt(:), wigner_seitz%irvec_pw90(:, ir)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (alpha == 0) then OO(:, :) = OO(:, :) + phase_fac*OO_R(:, :, ir) elseif (alpha == 1 .or. alpha == 2 .or. alpha == 3) then OO(:, :) = OO(:, :) + & cmplx_i*wigner_seitz%crvec_pw90(alpha, ir)*phase_fac*OO_R(:, :, ir) else stop 'wrong value of alpha in pw90common_fourier_R_to_k' end if end do end subroutine pw90common_fourier_R_to_k !================================================! subroutine pw90common_fourier_R_to_k_new(ws_region, wannier_data, ws_distance, wigner_seitz, & OO_R, kpt, real_lattice, mp_grid, num_wann, error, & comm, OO, OO_dx, OO_dy, OO_dz) !================================================! ! !! For OO: !! $$O_{ij}(k) = \sum_R e^{+ik.R}.O_{ij}(R)$$ !! For $$OO_{dx,dy,dz}$$: !! $$\sum_R [i.R_{dx,dy,dz}.e^{+ik.R}.O_{ij}(R)]$$ !! where R_{x,y,z} are the Cartesian components of R ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(wigner_seitz_type), intent(in) :: wigner_seitz type(ws_distance_type), intent(inout) :: ws_distance type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: num_wann integer, intent(in) :: mp_grid(3) real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :) complex(kind=dp), optional, intent(out) :: OO(:, :) complex(kind=dp), optional, intent(out) :: OO_dx(:, :) complex(kind=dp), optional, intent(out) :: OO_dy(:, :) complex(kind=dp), optional, intent(out) :: OO_dz(:, :) ! local variables integer :: ir, i, j, ideg real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac if (present(OO)) OO = cmplx_0 if (present(OO_dx)) OO_dx = cmplx_0 if (present(OO_dy)) OO_dy = cmplx_0 if (present(OO_dz)) OO_dz = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: rdotk = twopi*dot_product(kpt(:), wigner_seitz%irvec_pw90(:, ir)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO)) OO(:, :) = OO(:, :) + phase_fac*OO_R(:, :, ir) if (present(OO_dx)) OO_dx(:, :) = OO_dx(:, :) + & cmplx_i*wigner_seitz%crvec_pw90(1, ir)*phase_fac*OO_R(:, :, ir) if (present(OO_dy)) OO_dy(:, :) = OO_dy(:, :) + & cmplx_i*wigner_seitz%crvec_pw90(2, ir)*phase_fac*OO_R(:, :, ir) if (present(OO_dz)) OO_dz(:, :) = OO_dz(:, :) + & cmplx_i*wigner_seitz%crvec_pw90(3, ir)*phase_fac*OO_R(:, :, ir) end do end subroutine pw90common_fourier_R_to_k_new !================================================! subroutine pw90common_fourier_R_to_k_new_second_d(kpt, OO_R, num_wann, ws_region, wannier_data, & real_lattice, mp_grid, ws_distance, & wigner_seitz, error, comm, OO, OO_da, OO_dadb) !================================================! ! !! For OO: !! $$O_{ij}(k) = \sum_R e^{+ik.R}.O_{ij}(R)$$ !! For $$OO_{dx,dy,dz}$$: !! $$\sum_R [i.R_{dx,dy,dz}.e^{+ik.R}.O_{ij}(R)]$$ !! where R_{x,y,z} are the Cartesian components of R !! For $$OO_{dx1,dy1,dz1;dx2,dy2,dz2}$$: !! $$-\sum_R [R_{dx1,dy1,dz1}.R_{dx2,dy2,dz2}.e^{+ik.R}.O_{ij}(R)]$$ !! where R_{xi,yi,zi} are the Cartesian components of R ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(wigner_seitz_type), intent(in) :: wigner_seitz type(ws_distance_type), intent(inout) :: ws_distance type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: mp_grid(3) integer, intent(in) :: num_wann real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :) complex(kind=dp), optional, intent(out) :: OO(:, :) complex(kind=dp), optional, intent(out) :: OO_da(:, :, :) complex(kind=dp), optional, intent(out) :: OO_dadb(:, :, :, :) ! local variables integer :: ir, i, j, ideg, a, b real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac if (present(OO)) OO = cmplx_0 if (present(OO_da)) OO_da = cmplx_0 if (present(OO_dadb)) OO_dadb = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: rdotk = twopi*dot_product(kpt(:), wigner_seitz%irvec_pw90(:, ir)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO)) OO(:, :) = OO(:, :) + phase_fac*OO_R(:, :, ir) if (present(OO_da)) then do a = 1, 3 OO_da(:, :, a) = OO_da(:, :, a) + cmplx_i*wigner_seitz%crvec_pw90(a, ir)*phase_fac & *OO_R(:, :, ir) end do end if if (present(OO_dadb)) then do a = 1, 3 do b = 1, 3 OO_dadb(:, :, a, b) = OO_dadb(:, :, a, b) - & wigner_seitz%crvec_pw90(a, ir)*wigner_seitz%crvec_pw90(b, ir)*phase_fac & *OO_R(:, :, ir) end do end do end if end do end subroutine pw90common_fourier_R_to_k_new_second_d !================================================! subroutine pw90common_fourier_R_to_k_new_second_d_TB_conv(kpt, OO_R, oo_a_R, num_wann, & ws_region, wannier_data, real_lattice, & mp_grid, ws_distance, wigner_seitz, & error, comm, OO, OO_da, OO_dadb) !================================================! ! modified version of pw90common_fourier_R_to_k_new_second_d, includes wannier centres in ! the exponential inside the sum (so called TB convention) ! !! For OO: !! $$O_{ij}(k) = \sum_R e^{+ik.(R+tau_ij)}.O_{ij}(R)$$ !! For $$OO_{dx,dy,dz}$$: !! $$\sum_R [i.(R+tau_ij)_{dx,dy,dz}.e^{+ik.(R+tau_ij)}.O_{ij}(R)]$$ !! where R_{x,y,z} are the Cartesian components of R !! For $$OO_{dx1,dy1,dz1;dx2,dy2,dz2}$$: !! $$-\sum_R [(R+tau_ij)_{dx1,dy1,dz1}.(R+tau_ij)_{dx2,dy2,dz2}.e^{+ik.(R+tau_ij)}.O_{ij}(R)]$$ !! where {xi,yi,zi} denote the Cartesian components and ! tau_ij = tau_j - tau_i, being tau_i=<0i|r|0i> the individual wannier centres !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_utility, only: utility_cart_to_frac, utility_inverse_mat use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(wigner_seitz_type), intent(in) :: wigner_seitz type(ws_distance_type), intent(inout) :: ws_distance type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: mp_grid(3) integer, intent(in) :: num_wann real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: oo_a_R(:, :, :, :) complex(kind=dp), intent(in) :: OO_R(:, :, :) complex(kind=dp), optional, intent(out) :: OO(:, :) complex(kind=dp), optional, intent(out) :: OO_da(:, :, :) complex(kind=dp), optional, intent(out) :: OO_dadb(:, :, :, :) ! local variables real(kind=dp) :: inv_lattice(3, 3) integer :: ir, i, j, ideg, a, b real(kind=dp) :: rdotk real(kind=dp) :: wannier_centres_frac(3, num_wann) real(kind=dp) :: r_sum(3) complex(kind=dp) :: phase_fac r_sum = 0.d0 ! rotate wannier centres from cartesian to fractional coordinates wannier_centres_frac(:, :) = 0.d0 call utility_inverse_mat(real_lattice, inv_lattice) do ir = 1, num_wann call utility_cart_to_frac(wigner_seitz%wannier_centres_from_AA_R(:, ir), & wannier_centres_frac(:, ir), inv_lattice) end do if (present(OO)) OO = cmplx_0 if (present(OO_da)) OO_da = cmplx_0 if (present(OO_dadb)) OO_dadb = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: do j = 1, num_wann do i = 1, num_wann r_sum(:) = real(wigner_seitz%irvec_pw90(:, ir)) & + wannier_centres_frac(:, j) - wannier_centres_frac(:, i) rdotk = twopi*dot_product(kpt(:), r_sum(:)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO)) OO(i, j) = OO(i, j) + phase_fac*OO_R(i, j, ir) if (present(OO_da)) then do a = 1, 3 OO_da(i, j, a) = OO_da(i, j, a) + cmplx_i* & (wigner_seitz%crvec_pw90(a, ir) + wigner_seitz%wannier_centres_from_AA_R(a, j) - & wigner_seitz%wannier_centres_from_AA_R(a, i))*phase_fac*OO_R(i, j, ir) end do end if if (present(OO_dadb)) then do a = 1, 3 do b = 1, 3 OO_dadb(i, j, a, b) = & OO_dadb(i, j, a, b) - & (wigner_seitz%crvec_pw90(a, ir) + wigner_seitz%wannier_centres_from_AA_R(a, j) - & wigner_seitz%wannier_centres_from_AA_R(a, i))*(wigner_seitz%crvec_pw90(b, ir) + & wigner_seitz%wannier_centres_from_AA_R(b, j) - & wigner_seitz%wannier_centres_from_AA_R(b, i))* & phase_fac*OO_R(i, j, ir) end do end do end if end do end do end do end subroutine pw90common_fourier_R_to_k_new_second_d_TB_conv !================================================! subroutine pw90common_fourier_R_to_k_vec(ws_region, wannier_data, ws_distance, wigner_seitz, & OO_R, kpt, real_lattice, mp_grid, num_wann, error, & comm, OO_true, OO_pseudo) !================================================! ! !! For OO_true (true vector): !! $${\vec O}_{ij}(k) = \sum_R e^{+ik.R} {\vec O}_{ij}(R)$$ ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(ws_distance_type), intent(inout) :: ws_distance type(wigner_seitz_type), intent(in) :: wigner_seitz type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: num_wann integer, intent(in) :: mp_grid(3) real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :, :) complex(kind=dp), optional, intent(out) :: OO_true(:, :, :) complex(kind=dp), optional, intent(out) :: OO_pseudo(:, :, :) ! local variables integer :: ir, i, j, ideg real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac if (present(OO_true)) OO_true = cmplx_0 if (present(OO_pseudo)) OO_pseudo = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: rdotk = twopi*dot_product(kpt(:), wigner_seitz%irvec_pw90(:, ir)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO_true)) then OO_true(:, :, 1) = OO_true(:, :, 1) + phase_fac*OO_R(:, :, ir, 1) OO_true(:, :, 2) = OO_true(:, :, 2) + phase_fac*OO_R(:, :, ir, 2) OO_true(:, :, 3) = OO_true(:, :, 3) + phase_fac*OO_R(:, :, ir, 3) end if if (present(OO_pseudo)) then OO_pseudo(:, :, 1) = OO_pseudo(:, :, 1) & + cmplx_i*wigner_seitz%crvec_pw90(2, ir)*phase_fac*OO_R(:, :, ir, 3) & - cmplx_i*wigner_seitz%crvec_pw90(3, ir)*phase_fac*OO_R(:, :, ir, 2) OO_pseudo(:, :, 2) = OO_pseudo(:, :, 2) & + cmplx_i*wigner_seitz%crvec_pw90(3, ir)*phase_fac*OO_R(:, :, ir, 1) & - cmplx_i*wigner_seitz%crvec_pw90(1, ir)*phase_fac*OO_R(:, :, ir, 3) OO_pseudo(:, :, 3) = OO_pseudo(:, :, 3) & + cmplx_i*wigner_seitz%crvec_pw90(1, ir)*phase_fac*OO_R(:, :, ir, 2) & - cmplx_i*wigner_seitz%crvec_pw90(2, ir)*phase_fac*OO_R(:, :, ir, 1) end if end do end subroutine pw90common_fourier_R_to_k_vec !================================================! subroutine pw90common_fourier_R_to_k_vec_dadb(ws_region, wannier_data, ws_distance, & wigner_seitz, OO_R, kpt, real_lattice, mp_grid, & num_wann, error, comm, OO_da, OO_dadb) !================================================! ! !! For $$OO_{ij;dx,dy,dz}$$: !! $$O_{ij;dx,dy,dz}(k) = \sum_R e^{+ik.R} O_{ij;dx,dy,dz}(R)$$ !! For $$OO_{ij;dx1,dy1,dz1;dx2,dy2,dz2}$$: !! $$O_{ij;dx1,dy1,dz1;dx2,dy2,dz2}(k) = \sum_R e^{+ik.R} i.R_{dx2,dy2,dz2} !! .O_{ij;dx1,dy1,dz1}(R)$$ ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(ws_distance_type), intent(inout) :: ws_distance type(wigner_seitz_type), intent(in) :: wigner_seitz type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: num_wann integer, intent(in) :: mp_grid(3) real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :, :) complex(kind=dp), optional, intent(out) :: OO_da(:, :, :) complex(kind=dp), optional, intent(out) :: OO_dadb(:, :, :, :) ! local variables integer :: ir, i, j, ideg, a, b real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac if (present(OO_da)) OO_da = cmplx_0 if (present(OO_dadb)) OO_dadb = cmplx_0 do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: rdotk = twopi*dot_product(kpt(:), wigner_seitz%irvec_pw90(:, ir)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO_da)) then OO_da(:, :, 1) = OO_da(:, :, 1) + phase_fac*OO_R(:, :, ir, 1) OO_da(:, :, 2) = OO_da(:, :, 2) + phase_fac*OO_R(:, :, ir, 2) OO_da(:, :, 3) = OO_da(:, :, 3) + phase_fac*OO_R(:, :, ir, 3) end if if (present(OO_dadb)) then do a = 1, 3 do b = 1, 3 OO_dadb(:, :, a, b) = OO_dadb(:, :, a, b) & + cmplx_i*wigner_seitz%crvec_pw90(b, ir)*phase_fac*OO_R(:, :, ir, a) end do end do end if end do end subroutine pw90common_fourier_R_to_k_vec_dadb !================================================! subroutine pw90common_fourier_R_to_k_vec_dadb_TB_conv(ws_region, wannier_data, ws_distance, & wigner_seitz, OO_R, kpt, real_lattice, & mp_grid, num_wann, error, comm, OO_da, & OO_dadb) !================================================! ! ! modified version of pw90common_fourier_R_to_k_vec_dadb, includes wannier centres in ! the exponential inside the sum (so called TB convention) ! !! For $$OO_{ij;dx,dy,dz}$$: !! $$O_{ij;dx,dy,dz}(k) = \sum_R e^{+ik.(R+tau_ij)} O_{ij;dx,dy,dz}(R)$$ !! For $$OO_{ij;dx1,dy1,dz1;dx2,dy2,dz2}$$: !! $$O_{ij;dx1,dy1,dz1;dx2,dy2,dz2}(k) = \sum_R e^{+ik.(R+tau_ij)} i.(R+tau_ij)_{dx2,dy2,dz2} !! .O_{ij;dx1,dy1,dz1}(R)$$ ! with tau_ij = tau_j - tau_i, being tau_i=<0i|r|0i> the individual wannier centres ! !================================================! use w90_constants, only: dp, cmplx_0, cmplx_i, twopi use w90_types, only: ws_region_type, wannier_data_type, ws_distance_type use w90_ws_distance, only: ws_translate_dist use w90_utility, only: utility_cart_to_frac, utility_inverse_mat use w90_postw90_types, only: wigner_seitz_type use w90_comms, only: w90_comm_type implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(wannier_data_type), intent(in) :: wannier_data type(ws_distance_type), intent(inout) :: ws_distance type(wigner_seitz_type), intent(in) :: wigner_seitz type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm integer, intent(in) :: num_wann integer, intent(in) :: mp_grid(3) real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3) complex(kind=dp), intent(in) :: OO_R(:, :, :, :) complex(kind=dp), optional, intent(out) :: OO_da(:, :, :) complex(kind=dp), optional, intent(out) :: OO_dadb(:, :, :, :) ! local variables real(kind=dp) :: inv_lattice(3, 3) integer :: ir, i, j, ideg, a, b real(kind=dp) :: rdotk complex(kind=dp) :: phase_fac real(kind=dp) :: wannier_centres_frac(3, num_wann) real(kind=dp) :: r_sum(3) r_sum = 0.d0 if (present(OO_da)) OO_da = cmplx_0 if (present(OO_dadb)) OO_dadb = cmplx_0 ! rotate wannier centres from cartesian to fractional coordinates wannier_centres_frac(:, :) = 0.d0 call utility_inverse_mat(real_lattice, inv_lattice) do ir = 1, num_wann call utility_cart_to_frac(wigner_seitz%wannier_centres_from_AA_R(:, ir), & wannier_centres_frac(:, ir), inv_lattice) end do ! print *, 'wannier_centres_frac' ! do ir = 1,num_wann ! print *, wannier_centres_frac(:,ir) ! enddo ! stop ! ! print *, 'crvec' ! do ir = 1,nrpts ! print *, crvec(:,ir) ! enddo ! stop ! print *, 'wannier_centres' ! do ir = 1,num_wann ! print *, wannier_centres(:,ir) ! enddo ! stop do ir = 1, wigner_seitz%nrpts_pw90 ! [lp] Original code, without IJ-dependent shift: do j = 1, num_wann do i = 1, num_wann r_sum(:) = real(wigner_seitz%irvec_pw90(:, ir)) & + wannier_centres_frac(:, j) - wannier_centres_frac(:, i) rdotk = twopi*dot_product(kpt(:), r_sum(:)) phase_fac = cmplx(cos(rdotk), sin(rdotk), dp) if (present(OO_da)) then ! if we are at the origin and at the same band, then the ! matrix element is zero in this convention if ((wigner_seitz%irvec_pw90(1, ir) .eq. 0) .and. (wigner_seitz%irvec_pw90(2, ir) .eq. 0) .and. & (wigner_seitz%irvec_pw90(3, ir) .eq. 0) .and. (i .eq. j)) then OO_da(i, j, 1) = OO_da(i, j, 1) + phase_fac*(OO_R(i, j, ir, 1) & - wigner_seitz%wannier_centres_from_AA_R(1, j)) OO_da(i, j, 2) = OO_da(i, j, 2) + phase_fac*(OO_R(i, j, ir, 2) & - wigner_seitz%wannier_centres_from_AA_R(2, j)) OO_da(i, j, 3) = OO_da(i, j, 3) + phase_fac*(OO_R(i, j, ir, 3) & - wigner_seitz%wannier_centres_from_AA_R(3, j)) ! print *, 'OO_R(i,j,ir,1)', OO_R(i,j,ir,1) ! print *, 'local_wannier_centres(1,j)', local_wannier_centres(1,j) ! print *, 'OO_R(i,j,ir,2)', OO_R(i,j,ir,2) ! print *, 'local_wannier_centres(2,j)', local_wannier_centres(2,j) ! cycle else OO_da(i, j, 1) = OO_da(i, j, 1) + phase_fac*OO_R(i, j, ir, 1) OO_da(i, j, 2) = OO_da(i, j, 2) + phase_fac*OO_R(i, j, ir, 2) OO_da(i, j, 3) = OO_da(i, j, 3) + phase_fac*OO_R(i, j, ir, 3) end if end if if (present(OO_dadb)) then ! same skip as before if ((wigner_seitz%irvec_pw90(1, ir) .eq. 0) .and. (wigner_seitz%irvec_pw90(2, ir) .eq. 0) .and. & (wigner_seitz%irvec_pw90(3, ir) .eq. 0) .and. (i .eq. j)) then do a = 1, 3 do b = 1, 3 OO_dadb(i, j, a, b) = OO_dadb(i, j, a, b) + & cmplx_i*(wigner_seitz%crvec_pw90(b, ir) + & wigner_seitz%wannier_centres_from_AA_R(b, j) & - wigner_seitz%wannier_centres_from_AA_R(b, i))*phase_fac* & (OO_R(i, j, ir, a) - wigner_seitz%wannier_centres_from_AA_R(a, j)) end do end do ! cycle else do a = 1, 3 do b = 1, 3 OO_dadb(i, j, a, b) = OO_dadb(i, j, a, b) + & cmplx_i*(wigner_seitz%crvec_pw90(b, ir) + & wigner_seitz%wannier_centres_from_AA_R(b, j) & - wigner_seitz%wannier_centres_from_AA_R(b, i))* & phase_fac*OO_R(i, j, ir, a) end do end do end if end if end do end do end do end subroutine pw90common_fourier_R_to_k_vec_dadb_TB_conv !================================================! ! PRIVATE PROCEDURES !================================================! !================================================! subroutine wignerseitz(print_output, real_lattice, mp_grid, ws_region, wigner_seitz, stdout, & count_pts, timer, error, comm) !================================================! !! Calculates a grid of lattice vectors r that fall inside (and eventually !! on the surface of) the Wigner-Seitz supercell centered on the !! origin of the Bravais superlattice with primitive translations !! mp_grid(1)*a_1, mp_grid(2)*a_2, and mp_grid(3)*a_3 !================================================! use w90_constants, only: eps8, dp use w90_io, only: io_stopwatch_start, io_stopwatch_stop use w90_types, only: print_output_type, timer_list_type, ws_region_type use w90_utility, only: utility_metric use w90_comms, only: w90_comm_type, mpirank use w90_postw90_types, only: wigner_seitz_type ! irvec(i,irpt) The irpt-th Wigner-Seitz grid point has components ! irvec(1:3,irpt) in the basis of the lattice vectors ! ndegen(irpt) Weight of the irpt-th point is 1/ndegen(irpt) ! nrpts number of Wigner-Seitz grid points implicit none ! arguments type(ws_region_type), intent(in) :: ws_region type(print_output_type), intent(in) :: print_output type(timer_list_type), intent(inout) :: timer type(w90_comm_type), intent(in) :: comm type(wigner_seitz_type), intent(inout) :: wigner_seitz type(w90_error_type), allocatable, intent(out) :: error integer, intent(in) :: mp_grid(3) integer, intent(in) :: stdout logical, intent(in) :: count_pts real(kind=dp), intent(in) :: real_lattice(3, 3) ! local variables integer :: ierr, dist_dim integer :: n1, n2, n3, i1, i2, i3, icnt, i, j, ir integer :: ndiff(3) logical :: on_root = .false. real(kind=dp), allocatable :: dist(:) real(kind=dp) :: real_metric(3, 3) real(kind=dp) :: tot, dist_min if (mpirank(comm) == 0) on_root = .true. if (print_output%timing_level > 1 .and. on_root) & call io_stopwatch_start('postw90_common: wigner_seitz', timer) call utility_metric(real_lattice, real_metric) dist_dim = 1 do i = 1, 3 dist_dim = dist_dim*((ws_region%ws_search_size(i) + 1)*2 + 1) end do allocate (dist(dist_dim), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating dist in wigner_seitz', comm) return end if ! The Wannier functions live in a periodic supercell of the real space unit ! cell. This supercell is mp_grid(i) unit cells long along each primitive ! translation vector a_i of the unit cell ! ! We loop over grid points r on a unit cell that is (2*ws_search_size+1)**3 times ! larger than this primitive supercell. ! ! One of these points is in the W-S supercell if it is closer to R=0 than any ! of the other points R (where R are the translation vectors of the ! supercell). ! In the end, nrpts contains the total number of grid points that have been ! found in the Wigner-Seitz cell wigner_seitz%nrpts = 0 ! Loop over the lattice vectors of the primitive cell ! that live in a supercell which is (2*ws_search_size+1)**2 ! larger than the Born-von Karman supercell. ! We need to find which among these live in the Wigner-Seitz cell do n1 = -ws_region%ws_search_size(1)*mp_grid(1), ws_region%ws_search_size(1)*mp_grid(1) do n2 = -ws_region%ws_search_size(2)*mp_grid(2), ws_region%ws_search_size(2)*mp_grid(2) do n3 = -ws_region%ws_search_size(3)*mp_grid(3), ws_region%ws_search_size(3)*mp_grid(3) ! Loop over the lattice vectors R of the Born-von Karman supercell ! that contains all the points of the previous loop. ! There are (2*(ws_search_size+1)+1)**3 points R. R=0 corresponds to ! i1=i2=i3=0, or icnt=((2*(ws_search_size+1)+1)**3 + 1)/2 icnt = 0 do i1 = -ws_region%ws_search_size(1) - 1, ws_region%ws_search_size(1) + 1 do i2 = -ws_region%ws_search_size(2) - 1, ws_region%ws_search_size(2) + 1 do i3 = -ws_region%ws_search_size(3) - 1, ws_region%ws_search_size(3) + 1 icnt = icnt + 1 ! Calculate distance squared |r-R|^2 ndiff(1) = n1 - i1*mp_grid(1) ndiff(2) = n2 - i2*mp_grid(2) ndiff(3) = n3 - i3*mp_grid(3) dist(icnt) = 0.0_dp do i = 1, 3 do j = 1, 3 dist(icnt) = dist(icnt) + & real(ndiff(i), dp)*real_metric(i, j)*real(ndiff(j), dp) end do end do end do end do end do dist_min = minval(dist) if (abs(dist((dist_dim + 1)/2) - dist_min) .lt. ws_region%ws_distance_tol**2) then wigner_seitz%nrpts = wigner_seitz%nrpts + 1 if (.not. count_pts) then wigner_seitz%ndegen(wigner_seitz%nrpts) = 0 do i = 1, dist_dim if (abs(dist(i) - dist_min) .lt. ws_region%ws_distance_tol**2) & wigner_seitz%ndegen(wigner_seitz%nrpts) = & wigner_seitz%ndegen(wigner_seitz%nrpts) + 1 end do wigner_seitz%irvec(1, wigner_seitz%nrpts) = n1 wigner_seitz%irvec(2, wigner_seitz%nrpts) = n2 wigner_seitz%irvec(3, wigner_seitz%nrpts) = n3 ! Remember which grid point r is at the origin if (n1 == 0 .and. n2 == 0 .and. n3 == 0) wigner_seitz%rpt_origin = wigner_seitz%nrpts end if end if !n3 end do !n2 end do !n1 end do ! deallocate (dist, stat=ierr) if (ierr /= 0) then call set_error_dealloc(error, 'Error in deallocating dist wigner_seitz', comm) return end if if (count_pts) then if (print_output%timing_level > 1 .and. on_root) & call io_stopwatch_stop('postw90_common: wigner_seitz', timer) return end if if (print_output%iprint >= 3 .and. on_root) then write (stdout, '(1x,i4,a,/)') wigner_seitz%nrpts, & ' lattice points in Wigner-Seitz supercell:' do ir = 1, wigner_seitz%nrpts write (stdout, '(4x,a,3(i3,1x),a,i2)') ' vector ', wigner_seitz%irvec(1, ir), & wigner_seitz%irvec(2, ir), wigner_seitz%irvec(3, ir), ' degeneracy: ', & wigner_seitz%ndegen(ir) end do !write (stdout, '(1x,a,f12.3)') ' tot = ', tot !not yet computed write (stdout, '(1x,a,i12)') ' mp_grid product = ', mp_grid(1)*mp_grid(2)*mp_grid(3) end if ! Check the "sum rule" tot = 0.0_dp do ir = 1, wigner_seitz%nrpts ! ! Corrects weights in Fourier sums for R-vectors on the boundary of the ! W-S supercell ! tot = tot + 1.0_dp/real(wigner_seitz%ndegen(ir), dp) end do if (abs(tot - real(mp_grid(1)*mp_grid(2)*mp_grid(3), dp)) > eps8) then call set_error_fatal(error, 'ERROR in wigner_seitz: error in finding Wigner-Seitz points', comm) return end if if (print_output%timing_level > 1 .and. on_root) & call io_stopwatch_stop('postw90_common: wigner_seitz', timer) return end subroutine wignerseitz !=========================================================! subroutine wigner_seitz_opt_setup(ws_distance, ws_region, wigner_seitz, num_wann, real_lattice, & nrpts, error, comm) !========================================================================== ! !! If use_ws_distance = true, get list of all R vectors that appear in !! irdist_ws and write them to irvec_pw90. Set the index map ir_ind_ws_to_pw90 !! !! If use_ws_distance = false, copy irvec to irvec_pw90. ! !========================================================================== use w90_constants, only: dp use w90_types, only: ws_region_type, ws_distance_type use w90_comms, only: w90_comm_type use w90_postw90_types, only: wigner_seitz_type implicit none type(ws_distance_type), intent(inout) :: ws_distance type(ws_region_type), intent(in) :: ws_region type(wigner_seitz_type), intent(inout) :: wigner_seitz type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm real(kind=dp), intent(in) :: real_lattice(3, 3) integer, intent(in) :: num_wann integer, intent(in) :: nrpts logical :: found integer :: i, j, ideg, ir, jr, ierr, nrpts_found, max_ndeg, ivdum(3) integer, allocatable :: irvec_found(:, :), irvec_temp(:, :) if (ws_region%use_ws_distance) then max_ndeg = maxval(ws_distance%ndeg) allocate (wigner_seitz%ir_ind_ws_to_pw90(max_ndeg, num_wann, num_wann, nrpts), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating ir_ind_ws_to_pw90 in wigner_seitz_opt_setup', comm) return end if wigner_seitz%ir_ind_ws_to_pw90 = -1 ! find the set of R vectors from irdist_ws, removing duplicates nrpts_found = 0 do ir = 1, nrpts do j = 1, num_wann do i = 1, num_wann do ideg = 1, ws_distance%ndeg(i, j, ir) ivdum = ws_distance%irdist(:, ideg, i, j, ir) if (nrpts_found == 0) then nrpts_found = 1 allocate (irvec_found(3, 1), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec_found in wigner_seitz_opt_setup', comm) return end if irvec_found(:, 1) = ivdum wigner_seitz%ir_ind_ws_to_pw90(ideg, i, j, ir) = 1 else ! find if ivdum is already found found = .false. do jr = 1, nrpts_found if (all(ivdum == irvec_found(:, jr))) then wigner_seitz%ir_ind_ws_to_pw90(ideg, i, j, ir) = jr found = .true. exit end if end do ! if not found, add ivdum to irvec_found if (.not. found) then ! copy irvec_found to irvec_temp allocate (irvec_temp(3, nrpts_found), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec_temp in wigner_seitz_opt_setup', comm) return end if irvec_temp = irvec_found deallocate (irvec_found, stat=ierr) if (ierr /= 0) then call set_error_dealloc(error, 'Error in deallocating irvec_found in wigner_seitz_opt_setup', comm) return end if allocate (irvec_found(3, nrpts_found + 1), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec_found in wigner_seitz_opt_setup', comm) return end if irvec_found(:, 1:nrpts_found) = irvec_temp irvec_found(:, nrpts_found + 1) = ivdum deallocate (irvec_temp, stat=ierr) if (ierr /= 0) then call set_error_dealloc(error, 'Error in deallocating irvec_temp in wigner_seitz_opt_setup', comm) return end if nrpts_found = nrpts_found + 1 wigner_seitz%ir_ind_ws_to_pw90(ideg, i, j, ir) = nrpts_found end if ! found end if ! nrpts_found == 0 end do end do end do end do wigner_seitz%nrpts_pw90 = nrpts_found allocate (wigner_seitz%irvec_pw90(3, wigner_seitz%nrpts_pw90), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec_pw90 in wigner_seitz_opt_setup', comm) return end if allocate (wigner_seitz%crvec_pw90(3, wigner_seitz%nrpts_pw90), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating crvec_pw90 in wigner_seitz_opt_setup', comm) return end if wigner_seitz%irvec_pw90 = irvec_found do ir = 1, wigner_seitz%nrpts_pw90 wigner_seitz%crvec_pw90(:, ir) = matmul(transpose(real_lattice), real(wigner_seitz%irvec_pw90(:, ir), dp)) end do else ! .not. use_ws_distance allocate (wigner_seitz%irvec_pw90(3, wigner_seitz%nrpts_pw90), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating irvec_pw90 in wigner_seitz_opt_setup', comm) return end if allocate (wigner_seitz%crvec_pw90(3, wigner_seitz%nrpts_pw90), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error in allocating crvec_pw90 in wigner_seitz_opt_setup', comm) return end if wigner_seitz%irvec_pw90 = wigner_seitz%irvec wigner_seitz%crvec_pw90 = wigner_seitz%crvec wigner_seitz%nrpts_pw90 = wigner_seitz%nrpts end if ! use_ws_distance ! Check degeneracy factor ndegen for R = 0 is 1 if (wigner_seitz%ndegen(wigner_seitz%rpt_origin) /= 1) then call set_error_fatal(error, 'ndegen for R=0 is not 1.', comm) return end if ! Check degeneracy factor ws_distance%ndeg for a Wannier function with itself, ! i.e. R = 0 and i = j, is 1. if (ws_region%use_ws_distance) then do ir = 1, nrpts do i = 1, num_wann do ideg = 1, ws_distance%ndeg(i, i, ir) ivdum = ws_distance%irdist(:, ideg, i, i, ir) if (all(ivdum == 0)) then if (ws_distance%ndeg(i, i, ir) /= 1) then call set_error_fatal(error, 'ws_distance%ndeg for R=0 and i=j is not 1.', comm) return end if end if end do end do end do end if end subroutine wigner_seitz_opt_setup end module w90_postw90_common