1!===============================================================================
2!
3! Routines:
4!
5!
6!===============================================================================
7
8#include "f_defs.h"
9
10module read_rho_vxc_m
11
12  use global_m
13  use gmap_m
14  use wfn_rho_vxc_io_m
15
16  implicit none
17
18  private
19
20  public :: read_vxc, read_rho
21
22contains
23
24  !---------------------------------------------------------------------------------------------------
25  !> Read in the exchange-correlation potential and store in array sig%vxc
26  subroutine read_vxc(sig, gvec, kp, syms, crys, isrti, isrt, vxc_type)
27    type(siginfo), intent(inout) :: sig
28    type(gspace), intent(in) :: gvec
29    type(kpoints), intent(in) :: kp
30    type(symmetry), intent(in) :: syms
31    type(crystal), intent(in) :: crys
32    integer, intent(in) :: isrti(:) !< (gvec%ng)
33    integer, intent(in) :: isrt(:) !< (gvec%ng)
34    integer, intent(in) :: vxc_type
35
36    character*3 :: sheader
37    integer :: iflavor, ig, ispin
38    type(gspace) :: gvec_dummy
39    type(crystal) :: crys_dummy
40    type(symmetry) :: syms_dummy
41    type(kpoints) :: kp_dummy
42    real(DP) :: discrepancy
43
44    PUSH_SUB(read_vxc)
45
46    call logit('Reading VXC')
47    if (vxc_type .eq. 1) then
48      SAFE_ALLOCATE(sig%vxc, (gvec%ng,kp%nspin))
49    else if (vxc_type .eq. 2) then
50      SAFE_ALLOCATE(sig%vxc2, (gvec%ng,kp%nspin))
51    else if (vxc_type .eq. 3) then
52      SAFE_ALLOCATE(sig%dvxc, (gvec%ng,kp%nspin))
53    else
54      call die("VXC type can only be 1 (VXC), 2 (VXC2), and 3 (dVXC).")
55    endif
56
57    if (vxc_type .eq. 1) then
58      if(peinf%inode == 0) call open_file(96,file='VXC',form='unformatted',status='old')
59    else if (vxc_type .eq. 2) then
60      if(peinf%inode == 0) call open_file(96,file='VXC2',form='unformatted',status='old')
61    else if (vxc_type .eq. 3) then
62      if(peinf%inode == 0) call open_file(96,file='dVXC',form='unformatted',status='old')
63    endif
64
65    sheader = 'VXC'
66    iflavor = 0
67    call read_binary_header_type(96, sheader, iflavor, kp_dummy, gvec_dummy, syms_dummy, &
68         crys_dummy, warn = .false.)
69
70    call check_header('WFN_inner', kp, gvec, syms, crys, 'VXC', kp_dummy, gvec_dummy, &
71         syms_dummy, crys_dummy, is_wfn = .false.)
72
73    SAFE_ALLOCATE(gvec_dummy%components, (3, gvec_dummy%ng))
74    call read_binary_gvectors(96, gvec_dummy%ng, gvec_dummy%ng, gvec_dummy%components)
75    do ig = 1, gvec%ng
76      if(any(gvec_dummy%components(:,isrt(ig)) .ne. gvec%components(:,ig))) call die("gvec mismatch in VXC")
77    enddo
78    SAFE_DEALLOCATE_P(gvec_dummy%components)
79
80    if (vxc_type .eq. 1) then
81      call read_binary_data(96, gvec_dummy%ng, gvec_dummy%ng, kp%nspin, sig%vxc, gindex = isrti)
82    else if (vxc_type .eq. 2) then
83      call read_binary_data(96, gvec_dummy%ng, gvec_dummy%ng, kp%nspin, sig%vxc2, gindex = isrti)
84    else if (vxc_type .eq. 3) then
85      call read_binary_data(96, gvec_dummy%ng, gvec_dummy%ng, kp%nspin, sig%dvxc, gindex = isrti)
86    endif
87
88    if(peinf%inode == 0) then
89      call close_file(96)
90
91      if (vxc_type .eq. 1) then
92        do ispin = 1, kp%nspin
93          discrepancy = check_field_is_real(sig%vxc(:, ispin), gvec)
94          if(discrepancy > TOL_Zero) then
95            write(0,*) 'WARNING: VXC is not real in real space, with discrepancy ', discrepancy, ' for spin ', ispin
96          endif
97        enddo
98      else if (vxc_type .eq. 2) then
99        do ispin = 1, kp%nspin
100          discrepancy = check_field_is_real(sig%vxc2(:, ispin), gvec)
101          if(discrepancy > TOL_Zero) then
102            write(0,*) 'WARNING: VXC2 is not real in real space, with discrepancy ', discrepancy, ' for spin ', ispin
103          endif
104        enddo
105      else if (vxc_type .eq. 3) then
106        write (0,*) 'NOTE: dVXC is complex in general.'
107      endif
108    endif
109
110    call dealloc_header_type(sheader, crys_dummy, kp_dummy)
111
112    PUSH_SUB(read_vxc)
113
114  end subroutine read_vxc
115
116  !---------------------------------------------------------------------------------------------------
117  !> Read in the charge density and store in array wpg%rho (formerly known as CD95)
118  !! CD95 Ref: http://www.nature.com/nature/journal/v471/n7337/full/nature09897.html
119  subroutine read_rho(wpg, gvec, kp, syms, crys, isrti, isrt, check_filename)
120    type(wpgen), intent(inout) :: wpg
121    type(gspace), intent(in) :: gvec
122    type(kpoints), intent(in) :: kp
123    type(symmetry), intent(in) :: syms
124    type(crystal), intent(in) :: crys
125    integer, intent(in) :: isrti(:) !< (gvec%ng)
126    integer, intent(in) :: isrt(:) !< (gvec%ng)
127    character(len=*), intent(in) :: check_filename !< This is the file
128                                            !< against which the header is
129                                            !< checked
130
131
132    character*3 :: sheader
133    integer :: iflavor, ig, ispin
134    type(gspace) :: gvec_dummy
135    type(crystal) :: crys_dummy
136    type(symmetry) :: syms_dummy
137    type(kpoints) :: kp_dummy
138    real(DP) :: discrepancy
139
140    PUSH_SUB(read_rho)
141
142    call logit('Reading RHO')
143
144    SAFE_ALLOCATE(wpg%rho, (gvec%ng,kp%nspin))
145
146    if(peinf%inode == 0) call open_file(95,file='RHO',form='unformatted',status='old')
147
148    sheader = 'RHO'
149    iflavor = 0
150    call read_binary_header_type(95, sheader, iflavor, kp_dummy, gvec_dummy, &
151        syms_dummy, crys_dummy, warn = .false.)
152
153    call check_header(trim(check_filename), kp, gvec, syms, crys, 'RHO', kp_dummy, &
154        gvec_dummy, syms_dummy, crys_dummy, is_wfn = .false.)
155
156    SAFE_ALLOCATE(gvec_dummy%components, (3, gvec_dummy%ng))
157    call read_binary_gvectors(95, gvec_dummy%ng, gvec_dummy%ng, gvec_dummy%components)
158    do ig = 1, gvec%ng
159      if(any(gvec_dummy%components(:,isrt(ig)) .ne. gvec%components(:,ig))) call die("gvec mismatch in RHO")
160    enddo
161    SAFE_DEALLOCATE_P(gvec_dummy%components)
162
163    call read_binary_data(95, gvec_dummy%ng, gvec_dummy%ng, kp%nspin, wpg%rho, gindex = isrti)
164
165    if(peinf%inode == 0) call close_file(95)
166
167    ! otherwise if nspin == 1, the 2 component may be uninitialized to NaN
168    wpg%wpsq(1:2) = 0d0
169    wpg%nelec(1:2) = 0d0
170
171    ! since they are sorted, if G = 0 is present, it is the first one
172    if(any(gvec%components(1:3, 1) /= 0)) call die("gvectors for RHO must include G = 0")
173    ! otherwise, the code below will not do what we think it does
174
175#ifdef CPLX
176    if(any(abs(aimag(wpg%rho(1,:))) > TOL_Zero)) then
177      call die("Charge density in RHO has imaginary part for G=0", only_root_writes = .true.)
178    endif
179#endif
180
181    do ispin=1,kp%nspin
182      wpg%nelec(ispin)=dble(wpg%rho(1,ispin))
183      wpg%wpsq(ispin)=ryd*ryd*16.0d0*PI_D*wpg%nelec(ispin)/crys%celvol
184    enddo
185
186    ! This is unacceptable because it means the number of electrons is negative,
187    ! and the plasma frequency will be imaginary!
188    if(any(wpg%nelec(1:kp%nspin) < TOL_Zero)) then
189      write(0,*) wpg%nelec(:)
190      call die("Charge density in RHO has negative part for G=0", only_root_writes = .true.)
191    endif
192
193    if(peinf%inode == 0) then
194      do ispin = 1, kp%nspin
195        discrepancy = check_field_is_real(wpg%rho(:, ispin), gvec)
196        if(discrepancy > TOL_Zero) then
197          write(0,*) 'WARNING: RHO is not real in real space, with discrepancy ', discrepancy, ' for spin ', ispin
198        endif
199      enddo
200    endif
201
202    call dealloc_header_type(sheader, crys_dummy, kp_dummy)
203
204    POP_SUB(read_rho)
205  end subroutine read_rho
206
207  !---------------------------------------------------------------------------------------------------
208  !> RHO and VXC must be real in real space, i.e. c(G) - c(-G)* = 0
209  real(DP) function check_field_is_real(field, gvec)
210    SCALAR, intent(in) :: field(:)
211    type(gspace), intent(in) :: gvec
212
213    integer :: ig, umklapp(3)
214    SCALAR :: diff
215    type(symmetry) :: syms_inv
216    integer, allocatable :: ind(:), identity(:)
217    SCALAR, allocatable :: phase(:)
218
219    PUSH_SUB(check_field_is_real)
220
221    syms_inv%ntran = 1
222    syms_inv%mtrx = 0
223    syms_inv%tnp = 0d0
224    syms_inv%mtrx(1:3, 1:3, 1) = reshape((/-1, 0, 0, 0, -1, 0, 0, 0, -1/), shape(syms_inv%mtrx(:,:,1)))
225    umklapp(1:3) = 0
226
227    SAFE_ALLOCATE(phase, (gvec%ng))
228    SAFE_ALLOCATE(ind, (gvec%ng))
229    SAFE_ALLOCATE(identity, (gvec%ng))
230
231    do ig = 1, gvec%ng
232      identity(ig) = ig
233    enddo
234
235    call gmap(gvec, syms_inv, gvec%ng, 1, umklapp, identity, identity, ind, phase, die_outside_sphere = .true.)
236    if(any(abs(phase(1:gvec%ng) - 1d0) > TOL_Zero)) call die("non-unity phase in check_field_is_real")
237    if(any(ind(1:gvec%ng) == 0)) call die("ind array from gmap has a zero")
238
239    check_field_is_real = 0d0
240    do ig = 1, gvec%ng
241      diff = field(ig) - MYCONJG(field(ind(ig))*phase(ig))
242      check_field_is_real = max(check_field_is_real, abs(diff))
243    enddo
244
245    SAFE_DEALLOCATE(phase)
246    SAFE_DEALLOCATE(ind)
247    SAFE_DEALLOCATE(identity)
248
249    POP_SUB(check_field_is_real)
250  end function check_field_is_real
251
252end module read_rho_vxc_m
253