1!==============================================================================
2!
3! Module bse_init_m
4!
5! (1) bse_init()         Originally By FHJ       Last Modified  3/24/2012 (FHJ)
6!
7!     Performs initializations required by absorption and inteqp calculations.
8!
9!==============================================================================
10
11#include "f_defs.h"
12
13module bse_init_m
14
15  use global_m
16  use fullbz_m
17  use intpts_m
18  use wfn_rho_vxc_io_m
19
20  implicit none
21
22  private
23
24  public :: bse_init
25
26contains
27
28!> FHJ: Figures out the number of k-points in the coarse grid and the
29!!      dimensionality of the system for interpolation pourposes.
30!! FIXME: am I freeing all buffers correctly?
31subroutine bse_init(xct,flag)
32  type (xctinfo), intent(inout) :: xct
33  type (flags), intent(in) :: flag
34
35  character(len=3) :: sheader
36  integer :: iflavor
37  type (crystal) :: crys, crys_co
38  type (gspace) :: gvec, gvec_co
39  type (grid) :: kg, kg_co
40  type (kpoints) :: kp, kp_co
41  type (symmetry) :: syms, syms_co
42  logical :: skip_checkbz
43
44  logical :: is_periodic_old(3)
45  integer :: jdim, npts_intp_kernel
46
47  PUSH_SUB(bse_init)
48
49  ! GKA: No additional information is needed if the eigenvalues are already computed
50  if (flag%spec.eq.1) then
51    return
52  endif
53
54  if (flag%read_dtmat) then
55    call open_file(unit=13,file='dtmat',form='unformatted',status='old')
56    read(13) xct%idimensions, xct%is_periodic(1:3), npts_intp_kernel, xct%nkpt_co
57    if (peinf%inode==0) then
58      if (xct%npts_intp_kernel==-1 .and. npts_intp_kernel==1) then
59        write(0,*)
60        write(0,'(a)') 'WARNING: dtmat was not calculated with "kernel_k_interpolation",'
61        write(0,'(a)') 'so we are turning off "kernel_k_interpolation".'
62        write(0,*)
63      elseif (xct%npts_intp_kernel==1 .and. npts_intp_kernel>1) then
64        write(0,*)
65        write(0,'(a)') 'WARNING: dtmat was calculated with "kernel_k_interpolation",'
66        write(0,'(a)') 'so we are turning on "kernel_k_interpolation".'
67        write(0,*)
68      endif
69    endif
70    xct%npts_intp_kernel = npts_intp_kernel
71    call close_file(13)
72    return
73  endif
74
75  if (peinf%inode==0) write(6,*)
76
77  ! FHJ: Read the header of WFN_co
78  if (peinf%inode == 0) &
79    call open_file(unit=25, file='WFN_co', form='unformatted', status='old')
80  sheader = 'WFN'
81  iflavor = 0
82  call read_binary_header_type(25, sheader, iflavor, kp_co, gvec_co, &
83    syms_co, crys_co, warn=.false., dont_warn_kgrid=.true.)
84  if (peinf%inode == 0) call close_file(25)
85  kg_co%nr = kp_co%nrk
86  SAFE_ALLOCATE(kg_co%r, (3,kg_co%nr))
87  kg_co%r(:,:) = kp_co%rk(:,:)
88
89  ! FHJ: We just need to unfold the k-points fast to get the number of k-points
90  !      in each direction, so we won`t build the wigner seitz cell.
91  if (flag%bzc==1) then
92    call fullbz(crys_co, syms_co, kg_co, 1, skip_checkbz, &
93                wigner_seitz=.false., paranoid=.false., do_nothing=.true.)
94  else
95    call fullbz(crys_co, syms_co, kg_co, syms_co%ntran, skip_checkbz, &
96                wigner_seitz=.false., paranoid=.false.)
97  endif
98
99  xct%nkpt_co = kg_co%nf
100  if (peinf%inode==0) then
101    write(6,'(1x,a,i0)') 'Number of k-points in the coarse k-grid: ', xct%nkpt_co
102  endif
103
104  ! FHJ: find the number of periodic dimensions from the kpt sampling.
105  !      Don`t trust kp%kgrid!
106  call get_ndims(kg_co, xct)
107
108  SAFE_DEALLOCATE_P(kg_co%r)
109  SAFE_DEALLOCATE_P(kg_co%f)
110
111
112  if (xct%skipinterp) then
113    xct%npts_intp_kernel = 1
114  else
115    ! FHJ: Read the header of WFN, if we are interpolating
116    if (peinf%inode == 0) &
117      call open_file(25, file='WFN_fi', form='unformatted', status='old')
118    sheader = 'WFN'
119    iflavor = 0
120    call read_binary_header_type(25, sheader, iflavor, kp, gvec, syms, crys, &
121      dont_warn_kgrid=.true.)
122    if (peinf%inode == 0) call close_file(25)
123    kg%nr = kp%nrk
124    SAFE_ALLOCATE(kg%r, (3,kg%nr))
125    kg%r(:,:) = kp%rk(:,:)
126
127    if (flag%bz0==0.and.xct%is_absorption) then
128      call fullbz(crys, syms, kg, syms%ntran, skip_checkbz, &
129                  wigner_seitz=.false., paranoid=.false.)
130    else
131      call fullbz(crys, syms, kg, 1, skip_checkbz, &
132                  wigner_seitz=.false., paranoid=.false., do_nothing=.true.)
133    endif
134
135    is_periodic_old = xct%is_periodic
136    call get_ndims(kg, xct)
137    xct%is_periodic = xct%is_periodic .or. is_periodic_old
138    xct%idimensions = 0
139    do jdim = 1, 3
140      if (xct%is_periodic(jdim)) xct%idimensions = xct%idimensions + 1
141    enddo
142
143    SAFE_DEALLOCATE_P(kg%r)
144    SAFE_DEALLOCATE_P(kg%f)
145
146    if (xct%npts_intp_kernel==-1) xct%npts_intp_kernel = xct%idimensions + 1
147
148    if (peinf%inode==0) then
149      write(6,'(1x,a,i1,a)') 'A ',xct%idimensions,'-D interpolation algorithm will be employed.'
150    endif
151  endif ! xct%skipinterp
152
153  POP_SUB(bse_init)
154
155end subroutine bse_init
156
157end module bse_init_m
158
159