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