Octopus
lcao.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
19#include "global.h"
20
21module lcao_oct_m
22 use atom_oct_m
24 use batch_oct_m
25 use blacs_oct_m
28 use comm_oct_m
29 use debug_oct_m
31 use epot_oct_m
32 use global_oct_m
33 use grid_oct_m
36 use io_oct_m
38 use ions_oct_m
39 use, intrinsic :: iso_fortran_env
43 use lapack_oct_m
44 use loct_oct_m
46 use math_oct_m
47 use mesh_oct_m
50 use mpi_oct_m
53 use parser_oct_m
55 use ps_oct_m
59 use smear_oct_m
60 use space_oct_m
69 use unit_oct_m
71 use v_ks_oct_m
74 use xc_oct_m
75 use xc_sic_oct_m
76
77 implicit none
78
79 private
80 public :: &
81 lcao_t, &
82 lcao_init, &
84 lcao_wf, &
85 lcao_run, &
86 lcao_end, &
89
90 integer, parameter :: &
91 INITRHO_PARAMAGNETIC = 1, &
93 initrho_random = 3, &
95
96 type lcao_t
97 private
98 integer :: mode
99 logical :: complex_ylms
100 logical :: initialized
101 integer :: norbs
102 integer :: maxorbs
103 integer, allocatable :: atom(:)
104 integer, allocatable :: level(:)
105 integer, allocatable :: ddim(:)
106 logical :: alternative
107 logical :: derivative
108 integer, allocatable :: cst(:, :)
109 integer, allocatable :: ck(:, :)
110 real(4), allocatable :: dbuff_single(:, :, :, :)
111 complex(4), allocatable :: zbuff_single(:, :, :, :)
112 real(real64), allocatable :: dbuff(:, :, :, :)
113 complex(real64), allocatable :: zbuff(:, :, :, :)
114 logical :: save_memory
115 logical :: initialized_orbitals
116 real(real64) :: orbital_scale_factor
117
119 logical, allocatable :: is_empty(:)
120 ! !< This occurs in domain parallelization if the atom is not
121 ! !< in the local domain
122
124 logical :: keep_orb
125 real(real64), allocatable :: radius(:)
126 real(real64) :: lapdist
127 integer :: mult
128 integer :: maxorb
130 integer, allocatable :: basis_atom(:)
131 integer, allocatable :: basis_orb(:)
132 integer, allocatable :: atom_orb_basis(:, :)
133 integer, allocatable :: norb_atom(:)
134 logical :: parallel
135 integer :: lsize(1:2)
136 integer :: nproc(1:2)
137 integer :: myroc(1:2)
138 integer :: desc(1:BLACS_DLEN)
139 logical, allocatable :: calc_atom(:)
140 real(real64) :: diag_tol
141 type(submesh_t), allocatable :: sphere(:)
142 type(batch_t), allocatable :: orbitals(:)
143 logical, allocatable :: is_orbital_initialized(:)
144
145 integer :: gmd_opt
146 end type lcao_t
147
148
149
150 real(real64), parameter, public :: default_eigenval = 1.0e10_real64
151
152contains
153
154 ! ---------------------------------------------------------
155 subroutine lcao_init(this, namespace, space, gr, ions, st, st_start)
156 type(lcao_t), intent(out) :: this
157 type(namespace_t), intent(in) :: namespace
158 type(electron_space_t), intent(in) :: space
159 type(grid_t), intent(in) :: gr
160 type(ions_t), intent(in) :: ions
161 type(states_elec_t), intent(in) :: st
162 integer, intent(in) :: st_start
163
164 integer :: ia, n, iorb, jj, maxj, idim
165 integer :: ii, ll, mm, norbs, ii_tmp
166 integer :: mode_default
167 real(real64) :: max_orb_radius, maxradius
168 integer :: iunit_o
169
170 push_sub(lcao_init)
171
172 this%initialized = .true.
173 this%gmd_opt = initrho_userdef
174
175 ! The initial LCAO calculation is done by default if we have species representing atoms.
176 ! Otherwise, it is not the default value and has to be enforced in the input file.
177 mode_default = option__lcaostart__lcao_states
178 if (ions%only_user_def) mode_default = option__lcaostart__lcao_none
179
180 !%Variable LCAOStart
181 !%Type integer
182 !%Section SCF::LCAO
183 !%Description
184 !% Before starting a SCF calculation, <tt>Octopus</tt> can perform
185 !% a linear combination of atomic orbitals (LCAO) calculation.
186 !% These can provide <tt>Octopus</tt> with a good set
187 !% of initial wavefunctions and with a new guess for the density.
188 !% (Up to the current version, only a minimal basis set is used.)
189 !% The default is <tt>lcao_states</tt> if at least one species representing an atom is present.
190 !% The default is <tt>lcao_none</tt> if all species are <tt>species_user_defined</tt>,
191 !% <tt>species_charge_density</tt>, <tt>species_from_file</tt>, or <tt>species_jellium_slab</tt>.
192 !%
193 !% The initial guess densities for LCAO are taken from the atomic orbitals for pseudopotential species;
194 !% from the natural charge density for <tt>species_charge_density</tt>, <tt>species_point</tt>,
195 !% <tt>species_jellium</tt>, and <tt>species_jellium_slab</tt>;
196 !% or uniform for <tt>species_full_delta</tt>, <tt>species_full_gaussian</tt>,
197 !% <tt>species_user_defined</tt>, or <tt>species_from_file</tt>.
198 !% Pseudopotential species use the pseudo-wavefunctions as orbitals, full-potential atomic species
199 !% (<tt>species_full_delta</tt> and <tt>species_full_gaussian</tt>) use hydrogenic wavefunctions, and
200 !% others use harmonic-oscillator wavefunctions.
201 !%
202 !% The LCAO implementation currently does not simultaneously support k-point parallelism and states parallelism.
203 !% You must select one or the other, or alternatively use lcao_none (random wavefunctions) as an initial guess.
204 !%
205 !% Note: Some pseudopotential files (CPI, FHI for example) do not
206 !% contain full information about the orbitals. In this case,
207 !% Octopus generates the starting density from the normalized
208 !% square root of the local potential. If no orbitals are
209 !% available at all from the pseudopotential files, Octopus will
210 !% not be able to perform an LCAO and the initial states will be
211 !% randomized.
212 !%
213 !%Option lcao_none 0
214 !% Do not perform a LCAO calculation before the SCF cycle. Instead use random wavefunctions.
215 !%Option lcao_states 2
216 !% Do a LCAO calculation before the SCF cycle and use the resulting wavefunctions as
217 !% initial wavefunctions without changing the guess density.
218 !% This will speed up the convergence of the eigensolver during the first SCF iterations.
219 !%Option lcao_full 3
220 !% Do a LCAO calculation before the SCF cycle and use the LCAO wavefunctions to build a new
221 !% guess density and a new KS potential.
222 !% Using the LCAO density as a new guess density may improve the convergence, but can
223 !% also slow it down or yield wrong results (especially for spin-polarized calculations).
224 !%Option lcao_states_batch 4
225 !% (Experimental) Same as lcao_states, but faster for GPUs, large systems and parallel in states.
226 !% It is not working for spinors, however.
227 !% This is not used for unocc calculations at the moment.
228 !%End
229 call parse_variable(namespace, 'LCAOStart', mode_default, this%mode)
230 if (.not. varinfo_valid_option('LCAOStart', this%mode)) call messages_input_error(namespace, 'LCAOStart')
232 call messages_print_var_option('LCAOStart', this%mode, namespace=namespace)
234 ! LCAO alternative not implemented for st_start > 1
235 if (this%mode == option__lcaostart__lcao_states_batch .and. st_start > 1) then
236 message(1) = "LCAOStart = lcao_states_batch not compatible with this run."
237 message(2) = "Please use LCAOStart = lcao_states instead"
238 call messages_fatal(2, namespace=namespace)
239 end if
241 !TODO: Remove the alternative flag
242 this%alternative = this%mode == option__lcaostart__lcao_states_batch
243
244 if (this%mode == option__lcaostart__lcao_none) then
245 pop_sub(lcao_init)
246 return
247 end if
248
249 call messages_obsolete_variable(namespace, 'LCAOAlternative')
251
252 ! DAS: For spinors, you will always get magnetization in (1, 0, 0) direction, and the
253 ! eigenvalues will be incorrect. This is due to confusion between spins and spinors in the code.
254 if (st%d%ispin == spinors .and. this%alternative) then
255 message(1) = "LCAOStart = lcao_states_batch is not working for spinors."
256 call messages_fatal(1, namespace=namespace)
257 end if
258 if (space%is_periodic() .and. this%alternative) then
259 call messages_experimental("LCAOStart = lcao_states_batch in periodic systems", namespace=namespace)
260 ! specifically, if you get the message about submesh radius > box size, results will probably be totally wrong.
261 end if
262
263 !%Variable LCAOComplexYlms
264 !%Type logical
265 !%Default false
266 !%Section SCF::LCAO
267 !%Description
268 !% If set to true, and using complex states, complex spherical harmonics will be used, <i>i.e.</i>
269 !% with <math>e^{\pm i m \phi}</math>.
270 !% If false, real spherical harmonics with <math>\sin(m \phi)</math> or <math>\cos(m \phi)</math> are used.
271 !% This variable will make it more likely to get states that are eigenvectors of the <math>L_z</math>
272 !% operator, with a definite angular momentum.
273 !%End
274
275 if (states_are_complex(st)) then
276 call parse_variable(namespace, 'LCAOComplexYlms', .false., this%complex_ylms)
277 else
278 this%complex_ylms = .false.
279 end if
280
281 !%Variable LCAOSaveMemory
282 !%Type logical
283 !%Default false
284 !%Section SCF::LCAO
285 !%Description
286 !% If set to true, the LCAO will allocate extra memory needed in single precision instead of
287 !% double precision.
288 !%End
289 call parse_variable(namespace, 'LCAOSaveMemory', .false., this%save_memory)
290
291
292 if (debug%info .and. mpi_world%is_root()) then
293 call io_mkdir('debug/lcao', namespace)
294 iunit_o = io_open('debug/lcao/orbitals', namespace, action='write')
295 write(iunit_o,'(7a6)') 'iorb', 'atom', 'level', 'i', 'l', 'm', 'spin'
296 end if
297
298 if (.not. this%alternative) then
299
300 !%Variable LCAOScaleFactor
301 !%Type float
302 !%Default 1.0
303 !%Section SCF::LCAO
304 !%Description
305 !% The coordinates of the atomic orbitals used by the LCAO
306 !% procedure will be rescaled by the value of this variable. 1.0 means no rescaling.
307 !%End
308 call parse_variable(namespace, 'LCAOScaleFactor', m_one, this%orbital_scale_factor)
309 call messages_print_var_value('LCAOScaleFactor', this%orbital_scale_factor, namespace=namespace)
310
311 !%Variable LCAOMaximumOrbitalRadius
312 !%Type float
313 !%Default 20.0 a.u.
314 !%Section SCF::LCAO
315 !%Description
316 !% The LCAO procedure will ignore orbitals that have an
317 !% extent greater that this value.
318 !%End
319 call parse_variable(namespace, 'LCAOMaximumOrbitalRadius', 20.0_real64, max_orb_radius, unit = units_inp%length)
320 call messages_print_var_value('LCAOMaximumOrbitalRadius', max_orb_radius, units_out%length, namespace=namespace)
321
322 ! count the number of orbitals available
323 maxj = 0
324 this%maxorbs = 0
325 do ia = 1, ions%natoms
326 maxj = max(maxj, ions%atom(ia)%species%get_niwfs())
327 this%maxorbs = this%maxorbs + ions%atom(ia)%species%get_niwfs()
328 end do
329
330 this%maxorbs = this%maxorbs*st%d%dim
331
332 if (this%maxorbs == 0) then
333 call messages_write('The are no atomic orbitals available, cannot do LCAO.')
334 call messages_warning(namespace=namespace)
335 this%mode = option__lcaostart__lcao_none
336 pop_sub(lcao_init)
337 return
338 end if
339
340 ! generate tables to know which indices each atomic orbital has
341
342 safe_allocate( this%atom(1:this%maxorbs))
343 safe_allocate(this%level(1:this%maxorbs))
344 safe_allocate( this%ddim(1:this%maxorbs))
345
346 safe_allocate(this%is_empty(1:this%maxorbs))
347 this%is_empty = .false.
348
349 ! this is determined by the stencil we are using and the spacing
350 this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
351
352 ! calculate the radius of each orbital
353 safe_allocate(this%radius(1:ions%natoms))
354
355 do ia = 1, ions%natoms
356 norbs = ions%atom(ia)%species%get_niwfs()
357
358 maxradius = m_zero
359 do iorb = 1, norbs
360 call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
361 ! For all-electron species, we need to use the principal quantum number
362 if(ions%atom(ia)%species%is_full()) call ions%atom(ia)%species%get_iwf_n( iorb, 1, ii)
363 maxradius = max(maxradius, ions%atom(ia)%species%get_iwf_radius( ii, is = 1))
364 end do
365
366 maxradius = min(maxradius, m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
367
368 this%radius(ia) = maxradius
369 end do
370
371
372 ! Each atom provides niwfs pseudo-orbitals (this number is given in
373 ! ions%atom(ia)%species%niwfs for atom number ia). This number is
374 ! actually multiplied by two in case of spin-unrestricted or spinors
375 ! calculations.
376 !
377 ! The pseudo-orbitals are placed in order in the following way (Natoms
378 ! is the total number of atoms).
379 !
380 ! n = 1 => first orbital of atom 1,
381 ! n = 2 => first orbital of atom 2.
382 ! n = 3 => first orbital of atom 3.
383 ! ....
384 ! n = Natoms => first orbital of atom Natoms
385 ! n = Natoms + 1 = > second orbital of atom 1
386 ! ....
387 !
388 ! If at some point in this loop an atom pseudo cannot provide the corresponding
389 ! orbital (because the niws orbitals have been exhausted), it moves on to the following
390 ! atom.
391 !
392 ! In the spinors case, it changes a bit:
393 !
394 ! n = 1 => first spin-up orbital of atom 1, assigned to the spin-up component of the spinor.
395 ! n = 2 => first spin-down orbital of atom 1, assigned to the spin-down component of the spinor.
396 ! n = 3 => first spin-up orbital of atom 2, assigned to the spin-up component of the spinor.
397
398 iorb = 1
399 do jj = 1, maxj
400 do ia = 1, ions%natoms
401 do idim = 1,st%d%dim
402 if (jj > ions%atom(ia)%species%get_niwfs()) cycle
403 call ions%atom(ia)%species%get_iwf_ilm(jj, idim, ii, ll, mm)
404 ! For all-electron species, we need to use the principal quantum number
405 if(ions%atom(ia)%species%is_full()) then
406 ii_tmp = ii
407 call ions%atom(ia)%species%get_iwf_n( ii_tmp, 1, ii)
408 end if
409 if (this%orbital_scale_factor*ions%atom(ia)%species%get_iwf_radius( ii, is = 1) >= max_orb_radius) cycle
410
411 this%atom(iorb) = ia
412 this%level(iorb) = jj
413 this%ddim(iorb) = idim
414
415 if (debug%info .and. mpi_world%is_root()) then
416 write(iunit_o,'(7i6)') iorb, this%atom(iorb), this%level(iorb), ii, ll, mm, this%ddim(iorb)
417 end if
418
419 iorb = iorb + 1
420 end do
421 end do
422 end do
423
424 if (debug%info .and. mpi_world%is_root()) then
425 call io_close(iunit_o)
426 end if
427
428 ! some orbitals might have been removed because of their radii
429 if (this%maxorbs /= iorb - 1) then
430 call messages_write('Info: ')
431 call messages_write(this%maxorbs - iorb + 1)
432 call messages_write(' of ')
433 call messages_write(this%maxorbs)
434 call messages_write(' orbitals cannot be used for the LCAO calculation,')
435 call messages_new_line()
436 call messages_write(' their radii exceeds LCAOMaximumOrbitalRadius (')
437 call messages_write(max_orb_radius, units = units_out%length, fmt = '(f6.1)')
438 call messages_write(').')
439 call messages_warning(namespace=namespace)
440
441 this%maxorbs = iorb - 1
442 end if
443
444 if (this%maxorbs < st%nst) then
445 call messages_write('Cannot do LCAO for all states because there are not enough atomic orbitals.')
446 call messages_new_line()
447
448 call messages_write('Required: ')
449 call messages_write(st%nst)
450 call messages_write('. Available: ')
451 call messages_write(this%maxorbs)
452 call messages_write('. ')
453 call messages_write(st%nst - this%maxorbs)
454 call messages_write(' orbitals will be randomized.')
455 call messages_warning(namespace=namespace)
456 end if
457
458 !%Variable LCAODimension
459 !%Type integer
460 !%Section SCF::LCAO
461 !%Description
462 !% (Only applies if <tt>LCAOStart /= lcao_states_batch</tt>.)
463 !% Before starting the SCF cycle, an initial LCAO calculation can be performed
464 !% in order to obtain reasonable initial guesses for spin-orbitals and densities.
465 !% For this purpose, the code calculates a number of atomic orbitals.
466 !% The number available for a species described by a pseudopotential is all the
467 !% orbitals up the maximum angular momentum in the pseudopotential, minus any orbitals that
468 !% are found to be unbound. For non-pseudopotential species, the number is equal to
469 !% twice the valence charge.
470 !% The default dimension for the LCAO basis
471 !% set will be the sum of all these numbers, or twice the number of required orbitals
472 !% for the full calculation, whichever is less.
473 !%
474 !% This dimension however can be changed by making use of this
475 !% variable. Note that <tt>LCAODimension</tt> cannot be smaller than the
476 !% number of orbitals needed in the full calculation -- if
477 !% <tt>LCAODimension</tt> is smaller, it will be silently increased to meet
478 !% this requirement. In the same way, if <tt>LCAODimension</tt> is larger
479 !% than the available number of atomic orbitals, it will be
480 !% reduced. If you want to use the largest possible number, set
481 !% <tt>LCAODimension</tt> to a negative number.
482 !%End
483 call parse_variable(namespace, 'LCAODimension', 0, n)
484
485 if (n > 0 .and. n <= st%nst .and. st%nst <= this%maxorbs) then
486 ! n was made too small
487 this%norbs = st%nst
488 else if (n > st%nst .and. n <= this%maxorbs) then
489 ! n is a reasonable value
490 this%norbs = n
491 else if (n == 0) then
492 ! using the default
493 this%norbs = min(this%maxorbs, 2*st%nst)
494 else
495 ! n was negative, or greater than maxorbs
496 this%norbs = this%maxorbs
497 end if
498
499 assert(this%norbs <= this%maxorbs)
500
501 safe_allocate(this%cst(1:this%norbs, 1:st%d%spin_channels))
502 safe_allocate(this%ck(1:this%norbs, 1:st%d%spin_channels))
503 this%initialized_orbitals = .false.
504 else
505 call lcao2_init()
506 end if
507
508 pop_sub(lcao_init)
509
510 contains
511
512 subroutine lcao2_init()
513 integer :: iatom, iorb, norbs
514 real(real64) :: maxradius
515 integer :: ibasis
516#ifdef HAVE_SCALAPACK
517 integer :: jatom, jorb, jbasis, ilbasis, jlbasis, proc(1:2), info, nbl
518#endif
519 push_sub(lcao_init.lcao2_init)
520
521 call messages_write('Info: Using LCAO batched implementation.')
522 call messages_info(namespace=namespace)
523
524 call messages_experimental('LCAO alternative implementation', namespace=namespace)
525
526 !%Variable LCAOKeepOrbitals
527 !%Type logical
528 !%Default yes
529 !%Section SCF::LCAO
530 !%Description
531 !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>.
532 !% If set to yes (the default) Octopus keeps atomic orbitals in
533 !% memory during the LCAO procedure. If set to no, the orbitals
534 !% are generated each time that they are needed, increasing
535 !% computational time but saving memory.
536 !%
537 !% When set to yes, Octopus prints the amount of memory per node
538 !% that is required to store the orbitals.
539 !%
540 !%End
541 call parse_variable(namespace, 'LCAOKeepOrbitals', .true., this%keep_orb)
542
543 !%Variable LCAOExtraOrbitals
544 !%Type logical
545 !%Default false
546 !%Section SCF::LCAO
547 !%Description
548 !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>, and all species are pseudopotentials.
549 !% (experimental) If this variable is set to yes, the LCAO
550 !% procedure will add an extra set of numerical orbitals (by
551 !% using the derivative of the radial part of the original
552 !% orbitals). Note that this corresponds roughly to adding orbitals
553 !% with higher principal quantum numbers, but the same angular momentum.
554 !% This option may cause problems for unoccupied states since you may miss
555 !% some lower-lying states which correspond to higher angular momenta instead
556 !% of higher principal quantum number.
557 !%End
558 call parse_variable(namespace, 'LCAOExtraOrbitals', .false., this%derivative)
559
560 ! DAS: if you calculate the Na atom this way, spin-polarized, with just one unoccupied state,
561 ! you will obtain states (up and down) which are actually the 10th states if you start with
562 ! random wavefunctions! We really need to implement taking the derivative of the angular part
563 ! instead to be sure of getting decent results!
564
565 if (this%derivative) then
566 call messages_experimental('LCAO extra orbitals', namespace=namespace)
567
568 if (st%nst * st%smear%el_per_state > st%qtot) then
569 message(1) = "Lower-lying empty states may be missed with LCAOExtraOrbitals."
570 call messages_warning(1, namespace=namespace)
571 end if
572 end if
573
574 !%Variable LCAODiagTol
575 !%Type float
576 !%Default 1e-10
577 !%Section SCF::LCAO
578 !%Description
579 !% Only applies if <tt>LCAOStart = lcao_states_batch</tt>.
580 !% The tolerance for the diagonalization of the LCAO Hamiltonian.
581 !%End
582 call parse_variable(namespace, 'LCAODiagTol', 1e-10_real64, this%diag_tol)
583
584 if (this%derivative) then
585 this%mult = 2
586 else
587 this%mult = 1
588 end if
589
590 safe_allocate(this%sphere(1:ions%natoms))
591 safe_allocate(this%orbitals(1:ions%natoms))
592 safe_allocate(this%is_orbital_initialized(1:ions%natoms))
593 this%is_orbital_initialized = .false.
594
595 safe_allocate(this%norb_atom(1:ions%natoms))
596
597 this%maxorb = 0
598 this%norbs = 0
599 do iatom = 1, ions%natoms
600 this%norb_atom(iatom) = this%mult*ions%atom(iatom)%species%get_niwfs()
601 this%maxorb = max(this%maxorb, ions%atom(iatom)%species%get_niwfs())
602 this%norbs = this%norbs + ions%atom(iatom)%species%get_niwfs()
603 end do
604
605 this%maxorb = this%maxorb*this%mult
606 this%norbs = this%norbs*this%mult
608 safe_allocate(this%basis_atom(1:this%norbs))
609 safe_allocate(this%basis_orb(1:this%norbs))
610 safe_allocate(this%atom_orb_basis(1:ions%natoms, 1:this%maxorb))
611
612 ! Initialize the mapping between indices
613
614 ibasis = 0
615 do iatom = 1, ions%natoms
616 norbs = ions%atom(iatom)%species%get_niwfs()
617
618 do iorb = 1, this%mult*norbs
619 ibasis = ibasis + 1
620 this%atom_orb_basis(iatom, iorb) = ibasis
621 this%basis_atom(ibasis) = iatom
622 this%basis_orb(ibasis) = iorb
623
624 ! no stored spin index in alternative mode
625 if (debug%info .and. mpi_world%is_root()) then
626 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
627 write(iunit_o,'(7i6)') ibasis, iatom, iorb, ii, ll, mm, 1
628 end if
629 end do
630 end do
631
632 if (debug%info .and. mpi_world%is_root()) then
633 call io_close(iunit_o)
634 end if
635
636 ! this is determined by the stencil we are using and the spacing
637 this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
638
639 ! calculate the radius of each orbital
640 safe_allocate(this%radius(1:ions%natoms))
641
642 do iatom = 1, ions%natoms
643 norbs = ions%atom(iatom)%species%get_niwfs()
644
645 maxradius = m_zero
646 do iorb = 1, norbs
647 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
648 ! For all-electron species, we need to use the principal quantum number
649 if(ions%atom(iatom)%species%is_full()) call ions%atom(iatom)%species%get_iwf_n( iorb, 1, ii)
650 maxradius = max(maxradius, ions%atom(iatom)%species%get_iwf_radius( ii, is = 1))
651 end do
652
653 if (this%derivative) maxradius = maxradius + this%lapdist
654
655 maxradius = min(maxradius, m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
656
657 this%radius(iatom) = maxradius
658 end do
659
660 safe_allocate(this%calc_atom(1:ions%natoms))
661 this%calc_atom = .true.
662
663 ! initialize parallel data
664#ifndef HAVE_SCALAPACK
665 this%parallel = .false.
666#else
667 this%parallel = (st%parallel_in_states .or. gr%parallel_in_domains) &
668 .and. .not. blacs_proc_grid_null(st%dom_st_proc_grid)
669
670 if (this%parallel) then
671 nbl = min(16, this%norbs)
672
673 ! The size of the distributed matrix in each node
674 this%lsize(1) = max(1, numroc(this%norbs, nbl, st%dom_st_proc_grid%myrow, 0, st%dom_st_proc_grid%nprow))
675 this%lsize(2) = max(1, numroc(this%norbs, nbl, st%dom_st_proc_grid%mycol, 0, st%dom_st_proc_grid%npcol))
676
677 this%nproc(1) = st%dom_st_proc_grid%nprow
678 this%nproc(2) = st%dom_st_proc_grid%npcol
679 this%myroc(1) = st%dom_st_proc_grid%myrow
680 this%myroc(2) = st%dom_st_proc_grid%mycol
681
682 call descinit(this%desc(1), this%norbs, this%norbs, nbl, nbl, 0, 0, &
683 st%dom_st_proc_grid%context, this%lsize(1), info)
684
685 if (info /= 0) then
686 write(message(1), '(a,i6)') 'descinit for BLACS failed with error code ', info
687 call messages_fatal(1, namespace=namespace)
688 end if
689
690 this%calc_atom = .false.
691 do iatom = 1, ions%natoms
692 ibasis = this%atom_orb_basis(iatom, 1)
693
694 do jatom = 1, ions%natoms
695 jbasis = this%atom_orb_basis(jatom, 1)
696
697 do iorb = 1, this%norb_atom(iatom)
698 do jorb = 1, this%norb_atom(jatom)
699 call lcao_local_index(this, ibasis - 1 + iorb, jbasis - 1 + jorb, &
700 ilbasis, jlbasis, proc(1), proc(2))
701
702 this%calc_atom(this%basis_atom(jbasis)) = &
703 this%calc_atom(this%basis_atom(jbasis)) .or. proc(2) == this%myroc(2)
704
705 end do
706 end do
707
708 end do
709 end do
710
711 end if
712#endif
713
714 pop_sub(lcao_init.lcao2_init)
715 end subroutine lcao2_init
716
717 end subroutine lcao_init
718
719
720 ! ---------------------------------------------------------
721 subroutine lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r)
722 type(namespace_t), intent(in) :: namespace
723 type(electron_space_t), intent(in) :: space
724 type(grid_t), intent(in) :: gr
725 type(ions_t), intent(in) :: ions
726 type(partner_list_t), intent(in) :: ext_partners
727 type(states_elec_t), intent(inout) :: st
728 type(v_ks_t), intent(inout) :: ks
729 type(hamiltonian_elec_t), intent(inout) :: hm
730 integer, optional, intent(in) :: st_start
731 real(real64), optional, intent(in) :: lmm_r
732
733 type(lcao_t) :: lcao
734 integer :: st_start_random, required_min_nst
735 logical :: lcao_done
736 logical :: is_orbital_dependent
737
738 push_sub(lcao_run)
739
740 if (present(st_start)) then
741 ! If we are doing unocc calculation, do not mess with the correct eigenvalues
742 ! of the occupied states.
743 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
744 calc_eigenval=.not. present(st_start), calc_current=.false.)
745
746 assert(st_start >= 1)
747 if (st_start > st%nst) then ! nothing to be done in LCAO
748 pop_sub(lcao_run)
749 return
750 end if
751 end if
752
753 call profiling_in('LCAO_RUN')
754
755 call lcao_init(lcao, namespace, space, gr, ions, st, optional_default(st_start, 1))
756
757 call lcao_init_orbitals(lcao, namespace, st, gr, ions, start = st_start)
758
759 ! By default, we want to use vxc for the LCAO.
760 ! However, we can only do this if vxc depends on the density only
761 ! For cases like MGGA or hybrids, OEP, we cannot do this
762 is_orbital_dependent = (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
763 .or. (ks%theory_level == kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
764 .or. (ks%theory_level == generalized_kohn_sham_dft .and. xc_is_orbital_dependent(ks%xc)) &
765 .or. ks%sic%level == sic_pz_oep)
766
767 if (.not. present(st_start)) then
768 call lcao_guess_density(lcao, namespace, st, gr, hm, ions, st%qtot, st%d%ispin, st%rho)
770 if (st%d%ispin > unpolarized) then
771 assert(present(lmm_r))
772 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, lmm_r, namespace=namespace)
773 end if
774
775 ! set up Hamiltonian (we do not call v_ks_h_setup here because we do not want to
776 ! overwrite the guess density)
777 message(1) = 'Info: Setting up Hamiltonian.'
778 call messages_info(1, namespace=namespace)
779
780 ! get the effective potential (we don`t need the eigenvalues yet)
781 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.false., &
782 calc_current=.false., calc_energy=.false., force_semilocal=is_orbital_dependent)
783 ! eigenvalues have nevertheless to be initialized to something
784 ! This value must be larger that the highest eigenvalue from LCAO to get the correct occupations
785 st%eigenval = default_eigenval
786 end if
787
788 lcao_done = .false.
789
790 ! after initialized, can check that LCAO is possible
791 if (lcao_is_available(lcao)) then
792 lcao_done = .true.
793
794 if (present(st_start)) then
795 write(message(1),'(a,i8,a)') 'Performing LCAO for states ', st_start, ' and above'
796 call messages_info(1, namespace=namespace)
797 end if
798
799 call lcao_wf(lcao, st, gr, ions, hm, namespace, start = st_start)
800
801 ! In some rare cases, like bad pseudopotentials with unbound orbitals,
802 ! we might not have enough orbitals to go up to the Fermi energy
803 ! In this case, we cannot set the eigenvales to a huge value, as
804 ! the smearing will not converge. We set them to zero then
805 select case (st%d%ispin)
806 case (unpolarized)
807 required_min_nst = int(st%qtot/2)
808 case (spin_polarized)
809 required_min_nst = int(st%qtot/2)
810 case (spinors)
811 required_min_nst = int(st%qtot)
812 end select
813 if (st%smear%method /= smear_fixed_occ .and. st%smear%method /= smear_semiconductor) then
814 if (lcao%norbs <= required_min_nst .and. lcao%norbs < st%nst) then
815 st%eigenval(lcao%norbs+1:,:) = m_zero
816 end if
817 end if
818
819 if (.not. present(st_start)) then
820 call states_elec_fermi(st, namespace, gr)
821 call states_elec_write_eigenvalues(min(st%nst, lcao%norbs), st, space, hm%kpoints, namespace=namespace)
822
823 ! Update the density and the Hamiltonian
824 if (lcao%mode == option__lcaostart__lcao_full) then
825 call v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, &
826 calc_eigenval = .false., calc_current=.false.)
827 if (st%d%ispin > unpolarized) then
828 assert(present(lmm_r))
829 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, lmm_r, namespace=namespace)
830 end if
831 end if
832 end if
833 end if
834
835 if (.not. lcao_done .or. lcao%norbs < st%nst) then
836
837 if (lcao_done) then
838 st_start_random = lcao%norbs + 1
839 else
840 st_start_random = 1
841 end if
842 if (present(st_start)) st_start_random = max(st_start, st_start_random)
843
844 if (st_start_random > 1) then
845 write(message(1),'(a,i8,a)') 'Generating random wavefunctions for states ', st_start_random, ' and above'
846 call messages_info(1, namespace=namespace)
847 end if
848
849 ! Randomly generate the initial wavefunctions.
850 call states_elec_generate_random(st, gr, hm%kpoints, ist_start_ = st_start_random, normalized = .false.)
851
852 ! For spinors, we rotate the random states in the local frame of the magnetization
853 ! Unless the user specified the spin directions using InitialSpins
854 if (.not. st%fixed_spins) then
855 call rotate_random_states_to_local_frame(st, gr, st_start_random, lcao%gmd_opt)
856 end if
857
858 call messages_write('Orthogonalizing wavefunctions.')
859 call messages_info(namespace=namespace)
860 call states_elec_orthogonalize(st, namespace, gr)
861
862 if (.not. lcao_done) then
863 ! If we are doing unocc calculation, do not mess with the correct eigenvalues and occupations
864 ! of the occupied states.
865 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
866 calc_eigenval=.not. present(st_start), calc_current=.false., force_semilocal=is_orbital_dependent) ! get potentials
867 if (.not. present(st_start)) then
868 call states_elec_fermi(st, namespace, gr) ! occupations
869 end if
870
871 end if
872
873 else if (present(st_start)) then
874
875 if (st_start > 1) then
876 call messages_write('Orthogonalizing wavefunctions.')
877 call messages_info(namespace=namespace)
878 call states_elec_orthogonalize(st, namespace, gr)
879 end if
880
881 end if
882
883 call lcao_end(lcao)
884
885
886 call profiling_out('LCAO_RUN')
887 pop_sub(lcao_run)
888 end subroutine lcao_run
889
890 ! ---------------------------------------------------------
891 subroutine lcao_end(this)
892 type(lcao_t), intent(inout) :: this
893
894 push_sub(lcao_end)
895
896 safe_deallocate_a(this%calc_atom)
897 safe_deallocate_a(this%norb_atom)
898 safe_deallocate_a(this%basis_atom)
899 safe_deallocate_a(this%basis_orb)
900 safe_deallocate_a(this%atom_orb_basis)
901 safe_deallocate_a(this%radius)
902 safe_deallocate_a(this%sphere)
903 safe_deallocate_a(this%orbitals)
904
905 safe_deallocate_a(this%atom)
906 safe_deallocate_a(this%level)
907 safe_deallocate_a(this%ddim)
908 safe_deallocate_a(this%cst)
909 safe_deallocate_a(this%ck)
910 safe_deallocate_a(this%dbuff_single)
911 safe_deallocate_a(this%zbuff_single)
912 safe_deallocate_a(this%dbuff)
913 safe_deallocate_a(this%zbuff)
914
915 this%initialized = .false.
916 pop_sub(lcao_end)
917 end subroutine lcao_end
918
919
920 ! ---------------------------------------------------------
921 subroutine lcao_wf(this, st, gr, ions, hm, namespace, start)
922 type(lcao_t), intent(inout) :: this
923 type(states_elec_t), intent(inout) :: st
924 type(grid_t), intent(in) :: gr
925 type(ions_t), intent(in) :: ions
926 type(hamiltonian_elec_t), intent(in) :: hm
927 type(namespace_t), intent(in) :: namespace
928 integer, optional, intent(in) :: start
929
930 integer :: start_
931
932 assert(this%initialized)
933
934 call profiling_in("LCAO")
935 push_sub(lcao_wf)
936
937 start_ = 1
938 if (present(start)) start_ = start
940 if (this%alternative) then
941 if (states_are_real(st)) then
942 call dlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
943 else
944 call zlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
945 end if
946 else
947 if (states_are_real(st)) then
948 call dlcao_wf(this, st, gr, ions, hm, namespace, start_)
949 else
950 call zlcao_wf(this, st, gr, ions, hm, namespace, start_)
951 end if
952 end if
953 pop_sub(lcao_wf)
954 call profiling_out("LCAO")
955 end subroutine lcao_wf
956
957
958 ! ---------------------------------------------------------
960 logical function lcao_is_available(this) result(available)
961 type(lcao_t), intent(in) :: this
962
963 push_sub(lcao_is_available)
964
965 available = this%initialized .and. this%mode /= option__lcaostart__lcao_none &
966 .and. this%norbs > 0
967
968 pop_sub(lcao_is_available)
969 end function lcao_is_available
970
971
972 ! ---------------------------------------------------------
974 integer function lcao_num_orbitals(this) result(norbs)
975 type(lcao_t), intent(in) :: this
976
977 push_sub(lcao_num_orbitals)
978 norbs = this%norbs
979
980 pop_sub(lcao_num_orbitals)
981 end function lcao_num_orbitals
982
983 ! ---------------------------------------------------------
984
985 subroutine lcao_local_index(this, ig, jg, il, jl, prow, pcol)
986 type(lcao_t), intent(in) :: this
987 integer, intent(in) :: ig
988 integer, intent(in) :: jg
989 integer, intent(out) :: il
990 integer, intent(out) :: jl
991 integer, intent(out) :: prow
992 integer, intent(out) :: pcol
993
994 ! no PUSH_SUB, called too often
995#ifdef HAVE_SCALAPACK
996 call infog2l(ig, jg, this%desc(1), this%nproc(1), this%nproc(2), this%myroc(1), this%myroc(2), &
997 il, jl, prow, pcol)
998#else
999 il = ig
1000 jl = jg
1001 prow = 0
1002 pcol = 0
1003#endif
1004
1005 end subroutine lcao_local_index
1006
1007 ! ---------------------------------------------------------
1012 subroutine lcao_alt_end_orbital(this, iatom)
1013 type(lcao_t), intent(inout) :: this
1014 integer, intent(in) :: iatom
1015
1016 push_sub(lcao_alt_end_orbital)
1017
1018 if (this%is_orbital_initialized(iatom)) then
1019 call this%orbitals(iatom)%end()
1020 this%is_orbital_initialized(iatom) = .false.
1021 end if
1023 pop_sub(lcao_alt_end_orbital)
1024
1025 end subroutine lcao_alt_end_orbital
1026
1027 ! ---------------------------------------------------------
1028
1029 subroutine lcao_atom_density(this, st, mesh, ions, iatom, spin_channels, rho)
1030 type(lcao_t), intent(inout) :: this
1031 type(states_elec_t), intent(in) :: st
1032 class(mesh_t), intent(in) :: mesh
1033 type(ions_t), target, intent(in) :: ions
1034 integer, intent(in) :: iatom
1035 integer, intent(in) :: spin_channels
1036 real(real64), intent(inout) :: rho(:, :)
1037
1038 real(real64), allocatable :: dorbital(:, :)
1039 complex(real64), allocatable :: zorbital(:, :)
1040 real(real64), allocatable :: factors(:)
1041 real(real64) :: factor, aa
1042 integer :: iorb, ip, ii, ll, mm, ispin
1043 type(ps_t), pointer :: ps
1044 logical :: use_stored_orbitals
1045
1046 push_sub(lcao_atom_density)
1047
1048 rho = m_zero
1049
1050 use_stored_orbitals = ions%atom(iatom)%species%is_ps() &
1051 .and. states_are_real(st) .and. spin_channels == 1 .and. lcao_is_available(this) &
1052 .and. st%d%dim == 1 .and. .not. ions%space%is_periodic()
1053
1054
1055 ! we can use the orbitals we already have calculated
1056 if (use_stored_orbitals) then
1057 !There is no periodic copies here, so this will not work for periodic systems
1058 assert(.not. ions%space%is_periodic())
1059
1060 select type(spec=>ions%atom(iatom)%species)
1061 class is(pseudopotential_t)
1062 ps => spec%ps
1063 class default
1064 assert(.false.)
1065 end select
1066
1067 if (.not. this%alternative) then
1068
1069 if (states_are_real(st)) then
1070 safe_allocate(dorbital(1:mesh%np, 1:st%d%dim))
1071 else
1072 safe_allocate(zorbital(1:mesh%np, 1:st%d%dim))
1073 end if
1074
1075 do iorb = 1, this%norbs
1076 if (iatom /= this%atom(iorb)) cycle
1078 call ions%atom(iatom)%species%get_iwf_ilm(this%level(iorb), 1, ii, ll, mm)
1079 factor = ps%conf%occ(ii, 1)/(m_two*ll + m_one)
1080
1081 if (states_are_real(st)) then
1082 call dget_ao(this, st, mesh, ions, iorb, 1, dorbital, use_psi = .true.)
1083 !$omp parallel do
1084 do ip = 1, mesh%np
1085 rho(ip, 1) = rho(ip, 1) + factor*dorbital(ip, 1)**2
1086 end do
1087 else
1088 call zget_ao(this, st, mesh, ions, iorb, 1, zorbital, use_psi = .true.)
1089 !$omp parallel do
1090 do ip = 1, mesh%np
1091 rho(ip, 1) = rho(ip, 1) + factor*abs(zorbital(ip, 1))**2
1092 end do
1093 end if
1094
1095 end do
1096
1097 safe_deallocate_a(dorbital)
1098 safe_deallocate_a(zorbital)
1099
1100 else
1101
1102 ! for simplicity, always use real ones here.
1103 call dlcao_alt_get_orbital(this, this%sphere(iatom), ions, 1, iatom, this%norb_atom(iatom))
1104
1105 ! the extra orbitals with the derivative are not relevant here, hence we divide by this%mult
1106 safe_allocate(factors(1:this%norb_atom(iatom)/this%mult))
1107
1108 do iorb = 1, this%norb_atom(iatom)/this%mult
1109 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
1110 factors(iorb) = ps%conf%occ(ii, 1)/(m_two*ll + m_one)
1111 end do
1112
1113 !$omp parallel do private(ip, aa, iorb) if(.not. this%sphere(iatom)%overlap)
1114 do ip = 1, this%sphere(iatom)%np
1115 aa = m_zero
1116 do iorb = 1, this%norb_atom(iatom)/this%mult
1117 aa = aa + factors(iorb)*this%orbitals(iatom)%dff_linear(ip, iorb)**2
1118 end do
1119 rho(this%sphere(iatom)%map(ip), 1) = rho(this%sphere(iatom)%map(ip), 1) + aa
1120 end do
1121
1122 safe_deallocate_a(factors)
1123
1124 end if
1125
1126 else
1127 call species_atom_density(ions%atom(iatom)%species, ions%namespace, ions%space, ions%latt, &
1128 ions%pos(:, iatom), mesh, spin_channels, rho)
1129 end if
1130
1131 ! The above code can sometimes return negative values of the density. Here we avoid introducing
1132 ! them in the calculation of v_s, mixing, ...
1133 do ispin = 1, spin_channels
1134 !$omp parallel do simd
1135 do ip = 1, mesh%np
1136 rho(ip, ispin) = max(rho(ip, ispin), m_zero)
1137 end do
1138 end do
1139
1140 pop_sub(lcao_atom_density)
1141 end subroutine lcao_atom_density
1142
1143 ! ---------------------------------------------------------
1145 subroutine lcao_guess_density(this, namespace, st, gr, hm, ions, qtot, ispin, rho)
1146 type(lcao_t), intent(inout) :: this
1147 type(namespace_t), intent(in) :: namespace
1148 type(states_elec_t), intent(in) :: st
1149 type(grid_t), intent(in) :: gr
1150 type(hamiltonian_elec_t), intent(in) :: hm
1151 type(ions_t), intent(in) :: ions
1152 real(real64), intent(in) :: qtot
1153 integer, intent(in) :: ispin
1154 real(real64), contiguous, intent(out) :: rho(:, :)
1155
1156 integer :: ia, is, idir, ip, m_dim
1157 integer(int64), save :: iseed = splitmix64_321
1158 type(block_t) :: blk
1159 real(real64) :: rr, rnd, phi, theta, lmag, n1, n2
1160 real(real64), allocatable :: atom_rho(:,:), mag(:,:)
1161 real(real64), parameter :: tol_min_mag = 1.0e-20_real64
1162
1163 push_sub(lcao_guess_density)
1164
1165 if (st%d%spin_channels == 1) then
1166 this%gmd_opt = initrho_paramagnetic
1167 else
1168 !%Variable GuessMagnetDensity
1169 !%Type integer
1170 !%Default ferromagnetic
1171 !%Section SCF::LCAO
1172 !%Description
1173 !% The guess density for the SCF cycle is just the sum of all the atomic densities.
1174 !% When performing spin-polarized or non-collinear-spin calculations this option sets
1175 !% the guess magnetization density.
1176 !%
1177 !% For anti-ferromagnetic configurations, the <tt>user_defined</tt> option should be used.
1178 !%
1179 !% Note that if the <tt>paramagnetic</tt> option is used, the final ground state will also be
1180 !% paramagnetic, but the same is not true for the other options.
1181 !%Option paramagnetic 1
1182 !% Magnetization density is zero.
1183 !%Option ferromagnetic 2
1184 !% Magnetization density is the sum of the atomic magnetization densities.
1185 !%Option random 3
1186 !% Each atomic magnetization density is randomly rotated.
1187 !%Option user_defined 77
1188 !% The atomic magnetization densities are rotated so that the magnetization
1189 !% vector has the same direction as a vector provided by the user. In this case,
1190 !% the <tt>AtomsMagnetDirection</tt> block has to be set.
1191 !%End
1192 call parse_variable(namespace, 'GuessMagnetDensity', initrho_ferromagnetic, this%gmd_opt)
1193 if (.not. varinfo_valid_option('GuessMagnetDensity', this%gmd_opt)) call messages_input_error(namespace, 'GuessMagnetDensity')
1194 call messages_print_var_option('GuessMagnetDensity', this%gmd_opt, namespace=namespace)
1195 end if
1196
1197 if (parse_is_defined(namespace, 'GuessMagnetDensity') .and. (hm%theory_level == hartree_fock &
1198 .or. hm%theory_level == generalized_kohn_sham_dft) .and. st%d%spin_channels > 1) then
1199 message(1) = "GuessMagnetDensity cannot be used for Hartree-Fock and generalized Kohn-Sham calculation."
1200 message(2) = "Please perform a LDA or GGA calculation first and restart from this calculation."
1201 call messages_fatal(2, namespace=namespace)
1202 end if
1203
1204 if (this%gmd_opt == initrho_userdef) then
1205 !%Variable AtomsMagnetDirection
1206 !%Type block
1207 !%Section SCF::LCAO
1208 !%Description
1209 !% This option is only used when <tt>GuessMagnetDensity</tt> is
1210 !% set to <tt>user_defined</tt>. It provides a direction for the
1211 !% magnetization vector of each atom when building the guess
1212 !% density. In order to do that, the user should specify the
1213 !% coordinates of a vector that has the desired direction and
1214 !% norm. Note that it is necessary to maintain the ordering in
1215 !% which the species were defined in the coordinates
1216 !% specifications.
1217 !%
1218 !% For spin-polarized calculations, the vectors should have only
1219 !% one component; for non-collinear-spin calculations, they
1220 !% should have three components. If the norm of the vector is greater
1221 !% than the number of valence electrons in the atom, it will be rescaled
1222 !% to this number, which is the maximum possible magnetization.
1223 !%End
1224 if (parse_block(namespace, 'AtomsMagnetDirection', blk) < 0) then
1225 message(1) = "AtomsMagnetDirection block is not defined."
1226 call messages_fatal(1, namespace=namespace)
1227 end if
1228
1229 if (parse_block_n(blk) /= ions%natoms) then
1230 message(1) = "The number of rows in the AtomsMagnetDirection block does not equal the number of atoms."
1231 call messages_fatal(1, namespace=namespace)
1232 end if
1233
1234 if (ispin == spin_polarized) then
1235 m_dim = 1
1236 elseif(ispin == spinors) then
1237 m_dim = 3
1238 endif
1239
1240 safe_allocate(mag(1:m_dim, 1:ions%natoms))
1241 do ia = 1, ions%natoms
1242 !Read from AtomsMagnetDirection block
1243 do idir = 1, m_dim
1244 call parse_block_float(blk, ia-1, idir-1, mag(idir, ia))
1245 if (abs(mag(idir, ia)) < tol_min_mag) mag(idir, ia) = m_zero
1246 end do
1247 end do
1248 call parse_block_end(blk)
1249 end if
1250
1251 rho = m_zero
1252
1253 safe_allocate(atom_rho(1:gr%np, 1:st%d%spin_channels))
1254 select case (this%gmd_opt)
1255 case (initrho_paramagnetic)
1256 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1257 call lcao_atom_density(this, st, gr, ions, ia, st%d%spin_channels, atom_rho)
1258 call lalg_axpy(gr%np, st%d%spin_channels, m_one, atom_rho, rho)
1259 end do
1260
1261 if (st%d%spin_channels == 2) then
1262 !$omp parallel do
1263 do ip = 1, gr%np
1264 rho(ip, 1) = m_half*(rho(ip, 1) + rho(ip, 2))
1265 rho(ip, 2) = rho(ip, 1)
1266 end do
1267 end if
1268
1270 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1271 call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho(1:gr%np, 1:2))
1272 rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
1273 end do
1274
1275 case (initrho_random) ! Randomly oriented spins
1276 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1277 call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho)
1278
1279 ! For charge neutral atoms, the magnetization density will always be zero
1280 ! In order to still make a random spin structure, we then reassign the charge in
1281 ! one spin channel to have a 1 \nu_B for this atom, with a random direction
1282 !
1283 ! An example where this is needed is a Xe3 cluster where we put an excess charge
1284 ! Each atom is neutral but the full system has a net magnetiation
1285 n1 = dmf_integrate(gr, atom_rho(:, 1))
1286 n2 = dmf_integrate(gr, atom_rho(:, 2))
1287
1288 if (is_close(n1, n2)) then
1289 lmag = m_one
1290
1291 call lalg_axpy(gr%np, (lmag - n1 + n2)/m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
1292 call lalg_scal(gr%np, (n1 + n2 - lmag)/m_two/n2, atom_rho(:, 2))
1293 end if
1294
1295
1296 if (ispin == spin_polarized) then
1297 call quickrnd(iseed, rnd)
1298 rnd = rnd - m_half
1299 if (rnd > m_zero) then
1300 rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
1301 else
1302 rho(1:gr%np, 1) = rho(1:gr%np, 1) + atom_rho(1:gr%np, 2)
1303 rho(1:gr%np, 2) = rho(1:gr%np, 2) + atom_rho(1:gr%np, 1)
1304 end if
1305 elseif (ispin == spinors) then
1306 call quickrnd(iseed, phi)
1307 call quickrnd(iseed, theta)
1308 phi = phi*m_two*m_pi
1309 theta = theta*m_pi*m_half
1310
1311 call accumulate_rotated_density(gr, rho, atom_rho, theta, phi)
1312 end if
1313 end do
1314
1315 case (initrho_userdef) ! User-defined
1316 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1317 !Get atomic density
1318 call lcao_atom_density(this, st, gr, ions, ia, 2, atom_rho)
1319
1320 !Scale magnetization density
1321 n1 = dmf_integrate(gr, atom_rho(:, 1))
1322 n2 = dmf_integrate(gr, atom_rho(:, 2))
1323
1324 lmag = norm2(mag(:, ia))
1325 if (lmag > n1 + n2) then
1326 mag = mag*(n1 + n2)/lmag
1327 lmag = n1 + n2
1328 elseif (abs(lmag) <= m_epsilon) then
1329 if (abs(n1 - n2) <= m_epsilon) then
1330 call lalg_axpy(gr%np, 2, m_one, atom_rho, rho)
1331 else
1332 !$omp parallel do simd
1333 do ip = 1, gr%np
1334 atom_rho(ip, 1) = m_half*(atom_rho(ip, 1) + atom_rho(ip, 2))
1335 rho(ip, 1) = rho(ip, 1) + atom_rho(ip, 1)
1336 rho(ip, 2) = rho(ip, 2) + atom_rho(ip, 1)
1337 end do
1338 end if
1339 cycle
1340 end if
1341
1342 if (.not. is_close(n1 - n2, lmag) .and. abs(n2) > m_epsilon) then
1343 if (n1 - n2 < lmag) then
1344 call lalg_axpy(gr%np, (lmag - n1 + n2)/m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
1345 call lalg_scal(gr%np, (n1 + n2 - lmag)/m_two/n2, atom_rho(:, 2))
1346 elseif (n1 - n2 > lmag) then
1347 call lalg_axpy(gr%np, (n1 - n2 - lmag)/m_two/n1, atom_rho(:, 1), atom_rho(:, 2))
1348 call lalg_scal(gr%np, (n1 + n2 + lmag)/m_two/n1, atom_rho(:, 1))
1349 end if
1350 end if
1351
1352 !Rotate magnetization density
1353 if (ispin == spin_polarized) then
1354
1355 if (mag(1, ia) > m_zero) then
1356 call lalg_axpy(gr%np, 2, m_one, atom_rho, rho)
1357 else
1358 call lalg_axpy(gr%np, m_one, atom_rho(:,2), rho(:,1))
1359 call lalg_axpy(gr%np, m_one, atom_rho(:,1), rho(:,2))
1360 end if
1361
1362 elseif (ispin == spinors) then
1363
1364 call get_angles_from_magnetization(mag(:, ia), lmag, theta, phi)
1365 call accumulate_rotated_density(gr, rho, atom_rho, theta, phi)
1366 end if
1367 end do
1368
1369 end select
1370
1371
1372 if (ions%atoms_dist%parallel) then
1373 do is = 1, st%d%nspin
1374 call lalg_copy(gr%np, rho(:,is), atom_rho(:,1))
1375 call ions%atoms_dist%mpi_grp%allreduce(atom_rho(1, 1), rho(1, is), gr%np, mpi_double_precision, mpi_sum)
1376 end do
1377 end if
1378
1379 ! we now renormalize the density (necessary if we have a charged system)
1380 rr = integrated_charge_density(gr, st, rho)
1381 write(message(1),'(a,f13.6)')'Info: Unnormalized total charge = ', rr
1382 call messages_info(1, namespace=namespace)
1383
1384 ! We only renormalize if the density is not zero
1385 if (abs(rr) > m_epsilon) then
1386 call lalg_scal(gr%np, st%d%nspin, qtot/rr, rho)
1387 rr = integrated_charge_density(gr, st, rho)
1388 write(message(1),'(a,f13.6)')'Info: Renormalized total charge = ', rr
1389 call messages_info(1, namespace=namespace)
1390 end if
1391
1392 ! Symmetrize the density if needed
1393 if (st%symmetrize_density) then
1394 do is = 1, st%d%nspin
1395 call dgrid_symmetrize_scalar_field(gr, rho(:, is))
1396 end do
1397 end if
1398
1399 safe_deallocate_a(atom_rho)
1400 safe_deallocate_a(mag)
1401
1402 pop_sub(lcao_guess_density)
1403 end subroutine lcao_guess_density
1404
1406 real(real64) function integrated_charge_density(gr, st, rho) result(rr)
1407 type(grid_t), intent(in) :: gr
1408 type(states_elec_t), intent(in) :: st
1409 real(real64), intent(in) :: rho(:,:)
1410
1411 integer :: is
1412
1413 rr = m_zero
1414 do is = 1, st%d%spin_channels
1415 rr = rr + dmf_integrate(gr, rho(:, is), reduce = .false.)
1416 end do
1417
1418 call gr%allreduce(rr)
1419 end function integrated_charge_density
1420
1421 ! ---------------------------------------------------------
1422 subroutine accumulate_rotated_density(mesh, rho, atom_rho, theta, phi)
1423 class(mesh_t), intent(in) :: mesh
1424 real(real64), intent(inout) :: rho(:,:)
1425 real(real64), intent(in) :: atom_rho(:,:)
1426 real(real64), intent(in) :: theta, phi
1427
1428 integer :: ip
1429 real(real64) :: half_theta
1430
1432
1433 half_theta = m_half * theta
1434
1435 !$omp parallel do simd
1436 do ip = 1, mesh%np
1437 rho(ip, 1) = rho(ip, 1) + cos(half_theta)**2*atom_rho(ip, 1) + sin(half_theta)**2*atom_rho(ip, 2)
1438 rho(ip, 2) = rho(ip, 2) + sin(half_theta)**2*atom_rho(ip, 1) + cos(half_theta)**2*atom_rho(ip, 2)
1439 rho(ip, 3) = rho(ip, 3) + cos(half_theta)*sin(half_theta)*cos(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
1440 rho(ip, 4) = rho(ip, 4) - cos(half_theta)*sin(half_theta)*sin(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
1441 end do
1442
1444 end subroutine accumulate_rotated_density
1445
1446 ! ---------------------------------------------------------
1447
1448 subroutine lcao_init_orbitals(this, namespace, st, gr, ions, start)
1449 type(lcao_t), intent(inout) :: this
1450 type(namespace_t), intent(in) :: namespace
1451 type(states_elec_t), intent(inout) :: st
1452 type(grid_t), intent(in) :: gr
1453 type(ions_t), intent(in) :: ions
1454 integer, optional, intent(in) :: start
1455
1456 if (.not. lcao_is_available(this)) return
1457
1458 push_sub(lcao_init_orbitals)
1459
1460 if (.not. this%alternative) then
1461 if (states_are_real(st)) then
1462 call dinit_orbitals(this, namespace, st, gr, ions, start)
1463 else
1464 call zinit_orbitals(this, namespace, st, gr, ions, start)
1465 end if
1466 else
1467 if (states_are_real(st)) then
1468 call dlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
1469 else
1470 call zlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
1471 end if
1472
1473 end if
1474
1475 pop_sub(lcao_init_orbitals)
1476 end subroutine lcao_init_orbitals
1477
1479 subroutine get_angles_from_magnetization(mag, lmag, theta, phi)
1480 real(real64), intent(in) :: mag(3)
1481 real(real64), intent(in) :: lmag
1482 real(real64), intent(out) :: theta
1483 real(real64), intent(out) :: phi
1484
1485 assert(lmag > m_zero)
1486
1487 theta = acos(max(-1.0_real64, min(1.0_real64, mag(3)/lmag)))
1488
1489 if (abs(mag(1)) <= m_epsilon) then
1490 if (abs(mag(2)) <= m_epsilon) then
1491 phi = m_zero
1492 elseif (mag(2) < m_zero) then
1493 phi = m_pi*m_three*m_half
1494 else
1495 phi = m_pi*m_half
1496 end if
1497 else
1498
1499 phi = atan2(mag(2), mag(1))
1500 !note: atan2 returns the correct phase. The correction from wikipedia on spherical coordinates is not necessary.
1501
1502 end if
1503
1504 end subroutine get_angles_from_magnetization
1505
1506
1513 subroutine rotate_random_states_to_local_frame(st, gr, ist_start, gmd_opt)
1514 type(states_elec_t), intent(inout) :: st
1515 type(grid_t), intent(in) :: gr
1516 integer, intent(in) :: ist_start
1517 integer, intent(in) :: gmd_opt
1518
1519 integer :: ik, ist, ip
1520 real(real64), allocatable :: md(:,:), up(:)
1521 complex(real64), allocatable :: down(:), psi(:,:)
1522 real(real64) :: lmag, theta, phi, norm, mm(3)
1523
1524 if (st%d%ispin /= spinors) return
1525 if (gmd_opt == initrho_random) return
1526
1528
1529 safe_allocate(psi(1:gr%np, 1:st%d%dim))
1530
1531 select case (gmd_opt)
1533 ! Random spinors should have no net magnetization
1534 do ik = st%d%kpt%start, st%d%kpt%end
1535 do ist = ist_start, st%st_end
1536 call states_elec_get_state(st, gr, ist, ik, psi)
1537 !$omp parallel do private(norm)
1538 do ip = 1, gr%np
1539 psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
1540 psi(ip, 2) = psi(ip, 1)
1541 end do
1542 call zmf_normalize(gr, st%d%dim, psi)
1543 call states_elec_set_state(st, gr, ist, ik, psi)
1544 end do
1545 end do
1546
1548 ! Random spinors magnetization is aligned along the z axis
1549 do ik = st%d%kpt%start, st%d%kpt%end
1550 do ist = ist_start, st%st_end
1551 call states_elec_get_state(st, gr, ist, ik, psi)
1552 !$omp parallel do
1553 do ip = 1, gr%np
1554 psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
1555 psi(ip, 2) = m_zero
1556 end do
1557 call zmf_normalize(gr, st%d%dim, psi)
1558 call states_elec_set_state(st, gr, ist, ik, psi)
1559 end do
1560 end do
1562 case(initrho_userdef)
1563 ! Random spinors magnetization is aligned along the local frame
1564 safe_allocate(md(1:gr%np, 1:3))
1565 safe_allocate(up(1:gr%np))
1566 safe_allocate(down(1:gr%np))
1567 call magnetic_density(gr, st%d, st%rho, md)
1568
1569 mm(1) = dmf_integrate(gr, md(:, 1), reduce = .false.)
1570 mm(2) = dmf_integrate(gr, md(:, 2), reduce = .false.)
1571 mm(3) = dmf_integrate(gr, md(:, 3), reduce = .false.)
1572
1573 if (gr%parallel_in_domains) then
1574 call gr%allreduce(mm)
1575 end if
1576
1577 ! If the density has a net direction, we align globally the spinors to this direction
1578 ! This allows to generate random pure states
1579 lmag = norm2(mm)
1580 if (lmag > 1.0e-2_real64) then
1581 call get_angles_from_magnetization(mm, lmag, theta, phi)
1582 up(:) = cos(theta/m_two)
1583 down(:) = exp(-m_zi * phi) * sin(theta/m_two)
1584 else
1585 ! We do a rotation in the local frame
1586 ! Note that this produces mixed states, so the norm of the magnetization is smaller than 1
1587 !
1588 !$omp parallel do private(lmag, theta, phi)
1589 do ip = 1, gr%np
1590 lmag = norm2(md(ip, 1:3))
1591 call get_angles_from_magnetization(md(ip,:), lmag, theta, phi)
1592 up(ip) = cos(theta/m_two)
1593 down(ip) = exp(-m_zi * phi) * sin(theta/m_two)
1594 end do
1595 safe_deallocate_a(md)
1596 end if
1597
1598 ! Rotate the spinors
1599 do ik = st%d%kpt%start, st%d%kpt%end
1600 do ist = ist_start, st%st_end
1601 call states_elec_get_state(st, gr, ist, ik, psi)
1602 call zmf_normalize(gr, st%d%dim, psi)
1603 !$omp parallel do private(norm)
1604 do ip = 1, gr%np
1605 norm = sqrt(real(psi(ip, 1)*conjg(psi(ip, 1)) + psi(ip, 2)*conjg(psi(ip, 2))))
1606 psi(ip, 1) = norm * up(ip)
1607 psi(ip, 2) = norm * down(ip)
1608 end do
1609 call states_elec_set_state(st, gr, ist, ik, psi)
1610 end do
1611 end do
1612 safe_deallocate_a(up)
1613 safe_deallocate_a(down)
1614
1615 end select
1616 safe_deallocate_a(psi)
1617
1620
1621#include "undef.F90"
1622#include "real.F90"
1623#include "lcao_inc.F90"
1624
1625#include "undef.F90"
1626#include "complex.F90"
1627#include "lcao_inc.F90"
1628
1629
1630end module lcao_oct_m
1631
1632!! Local Variables:
1633!! mode: f90
1634!! coding: utf-8
1635!! End:
subroutine info()
Definition: em_resp.F90:1098
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
scales a vector by a constant
Definition: lalg_basic.F90:159
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
double acos(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
subroutine lcao2_init()
Definition: lcao.F90:608
This module implements batches of mesh functions.
Definition: batch.F90:135
This module contains interfaces for BLACS routines Interfaces are from http:
Definition: blacs.F90:27
This module provides the BLACS processor grid.
logical pure function, public blacs_proc_grid_null(this)
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
type(debug_t), save, public debug
Definition: debug.F90:158
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
Definition: global.F90:193
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
integer, parameter, public hartree_fock
Definition: global.F90:237
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:237
integer, parameter, public kohn_sham_dft
Definition: global.F90:237
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
integer, parameter, public hartree
Definition: global.F90:237
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public dgrid_symmetrize_scalar_field(gr, field, suppress_warning)
Definition: grid.F90:670
This module defines classes and functions for interaction partners.
Definition: io.F90:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:360
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
A module to handle KS potential, without the external potential.
This module contains interfaces for LAPACK routines.
Definition: lapack.F90:120
subroutine, public lcao_init_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:1497
subroutine, public lcao_wf(this, st, gr, ions, hm, namespace, start)
Definition: lcao.F90:970
subroutine zlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:3333
subroutine zget_ao(this, st, mesh, ions, iorb, ispin, ao, use_psi)
Definition: lcao.F90:3301
subroutine, public lcao_end(this)
Definition: lcao.F90:940
subroutine, public lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r)
Definition: lcao.F90:770
subroutine zlcao_alt_wf(this, st, gr, ions, hm, namespace, start)
The alternative implementation.
Definition: lcao.F90:3386
subroutine dlcao_wf(this, st, gr, ions, hm, namespace, start)
Definition: lcao.F90:1810
subroutine dlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:2211
subroutine zlcao_wf(this, st, gr, ions, hm, namespace, start)
Definition: lcao.F90:2924
subroutine rotate_random_states_to_local_frame(st, gr, ist_start, gmd_opt)
Rotate the spinors with band index >= ist_start to the local frame of the magnetization.
Definition: lcao.F90:1562
subroutine dlcao_alt_wf(this, st, gr, ions, hm, namespace, start)
The alternative implementation.
Definition: lcao.F90:2264
integer, parameter initrho_userdef
Definition: lcao.F90:185
subroutine get_angles_from_magnetization(mag, lmag, theta, phi)
Given a magnetization vector, and its norm, this returns the corresponding inclination and azimuthal ...
Definition: lcao.F90:1528
integer, parameter initrho_random
Definition: lcao.F90:185
subroutine lcao_alt_end_orbital(this, iatom)
This function deallocates a set of an atomic orbitals for an atom. It can be called when the batch is...
Definition: lcao.F90:1061
subroutine lcao_local_index(this, ig, jg, il, jl, prow, pcol)
Definition: lcao.F90:1034
integer function, public lcao_num_orbitals(this)
Returns the number of LCAO orbitas.
Definition: lcao.F90:1023
subroutine lcao_guess_density(this, namespace, st, gr, hm, ions, qtot, ispin, rho)
builds a density which is the sum of the atomic densities
Definition: lcao.F90:1194
subroutine dinit_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:2070
real(real64) function integrated_charge_density(gr, st, rho)
Computes the integral of rho, summed over spin channels.
Definition: lcao.F90:1455
integer, parameter initrho_paramagnetic
Definition: lcao.F90:185
subroutine dget_ao(this, st, mesh, ions, iorb, ispin, ao, use_psi)
Definition: lcao.F90:2179
subroutine dlcao_alt_get_orbital(this, sphere, ions, ispin, iatom, norbs)
This function generates the set of an atomic orbitals for an atom and stores it in the batch orbitalb...
Definition: lcao.F90:2728
subroutine zinit_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:3192
subroutine lcao_atom_density(this, st, mesh, ions, iatom, spin_channels, rho)
Definition: lcao.F90:1078
subroutine, public lcao_init(this, namespace, space, gr, ions, st, st_start)
Definition: lcao.F90:251
subroutine accumulate_rotated_density(mesh, rho, atom_rho, theta, phi)
Definition: lcao.F90:1471
integer, parameter initrho_ferromagnetic
Definition: lcao.F90:185
logical function, public lcao_is_available(this)
Returns true if LCAO can be done.
Definition: lcao.F90:1009
subroutine, public write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
Definition: magnetic.F90:207
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
subroutine, public messages_new_line()
Definition: messages.F90:1112
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:455
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
Definition: ps.F90:116
integer(int64), parameter, public splitmix64_321
Definition: quickrnd.F90:137
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
Definition: scalapack.F90:133
integer, parameter, public smear_semiconductor
Definition: smear.F90:173
integer, parameter, public smear_fixed_occ
Definition: smear.F90:173
subroutine, public species_atom_density(species, namespace, space, latt, pos, mesh, spin_channels, rho)
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
subroutine, public states_elec_orthogonalize(st, namespace, mesh)
Orthonormalizes nst orbitals in mesh (honours state parallelization).
This module defines routines to write information about states.
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
subroutine, public states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
randomize states
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:749
subroutine, public v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
Definition: v_ks.F90:695
Definition: xc.F90:116
logical pure function, public xc_is_orbital_dependent(xcs)
Is the xc family orbital dependent.
Definition: xc.F90:590
integer, parameter, public sic_pz_oep
Perdew-Zunger SIC (OEP way)
Definition: xc_sic.F90:151
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:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
A type storing the information and data about a pseudopotential.
Definition: ps.F90:188
The states_elec_t class contains all electronic wave functions.
int true(void)