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