1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module unocc_oct_m 22 use density_oct_m 23 use eigensolver_oct_m 24 use energy_calc_oct_m 25 use global_oct_m 26 use output_oct_m 27 use hamiltonian_elec_oct_m 28 use io_oct_m 29 use kpoints_oct_m 30 use lcao_oct_m 31 use lda_u_oct_m 32 use lda_u_io_oct_m 33 use loct_oct_m 34 use mesh_oct_m 35 use messages_oct_m 36 use mpi_oct_m 37 use namespace_oct_m 38 use parser_oct_m 39 use profiling_oct_m 40 use restart_oct_m 41 use scf_oct_m 42 use simul_box_oct_m 43 use states_abst_oct_m 44 use states_elec_oct_m 45 use states_elec_io_oct_m 46 use states_elec_restart_oct_m 47 use system_oct_m 48 use v_ks_oct_m 49 use xc_oct_m 50 51 implicit none 52 53 private 54 public :: & 55 unocc_run 56 57 58contains 59 60 ! --------------------------------------------------------- 61 subroutine unocc_run(sys, fromscratch) 62 type(system_t), intent(inout) :: sys 63 logical, intent(inout) :: fromscratch 64 65 type(eigensolver_t) :: eigens 66 integer :: iunit, ierr, iter, ierr_rho, ik 67 logical :: read_gs, converged, forced_finish, showoccstates, is_orbital_dependent, occ_missing 68 integer :: max_iter, nst_calculated, showstart 69 integer :: n_filled, n_partially_filled, n_half_filled 70 integer, allocatable :: lowest_missing(:, :), occ_states(:) 71 character(len=10) :: dirname 72 type(restart_t) :: restart_load_unocc, restart_load_gs, restart_dump 73 logical :: write_density, bandstructure_mode, read_td_states 74 75 PUSH_SUB(unocc_run) 76 77 if (sys%hm%pcm%run_pcm) then 78 call messages_not_implemented("PCM for CalculationMode /= gs or td") 79 end if 80 81 !%Variable MaximumIter 82 !%Type integer 83 !%Default 50 84 !%Section Calculation Modes::Unoccupied States 85 !%Description 86 !% Maximum number of eigensolver iterations. The code will stop even if convergence 87 !% has not been achieved. -1 means unlimited. 0 means just do LCAO or read from 88 !% restart, and stop. 89 !%End 90 call parse_variable(sys%namespace, 'MaximumIter', 50, max_iter) 91 call messages_obsolete_variable(sys%namespace, 'UnoccMaximumIter', 'MaximumIter') 92 if(max_iter < 0) max_iter = huge(max_iter) 93 94 !%Variable UnoccShowOccStates 95 !%Type logical 96 !%Default false 97 !%Section Calculation Modes::Unoccupied States 98 !%Description 99 !% If true, the convergence for the occupied states will be shown too in the output. 100 !% This is useful for testing, or if the occupied states fail to converge. 101 !% It will be enabled automatically if only occupied states are being calculated. 102 !%End 103 call parse_variable(sys%namespace, 'UnoccShowOccStates', .false., showoccstates) 104 105 bandstructure_mode = .false. 106 if(simul_box_is_periodic(sys%gr%sb) .and. kpoints_get_kpoint_method(sys%gr%sb%kpoints) == KPOINTS_PATH) then 107 bandstructure_mode = .true. 108 end if 109 110 call init_(sys%gr%mesh, sys%st) 111 converged = .false. 112 113 read_td_states = .false. 114 if(bandstructure_mode) then 115 !%Variable UnoccUseTD 116 !%Type logical 117 !%Default no 118 !%Section Calculation Modes::Unoccupied States 119 !%Description 120 !% If true, Octopus will use the density and states from the restart/td folder to compute 121 !% the bandstructure, instead of the restart/gs ones. 122 !%End 123 call parse_variable(sys%namespace, 'UnoccUseTD', .false., read_td_states) 124 end if 125 126 127 SAFE_ALLOCATE(lowest_missing(1:sys%st%d%dim, 1:sys%st%d%nik)) 128 ! if there is no restart info to read, this will not get set otherwise 129 ! setting to zero means everything is missing. 130 lowest_missing(:,:) = 0 131 132 read_gs = .true. 133 if (.not. fromScratch) then 134 call restart_init(restart_load_unocc, sys%namespace, RESTART_UNOCC, RESTART_TYPE_LOAD, sys%mc, ierr, & 135 mesh = sys%gr%mesh, exact = .true.) 136 137 if(ierr == 0) then 138 call states_elec_load(restart_load_unocc, sys%namespace, sys%st, sys%gr, ierr, lowest_missing = lowest_missing) 139 call restart_end(restart_load_unocc) 140 end if 141 142 ! If successfully read states from unocc, do not read from gs. 143 ! If RESTART_GS and RESTART_UNOCC have the same directory (the default), and we tried RESTART_UNOCC 144 ! already and failed, it is a waste of time to try to read again. 145 if(ierr == 0 .or. restart_are_basedirs_equal(RESTART_GS, RESTART_UNOCC)) & 146 read_gs = .false. 147 end if 148 149 if(read_td_states) then 150 call restart_init(restart_load_gs, sys%namespace, RESTART_TD, RESTART_TYPE_LOAD, sys%mc, ierr_rho, mesh=sys%gr%mesh, & 151 exact=.true.) 152 else 153 call restart_init(restart_load_gs, sys%namespace, RESTART_GS, RESTART_TYPE_LOAD, sys%mc, ierr_rho, mesh=sys%gr%mesh, & 154 exact=.true.) 155 end if 156 157 if(ierr_rho == 0) then 158 if (read_gs) then 159 call states_elec_load(restart_load_gs, sys%namespace, sys%st, sys%gr, ierr, lowest_missing = lowest_missing) 160 end if 161 if(sys%hm%lda_u_level /= DFT_U_NONE) & 162 call lda_u_load(restart_load_gs, sys%hm%lda_u, sys%st, sys%hm%energy%dft_u, ierr) 163 call states_elec_load_rho(restart_load_gs, sys%st, sys%gr, ierr_rho) 164 write_density = restart_has_map(restart_load_gs) 165 call restart_end(restart_load_gs) 166 else 167 write_density = .true. 168 end if 169 170 SAFE_ALLOCATE(occ_states(1:sys%st%d%nik)) 171 do ik = 1, sys%st%d%nik 172 call occupied_states(sys%st, sys%namespace, ik, n_filled, n_partially_filled, n_half_filled) 173 occ_states(ik) = n_filled + n_partially_filled + n_half_filled 174 end do 175 176 is_orbital_dependent = (sys%ks%theory_level == HARTREE .or. sys%ks%theory_level == HARTREE_FOCK .or. & 177 (sys%ks%theory_level == KOHN_SHAM_DFT .and. xc_is_orbital_dependent(sys%ks%xc))) 178 179 if(is_orbital_dependent) then 180 message(1) = "Be sure your gs run is well converged since you have an orbital-dependent functional." 181 message(2) = "Otherwise, the occupied states may change in CalculationMode = unocc, and your" 182 message(3) = "unoccupied states will not be consistent with the gs run." 183 call messages_warning(3) 184 end if 185 186 if(ierr_rho /= 0 .or. is_orbital_dependent) then 187 occ_missing = .false. 188 do ik = 1, sys%st%d%nik 189 if(any(lowest_missing(1:sys%st%d%dim, ik) <= occ_states(ik))) then 190 occ_missing = .true. 191 end if 192 end do 193 194 if(occ_missing) then 195 if(is_orbital_dependent) then 196 message(1) = "For an orbital-dependent functional, all occupied orbitals must be provided." 197 else if(ierr_rho /= 0) then 198 message(1) = "Since density could not be read, all occupied orbitals must be provided." 199 end if 200 201 message(2) = "Not all the occupied orbitals could be read." 202 message(3) = "Please run a ground-state calculation first!" 203 call messages_fatal(3, only_root_writes = .true.) 204 end if 205 206 message(1) = "Unable to read density: Building density from wavefunctions." 207 call messages_info(1) 208 209 call density_calc(sys%st, sys%gr, sys%st%rho) 210 end if 211 212 call scf_state_info(sys%st) 213 214 if(fromScratch .or. ierr /= 0) then 215 if(fromScratch) then 216 ! do not use previously calculated occupied states 217 nst_calculated = min(maxval(occ_states), minval(lowest_missing) - 1) 218 else 219 ! or, use as many states as have been calculated 220 nst_calculated = minval(lowest_missing) - 1 221 end if 222 showstart = max(nst_calculated + 1, 1) 223 call lcao_run(sys, st_start = showstart) 224 else 225 ! we successfully read all the states and are planning to use them, no need for LCAO 226 call v_ks_calc(sys%ks, sys%namespace, sys%hm, sys%st, sys%geo, calc_eigenval = .false.) 227 showstart = minval(occ_states(:)) + 1 228 end if 229 230 231 232 ! it is strange and useless to see no eigenvalues written if you are only calculating 233 ! occupied states, on a different k-point. 234 if(showstart > sys%st%nst) showstart = 1 235 236 SAFE_DEALLOCATE_A(lowest_missing) 237 238 if(showoccstates) showstart = 1 239 240 ! In the case of someone using KPointsPath, the code assume that this is only for plotting a 241 ! bandstructure. This mode ensure that no restart information will be written for the new grid 242 if(bandstructure_mode) then 243 message(1) = "Info: The code will run in band structure mode." 244 message(2) = " No restart information will be printed." 245 call messages_info(2) 246 end if 247 248 if(.not. bandstructure_mode) then 249 ! Restart dump should be initialized after restart_load, as the mesh might have changed 250 call restart_init(restart_dump, sys%namespace, RESTART_UNOCC, RESTART_TYPE_DUMP, sys%mc, ierr, mesh=sys%gr%mesh) 251 252 ! make sure the density is defined on the same mesh as the wavefunctions that will be written 253 if(write_density) & 254 call states_elec_dump_rho(restart_dump, sys%st, sys%gr, ierr_rho) 255 end if 256 257 message(1) = "Info: Starting calculation of unoccupied states." 258 call messages_info(1) 259 260 ! reset this variable, so that the eigensolver passes through all states 261 eigens%converged(:) = 0 262 263 ! If not all gs wavefunctions were read when starting, in particular for nscf with different k-points, 264 ! the occupations must be recalculated each time, though they do not affect the result of course. 265 ! FIXME: This is wrong for metals where we must use the Fermi level from the original calculation! 266 call states_elec_fermi(sys%st, sys%namespace, sys%gr%mesh) 267 268 if(sys%st%d%pack_states .and. hamiltonian_elec_apply_packed(sys%hm)) call sys%st%pack() 269 270 do iter = 1, max_iter 271 call eigensolver_run(eigens, sys%namespace, sys%gr, sys%st, sys%hm, 1, converged, sys%st%nst_conv) 272 273 ! If not all gs wavefunctions were read when starting, in particular for nscf with different k-points, 274 ! the occupations must be recalculated each time, though they do not affect the result of course. 275 ! FIXME: This is wrong for metals where we must use the Fermi level from the original calculation! 276 call states_elec_fermi(sys%st, sys%namespace, sys%gr%mesh) 277 278 call write_iter_(sys%st) 279 280 ! write output file 281 if(mpi_grp_is_root(mpi_world)) then 282 call io_mkdir(STATIC_DIR, sys%namespace) 283 iunit = io_open(STATIC_DIR//'/eigenvalues', sys%namespace, action='write') 284 285 if(converged) then 286 write(iunit,'(a)') 'All states converged.' 287 else 288 write(iunit,'(a)') 'Some of the states are not fully converged!' 289 end if 290 write(iunit,'(a, e17.6)') 'Criterion = ', eigens%tolerance 291 write(iunit,'(1x)') 292 call states_elec_write_eigenvalues(iunit, sys%st%nst, sys%st, sys%gr%sb, eigens%diff) 293 call io_close(iunit) 294 end if 295 296 forced_finish = clean_stop(sys%mc%master_comm) 297 298 if(.not. bandstructure_mode) then 299 ! write restart information. 300 if(converged .or. (modulo(iter, sys%outp%restart_write_interval) == 0) & 301 .or. iter == max_iter .or. forced_finish) then 302 call states_elec_dump(restart_dump, sys%st, sys%gr, ierr, iter=iter) 303 if(ierr /= 0) then 304 message(1) = "Unable to write states wavefunctions." 305 call messages_warning(1) 306 end if 307 end if 308 end if 309 310 if(sys%outp%output_interval /= 0 .and. mod(iter, sys%outp%output_interval) == 0 & 311 .and. sys%outp%duringscf) then 312 write(dirname,'(a,i4.4)') "unocc.",iter 313 call output_all(sys%outp, sys%namespace, dirname, sys%gr, sys%geo, sys%st, sys%hm, sys%ks) 314 end if 315 316 if(converged .or. forced_finish) exit 317 318 end do 319 320 if(.not. bandstructure_mode) call restart_end(restart_dump) 321 322 if(sys%st%d%pack_states .and. hamiltonian_elec_apply_packed(sys%hm)) & 323 call sys%st%unpack() 324 325 if(any(eigens%converged(:) < occ_states(:))) then 326 write(message(1),'(a)') 'Some of the occupied states are not fully converged!' 327 call messages_warning(1) 328 end if 329 330 SAFE_DEALLOCATE_A(occ_states) 331 332 if(.not. converged) then 333 write(message(1),'(a)') 'Some of the unoccupied states are not fully converged!' 334 call messages_warning(1) 335 end if 336 337 if(simul_box_is_periodic(sys%gr%sb).and. sys%st%d%nik > sys%st%d%nspin) then 338 if(bitand(sys%gr%sb%kpoints%method, KPOINTS_PATH) /= 0) then 339 call states_elec_write_bandstructure(STATIC_DIR, sys%namespace, sys%st%nst, sys%st, & 340 sys%gr%sb, sys%geo, sys%gr%mesh, & 341 sys%hm%hm_base%phase, vec_pot = sys%hm%hm_base%uniform_vector_potential, & 342 vec_pot_var = sys%hm%hm_base%vector_potential) 343 end if 344 end if 345 346 347 call output_all(sys%outp, sys%namespace, STATIC_DIR, sys%gr, sys%geo, sys%st, sys%hm, sys%ks) 348 349 call end_() 350 POP_SUB(unocc_run) 351 352 contains 353 354 ! --------------------------------------------------------- 355 subroutine init_(mesh, st) 356 type(mesh_t), intent(in) :: mesh 357 type(states_elec_t), intent(inout) :: st 358 359 PUSH_SUB(unocc_run.init_) 360 361 call messages_obsolete_variable(sys%namespace, "NumberUnoccStates", "ExtraStates") 362 363 call states_elec_allocate_wfns(st, mesh, packed=.true.) 364 365 ! now the eigensolver stuff 366 call eigensolver_init(eigens, sys%namespace, sys%gr, st, sys%geo, sys%mc) 367 368 if(eigens%es_type == RS_RMMDIIS) then 369 message(1) = "With the RMMDIIS eigensolver for unocc, you will need to stop the calculation" 370 message(2) = "by hand, since the highest states will probably never converge." 371 call messages_warning(2) 372 end if 373 374 POP_SUB(unocc_run.init_) 375 end subroutine init_ 376 377 378 ! --------------------------------------------------------- 379 subroutine end_() 380 PUSH_SUB(unocc_run.end_) 381 382 call eigensolver_end(eigens, sys%gr) 383 384 POP_SUB(unocc_run.end_) 385 end subroutine end_ 386 387 ! --------------------------------------------------------- 388 subroutine write_iter_(st) 389 type(states_elec_t), intent(in) :: st 390 391 character(len=50) :: str 392 393 PUSH_SUB(unocc_run.write_iter_) 394 395 write(str, '(a,i5)') 'Unoccupied states iteration #', iter 396 call messages_print_stress(stdout, trim(str)) 397 398 write(message(1),'(a,i6,a,i6)') 'Converged states: ', minval(eigens%converged(1:st%d%nik)) 399 call messages_info(1) 400 401 call states_elec_write_eigenvalues(stdout, sys%st%nst, sys%st, sys%gr%sb, eigens%diff, st_start = showstart, compact = .true.) 402 403 call scf_print_mem_use() 404 405 call messages_print_stress(stdout) 406 407 POP_SUB(unocc_run.write_iter_) 408 end subroutine write_iter_ 409 410 411 end subroutine unocc_run 412 413 414end module unocc_oct_m 415 416!! Local Variables: 417!! mode: f90 418!! coding: utf-8 419!! End: 420