1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! Copyright (C) 2012-2013 M. Gruning, P. Melo, M. Oliveira 3!! 4!! This program is free software; you can redistribute it and/or modify 5!! it under the terms of the GNU General Public License as published by 6!! the Free Software Foundation; either version 2, or (at your option) 7!! any later version. 8!! 9!! This program is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU General Public License for more details. 13!! 14!! You should have received a copy of the GNU General Public License 15!! along with this program; if not, write to the Free Software 16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 17!! 02110-1301, USA. 18!! 19 20#include "global.h" 21 22module xc_oep_oct_m 23 use comm_oct_m 24 use derivatives_oct_m 25 use exchange_operator_oct_m 26 use geometry_oct_m 27 use global_oct_m 28 use grid_oct_m 29 use hamiltonian_elec_oct_m 30 use lalg_basic_oct_m 31 use lalg_adv_oct_m 32 use linear_response_oct_m 33 use linear_solver_oct_m 34 use mesh_function_oct_m 35 use mesh_oct_m 36 use messages_oct_m 37 use mpi_oct_m 38 use multicomm_oct_m 39 use namespace_oct_m 40 use parser_oct_m 41 use photon_mode_oct_m 42 use poisson_oct_m 43 use profiling_oct_m 44 use states_abst_oct_m 45 use states_elec_oct_m 46 use states_elec_calc_oct_m 47 use states_elec_dim_oct_m 48 use scf_tol_oct_m 49 use varinfo_oct_m 50 use xc_oct_m 51 use XC_F90(lib_m) 52 use xc_functl_oct_m 53 54 implicit none 55 56 private 57 public :: & 58 xc_oep_t, & 59 xc_oep_init, & 60 xc_oep_end, & 61 xc_oep_write_info, & 62 dxc_oep_calc, & 63 zxc_oep_calc, & 64 doep_x, & 65 zoep_x 66 67 !> the OEP levels 68 integer, public, parameter :: & 69 XC_OEP_NONE = 1, & 70 XC_OEP_KLI = 3, & 71 XC_OEP_FULL = 5, & 72 OEP_MIXING_SCHEME_CONST = 1, & 73 OEP_MIXING_SCHEME_BB = 2, & 74 OEP_MIXING_SCHEME_DENS = 3 75 76 type xc_oep_t 77 private 78 integer, public :: level !< 0 = no oep, 1 = Slater, 2 = KLI, 4 = full OEP 79 FLOAT :: mixing !< how much of the function S(r) to add to vxc in every iteration 80 type(lr_t) :: lr !< to solve the equation H psi = b 81 type(linear_solver_t) :: solver 82 type(scf_tol_t) :: scftol 83 integer :: eigen_n 84 integer, pointer :: eigen_type(:), eigen_index(:) 85 FLOAT :: socc, sfact 86 FLOAT, pointer, public :: vxc(:,:), uxc_bar(:,:) 87 FLOAT, pointer :: dlxc(:, :, :) 88 CMPLX, pointer :: zlxc(:, :, :) 89 integer :: mixing_scheme 90 logical, public :: has_photons ! one-photon OEP 91 type(photon_mode_t), public :: pt 92 type(lr_t) :: photon_lr !< to solve the equation H psi = b 93 FLOAT, public :: norm2ss 94 FLOAT, pointer :: vxc_old(:,:), ss_old(:,:) 95 integer :: noccst 96 logical :: coc_translation 97 end type xc_oep_t 98 99 type(profile_t), save :: & 100 C_PROFILING_XC_OEP, & 101 C_PROFILING_XC_EXX, & 102 C_PROFILING_XC_SIC, & 103 C_PROFILING_XC_OEP_FULL, & 104 C_PROFILING_XC_KLI 105 106contains 107 108 ! --------------------------------------------------------- 109 subroutine xc_oep_init(oep, namespace, family, gr, st, geo, mc) 110 type(xc_oep_t), intent(out) :: oep 111 type(namespace_t), intent(in) :: namespace 112 integer, intent(in) :: family 113 type(grid_t), intent(inout) :: gr 114 type(states_elec_t), intent(in) :: st 115 type(geometry_t), intent(in) :: geo 116 type(multicomm_t), intent(in) :: mc 117 118 PUSH_SUB(xc_oep_init) 119 120 if(bitand(family, XC_FAMILY_OEP) == 0) then 121 oep%level = XC_OEP_NONE 122 POP_SUB(xc_oep_init) 123 return 124 end if 125 126 !%Variable OEPLevel 127 !%Type integer 128 !%Default oep_kli 129 !%Section Hamiltonian::XC 130 !%Description 131 !% At what level shall <tt>Octopus</tt> handle the optimized effective potential (OEP) equation. 132 !%Option oep_none 1 133 !% Do not solve OEP equation. 134 !%Option oep_kli 3 135 !% Krieger-Li-Iafrate (KLI) approximation. For spinors, the iterative solution is controlled by the variables 136 !% in section <tt>Linear Response::Solver</tt>, and the default for <tt>LRMaximumIter</tt> is set to 50. 137 !% Ref: JB Krieger, Y Li, GJ Iafrate, <i>Phys. Lett. A</i> <b>146</b>, 256 (1990). 138 !%Option oep_full 5 139 !% (Experimental) Full solution of OEP equation using the Sternheimer approach. 140 !% The linear solver will be controlled by the variables in section <tt>Linear Response::Solver</tt>, 141 !% and the iterations for OEP by <tt>Linear Response::SCF in LR calculations</tt> and variable 142 !% <tt>OEPMixing</tt>. Note that default for <tt>LRMaximumIter</tt> is set to 10. 143 !% Ref: S. Kuemmel and J. Perdew, <i>Phys. Rev. Lett.</i> <b>90</b>, 043004 (2003). 144 !%End 145 call messages_obsolete_variable(namespace, 'OEP_Level', 'OEPLevel') 146 call parse_variable(namespace, 'OEPLevel', XC_OEP_KLI, oep%level) 147 if(.not. varinfo_valid_option('OEPLevel', oep%level)) call messages_input_error(namespace, 'OEPLevel') 148 149 if(oep%level /= XC_OEP_NONE) then 150 151 !%Variable Photons 152 !%Type logical 153 !%Default .false. 154 !%Section Hamiltonian::XC 155 !%Description 156 !% Activate the one-photon OEP 157 !%End 158 call messages_obsolete_variable(namespace, 'OEPPtX', 'Photons') 159 call parse_variable(namespace, 'Photons', .false., oep%has_photons) 160 if (oep%has_photons) then 161 call messages_experimental("Photons = yes") 162 call photon_mode_init(oep%pt, namespace, gr%mesh, gr%sb%dim, st%qtot) 163 if (oep%pt%nmodes > 1) then 164 call messages_not_implemented('Photon OEP for more than one photon mode.') 165 end if 166 if (oep%level == XC_OEP_FULL .and. st%d%nspin /= UNPOLARIZED) then 167 call messages_not_implemented('Spin-polarized calculations with photon OEP.') 168 end if 169 end if 170 171 if(oep%level == XC_OEP_FULL) then 172 173 if(st%d%nspin == SPINORS) & 174 call messages_not_implemented("Full OEP with spinors", namespace=namespace) 175 176 call messages_experimental("Full OEP") 177 !%Variable OEPMixing 178 !%Type float 179 !%Default 1.0 180 !%Section Hamiltonian::XC 181 !%Description 182 !% The linear mixing factor used to solve the Sternheimer 183 !% equation in the full OEP procedure. 184 !%End 185 call messages_obsolete_variable(namespace, 'OEP_Mixing', 'OEPMixing') 186 call parse_variable(namespace, 'OEPMixing', M_ONE, oep%mixing) 187 188 !%Variable OEPMixingScheme 189 !%Type integer 190 !%Default 1.0 191 !%Section Hamiltonian::XC 192 !%Description 193 !%Different Mixing Schemes are possible 194 !%Option OEP_MIXING_SCHEME_CONST 1 195 !%Use a constant 196 !%Reference: S. Kuemmel and J. Perdew, <i>Phys. Rev. Lett.</i> <b>90</b>, 4, 043004 (2003) 197 !%Option OEP_MIXING_SCHEME_BB 2 198 !%Use the Barzilai-Borwein (BB) Method 199 !%Reference: T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos, 200 !%<i>Phys. Rev. B</i> <b>85<\b>, 235126 (2012) 201 !%Option OEP_MIXING_SCHEME_DENS 3 202 !%Use the inverse of the electron density 203 !%Reference: S. Kuemmel and J. Perdew, <i>Phys. Rev. B</i> <b>68</b>, 035103 (2003) 204 !%End 205 call parse_variable(namespace, 'OEPMixingScheme', OEP_MIXING_SCHEME_CONST, oep%mixing_scheme) 206 207 if (oep%mixing_scheme == OEP_MIXING_SCHEME_BB) then 208 SAFE_ALLOCATE(oep%vxc_old(1:gr%mesh%np,st%d%ispin)) 209 SAFE_ALLOCATE(oep%ss_old(1:gr%mesh%np,st%d%ispin)) 210 oep%vxc_old = M_ZERO 211 oep%ss_old = M_ZERO 212 end if 213 214 oep%norm2ss = M_ZERO 215 end if 216 217 ! this routine is only prepared for finite systems. (Why not?) 218 if(st%d%nik > st%d%ispin) & 219 call messages_not_implemented("OEP for periodic systems", namespace=namespace) 220 221 ! obtain the spin factors 222 call xc_oep_SpinFactor(oep, st%d%nspin) 223 224 ! This variable will keep vxc across iterations 225 if ((st%d%ispin==3) .or. oep%level == XC_OEP_FULL) then 226 SAFE_ALLOCATE(oep%vxc(1:gr%mesh%np,st%d%nspin)) 227 else 228 SAFE_ALLOCATE(oep%vxc(1:gr%mesh%np,1:min(st%d%nspin, 2))) 229 end if 230 oep%vxc = M_ZERO 231 232 !%Variable KLIPhotonCOC 233 !%Type logical 234 !%Default .false. 235 !%Section Hamiltonian::XC 236 !%Description 237 !% Activate the center of charge translation of the electric dipole operator which should avoid the dependence of the photon KLI on an permanent dipole. 238 !%End 239 240 ! when performing full OEP, we need to solve a linear equation 241 if((oep%level == XC_OEP_FULL).or.(oep%has_photons)) then 242 call scf_tol_init(oep%scftol, namespace, st%qtot, def_maximumiter=10) 243 call linear_solver_init(oep%solver, namespace, gr, states_are_real(st), geo, mc) 244 call lr_init(oep%lr) 245 if(oep%has_photons) then 246 call lr_init(oep%photon_lr) 247 call parse_variable(namespace, 'KLIPhotonCOC', .false., oep%coc_translation) 248 end if 249 end if 250 251 ! the linear equation has to be more converged if we are to attain the required precision 252 !oep%lr%conv_abs_dens = oep%lr%conv_abs_dens / (oep%mixing) 253 254 if(st%d%nspin == SPINORS) & 255 call messages_experimental("OEP with spinors") 256 257 if(st%d%kpt%parallel) & 258 call messages_not_implemented("OEP parallel in spin/k-points", namespace=namespace) 259 260 end if 261 262 POP_SUB(xc_oep_init) 263 end subroutine xc_oep_init 264 265 266 ! --------------------------------------------------------- 267 subroutine xc_oep_end(oep) 268 type(xc_oep_t), intent(inout) :: oep 269 270 PUSH_SUB(xc_oep_end) 271 272 if(oep%level /= XC_OEP_NONE) then 273 SAFE_DEALLOCATE_P(oep%vxc) 274 if (oep%level == XC_OEP_FULL .or. oep%has_photons) then 275 call lr_dealloc(oep%lr) 276 call linear_solver_end(oep%solver) 277 end if 278 if (oep%has_photons) then 279 call lr_dealloc(oep%photon_lr) 280 call photon_mode_end(oep%pt) 281 end if 282 if (oep%level == XC_OEP_FULL .and. oep%mixing_scheme == OEP_MIXING_SCHEME_BB) then 283 SAFE_DEALLOCATE_P(oep%vxc_old) 284 SAFE_DEALLOCATE_P(oep%ss_old) 285 end if 286 end if 287 288 POP_SUB(xc_oep_end) 289 end subroutine xc_oep_end 290 291 292 ! --------------------------------------------------------- 293 subroutine xc_oep_write_info(oep, iunit) 294 type(xc_oep_t), intent(in) :: oep 295 integer, intent(in) :: iunit 296 297 if(oep%level == XC_OEP_NONE) return 298 299 PUSH_SUB(xc_oep_write_info) 300 call messages_print_var_option(iunit, 'OEPLevel', oep%level) 301 302 POP_SUB(xc_oep_write_info) 303 end subroutine xc_oep_write_info 304 305 306 ! --------------------------------------------------------- 307 !> A couple of auxiliary functions for oep 308 ! --------------------------------------------------------- 309 subroutine xc_oep_SpinFactor(oep, nspin) 310 type(xc_oep_t), intent(inout) :: oep 311 integer, intent(in) :: nspin 312 313 PUSH_SUB(xc_oep_SpinFactor) 314 315 select case(nspin) 316 case(1) ! we need to correct or the spin occupancies 317 oep%socc = M_HALF 318 oep%sfact = M_TWO 319 case(2, 4) 320 oep%socc = M_ONE 321 oep%sfact = M_ONE 322 case default ! cannot handle any other case 323 ASSERT(.false.) 324 end select 325 326 POP_SUB(xc_oep_SpinFactor) 327 end subroutine xc_oep_SpinFactor 328 329 330 ! --------------------------------------------------------- 331 subroutine xc_oep_AnalyzeEigen(oep, st, is) 332 type(xc_oep_t), intent(inout) :: oep 333 type(states_elec_t), intent(in) :: st 334 integer, intent(in) :: is 335 336 integer :: ist 337 FLOAT :: max_eigen 338 FLOAT, allocatable :: eigenval(:), occ(:) 339 340 PUSH_SUB(xc_oep_AnalyzeEigen) 341 342 SAFE_ALLOCATE(eigenval(1:st%nst)) 343 SAFE_ALLOCATE (occ(1:st%nst)) 344 eigenval = M_ZERO 345 occ = M_ZERO 346 347 do ist = st%st_start, st%st_end 348 eigenval(ist) = st%eigenval(ist, is) 349 occ(ist) = st%occ(ist, is) 350 end do 351 352#if defined(HAVE_MPI) 353 if(st%parallel_in_states) then 354 call MPI_Barrier(st%mpi_grp%comm, mpi_err) 355 do ist = 1, st%nst 356 call MPI_Bcast(eigenval(ist), 1, MPI_FLOAT, st%node(ist), st%mpi_grp%comm, mpi_err) 357 call MPI_Bcast(occ(ist), 1, MPI_FLOAT, st%node(ist), st%mpi_grp%comm, mpi_err) 358 end do 359 end if 360#endif 361 362 ! find out the top occupied state, to correct for the asymptotics 363 ! of the potential 364 max_eigen = CNST(-1e30) 365 do ist = 1, st%nst 366 if((occ(ist) > M_EPSILON).and.(eigenval(ist) > max_eigen)) then 367 max_eigen = eigenval(ist) 368 end if 369 end do 370 371 oep%eigen_n = 1 372 do ist = 1, st%nst 373 if(occ(ist) > M_EPSILON) then 374 ! criterion for degeneracy 375 if(abs(eigenval(ist)-max_eigen) <= CNST(1e-3)) then 376 oep%eigen_type(ist) = 2 377 else 378 oep%eigen_type(ist) = 1 379 oep%eigen_index(oep%eigen_n) = ist 380 oep%eigen_n = oep%eigen_n + 1 381 end if 382 else 383 oep%eigen_type(ist) = 0 384 end if 385 end do 386 oep%eigen_n = oep%eigen_n - 1 387 388 ! find how many states are occupied. 389 oep%noccst = 0 390 do ist = 1, st%nst 391 if(st%occ(ist, is) > M_EPSILON) oep%noccst = ist 392 end do 393 394 SAFE_DEALLOCATE_A(eigenval) 395 SAFE_DEALLOCATE_A(occ) 396 POP_SUB(xc_oep_AnalyzeEigen) 397 end subroutine xc_oep_AnalyzeEigen 398 399 400#include "xc_kli_pauli_inc.F90" 401 402#include "undef.F90" 403#include "real.F90" 404#include "xc_kli_inc.F90" 405#include "xc_oep_x_inc.F90" 406#include "xc_oep_sic_inc.F90" 407#include "xc_oep_inc.F90" 408#include "xc_oep_qed_inc.F90" 409 410#include "undef.F90" 411#include "complex.F90" 412#include "xc_kli_inc.F90" 413#include "xc_oep_x_inc.F90" 414#include "xc_oep_sic_inc.F90" 415#include "xc_oep_inc.F90" 416#include "xc_oep_qed_inc.F90" 417 418end module xc_oep_oct_m 419 420!! Local Variables: 421!! mode: f90 422!! coding: utf-8 423!! End: 424