39 use,
intrinsic :: iso_fortran_env
90 integer,
parameter :: &
91 INITRHO_PARAMAGNETIC = 1, &
99 logical :: complex_ylms
100 logical :: initialized
103 integer,
allocatable :: atom(:)
104 integer,
allocatable :: level(:)
105 integer,
allocatable :: ddim(:)
106 logical :: alternative
107 logical :: derivative
108 integer,
allocatable :: cst(:, :)
109 integer,
allocatable :: ck(:, :)
110 real(4),
allocatable :: dbuff_single(:, :, :, :)
111 complex(4),
allocatable :: zbuff_single(:, :, :, :)
112 real(real64),
allocatable :: dbuff(:, :, :, :)
113 complex(real64),
allocatable :: zbuff(:, :, :, :)
114 logical :: save_memory
115 logical :: initialized_orbitals
116 real(real64) :: orbital_scale_factor
119 logical,
allocatable :: is_empty(:)
125 real(real64),
allocatable :: radius(:)
126 real(real64) :: lapdist
130 integer,
allocatable :: basis_atom(:)
131 integer,
allocatable :: basis_orb(:)
132 integer,
allocatable :: atom_orb_basis(:, :)
133 integer,
allocatable :: norb_atom(:)
135 integer :: lsize(1:2)
136 integer :: nproc(1:2)
137 integer :: myroc(1:2)
138 integer :: desc(1:BLACS_DLEN)
139 logical,
allocatable :: calc_atom(:)
140 real(real64) :: diag_tol
141 type(submesh_t),
allocatable :: sphere(:)
142 type(batch_t),
allocatable :: orbitals(:)
143 logical,
allocatable :: is_orbital_initialized(:)
150 real(real64),
parameter,
public :: default_eigenval = 1.0e10_real64
155 subroutine lcao_init(this, namespace, space, gr, ions, st, st_start)
156 type(lcao_t),
intent(out) :: this
157 type(namespace_t),
intent(in) :: namespace
158 type(electron_space_t),
intent(in) :: space
159 type(grid_t),
intent(in) :: gr
160 type(ions_t),
intent(in) :: ions
161 type(states_elec_t),
intent(in) :: st
162 integer,
intent(in) :: st_start
164 integer :: ia, n, iorb, jj, maxj, idim
165 integer :: ii, ll, mm, norbs, ii_tmp
166 integer :: mode_default
167 real(real64) :: max_orb_radius, maxradius
172 this%initialized = .
true.
177 mode_default = option__lcaostart__lcao_states
178 if (ions%only_user_def) mode_default = option__lcaostart__lcao_none
235 if (this%mode == option__lcaostart__lcao_states_batch .and. st_start > 1)
then
236 message(1) =
"LCAOStart = lcao_states_batch not compatible with this run."
237 message(2) =
"Please use LCAOStart = lcao_states instead"
242 this%alternative = this%mode == option__lcaostart__lcao_states_batch
244 if (this%mode == option__lcaostart__lcao_none)
then
254 if (st%d%ispin ==
spinors .and. this%alternative)
then
255 message(1) =
"LCAOStart = lcao_states_batch is not working for spinors."
258 if (space%is_periodic() .and. this%alternative)
then
276 call parse_variable(namespace,
'LCAOComplexYlms', .false., this%complex_ylms)
278 this%complex_ylms = .false.
289 call parse_variable(namespace,
'LCAOSaveMemory', .false., this%save_memory)
293 call io_mkdir(
'debug/lcao', namespace)
294 iunit_o =
io_open(
'debug/lcao/orbitals', namespace, action=
'write')
295 write(iunit_o,
'(7a6)')
'iorb',
'atom',
'level',
'i',
'l',
'm',
'spin'
298 if (.not. this%alternative)
then
319 call parse_variable(namespace,
'LCAOMaximumOrbitalRadius', 20.0_real64, max_orb_radius, unit =
units_inp%length)
325 do ia = 1, ions%natoms
326 maxj = max(maxj, ions%atom(ia)%species%get_niwfs())
327 this%maxorbs = this%maxorbs + ions%atom(ia)%species%get_niwfs()
330 this%maxorbs = this%maxorbs*st%d%dim
332 if (this%maxorbs == 0)
then
333 call messages_write(
'The are no atomic orbitals available, cannot do LCAO.')
335 this%mode = option__lcaostart__lcao_none
342 safe_allocate( this%atom(1:this%maxorbs))
343 safe_allocate(this%level(1:this%maxorbs))
344 safe_allocate( this%ddim(1:this%maxorbs))
346 safe_allocate(this%is_empty(1:this%maxorbs))
347 this%is_empty = .false.
350 this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
353 safe_allocate(this%radius(1:ions%natoms))
355 do ia = 1, ions%natoms
356 norbs = ions%atom(ia)%species%get_niwfs()
360 call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
362 if(ions%atom(ia)%species%is_full())
call ions%atom(ia)%species%get_iwf_n( iorb, 1, ii)
363 maxradius = max(maxradius, ions%atom(ia)%species%get_iwf_radius( ii, is = 1))
366 maxradius = min(maxradius,
m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
368 this%radius(ia) = maxradius
400 do ia = 1, ions%natoms
402 if (jj > ions%atom(ia)%species%get_niwfs()) cycle
403 call ions%atom(ia)%species%get_iwf_ilm(jj, idim, ii, ll, mm)
405 if(ions%atom(ia)%species%is_full())
then
407 call ions%atom(ia)%species%get_iwf_n( ii_tmp, 1, ii)
409 if (this%orbital_scale_factor*ions%atom(ia)%species%get_iwf_radius( ii, is = 1) >= max_orb_radius) cycle
412 this%level(iorb) = jj
413 this%ddim(iorb) = idim
416 write(iunit_o,
'(7i6)') iorb, this%atom(iorb), this%level(iorb), ii, ll, mm, this%ddim(iorb)
429 if (this%maxorbs /= iorb - 1)
then
434 call messages_write(
' orbitals cannot be used for the LCAO calculation,')
436 call messages_write(
' their radii exceeds LCAOMaximumOrbitalRadius (')
441 this%maxorbs = iorb - 1
444 if (this%maxorbs < st%nst)
then
445 call messages_write(
'Cannot do LCAO for all states because there are not enough atomic orbitals.')
485 if (n > 0 .and. n <= st%nst .and. st%nst <= this%maxorbs)
then
488 else if (n > st%nst .and. n <= this%maxorbs)
then
491 else if (n == 0)
then
493 this%norbs = min(this%maxorbs, 2*st%nst)
496 this%norbs = this%maxorbs
499 assert(this%norbs <= this%maxorbs)
501 safe_allocate(this%cst(1:this%norbs, 1:st%d%spin_channels))
502 safe_allocate(this%ck(1:this%norbs, 1:st%d%spin_channels))
503 this%initialized_orbitals = .false.
513 integer :: iatom, iorb, norbs
514 real(real64) :: maxradius
517 integer :: jatom, jorb, jbasis, ilbasis, jlbasis, proc(1:2),
info, nbl
558 call parse_variable(namespace,
'LCAOExtraOrbitals', .false., this%derivative)
565 if (this%derivative)
then
568 if (st%nst * st%smear%el_per_state > st%qtot)
then
569 message(1) =
"Lower-lying empty states may be missed with LCAOExtraOrbitals."
582 call parse_variable(namespace,
'LCAODiagTol', 1e-10_real64, this%diag_tol)
584 if (this%derivative)
then
590 safe_allocate(this%sphere(1:ions%natoms))
591 safe_allocate(this%orbitals(1:ions%natoms))
592 safe_allocate(this%is_orbital_initialized(1:ions%natoms))
593 this%is_orbital_initialized = .false.
595 safe_allocate(this%norb_atom(1:ions%natoms))
599 do iatom = 1, ions%natoms
600 this%norb_atom(iatom) = this%mult*ions%atom(iatom)%species%get_niwfs()
601 this%maxorb = max(this%maxorb, ions%atom(iatom)%species%get_niwfs())
602 this%norbs = this%norbs + ions%atom(iatom)%species%get_niwfs()
605 this%maxorb = this%maxorb*this%mult
606 this%norbs = this%norbs*this%mult
608 safe_allocate(this%basis_atom(1:this%norbs))
609 safe_allocate(this%basis_orb(1:this%norbs))
610 safe_allocate(this%atom_orb_basis(1:ions%natoms, 1:this%maxorb))
615 do iatom = 1, ions%natoms
616 norbs = ions%atom(iatom)%species%get_niwfs()
618 do iorb = 1, this%mult*norbs
620 this%atom_orb_basis(iatom, iorb) = ibasis
621 this%basis_atom(ibasis) = iatom
622 this%basis_orb(ibasis) = iorb
626 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
627 write(iunit_o,
'(7i6)') ibasis, iatom, iorb, ii, ll, mm, 1
637 this%lapdist = maxval(abs(gr%idx%enlarge)*gr%spacing)
640 safe_allocate(this%radius(1:ions%natoms))
642 do iatom = 1, ions%natoms
643 norbs = ions%atom(iatom)%species%get_niwfs()
647 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
649 if(ions%atom(iatom)%species%is_full())
call ions%atom(iatom)%species%get_iwf_n( iorb, 1, ii)
650 maxradius = max(maxradius, ions%atom(iatom)%species%get_iwf_radius( ii, is = 1))
653 if (this%derivative) maxradius = maxradius + this%lapdist
655 maxradius = min(maxradius,
m_two*maxval(gr%box%bounding_box_l(1:space%dim)))
657 this%radius(iatom) = maxradius
660 safe_allocate(this%calc_atom(1:ions%natoms))
661 this%calc_atom = .
true.
664#ifndef HAVE_SCALAPACK
665 this%parallel = .false.
667 this%parallel = (st%parallel_in_states .or. gr%parallel_in_domains) &
670 if (this%parallel)
then
671 nbl = min(16, this%norbs)
674 this%lsize(1) = max(1,
numroc(this%norbs, nbl, st%dom_st_proc_grid%myrow, 0, st%dom_st_proc_grid%nprow))
675 this%lsize(2) = max(1,
numroc(this%norbs, nbl, st%dom_st_proc_grid%mycol, 0, st%dom_st_proc_grid%npcol))
677 this%nproc(1) = st%dom_st_proc_grid%nprow
678 this%nproc(2) = st%dom_st_proc_grid%npcol
679 this%myroc(1) = st%dom_st_proc_grid%myrow
680 this%myroc(2) = st%dom_st_proc_grid%mycol
682 call descinit(this%desc(1), this%norbs, this%norbs, nbl, nbl, 0, 0, &
683 st%dom_st_proc_grid%context, this%lsize(1),
info)
686 write(
message(1),
'(a,i6)')
'descinit for BLACS failed with error code ',
info
690 this%calc_atom = .false.
691 do iatom = 1, ions%natoms
692 ibasis = this%atom_orb_basis(iatom, 1)
694 do jatom = 1, ions%natoms
695 jbasis = this%atom_orb_basis(jatom, 1)
697 do iorb = 1, this%norb_atom(iatom)
698 do jorb = 1, this%norb_atom(jatom)
700 ilbasis, jlbasis, proc(1), proc(2))
702 this%calc_atom(this%basis_atom(jbasis)) = &
703 this%calc_atom(this%basis_atom(jbasis)) .or. proc(2) == this%myroc(2)
721 subroutine lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r)
724 type(
grid_t),
intent(in) :: gr
725 type(
ions_t),
intent(in) :: ions
728 type(
v_ks_t),
intent(inout) :: ks
730 integer,
optional,
intent(in) :: st_start
731 real(real64),
optional,
intent(in) :: lmm_r
734 integer :: st_start_random, required_min_nst
736 logical :: is_orbital_dependent
740 if (
present(st_start))
then
743 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
744 calc_eigenval=.not.
present(st_start), calc_current=.false.)
746 assert(st_start >= 1)
747 if (st_start > st%nst)
then
767 if (.not.
present(st_start))
then
768 call lcao_guess_density(lcao, namespace, st, gr, hm, ions, st%qtot, st%d%ispin, st%rho)
771 assert(
present(lmm_r))
777 message(1) =
'Info: Setting up Hamiltonian.'
781 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.false., &
782 calc_current=.false., calc_energy=.false., force_semilocal=is_orbital_dependent)
785 st%eigenval = default_eigenval
794 if (
present(st_start))
then
795 write(
message(1),
'(a,i8,a)')
'Performing LCAO for states ', st_start,
' and above'
799 call lcao_wf(lcao, st, gr, ions, hm, namespace, start = st_start)
805 select case (st%d%ispin)
807 required_min_nst = int(st%qtot/2)
809 required_min_nst = int(st%qtot/2)
811 required_min_nst = int(st%qtot)
814 if (lcao%norbs <= required_min_nst .and. lcao%norbs < st%nst)
then
815 st%eigenval(lcao%norbs+1:,:) =
m_zero
819 if (.not.
present(st_start))
then
824 if (lcao%mode == option__lcaostart__lcao_full)
then
825 call v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, &
826 calc_eigenval = .false., calc_current=.false.)
828 assert(
present(lmm_r))
835 if (.not. lcao_done .or. lcao%norbs < st%nst)
then
838 st_start_random = lcao%norbs + 1
842 if (
present(st_start)) st_start_random = max(st_start, st_start_random)
844 if (st_start_random > 1)
then
845 write(
message(1),
'(a,i8,a)')
'Generating random wavefunctions for states ', st_start_random,
' and above'
854 if (.not. st%fixed_spins)
then
862 if (.not. lcao_done)
then
865 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
866 calc_eigenval=.not.
present(st_start), calc_current=.false., force_semilocal=is_orbital_dependent)
867 if (.not.
present(st_start))
then
873 else if (
present(st_start))
then
875 if (st_start > 1)
then
892 type(
lcao_t),
intent(inout) :: this
896 safe_deallocate_a(this%calc_atom)
897 safe_deallocate_a(this%norb_atom)
898 safe_deallocate_a(this%basis_atom)
899 safe_deallocate_a(this%basis_orb)
900 safe_deallocate_a(this%atom_orb_basis)
901 safe_deallocate_a(this%radius)
902 safe_deallocate_a(this%sphere)
903 safe_deallocate_a(this%orbitals)
905 safe_deallocate_a(this%atom)
906 safe_deallocate_a(this%level)
907 safe_deallocate_a(this%ddim)
908 safe_deallocate_a(this%cst)
909 safe_deallocate_a(this%ck)
910 safe_deallocate_a(this%dbuff_single)
911 safe_deallocate_a(this%zbuff_single)
912 safe_deallocate_a(this%dbuff)
913 safe_deallocate_a(this%zbuff)
915 this%initialized = .false.
921 subroutine lcao_wf(this, st, gr, ions, hm, namespace, start)
922 type(
lcao_t),
intent(inout) :: this
924 type(
grid_t),
intent(in) :: gr
925 type(
ions_t),
intent(in) :: ions
928 integer,
optional,
intent(in) :: start
932 assert(this%initialized)
938 if (
present(start)) start_ = start
940 if (this%alternative)
then
942 call dlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
944 call zlcao_alt_wf(this, st, gr, ions, hm, namespace, start_)
948 call dlcao_wf(this, st, gr, ions, hm, namespace, start_)
950 call zlcao_wf(this, st, gr, ions, hm, namespace, start_)
961 type(
lcao_t),
intent(in) :: this
965 available = this%initialized .and. this%mode /= option__lcaostart__lcao_none &
975 type(
lcao_t),
intent(in) :: this
986 type(
lcao_t),
intent(in) :: this
987 integer,
intent(in) :: ig
988 integer,
intent(in) :: jg
989 integer,
intent(out) :: il
990 integer,
intent(out) :: jl
991 integer,
intent(out) :: prow
992 integer,
intent(out) :: pcol
996 call infog2l(ig, jg, this%desc(1), this%nproc(1), this%nproc(2), this%myroc(1), this%myroc(2), &
1013 type(
lcao_t),
intent(inout) :: this
1014 integer,
intent(in) :: iatom
1018 if (this%is_orbital_initialized(iatom))
then
1019 call this%orbitals(iatom)%end()
1020 this%is_orbital_initialized(iatom) = .false.
1030 type(
lcao_t),
intent(inout) :: this
1032 class(
mesh_t),
intent(in) :: mesh
1034 integer,
intent(in) :: iatom
1035 integer,
intent(in) :: spin_channels
1036 real(real64),
intent(inout) :: rho(:, :)
1038 real(real64),
allocatable :: dorbital(:, :)
1039 complex(real64),
allocatable :: zorbital(:, :)
1040 real(real64),
allocatable :: factors(:)
1041 real(real64) :: factor, aa
1042 integer :: iorb, ip, ii, ll, mm, ispin
1043 type(
ps_t),
pointer :: ps
1044 logical :: use_stored_orbitals
1050 use_stored_orbitals = ions%atom(iatom)%species%is_ps() &
1052 .and. st%d%dim == 1 .and. .not. ions%space%is_periodic()
1056 if (use_stored_orbitals)
then
1058 assert(.not. ions%space%is_periodic())
1060 select type(spec=>ions%atom(iatom)%species)
1067 if (.not. this%alternative)
then
1070 safe_allocate(dorbital(1:mesh%np, 1:st%d%dim))
1072 safe_allocate(zorbital(1:mesh%np, 1:st%d%dim))
1075 do iorb = 1, this%norbs
1076 if (iatom /= this%atom(iorb)) cycle
1078 call ions%atom(iatom)%species%get_iwf_ilm(this%level(iorb), 1, ii, ll, mm)
1079 factor = ps%conf%occ(ii, 1)/(
m_two*ll +
m_one)
1082 call dget_ao(this, st, mesh, ions, iorb, 1, dorbital, use_psi = .
true.)
1085 rho(ip, 1) = rho(ip, 1) + factor*dorbital(ip, 1)**2
1088 call zget_ao(this, st, mesh, ions, iorb, 1, zorbital, use_psi = .
true.)
1091 rho(ip, 1) = rho(ip, 1) + factor*abs(zorbital(ip, 1))**2
1097 safe_deallocate_a(dorbital)
1098 safe_deallocate_a(zorbital)
1106 safe_allocate(factors(1:this%norb_atom(iatom)/this%mult))
1108 do iorb = 1, this%norb_atom(iatom)/this%mult
1109 call ions%atom(iatom)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
1110 factors(iorb) = ps%conf%occ(ii, 1)/(
m_two*ll +
m_one)
1114 do ip = 1, this%sphere(iatom)%np
1116 do iorb = 1, this%norb_atom(iatom)/this%mult
1117 aa = aa + factors(iorb)*this%orbitals(iatom)%dff_linear(ip, iorb)**2
1119 rho(this%sphere(iatom)%map(ip), 1) = rho(this%sphere(iatom)%map(ip), 1) + aa
1122 safe_deallocate_a(factors)
1128 ions%pos(:, iatom), mesh, spin_channels, rho)
1133 do ispin = 1, spin_channels
1136 rho(ip, ispin) = max(rho(ip, ispin),
m_zero)
1146 type(
lcao_t),
intent(inout) :: this
1149 type(
grid_t),
intent(in) :: gr
1151 type(
ions_t),
intent(in) :: ions
1152 real(real64),
intent(in) :: qtot
1153 integer,
intent(in) :: ispin
1154 real(real64),
contiguous,
intent(out) :: rho(:, :)
1156 integer :: ia, is, idir, ip, m_dim
1159 real(real64) :: rr, rnd, phi, theta, lmag, n1, n2
1160 real(real64),
allocatable :: atom_rho(:,:), mag(:,:)
1161 real(real64),
parameter :: tol_min_mag = 1.0e-20_real64
1165 if (st%d%spin_channels == 1)
then
1166 this%gmd_opt = initrho_paramagnetic
1199 message(1) =
"GuessMagnetDensity cannot be used for Hartree-Fock and generalized Kohn-Sham calculation."
1200 message(2) =
"Please perform a LDA or GGA calculation first and restart from this calculation."
1224 if (
parse_block(namespace,
'AtomsMagnetDirection', blk) < 0)
then
1225 message(1) =
"AtomsMagnetDirection block is not defined."
1230 message(1) =
"The number of rows in the AtomsMagnetDirection block does not equal the number of atoms."
1240 safe_allocate(mag(1:m_dim, 1:ions%natoms))
1241 do ia = 1, ions%natoms
1245 if (abs(mag(idir, ia)) < tol_min_mag) mag(idir, ia) =
m_zero
1253 safe_allocate(atom_rho(1:gr%np, 1:st%d%spin_channels))
1254 select case (this%gmd_opt)
1255 case (initrho_paramagnetic)
1256 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1261 if (st%d%spin_channels == 2)
then
1264 rho(ip, 1) =
m_half*(rho(ip, 1) + rho(ip, 2))
1265 rho(ip, 2) = rho(ip, 1)
1270 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1272 rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
1276 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1291 call lalg_axpy(gr%np, (lmag - n1 + n2)/
m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
1300 rho(1:gr%np, 1:2) = rho(1:gr%np, 1:2) + atom_rho(1:gr%np, 1:2)
1302 rho(1:gr%np, 1) = rho(1:gr%np, 1) + atom_rho(1:gr%np, 2)
1303 rho(1:gr%np, 2) = rho(1:gr%np, 2) + atom_rho(1:gr%np, 1)
1305 elseif (ispin ==
spinors)
then
1316 do ia = ions%atoms_dist%start, ions%atoms_dist%end
1324 lmag = norm2(mag(:, ia))
1325 if (lmag > n1 + n2)
then
1326 mag = mag*(n1 + n2)/lmag
1334 atom_rho(ip, 1) =
m_half*(atom_rho(ip, 1) + atom_rho(ip, 2))
1335 rho(ip, 1) = rho(ip, 1) + atom_rho(ip, 1)
1336 rho(ip, 2) = rho(ip, 2) + atom_rho(ip, 1)
1343 if (n1 - n2 < lmag)
then
1344 call lalg_axpy(gr%np, (lmag - n1 + n2)/
m_two/n2, atom_rho(:, 2), atom_rho(:, 1))
1346 elseif (n1 - n2 > lmag)
then
1347 call lalg_axpy(gr%np, (n1 - n2 - lmag)/
m_two/n1, atom_rho(:, 1), atom_rho(:, 2))
1355 if (mag(1, ia) >
m_zero)
then
1362 elseif (ispin ==
spinors)
then
1372 if (ions%atoms_dist%parallel)
then
1373 do is = 1, st%d%nspin
1374 call lalg_copy(gr%np, rho(:,is), atom_rho(:,1))
1375 call ions%atoms_dist%mpi_grp%allreduce(atom_rho(1, 1), rho(1, is), gr%np, mpi_double_precision, mpi_sum)
1381 write(
message(1),
'(a,f13.6)')
'Info: Unnormalized total charge = ', rr
1386 call lalg_scal(gr%np, st%d%nspin, qtot/rr, rho)
1388 write(
message(1),
'(a,f13.6)')
'Info: Renormalized total charge = ', rr
1393 if (st%symmetrize_density)
then
1394 do is = 1, st%d%nspin
1399 safe_deallocate_a(atom_rho)
1400 safe_deallocate_a(mag)
1406 real(real64) function integrated_charge_density(gr, st, rho) result(rr)
1407 type(
grid_t),
intent(in) :: gr
1409 real(real64),
intent(in) :: rho(:,:)
1414 do is = 1, st%d%spin_channels
1418 call gr%allreduce(rr)
1423 class(mesh_t),
intent(in) :: mesh
1424 real(real64),
intent(inout) :: rho(:,:)
1425 real(real64),
intent(in) :: atom_rho(:,:)
1426 real(real64),
intent(in) :: theta, phi
1429 real(real64) :: half_theta
1433 half_theta = m_half * theta
1437 rho(ip, 1) = rho(ip, 1) +
cos(half_theta)**2*atom_rho(ip, 1) +
sin(half_theta)**2*atom_rho(ip, 2)
1438 rho(ip, 2) = rho(ip, 2) +
sin(half_theta)**2*atom_rho(ip, 1) +
cos(half_theta)**2*atom_rho(ip, 2)
1439 rho(ip, 3) = rho(ip, 3) +
cos(half_theta)*
sin(half_theta)*
cos(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
1440 rho(ip, 4) = rho(ip, 4) -
cos(half_theta)*
sin(half_theta)*
sin(phi)*(atom_rho(ip, 1)-atom_rho(ip, 2))
1449 type(
lcao_t),
intent(inout) :: this
1450 type(namespace_t),
intent(in) :: namespace
1451 type(states_elec_t),
intent(inout) :: st
1452 type(grid_t),
intent(in) :: gr
1453 type(ions_t),
intent(in) :: ions
1454 integer,
optional,
intent(in) :: start
1460 if (.not. this%alternative)
then
1461 if (states_are_real(st))
then
1467 if (states_are_real(st))
then
1480 real(real64),
intent(in) :: mag(3)
1481 real(real64),
intent(in) :: lmag
1482 real(real64),
intent(out) :: theta
1483 real(real64),
intent(out) :: phi
1485 assert(lmag > m_zero)
1487 theta =
acos(max(-1.0_real64, min(1.0_real64, mag(3)/lmag)))
1489 if (abs(mag(1)) <= m_epsilon)
then
1490 if (abs(mag(2)) <= m_epsilon)
then
1492 elseif (mag(2) < m_zero)
then
1493 phi = m_pi*m_three*m_half
1499 phi =
atan2(mag(2), mag(1))
1514 type(states_elec_t),
intent(inout) :: st
1515 type(grid_t),
intent(in) :: gr
1516 integer,
intent(in) :: ist_start
1517 integer,
intent(in) :: gmd_opt
1519 integer :: ik, ist, ip
1520 real(real64),
allocatable :: md(:,:), up(:)
1521 complex(real64),
allocatable :: down(:), psi(:,:)
1522 real(real64) :: lmag, theta, phi, norm, mm(3)
1524 if (st%d%ispin /= spinors)
return
1529 safe_allocate(psi(1:gr%np, 1:st%d%dim))
1531 select case (gmd_opt)
1534 do ik = st%d%kpt%start, st%d%kpt%end
1535 do ist = ist_start, st%st_end
1536 call states_elec_get_state(st, gr, ist, ik, psi)
1539 psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
1540 psi(ip, 2) = psi(ip, 1)
1542 call zmf_normalize(gr, st%d%dim, psi)
1543 call states_elec_set_state(st, gr, ist, ik, psi)
1549 do ik = st%d%kpt%start, st%d%kpt%end
1550 do ist = ist_start, st%st_end
1551 call states_elec_get_state(st, gr, ist, ik, psi)
1554 psi(ip, 1) = psi(ip, 1) + psi(ip, 2)
1557 call zmf_normalize(gr, st%d%dim, psi)
1558 call states_elec_set_state(st, gr, ist, ik, psi)
1564 safe_allocate(md(1:gr%np, 1:3))
1565 safe_allocate(up(1:gr%np))
1566 safe_allocate(down(1:gr%np))
1567 call magnetic_density(gr, st%d, st%rho, md)
1569 mm(1) = dmf_integrate(gr, md(:, 1), reduce = .false.)
1570 mm(2) = dmf_integrate(gr, md(:, 2), reduce = .false.)
1571 mm(3) = dmf_integrate(gr, md(:, 3), reduce = .false.)
1573 if (gr%parallel_in_domains)
then
1574 call gr%allreduce(mm)
1580 if (lmag > 1.0e-2_real64)
then
1582 up(:) =
cos(theta/m_two)
1583 down(:) =
exp(-m_zi * phi) *
sin(theta/m_two)
1590 lmag = norm2(md(ip, 1:3))
1592 up(ip) =
cos(theta/m_two)
1593 down(ip) =
exp(-m_zi * phi) *
sin(theta/m_two)
1595 safe_deallocate_a(md)
1599 do ik = st%d%kpt%start, st%d%kpt%end
1600 do ist = ist_start, st%st_end
1601 call states_elec_get_state(st, gr, ist, ik, psi)
1602 call zmf_normalize(gr, st%d%dim, psi)
1605 norm =
sqrt(real(psi(ip, 1)*conjg(psi(ip, 1)) + psi(ip, 2)*conjg(psi(ip, 2))))
1606 psi(ip, 1) = norm * up(ip)
1607 psi(ip, 2) = norm * down(ip)
1609 call states_elec_set_state(st, gr, ist, ik, psi)
1612 safe_deallocate_a(up)
1613 safe_deallocate_a(down)
1616 safe_deallocate_a(psi)
1623#include "lcao_inc.F90"
1626#include "complex.F90"
1627#include "lcao_inc.F90"
constant times a vector plus a vector
Copies a vector x, to a vector y.
scales a vector by a constant
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
double acos(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
This module implements batches of mesh functions.
This module contains interfaces for BLACS routines Interfaces are from http:
This module provides the BLACS processor grid.
logical pure function, public blacs_proc_grid_null(this)
Module implementing boundary conditions in Octopus.
type(debug_t), save, public debug
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_pi
some mathematical constants
integer, parameter, public hartree_fock
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
integer, parameter, public hartree
This module implements the underlying real-space grid.
subroutine, public dgrid_symmetrize_scalar_field(gr, field, suppress_warning)
This module defines classes and functions for interaction partners.
subroutine, public io_close(iunit, grp)
subroutine, public io_mkdir(fname, namespace, parents)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
A module to handle KS potential, without the external potential.
This module contains interfaces for LAPACK routines.
subroutine, public lcao_init_orbitals(this, namespace, st, gr, ions, start)
subroutine, public lcao_wf(this, st, gr, ions, hm, namespace, start)
subroutine zlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
subroutine zget_ao(this, st, mesh, ions, iorb, ispin, ao, use_psi)
subroutine, public lcao_end(this)
subroutine, public lcao_run(namespace, space, gr, ions, ext_partners, st, ks, hm, st_start, lmm_r)
subroutine zlcao_alt_wf(this, st, gr, ions, hm, namespace, start)
The alternative implementation.
subroutine dlcao_wf(this, st, gr, ions, hm, namespace, start)
subroutine dlcao_alt_init_orbitals(this, namespace, st, gr, ions, start)
subroutine zlcao_wf(this, st, gr, ions, hm, namespace, start)
subroutine rotate_random_states_to_local_frame(st, gr, ist_start, gmd_opt)
Rotate the spinors with band index >= ist_start to the local frame of the magnetization.
subroutine dlcao_alt_wf(this, st, gr, ions, hm, namespace, start)
The alternative implementation.
integer, parameter initrho_userdef
subroutine get_angles_from_magnetization(mag, lmag, theta, phi)
Given a magnetization vector, and its norm, this returns the corresponding inclination and azimuthal ...
integer, parameter initrho_random
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...
subroutine lcao_local_index(this, ig, jg, il, jl, prow, pcol)
integer function, public lcao_num_orbitals(this)
Returns the number of LCAO orbitas.
subroutine lcao_guess_density(this, namespace, st, gr, hm, ions, qtot, ispin, rho)
builds a density which is the sum of the atomic densities
subroutine dinit_orbitals(this, namespace, st, gr, ions, start)
real(real64) function integrated_charge_density(gr, st, rho)
Computes the integral of rho, summed over spin channels.
integer, parameter initrho_paramagnetic
subroutine dget_ao(this, st, mesh, ions, iorb, ispin, ao, use_psi)
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...
subroutine zinit_orbitals(this, namespace, st, gr, ions, start)
subroutine lcao_atom_density(this, st, mesh, ions, iatom, spin_channels, rho)
subroutine, public lcao_init(this, namespace, space, gr, ions, st, st_start)
subroutine accumulate_rotated_density(mesh, rho, atom_rho, theta, phi)
integer, parameter initrho_ferromagnetic
logical function, public lcao_is_available(this)
Returns true if LCAO can be done.
subroutine, public write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
This module is intended to contain "only mathematical" functions and procedures.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
subroutine, public messages_new_line()
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
type(mpi_grp_t), public mpi_world
logical function, public parse_is_defined(namespace, name)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
integer(int64), parameter, public splitmix64_321
This module contains interfaces for ScaLAPACK routines Interfaces are from http:
integer, parameter, public smear_semiconductor
integer, parameter, public smear_fixed_occ
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.
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)
subroutine, public v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
logical pure function, public xc_is_orbital_dependent(xcs)
Is the xc family orbital dependent.
integer, parameter, public sic_pz_oep
Perdew-Zunger SIC (OEP way)
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Describes mesh distribution to nodes.
A type storing the information and data about a pseudopotential.
The states_elec_t class contains all electronic wave functions.