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