Octopus
states_elec.F90
Go to the documentation of this file.
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#include "global.h"
19
21 use accel_oct_m
22 use batch_oct_m
26 use batch_oct_m
29 use comm_oct_m
30 use debug_oct_m
34 use global_oct_m
35 use grid_oct_m
36 use io_oct_m
37 use, intrinsic :: iso_fortran_env
40 use loct_oct_m
41 use math_oct_m
42 use mesh_oct_m
47 use mpi_oct_m
50#ifdef HAVE_OPENMP
51 use omp_lib
52#endif
53 use parser_oct_m
57 use smear_oct_m
58 use space_oct_m
62 use types_oct_m
63 use unit_oct_m
67
68 implicit none
69
70 private
71
72 public :: &
103 stress_t, &
105
106
107 ! this type must be moved to stress module but due to circular dependency it is not possible now
108 type stress_t
109 real(real64) :: total(3,3) = m_zero
110 real(real64) :: kinetic(3,3) = m_zero
111 real(real64) :: Hartree(3,3) = m_zero
112 real(real64) :: xc(3,3) = m_zero
113 real(real64) :: ps_local(3,3) = m_zero
114 real(real64) :: ps_nl(3,3) = m_zero
115 real(real64) :: ion_ion(3,3) = m_zero
116 real(real64) :: vdw(3,3) = m_zero
117 real(real64) :: hubbard(3,3) = m_zero
118
119 real(real64) :: kinetic_sumrule
120 real(real64) :: hartree_sumrule
121 end type stress_t
122
123 ! TODO(Alex) Issue #672 Decouple k-point info from `states_elec_dim_t`
124
132 !
133 type, extends(states_abst_t) :: states_elec_t
134 ! Components are public by default
135 type(states_elec_dim_t) :: d
136 integer :: nst_conv
137
138 logical :: only_userdef_istates
139
140 type(states_elec_group_t) :: group
141 integer :: block_size
145 logical :: pack_states
149
150
151 character(len=1024), allocatable :: user_def_states(:,:,:)
154
155 ! TODO(Alex) Issue #821. Collate current quantities into an object.
156 ! the densities and currents (after all we are doing DFT :)
157 real(real64), allocatable :: rho(:,:)
158 real(real64), allocatable :: rho_core(:)
159
160 real(real64), allocatable :: current(:, :, :)
161 real(real64), allocatable :: current_para(:, :, :)
162 real(real64), allocatable :: current_dia(:, :, :)
163 real(real64), allocatable :: current_mag(:, :, :)
164 real(real64), allocatable :: current_kpt(:,:,:)
165
166 ! TODO(Alex) Issue #673. Create frozen density class and replace in states_elec_t
167 ! It may be required to "freeze" the deepest orbitals during the evolution; the density
168 ! of these orbitals is kept in frozen_rho. It is different from rho_core.
169 real(real64), allocatable :: frozen_rho(:, :)
170 real(real64), allocatable :: frozen_tau(:, :)
171 real(real64), allocatable :: frozen_gdens(:,:,:)
172 real(real64), allocatable :: frozen_ldens(:,:)
173
174 logical :: uniform_occ
175
176 real(real64), allocatable :: eigenval(:,:)
177 logical :: fixed_occ
178 logical :: restart_fixed_occ
179 logical :: restart_reorder_occs
180 real(real64), allocatable :: occ(:,:)
181 real(real64), allocatable :: kweights(:)
182 integer :: nik
183
184 logical, private :: fixed_spins
186 real(real64), allocatable :: spin(:, :, :)
187
188 real(real64) :: qtot
189 real(real64) :: val_charge
190
191 type(stress_t) :: stress_tensors
192
193 logical :: fromScratch
194 type(smear_t) :: smear
195
196 ! TODO(Alex) Issue #823 Move modelmbparticles out of states_elec_t
197 type(modelmb_particle_t) :: modelmbparticles
198 integer, allocatable :: mmb_nspindown(:,:)
199 integer, allocatable :: mmb_iyoung(:,:)
200 real(real64), allocatable :: mmb_proj(:)
202 logical :: parallel_in_states = .false.
204 ! TODO(Alex) Issue #824. Package the communicators in a single instance prior to removing
205 ! or consider creating a distributed_t instance for each (as distributed_t contains the an instance of mpi_grp_t)
206 type(mpi_grp_t) :: mpi_grp
207 type(mpi_grp_t) :: dom_st_mpi_grp
208 type(mpi_grp_t) :: st_kpt_mpi_grp
209 type(mpi_grp_t) :: dom_st_kpt_mpi_grp
210 type(blacs_proc_grid_t) :: dom_st_proc_grid
211
212 type(distributed_t) :: dist
213 logical :: scalapack_compatible
214
215 ! TODO(Alex) Issue #820. Remove lnst, st_start, st_end and node, as they are all contained within dist
216 integer :: lnst
217 integer :: st_start, st_end
218 integer, allocatable :: node(:)
219
220 ! TODO(Alex) Issue #824. Either change from data to a method, or package with `st_kpt_mpi_grp`
221 integer, allocatable :: st_kpt_task(:,:)
222
223 type(multicomm_all_pairs_t), private :: ap
224 logical :: symmetrize_density
225 integer :: randomization
226 integer :: orth_method = 0
227
228 real(real64) :: cl_states_mem
230 contains
231 procedure :: nullify => states_elec_null
232 procedure :: write_info => states_elec_write_info
233 procedure :: pack => states_elec_pack
234 procedure :: unpack => states_elec_unpack
235 procedure :: set_zero => states_elec_set_zero
236 procedure :: dipole => states_elec_calculate_dipole
237 end type states_elec_t
240 integer, public, parameter :: &
241 par_independent = 1, &
242 par_dependent = 2
243
245 interface states_elec_get_state
248 end interface states_elec_get_state
249
258
261 end interface
263contains
265 ! TODO(Alex): Issue #826. Rename to something like "states_elec_default_wfs_type", or remove
266 subroutine states_elec_null(st)
267 class(states_elec_t), intent(inout) :: st
268
271 st%wfs_type = type_float ! By default, calculations use real wavefunctions
273 st%packed = .false.
276 end subroutine states_elec_null
278
280 subroutine states_elec_init(st, namespace, space, valence_charge, kpoints)
281 type(states_elec_t), target, intent(inout) :: st
282 type(namespace_t), intent(in) :: namespace
283 type(electron_space_t), intent(in) :: space
284 real(real64), intent(in) :: valence_charge
285 type(kpoints_t), intent(in) :: kpoints
287 real(real64) :: excess_charge, nempty_percent
288 integer :: nempty, ntot, default
289 integer :: nempty_conv
290 logical :: force
291 real(real64), parameter :: tol = 1e-13_real64
293 push_sub_with_profile(states_elec_init)
294
295 st%fromScratch = .true. ! this will be reset if restart_read is called
296 call states_elec_null(st)
297
298 ! We get the spin dimension from the electronic space
299 ! TODO: Remove spin space information from states_elec_dim
300 st%d%ispin = space%ispin
302 ! Use of spinors requires complex wavefunctions.
303 if (st%d%ispin == spinors) call states_set_complex(st)
304
305 if (st%d%ispin /= unpolarized .and. kpoints%use_time_reversal) then
306 message(1) = "Time reversal symmetry is only implemented for unpolarized spins."
307 message(2) = "Use KPointsUseTimeReversal = no."
308 call messages_fatal(2, namespace=namespace)
309 end if
312 !%Variable ExcessCharge
313 !%Type float
314 !%Default 0.0
315 !%Section States
316 !%Description
317 !% The net charge of the system. A negative value means that we are adding
318 !% electrons, while a positive value means we are taking electrons
319 !% from the system.
320 !%End
321 call parse_variable(namespace, 'ExcessCharge', m_zero, excess_charge)
322
323 !%Variable TotalStates
324 !%Type integer
325 !%Default 0
326 !%Section States
327 !%Description
328 !% This variable sets the total number of states that Octopus will
329 !% use. This is normally not necessary since by default Octopus
330 !% sets the number of states to the minimum necessary to hold the
331 !% electrons present in the system. (This default behavior is
332 !% obtained by setting <tt>TotalStates</tt> to 0).
333 !%
334 !% If you want to add some unoccupied states, probably it is more convenient to use the variable
335 !% <tt>ExtraStates</tt>.
336 !%End
337 call parse_variable(namespace, 'TotalStates', 0, ntot)
338 if (ntot < 0) then
339 write(message(1), '(a,i5,a)') "Input: '", ntot, "' is not a valid value for TotalStates."
340 call messages_fatal(1, namespace=namespace)
341 end if
342
343 !%Variable ExtraStates
344 !%Type integer
345 !%Default 0
346 !%Section States
347 !%Description
348 !% The number of states is in principle calculated considering the minimum
349 !% numbers of states necessary to hold the electrons present in the system.
350 !% The number of electrons is
351 !% in turn calculated considering the nature of the species supplied in the
352 !% <tt>Species</tt> block, and the value of the <tt>ExcessCharge</tt> variable.
353 !% However, one may command <tt>Octopus</tt> to use more states, which is necessary if one wants to
354 !% use fractional occupational numbers, either fixed from the beginning through
355 !% the <tt>Occupations</tt> block or by prescribing
356 !% an electronic temperature with <tt>Smearing</tt>, or in order to calculate
357 !% excited states (including with <tt>CalculationMode = unocc</tt>).
358 !%End
359 call parse_variable(namespace, 'ExtraStates', 0, nempty)
360 if (nempty < 0) then
361 write(message(1), '(a,i5,a)') "Input: '", nempty, "' is not a valid value for ExtraStates."
362 message(2) = '(0 <= ExtraStates)'
363 call messages_fatal(2, namespace=namespace)
364 end if
365
366 if (ntot > 0 .and. nempty > 0) then
367 message(1) = 'You cannot set TotalStates and ExtraStates at the same time.'
368 call messages_fatal(1, namespace=namespace)
369 end if
370
371 !%Variable ExtraStatesInPercent
372 !%Type float
373 !%Default 0
374 !%Section States
375 !%Description
376 !% This variable allows to set the number of extra/empty states as percentage of the
377 !% used occupied states. For example, a value 35 for ExtraStatesInPercent would amount
378 !% to ceiling(35/100 * nstates) extra states, where nstates denotes the amount of occupied
379 !% states Octopus is using for the system at hand.
380 !%End
381 call parse_variable(namespace, 'ExtraStatesInPercent', m_zero, nempty_percent)
382 if (nempty_percent < 0) then
383 write(message(1), '(a,f8.6,a)') "Input: '", nempty_percent, &
384 "' should be a percentage value x (where x is parts in hundred) larger or equal 0"
385 call messages_fatal(1, namespace=namespace)
386 end if
387
388 if (nempty > 0 .and. nempty_percent > 0) then
389 message(1) = 'You cannot set ExtraStates and ExtraStatesInPercent at the same time.'
390 call messages_fatal(1, namespace=namespace)
391 end if
392
393 !%Variable ExtraStatesToConverge
394 !%Type integer
395 !%Default <tt>ExtraStates</tt> (Default 0)
396 !%Section States
397 !%Description
398 !% For <tt>gs</tt> and <tt>unocc</tt> calculations.
399 !% (For the <tt>gs</tt> calculation one needs to set <tt>ConvEigenError=yes</tt>)
400 !% Specifies the number of extra states that will be considered for reaching the convergence.
401 !% The calculation will consider the number off occupied states plus
402 !% <tt>ExtraStatesToConverge</tt> for the convergence criteria.
403 !% By default, all extra states need to be converged (For <tt>gs</tt> calculations only with <tt>ConvEigenError=yes</tt>).
404 !% Thus, together with <tt>ExtraStates</tt>, one can have some more states which will not be
405 !% considered for the convergence criteria, thus making the convergence of the
406 !% unocc calculation faster.
407 !%End
408 call parse_variable(namespace, 'ExtraStatesToConverge', nempty, nempty_conv)
409 if (nempty < 0) then
410 write(message(1), '(a,i5,a)') "Input: '", nempty_conv, "' is not a valid value for ExtraStatesToConverge."
411 message(2) = '(0 <= ExtraStatesToConverge)'
412 call messages_fatal(2, namespace=namespace)
413 end if
414
415 if (nempty_conv > nempty) then
416 message(1) = 'You cannot set ExtraStatesToConverge to a higher value than ExtraStates.'
417 call messages_fatal(1, namespace=namespace)
418 end if
419
420 ! For non-periodic systems this should just return the Gamma point
421 call states_elec_choose_kpoints(st, kpoints, namespace)
422
423 st%val_charge = valence_charge
424
425 st%qtot = -(st%val_charge + excess_charge)
426
427 if (st%qtot < -m_epsilon) then
428 write(message(1),'(a,f12.6,a)') 'Total charge = ', st%qtot, ' < 0'
429 message(2) = 'Check Species and ExcessCharge.'
430 call messages_fatal(2, only_root_writes = .true., namespace=namespace)
431 end if
432
433 select case (st%d%ispin)
434 case (unpolarized)
435 st%d%dim = 1
436 st%nst = nint(st%qtot/2)
437 if (st%nst*2 - st%qtot < -tol) st%nst = st%nst + 1
438 st%d%nspin = 1
439 st%d%spin_channels = 1
440 case (spin_polarized)
441 st%d%dim = 1
442 st%nst = nint(st%qtot/2)
443 if (st%nst*2 - st%qtot < -tol) st%nst = st%nst + 1
444 st%d%nspin = 2
445 st%d%spin_channels = 2
446 case (spinors)
447 st%d%dim = 2
448 st%nst = nint(st%qtot)
449 if (st%nst - st%qtot < -tol) st%nst = st%nst + 1
450 st%d%nspin = 4
451 st%d%spin_channels = 2
452 end select
453
454 if (ntot > 0) then
455 if (ntot < st%nst) then
456 message(1) = 'TotalStates is smaller than the number of states required by the system.'
457 call messages_fatal(1, namespace=namespace)
458 end if
459
460 st%nst = ntot
461 end if
462
463 if (nempty_percent > 0) then
464 nempty = ceiling(nempty_percent * st%nst / 100)
465 end if
466
467 st%nst_conv = st%nst + nempty_conv
468 st%nst = st%nst + nempty
469 if (st%nst == 0) then
470 message(1) = "Cannot run with number of states = zero."
471 call messages_fatal(1, namespace=namespace)
472 end if
473
474 !%Variable StatesBlockSize
475 !%Type integer
476 !%Section Execution::Optimization
477 !%Description
478 !% Some routines work over blocks of eigenfunctions, which
479 !% generally improves performance at the expense of increased
480 !% memory consumption. This variable selects the size of the
481 !% blocks to be used. If GPUs are used, the default is the
482 !% warp size (32 for NVIDIA, 32 or 64 for AMD);
483 !% otherwise it is 4.
484 !%End
485
486 if (accel_is_enabled()) then
487 ! Some AMD GPUs have a warp size of 64. When OpenCL is used
488 ! accel%warp_size = 1 which is why we use max(accel%warp_size, 32)
489 ! here so that StatesBlockSize is at least 32.
490 default = max(accel%warp_size, 32)
491 else
492 default = 4
493 end if
494
495 if (default > pad_pow2(st%nst)) default = pad_pow2(st%nst)
496
497 assert(default > 0)
498
499 call parse_variable(namespace, 'StatesBlockSize', default, st%block_size)
500 if (st%block_size < 1) then
501 call messages_write("The variable 'StatesBlockSize' must be greater than 0.")
502 call messages_fatal(namespace=namespace)
503 end if
504
505 st%block_size = min(st%block_size, st%nst)
506 conf%target_states_block_size = st%block_size
507
508 safe_allocate(st%eigenval(1:st%nst, 1:st%nik))
509 st%eigenval = huge(st%eigenval)
510
511 ! Periodic systems require complex wavefunctions
512 ! but not if it is Gamma-point only
513 if (.not. kpoints%gamma_only()) then
514 call states_set_complex(st)
515 end if
516
517 !%Variable OnlyUserDefinedInitialStates
518 !%Type logical
519 !%Default no
520 !%Section States
521 !%Description
522 !% If true, then only user-defined states from the block <tt>UserDefinedStates</tt>
523 !% will be used as initial states for a time-propagation. No attempt is made
524 !% to load ground-state orbitals from a previous ground-state run.
525 !%End
526 call parse_variable(namespace, 'OnlyUserDefinedInitialStates', .false., st%only_userdef_istates)
527
528 ! we now allocate some arrays
529 safe_allocate(st%occ (1:st%nst, 1:st%nik))
530 st%occ = m_zero
531 ! allocate space for formula strings that define user-defined states
532 if (parse_is_defined(namespace, 'UserDefinedStates') .or. parse_is_defined(namespace, 'OCTInitialUserdefined') &
533 .or. parse_is_defined(namespace, 'OCTTargetUserdefined')) then
534 safe_allocate(st%user_def_states(1:st%d%dim, 1:st%nst, 1:st%nik))
535 ! initially we mark all 'formulas' as undefined
536 st%user_def_states(1:st%d%dim, 1:st%nst, 1:st%nik) = 'undefined'
537 end if
538
539 if (st%d%ispin == spinors) then
540 safe_allocate(st%spin(1:3, 1:st%nst, 1:st%nik))
541 end if
542
543 !%Variable StatesRandomization
544 !%Type integer
545 !%Default par_independent
546 !%Section States
547 !%Description
548 !% The randomization of states can be done in two ways:
549 !% i) a parallelisation independent way (default), where the random states are identical,
550 !% irrespectively of the number of tasks and
551 !% ii) a parallelisation dependent way, which can prevent linear dependency
552 !% to occur for large systems.
553 !%Option par_independent 1
554 !% Parallelisation-independent randomization of states.
555 !%Option par_dependent 2
556 !% The randomization depends on the number of taks used in the calculation.
557 !%End
558 call parse_variable(namespace, 'StatesRandomization', par_independent, st%randomization)
559
560
561 call states_elec_read_initial_occs(st, namespace, excess_charge, kpoints)
562 call states_elec_read_initial_spins(st, namespace)
563
564 ! This test can only be done here, as smear_init is called inside states_elec_read_initial_occs, and
565 ! only there smear%photodop is set.
566
567 if (st%smear%photodop) then
568 if (nempty == 0) then
569 write(message(1), '(a,i5,a)') "PhotoDoping requires to specify ExtraStates."
570 message(2) = '(0 == ExtraStates)'
571 call messages_fatal(2, namespace=namespace)
572 end if
573 end if
574
575 st%st_start = 1
576 st%st_end = st%nst
577 st%lnst = st%nst
578 safe_allocate(st%node(1:st%nst))
579 st%node(1:st%nst) = 0
580
581 call mpi_grp_init(st%mpi_grp, mpi_comm_undefined)
582 st%parallel_in_states = .false.
583
584 call distributed_nullify(st%d%kpt, st%nik)
585
586 call modelmb_particles_init(st%modelmbparticles, namespace, space)
587 if (st%modelmbparticles%nparticle > 0) then
588 ! FIXME: check why this is not initialized properly in the test, or why it is written out when not initialized
589 safe_allocate(st%mmb_nspindown(1:st%modelmbparticles%ntype_of_particle, 1:st%nst))
590 st%mmb_nspindown(:,:) = -1
591 safe_allocate(st%mmb_iyoung(1:st%modelmbparticles%ntype_of_particle, 1:st%nst))
592 st%mmb_iyoung(:,:) = -1
593 safe_allocate(st%mmb_proj(1:st%nst))
594 st%mmb_proj(:) = m_zero
595 end if
596
597 !%Variable SymmetrizeDensity
598 !%Type logical
599 !%Default no
600 !%Section States
601 !%Description
602 !% When enabled the density is symmetrized. Currently, this can
603 !% only be done for periodic systems. (Experimental.)
604 !%End
605 call parse_variable(namespace, 'SymmetrizeDensity', kpoints%use_symmetries, st%symmetrize_density)
606 call messages_print_var_value('SymmetrizeDensity', st%symmetrize_density, namespace=namespace)
607
608 !%Variable ForceComplex
609 !%Type logical
610 !%Default no
611 !%Section Execution::Debug
612 !%Description
613 !% Normally <tt>Octopus</tt> determines automatically the type necessary
614 !% for the wavefunctions. When set to yes this variable will
615 !% force the use of complex wavefunctions.
616 !%
617 !% Warning: This variable is designed for testing and
618 !% benchmarking and normal users need not use it.
619 !%End
620 call parse_variable(namespace, 'ForceComplex', .false., force)
621
622 if (force) call states_set_complex(st)
623
624 st%packed = .false.
625
626 pop_sub_with_profile(states_elec_init)
627 end subroutine states_elec_init
628
629 ! ---------------------------------------------------------
632 !
633 subroutine states_elec_look(restart, nik, dim, nst, ierr)
634 type(restart_t), intent(in) :: restart
635 integer, intent(out) :: nik
636 integer, intent(out) :: dim
637 integer, intent(out) :: nst
638 integer, intent(out) :: ierr
639
640 character(len=256) :: lines(3)
641 character(len=20) :: char
642 integer :: iunit
643
644 push_sub(states_elec_look)
645
646 ierr = 0
647
648 iunit = restart_open(restart, 'states')
649 call restart_read(restart, iunit, lines, 3, ierr)
650 if (ierr == 0) then
651 read(lines(1), *) char, nst
652 read(lines(2), *) char, dim
653 read(lines(3), *) char, nik
654 end if
655 call restart_close(restart, iunit)
656
657 pop_sub(states_elec_look)
658 end subroutine states_elec_look
659
660 ! ---------------------------------------------------------
669 !
670 subroutine states_elec_read_initial_occs(st, namespace, excess_charge, kpoints)
671 type(states_elec_t), intent(inout) :: st
672 type(namespace_t), intent(in) :: namespace
673 real(real64), intent(in) :: excess_charge
674 type(kpoints_t), intent(in) :: kpoints
675
676 integer :: ik, ist, ispin, nspin, ncols, nrows, el_per_state, icol, start_pos, spin_n
677 type(block_t) :: blk
678 real(real64) :: rr, charge
679 logical :: integral_occs, unoccupied_states
680 real(real64), allocatable :: read_occs(:, :)
681 real(real64) :: charge_in_block
682
684
685 !%Variable RestartFixedOccupations
686 !%Type logical
687 !%Default yes
688 !%Section States
689 !%Description
690 !% Setting this variable will make the restart proceed as
691 !% if the occupations from the previous calculation had been set via the <tt>Occupations</tt> block,
692 !% <i>i.e.</i> fixed. Otherwise, occupations will be determined by smearing.
693 !%End
694 call parse_variable(namespace, 'RestartFixedOccupations', .true., st%restart_fixed_occ)
695 ! we will turn on st%fixed_occ if restart_read is ever called
696
697 !%Variable Occupations
698 !%Type block
699 !%Section States
700 !%Description
701 !% The occupation numbers of the orbitals can be fixed through the use of this
702 !% variable. For example:
703 !%
704 !% <tt>%Occupations
705 !% <br>&nbsp;&nbsp;2 | 2 | 2 | 2 | 2
706 !% <br>%</tt>
707 !%
708 !% would fix the occupations of the five states to 2. There can be
709 !% at most as many columns as states in the calculation. If there are fewer columns
710 !% than states, then the code will assume that the user is indicating the occupations
711 !% of the uppermost states where all lower states have full occupation (i.e. 2 for spin-unpolarized
712 !% calculations, 1 otherwise) and all higher states have zero occupation. The first column
713 !% will be taken to refer to the lowest state such that the occupations would be consistent
714 !% with the correct total charge. For example, if there are 8 electrons and 10 states (from
715 !% <tt>ExtraStates = 6</tt>), then an abbreviated specification
716 !%
717 !% <tt>%Occupations
718 !% <br>&nbsp;&nbsp;1 | 0 | 1
719 !% <br>%</tt>
720 !%
721 !% would be equivalent to a full specification
722 !%
723 !% <tt>%Occupations
724 !% <br>&nbsp;&nbsp;2 | 2 | 2 | 1 | 0 | 1 | 0 | 0 | 0 | 0
725 !% <br>%</tt>
726 !%
727 !% This is an example of use for constrained density-functional theory,
728 !% crudely emulating a HOMO->LUMO+1 optical excitation.
729 !% The number of rows should be equal
730 !% to the number of k-points times the number of spins. For example, for a finite system
731 !% with <tt>SpinComponents == spin_polarized</tt>,
732 !% this block should contain two lines, one for each spin channel.
733 !% All rows must have the same number of columns.
734 !%
735 !% The <tt>Occupations</tt> block is useful for the ground state of highly symmetric
736 !% small systems (like an open-shell atom), to fix the occupation numbers
737 !% of degenerate states in order to help <tt>octopus</tt> to converge. This is to
738 !% be used in conjuction with <tt>ExtraStates</tt>. For example, to calculate the
739 !% carbon atom, one would do:
740 !%
741 !% <tt>ExtraStates = 2
742 !% <br>%Occupations
743 !% <br>&nbsp;&nbsp;2 | 2/3 | 2/3 | 2/3
744 !% <br>%</tt>
745 !%
746 !% If you want the calculation to be spin-polarized (which makes more sense), you could do:
747 !%
748 !% <tt>ExtraStates = 2
749 !% <br>%Occupations
750 !% <br>&nbsp;&nbsp; 2/3 | 2/3 | 2/3
751 !% <br>&nbsp;&nbsp; 0 | 0 | 0
752 !% <br>%</tt>
753 !%
754 !% Note that in this case the first state is absent, the code will calculate four states
755 !% (two because there are four electrons, plus two because <tt>ExtraStates</tt> = 2), and since
756 !% it finds only three columns, it will occupy the first state with one electron for each
757 !% of the spin options.
758 !%
759 !% If the sum of occupations is not equal to the total charge set by <tt>ExcessCharge</tt>,
760 !% an error message is printed.
761 !% If <tt>FromScratch = no</tt> and <tt>RestartFixedOccupations = yes</tt>,
762 !% this block will be ignored.
763 !%End
764
765 integral_occs = .true.
766
767 occ_fix: if (parse_block(namespace, 'Occupations', blk) == 0) then
768 ! read in occupations
769 st%fixed_occ = .true.
770
771 ncols = parse_block_cols(blk, 0)
772 if (ncols > st%nst) then
773 call messages_input_error(namespace, "Occupations", "Too many columns in block Occupations.")
774 end if
775
776 nrows = parse_block_n(blk)
777 if (nrows /= st%nik) then
778 call messages_input_error(namespace, "Occupations", "Wrong number of rows in block Occupations.")
779 end if
780
781 do ik = 1, st%nik - 1
782 if (parse_block_cols(blk, ik) /= ncols) then
783 call messages_input_error(namespace, "Occupations", &
784 "All rows in block Occupations must have the same number of columns.")
785 end if
786 end do
787
788 ! Now we fill all the "missing" states with the maximum occupation.
789 if (st%d%ispin == unpolarized) then
790 el_per_state = 2
791 else
792 el_per_state = 1
793 end if
794
795 safe_allocate(read_occs(1:ncols, 1:st%nik))
796
797 charge_in_block = m_zero
798 do ik = 1, st%nik
799 do icol = 1, ncols
800 call parse_block_float(blk, ik - 1, icol - 1, read_occs(icol, ik))
801 charge_in_block = charge_in_block + read_occs(icol, ik) * st%kweights(ik)
802 end do
803 end do
804
805 spin_n = 2
806 select case (st%d%ispin)
807 case (unpolarized)
808 spin_n = 2
809 case (spin_polarized)
810 spin_n = 2
811 case (spinors)
812 spin_n = 1
813 end select
814
815 start_pos = nint((st%qtot - charge_in_block)/spin_n)
816
817 if (start_pos + ncols > st%nst) then
818 message(1) = "To balance charge, the first column in block Occupations is taken to refer to state"
819 write(message(2),'(a,i6,a)') "number ", start_pos, " but there are too many columns for the number of states."
820 write(message(3),'(a,i6,a)') "Solution: set ExtraStates = ", start_pos + ncols - st%nst
821 call messages_fatal(3, namespace=namespace)
822 end if
823
824 do ik = 1, st%nik
825 do ist = 1, start_pos
826 st%occ(ist, ik) = el_per_state
827 end do
828 end do
829
830 do ik = 1, st%nik
831 do ist = start_pos + 1, start_pos + ncols
832 st%occ(ist, ik) = read_occs(ist - start_pos, ik)
833 integral_occs = integral_occs .and. &
834 abs((st%occ(ist, ik) - el_per_state) * st%occ(ist, ik)) <= m_epsilon
835 end do
836 end do
837
838 do ik = 1, st%nik
839 do ist = start_pos + ncols + 1, st%nst
840 st%occ(ist, ik) = m_zero
841 end do
842 end do
843
844 call parse_block_end(blk)
845
846 safe_deallocate_a(read_occs)
847
848 else
849 st%fixed_occ = .false.
850 integral_occs = .false.
851
852 ! first guess for occupation...paramagnetic configuration
853 rr = m_one
854 if (st%d%ispin == unpolarized) rr = m_two
855
856 st%occ = m_zero
857 st%qtot = -(st%val_charge + excess_charge)
858
859 nspin = 1
860 if (st%d%nspin == 2) nspin = 2
861
862 do ik = 1, st%nik, nspin
863 charge = m_zero
864 do ist = 1, st%nst
865 do ispin = ik, ik + nspin - 1
866 st%occ(ist, ispin) = min(rr, -(st%val_charge + excess_charge) - charge)
867 charge = charge + st%occ(ist, ispin)
868 end do
869 end do
870 end do
871
872 end if occ_fix
873
874 !%Variable RestartReorderOccs
875 !%Type logical
876 !%Default no
877 !%Section States
878 !%Description
879 !% Consider doing a ground-state calculation, and then restarting with new occupations set
880 !% with the <tt>Occupations</tt> block, in an attempt to populate the orbitals of the original
881 !% calculation. However, the eigenvalues may reorder as the density changes, in which case the
882 !% occupations will now be referring to different orbitals. Setting this variable to yes will
883 !% try to solve this issue when the restart data is being read, by reordering the occupations
884 !% according to the order of the expectation values of the restart wavefunctions.
885 !%End
886 if (st%fixed_occ) then
887 call parse_variable(namespace, 'RestartReorderOccs', .false., st%restart_reorder_occs)
888 else
889 st%restart_reorder_occs = .false.
890 end if
891
892 call smear_init(st%smear, namespace, st%d%ispin, st%fixed_occ, integral_occs, kpoints)
893
894 unoccupied_states = (st%d%ispin /= spinors .and. st%nst*2 > st%qtot) .or. (st%d%ispin == spinors .and. st%nst > st%qtot)
895
896 if (.not. smear_is_semiconducting(st%smear) .and. .not. st%smear%method == smear_fixed_occ) then
897 if (.not. unoccupied_states) then
898 call messages_write('Smearing needs unoccupied states (via ExtraStates or TotalStates) to be useful.')
899 call messages_warning(namespace=namespace)
900 end if
901 end if
902
903 ! sanity check
904 charge = m_zero
905 do ist = 1, st%nst
906 charge = charge + sum(st%occ(ist, 1:st%nik) * st%kweights(1:st%nik))
907 end do
908 if (abs(charge - st%qtot) > 1e-6_real64) then
909 message(1) = "Initial occupations do not integrate to total charge."
910 write(message(2), '(6x,f12.6,a,f12.6)') charge, ' != ', st%qtot
911 call messages_fatal(2, only_root_writes = .true., namespace=namespace)
912 end if
913
914 st%uniform_occ = smear_is_semiconducting(st%smear) .and. .not. unoccupied_states
915
917 end subroutine states_elec_read_initial_occs
918
919
920 ! ---------------------------------------------------------
927 !
928 subroutine states_elec_read_initial_spins(st, namespace)
929 type(states_elec_t), intent(inout) :: st
930 type(namespace_t), intent(in) :: namespace
931
932 integer :: i, j
933 type(block_t) :: blk
934
936
937 st%fixed_spins = .false.
938 if (st%d%ispin /= spinors) then
940 return
941 end if
942
943 !%Variable InitialSpins
944 !%Type block
945 !%Section States
946 !%Description
947 !% The spin character of the initial random guesses for the spinors can
948 !% be fixed by making use of this block. Note that this will not "fix" the
949 !% the spins during the calculation (this cannot be done in spinors mode, in
950 !% being able to change the spins is why the spinors mode exists in the first
951 !% place).
952 !%
953 !% This block is meaningless and ignored if the run is not in spinors mode
954 !% (<tt>SpinComponents = spinors</tt>).
955 !%
956 !% The structure of the block is very simple: each column contains the desired
957 !% <math>\left< S_x \right>, \left< S_y \right>, \left< S_z \right> </math> for each spinor.
958 !% If the calculation is for a periodic system
959 !% and there is more than one <i>k</i>-point, the spins of all the <i>k</i>-points are
960 !% the same.
961 !%
962 !% For example, if we have two spinors, and we want one in the <math>S_x</math> "down" state,
963 !% and another one in the <math>S_x</math> "up" state:
964 !%
965 !% <tt>%InitialSpins
966 !% <br>&nbsp;&nbsp;&nbsp; 0.5 | 0.0 | 0.0
967 !% <br>&nbsp;&nbsp; -0.5 | 0.0 | 0.0
968 !% <br>%</tt>
969 !%
970 !% WARNING: if the calculation is for a system described by pseudopotentials (as
971 !% opposed to user-defined potentials or model systems), this option is
972 !% meaningless since the random spinors are overwritten by the atomic orbitals.
973 !%
974 !% This constraint must be fulfilled:
975 !% <br><math> \left< S_x \right>^2 + \left< S_y \right>^2 + \left< S_z \right>^2 = \frac{1}{4} </math>
976 !%End
977 spin_fix: if (parse_block(namespace, 'InitialSpins', blk) == 0) then
978 do i = 1, st%nst
979 do j = 1, 3
980 call parse_block_float(blk, i-1, j-1, st%spin(j, i, 1))
981 end do
982 if (abs(sum(st%spin(1:3, i, 1)**2) - m_fourth) > 1.0e-6_real64) call messages_input_error(namespace, 'InitialSpins')
983 end do
984 call parse_block_end(blk)
985 st%fixed_spins = .true.
986 do i = 2, st%nik
987 st%spin(:, :, i) = st%spin(:, :, 1)
988 end do
989 end if spin_fix
990
992 end subroutine states_elec_read_initial_spins
993
994
995 ! ---------------------------------------------------------
997 !
998 subroutine states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
999 type(states_elec_t), intent(inout) :: st
1000 class(mesh_t), intent(in) :: mesh
1001 type(type_t), optional, intent(in) :: wfs_type
1002 logical, optional, intent(in) :: skip(:)
1003 logical, optional, intent(in) :: packed
1004
1006
1007 if (present(wfs_type)) then
1008 assert(wfs_type == type_float .or. wfs_type == type_cmplx)
1009 st%wfs_type = wfs_type
1010 end if
1011
1012 call states_elec_init_block(st, mesh, skip = skip, packed=packed)
1013 call states_elec_set_zero(st)
1014
1016 end subroutine states_elec_allocate_wfns
1017
1018 !---------------------------------------------------------------------
1035 subroutine states_elec_init_block(st, mesh, verbose, skip, packed)
1036 type(states_elec_t), intent(inout) :: st
1037 type(mesh_t), intent(in) :: mesh
1038 logical, optional, intent(in) :: verbose
1039 logical, optional, intent(in) :: skip(:)
1040 logical, optional, intent(in) :: packed
1041
1042 integer :: ib, iqn, ist, istmin, istmax
1043 logical :: same_node, verbose_, packed_
1044 integer, allocatable :: bstart(:), bend(:)
1045
1046 push_sub(states_elec_init_block)
1047
1048 safe_allocate(bstart(1:st%nst))
1049 safe_allocate(bend(1:st%nst))
1050 safe_allocate(st%group%iblock(1:st%nst))
1051
1052 st%group%iblock = 0
1053
1054 verbose_ = optional_default(verbose, .true.)
1055 packed_ = optional_default(packed, .false.)
1056
1057 !In case we have a list of state to skip, we do not allocate them
1058 istmin = 1
1059 if (present(skip)) then
1060 do ist = 1, st%nst
1061 if (.not. skip(ist)) then
1062 istmin = ist
1063 exit
1064 end if
1065 end do
1066 end if
1067
1068 istmax = st%nst
1069 if (present(skip)) then
1070 do ist = st%nst, istmin, -1
1071 if (.not. skip(ist)) then
1072 istmax = ist
1073 exit
1074 end if
1075 end do
1076 end if
1077
1078 if (present(skip) .and. verbose_) then
1079 call messages_write('Info: Allocating states from ')
1080 call messages_write(istmin, fmt = 'i8')
1081 call messages_write(' to ')
1082 call messages_write(istmax, fmt = 'i8')
1083 call messages_info()
1084 end if
1085
1086 ! count and assign blocks
1087 ib = 0
1088 st%group%nblocks = 0
1089 bstart(1) = istmin
1090 do ist = istmin, istmax
1091 ib = ib + 1
1092
1093 st%group%iblock(ist) = st%group%nblocks + 1
1094
1095 same_node = .true.
1096 if (st%parallel_in_states .and. ist /= istmax) then
1097 ! We have to avoid that states that are in different nodes end
1098 ! up in the same block
1099 same_node = (st%node(ist + 1) == st%node(ist))
1100 end if
1101
1102 if (ib == st%block_size .or. ist == istmax .or. .not. same_node) then
1103 ib = 0
1104 st%group%nblocks = st%group%nblocks + 1
1105 bend(st%group%nblocks) = ist
1106 if (ist /= istmax) bstart(st%group%nblocks + 1) = ist + 1
1107 end if
1108 end do
1109
1110 safe_allocate(st%group%psib(1:st%group%nblocks, st%d%kpt%start:st%d%kpt%end))
1111
1112 safe_allocate(st%group%block_is_local(1:st%group%nblocks, 1:st%nik))
1113 st%group%block_is_local = .false.
1114 st%group%block_start = -1
1115 st%group%block_end = -2 ! this will make that loops block_start:block_end do not run if not initialized
1116
1117 do ib = 1, st%group%nblocks
1118 if (bstart(ib) >= st%st_start .and. bend(ib) <= st%st_end) then
1119 if (st%group%block_start == -1) st%group%block_start = ib
1120 st%group%block_end = ib
1121 do iqn = st%d%kpt%start, st%d%kpt%end
1122 st%group%block_is_local(ib, iqn) = .true.
1123
1124 if (states_are_real(st)) then
1125 call dwfs_elec_init(st%group%psib(ib, iqn), st%d%dim, bstart(ib), bend(ib), mesh%np_part, iqn, &
1126 special=.true., packed=packed_)
1127 else
1128 call zwfs_elec_init(st%group%psib(ib, iqn), st%d%dim, bstart(ib), bend(ib), mesh%np_part, iqn, &
1129 special=.true., packed=packed_)
1130 end if
1131
1132 end do
1133 end if
1134 end do
1135
1136 safe_allocate(st%group%block_range(1:st%group%nblocks, 1:2))
1137 safe_allocate(st%group%block_size(1:st%group%nblocks))
1138
1139 st%group%block_range(1:st%group%nblocks, 1) = bstart(1:st%group%nblocks)
1140 st%group%block_range(1:st%group%nblocks, 2) = bend(1:st%group%nblocks)
1141 st%group%block_size(1:st%group%nblocks) = bend(1:st%group%nblocks) - bstart(1:st%group%nblocks) + 1
1142
1143 st%group%block_initialized = .true.
1144
1145 safe_allocate(st%group%block_node(1:st%group%nblocks, 1:st%nik))
1146 st%group%block_node = 0
1147
1148 assert(allocated(st%node))
1149 assert(all(st%node >= 0) .and. all(st%node < st%mpi_grp%size))
1150
1151 do iqn = st%d%kpt%start, st%d%kpt%end
1152 do ib = st%group%block_start, st%group%block_end
1153 st%group%block_node(ib, iqn) = st%st_kpt_mpi_grp%rank
1154 end do
1155 end do
1156
1157 call comm_allreduce(st%st_kpt_mpi_grp, st%group%block_node)
1158
1159 if (verbose_) then
1160 call messages_write('Info: Blocks of states')
1161 call messages_info()
1162 do ib = 1, st%group%nblocks
1163 call messages_write(' Block ')
1164 call messages_write(ib, fmt = 'i8')
1165 call messages_write(' contains ')
1166 call messages_write(st%group%block_size(ib), fmt = 'i8')
1167 call messages_write(' states')
1168 if (st%group%block_size(ib) > 0) then
1169 call messages_write(':')
1170 call messages_write(st%group%block_range(ib, 1), fmt = 'i8')
1171 call messages_write(' - ')
1172 call messages_write(st%group%block_range(ib, 2), fmt = 'i8')
1173 end if
1174 call messages_info()
1175 end do
1176 end if
1177
1178!!$!!!!DEBUG
1179!!$ ! some debug output that I will keep here for the moment
1180!!$ if (mpi_grp_is_root(mpi_world)) then
1181!!$ print*, "NST ", st%nst
1182!!$ print*, "BLOCKSIZE ", st%block_size
1183!!$ print*, "NBLOCKS ", st%group%nblocks
1184!!$
1185!!$ print*, "==============="
1186!!$ do ist = 1, st%nst
1187!!$ print*, st%node(ist), ist, st%group%iblock(ist)
1188!!$ end do
1189!!$ print*, "==============="
1190!!$
1191!!$ do ib = 1, st%group%nblocks
1192!!$ print*, ib, bstart(ib), bend(ib)
1193!!$ end do
1194!!$
1195!!$ end if
1196!!$!!!!ENDOFDEBUG
1197
1198 safe_deallocate_a(bstart)
1199 safe_deallocate_a(bend)
1200 pop_sub(states_elec_init_block)
1201 end subroutine states_elec_init_block
1202
1203
1204 ! ---------------------------------------------------------
1206 subroutine states_elec_deallocate_wfns(st)
1207 type(states_elec_t), intent(inout) :: st
1208
1210
1211 call states_elec_group_end(st%group, st%d)
1212
1214 end subroutine states_elec_deallocate_wfns
1215
1216
1217 ! ---------------------------------------------------------
1218 subroutine states_elec_densities_init(st, gr)
1219 type(states_elec_t), target, intent(inout) :: st
1220 type(grid_t), intent(in) :: gr
1221
1222 real(real64) :: fsize
1223
1225
1226 safe_allocate(st%rho(1:gr%np_part, 1:st%d%nspin))
1227 st%rho = m_zero
1228
1229 fsize = gr%np_part*8.0_real64*st%block_size
1230
1231 call messages_write('Info: states-block size = ')
1232 call messages_write(fsize, fmt = '(f10.1)', align_left = .true., units = unit_megabytes, print_units = .true.)
1233 call messages_info()
1234
1236 end subroutine states_elec_densities_init
1237
1238 !---------------------------------------------------------------------
1239 subroutine states_elec_allocate_current(st, space, mesh)
1240 type(states_elec_t), intent(inout) :: st
1241 class(space_t), intent(in) :: space
1242 class(mesh_t), intent(in) :: mesh
1243
1245
1246 if (.not. allocated(st%current)) then
1247 safe_allocate(st%current(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1248 st%current = m_zero
1249 end if
1250
1251 if (.not. allocated(st%current_para)) then
1252 safe_allocate(st%current_para(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1253 st%current_para = m_zero
1254 end if
1255
1256 if (.not. allocated(st%current_dia)) then
1257 safe_allocate(st%current_dia(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1258 st%current_dia= m_zero
1259 end if
1260
1261 if (.not. allocated(st%current_mag)) then
1262 safe_allocate(st%current_mag(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1263 st%current_mag= m_zero
1264 end if
1265
1266 if (.not. allocated(st%current_kpt)) then
1267 safe_allocate(st%current_kpt(1:mesh%np,1:space%dim,st%d%kpt%start:st%d%kpt%end))
1268 st%current_kpt = m_zero
1269 end if
1270
1272 end subroutine states_elec_allocate_current
1273
1274 !---------------------------------------------------------------------
1282 subroutine states_elec_exec_init(st, namespace, mc)
1283 type(states_elec_t), intent(inout) :: st
1284 type(namespace_t), intent(in) :: namespace
1285 type(multicomm_t), intent(in) :: mc
1286
1287 integer :: default
1288
1289 push_sub(states_elec_exec_init)
1290
1291 !%Variable StatesPack
1292 !%Type logical
1293 !%Section Execution::Optimization
1294 !%Description
1295 !% When set to yes, states are stored in packed mode, which improves
1296 !% performance considerably. Not all parts of the code will profit from
1297 !% this, but should nevertheless work regardless of how the states are
1298 !% stored.
1300 !% If GPUs are used and this variable is set to yes, Octopus
1301 !% will store the wave-functions in device (GPU) memory. If
1302 !% there is not enough memory to store all the wave-functions,
1303 !% execution will stop with an error.
1304 !%
1305 !% See also the related <tt>HamiltonianApplyPacked</tt> variable.
1306 !%
1307 !% The default is yes.
1308 !%End
1309
1310 call parse_variable(namespace, 'StatesPack', .true., st%pack_states)
1312 call messages_print_var_value('StatesPack', st%pack_states, namespace=namespace)
1313
1314 call messages_obsolete_variable(namespace, 'StatesMirror')
1315
1316 !%Variable StatesOrthogonalization
1317 !%Type integer
1318 !%Section SCF::Eigensolver
1319 !%Description
1320 !% The full orthogonalization method used by some
1321 !% eigensolvers. The default is <tt>cholesky_serial</tt>, except with state
1322 !% parallelization, the default is <tt>cholesky_parallel</tt>.
1323 !%Option cholesky_serial 1
1324 !% Cholesky decomposition implemented using
1325 !% BLAS/LAPACK. Can be used with domain parallelization but not
1326 !% state parallelization.
1327 !%Option cholesky_parallel 2
1328 !% Cholesky decomposition implemented using
1329 !% ScaLAPACK. Compatible with states parallelization.
1330 !%Option cgs 3
1331 !% Classical Gram-Schmidt (CGS) orthogonalization.
1332 !% Can be used with domain parallelization but not state parallelization.
1333 !% The algorithm is defined in Giraud et al., Computers and Mathematics with Applications 50, 1069 (2005).
1334 !%Option mgs 4
1335 !% Modified Gram-Schmidt (MGS) orthogonalization.
1336 !% Can be used with domain parallelization but not state parallelization.
1337 !% The algorithm is defined in Giraud et al., Computers and Mathematics with Applications 50, 1069 (2005).
1338 !%Option drcgs 5
1339 !% Classical Gram-Schmidt orthogonalization with double-step reorthogonalization.
1340 !% Can be used with domain parallelization but not state parallelization.
1341 !% The algorithm is taken from Giraud et al., Computers and Mathematics with Applications 50, 1069 (2005).
1342 !% According to this reference, this is much more precise than CGS or MGS algorithms.
1343 !% The MGS version seems not to improve much the stability and would require more communications over the domains.
1344 !%End
1345
1346 default = option__statesorthogonalization__cholesky_serial
1347#ifdef HAVE_SCALAPACK
1349 default = option__statesorthogonalization__cholesky_parallel
1350 end if
1351#endif
1352
1353 call parse_variable(namespace, 'StatesOrthogonalization', default, st%orth_method)
1354
1355 if (.not. varinfo_valid_option('StatesOrthogonalization', st%orth_method)) then
1356 call messages_input_error(namespace, 'StatesOrthogonalization')
1357 end if
1358 call messages_print_var_option('StatesOrthogonalization', st%orth_method, namespace=namespace)
1359
1360 !%Variable StatesCLDeviceMemory
1361 !%Type float
1362 !%Section Execution::Optimization
1363 !%Default -512
1364 !%Description
1365 !% This variable selects the amount of OpenCL device memory that
1366 !% will be used by Octopus to store the states.
1367 !%
1368 !% A positive number smaller than 1 indicates a fraction of the total
1369 !% device memory. A number larger than one indicates an absolute
1370 !% amount of memory in megabytes. A negative number indicates an
1371 !% amount of memory in megabytes that would be subtracted from
1372 !% the total device memory.
1373 !%End
1374 call parse_variable(namespace, 'StatesCLDeviceMemory', -512.0_real64, st%cl_states_mem)
1376 pop_sub(states_elec_exec_init)
1377 end subroutine states_elec_exec_init
1378
1379
1380 ! ---------------------------------------------------------
1382 !
1383 subroutine states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
1384 type(states_elec_t), target, intent(inout) :: stout
1385 type(states_elec_t), intent(in) :: stin
1386 logical, optional, intent(in) :: exclude_wfns
1387 logical, optional, intent(in) :: exclude_eigenval
1388 logical, optional, intent(in) :: special
1389
1390 logical :: exclude_wfns_
1391
1392 push_sub(states_elec_copy)
1393
1394 exclude_wfns_ = optional_default(exclude_wfns, .false.)
1395
1396 call states_elec_null(stout)
1397
1398 call states_elec_dim_copy(stout%d, stin%d)
1399 safe_allocate_source_a(stout%kweights, stin%kweights)
1400 stout%nik = stin%nik
1401
1402 call modelmb_particles_copy(stout%modelmbparticles, stin%modelmbparticles)
1403 if (stin%modelmbparticles%nparticle > 0) then
1404 safe_allocate_source_a(stout%mmb_nspindown, stin%mmb_nspindown)
1405 safe_allocate_source_a(stout%mmb_iyoung, stin%mmb_iyoung)
1406 safe_allocate_source_a(stout%mmb_proj, stin%mmb_proj)
1407 end if
1408
1409 stout%wfs_type = stin%wfs_type
1410 stout%nst = stin%nst
1411 stout%block_size = stin%block_size
1412 stout%orth_method = stin%orth_method
1413
1414 stout%cl_states_mem = stin%cl_states_mem
1415 stout%pack_states = stin%pack_states
1416
1417
1418 stout%only_userdef_istates = stin%only_userdef_istates
1419
1420 if (.not. exclude_wfns_) then
1421 safe_allocate_source_a(stout%rho, stin%rho)
1422 end if
1423
1424 stout%uniform_occ = stin%uniform_occ
1425
1426 if (.not. optional_default(exclude_eigenval, .false.)) then
1427 safe_allocate_source_a(stout%eigenval, stin%eigenval)
1428 safe_allocate_source_a(stout%occ, stin%occ)
1429 safe_allocate_source_a(stout%spin, stin%spin)
1430 end if
1431
1432 ! the call to init_block is done at the end of this subroutine
1433 ! it allocates iblock, psib, block_is_local
1434 stout%group%nblocks = stin%group%nblocks
1435
1436 safe_allocate_source_a(stout%user_def_states, stin%user_def_states)
1437
1438 safe_allocate_source_a(stout%current, stin%current)
1439 safe_allocate_source_a(stout%current_kpt, stin%current_kpt)
1440 safe_allocate_source_a(stout%rho_core, stin%rho_core)
1441 safe_allocate_source_a(stout%frozen_rho, stin%frozen_rho)
1442 safe_allocate_source_a(stout%frozen_tau, stin%frozen_tau)
1443 safe_allocate_source_a(stout%frozen_gdens, stin%frozen_gdens)
1444 safe_allocate_source_a(stout%frozen_ldens, stin%frozen_ldens)
1445
1446 stout%fixed_occ = stin%fixed_occ
1447 stout%restart_fixed_occ = stin%restart_fixed_occ
1448
1449 stout%fixed_spins = stin%fixed_spins
1450
1451 stout%qtot = stin%qtot
1452 stout%val_charge = stin%val_charge
1453
1454 call smear_copy(stout%smear, stin%smear)
1455
1456 stout%parallel_in_states = stin%parallel_in_states
1457 call mpi_grp_copy(stout%mpi_grp, stin%mpi_grp)
1458 call mpi_grp_copy(stout%dom_st_kpt_mpi_grp, stin%dom_st_kpt_mpi_grp)
1459 call mpi_grp_copy(stout%st_kpt_mpi_grp, stin%st_kpt_mpi_grp)
1460 call mpi_grp_copy(stout%dom_st_mpi_grp, stin%dom_st_mpi_grp)
1461 safe_allocate_source_a(stout%node, stin%node)
1462 safe_allocate_source_a(stout%st_kpt_task, stin%st_kpt_task)
1463
1464#ifdef HAVE_SCALAPACK
1465 call blacs_proc_grid_copy(stin%dom_st_proc_grid, stout%dom_st_proc_grid)
1466#endif
1467
1468 call distributed_copy(stin%dist, stout%dist)
1469
1470 stout%scalapack_compatible = stin%scalapack_compatible
1471
1472 stout%lnst = stin%lnst
1473 stout%st_start = stin%st_start
1474 stout%st_end = stin%st_end
1475
1476 if (stin%parallel_in_states) call multicomm_all_pairs_copy(stout%ap, stin%ap)
1477
1478 stout%symmetrize_density = stin%symmetrize_density
1479
1480 if (.not. exclude_wfns_) call states_elec_group_copy(stin%d, stin%group, stout%group, special=special)
1481
1482 stout%packed = stin%packed
1483
1484 stout%randomization = stin%randomization
1485
1486 pop_sub(states_elec_copy)
1487 end subroutine states_elec_copy
1488
1489
1490 ! ---------------------------------------------------------
1492 !
1493 subroutine states_elec_end(st)
1494 type(states_elec_t), intent(inout) :: st
1495
1496 push_sub(states_elec_end)
1497
1498 call states_elec_dim_end(st%d)
1499
1500 if (st%modelmbparticles%nparticle > 0) then
1501 safe_deallocate_a(st%mmb_nspindown)
1502 safe_deallocate_a(st%mmb_iyoung)
1503 safe_deallocate_a(st%mmb_proj)
1504 end if
1505 call modelmb_particles_end(st%modelmbparticles)
1506
1507 ! this deallocates dpsi, zpsi, psib, iblock, iblock
1509
1510 safe_deallocate_a(st%user_def_states)
1511
1512 safe_deallocate_a(st%rho)
1513 safe_deallocate_a(st%eigenval)
1514
1515 safe_deallocate_a(st%current)
1516 safe_deallocate_a(st%current_para)
1517 safe_deallocate_a(st%current_dia)
1518 safe_deallocate_a(st%current_mag)
1519 safe_deallocate_a(st%current_kpt)
1520 safe_deallocate_a(st%rho_core)
1521 safe_deallocate_a(st%frozen_rho)
1522 safe_deallocate_a(st%frozen_tau)
1523 safe_deallocate_a(st%frozen_gdens)
1524 safe_deallocate_a(st%frozen_ldens)
1525 safe_deallocate_a(st%occ)
1526 safe_deallocate_a(st%spin)
1527 safe_deallocate_a(st%kweights)
1528
1529
1530#ifdef HAVE_SCALAPACK
1531 call blacs_proc_grid_end(st%dom_st_proc_grid)
1532#endif
1533
1534 call distributed_end(st%dist)
1535
1536 safe_deallocate_a(st%node)
1537 safe_deallocate_a(st%st_kpt_task)
1538
1539 if (st%parallel_in_states) then
1540 safe_deallocate_a(st%ap%schedule)
1541 end if
1542
1543 pop_sub(states_elec_end)
1544 end subroutine states_elec_end
1545
1546
1547 ! TODO(Alex) Issue #684. Abstract duplicate code in states_elec_generate_random, to get it to
1548 ! a point where it can be refactored.
1550 subroutine states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
1551 type(states_elec_t), intent(inout) :: st
1552 class(mesh_t), intent(in) :: mesh
1553 type(kpoints_t), intent(in) :: kpoints
1554 integer, optional, intent(in) :: ist_start_
1555 integer, optional, intent(in) :: ist_end_
1556 integer, optional, intent(in) :: ikpt_start_
1557 integer, optional, intent(in) :: ikpt_end_
1558 logical, optional, intent(in) :: normalized
1560
1561 integer :: ist, ik, id, ist_start, ist_end, jst, ikpt_start, ikpt_end
1562 complex(real64) :: alpha, beta
1563 real(real64), allocatable :: dpsi(:, :)
1564 complex(real64), allocatable :: zpsi(:, :), zpsi2(:)
1565 integer :: ikpoint, ip
1566 type(batch_t) :: ffb
1567
1568 logical :: normalized_
1569
1570 normalized_ = optional_default(normalized, .true.)
1571
1573
1574 ist_start = optional_default(ist_start_, 1)
1575 ist_end = optional_default(ist_end_, st%nst)
1576 ikpt_start = optional_default(ikpt_start_, 1)
1577 ikpt_end = optional_default(ikpt_end_, st%nik)
1578
1579 safe_allocate(dpsi(1:mesh%np, 1:st%d%dim))
1580 if (states_are_complex(st)) then
1581 safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
1582 end if
1583
1584 select case (st%d%ispin)
1587 do ik = ikpt_start, ikpt_end
1588 ikpoint = st%d%get_kpoint_index(ik)
1589 do ist = ist_start, ist_end
1590 if (states_are_real(st).or.kpoints_point_is_gamma(kpoints, ikpoint)) then
1591 if (st%randomization == par_independent) then
1592 call dmf_random(mesh, dpsi(:, 1), &
1593 pre_shift = mesh%pv%xlocal-1, &
1594 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1595 normalized = normalized)
1596 ! Ensures that the grid points are properly distributed in the domain parallel case
1597 if(mesh%parallel_in_domains) then
1598 call batch_init(ffb, dpsi(:,1))
1599 call dmesh_batch_exchange_points(mesh, ffb, backward_map = .true.)
1600 call ffb%end()
1601 end if
1602 else
1603 call dmf_random(mesh, dpsi(:, 1), normalized = normalized)
1604 end if
1605 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1606 if (states_are_complex(st)) then !Gamma point
1607 do ip = 1, mesh%np
1608 zpsi(ip,1) = cmplx(dpsi(ip,1), m_zero, real64)
1609 end do
1610 call states_elec_set_state(st, mesh, ist, ik, zpsi)
1611 else
1612 call states_elec_set_state(st, mesh, ist, ik, dpsi)
1613 end if
1614 else
1615 if (st%randomization == par_independent) then
1616 call zmf_random(mesh, zpsi(:, 1), &
1617 pre_shift = mesh%pv%xlocal-1, &
1618 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1619 normalized = normalized)
1620 ! Ensures that the grid points are properly distributed in the domain parallel case
1621 if(mesh%parallel_in_domains) then
1622 call batch_init(ffb, zpsi(:,1))
1623 call zmesh_batch_exchange_points(mesh, ffb, backward_map = .true.)
1624 call ffb%end()
1625 end if
1626 else
1627 call zmf_random(mesh, zpsi(:, 1), normalized = normalized)
1628 end if
1629 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1630 call states_elec_set_state(st, mesh, ist, ik, zpsi)
1631 end if
1632 end do
1633 end do
1634
1635 case (spinors)
1636
1637 assert(states_are_complex(st))
1638
1639 if (st%fixed_spins) then
1640
1641 do ik = ikpt_start, ikpt_end
1642 ikpoint = st%d%get_kpoint_index(ik)
1643 do ist = ist_start, ist_end
1644 if (kpoints_point_is_gamma(kpoints, ikpoint)) then
1645 if (st%randomization == par_independent) then
1646 call dmf_random(mesh, dpsi(:, 1), &
1647 pre_shift = mesh%pv%xlocal-1, &
1648 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1649 normalized = normalized)
1650 ! Ensures that the grid points are properly distributed in the domain parallel case
1651 if(mesh%parallel_in_domains) then
1652 call batch_init(ffb, dpsi(:,1))
1653 call dmesh_batch_exchange_points(mesh, ffb, backward_map = .true.)
1654 call ffb%end()
1655 end if
1656 else
1657 call dmf_random(mesh, dpsi(:, 1), normalized = normalized)
1658 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1659 end if
1660 do ip = 1, mesh%np
1661 zpsi(ip,1) = cmplx(dpsi(ip,1), m_zero, real64)
1662 end do
1663 call states_elec_set_state(st, mesh, ist, ik, zpsi)
1664 else
1665 if (st%randomization == par_independent) then
1666 call zmf_random(mesh, zpsi(:, 1), &
1667 pre_shift = mesh%pv%xlocal-1, &
1668 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1669 normalized = normalized)
1670 ! Ensures that the grid points are properly distributed in the domain parallel case
1671 if(mesh%parallel_in_domains) then
1672 call batch_init(ffb, zpsi(:,1))
1673 call zmesh_batch_exchange_points(mesh, ffb, backward_map = .true.)
1674 call ffb%end()
1675 end if
1676 else
1677 call zmf_random(mesh, zpsi(:, 1), normalized = normalized)
1678 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1679 end if
1680 end if
1681 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1682 ! In this case, the spinors are made of a spatial part times a vector [alpha beta]^T in
1683 ! spin space (i.e., same spatial part for each spin component). So (alpha, beta)
1684 ! determines the spin values. The values of (alpha, beta) can be be obtained
1685 ! with simple formulae from <Sx>, <Sy>, <Sz>.
1686 !
1687 ! Note that here we orthonormalize the orbital part. This ensures that the spinors
1688 ! are untouched later in the general orthonormalization, and therefore the spin values
1689 ! of each spinor remain the same.
1690 safe_allocate(zpsi2(1:mesh%np))
1691 do jst = ist_start, ist - 1
1692 call states_elec_get_state(st, mesh, 1, jst, ik, zpsi2)
1693 zpsi(1:mesh%np, 1) = zpsi(1:mesh%np, 1) - zmf_dotp(mesh, zpsi(:, 1), zpsi2)*zpsi2(1:mesh%np)
1694 end do
1695 safe_deallocate_a(zpsi2)
1696
1697 call zmf_normalize(mesh, 1, zpsi)
1698 zpsi(1:mesh%np, 2) = zpsi(1:mesh%np, 1)
1699
1700 alpha = cmplx(sqrt(m_half + st%spin(3, ist, ik)), m_zero, real64)
1701 beta = cmplx(sqrt(m_one - abs(alpha)**2), m_zero, real64)
1702 if (abs(alpha) > m_zero) then
1703 beta = cmplx(st%spin(1, ist, ik) / abs(alpha), st%spin(2, ist, ik) / abs(alpha), real64)
1704 end if
1705 zpsi(1:mesh%np, 1) = alpha*zpsi(1:mesh%np, 1)
1706 zpsi(1:mesh%np, 2) = beta*zpsi(1:mesh%np, 2)
1707 call states_elec_set_state(st, mesh, ist, ik, zpsi)
1708 end do
1709 end do
1710 else
1711 do ik = ikpt_start, ikpt_end
1712 do ist = ist_start, ist_end
1713 do id = 1, st%d%dim
1714 if (st%randomization == par_independent) then
1715 call zmf_random(mesh, zpsi(:, id), &
1716 pre_shift = mesh%pv%xlocal-1, &
1717 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1718 normalized = .false.)
1719 ! Ensures that the grid points are properly distributed in the domain parallel case
1720 if(mesh%parallel_in_domains) then
1721 call batch_init(ffb, zpsi(:, id))
1722 call zmesh_batch_exchange_points(mesh, ffb, backward_map = .true.)
1723 call ffb%end()
1724 end if
1725 else
1726 call zmf_random(mesh, zpsi(:, id), normalized = .false.)
1727 end if
1728 end do
1729 ! We need to generate the wave functions even if not using them in order to be consistent
1730 ! with the random seed in parallel runs.
1731 if (.not. state_kpt_is_local(st, ist, ik)) cycle
1732 ! Note that mf_random normalizes each spin channel independently to 1.
1733 ! Therefore we need to renormalize the spinor:
1734 if (normalized_) call zmf_normalize(mesh, st%d%dim, zpsi)
1735 call states_elec_set_state(st, mesh, ist, ik, zpsi)
1736 end do
1737 end do
1738 end if
1739
1740 end select
1741
1742 safe_deallocate_a(dpsi)
1743 safe_deallocate_a(zpsi)
1744
1746 end subroutine states_elec_generate_random
1747
1748 ! ---------------------------------------------------------
1750 !
1751 subroutine states_elec_fermi(st, namespace, mesh, compute_spin)
1752 type(states_elec_t), intent(inout) :: st
1753 type(namespace_t), intent(in) :: namespace
1754 class(mesh_t), intent(in) :: mesh
1755 logical, optional, intent(in) :: compute_spin
1756
1758 integer :: ist, ik
1759 real(real64) :: charge
1760 complex(real64), allocatable :: zpsi(:, :)
1761
1762 push_sub(states_elec_fermi)
1763
1764 call smear_find_fermi_energy(st%smear, namespace, st%eigenval, st%occ, st%qtot, &
1765 st%nik, st%nst, st%kweights)
1766
1767 call smear_fill_occupations(st%smear, st%eigenval, st%occ, st%nik, st%nst)
1768
1769 ! check if everything is OK
1770 charge = m_zero
1771 do ist = 1, st%nst
1772 charge = charge + sum(st%occ(ist, 1:st%nik) * st%kweights(1:st%nik))
1773 end do
1774 if (abs(charge-st%qtot) > 1e-6_real64) then
1775 message(1) = 'Occupations do not integrate to total charge.'
1776 write(message(2), '(6x,f12.8,a,f12.8)') charge, ' != ', st%qtot
1777 call messages_warning(2, namespace=namespace)
1778 if (charge < m_epsilon) then
1779 message(1) = "There don't seem to be any electrons at all!"
1780 call messages_fatal(1, namespace=namespace)
1781 end if
1782 end if
1783
1784 if (st%d%ispin == spinors .and. optional_default(compute_spin,.true.)) then
1785 assert(states_are_complex(st))
1786
1787 st%spin(:,:,:) = m_zero
1788
1789 safe_allocate(zpsi(1:mesh%np, st%d%dim))
1790 do ik = st%d%kpt%start, st%d%kpt%end
1791 do ist = st%st_start, st%st_end
1792 call states_elec_get_state(st, mesh, ist, ik, zpsi)
1793 st%spin(1:3, ist, ik) = state_spin(mesh, zpsi)
1794 end do
1795 end do
1796 safe_deallocate_a(zpsi)
1797
1798 if (st%parallel_in_states .or. st%d%kpt%parallel) then
1799 call comm_allreduce(st%st_kpt_mpi_grp, st%spin)
1800 end if
1801
1802 end if
1803
1804 pop_sub(states_elec_fermi)
1805 end subroutine states_elec_fermi
1806
1807
1808 ! ---------------------------------------------------------
1810 !
1811 function states_elec_eigenvalues_sum(st, alt_eig) result(tot)
1812 type(states_elec_t), intent(in) :: st
1813 real(real64), optional, intent(in) :: alt_eig(st%st_start:, st%d%kpt%start:)
1815 real(real64) :: tot
1816
1817 integer :: ik
1818
1820
1821 tot = m_zero
1822 do ik = st%d%kpt%start, st%d%kpt%end
1823 if (present(alt_eig)) then
1824 tot = tot + st%kweights(ik) * sum(st%occ(st%st_start:st%st_end, ik) * &
1825 alt_eig(st%st_start:st%st_end, ik))
1826 else
1827 tot = tot + st%kweights(ik) * sum(st%occ(st%st_start:st%st_end, ik) * &
1828 st%eigenval(st%st_start:st%st_end, ik))
1829 end if
1830 end do
1831
1832 if (st%parallel_in_states .or. st%d%kpt%parallel) call comm_allreduce(st%st_kpt_mpi_grp, tot)
1833
1835 end function states_elec_eigenvalues_sum
1836
1837
1839 subroutine states_elec_distribute_nodes(st, namespace, mc)
1840 type(states_elec_t), intent(inout) :: st
1841 type(namespace_t), intent(in) :: namespace
1842 type(multicomm_t), intent(in) :: mc
1843
1844 logical :: default_scalapack_compatible
1845
1847
1848 ! TODO(Alex) Issue #820. This is superflous. These defaults are set in initialisation of
1849 ! states, and in the state distribution instance
1850 ! Defaults.
1851 st%node(:) = 0
1852 st%st_start = 1
1853 st%st_end = st%nst
1854 st%lnst = st%nst
1855 st%parallel_in_states = .false.
1856
1857 call mpi_grp_init(st%mpi_grp, mc%group_comm(p_strategy_states))
1858 call mpi_grp_init(st%dom_st_kpt_mpi_grp, mc%dom_st_kpt_comm)
1859 call mpi_grp_init(st%dom_st_mpi_grp, mc%dom_st_comm)
1860 call mpi_grp_init(st%st_kpt_mpi_grp, mc%st_kpt_comm)
1861
1862 default_scalapack_compatible = calc_mode_par%scalapack_compat() .and. .not. st%d%kpt%parallel
1863
1864 !%Variable ScaLAPACKCompatible
1865 !%Type logical
1866 !%Section Execution::Parallelization
1867 !%Description
1868 !% Whether to use a layout for states parallelization which is compatible with ScaLAPACK.
1869 !% The default is yes for <tt>CalculationMode = gs, unocc, go</tt> without k-point parallelization,
1870 !% and no otherwise. (Setting to other than default is experimental.)
1871 !% The value must be yes if any ScaLAPACK routines are called in the course of the run;
1872 !% it must be set by hand for <tt>td</tt> with <tt>TDDynamics = bo</tt>.
1873 !% This variable has no effect unless you are using states parallelization and have linked ScaLAPACK.
1874 !% Note: currently, use of ScaLAPACK is not compatible with task parallelization (<i>i.e.</i> slaves).
1875 !%End
1876 call parse_variable(namespace, 'ScaLAPACKCompatible', default_scalapack_compatible, st%scalapack_compatible)
1877
1878#ifdef HAVE_SCALAPACK
1879 if (default_scalapack_compatible .neqv. st%scalapack_compatible) then
1880 call messages_experimental('Setting ScaLAPACKCompatible to other than default', namespace=namespace)
1881 end if
1882
1883 if (st%scalapack_compatible) then
1884 if (multicomm_have_slaves(mc)) then
1885 call messages_not_implemented("ScaLAPACK usage with task parallelization (slaves)", namespace=namespace)
1886 end if
1887 call blacs_proc_grid_init(st%dom_st_proc_grid, st%dom_st_mpi_grp)
1888 end if
1889#else
1890 st%scalapack_compatible = .false.
1891#endif
1892
1894
1895#ifdef HAVE_MPI
1896 call multicomm_create_all_pairs(st%mpi_grp, st%ap)
1897#endif
1898
1899 if (st%nst < st%mpi_grp%size) then
1900 message(1) = "Have more processors than necessary"
1901 write(message(2),'(i4,a,i4,a)') st%mpi_grp%size, " processors and ", st%nst, " states."
1902 call messages_fatal(2, namespace=namespace)
1903 end if
1905 call distributed_init(st%dist, st%nst, st%mpi_grp%comm, "states", scalapack_compat = st%scalapack_compatible)
1906
1907 st%parallel_in_states = st%dist%parallel
1908
1909 ! TODO(Alex) Issue #820. Remove lnst, st_start, st_end and node, as they are all contained within dist
1910 st%st_start = st%dist%start
1911 st%st_end = st%dist%end
1912 st%lnst = st%dist%nlocal
1913 st%node(1:st%nst) = st%dist%node(1:st%nst)
1914
1915 end if
1916
1918
1920 end subroutine states_elec_distribute_nodes
1921
1922
1928 !
1929 subroutine states_elec_calc_quantities(gr, st, kpoints, nlcc, &
1930 kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, &
1931 gi_kinetic_energy_density, st_end)
1932 type(grid_t), intent(in) :: gr
1933 type(states_elec_t), intent(in) :: st
1934 type(kpoints_t), intent(in) :: kpoints
1935 logical, intent(in) :: nlcc
1936 real(real64), contiguous, optional, target, intent(out) :: kinetic_energy_density(:,:)
1937 real(real64), contiguous, optional, target, intent(out) :: paramagnetic_current(:,:,:)
1938 real(real64), contiguous, optional, intent(out) :: density_gradient(:,:,:)
1939 real(real64), contiguous, optional, intent(out) :: density_laplacian(:,:)
1940 real(real64), contiguous, optional, intent(out) :: gi_kinetic_energy_density(:,:)
1941 integer, optional, intent(in) :: st_end
1942
1943 real(real64), pointer, contiguous :: jp(:, :, :)
1944 real(real64), pointer, contiguous :: tau(:, :)
1945 complex(real64), allocatable :: wf_psi(:,:), gwf_psi(:,:,:), wf_psi_conj(:,:), lwf_psi(:,:)
1946 real(real64), allocatable :: abs_wf_psi(:,:), abs_gwf_psi(:,:)
1947 complex(real64), allocatable :: psi_gpsi(:,:)
1948 complex(real64) :: c_tmp
1949 integer :: is, ik, ist, i_dim, st_dim, ii, st_end_
1950 real(real64) :: ww, kpoint(gr%der%dim)
1951 logical :: something_to_do
1952
1953 call profiling_in("STATES_CALC_QUANTITIES")
1954
1956
1957 st_end_ = min(st%st_end, optional_default(st_end, st%st_end))
1958
1959 something_to_do = present(kinetic_energy_density) .or. present(gi_kinetic_energy_density) .or. &
1960 present(paramagnetic_current) .or. present(density_gradient) .or. present(density_laplacian)
1961 assert(something_to_do)
1962
1963 safe_allocate( wf_psi(1:gr%np_part, 1:st%d%dim))
1964 safe_allocate( wf_psi_conj(1:gr%np_part, 1:st%d%dim))
1965 safe_allocate(gwf_psi(1:gr%np, 1:gr%der%dim, 1:st%d%dim))
1966 safe_allocate(abs_wf_psi(1:gr%np, 1:st%d%dim))
1967 safe_allocate(abs_gwf_psi(1:gr%np, 1:st%d%dim))
1968 safe_allocate(psi_gpsi(1:gr%np, 1:st%d%dim))
1969 if (present(density_laplacian)) then
1970 safe_allocate(lwf_psi(1:gr%np, 1:st%d%dim))
1971 end if
1972
1973 nullify(tau)
1974 if (present(kinetic_energy_density)) tau => kinetic_energy_density
1975
1976 nullify(jp)
1977 if (present(paramagnetic_current)) jp => paramagnetic_current
1978
1979 ! for the gauge-invariant kinetic energy density we need the
1980 ! current and the kinetic energy density
1981 if (present(gi_kinetic_energy_density)) then
1982 if (.not. present(paramagnetic_current) .and. states_are_complex(st)) then
1983 safe_allocate(jp(1:gr%np, 1:gr%der%dim, 1:st%d%nspin))
1984 end if
1985 if (.not. present(kinetic_energy_density)) then
1986 safe_allocate(tau(1:gr%np, 1:st%d%nspin))
1987 end if
1988 end if
1989
1990 if (associated(tau)) tau = m_zero
1991 if (associated(jp)) jp = m_zero
1992 if (present(density_gradient)) density_gradient(:,:,:) = m_zero
1993 if (present(density_laplacian)) density_laplacian(:,:) = m_zero
1994 if (present(gi_kinetic_energy_density)) gi_kinetic_energy_density = m_zero
1995
1996 do ik = st%d%kpt%start, st%d%kpt%end
1997
1998 kpoint(1:gr%der%dim) = kpoints%get_point(st%d%get_kpoint_index(ik))
1999 is = st%d%get_spin_index(ik)
2000
2001 do ist = st%st_start, st_end_
2002 ww = st%kweights(ik)*st%occ(ist, ik)
2003 if (abs(ww) <= m_epsilon) cycle
2004
2005 ! all calculations will be done with complex wavefunctions
2006 call states_elec_get_state(st, gr, ist, ik, wf_psi)
2007
2008 do st_dim = 1, st%d%dim
2009 call boundaries_set(gr%der%boundaries, gr, wf_psi(:, st_dim))
2010 end do
2011
2012 ! calculate gradient of the wavefunction
2013 do st_dim = 1, st%d%dim
2014 call zderivatives_grad(gr%der, wf_psi(:,st_dim), gwf_psi(:,:,st_dim), set_bc = .false.)
2015 end do
2016
2017 ! calculate the Laplacian of the wavefunction
2018 if (present(density_laplacian)) then
2019 do st_dim = 1, st%d%dim
2020 call zderivatives_lapl(gr%der, wf_psi(:,st_dim), lwf_psi(:,st_dim), ghost_update = .false., set_bc = .false.)
2021 end do
2022 end if
2023
2024 ! We precompute some quantites, to avoid to compute it many times
2025 wf_psi_conj(1:gr%np, 1:st%d%dim) = conjg(wf_psi(1:gr%np,1:st%d%dim))
2026 abs_wf_psi(1:gr%np, 1:st%d%dim) = real(wf_psi_conj(1:gr%np, 1:st%d%dim) * wf_psi(1:gr%np, 1:st%d%dim), real64)
2027
2028 if (present(density_laplacian)) then
2029 density_laplacian(1:gr%np, is) = density_laplacian(1:gr%np, is) + &
2030 ww * m_two*real(wf_psi_conj(1:gr%np, 1) * lwf_psi(1:gr%np, 1), real64)
2031 if (st%d%ispin == spinors) then
2032 density_laplacian(1:gr%np, 2) = density_laplacian(1:gr%np, 2) + &
2033 ww * m_two*real(wf_psi_conj(1:gr%np, 2) * lwf_psi(1:gr%np, 2), real64)
2034 !$omp parallel do private(c_tmp)
2035 do ii = 1, gr%np
2036 c_tmp = ww*(lwf_psi(ii, 1) * wf_psi_conj(ii, 2) + wf_psi(ii, 1) * conjg(lwf_psi(ii, 2)))
2037 density_laplacian(ii, 3) = density_laplacian(ii, 3) + real(c_tmp, real64)
2038 density_laplacian(ii, 4) = density_laplacian(ii, 4) + aimag(c_tmp)
2039 end do
2040 end if
2041 end if
2042
2043 if (associated(tau)) then
2044 tau(1:gr%np, is) = tau(1:gr%np, is) &
2045 + ww * sum(kpoint(1:gr%der%dim)**2) * abs_wf_psi(1:gr%np, 1)
2046 if (st%d%ispin == spinors) then
2047 tau(1:gr%np, 2) = tau(1:gr%np, 2) &
2048 + ww * sum(kpoint(1:gr%der%dim)**2) * abs_wf_psi(1:gr%np, 2)
2049
2050 !$omp parallel do private(c_tmp)
2051 do ii = 1, gr%np
2052 c_tmp = ww * sum(kpoint(1:gr%der%dim)**2) * wf_psi(ii, 1) * wf_psi_conj(ii, 2)
2053 tau(ii, 3) = tau(ii, 3) + real(c_tmp, real64)
2054 tau(ii, 4) = tau(ii, 4) + aimag(c_tmp)
2055 end do
2056 end if
2057 end if
2058
2059 do i_dim = 1, gr%der%dim
2060
2061 ! We precompute some quantites, to avoid to compute them many times
2062 psi_gpsi(1:gr%np, 1:st%d%dim) = wf_psi_conj(1:gr%np, 1:st%d%dim) &
2063 * gwf_psi(1:gr%np,i_dim,1:st%d%dim)
2064 abs_gwf_psi(1:gr%np, 1:st%d%dim) = real(conjg(gwf_psi(1:gr%np, i_dim, 1:st%d%dim)) &
2065 * gwf_psi(1:gr%np, i_dim, 1:st%d%dim), real64)
2066
2067 if (present(density_gradient)) then
2068 density_gradient(1:gr%np, i_dim, is) = density_gradient(1:gr%np, i_dim, is) &
2069 + ww * m_two * real(psi_gpsi(1:gr%np, 1), real64)
2070 if (st%d%ispin == spinors) then
2071 density_gradient(1:gr%np, i_dim, 2) = density_gradient(1:gr%np, i_dim, 2) &
2072 + ww * m_two*real(psi_gpsi(1:gr%np, 2), real64)
2073 !$omp parallel do private(c_tmp)
2074 do ii = 1, gr%np
2075 c_tmp = ww * (gwf_psi(ii, i_dim, 1) * wf_psi_conj(ii, 2) + wf_psi(ii, 1) * conjg(gwf_psi(ii, i_dim, 2)))
2076 density_gradient(ii, i_dim, 3) = density_gradient(ii, i_dim, 3) + real(c_tmp, real64)
2077 density_gradient(ii, i_dim, 4) = density_gradient(ii, i_dim, 4) + aimag(c_tmp)
2078 end do
2079 end if
2080 end if
2081
2082 if (present(density_laplacian)) then
2083 call lalg_axpy(gr%np, ww*m_two, abs_gwf_psi(:,1), density_laplacian(:,is))
2084 if (st%d%ispin == spinors) then
2085 call lalg_axpy(gr%np, ww*m_two, abs_gwf_psi(:,2), density_laplacian(:,2))
2086 !$omp parallel do private(c_tmp)
2087 do ii = 1, gr%np
2088 c_tmp = m_two * ww * gwf_psi(ii, i_dim, 1) * conjg(gwf_psi(ii, i_dim, 2))
2089 density_laplacian(ii, 3) = density_laplacian(ii, 3) + real(c_tmp, real64)
2090 density_laplacian(ii, 4) = density_laplacian(ii, 4) + aimag(c_tmp)
2091 end do
2092 end if
2093 end if
2094
2095 ! the expression for the paramagnetic current with spinors is
2096 ! j = ( jp(1) jp(3) + i jp(4) )
2097 ! (-jp(3) + i jp(4) jp(2) )
2098 if (associated(jp)) then
2099 if (.not.(states_are_real(st))) then
2100 jp(1:gr%np, i_dim, is) = jp(1:gr%np, i_dim, is) + &
2101 ww*(aimag(psi_gpsi(1:gr%np, 1)) - abs_wf_psi(1:gr%np, 1)*kpoint(i_dim))
2102 if (st%d%ispin == spinors) then
2103 jp(1:gr%np, i_dim, 2) = jp(1:gr%np, i_dim, 2) + &
2104 ww*( aimag(psi_gpsi(1:gr%np, 2)) - abs_wf_psi(1:gr%np, 2)*kpoint(i_dim))
2105 !$omp parallel do private(c_tmp)
2106 do ii = 1, gr%np
2107 c_tmp = -ww*m_half*m_zi*(gwf_psi(ii, i_dim, 1)*wf_psi_conj(ii, 2) - wf_psi(ii, 1)*conjg(gwf_psi(ii, i_dim, 2)) &
2108 - m_two * m_zi*wf_psi(ii, 1)*wf_psi_conj(ii, 2)*kpoint(i_dim))
2109 jp(ii, i_dim, 3) = jp(ii, i_dim, 3) + real(c_tmp, real64)
2110 jp(ii, i_dim, 4) = jp(ii, i_dim, 4) + aimag(c_tmp)
2111 end do
2112 end if
2113 end if
2114 end if
2115
2116 ! the expression for the paramagnetic current with spinors is
2117 ! t = ( tau(1) tau(3) + i tau(4) )
2118 ! ( tau(3) - i tau(4) tau(2) )
2119 if (associated(tau)) then
2120 tau(1:gr%np, is) = tau(1:gr%np, is) + ww*(abs_gwf_psi(1:gr%np,1) &
2121 - m_two*aimag(psi_gpsi(1:gr%np, 1))*kpoint(i_dim))
2122 if (st%d%ispin == spinors) then
2123 tau(1:gr%np, 2) = tau(1:gr%np, 2) + ww*(abs_gwf_psi(1:gr%np, 2) &
2124 - m_two*aimag(psi_gpsi(1:gr%np, 2))*kpoint(i_dim))
2125 !$omp parallel do private(c_tmp)
2126 do ii = 1, gr%np
2127 c_tmp = ww * ( gwf_psi(ii, i_dim, 1)*conjg(gwf_psi(ii, i_dim, 2)) &
2128 + m_zi * (gwf_psi(ii,i_dim,1)*wf_psi_conj(ii, 2) &
2129 - wf_psi(ii, 1)*conjg(gwf_psi(ii,i_dim,2)))*kpoint(i_dim))
2130 tau(ii, 3) = tau(ii, 3) + real(c_tmp, real64)
2131 tau(ii, 4) = tau(ii, 4) + aimag(c_tmp)
2132 end do
2133 end if
2134 end if
2135
2136 end do
2137
2138 end do
2139 end do
2140
2141 safe_deallocate_a(wf_psi_conj)
2142 safe_deallocate_a(abs_wf_psi)
2143 safe_deallocate_a(abs_gwf_psi)
2144 safe_deallocate_a(psi_gpsi)
2145
2146 if (.not. present(gi_kinetic_energy_density)) then
2147 if (.not. present(paramagnetic_current)) then
2148 safe_deallocate_p(jp)
2149 end if
2150 if (.not. present(kinetic_energy_density)) then
2151 safe_deallocate_p(tau)
2152 end if
2153 end if
2154
2155 if (st%parallel_in_states .or. st%d%kpt%parallel) call reduce_all(st%st_kpt_mpi_grp)
2156
2157 ! We have to symmetrize everything as they are calculated from the
2158 ! wavefunctions.
2159 ! This must be done before compute the gauge-invariant kinetic energy density
2160 if (st%symmetrize_density) then
2161 do is = 1, st%d%nspin
2162 if (associated(tau)) then
2163 call dgrid_symmetrize_scalar_field(gr, tau(:, is), suppress_warning = .true.)
2164 end if
2165
2166 if (present(density_laplacian)) then
2167 call dgrid_symmetrize_scalar_field(gr, density_laplacian(:, is), suppress_warning = .true.)
2168 end if
2169
2170 if (associated(jp)) then
2171 call dgrid_symmetrize_vector_field(gr, jp(:, :, is), suppress_warning = .true.)
2172 end if
2173
2174 if (present(density_gradient)) then
2175 call dgrid_symmetrize_vector_field(gr, density_gradient(:, :, is), suppress_warning = .true.)
2176 end if
2177 end do
2178 end if
2179
2180
2181 if (allocated(st%rho_core) .and. nlcc .and. (present(density_laplacian) .or. present(density_gradient))) then
2182 do ii = 1, gr%np
2183 wf_psi(ii, 1) = st%rho_core(ii)/st%d%spin_channels
2184 end do
2185
2186 call boundaries_set(gr%der%boundaries, gr, wf_psi(:, 1))
2187
2188 if (present(density_gradient)) then
2189 ! calculate gradient of the NLCC
2190 call zderivatives_grad(gr%der, wf_psi(:,1), gwf_psi(:,:,1), set_bc = .false.)
2191 do is = 1, st%d%spin_channels
2192 density_gradient(1:gr%np, 1:gr%der%dim, is) = density_gradient(1:gr%np, 1:gr%der%dim, is) + &
2193 real(gwf_psi(1:gr%np, 1:gr%der%dim,1))
2194 end do
2195 end if
2196
2197 ! calculate the Laplacian of the wavefunction
2198 if (present(density_laplacian)) then
2199 call zderivatives_lapl(gr%der, wf_psi(:,1), lwf_psi(:,1), set_bc = .false.)
2200
2201 do is = 1, st%d%spin_channels
2202 density_laplacian(1:gr%np, is) = density_laplacian(1:gr%np, is) + real(lwf_psi(1:gr%np, 1))
2203 end do
2204 end if
2205 end if
2206
2207 !If we freeze some of the orbitals, we need to had the contributions here
2208 !Only in the case we are not computing it
2209 if (allocated(st%frozen_tau) .and. .not. present(st_end)) then
2210 call lalg_axpy(gr%np, st%d%nspin, m_one, st%frozen_tau, tau)
2211 end if
2212 if (allocated(st%frozen_gdens) .and. .not. present(st_end)) then
2213 do is = 1, st%d%nspin
2214 call lalg_axpy(gr%np, gr%der%dim, m_one, st%frozen_gdens(:,:,is), density_gradient(:,:,is))
2215 end do
2216 end if
2217 if (allocated(st%frozen_tau) .and. .not. present(st_end)) then
2218 call lalg_axpy(gr%np, st%d%nspin, m_one, st%frozen_ldens, density_laplacian)
2219 end if
2220
2221 safe_deallocate_a(wf_psi)
2222 safe_deallocate_a(gwf_psi)
2223 safe_deallocate_a(lwf_psi)
2224
2225
2226 !We compute the gauge-invariant kinetic energy density
2227 if (present(gi_kinetic_energy_density)) then
2228 do is = 1, st%d%nspin
2229 assert(associated(tau))
2230 gi_kinetic_energy_density(1:gr%np, is) = tau(1:gr%np, is)
2231 end do
2232 if (states_are_complex(st)) then
2233 assert(associated(jp))
2234 if (st%d%ispin /= spinors) then
2235 do is = 1, st%d%nspin
2236 !$omp parallel do
2237 do ii = 1, gr%np
2238 if (st%rho(ii, is) < 1.0e-7_real64) cycle
2239 gi_kinetic_energy_density(ii, is) = &
2240 gi_kinetic_energy_density(ii, is) - sum(jp(ii,1:gr%der%dim, is)**2)/st%rho(ii, is)
2241 end do
2242 end do
2243 else ! Note that this is only the U(1) part of the gauge-invariant KED
2244 !$omp parallel do
2245 do ii = 1, gr%np
2246 gi_kinetic_energy_density(ii, 1) = gi_kinetic_energy_density(ii, 1) &
2247 - sum(jp(ii,1:gr%der%dim, 1)**2)/(safe_tol(st%rho(ii, 1),m_epsilon))
2248 gi_kinetic_energy_density(ii, 2) = gi_kinetic_energy_density(ii, 2) &
2249 - sum(jp(ii,1:gr%der%dim, 2)**2)/(safe_tol(st%rho(ii, 2),m_epsilon))
2250 gi_kinetic_energy_density(ii, 3) = &
2251 gi_kinetic_energy_density(ii, 3) - sum(jp(ii,1:gr%der%dim, 3)**2 + jp(ii,1:gr%der%dim, 4)**2) &
2252 /(safe_tol((st%rho(ii, 3)**2 + st%rho(ii, 4)**2), m_epsilon))*st%rho(ii, 3)
2253 gi_kinetic_energy_density(ii, 4) = &
2254 gi_kinetic_energy_density(ii, 4) + sum(jp(ii,1:gr%der%dim, 3)**2 + jp(ii,1:gr%der%dim, 4)**2) &
2255 /(safe_tol((st%rho(ii, 3)**2 + st%rho(ii, 4)**2), m_epsilon))*st%rho(ii, 4)
2256 end do
2257 end if
2258 end if
2259 end if
2260
2261 if (.not. present(kinetic_energy_density)) then
2262 safe_deallocate_p(tau)
2263 end if
2264 if (.not. present(paramagnetic_current)) then
2265 safe_deallocate_p(jp)
2266 end if
2267
2268
2270
2271 call profiling_out("STATES_CALC_QUANTITIES")
2272
2273 contains
2274
2275 subroutine reduce_all(grp)
2276 type(mpi_grp_t), intent(in) :: grp
2277
2279
2280 if (associated(tau)) call comm_allreduce(grp, tau, dim = (/gr%np, st%d%nspin/))
2281
2282 if (present(density_laplacian)) call comm_allreduce(grp, density_laplacian, dim = (/gr%np, st%d%nspin/))
2283
2284 do is = 1, st%d%nspin
2285 if (associated(jp)) call comm_allreduce(grp, jp(:, :, is), dim = (/gr%np, gr%der%dim/))
2286
2287 if (present(density_gradient)) then
2288 call comm_allreduce(grp, density_gradient(:, :, is), dim = (/gr%np, gr%der%dim/))
2289 end if
2290 end do
2291
2293 end subroutine reduce_all
2294
2295 end subroutine states_elec_calc_quantities
2296
2297
2298 ! ---------------------------------------------------------
2300 !
2301 function state_spin(mesh, f1) result(spin)
2302 type(mesh_t), intent(in) :: mesh
2303 complex(real64), intent(in) :: f1(:, :)
2304 real(real64) :: spin(1:3)
2305
2306 complex(real64) :: z
2307
2308 push_sub(state_spin)
2309
2310 z = zmf_dotp(mesh, f1(:, 1) , f1(:, 2))
2311
2312 spin(1) = m_two*real(z, real64)
2313 spin(2) = m_two*aimag(z)
2314 spin(3) = zmf_nrm2(mesh, f1(:, 1))**2 - zmf_nrm2(mesh, f1(:, 2))**2
2315 spin = m_half*spin ! spin is half the sigma matrix.
2316
2317 pop_sub(state_spin)
2318 end function state_spin
2319
2320 ! ---------------------------------------------------------
2322 !
2323 logical function state_is_local(st, ist)
2324 type(states_elec_t), intent(in) :: st
2325 integer, intent(in) :: ist
2326
2327 push_sub(state_is_local)
2328
2329 state_is_local = ist >= st%st_start.and.ist <= st%st_end
2330
2331 pop_sub(state_is_local)
2332 end function state_is_local
2333
2334 ! ---------------------------------------------------------
2336 !
2337 logical function state_kpt_is_local(st, ist, ik)
2338 type(states_elec_t), intent(in) :: st
2339 integer, intent(in) :: ist
2340 integer, intent(in) :: ik
2341
2342 push_sub(state_kpt_is_local)
2343
2344 state_kpt_is_local = ist >= st%st_start .and. ist <= st%st_end .and. &
2345 ik >= st%d%kpt%start .and. ik <= st%d%kpt%end
2346
2347 pop_sub(state_kpt_is_local)
2348 end function state_kpt_is_local
2349
2350
2351 ! ---------------------------------------------------------
2353 real(real64) function states_elec_wfns_memory(st, mesh) result(memory)
2354 type(states_elec_t), intent(in) :: st
2355 class(mesh_t), intent(in) :: mesh
2356
2357 push_sub(states_elec_wfns_memory)
2358 memory = m_zero
2359
2360 ! orbitals
2361 memory = memory + real64*real(mesh%np_part_global, real64) *st%d%dim*real(st%nst, real64) *st%d%kpt%nglobal
2362
2364 end function states_elec_wfns_memory
2365
2366 ! ---------------------------------------------------------
2369 subroutine states_elec_pack(st, copy)
2370 class(states_elec_t), intent(inout) :: st
2371 logical, optional, intent(in) :: copy
2372
2373 integer :: iqn, ib
2374 integer(int64) :: max_mem, mem
2375
2376 push_sub(states_elec_pack)
2377
2378 ! nothing to do, already packed
2379 if (st%packed) then
2380 pop_sub(states_elec_pack)
2381 return
2382 end if
2383
2384 st%packed = .true.
2385
2386 if (accel_is_enabled()) then
2387 max_mem = accel_global_memory_size()
2388
2389 if (st%cl_states_mem > m_one) then
2390 max_mem = int(st%cl_states_mem, int64)*(1024_8)**2
2391 else if (st%cl_states_mem < 0.0_real64) then
2392 max_mem = max_mem + int(st%cl_states_mem, int64)*(1024_8)**2
2393 else
2394 max_mem = int(st%cl_states_mem*real(max_mem, real64) , int64)
2395 end if
2396 else
2397 max_mem = huge(max_mem)
2398 end if
2399
2400 mem = 0
2401 qnloop: do iqn = st%d%kpt%start, st%d%kpt%end
2402 do ib = st%group%block_start, st%group%block_end
2403
2404 mem = mem + st%group%psib(ib, iqn)%pack_total_size()
2405
2406 if (mem > max_mem) then
2407 call messages_write('Not enough CL device memory to store all states simultaneously.', new_line = .true.)
2408 call messages_write('Only ')
2409 call messages_write(ib - st%group%block_start)
2410 call messages_write(' of ')
2411 call messages_write(st%group%block_end - st%group%block_start + 1)
2412 call messages_write(' blocks will be stored in device memory.', new_line = .true.)
2413 call messages_warning()
2414 exit qnloop
2415 end if
2417 call st%group%psib(ib, iqn)%do_pack(copy)
2418 end do
2419 end do qnloop
2420
2421 pop_sub(states_elec_pack)
2422 end subroutine states_elec_pack
2423
2424 ! ------------------------------------------------------------
2426 !
2427 subroutine states_elec_unpack(st, copy)
2428 class(states_elec_t), intent(inout) :: st
2429 logical, optional, intent(in) :: copy
2431 integer :: iqn, ib
2432
2433 push_sub(states_elec_unpack)
2434
2435 if (st%packed) then
2436 st%packed = .false.
2437
2438 do iqn = st%d%kpt%start, st%d%kpt%end
2439 do ib = st%group%block_start, st%group%block_end
2440 if (st%group%psib(ib, iqn)%is_packed()) call st%group%psib(ib, iqn)%do_unpack(copy)
2441 end do
2442 end do
2443 end if
2444
2445 pop_sub(states_elec_unpack)
2446 end subroutine states_elec_unpack
2447
2448 ! -----------------------------------------------------------
2450 !
2451 subroutine states_elec_write_info(st, namespace)
2452 class(states_elec_t), intent(in) :: st
2453 type(namespace_t), intent(in) :: namespace
2454
2455 push_sub(states_elec_write_info)
2456
2457 call messages_print_with_emphasis(msg="States", namespace=namespace)
2458
2459 write(message(1), '(a,f12.3)') 'Total electronic charge = ', st%qtot
2460 write(message(2), '(a,i8)') 'Number of states = ', st%nst
2461 write(message(3), '(a,i8)') 'States block-size = ', st%block_size
2462 call messages_info(3, namespace=namespace)
2463
2464 call messages_print_with_emphasis(namespace=namespace)
2465
2466 pop_sub(states_elec_write_info)
2467 end subroutine states_elec_write_info
2468
2469 ! ------------------------------------------------------------
2471 !
2472 subroutine states_elec_set_zero(st)
2473 class(states_elec_t), intent(inout) :: st
2474
2475 integer :: iqn, ib
2476
2477 push_sub(states_elec_set_zero)
2478
2479 do iqn = st%d%kpt%start, st%d%kpt%end
2480 do ib = st%group%block_start, st%group%block_end
2481 call batch_set_zero(st%group%psib(ib, iqn))
2482 end do
2483 end do
2484
2485 pop_sub(states_elec_set_zero)
2486 end subroutine states_elec_set_zero
2487
2488 ! ------------------------------------------------------------
2490 !
2491 integer pure function states_elec_block_min(st, ib) result(range)
2492 type(states_elec_t), intent(in) :: st
2493 integer, intent(in) :: ib
2494
2495 range = st%group%block_range(ib, 1)
2496 end function states_elec_block_min
2497
2498 ! ------------------------------------------------------------
2500 !
2501 integer pure function states_elec_block_max(st, ib) result(range)
2502 type(states_elec_t), intent(in) :: st
2503 integer, intent(in) :: ib
2504
2505 range = st%group%block_range(ib, 2)
2506 end function states_elec_block_max
2507
2508 ! ------------------------------------------------------------
2510 !
2511 integer pure function states_elec_block_size(st, ib) result(size)
2512 type(states_elec_t), intent(in) :: st
2513 integer, intent(in) :: ib
2514
2515 size = st%group%block_size(ib)
2516 end function states_elec_block_size
2517
2518 ! ---------------------------------------------------------
2521 subroutine states_elec_count_pairs(st, namespace, n_pairs, n_occ, n_unocc, is_included, is_frac_occ)
2522 type(states_elec_t), intent(in) :: st
2523 type(namespace_t), intent(in) :: namespace
2524 integer, intent(out) :: n_pairs
2525 integer, intent(out) :: n_occ(:)
2526 integer, intent(out) :: n_unocc(:)
2527 logical, allocatable, intent(out) :: is_included(:,:,:)
2529 logical, intent(out) :: is_frac_occ
2530
2531 integer :: ik, ist, ast, n_filled, n_partially_filled, n_half_filled
2532 character(len=80) :: nst_string, default, wfn_list
2533 real(real64) :: energy_window
2534
2535 push_sub(states_elec_count_pairs)
2536
2537 is_frac_occ = .false.
2538 do ik = 1, st%nik
2539 call occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled)
2540 if (n_partially_filled > 0 .or. n_half_filled > 0) is_frac_occ = .true.
2541 n_occ(ik) = n_filled + n_partially_filled + n_half_filled
2542 n_unocc(ik) = st%nst - n_filled
2543 ! when we implement occupations, partially occupied levels need to be counted as both occ and unocc.
2544 end do
2545
2546 !%Variable CasidaKSEnergyWindow
2547 !%Type float
2548 !%Section Linear Response::Casida
2549 !%Description
2550 !% An alternative to <tt>CasidaKohnShamStates</tt> for specifying which occupied-unoccupied
2551 !% transitions will be used: all those whose eigenvalue differences are less than this
2552 !% number will be included. If a value less than 0 is supplied, this criterion will not be used.
2553 !%End
2554
2555 call parse_variable(namespace, 'CasidaKSEnergyWindow', -m_one, energy_window, units_inp%energy)
2556
2557 !%Variable CasidaKohnShamStates
2558 !%Type string
2559 !%Section Linear Response::Casida
2560 !%Default all states
2561 !%Description
2562 !% The calculation of the excitation spectrum of a system in the Casida frequency-domain
2563 !% formulation of linear-response time-dependent density functional theory (TDDFT)
2564 !% implies the use of a basis set of occupied/unoccupied Kohn-Sham orbitals. This
2565 !% basis set should, in principle, include all pairs formed by all occupied states,
2566 !% and an infinite number of unoccupied states. In practice, one has to truncate this
2567 !% basis set, selecting a number of occupied and unoccupied states that will form the
2568 !% pairs. These states are specified with this variable. If there are, say, 15 occupied
2569 !% states, and one sets this variable to the value "10-18", this means that occupied
2570 !% states from 10 to 15, and unoccupied states from 16 to 18 will be considered.
2571 !%
2572 !% This variable is a string in list form, <i>i.e.</i> expressions such as "1,2-5,8-15" are
2573 !% valid. You should include a non-zero number of unoccupied states and a non-zero number
2574 !% of occupied states.
2575 !%End
2576
2577 n_pairs = 0
2578 safe_allocate(is_included(maxval(n_occ), minval(n_occ) + 1:st%nst , st%nik))
2579 is_included(:,:,:) = .false.
2580
2581 if (energy_window < m_zero) then
2582 write(nst_string,'(i6)') st%nst
2583 write(default,'(a,a)') "1-", trim(adjustl(nst_string))
2584 call parse_variable(namespace, 'CasidaKohnShamStates', default, wfn_list)
2585
2586 write(message(1),'(a,a)') "Info: States that form the basis: ", trim(wfn_list)
2587 call messages_info(1, namespace=namespace)
2588
2589 ! count pairs
2590 n_pairs = 0
2591 do ik = 1, st%nik
2592 do ast = n_occ(ik) + 1, st%nst
2593 if (loct_isinstringlist(ast, wfn_list)) then
2594 do ist = 1, n_occ(ik)
2595 if (loct_isinstringlist(ist, wfn_list)) then
2596 n_pairs = n_pairs + 1
2597 is_included(ist, ast, ik) = .true.
2598 end if
2599 end do
2600 end if
2601 end do
2602 end do
2603
2604 else ! using CasidaKSEnergyWindow
2605
2606 write(message(1),'(a,f12.6,a)') "Info: including transitions with energy < ", &
2607 units_from_atomic(units_out%energy, energy_window), trim(units_abbrev(units_out%energy))
2608 call messages_info(1, namespace=namespace)
2609
2610 ! count pairs
2611 n_pairs = 0
2612 do ik = 1, st%nik
2613 do ast = n_occ(ik) + 1, st%nst
2614 do ist = 1, n_occ(ik)
2615 if (st%eigenval(ast, ik) - st%eigenval(ist, ik) < energy_window) then
2616 n_pairs = n_pairs + 1
2617 is_included(ist, ast, ik) = .true.
2618 end if
2619 end do
2620 end do
2621 end do
2622
2623 end if
2624
2626 end subroutine states_elec_count_pairs
2627
2628
2644 !
2645 subroutine occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled, &
2646 filled, partially_filled, half_filled)
2647 type(states_elec_t), intent(in) :: st
2648 type(namespace_t), intent(in) :: namespace
2649 integer, intent(in) :: ik
2650 integer, intent(out) :: n_filled, n_partially_filled, n_half_filled
2651 integer, optional, intent(out) :: filled(:), partially_filled(:), half_filled(:)
2652
2653 integer :: ist
2654
2655 push_sub(occupied_states)
2656
2657 if (present(filled)) filled(:) = 0
2658 if (present(partially_filled)) partially_filled(:) = 0
2659 if (present(half_filled)) half_filled(:) = 0
2660 n_filled = 0
2661 n_partially_filled = 0
2662 n_half_filled = 0
2663
2664 select case (st%d%ispin)
2665 case (unpolarized)
2666 do ist = 1, st%nst
2667 if (abs(st%occ(ist, ik) - m_two) < m_min_occ) then
2668 n_filled = n_filled + 1
2669 if (present(filled)) filled(n_filled) = ist
2670 else if (abs(st%occ(ist, ik) - m_one) < m_min_occ) then
2671 n_half_filled = n_half_filled + 1
2672 if (present(half_filled)) half_filled(n_half_filled) = ist
2673 else if (st%occ(ist, ik) > m_min_occ) then
2674 n_partially_filled = n_partially_filled + 1
2675 if (present(partially_filled)) partially_filled(n_partially_filled) = ist
2676 else if (abs(st%occ(ist, ik)) > m_min_occ) then
2677 write(message(1),*) 'Internal error in occupied_states: Illegal occupation value ', st%occ(ist, ik)
2678 call messages_fatal(1, namespace=namespace)
2679 end if
2680 end do
2681 case (spin_polarized, spinors)
2682 do ist = 1, st%nst
2683 if (abs(st%occ(ist, ik)-m_one) < m_min_occ) then
2684 n_filled = n_filled + 1
2685 if (present(filled)) filled(n_filled) = ist
2686 else if (st%occ(ist, ik) > m_min_occ) then
2687 n_partially_filled = n_partially_filled + 1
2688 if (present(partially_filled)) partially_filled(n_partially_filled) = ist
2689 else if (abs(st%occ(ist, ik)) > m_min_occ) then
2690 write(message(1),*) 'Internal error in occupied_states: Illegal occupation value ', st%occ(ist, ik)
2691 call messages_fatal(1, namespace=namespace)
2692 end if
2693 end do
2694 end select
2695
2696 pop_sub(occupied_states)
2697 end subroutine occupied_states
2698
2699
2701 subroutine kpoints_distribute(this, mc)
2702 type(states_elec_t), intent(inout) :: this
2703 type(multicomm_t), intent(in) :: mc
2704
2705 push_sub(kpoints_distribute)
2706 call distributed_init(this%d%kpt, this%nik, mc%group_comm(p_strategy_kpoints), "k-points")
2707
2708 pop_sub(kpoints_distribute)
2709 end subroutine kpoints_distribute
2710
2711
2712 ! TODO(Alex) Issue #824. Consider converting this to a function to returns `st_kpt_task`
2713 ! as this is only called in a couple of places, or package with the `st_kpt_mpi_grp` when split
2714 ! from st instance
2717 type(states_elec_t), intent(inout) :: st
2718
2720
2721 if (.not. allocated(st%st_kpt_task)) then
2722 safe_allocate(st%st_kpt_task(0:st%st_kpt_mpi_grp%size-1, 1:4))
2723 end if
2724
2725 st%st_kpt_task(0:st%st_kpt_mpi_grp%size-1, :) = 0
2726 st%st_kpt_task(st%st_kpt_mpi_grp%rank, :) = [st%st_start, st%st_end, st%d%kpt%start, st%d%kpt%end]
2727
2728 if (st%parallel_in_states .or. st%d%kpt%parallel) then
2729 call comm_allreduce(st%st_kpt_mpi_grp, st%st_kpt_task)
2730 end if
2731
2734
2735 ! ---------------------------------------------------------
2737 !
2738 subroutine states_elec_choose_kpoints(st, kpoints, namespace)
2739 type(states_elec_t), target, intent(inout) :: st
2740 type(kpoints_t), intent(in) :: kpoints
2741 type(namespace_t), intent(in) :: namespace
2742
2743 integer :: ik, iq
2744 type(states_elec_dim_t), pointer :: dd
2745
2746
2748
2749 dd => st%d
2750
2751 st%nik = kpoints_number(kpoints)
2752
2753 if (dd%ispin == spin_polarized) st%nik = 2*st%nik
2754
2755 safe_allocate(st%kweights(1:st%nik))
2756
2757 do iq = 1, st%nik
2758 ik = dd%get_kpoint_index(iq)
2759 st%kweights(iq) = kpoints%get_weight(ik)
2760 end do
2761
2762 if (debug%info) call print_kpoints_debug
2764
2765 contains
2766 subroutine print_kpoints_debug
2767 integer :: iunit
2768
2770
2771 call io_mkdir('debug/', namespace)
2772 iunit = io_open('debug/kpoints', namespace, action = 'write')
2773 call kpoints%write_info(iunit=iunit, absolute_coordinates = .true.)
2774 call io_close(iunit)
2775
2777 end subroutine print_kpoints_debug
2778
2779 end subroutine states_elec_choose_kpoints
2780
2785 !
2786 function states_elec_calculate_dipole(this, gr) result(dipole)
2787 class(states_elec_t), intent(in) :: this
2788 class(mesh_t), intent(in) :: gr
2789 real(real64) :: dipole(1:gr%box%dim)
2790
2791 integer :: ispin
2792 real(real64) :: e_dip(1:gr%box%dim, this%d%spin_channels)
2793
2795
2796 do ispin = 1, this%d%spin_channels
2797 call dmf_dipole(gr, this%rho(:, ispin), e_dip(:, ispin))
2798 end do
2799
2800 dipole(:) = sum(e_dip(:, 1:this%d%spin_channels), 2) ! dipole moment <mu_el> = \sum_i -e <x_i>
2801
2803 end function states_elec_calculate_dipole
2804
2805
2806#include "undef.F90"
2807#include "real.F90"
2808#include "states_elec_inc.F90"
2810#include "undef.F90"
2811#include "complex.F90"
2812#include "states_elec_inc.F90"
2813#include "undef.F90"
2814
2815end module states_elec_oct_m
2816
2817
2818!! Local Variables:
2819!! mode: f90
2820!! coding: utf-8
2821!! End:
initialize a batch with existing memory
Definition: batch.F90:273
constant times a vector plus a vector
Definition: lalg_basic.F90:171
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:180
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
type(accel_t), public accel
Definition: accel.F90:270
This module implements batches of mesh functions.
Definition: batch.F90:133
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:116
This module provides the BLACS processor grid.
subroutine, public blacs_proc_grid_init(this, mpi_grp, procdim)
Initializes a blacs context from an MPI communicator with topological information.
subroutine, public blacs_proc_grid_end(this)
subroutine, public blacs_proc_grid_copy(cin, cout)
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:122
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public zderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public zderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
subroutine, public distributed_end(this)
subroutine, public distributed_nullify(this, total)
subroutine, public distributed_init(this, total, comm, tag, scalapack_compat)
Distribute N instances across M processes of communicator comm
subroutine, public distributed_copy(in, out)
@Brief Create a copy of a distributed instance
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_fourth
Definition: global.F90:197
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_epsilon
Definition: global.F90:204
type(conf_t), public conf
Global instance of Octopus configuration.
Definition: global.F90:178
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public dgrid_symmetrize_vector_field(gr, field, suppress_warning)
Definition: grid.F90:692
subroutine, public dgrid_symmetrize_scalar_field(gr, field, suppress_warning)
Definition: grid.F90:666
Definition: io.F90:114
logical pure function, public kpoints_point_is_gamma(this, ik)
Definition: kpoints.F90:1549
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
Definition: math.F90:889
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:116
subroutine, public zmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
subroutine, public dmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
This module defines various routines, operating on mesh functions.
subroutine, public zmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
subroutine, public dmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
subroutine, public zmf_normalize(mesh, dim, psi, norm)
Normalize a mesh function psi.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
general module for modelmb particles
subroutine, public modelmb_particles_end(this)
subroutine, public modelmb_particles_init(this, namespace, space)
==============================================================
subroutine, public modelmb_particles_copy(modelmb_out, modelmb_in)
subroutine mpi_grp_copy(mpi_grp_out, mpi_grp_in)
Definition: mpi.F90:407
type(mpi_comm), parameter, public mpi_comm_undefined
used to indicate a communicator has not been initialized
Definition: mpi.F90:136
subroutine mpi_grp_init(grp, comm)
Initialize MPI group instance.
Definition: mpi.F90:350
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
logical pure function, public multicomm_strategy_is_parallel(mc, level)
Definition: multicomm.F90:925
subroutine, public multicomm_all_pairs_copy(apout, apin)
Definition: multicomm.F90:246
logical pure function, public multicomm_have_slaves(this)
Definition: multicomm.F90:1141
subroutine, public multicomm_create_all_pairs(mpi_grp, ap)
This routine uses the one-factorization (or near-one-factorization of a complete graph to construct a...
Definition: multicomm.F90:943
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
subroutine, public restart_read(restart, iunit, lines, nlines, ierr)
Definition: restart.F90:936
subroutine, public restart_close(restart, iunit)
Close a file previously opened with restart_open.
Definition: restart.F90:953
integer function, public restart_open(restart, filename, status, position, silent)
Open file "filename" found inside the current restart directory. Depending on the type of restart,...
Definition: restart.F90:862
subroutine, public smear_fill_occupations(this, eigenvalues, occupations, nik, nst)
Definition: smear.F90:529
subroutine, public smear_find_fermi_energy(this, namespace, eigenvalues, occupations, qtot, nik, nst, kweights)
Definition: smear.F90:344
subroutine, public smear_copy(to, from)
Definition: smear.F90:323
integer, parameter, public smear_fixed_occ
Definition: smear.F90:171
subroutine, public smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
Definition: smear.F90:185
logical pure function, public smear_is_semiconducting(this)
Definition: smear.F90:833
subroutine, public states_set_complex(st)
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_dim_copy(dout, din)
subroutine, public states_elec_dim_end(dim)
This module handles groups of electronic batches and their parallel distribution.
subroutine, public states_elec_group_copy(d, group_in, group_out, copy_data, special)
make a copy of a group
subroutine, public states_elec_group_end(this, d)
finalize the local blocks of wave functions and release local arrays
real(real64) function, public states_elec_wfns_memory(st, mesh)
return the memory usage of a states_elec_t object
integer pure function, public states_elec_block_max(st, ib)
return index of last state in block ib
subroutine zstates_elec_generate_random_vector(mesh, st, vector, normalized, reset_seed)
Generate a random vector.
subroutine, public states_elec_null(st)
subroutine, public states_elec_distribute_nodes(st, namespace, mc)
@Brief. Distribute states over the processes for states parallelization
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
logical function, public state_kpt_is_local(st, ist, ik)
check whether a given state (ist, ik) is on the local node
subroutine, public states_elec_calc_quantities(gr, st, kpoints, nlcc, kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, gi_kinetic_energy_density, st_end)
calculated selected quantities
real(real64) function, dimension(1:gr%box%dim) states_elec_calculate_dipole(this, gr)
calculate the expectation value of the dipole moment of electrons
subroutine, public states_elec_densities_init(st, gr)
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine zstates_elec_get_state1(st, mesh, idim, ist, iqn, psi)
Write one component (dim) of a wave function into a state_elec_t object.
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine zstates_elec_get_points1(st, start_point, end_point, iqn, psi)
Returns the value of all the states for given k-point in the range of points [start_point:end_point].
subroutine states_elec_read_initial_occs(st, namespace, excess_charge, kpoints)
Reads from the input file the initial occupations, if the block "Occupations" is present....
subroutine, public states_elec_init(st, namespace, space, valence_charge, kpoints)
Initialize a new states_elec_t object.
subroutine zstates_elec_set_state1(st, mesh, idim, ist, iqn, psi)
get one dimension of local wave function for given k-point and states index from a states_elec_t obje...
real(real64) function, public states_elec_eigenvalues_sum(st, alt_eig)
function to calculate the eigenvalues sum using occupations as weights
integer, parameter, public par_independent
Method used to generate random states.
subroutine dstates_elec_generate_random_vector(mesh, st, vector, normalized, reset_seed)
Generate a random vector.
subroutine, public occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled, filled, partially_filled, half_filled)
return information about occupied orbitals in many-body state
subroutine dstates_elec_get_state2(st, mesh, ist, iqn, psi)
Write a wave function into a state_elec_t object.
subroutine zstates_elec_set_state2(st, mesh, ist, iqn, psi)
get local wave function for given k-point and states index from a states_elec_t object
subroutine zstates_elec_get_state2(st, mesh, ist, iqn, psi)
Write a wave function into a state_elec_t object.
subroutine dstates_elec_get_points1(st, start_point, end_point, iqn, psi)
Returns the value of all the states for given k-point in the range of points [start_point:end_point].
integer pure function, public states_elec_block_size(st, ib)
return number of states in block ib
subroutine, public kpoints_distribute(this, mc)
distribute k-points over the nodes in the corresponding communicator
integer, parameter, public par_dependent
subroutine states_elec_pack(st, copy)
pack the batches in this states object
subroutine, public states_elec_choose_kpoints(st, kpoints, namespace)
double up k-points for SPIN_POLARIZED calculations
subroutine zstates_elec_get_state3(st, mesh, iqn, psi)
get local wave functions for given k-point from a states_elec_t object
logical function, public state_is_local(st, ist)
check whether a given state (ist) is on the local node
subroutine dstates_elec_set_state1(st, mesh, idim, ist, iqn, psi)
get one dimension of local wave function for given k-point and states index from a states_elec_t obje...
subroutine, public states_elec_exec_init(st, namespace, mc)
Further initializations.
subroutine dstates_elec_set_state4(st, mesh, psi)
set all local wave functions in a states_elec_t object
subroutine dstates_elec_set_state3(st, mesh, iqn, psi)
set local wave functions for given k-point in a states_elec_t object
subroutine states_elec_write_info(st, namespace)
write information about the states object
subroutine states_elec_init_block(st, mesh, verbose, skip, packed)
Initializes the data components in st that describe how the states are distributed in blocks:
real(real64) function, dimension(1:3) state_spin(mesh, f1)
calculate the spin vector for a spinor wave function f1
subroutine states_elec_kpoints_distribution(st)
Assign the start and end indices for states and kpoints, for "st_kpt_mpi_grp" communicator.
subroutine zstates_elec_get_state4(st, mesh, psi)
get all local wave functions from a states_elec_t object
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
randomize states
subroutine dstates_elec_get_state1(st, mesh, idim, ist, iqn, psi)
Write one component (dim) of a wave function into a state_elec_t object.
subroutine dstates_elec_set_state2(st, mesh, ist, iqn, psi)
get local wave function for given k-point and states index from a states_elec_t object
subroutine zstates_elec_set_state4(st, mesh, psi)
set all local wave functions in a states_elec_t object
subroutine states_elec_read_initial_spins(st, namespace)
Reads, if present, the "InitialSpins" block.
subroutine dstates_elec_get_state3(st, mesh, iqn, psi)
get local wave functions for given k-point from a states_elec_t object
subroutine dstates_elec_get_points2(st, start_point, end_point, psi)
Returns the value of all the states in the range of points [start_point:end_point].
integer pure function, public states_elec_block_min(st, ib)
return index of first state in block ib
subroutine states_elec_unpack(st, copy)
unpack the batches in this states object
subroutine, public states_elec_look(restart, nik, dim, nst, ierr)
Reads the 'states' file in the restart directory, and finds out the nik, dim, and nst contained in it...
subroutine dstates_elec_get_state4(st, mesh, psi)
get all local wave functions from a states_elec_t object
subroutine, public states_elec_count_pairs(st, namespace, n_pairs, n_occ, n_unocc, is_included, is_frac_occ)
number of occupied-unoccupied pairs for Casida
subroutine, public states_elec_allocate_current(st, space, mesh)
subroutine, public states_elec_set_zero(st)
expclicitely set all wave functions in the states to zero
subroutine zstates_elec_set_state3(st, mesh, iqn, psi)
set local wave functions for given k-point in a states_elec_t object
subroutine zstates_elec_get_points2(st, start_point, end_point, psi)
Returns the value of all the states in the range of points [start_point:end_point].
type(type_t), public type_float
Definition: types.F90:133
type(type_t), public type_cmplx
Definition: types.F90:134
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
type(unit_t), public unit_megabytes
For large amounts of data (natural code units are bytes)
subroutine, public zwfs_elec_init(this, dim, st_start, st_end, np, ik, special, packed)
initialize an empty wfs_elec_t object
Definition: wfs_elec.F90:548
subroutine, public dwfs_elec_init(this, dim, st_start, st_end, np, ik, special, packed)
initialize an empty wfs_elec_t object
Definition: wfs_elec.F90:404
subroutine reduce_all(grp)
subroutine print_kpoints_debug
Class defining batches of mesh functions.
Definition: batch.F90:159
Distribution of N instances over mpi_grpsize processes, for the local rank mpi_grprank....
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:169
Describes mesh distribution to nodes.
Definition: mesh.F90:186
This is defined even when running serial.
Definition: mpi.F90:142
An all-pairs communication schedule for a given group.
Definition: multicomm.F90:235
Stores all communicators and groups.
Definition: multicomm.F90:206
abstract class for states
The states_elec_t class contains all electronic wave functions.
int true(void)