37 use,
intrinsic :: iso_fortran_env
109 real(real64) :: total(3,3) =
m_zero
110 real(real64) :: kinetic(3,3) =
m_zero
111 real(real64) :: Hartree(3,3) =
m_zero
112 real(real64) :: xc(3,3) =
m_zero
114 real(real64) :: ps_nl(3,3) =
m_zero
115 real(real64) :: ion_ion(3,3) =
m_zero
116 real(real64) :: vdw(3,3) =
m_zero
117 real(real64) :: hubbard(3,3) =
m_zero
119 real(real64) :: kinetic_sumrule
120 real(real64) :: hartree_sumrule
135 type(states_elec_dim_t) :: d
138 logical :: only_userdef_istates
140 type(states_elec_group_t) :: group
141 integer :: block_size
145 logical :: pack_states
151 character(len=1024),
allocatable :: user_def_states(:,:,:)
157 real(real64),
allocatable :: rho(:,:)
158 real(real64),
allocatable :: rho_core(:)
160 real(real64),
allocatable :: current(:, :, :)
161 real(real64),
allocatable :: current_para(:, :, :)
162 real(real64),
allocatable :: current_dia(:, :, :)
163 real(real64),
allocatable :: current_mag(:, :, :)
164 real(real64),
allocatable :: current_kpt(:,:,:)
169 real(real64),
allocatable :: frozen_rho(:, :)
170 real(real64),
allocatable :: frozen_tau(:, :)
171 real(real64),
allocatable :: frozen_gdens(:,:,:)
172 real(real64),
allocatable :: frozen_ldens(:,:)
174 logical :: uniform_occ
176 real(real64),
allocatable :: eigenval(:,:)
178 logical :: restart_fixed_occ
179 logical :: restart_reorder_occs
180 real(real64),
allocatable :: occ(:,:)
181 real(real64),
allocatable :: kweights(:)
184 logical,
private :: fixed_spins
186 real(real64),
allocatable :: spin(:, :, :)
189 real(real64) :: val_charge
191 type(stress_t) :: stress_tensors
193 logical :: fromScratch
194 type(smear_t) :: smear
197 type(modelmb_particle_t) :: modelmbparticles
198 integer,
allocatable :: mmb_nspindown(:,:)
199 integer,
allocatable :: mmb_iyoung(:,:)
200 real(real64),
allocatable :: mmb_proj(:)
202 logical :: parallel_in_states = .false.
213 logical :: scalapack_compatible
217 integer :: st_start, st_end
218 integer,
allocatable :: node(:)
221 integer,
allocatable :: st_kpt_task(:,:)
224 logical :: symmetrize_density
225 integer :: randomization
226 integer :: orth_method = 0
228 real(real64) :: cl_states_mem
240 integer,
public,
parameter :: &
284 real(real64),
intent(in) :: valence_charge
287 real(real64) :: excess_charge, nempty_percent
288 integer :: nempty, ntot, default
289 integer :: nempty_conv
291 real(real64),
parameter :: tol = 1e-13_real64
300 st%d%ispin = space%ispin
305 if (st%d%ispin /=
unpolarized .and. kpoints%use_time_reversal)
then
306 message(1) =
"Time reversal symmetry is only implemented for unpolarized spins."
307 message(2) =
"Use KPointsUseTimeReversal = no."
339 write(
message(1),
'(a,i5,a)')
"Input: '", ntot,
"' is not a valid value for TotalStates."
361 write(
message(1),
'(a,i5,a)')
"Input: '", nempty,
"' is not a valid value for ExtraStates."
362 message(2) =
'(0 <= ExtraStates)'
366 if (ntot > 0 .and. nempty > 0)
then
367 message(1) =
'You cannot set TotalStates and ExtraStates at the same time.'
382 if (nempty_percent < 0)
then
383 write(
message(1),
'(a,f8.6,a)')
"Input: '", nempty_percent, &
384 "' should be a percentage value x (where x is parts in hundred) larger or equal 0"
388 if (nempty > 0 .and. nempty_percent > 0)
then
389 message(1) =
'You cannot set ExtraStates and ExtraStatesInPercent at the same time.'
408 call parse_variable(namespace,
'ExtraStatesToConverge', nempty, nempty_conv)
410 write(
message(1),
'(a,i5,a)')
"Input: '", nempty_conv,
"' is not a valid value for ExtraStatesToConverge."
411 message(2) =
'(0 <= ExtraStatesToConverge)'
415 if (nempty_conv > nempty)
then
416 message(1) =
'You cannot set ExtraStatesToConverge to a higher value than ExtraStates.'
423 st%val_charge = valence_charge
425 st%qtot = -(st%val_charge + excess_charge)
428 write(
message(1),
'(a,f12.6,a)')
'Total charge = ', st%qtot,
' < 0'
429 message(2) =
'Check Species and ExcessCharge.'
433 select case (st%d%ispin)
436 st%nst = nint(st%qtot/2)
437 if (st%nst*2 - st%qtot < -tol) st%nst = st%nst + 1
439 st%d%spin_channels = 1
442 st%nst = nint(st%qtot/2)
443 if (st%nst*2 - st%qtot < -tol) st%nst = st%nst + 1
445 st%d%spin_channels = 2
448 st%nst = nint(st%qtot)
449 if (st%nst - st%qtot < -tol) st%nst = st%nst + 1
451 st%d%spin_channels = 2
455 if (ntot < st%nst)
then
456 message(1) =
'TotalStates is smaller than the number of states required by the system.'
463 if (nempty_percent > 0)
then
464 nempty = ceiling(nempty_percent * st%nst / 100)
467 st%nst_conv = st%nst + nempty_conv
468 st%nst = st%nst + nempty
469 if (st%nst == 0)
then
470 message(1) =
"Cannot run with number of states = zero."
490 default = max(
accel%warp_size, 32)
499 call parse_variable(namespace,
'StatesBlockSize', default, st%block_size)
500 if (st%block_size < 1)
then
501 call messages_write(
"The variable 'StatesBlockSize' must be greater than 0.")
505 st%block_size = min(st%block_size, st%nst)
506 conf%target_states_block_size = st%block_size
508 safe_allocate(st%eigenval(1:st%nst, 1:st%nik))
509 st%eigenval = huge(st%eigenval)
513 if (.not. kpoints%gamma_only())
then
526 call parse_variable(namespace,
'OnlyUserDefinedInitialStates', .false., st%only_userdef_istates)
529 safe_allocate(st%occ (1:st%nst, 1:st%nik))
534 safe_allocate(st%user_def_states(1:st%d%dim, 1:st%nst, 1:st%nik))
536 st%user_def_states(1:st%d%dim, 1:st%nst, 1:st%nik) =
'undefined'
539 if (st%d%ispin ==
spinors)
then
540 safe_allocate(st%spin(1:3, 1:st%nst, 1:st%nik))
567 if (st%smear%photodop)
then
568 if (nempty == 0)
then
569 write(
message(1),
'(a,i5,a)')
"PhotoDoping requires to specify ExtraStates."
570 message(2) =
'(0 == ExtraStates)'
578 safe_allocate(st%node(1:st%nst))
579 st%node(1:st%nst) = 0
582 st%parallel_in_states = .false.
587 if (st%modelmbparticles%nparticle > 0)
then
589 safe_allocate(st%mmb_nspindown(1:st%modelmbparticles%ntype_of_particle, 1:st%nst))
590 st%mmb_nspindown(:,:) = -1
591 safe_allocate(st%mmb_iyoung(1:st%modelmbparticles%ntype_of_particle, 1:st%nst))
592 st%mmb_iyoung(:,:) = -1
593 safe_allocate(st%mmb_proj(1:st%nst))
605 call parse_variable(namespace,
'SymmetrizeDensity', kpoints%use_symmetries, st%symmetrize_density)
635 integer,
intent(out) :: nik
636 integer,
intent(out) :: dim
637 integer,
intent(out) :: nst
638 integer,
intent(out) :: ierr
640 character(len=256) :: lines(3)
641 character(len=20) :: char
651 read(lines(1), *) char, nst
652 read(lines(2), *) char, dim
653 read(lines(3), *) char, nik
673 real(real64),
intent(in) :: excess_charge
676 integer :: ik, ist, ispin, nspin, ncols, nrows, el_per_state, icol, start_pos, spin_n
678 real(real64) :: rr, charge
679 logical :: integral_occs, unoccupied_states
680 real(real64),
allocatable :: read_occs(:, :)
681 real(real64) :: charge_in_block
694 call parse_variable(namespace,
'RestartFixedOccupations', .
true., st%restart_fixed_occ)
765 integral_occs = .
true.
767 occ_fix:
if (
parse_block(namespace,
'Occupations', blk) == 0)
then
769 st%fixed_occ = .
true.
772 if (ncols > st%nst)
then
777 if (nrows /= st%nik)
then
778 call messages_input_error(namespace,
"Occupations",
"Wrong number of rows in block Occupations.")
781 do ik = 1, st%nik - 1
784 "All rows in block Occupations must have the same number of columns.")
795 safe_allocate(read_occs(1:ncols, 1:st%nik))
801 charge_in_block = charge_in_block + read_occs(icol, ik) * st%kweights(ik)
806 select case (st%d%ispin)
815 start_pos = nint((st%qtot - charge_in_block)/spin_n)
817 if (start_pos + ncols > st%nst)
then
818 message(1) =
"To balance charge, the first column in block Occupations is taken to refer to state"
819 write(
message(2),
'(a,i6,a)')
"number ", start_pos,
" but there are too many columns for the number of states."
820 write(
message(3),
'(a,i6,a)')
"Solution: set ExtraStates = ", start_pos + ncols - st%nst
825 do ist = 1, start_pos
826 st%occ(ist, ik) = el_per_state
831 do ist = start_pos + 1, start_pos + ncols
832 st%occ(ist, ik) = read_occs(ist - start_pos, ik)
833 integral_occs = integral_occs .and. &
834 abs((st%occ(ist, ik) - el_per_state) * st%occ(ist, ik)) <=
m_epsilon
839 do ist = start_pos + ncols + 1, st%nst
846 safe_deallocate_a(read_occs)
849 st%fixed_occ = .false.
850 integral_occs = .false.
857 st%qtot = -(st%val_charge + excess_charge)
860 if (st%d%nspin == 2) nspin = 2
862 do ik = 1, st%nik, nspin
865 do ispin = ik, ik + nspin - 1
866 st%occ(ist, ispin) = min(rr, -(st%val_charge + excess_charge) - charge)
867 charge = charge + st%occ(ist, ispin)
886 if (st%fixed_occ)
then
887 call parse_variable(namespace,
'RestartReorderOccs', .false., st%restart_reorder_occs)
889 st%restart_reorder_occs = .false.
892 call smear_init(st%smear, namespace, st%d%ispin, st%fixed_occ, integral_occs, kpoints)
894 unoccupied_states = (st%d%ispin /=
spinors .and. st%nst*2 > st%qtot) .or. (st%d%ispin ==
spinors .and. st%nst > st%qtot)
897 if (.not. unoccupied_states)
then
898 call messages_write(
'Smearing needs unoccupied states (via ExtraStates or TotalStates) to be useful.')
906 charge = charge + sum(st%occ(ist, 1:st%nik) * st%kweights(1:st%nik))
908 if (abs(charge - st%qtot) > 1e-6_real64)
then
909 message(1) =
"Initial occupations do not integrate to total charge."
910 write(
message(2),
'(6x,f12.6,a,f12.6)') charge,
' != ', st%qtot
937 st%fixed_spins = .false.
938 if (st%d%ispin /=
spinors)
then
977 spin_fix:
if (
parse_block(namespace,
'InitialSpins', blk) == 0)
then
985 st%fixed_spins = .
true.
987 st%spin(:, :, i) = st%spin(:, :, 1)
1000 class(
mesh_t),
intent(in) :: mesh
1001 type(
type_t),
optional,
intent(in) :: wfs_type
1002 logical,
optional,
intent(in) :: skip(:)
1003 logical,
optional,
intent(in) :: packed
1007 if (
present(wfs_type))
then
1009 st%wfs_type = wfs_type
1037 type(
mesh_t),
intent(in) :: mesh
1038 logical,
optional,
intent(in) :: verbose
1039 logical,
optional,
intent(in) :: skip(:)
1040 logical,
optional,
intent(in) :: packed
1042 integer :: ib, iqn, ist, istmin, istmax
1043 logical :: same_node, verbose_, packed_
1044 integer,
allocatable :: bstart(:), bend(:)
1048 safe_allocate(bstart(1:st%nst))
1049 safe_allocate(bend(1:st%nst))
1050 safe_allocate(st%group%iblock(1:st%nst))
1059 if (
present(skip))
then
1061 if (.not. skip(ist))
then
1069 if (
present(skip))
then
1070 do ist = st%nst, istmin, -1
1071 if (.not. skip(ist))
then
1078 if (
present(skip) .and. verbose_)
then
1088 st%group%nblocks = 0
1090 do ist = istmin, istmax
1093 st%group%iblock(ist) = st%group%nblocks + 1
1096 if (st%parallel_in_states .and. ist /= istmax)
then
1099 same_node = (st%node(ist + 1) == st%node(ist))
1102 if (ib == st%block_size .or. ist == istmax .or. .not. same_node)
then
1104 st%group%nblocks = st%group%nblocks + 1
1105 bend(st%group%nblocks) = ist
1106 if (ist /= istmax) bstart(st%group%nblocks + 1) = ist + 1
1110 safe_allocate(st%group%psib(1:st%group%nblocks, st%d%kpt%start:st%d%kpt%end))
1112 safe_allocate(st%group%block_is_local(1:st%group%nblocks, 1:st%nik))
1113 st%group%block_is_local = .false.
1114 st%group%block_start = -1
1115 st%group%block_end = -2
1117 do ib = 1, st%group%nblocks
1118 if (bstart(ib) >= st%st_start .and. bend(ib) <= st%st_end)
then
1119 if (st%group%block_start == -1) st%group%block_start = ib
1120 st%group%block_end = ib
1121 do iqn = st%d%kpt%start, st%d%kpt%end
1122 st%group%block_is_local(ib, iqn) = .
true.
1125 call dwfs_elec_init(st%group%psib(ib, iqn), st%d%dim, bstart(ib), bend(ib), mesh%np_part, iqn, &
1126 special=.
true., packed=packed_)
1128 call zwfs_elec_init(st%group%psib(ib, iqn), st%d%dim, bstart(ib), bend(ib), mesh%np_part, iqn, &
1129 special=.
true., packed=packed_)
1136 safe_allocate(st%group%block_range(1:st%group%nblocks, 1:2))
1137 safe_allocate(st%group%block_size(1:st%group%nblocks))
1139 st%group%block_range(1:st%group%nblocks, 1) = bstart(1:st%group%nblocks)
1140 st%group%block_range(1:st%group%nblocks, 2) = bend(1:st%group%nblocks)
1141 st%group%block_size(1:st%group%nblocks) = bend(1:st%group%nblocks) - bstart(1:st%group%nblocks) + 1
1143 st%group%block_initialized = .
true.
1145 safe_allocate(st%group%block_node(1:st%group%nblocks, 1:st%nik))
1146 st%group%block_node = 0
1148 assert(
allocated(st%node))
1149 assert(all(st%node >= 0) .and. all(st%node < st%mpi_grp%size))
1151 do iqn = st%d%kpt%start, st%d%kpt%end
1152 do ib = st%group%block_start, st%group%block_end
1153 st%group%block_node(ib, iqn) = st%st_kpt_mpi_grp%rank
1162 do ib = 1, st%group%nblocks
1168 if (st%group%block_size(ib) > 0)
then
1198 safe_deallocate_a(bstart)
1199 safe_deallocate_a(bend)
1220 type(
grid_t),
intent(in) :: gr
1222 real(real64) :: fsize
1226 safe_allocate(st%rho(1:gr%np_part, 1:st%d%nspin))
1229 fsize = gr%np_part*8.0_real64*st%block_size
1241 class(
space_t),
intent(in) :: space
1242 class(
mesh_t),
intent(in) :: mesh
1246 if (.not.
allocated(st%current))
then
1247 safe_allocate(st%current(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1251 if (.not.
allocated(st%current_para))
then
1252 safe_allocate(st%current_para(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1256 if (.not.
allocated(st%current_dia))
then
1257 safe_allocate(st%current_dia(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1261 if (.not.
allocated(st%current_mag))
then
1262 safe_allocate(st%current_mag(1:mesh%np_part, 1:space%dim, 1:st%d%nspin))
1266 if (.not.
allocated(st%current_kpt))
then
1267 safe_allocate(st%current_kpt(1:mesh%np,1:space%dim,st%d%kpt%start:st%d%kpt%end))
1346 default = option__statesorthogonalization__cholesky_serial
1347#ifdef HAVE_SCALAPACK
1349 default = option__statesorthogonalization__cholesky_parallel
1353 call parse_variable(namespace,
'StatesOrthogonalization', default, st%orth_method)
1374 call parse_variable(namespace,
'StatesCLDeviceMemory', -512.0_real64, st%cl_states_mem)
1383 subroutine states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
1386 logical,
optional,
intent(in) :: exclude_wfns
1387 logical,
optional,
intent(in) :: exclude_eigenval
1388 logical,
optional,
intent(in) :: special
1390 logical :: exclude_wfns_
1399 safe_allocate_source_a(stout%kweights, stin%kweights)
1400 stout%nik = stin%nik
1403 if (stin%modelmbparticles%nparticle > 0)
then
1404 safe_allocate_source_a(stout%mmb_nspindown, stin%mmb_nspindown)
1405 safe_allocate_source_a(stout%mmb_iyoung, stin%mmb_iyoung)
1406 safe_allocate_source_a(stout%mmb_proj, stin%mmb_proj)
1409 stout%wfs_type = stin%wfs_type
1410 stout%nst = stin%nst
1411 stout%block_size = stin%block_size
1412 stout%orth_method = stin%orth_method
1414 stout%cl_states_mem = stin%cl_states_mem
1415 stout%pack_states = stin%pack_states
1418 stout%only_userdef_istates = stin%only_userdef_istates
1420 if (.not. exclude_wfns_)
then
1421 safe_allocate_source_a(stout%rho, stin%rho)
1424 stout%uniform_occ = stin%uniform_occ
1427 safe_allocate_source_a(stout%eigenval, stin%eigenval)
1428 safe_allocate_source_a(stout%occ, stin%occ)
1429 safe_allocate_source_a(stout%spin, stin%spin)
1434 stout%group%nblocks = stin%group%nblocks
1436 safe_allocate_source_a(stout%user_def_states, stin%user_def_states)
1438 safe_allocate_source_a(stout%current, stin%current)
1439 safe_allocate_source_a(stout%current_kpt, stin%current_kpt)
1440 safe_allocate_source_a(stout%rho_core, stin%rho_core)
1441 safe_allocate_source_a(stout%frozen_rho, stin%frozen_rho)
1442 safe_allocate_source_a(stout%frozen_tau, stin%frozen_tau)
1443 safe_allocate_source_a(stout%frozen_gdens, stin%frozen_gdens)
1444 safe_allocate_source_a(stout%frozen_ldens, stin%frozen_ldens)
1446 stout%fixed_occ = stin%fixed_occ
1447 stout%restart_fixed_occ = stin%restart_fixed_occ
1449 stout%fixed_spins = stin%fixed_spins
1451 stout%qtot = stin%qtot
1452 stout%val_charge = stin%val_charge
1456 stout%parallel_in_states = stin%parallel_in_states
1458 call mpi_grp_copy(stout%dom_st_kpt_mpi_grp, stin%dom_st_kpt_mpi_grp)
1459 call mpi_grp_copy(stout%st_kpt_mpi_grp, stin%st_kpt_mpi_grp)
1460 call mpi_grp_copy(stout%dom_st_mpi_grp, stin%dom_st_mpi_grp)
1461 safe_allocate_source_a(stout%node, stin%node)
1462 safe_allocate_source_a(stout%st_kpt_task, stin%st_kpt_task)
1464#ifdef HAVE_SCALAPACK
1470 stout%scalapack_compatible = stin%scalapack_compatible
1472 stout%lnst = stin%lnst
1473 stout%st_start = stin%st_start
1474 stout%st_end = stin%st_end
1478 stout%symmetrize_density = stin%symmetrize_density
1482 stout%packed = stin%packed
1484 stout%randomization = stin%randomization
1500 if (st%modelmbparticles%nparticle > 0)
then
1501 safe_deallocate_a(st%mmb_nspindown)
1502 safe_deallocate_a(st%mmb_iyoung)
1503 safe_deallocate_a(st%mmb_proj)
1510 safe_deallocate_a(st%user_def_states)
1512 safe_deallocate_a(st%rho)
1513 safe_deallocate_a(st%eigenval)
1515 safe_deallocate_a(st%current)
1516 safe_deallocate_a(st%current_para)
1517 safe_deallocate_a(st%current_dia)
1518 safe_deallocate_a(st%current_mag)
1519 safe_deallocate_a(st%current_kpt)
1520 safe_deallocate_a(st%rho_core)
1521 safe_deallocate_a(st%frozen_rho)
1522 safe_deallocate_a(st%frozen_tau)
1523 safe_deallocate_a(st%frozen_gdens)
1524 safe_deallocate_a(st%frozen_ldens)
1525 safe_deallocate_a(st%occ)
1526 safe_deallocate_a(st%spin)
1527 safe_deallocate_a(st%kweights)
1530#ifdef HAVE_SCALAPACK
1536 safe_deallocate_a(st%node)
1537 safe_deallocate_a(st%st_kpt_task)
1539 if (st%parallel_in_states)
then
1540 safe_deallocate_a(st%ap%schedule)
1552 class(
mesh_t),
intent(in) :: mesh
1554 integer,
optional,
intent(in) :: ist_start_
1555 integer,
optional,
intent(in) :: ist_end_
1556 integer,
optional,
intent(in) :: ikpt_start_
1557 integer,
optional,
intent(in) :: ikpt_end_
1558 logical,
optional,
intent(in) :: normalized
1561 integer :: ist, ik, id, ist_start, ist_end, jst, ikpt_start, ikpt_end
1562 complex(real64) :: alpha, beta
1563 real(real64),
allocatable :: dpsi(:, :)
1564 complex(real64),
allocatable :: zpsi(:, :), zpsi2(:)
1565 integer :: ikpoint, ip
1568 logical :: normalized_
1579 safe_allocate(dpsi(1:mesh%np, 1:st%d%dim))
1581 safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
1584 select case (st%d%ispin)
1587 do ik = ikpt_start, ikpt_end
1588 ikpoint = st%d%get_kpoint_index(ik)
1589 do ist = ist_start, ist_end
1593 pre_shift = mesh%pv%xlocal-1, &
1594 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1595 normalized = normalized)
1597 if(mesh%parallel_in_domains)
then
1603 call dmf_random(mesh, dpsi(:, 1), normalized = normalized)
1608 zpsi(ip,1) = cmplx(dpsi(ip,1),
m_zero, real64)
1617 pre_shift = mesh%pv%xlocal-1, &
1618 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1619 normalized = normalized)
1621 if(mesh%parallel_in_domains)
then
1627 call zmf_random(mesh, zpsi(:, 1), normalized = normalized)
1639 if (st%fixed_spins)
then
1641 do ik = ikpt_start, ikpt_end
1642 ikpoint = st%d%get_kpoint_index(ik)
1643 do ist = ist_start, ist_end
1647 pre_shift = mesh%pv%xlocal-1, &
1648 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1649 normalized = normalized)
1651 if(mesh%parallel_in_domains)
then
1657 call dmf_random(mesh, dpsi(:, 1), normalized = normalized)
1661 zpsi(ip,1) = cmplx(dpsi(ip,1),
m_zero, real64)
1667 pre_shift = mesh%pv%xlocal-1, &
1668 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1669 normalized = normalized)
1671 if(mesh%parallel_in_domains)
then
1677 call zmf_random(mesh, zpsi(:, 1), normalized = normalized)
1690 safe_allocate(zpsi2(1:mesh%np))
1691 do jst = ist_start, ist - 1
1693 zpsi(1:mesh%np, 1) = zpsi(1:mesh%np, 1) -
zmf_dotp(mesh, zpsi(:, 1), zpsi2)*zpsi2(1:mesh%np)
1695 safe_deallocate_a(zpsi2)
1698 zpsi(1:mesh%np, 2) = zpsi(1:mesh%np, 1)
1702 if (abs(alpha) >
m_zero)
then
1703 beta = cmplx(st%spin(1, ist, ik) / abs(alpha), st%spin(2, ist, ik) / abs(alpha), real64)
1705 zpsi(1:mesh%np, 1) = alpha*zpsi(1:mesh%np, 1)
1706 zpsi(1:mesh%np, 2) = beta*zpsi(1:mesh%np, 2)
1711 do ik = ikpt_start, ikpt_end
1712 do ist = ist_start, ist_end
1716 pre_shift = mesh%pv%xlocal-1, &
1717 post_shift = mesh%pv%np_global - mesh%pv%xlocal - mesh%np + 1, &
1718 normalized = .false.)
1720 if(mesh%parallel_in_domains)
then
1726 call zmf_random(mesh, zpsi(:, id), normalized = .false.)
1742 safe_deallocate_a(dpsi)
1743 safe_deallocate_a(zpsi)
1754 class(
mesh_t),
intent(in) :: mesh
1755 logical,
optional,
intent(in) :: compute_spin
1759 real(real64) :: charge
1760 complex(real64),
allocatable :: zpsi(:, :)
1765 st%nik, st%nst, st%kweights)
1772 charge = charge + sum(st%occ(ist, 1:st%nik) * st%kweights(1:st%nik))
1774 if (abs(charge-st%qtot) > 1e-6_real64)
then
1775 message(1) =
'Occupations do not integrate to total charge.'
1776 write(
message(2),
'(6x,f12.8,a,f12.8)') charge,
' != ', st%qtot
1779 message(1) =
"There don't seem to be any electrons at all!"
1789 safe_allocate(zpsi(1:mesh%np, st%d%dim))
1790 do ik = st%d%kpt%start, st%d%kpt%end
1791 do ist = st%st_start, st%st_end
1793 st%spin(1:3, ist, ik) =
state_spin(mesh, zpsi)
1796 safe_deallocate_a(zpsi)
1798 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
1813 real(real64),
optional,
intent(in) :: alt_eig(st%st_start:, st%d%kpt%start:)
1822 do ik = st%d%kpt%start, st%d%kpt%end
1823 if (
present(alt_eig))
then
1824 tot = tot + st%kweights(ik) * sum(st%occ(st%st_start:st%st_end, ik) * &
1825 alt_eig(st%st_start:st%st_end, ik))
1827 tot = tot + st%kweights(ik) * sum(st%occ(st%st_start:st%st_end, ik) * &
1828 st%eigenval(st%st_start:st%st_end, ik))
1832 if (st%parallel_in_states .or. st%d%kpt%parallel)
call comm_allreduce(st%st_kpt_mpi_grp, tot)
1844 logical :: default_scalapack_compatible
1855 st%parallel_in_states = .false.
1858 call mpi_grp_init(st%dom_st_kpt_mpi_grp, mc%dom_st_kpt_comm)
1862 default_scalapack_compatible =
calc_mode_par%scalapack_compat() .and. .not. st%d%kpt%parallel
1876 call parse_variable(namespace,
'ScaLAPACKCompatible', default_scalapack_compatible, st%scalapack_compatible)
1878#ifdef HAVE_SCALAPACK
1879 if (default_scalapack_compatible .neqv. st%scalapack_compatible)
then
1880 call messages_experimental(
'Setting ScaLAPACKCompatible to other than default', namespace=namespace)
1883 if (st%scalapack_compatible)
then
1890 st%scalapack_compatible = .false.
1899 if (st%nst < st%mpi_grp%size)
then
1900 message(1) =
"Have more processors than necessary"
1901 write(
message(2),
'(i4,a,i4,a)') st%mpi_grp%size,
" processors and ", st%nst,
" states."
1905 call distributed_init(st%dist, st%nst, st%mpi_grp%comm,
"states", scalapack_compat = st%scalapack_compatible)
1907 st%parallel_in_states = st%dist%parallel
1910 st%st_start = st%dist%start
1911 st%st_end = st%dist%end
1912 st%lnst = st%dist%nlocal
1913 st%node(1:st%nst) = st%dist%node(1:st%nst)
1930 kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, &
1931 gi_kinetic_energy_density, st_end)
1935 logical,
intent(in) :: nlcc
1936 real(real64),
contiguous,
optional,
target,
intent(out) :: kinetic_energy_density(:,:)
1937 real(real64),
contiguous,
optional,
target,
intent(out) :: paramagnetic_current(:,:,:)
1938 real(real64),
contiguous,
optional,
intent(out) :: density_gradient(:,:,:)
1939 real(real64),
contiguous,
optional,
intent(out) :: density_laplacian(:,:)
1940 real(real64),
contiguous,
optional,
intent(out) :: gi_kinetic_energy_density(:,:)
1941 integer,
optional,
intent(in) :: st_end
1943 real(real64),
pointer,
contiguous :: jp(:, :, :)
1944 real(real64),
pointer,
contiguous :: tau(:, :)
1945 complex(real64),
allocatable :: wf_psi(:,:), gwf_psi(:,:,:), wf_psi_conj(:,:), lwf_psi(:,:)
1946 real(real64),
allocatable :: abs_wf_psi(:,:), abs_gwf_psi(:,:)
1947 complex(real64),
allocatable :: psi_gpsi(:,:)
1948 complex(real64) :: c_tmp
1949 integer :: is, ik, ist, i_dim, st_dim, ii, st_end_
1950 real(real64) :: ww, kpoint(gr%der%dim)
1951 logical :: something_to_do
1959 something_to_do =
present(kinetic_energy_density) .or.
present(gi_kinetic_energy_density) .or. &
1960 present(paramagnetic_current) .or.
present(density_gradient) .or.
present(density_laplacian)
1961 assert(something_to_do)
1963 safe_allocate( wf_psi(1:gr%np_part, 1:st%d%dim))
1964 safe_allocate( wf_psi_conj(1:gr%np_part, 1:st%d%dim))
1965 safe_allocate(gwf_psi(1:gr%np, 1:gr%der%dim, 1:st%d%dim))
1966 safe_allocate(abs_wf_psi(1:gr%np, 1:st%d%dim))
1967 safe_allocate(abs_gwf_psi(1:gr%np, 1:st%d%dim))
1968 safe_allocate(psi_gpsi(1:gr%np, 1:st%d%dim))
1969 if (
present(density_laplacian))
then
1970 safe_allocate(lwf_psi(1:gr%np, 1:st%d%dim))
1974 if (
present(kinetic_energy_density)) tau => kinetic_energy_density
1977 if (
present(paramagnetic_current)) jp => paramagnetic_current
1981 if (
present(gi_kinetic_energy_density))
then
1983 safe_allocate(jp(1:gr%np, 1:gr%der%dim, 1:st%d%nspin))
1985 if (.not.
present(kinetic_energy_density))
then
1986 safe_allocate(tau(1:gr%np, 1:st%d%nspin))
1990 if (
associated(tau)) tau =
m_zero
1991 if (
associated(jp)) jp =
m_zero
1992 if (
present(density_gradient)) density_gradient(:,:,:) =
m_zero
1993 if (
present(density_laplacian)) density_laplacian(:,:) =
m_zero
1994 if (
present(gi_kinetic_energy_density)) gi_kinetic_energy_density =
m_zero
1996 do ik = st%d%kpt%start, st%d%kpt%end
1998 kpoint(1:gr%der%dim) = kpoints%get_point(st%d%get_kpoint_index(ik))
1999 is = st%d%get_spin_index(ik)
2001 do ist = st%st_start, st_end_
2002 ww = st%kweights(ik)*st%occ(ist, ik)
2008 do st_dim = 1, st%d%dim
2013 do st_dim = 1, st%d%dim
2014 call zderivatives_grad(gr%der, wf_psi(:,st_dim), gwf_psi(:,:,st_dim), set_bc = .false.)
2018 if (
present(density_laplacian))
then
2019 do st_dim = 1, st%d%dim
2020 call zderivatives_lapl(gr%der, wf_psi(:,st_dim), lwf_psi(:,st_dim), ghost_update = .false., set_bc = .false.)
2025 wf_psi_conj(1:gr%np, 1:st%d%dim) = conjg(wf_psi(1:gr%np,1:st%d%dim))
2026 abs_wf_psi(1:gr%np, 1:st%d%dim) = real(wf_psi_conj(1:gr%np, 1:st%d%dim) * wf_psi(1:gr%np, 1:st%d%dim), real64)
2028 if (
present(density_laplacian))
then
2029 density_laplacian(1:gr%np, is) = density_laplacian(1:gr%np, is) + &
2030 ww *
m_two*real(wf_psi_conj(1:gr%np, 1) * lwf_psi(1:gr%np, 1), real64)
2031 if (st%d%ispin ==
spinors)
then
2032 density_laplacian(1:gr%np, 2) = density_laplacian(1:gr%np, 2) + &
2033 ww *
m_two*real(wf_psi_conj(1:gr%np, 2) * lwf_psi(1:gr%np, 2), real64)
2036 c_tmp = ww*(lwf_psi(ii, 1) * wf_psi_conj(ii, 2) + wf_psi(ii, 1) * conjg(lwf_psi(ii, 2)))
2037 density_laplacian(ii, 3) = density_laplacian(ii, 3) + real(c_tmp, real64)
2038 density_laplacian(ii, 4) = density_laplacian(ii, 4) + aimag(c_tmp)
2043 if (
associated(tau))
then
2044 tau(1:gr%np, is) = tau(1:gr%np, is) &
2045 + ww * sum(kpoint(1:gr%der%dim)**2) * abs_wf_psi(1:gr%np, 1)
2046 if (st%d%ispin ==
spinors)
then
2047 tau(1:gr%np, 2) = tau(1:gr%np, 2) &
2048 + ww * sum(kpoint(1:gr%der%dim)**2) * abs_wf_psi(1:gr%np, 2)
2052 c_tmp = ww * sum(kpoint(1:gr%der%dim)**2) * wf_psi(ii, 1) * wf_psi_conj(ii, 2)
2053 tau(ii, 3) = tau(ii, 3) + real(c_tmp, real64)
2054 tau(ii, 4) = tau(ii, 4) + aimag(c_tmp)
2059 do i_dim = 1, gr%der%dim
2062 psi_gpsi(1:gr%np, 1:st%d%dim) = wf_psi_conj(1:gr%np, 1:st%d%dim) &
2063 * gwf_psi(1:gr%np,i_dim,1:st%d%dim)
2064 abs_gwf_psi(1:gr%np, 1:st%d%dim) = real(conjg(gwf_psi(1:gr%np, i_dim, 1:st%d%dim)) &
2065 * gwf_psi(1:gr%np, i_dim, 1:st%d%dim), real64)
2067 if (
present(density_gradient))
then
2068 density_gradient(1:gr%np, i_dim, is) = density_gradient(1:gr%np, i_dim, is) &
2069 + ww *
m_two * real(psi_gpsi(1:gr%np, 1), real64)
2070 if (st%d%ispin ==
spinors)
then
2071 density_gradient(1:gr%np, i_dim, 2) = density_gradient(1:gr%np, i_dim, 2) &
2072 + ww *
m_two*real(psi_gpsi(1:gr%np, 2), real64)
2075 c_tmp = ww * (gwf_psi(ii, i_dim, 1) * wf_psi_conj(ii, 2) + wf_psi(ii, 1) * conjg(gwf_psi(ii, i_dim, 2)))
2076 density_gradient(ii, i_dim, 3) = density_gradient(ii, i_dim, 3) + real(c_tmp, real64)
2077 density_gradient(ii, i_dim, 4) = density_gradient(ii, i_dim, 4) + aimag(c_tmp)
2082 if (
present(density_laplacian))
then
2083 call lalg_axpy(gr%np, ww*
m_two, abs_gwf_psi(:,1), density_laplacian(:,is))
2084 if (st%d%ispin ==
spinors)
then
2085 call lalg_axpy(gr%np, ww*
m_two, abs_gwf_psi(:,2), density_laplacian(:,2))
2088 c_tmp =
m_two * ww * gwf_psi(ii, i_dim, 1) * conjg(gwf_psi(ii, i_dim, 2))
2089 density_laplacian(ii, 3) = density_laplacian(ii, 3) + real(c_tmp, real64)
2090 density_laplacian(ii, 4) = density_laplacian(ii, 4) + aimag(c_tmp)
2098 if (
associated(jp))
then
2100 jp(1:gr%np, i_dim, is) = jp(1:gr%np, i_dim, is) + &
2101 ww*(aimag(psi_gpsi(1:gr%np, 1)) - abs_wf_psi(1:gr%np, 1)*kpoint(i_dim))
2102 if (st%d%ispin ==
spinors)
then
2103 jp(1:gr%np, i_dim, 2) = jp(1:gr%np, i_dim, 2) + &
2104 ww*( aimag(psi_gpsi(1:gr%np, 2)) - abs_wf_psi(1:gr%np, 2)*kpoint(i_dim))
2107 c_tmp = -ww*
m_half*
m_zi*(gwf_psi(ii, i_dim, 1)*wf_psi_conj(ii, 2) - wf_psi(ii, 1)*conjg(gwf_psi(ii, i_dim, 2)) &
2108 -
m_two *
m_zi*wf_psi(ii, 1)*wf_psi_conj(ii, 2)*kpoint(i_dim))
2109 jp(ii, i_dim, 3) = jp(ii, i_dim, 3) + real(c_tmp, real64)
2110 jp(ii, i_dim, 4) = jp(ii, i_dim, 4) + aimag(c_tmp)
2119 if (
associated(tau))
then
2120 tau(1:gr%np, is) = tau(1:gr%np, is) + ww*(abs_gwf_psi(1:gr%np,1) &
2121 -
m_two*aimag(psi_gpsi(1:gr%np, 1))*kpoint(i_dim))
2122 if (st%d%ispin ==
spinors)
then
2123 tau(1:gr%np, 2) = tau(1:gr%np, 2) + ww*(abs_gwf_psi(1:gr%np, 2) &
2124 -
m_two*aimag(psi_gpsi(1:gr%np, 2))*kpoint(i_dim))
2127 c_tmp = ww * ( gwf_psi(ii, i_dim, 1)*conjg(gwf_psi(ii, i_dim, 2)) &
2128 +
m_zi * (gwf_psi(ii,i_dim,1)*wf_psi_conj(ii, 2) &
2129 - wf_psi(ii, 1)*conjg(gwf_psi(ii,i_dim,2)))*kpoint(i_dim))
2130 tau(ii, 3) = tau(ii, 3) + real(c_tmp, real64)
2131 tau(ii, 4) = tau(ii, 4) + aimag(c_tmp)
2141 safe_deallocate_a(wf_psi_conj)
2142 safe_deallocate_a(abs_wf_psi)
2143 safe_deallocate_a(abs_gwf_psi)
2144 safe_deallocate_a(psi_gpsi)
2146 if (.not.
present(gi_kinetic_energy_density))
then
2147 if (.not.
present(paramagnetic_current))
then
2148 safe_deallocate_p(jp)
2150 if (.not.
present(kinetic_energy_density))
then
2151 safe_deallocate_p(tau)
2155 if (st%parallel_in_states .or. st%d%kpt%parallel)
call reduce_all(st%st_kpt_mpi_grp)
2160 if (st%symmetrize_density)
then
2161 do is = 1, st%d%nspin
2162 if (
associated(tau))
then
2166 if (
present(density_laplacian))
then
2170 if (
associated(jp))
then
2174 if (
present(density_gradient))
then
2181 if (
allocated(st%rho_core) .and. nlcc .and. (
present(density_laplacian) .or.
present(density_gradient)))
then
2183 wf_psi(ii, 1) = st%rho_core(ii)/st%d%spin_channels
2188 if (
present(density_gradient))
then
2191 do is = 1, st%d%spin_channels
2192 density_gradient(1:gr%np, 1:gr%der%dim, is) = density_gradient(1:gr%np, 1:gr%der%dim, is) + &
2193 real(gwf_psi(1:gr%np, 1:gr%der%dim,1))
2198 if (
present(density_laplacian))
then
2201 do is = 1, st%d%spin_channels
2202 density_laplacian(1:gr%np, is) = density_laplacian(1:gr%np, is) + real(lwf_psi(1:gr%np, 1))
2209 if (
allocated(st%frozen_tau) .and. .not.
present(st_end))
then
2212 if (
allocated(st%frozen_gdens) .and. .not.
present(st_end))
then
2213 do is = 1, st%d%nspin
2214 call lalg_axpy(gr%np, gr%der%dim,
m_one, st%frozen_gdens(:,:,is), density_gradient(:,:,is))
2217 if (
allocated(st%frozen_tau) .and. .not.
present(st_end))
then
2218 call lalg_axpy(gr%np, st%d%nspin,
m_one, st%frozen_ldens, density_laplacian)
2221 safe_deallocate_a(wf_psi)
2222 safe_deallocate_a(gwf_psi)
2223 safe_deallocate_a(lwf_psi)
2227 if (
present(gi_kinetic_energy_density))
then
2228 do is = 1, st%d%nspin
2229 assert(
associated(tau))
2230 gi_kinetic_energy_density(1:gr%np, is) = tau(1:gr%np, is)
2233 assert(
associated(jp))
2234 if (st%d%ispin /=
spinors)
then
2235 do is = 1, st%d%nspin
2238 if (st%rho(ii, is) < 1.0e-7_real64) cycle
2239 gi_kinetic_energy_density(ii, is) = &
2240 gi_kinetic_energy_density(ii, is) - sum(jp(ii,1:gr%der%dim, is)**2)/st%rho(ii, is)
2246 gi_kinetic_energy_density(ii, 1) = gi_kinetic_energy_density(ii, 1) &
2247 - sum(jp(ii,1:gr%der%dim, 1)**2)/(safe_tol(st%rho(ii, 1),
m_epsilon))
2248 gi_kinetic_energy_density(ii, 2) = gi_kinetic_energy_density(ii, 2) &
2249 - sum(jp(ii,1:gr%der%dim, 2)**2)/(safe_tol(st%rho(ii, 2),
m_epsilon))
2250 gi_kinetic_energy_density(ii, 3) = &
2251 gi_kinetic_energy_density(ii, 3) - sum(jp(ii,1:gr%der%dim, 3)**2 + jp(ii,1:gr%der%dim, 4)**2) &
2252 /(safe_tol((st%rho(ii, 3)**2 + st%rho(ii, 4)**2),
m_epsilon))*st%rho(ii, 3)
2253 gi_kinetic_energy_density(ii, 4) = &
2254 gi_kinetic_energy_density(ii, 4) + sum(jp(ii,1:gr%der%dim, 3)**2 + jp(ii,1:gr%der%dim, 4)**2) &
2255 /(safe_tol((st%rho(ii, 3)**2 + st%rho(ii, 4)**2),
m_epsilon))*st%rho(ii, 4)
2261 if (.not.
present(kinetic_energy_density))
then
2262 safe_deallocate_p(tau)
2264 if (.not.
present(paramagnetic_current))
then
2265 safe_deallocate_p(jp)
2280 if (
associated(tau))
call comm_allreduce(grp, tau, dim = (/gr%np, st%d%nspin/))
2282 if (
present(density_laplacian))
call comm_allreduce(grp, density_laplacian, dim = (/gr%np, st%d%nspin/))
2284 do is = 1, st%d%nspin
2285 if (
associated(jp))
call comm_allreduce(grp, jp(:, :, is), dim = (/gr%np, gr%der%dim/))
2287 if (
present(density_gradient))
then
2288 call comm_allreduce(grp, density_gradient(:, :, is), dim = (/gr%np, gr%der%dim/))
2302 type(
mesh_t),
intent(in) :: mesh
2303 complex(real64),
intent(in) :: f1(:, :)
2304 real(real64) :: spin(1:3)
2306 complex(real64) :: z
2310 z =
zmf_dotp(mesh, f1(:, 1) , f1(:, 2))
2312 spin(1) =
m_two*real(z, real64)
2313 spin(2) =
m_two*aimag(z)
2325 integer,
intent(in) :: ist
2339 integer,
intent(in) :: ist
2340 integer,
intent(in) :: ik
2345 ik >= st%d%kpt%start .and. ik <= st%d%kpt%end
2355 class(
mesh_t),
intent(in) :: mesh
2361 memory = memory + real64*real(mesh%np_part_global, real64) *st%d%dim*real(st%nst, real64) *st%d%kpt%nglobal
2371 logical,
optional,
intent(in) :: copy
2374 integer(int64) :: max_mem, mem
2386 if (accel_is_enabled())
then
2387 max_mem = accel_global_memory_size()
2389 if (st%cl_states_mem > m_one)
then
2390 max_mem = int(st%cl_states_mem, int64)*(1024_8)**2
2391 else if (st%cl_states_mem < 0.0_real64)
then
2392 max_mem = max_mem + int(st%cl_states_mem, int64)*(1024_8)**2
2394 max_mem = int(st%cl_states_mem*real(max_mem, real64) , int64)
2397 max_mem = huge(max_mem)
2401 qnloop:
do iqn = st%d%kpt%start, st%d%kpt%end
2402 do ib = st%group%block_start, st%group%block_end
2404 mem = mem + st%group%psib(ib, iqn)%pack_total_size()
2406 if (mem > max_mem)
then
2407 call messages_write(
'Not enough CL device memory to store all states simultaneously.', new_line = .
true.)
2408 call messages_write(
'Only ')
2409 call messages_write(ib - st%group%block_start)
2410 call messages_write(
' of ')
2411 call messages_write(st%group%block_end - st%group%block_start + 1)
2412 call messages_write(
' blocks will be stored in device memory.', new_line = .
true.)
2413 call messages_warning()
2417 call st%group%psib(ib, iqn)%do_pack(copy)
2429 logical,
optional,
intent(in) :: copy
2438 do iqn = st%d%kpt%start, st%d%kpt%end
2439 do ib = st%group%block_start, st%group%block_end
2440 if (st%group%psib(ib, iqn)%is_packed())
call st%group%psib(ib, iqn)%do_unpack(copy)
2453 type(namespace_t),
intent(in) :: namespace
2457 call messages_print_with_emphasis(msg=
"States", namespace=namespace)
2459 write(message(1),
'(a,f12.3)')
'Total electronic charge = ', st%qtot
2460 write(message(2),
'(a,i8)')
'Number of states = ', st%nst
2461 write(message(3),
'(a,i8)')
'States block-size = ', st%block_size
2462 call messages_info(3, namespace=namespace)
2464 call messages_print_with_emphasis(namespace=namespace)
2479 do iqn = st%d%kpt%start, st%d%kpt%end
2480 do ib = st%group%block_start, st%group%block_end
2481 call batch_set_zero(st%group%psib(ib, iqn))
2491 integer pure function states_elec_block_min(st, ib) result(range)
2493 integer,
intent(in) :: ib
2495 range = st%group%block_range(ib, 1)
2501 integer pure function states_elec_block_max(st, ib) result(range)
2503 integer,
intent(in) :: ib
2505 range = st%group%block_range(ib, 2)
2511 integer pure function states_elec_block_size(st, ib) result(size)
2513 integer,
intent(in) :: ib
2515 size = st%group%block_size(ib)
2523 type(namespace_t),
intent(in) :: namespace
2524 integer,
intent(out) :: n_pairs
2525 integer,
intent(out) :: n_occ(:)
2526 integer,
intent(out) :: n_unocc(:)
2527 logical,
allocatable,
intent(out) :: is_included(:,:,:)
2529 logical,
intent(out) :: is_frac_occ
2531 integer :: ik, ist, ast, n_filled, n_partially_filled, n_half_filled
2532 character(len=80) :: nst_string, default, wfn_list
2533 real(real64) :: energy_window
2537 is_frac_occ = .false.
2539 call occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled)
2540 if (n_partially_filled > 0 .or. n_half_filled > 0) is_frac_occ = .
true.
2541 n_occ(ik) = n_filled + n_partially_filled + n_half_filled
2542 n_unocc(ik) = st%nst - n_filled
2555 call parse_variable(namespace,
'CasidaKSEnergyWindow', -m_one, energy_window, units_inp%energy)
2578 safe_allocate(is_included(maxval(n_occ), minval(n_occ) + 1:st%nst , st%nik))
2579 is_included(:,:,:) = .false.
2581 if (energy_window < m_zero)
then
2582 write(nst_string,
'(i6)') st%nst
2583 write(default,
'(a,a)')
"1-", trim(adjustl(nst_string))
2584 call parse_variable(namespace,
'CasidaKohnShamStates', default, wfn_list)
2586 write(message(1),
'(a,a)')
"Info: States that form the basis: ", trim(wfn_list)
2587 call messages_info(1, namespace=namespace)
2592 do ast = n_occ(ik) + 1, st%nst
2593 if (loct_isinstringlist(ast, wfn_list))
then
2594 do ist = 1, n_occ(ik)
2595 if (loct_isinstringlist(ist, wfn_list))
then
2596 n_pairs = n_pairs + 1
2597 is_included(ist, ast, ik) = .
true.
2606 write(message(1),
'(a,f12.6,a)')
"Info: including transitions with energy < ", &
2607 units_from_atomic(units_out%energy, energy_window), trim(units_abbrev(units_out%energy))
2608 call messages_info(1, namespace=namespace)
2613 do ast = n_occ(ik) + 1, st%nst
2614 do ist = 1, n_occ(ik)
2615 if (st%eigenval(ast, ik) - st%eigenval(ist, ik) < energy_window)
then
2616 n_pairs = n_pairs + 1
2617 is_included(ist, ast, ik) = .
true.
2645 subroutine occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled, &
2646 filled, partially_filled, half_filled)
2648 type(namespace_t),
intent(in) :: namespace
2649 integer,
intent(in) :: ik
2650 integer,
intent(out) :: n_filled, n_partially_filled, n_half_filled
2651 integer,
optional,
intent(out) :: filled(:), partially_filled(:), half_filled(:)
2657 if (
present(filled)) filled(:) = 0
2658 if (
present(partially_filled)) partially_filled(:) = 0
2659 if (
present(half_filled)) half_filled(:) = 0
2661 n_partially_filled = 0
2664 select case (st%d%ispin)
2667 if (abs(st%occ(ist, ik) - m_two) < m_min_occ)
then
2668 n_filled = n_filled + 1
2669 if (
present(filled)) filled(n_filled) = ist
2670 else if (abs(st%occ(ist, ik) - m_one) < m_min_occ)
then
2671 n_half_filled = n_half_filled + 1
2672 if (
present(half_filled)) half_filled(n_half_filled) = ist
2673 else if (st%occ(ist, ik) > m_min_occ)
then
2674 n_partially_filled = n_partially_filled + 1
2675 if (
present(partially_filled)) partially_filled(n_partially_filled) = ist
2676 else if (abs(st%occ(ist, ik)) > m_min_occ)
then
2677 write(message(1),*)
'Internal error in occupied_states: Illegal occupation value ', st%occ(ist, ik)
2678 call messages_fatal(1, namespace=namespace)
2681 case (spin_polarized, spinors)
2683 if (abs(st%occ(ist, ik)-m_one) < m_min_occ)
then
2684 n_filled = n_filled + 1
2685 if (
present(filled)) filled(n_filled) = ist
2686 else if (st%occ(ist, ik) > m_min_occ)
then
2687 n_partially_filled = n_partially_filled + 1
2688 if (
present(partially_filled)) partially_filled(n_partially_filled) = ist
2689 else if (abs(st%occ(ist, ik)) > m_min_occ)
then
2690 write(message(1),*)
'Internal error in occupied_states: Illegal occupation value ', st%occ(ist, ik)
2691 call messages_fatal(1, namespace=namespace)
2703 type(multicomm_t),
intent(in) :: mc
2706 call distributed_init(this%d%kpt, this%nik, mc%group_comm(p_strategy_kpoints),
"k-points")
2721 if (.not.
allocated(st%st_kpt_task))
then
2722 safe_allocate(st%st_kpt_task(0:st%st_kpt_mpi_grp%size-1, 1:4))
2725 st%st_kpt_task(0:st%st_kpt_mpi_grp%size-1, :) = 0
2726 st%st_kpt_task(st%st_kpt_mpi_grp%rank, :) = [st%st_start, st%st_end, st%d%kpt%start, st%d%kpt%end]
2728 if (st%parallel_in_states .or. st%d%kpt%parallel)
then
2729 call comm_allreduce(st%st_kpt_mpi_grp, st%st_kpt_task)
2740 type(kpoints_t),
intent(in) :: kpoints
2741 type(namespace_t),
intent(in) :: namespace
2744 type(states_elec_dim_t),
pointer :: dd
2751 st%nik = kpoints_number(kpoints)
2753 if (dd%ispin == spin_polarized) st%nik = 2*st%nik
2755 safe_allocate(st%kweights(1:st%nik))
2758 ik = dd%get_kpoint_index(iq)
2759 st%kweights(iq) = kpoints%get_weight(ik)
2771 call io_mkdir(
'debug/', namespace)
2772 iunit = io_open(
'debug/kpoints', namespace, action =
'write')
2773 call kpoints%write_info(iunit=iunit, absolute_coordinates = .
true.)
2774 call io_close(iunit)
2788 class(mesh_t),
intent(in) :: gr
2789 real(real64) :: dipole(1:gr%box%dim)
2792 real(real64) :: e_dip(1:gr%box%dim, this%d%spin_channels)
2796 do ispin = 1, this%d%spin_channels
2797 call dmf_dipole(gr, this%rho(:, ispin), e_dip(:, ispin))
2800 dipole(:) = sum(e_dip(:, 1:this%d%spin_channels), 2)
2808#include "states_elec_inc.F90"
2811#include "complex.F90"
2812#include "states_elec_inc.F90"
initialize a batch with existing memory
constant times a vector plus a vector
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
pure logical function, public accel_is_enabled()
type(accel_t), public accel
This module implements batches of mesh functions.
This module implements common operations on batches of mesh functions.
This module provides the BLACS processor grid.
subroutine, public blacs_proc_grid_init(this, mpi_grp, procdim)
Initializes a blacs context from an MPI communicator with topological information.
subroutine, public blacs_proc_grid_end(this)
subroutine, public blacs_proc_grid_copy(cin, cout)
Module implementing boundary conditions in Octopus.
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public zderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public zderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
subroutine, public distributed_end(this)
subroutine, public distributed_nullify(this, total)
subroutine, public distributed_init(this, total, comm, tag, scalapack_compat)
Distribute N instances across M processes of communicator comm
subroutine, public distributed_copy(in, out)
@Brief Create a copy of a distributed instance
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_fourth
complex(real64), parameter, public m_zi
real(real64), parameter, public m_epsilon
type(conf_t), public conf
Global instance of Octopus configuration.
real(real64), parameter, public m_half
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
subroutine, public dgrid_symmetrize_vector_field(gr, field, suppress_warning)
subroutine, public dgrid_symmetrize_scalar_field(gr, field, suppress_warning)
logical pure function, public kpoints_point_is_gamma(this, ik)
This module is intended to contain "only mathematical" functions and procedures.
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
This module defines functions over batches of mesh functions.
subroutine, public zmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
subroutine, public dmesh_batch_exchange_points(mesh, aa, forward_map, backward_map)
This functions exchanges points of a mesh according to a certain map. Two possible maps can be given....
This module defines various routines, operating on mesh functions.
subroutine, public zmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
subroutine, public dmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
subroutine, public zmf_normalize(mesh, dim, psi, norm)
Normalize a mesh function psi.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_not_implemented(feature, namespace)
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
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)
general module for modelmb particles
subroutine, public modelmb_particles_end(this)
subroutine, public modelmb_particles_init(this, namespace, space)
==============================================================
subroutine, public modelmb_particles_copy(modelmb_out, modelmb_in)
subroutine mpi_grp_copy(mpi_grp_out, mpi_grp_in)
type(mpi_comm), parameter, public mpi_comm_undefined
used to indicate a communicator has not been initialized
subroutine mpi_grp_init(grp, comm)
Initialize MPI group instance.
This module handles the communicators for the various parallelization strategies.
logical pure function, public multicomm_strategy_is_parallel(mc, level)
subroutine, public multicomm_all_pairs_copy(apout, apin)
logical pure function, public multicomm_have_slaves(this)
subroutine, public multicomm_create_all_pairs(mpi_grp, ap)
This routine uses the one-factorization (or near-one-factorization of a complete graph to construct a...
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.
subroutine, public restart_read(restart, iunit, lines, nlines, ierr)
subroutine, public restart_close(restart, iunit)
Close a file previously opened with restart_open.
integer function, public restart_open(restart, filename, status, position, silent)
Open file "filename" found inside the current restart directory. Depending on the type of restart,...
subroutine, public smear_fill_occupations(this, eigenvalues, occupations, nik, nst)
subroutine, public smear_find_fermi_energy(this, namespace, eigenvalues, occupations, qtot, nik, nst, kweights)
subroutine, public smear_copy(to, from)
integer, parameter, public smear_fixed_occ
subroutine, public smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
logical pure function, public smear_is_semiconducting(this)
subroutine, public states_set_complex(st)
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_dim_copy(dout, din)
subroutine, public states_elec_dim_end(dim)
This module handles groups of electronic batches and their parallel distribution.
subroutine, public states_elec_group_copy(d, group_in, group_out, copy_data, special)
make a copy of a group
subroutine, public states_elec_group_end(this, d)
finalize the local blocks of wave functions and release local arrays
real(real64) function, public states_elec_wfns_memory(st, mesh)
return the memory usage of a states_elec_t object
integer pure function, public states_elec_block_max(st, ib)
return index of last state in block ib
subroutine zstates_elec_generate_random_vector(mesh, st, vector, normalized, reset_seed)
Generate a random vector.
subroutine, public states_elec_null(st)
subroutine, public states_elec_distribute_nodes(st, namespace, mc)
@Brief. Distribute states over the processes for states parallelization
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
logical function, public state_kpt_is_local(st, ist, ik)
check whether a given state (ist, ik) is on the local node
subroutine, public states_elec_calc_quantities(gr, st, kpoints, nlcc, kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, gi_kinetic_energy_density, st_end)
calculated selected quantities
real(real64) function, dimension(1:gr%box%dim) states_elec_calculate_dipole(this, gr)
calculate the expectation value of the dipole moment of electrons
subroutine, public states_elec_densities_init(st, gr)
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine zstates_elec_get_state1(st, mesh, idim, ist, iqn, psi)
Write one component (dim) of a wave function into a state_elec_t object.
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine zstates_elec_get_points1(st, start_point, end_point, iqn, psi)
Returns the value of all the states for given k-point in the range of points [start_point:end_point].
subroutine states_elec_read_initial_occs(st, namespace, excess_charge, kpoints)
Reads from the input file the initial occupations, if the block "Occupations" is present....
subroutine, public states_elec_init(st, namespace, space, valence_charge, kpoints)
Initialize a new states_elec_t object.
subroutine zstates_elec_set_state1(st, mesh, idim, ist, iqn, psi)
get one dimension of local wave function for given k-point and states index from a states_elec_t obje...
real(real64) function, public states_elec_eigenvalues_sum(st, alt_eig)
function to calculate the eigenvalues sum using occupations as weights
integer, parameter, public par_independent
Method used to generate random states.
subroutine dstates_elec_generate_random_vector(mesh, st, vector, normalized, reset_seed)
Generate a random vector.
subroutine, public occupied_states(st, namespace, ik, n_filled, n_partially_filled, n_half_filled, filled, partially_filled, half_filled)
return information about occupied orbitals in many-body state
subroutine dstates_elec_get_state2(st, mesh, ist, iqn, psi)
Write a wave function into a state_elec_t object.
subroutine zstates_elec_set_state2(st, mesh, ist, iqn, psi)
get local wave function for given k-point and states index from a states_elec_t object
subroutine zstates_elec_get_state2(st, mesh, ist, iqn, psi)
Write a wave function into a state_elec_t object.
subroutine dstates_elec_get_points1(st, start_point, end_point, iqn, psi)
Returns the value of all the states for given k-point in the range of points [start_point:end_point].
integer pure function, public states_elec_block_size(st, ib)
return number of states in block ib
subroutine, public kpoints_distribute(this, mc)
distribute k-points over the nodes in the corresponding communicator
integer, parameter, public par_dependent
subroutine states_elec_pack(st, copy)
pack the batches in this states object
subroutine, public states_elec_choose_kpoints(st, kpoints, namespace)
double up k-points for SPIN_POLARIZED calculations
subroutine zstates_elec_get_state3(st, mesh, iqn, psi)
get local wave functions for given k-point from a states_elec_t object
logical function, public state_is_local(st, ist)
check whether a given state (ist) is on the local node
subroutine dstates_elec_set_state1(st, mesh, idim, ist, iqn, psi)
get one dimension of local wave function for given k-point and states index from a states_elec_t obje...
subroutine, public states_elec_exec_init(st, namespace, mc)
Further initializations.
subroutine dstates_elec_set_state4(st, mesh, psi)
set all local wave functions in a states_elec_t object
subroutine dstates_elec_set_state3(st, mesh, iqn, psi)
set local wave functions for given k-point in a states_elec_t object
subroutine states_elec_write_info(st, namespace)
write information about the states object
subroutine states_elec_init_block(st, mesh, verbose, skip, packed)
Initializes the data components in st that describe how the states are distributed in blocks:
real(real64) function, dimension(1:3) state_spin(mesh, f1)
calculate the spin vector for a spinor wave function f1
subroutine states_elec_kpoints_distribution(st)
Assign the start and end indices for states and kpoints, for "st_kpt_mpi_grp" communicator.
subroutine zstates_elec_get_state4(st, mesh, psi)
get all local wave functions from a states_elec_t object
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
randomize states
subroutine dstates_elec_get_state1(st, mesh, idim, ist, iqn, psi)
Write one component (dim) of a wave function into a state_elec_t object.
subroutine dstates_elec_set_state2(st, mesh, ist, iqn, psi)
get local wave function for given k-point and states index from a states_elec_t object
subroutine zstates_elec_set_state4(st, mesh, psi)
set all local wave functions in a states_elec_t object
subroutine states_elec_read_initial_spins(st, namespace)
Reads, if present, the "InitialSpins" block.
subroutine dstates_elec_get_state3(st, mesh, iqn, psi)
get local wave functions for given k-point from a states_elec_t object
subroutine dstates_elec_get_points2(st, start_point, end_point, psi)
Returns the value of all the states in the range of points [start_point:end_point].
integer pure function, public states_elec_block_min(st, ib)
return index of first state in block ib
subroutine states_elec_unpack(st, copy)
unpack the batches in this states object
subroutine, public states_elec_look(restart, nik, dim, nst, ierr)
Reads the 'states' file in the restart directory, and finds out the nik, dim, and nst contained in it...
subroutine dstates_elec_get_state4(st, mesh, psi)
get all local wave functions from a states_elec_t object
subroutine, public states_elec_count_pairs(st, namespace, n_pairs, n_occ, n_unocc, is_included, is_frac_occ)
number of occupied-unoccupied pairs for Casida
subroutine, public states_elec_allocate_current(st, space, mesh)
subroutine, public states_elec_set_zero(st)
expclicitely set all wave functions in the states to zero
subroutine zstates_elec_set_state3(st, mesh, iqn, psi)
set local wave functions for given k-point in a states_elec_t object
subroutine zstates_elec_get_points2(st, start_point, end_point, psi)
Returns the value of all the states in the range of points [start_point:end_point].
type(type_t), public type_float
type(type_t), public type_cmplx
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_t), public unit_megabytes
For large amounts of data (natural code units are bytes)
subroutine, public zwfs_elec_init(this, dim, st_start, st_end, np, ik, special, packed)
initialize an empty wfs_elec_t object
subroutine, public dwfs_elec_init(this, dim, st_start, st_end, np, ik, special, packed)
initialize an empty wfs_elec_t object
subroutine reduce_all(grp)
subroutine print_kpoints_debug
Class defining batches of mesh functions.
Distribution of N instances over mpi_grpsize processes, for the local rank mpi_grprank....
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Describes mesh distribution to nodes.
This is defined even when running serial.
An all-pairs communication schedule for a given group.
Stores all communicators and groups.
abstract class for states
The states_elec_t class contains all electronic wave functions.