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