52 USE dbcsr_api,
ONLY: dbcsr_add,&
200#include "./base/base_uses.f90"
206 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_scf_post_gpw'
239 CHARACTER(6),
OPTIONAL :: wf_type
240 LOGICAL,
OPTIONAL :: do_mp2
242 CHARACTER(len=*),
PARAMETER :: routinen =
'scf_post_calculation_gpw'
244 INTEGER :: handle, homo, ispin, min_lumos, n_rep, nchk_nmoloc, nhomo, nlumo, nlumo_stm, &
245 nlumo_tddft, nlumos, nmo, nspins, output_unit, unit_nr
246 INTEGER,
DIMENSION(:, :, :),
POINTER :: marked_states
247 LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints,
do_mixed, do_mo_cubes, do_stm, &
248 do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_do_mp2, &
249 my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo, p_loc_mixed
251 REAL(kind=
dp) :: gap, homo_lumo(2, 2), total_zeff_corr
252 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_eigenvalues
255 TYPE(
cp_1d_r_p_type),
DIMENSION(:),
POINTER :: mixed_evals, occupied_evals, &
256 unoccupied_evals, unoccupied_evals_stm
257 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: mixed_orbs, occupied_orbs
258 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
259 TARGET :: homo_localized, lumo_localized, &
261 TYPE(
cp_fm_type),
DIMENSION(:),
POINTER :: lumo_ptr, mo_loc_history, &
262 unoccupied_orbs, unoccupied_orbs_stm
265 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_rmpv, matrix_p_mp2, matrix_s, &
267 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: kinetic_m, rho_ao
279 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
287 localize_section, print_key, &
290 CALL timeset(routinen, handle)
297 IF (
PRESENT(do_mp2)) my_do_mp2 = do_mp2
298 IF (
PRESENT(wf_type))
THEN
299 IF (output_unit > 0)
THEN
300 WRITE (unit=output_unit, fmt=
'(/,(T1,A))') repeat(
"-", 40)
301 WRITE (unit=output_unit, fmt=
'(/,(T3,A,T19,A,T25,A))')
"Properties from ", wf_type,
" density"
302 WRITE (unit=output_unit, fmt=
'(/,(T1,A))') repeat(
"-", 40)
308 CALL write_available_results(qs_env, scf_env)
310 my_localized_wfn = .false.
311 NULLIFY (admm_env, dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, &
312 mo_coeff, ks_rmpv, matrix_s, qs_loc_env_homo, qs_loc_env_lumo, scf_control, &
313 unoccupied_orbs, mo_eigenvalues, unoccupied_evals, &
314 unoccupied_evals_stm, molecule_set, mo_derivs, &
315 subsys, particles, input, print_key, kinetic_m, marked_states, &
316 mixed_evals, qs_loc_env_mixed)
317 NULLIFY (lumo_ptr, rho_ao)
324 p_loc_mixed = .false.
326 cpassert(
ASSOCIATED(scf_env))
327 cpassert(
ASSOCIATED(qs_env))
330 dft_control=dft_control, &
331 molecule_set=molecule_set, &
332 scf_control=scf_control, &
333 do_kpoints=do_kpoints, &
338 particle_set=particle_set, &
339 atomic_kind_set=atomic_kind_set, &
340 qs_kind_set=qs_kind_set)
347 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
348 DO ispin = 1, dft_control%nspins
349 CALL dbcsr_add(rho_ao(ispin, 1)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
354 CALL update_hartree_with_mp2(rho, qs_env)
359 "DFT%PRINT%KINETIC_ENERGY") /= 0)
THEN
361 cpassert(
ASSOCIATED(kinetic_m))
362 cpassert(
ASSOCIATED(kinetic_m(1, 1)%matrix))
366 IF (unit_nr > 0)
THEN
367 WRITE (unit_nr,
'(T3,A,T55,F25.14)')
"Electronic kinetic energy:", e_kin
370 "DFT%PRINT%KINETIC_ENERGY")
374 CALL qs_scf_post_charges(input, logger, qs_env)
387 IF (loc_print_explicit)
THEN
405 IF (loc_explicit)
THEN
415 p_loc_mixed = .false.
419 IF (n_rep == 0 .AND. p_loc_lumo)
THEN
420 CALL cp_abort(__location__,
"No LIST_UNOCCUPIED was specified, "// &
421 "therefore localization of unoccupied states will be skipped!")
436 IF (loc_print_explicit)
THEN
440 do_wannier_cubes = .false.
445 IF (dft_control%do_tddfpt_calculation)
THEN
450 IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc)
THEN
451 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
453 CALL auxbas_pw_pool%create_pw(wf_r)
454 CALL auxbas_pw_pool%create_pw(wf_g)
457 IF (dft_control%restricted)
THEN
461 nspins = dft_control%nspins
464 IF (dft_control%restricted .AND. (do_mo_cubes .OR. p_loc_homo))
THEN
465 CALL cp_abort(__location__,
"Unclear how we define MOs / localization in the restricted case ... ")
470 IF (do_mo_cubes)
THEN
471 cpwarn(
"Print MO cubes not implemented for k-point calculations")
477 IF ((do_mo_cubes .AND. nhomo /= 0) .OR. do_stm .OR. dft_control%do_tddfpt_calculation)
THEN
479 IF (dft_control%do_admm)
THEN
481 CALL make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env=admm_env)
483 CALL make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs)
485 DO ispin = 1, dft_control%nspins
486 CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=homo)
487 homo_lumo(ispin, 1) = mo_eigenvalues(homo)
491 IF (do_mo_cubes .AND. nhomo /= 0)
THEN
494 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
495 eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
496 CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, &
497 mo_coeff, wf_g, wf_r, particles, homo, ispin)
508 cpwarn(
"Localization not implemented for k-point calculations!!")
509 ELSEIF (dft_control%restricted &
512 cpabort(
"ROKS works only with LOCALIZE METHOD NONE or JACOBI")
514 ALLOCATE (occupied_orbs(dft_control%nspins))
515 ALLOCATE (occupied_evals(dft_control%nspins))
516 ALLOCATE (homo_localized(dft_control%nspins))
517 DO ispin = 1, dft_control%nspins
518 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
519 eigenvalues=mo_eigenvalues)
520 occupied_orbs(ispin) = mo_coeff
521 occupied_evals(ispin)%array => mo_eigenvalues
522 CALL cp_fm_create(homo_localized(ispin), occupied_orbs(ispin)%matrix_struct)
523 CALL cp_fm_to_fm(occupied_orbs(ispin), homo_localized(ispin))
526 CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history)
529 ALLOCATE (qs_loc_env_homo)
532 CALL qs_loc_init(qs_env, qs_loc_env_homo, localize_section, homo_localized, do_homo, &
533 do_mo_cubes, mo_loc_history=mo_loc_history)
535 wf_r, wf_g, particles, occupied_orbs, occupied_evals, marked_states)
538 IF (qs_loc_env_homo%localized_wfn_control%use_history)
THEN
540 CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history)
545 homo_localized, do_homo)
547 DEALLOCATE (occupied_orbs)
548 DEALLOCATE (occupied_evals)
550 IF (qs_loc_env_homo%do_localize)
THEN
551 CALL loc_dipole(input, dft_control, qs_loc_env_homo, logger, qs_env)
558 IF (do_mo_cubes .OR. p_loc_lumo)
THEN
560 cpwarn(
"Localization and MO related output not implemented for k-point calculations!")
563 IF (nlumo .GT. -1)
THEN
564 nlumo = max(nlumo, nlumo_tddft)
566 compute_lumos = (do_mo_cubes .OR. dft_control%do_tddfpt_calculation) .AND. nlumo .NE. 0
567 compute_lumos = compute_lumos .OR. p_loc_lumo
569 DO ispin = 1, dft_control%nspins
570 CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo)
571 compute_lumos = compute_lumos .AND. homo == nmo
574 IF (do_mo_cubes .AND. .NOT. compute_lumos)
THEN
577 DO ispin = 1, dft_control%nspins
579 CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo, eigenvalues=mo_eigenvalues)
580 IF (nlumo > nmo - homo)
THEN
583 IF (nlumo .EQ. -1)
THEN
586 IF (output_unit > 0)
WRITE (output_unit, *)
" "
587 IF (output_unit > 0)
WRITE (output_unit, *)
" Lowest eigenvalues of the unoccupied subspace spin ", ispin
588 IF (output_unit > 0)
WRITE (output_unit, *)
"---------------------------------------------"
589 IF (output_unit > 0)
WRITE (output_unit,
'(4(1X,1F16.8))') mo_eigenvalues(homo + 1:homo + nlumo)
592 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
593 CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
594 mo_coeff, wf_g, wf_r, particles, nlumo, homo, ispin, lumo=homo + 1)
600 IF (compute_lumos)
THEN
603 IF (nlumo == 0) check_write = .false.
606 ALLOCATE (qs_loc_env_lumo)
609 min_lumos = max(maxval(qs_loc_env_lumo%localized_wfn_control%loc_states(:, :)), nlumo)
612 ALLOCATE (unoccupied_orbs(dft_control%nspins))
613 ALLOCATE (unoccupied_evals(dft_control%nspins))
614 CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, min_lumos, nlumos)
615 lumo_ptr => unoccupied_orbs
616 DO ispin = 1, dft_control%nspins
618 homo_lumo(ispin, 2) = unoccupied_evals(ispin)%array(1)
620 IF (check_write)
THEN
621 IF (p_loc_lumo .AND. nlumo .NE. -1) nlumos = min(nlumo, nlumos)
623 CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
624 unoccupied_orbs(ispin), wf_g, wf_r, particles, nlumos, homo, ispin)
629 IF (dft_control%do_tddfpt_calculation)
THEN
630 ALLOCATE (dft_control%tddfpt_control%lumos_eigenvalues(nlumos, dft_control%nspins))
631 DO ispin = 1, dft_control%nspins
632 dft_control%tddfpt_control%lumos_eigenvalues(1:nlumos, ispin) = &
633 unoccupied_evals(ispin)%array(1:nlumos)
635 dft_control%tddfpt_control%lumos => unoccupied_orbs
639 ALLOCATE (lumo_localized(dft_control%nspins))
640 DO ispin = 1, dft_control%nspins
641 CALL cp_fm_create(lumo_localized(ispin), unoccupied_orbs(ispin)%matrix_struct)
642 CALL cp_fm_to_fm(unoccupied_orbs(ispin), lumo_localized(ispin))
644 CALL qs_loc_init(qs_env, qs_loc_env_lumo, localize_section, lumo_localized, do_homo, do_mo_cubes, &
645 evals=unoccupied_evals)
646 CALL qs_loc_env_init(qs_loc_env_lumo, qs_loc_env_lumo%localized_wfn_control, qs_env, &
647 loc_coeff=unoccupied_orbs)
649 lumo_localized, wf_r, wf_g, particles, &
650 unoccupied_orbs, unoccupied_evals, marked_states)
651 CALL loc_write_restart(qs_loc_env_lumo, loc_print_section, mos, homo_localized, do_homo, &
652 evals=unoccupied_evals)
653 lumo_ptr => lumo_localized
657 IF (has_homo .AND. has_lumo)
THEN
658 IF (output_unit > 0)
WRITE (output_unit, *)
" "
659 DO ispin = 1, dft_control%nspins
660 IF (.NOT. scf_control%smear%do_smear)
THEN
661 gap = homo_lumo(ispin, 2) - homo_lumo(ispin, 1)
662 IF (output_unit > 0)
WRITE (output_unit,
'(T2,A,F12.6)') &
663 "HOMO - LUMO gap [eV] :", gap*
evolt
669 IF (p_loc_mixed)
THEN
671 cpwarn(
"Localization not implemented for k-point calculations!!")
672 ELSEIF (dft_control%restricted)
THEN
673 IF (output_unit > 0)
WRITE (output_unit, *) &
674 " Unclear how we define MOs / localization in the restricted case... skipping"
677 ALLOCATE (mixed_orbs(dft_control%nspins))
678 ALLOCATE (mixed_evals(dft_control%nspins))
679 ALLOCATE (mixed_localized(dft_control%nspins))
680 DO ispin = 1, dft_control%nspins
681 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
682 eigenvalues=mo_eigenvalues)
683 mixed_orbs(ispin) = mo_coeff
684 mixed_evals(ispin)%array => mo_eigenvalues
685 CALL cp_fm_create(mixed_localized(ispin), mixed_orbs(ispin)%matrix_struct)
686 CALL cp_fm_to_fm(mixed_orbs(ispin), mixed_localized(ispin))
689 CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history)
692 total_zeff_corr = scf_env%sum_zeff_corr
693 ALLOCATE (qs_loc_env_mixed)
696 CALL qs_loc_init(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, do_homo, &
697 do_mo_cubes, mo_loc_history=mo_loc_history, tot_zeff_corr=total_zeff_corr, &
700 DO ispin = 1, dft_control%nspins
701 CALL cp_fm_get_info(mixed_localized(ispin), ncol_global=nchk_nmoloc)
705 wf_r, wf_g, particles, mixed_orbs, mixed_evals, marked_states)
708 IF (qs_loc_env_mixed%localized_wfn_control%use_history)
THEN
710 CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history)
717 DEALLOCATE (mixed_orbs)
718 DEALLOCATE (mixed_evals)
728 IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc)
THEN
729 CALL auxbas_pw_pool%give_back_pw(wf_r)
730 CALL auxbas_pw_pool%give_back_pw(wf_g)
734 IF (.NOT. do_kpoints)
THEN
737 DEALLOCATE (qs_loc_env_homo)
741 DEALLOCATE (qs_loc_env_lumo)
743 IF (p_loc_mixed)
THEN
745 DEALLOCATE (qs_loc_env_mixed)
753 CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env)
754 CALL wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, &
755 output_unit, unoccupied_orbs=lumo_ptr, scf_env=scf_env, &
756 matrix_s=matrix_s, marked_states=marked_states)
760 IF (
ASSOCIATED(marked_states))
THEN
761 DEALLOCATE (marked_states)
765 IF (.NOT. do_kpoints)
THEN
766 IF (compute_lumos)
THEN
767 DO ispin = 1, dft_control%nspins
768 DEALLOCATE (unoccupied_evals(ispin)%array)
769 IF (.NOT. dft_control%do_tddfpt_calculation)
THEN
773 DEALLOCATE (unoccupied_evals)
774 IF (.NOT. dft_control%do_tddfpt_calculation)
THEN
775 DEALLOCATE (unoccupied_orbs)
783 cpwarn(
"STM not implemented for k-point calculations!")
785 NULLIFY (unoccupied_orbs_stm, unoccupied_evals_stm)
786 IF (nlumo_stm > 0)
THEN
787 ALLOCATE (unoccupied_orbs_stm(dft_control%nspins))
788 ALLOCATE (unoccupied_evals_stm(dft_control%nspins))
789 CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs_stm, unoccupied_evals_stm, &
793 CALL th_stm_image(qs_env, stm_section, particles, unoccupied_orbs_stm, &
794 unoccupied_evals_stm)
796 IF (nlumo_stm > 0)
THEN
797 DO ispin = 1, dft_control%nspins
798 DEALLOCATE (unoccupied_evals_stm(ispin)%array)
800 DEALLOCATE (unoccupied_evals_stm)
807 CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit)
810 CALL qs_scf_post_efg(input, logger, qs_env)
813 CALL qs_scf_post_et(input, qs_env, dft_control)
816 CALL qs_scf_post_epr(input, logger, qs_env)
819 CALL qs_scf_post_molopt(input, logger, qs_env)
822 CALL qs_scf_post_elf(input, logger, qs_env)
829 DO ispin = 1, dft_control%nspins
830 CALL dbcsr_add(rho_ao(ispin, 1)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
836 CALL timestop(handle)
849 SUBROUTINE make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
853 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(INOUT) :: unoccupied_orbs
855 INTEGER,
INTENT(IN) :: nlumo
856 INTEGER,
INTENT(OUT) :: nlumos
858 CHARACTER(len=*),
PARAMETER :: routinen =
'make_lumo_gpw'
860 INTEGER :: handle, homo, ispin, n, nao, nmo, &
867 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_rmpv, matrix_s
874 CALL timeset(routinen, handle)
876 NULLIFY (mos, ks_rmpv, scf_control, dft_control, admm_env, para_env, blacs_env)
880 scf_control=scf_control, &
881 dft_control=dft_control, &
890 DO ispin = 1, dft_control%nspins
891 NULLIFY (unoccupied_evals(ispin)%array)
893 IF (output_unit > 0)
WRITE (output_unit, *)
" "
894 IF (output_unit > 0)
WRITE (output_unit, *)
" Lowest Eigenvalues of the unoccupied subspace spin ", ispin
895 IF (output_unit > 0)
WRITE (output_unit, fmt=
'(1X,A)')
"-----------------------------------------------------"
896 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo)
898 nlumos = max(1, min(nlumo, nao - nmo))
899 IF (nlumo == -1) nlumos = nao - nmo
900 ALLOCATE (unoccupied_evals(ispin)%array(nlumos))
902 nrow_global=n, ncol_global=nlumos)
903 CALL cp_fm_create(unoccupied_orbs(ispin), fm_struct_tmp, name=
"lumos")
908 NULLIFY (local_preconditioner)
909 IF (
ASSOCIATED(scf_env%ot_preconditioner))
THEN
910 local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
913 NULLIFY (local_preconditioner)
918 IF (dft_control%do_admm)
THEN
922 CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
923 matrix_c_fm=unoccupied_orbs(ispin), &
924 matrix_orthogonal_space_fm=mo_coeff, &
925 eps_gradient=scf_control%eps_lumos, &
927 iter_max=scf_control%max_iter_lumos, &
928 size_ortho_space=nmo)
931 unoccupied_evals(ispin)%array, scr=output_unit, &
932 ionode=output_unit > 0)
935 IF (dft_control%do_admm)
THEN
941 CALL timestop(handle)
950 SUBROUTINE qs_scf_post_charges(input, logger, qs_env)
955 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_charges'
957 INTEGER :: handle, print_level, unit_nr
958 LOGICAL :: do_kpoints, print_it
961 CALL timeset(routinen, handle)
963 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
971 log_filename=.false.)
974 IF (print_it) print_level = 2
976 IF (print_it) print_level = 3
978 cpwarn(
"Lowdin charges not implemented for k-point calculations!")
991 unit_nr =
cp_print_key_unit_nr(logger, input,
"PROPERTIES%FIT_CHARGE", extension=
".Fitcharge", &
992 log_filename=.false.)
994 CALL get_ddapc(qs_env, .false., density_fit_section, iwc=unit_nr)
998 CALL timestop(handle)
1000 END SUBROUTINE qs_scf_post_charges
1016 SUBROUTINE qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, &
1017 mo_coeff, wf_g, wf_r, particles, homo, ispin)
1026 INTEGER,
INTENT(IN) :: homo, ispin
1028 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_occ_cubes'
1030 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1031 INTEGER :: handle, i, ir, ivector, n_rep, nhomo, &
1033 INTEGER,
DIMENSION(:),
POINTER ::
list, list_index
1034 LOGICAL :: append_cube, mpi_io
1039 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1041 CALL timeset(routinen, handle)
1043 NULLIFY (list_index)
1049 my_pos_cube =
"REWIND"
1050 IF (append_cube)
THEN
1051 my_pos_cube =
"APPEND"
1060 IF (
ASSOCIATED(
list))
THEN
1062 DO i = 1,
SIZE(
list)
1063 list_index(i + nlist) =
list(i)
1065 nlist = nlist +
SIZE(
list)
1070 IF (nhomo == -1) nhomo = homo
1071 nlist = homo - max(1, homo - nhomo + 1) + 1
1072 ALLOCATE (list_index(nlist))
1074 list_index(i) = max(1, homo - nhomo + 1) + i - 1
1078 ivector = list_index(i)
1080 atomic_kind_set=atomic_kind_set, &
1081 qs_kind_set=qs_kind_set, &
1083 particle_set=particle_set, &
1086 cell, dft_control, particle_set, pw_env)
1087 WRITE (filename,
'(a4,I5.5,a1,I1.1)')
"WFN_", ivector,
"_", ispin
1090 middle_name=trim(filename), file_position=my_pos_cube, log_filename=.false., &
1092 WRITE (title, *)
"WAVEFUNCTION ", ivector,
" spin ", ispin,
" i.e. HOMO - ", ivector - homo
1093 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, &
1097 IF (
ASSOCIATED(list_index))
DEALLOCATE (list_index)
1100 CALL timestop(handle)
1102 END SUBROUTINE qs_scf_post_occ_cubes
1120 SUBROUTINE qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
1121 unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, lumo)
1127 TYPE(
cp_fm_type),
INTENT(IN) :: unoccupied_orbs
1131 INTEGER,
INTENT(IN) :: nlumos, homo, ispin
1132 INTEGER,
INTENT(IN),
OPTIONAL :: lumo
1134 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_unocc_cubes'
1136 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1137 INTEGER :: handle, ifirst, index_mo, ivector, &
1139 LOGICAL :: append_cube, mpi_io
1144 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1146 CALL timeset(routinen, handle)
1150 NULLIFY (qs_kind_set, particle_set, pw_env, cell)
1152 my_pos_cube =
"REWIND"
1153 IF (append_cube)
THEN
1154 my_pos_cube =
"APPEND"
1157 IF (
PRESENT(lumo)) ifirst = lumo
1158 DO ivector = ifirst, ifirst + nlumos - 1
1160 atomic_kind_set=atomic_kind_set, &
1161 qs_kind_set=qs_kind_set, &
1163 particle_set=particle_set, &
1166 qs_kind_set, cell, dft_control, particle_set, pw_env)
1168 IF (ifirst == 1)
THEN
1169 index_mo = homo + ivector
1173 WRITE (filename,
'(a4,I5.5,a1,I1.1)')
"WFN_", index_mo,
"_", ispin
1177 middle_name=trim(filename), file_position=my_pos_cube, log_filename=.false., &
1179 WRITE (title, *)
"WAVEFUNCTION ", index_mo,
" spin ", ispin,
" i.e. LUMO + ", ifirst + ivector - 2
1180 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, &
1186 CALL timestop(handle)
1188 END SUBROUTINE qs_scf_post_unocc_cubes
1201 INTEGER,
INTENT(IN) :: output_unit
1203 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_moments'
1205 CHARACTER(LEN=default_path_length) :: filename
1206 INTEGER :: handle, maxmom, reference, unit_nr
1207 LOGICAL :: com_nl, do_kpoints, magnetic, periodic, &
1208 second_ref_point, vel_reprs
1209 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ref_point
1212 CALL timeset(routinen, handle)
1215 subsection_name=
"DFT%PRINT%MOMENTS")
1220 keyword_name=
"DFT%PRINT%MOMENTS%MAX_MOMENT")
1222 keyword_name=
"DFT%PRINT%MOMENTS%PERIODIC")
1224 keyword_name=
"DFT%PRINT%MOMENTS%REFERENCE")
1226 keyword_name=
"DFT%PRINT%MOMENTS%MAGNETIC")
1228 keyword_name=
"DFT%PRINT%MOMENTS%VEL_REPRS")
1230 keyword_name=
"DFT%PRINT%MOMENTS%COM_NL")
1232 keyword_name=
"DFT%PRINT%MOMENTS%SECOND_REFERENCE_POINT")
1237 print_key_path=
"DFT%PRINT%MOMENTS", extension=
".dat", &
1238 middle_name=
"moments", log_filename=.false.)
1240 IF (output_unit > 0)
THEN
1241 IF (unit_nr /= output_unit)
THEN
1242 INQUIRE (unit=unit_nr, name=filename)
1243 WRITE (unit=output_unit, fmt=
"(/,T2,A,2(/,T3,A),/)") &
1244 "MOMENTS",
"The electric/magnetic moments are written to file:", &
1247 WRITE (unit=output_unit, fmt=
"(/,T2,A)")
"ELECTRIC/MAGNETIC MOMENTS"
1251 CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
1253 IF (do_kpoints)
THEN
1254 cpwarn(
"Moments not implemented for k-point calculations!")
1259 CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr, vel_reprs, com_nl)
1264 basis_section=input, print_key_path=
"DFT%PRINT%MOMENTS")
1266 IF (second_ref_point)
THEN
1268 keyword_name=
"DFT%PRINT%MOMENTS%REFERENCE_2")
1273 print_key_path=
"DFT%PRINT%MOMENTS", extension=
".dat", &
1274 middle_name=
"moments_refpoint_2", log_filename=.false.)
1276 IF (output_unit > 0)
THEN
1277 IF (unit_nr /= output_unit)
THEN
1278 INQUIRE (unit=unit_nr, name=filename)
1279 WRITE (unit=output_unit, fmt=
"(/,T2,A,2(/,T3,A),/)") &
1280 "MOMENTS",
"The electric/magnetic moments for the second reference point are written to file:", &
1283 WRITE (unit=output_unit, fmt=
"(/,T2,A)")
"ELECTRIC/MAGNETIC MOMENTS"
1286 IF (do_kpoints)
THEN
1287 cpwarn(
"Moments not implemented for k-point calculations!")
1292 CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr, vel_reprs, com_nl)
1296 basis_section=input, print_key_path=
"DFT%PRINT%MOMENTS")
1301 CALL timestop(handle)
1313 SUBROUTINE qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit)
1318 INTEGER,
INTENT(IN) :: output_unit
1320 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_post_xray'
1322 CHARACTER(LEN=default_path_length) :: filename
1323 INTEGER :: handle, unit_nr
1324 REAL(kind=
dp) :: q_max
1327 CALL timeset(routinen, handle)
1330 subsection_name=
"DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
1334 keyword_name=
"PRINT%XRAY_DIFFRACTION_SPECTRUM%Q_MAX")
1336 basis_section=input, &
1337 print_key_path=
"DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM", &
1339 middle_name=
"xrd", &
1340 log_filename=.false.)
1341 IF (output_unit > 0)
THEN
1342 INQUIRE (unit=unit_nr, name=filename)
1343 WRITE (unit=output_unit, fmt=
"(/,/,T2,A)") &
1344 "X-RAY DIFFRACTION SPECTRUM"
1345 IF (unit_nr /= output_unit)
THEN
1346 WRITE (unit=output_unit, fmt=
"(/,T3,A,/,/,T3,A,/)") &
1347 "The coherent X-ray diffraction spectrum is written to the file:", &
1352 unit_number=unit_nr, &
1356 basis_section=input, &
1357 print_key_path=
"DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
1360 CALL timestop(handle)
1362 END SUBROUTINE qs_scf_post_xray
1370 SUBROUTINE qs_scf_post_efg(input, logger, qs_env)
1375 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_efg'
1380 CALL timeset(routinen, handle)
1383 subsection_name=
"DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
1389 CALL timestop(handle)
1391 END SUBROUTINE qs_scf_post_efg
1399 SUBROUTINE qs_scf_post_et(input, qs_env, dft_control)
1404 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_et'
1406 INTEGER :: handle, ispin
1408 TYPE(
cp_fm_type),
DIMENSION(:),
POINTER :: my_mos
1411 CALL timeset(routinen, handle)
1417 IF (qs_env%et_coupling%first_run)
THEN
1419 ALLOCATE (my_mos(dft_control%nspins))
1420 ALLOCATE (qs_env%et_coupling%et_mo_coeff(dft_control%nspins))
1421 DO ispin = 1, dft_control%nspins
1423 matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
1424 name=
"FIRST_RUN_COEFF"//trim(adjustl(
cp_to_string(ispin)))//
"MATRIX")
1433 CALL timestop(handle)
1435 END SUBROUTINE qs_scf_post_et
1446 SUBROUTINE qs_scf_post_elf(input, logger, qs_env)
1451 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_elf'
1453 CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube, &
1455 INTEGER :: handle, ispin, output_unit, unit_nr
1456 LOGICAL :: append_cube, gapw, mpi_io
1457 REAL(
dp) :: rho_cutoff
1467 CALL timeset(routinen, handle)
1474 NULLIFY (dft_control, pw_env, auxbas_pw_pool, pw_pools, particles, subsys)
1475 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, subsys=subsys)
1478 gapw = dft_control%qs_control%gapw
1479 IF (.NOT. gapw)
THEN
1481 ALLOCATE (elf_r(dft_control%nspins))
1482 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1484 DO ispin = 1, dft_control%nspins
1485 CALL auxbas_pw_pool%create_pw(elf_r(ispin))
1489 IF (output_unit > 0)
THEN
1490 WRITE (unit=output_unit, fmt=
"(/,T15,A,/)") &
1491 " ----- ELF is computed on the real space grid -----"
1498 my_pos_cube =
"REWIND"
1499 IF (append_cube)
THEN
1500 my_pos_cube =
"APPEND"
1503 DO ispin = 1, dft_control%nspins
1504 WRITE (filename,
'(a5,I1.1)')
"ELF_S", ispin
1505 WRITE (title, *)
"ELF spin ", ispin
1508 middle_name=trim(filename), file_position=my_pos_cube, log_filename=.false., &
1509 mpi_io=mpi_io, fout=mpi_filename)
1510 IF (output_unit > 0)
THEN
1511 IF (.NOT. mpi_io)
THEN
1512 INQUIRE (unit=unit_nr, name=filename)
1514 filename = mpi_filename
1516 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
1517 "ELF is written in cube file format to the file:", &
1521 CALL cp_pw_to_cube(elf_r(ispin), unit_nr, title, particles=particles, &
1525 CALL auxbas_pw_pool%give_back_pw(elf_r(ispin))
1533 cpwarn(
"ELF not implemented for GAPW calculations!!")
1539 CALL timestop(handle)
1541 END SUBROUTINE qs_scf_post_elf
1553 SUBROUTINE qs_scf_post_molopt(input, logger, qs_env)
1558 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_molopt'
1560 INTEGER :: handle, nao, unit_nr
1561 REAL(kind=
dp) :: s_cond_number
1562 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues
1566 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
1571 CALL timeset(routinen, handle)
1574 subsection_name=
"DFT%PRINT%BASIS_MOLOPT_QUANTITIES")
1578 CALL get_qs_env(qs_env, energy=energy, matrix_s=matrix_s, mos=mos)
1581 CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, nao=nao)
1583 nrow_global=nao, ncol_global=nao, &
1584 template_fmstruct=mo_coeff%matrix_struct)
1585 CALL cp_fm_create(fm_s, matrix_struct=ao_ao_fmstruct, &
1587 CALL cp_fm_create(fm_work, matrix_struct=ao_ao_fmstruct, &
1590 ALLOCATE (eigenvalues(nao))
1598 s_cond_number = maxval(abs(eigenvalues))/max(minval(abs(eigenvalues)), epsilon(0.0_dp))
1601 extension=
".molopt")
1603 IF (unit_nr > 0)
THEN
1606 WRITE (unit_nr,
'(T2,A28,2A25)')
"",
"Tot. Ener.",
"S Cond. Numb."
1607 WRITE (unit_nr,
'(T2,A28,2E25.17)')
"BASIS_MOLOPT_QUANTITIES", energy%total, s_cond_number
1611 "DFT%PRINT%BASIS_MOLOPT_QUANTITIES")
1615 CALL timestop(handle)
1617 END SUBROUTINE qs_scf_post_molopt
1625 SUBROUTINE qs_scf_post_epr(input, logger, qs_env)
1630 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_epr'
1635 CALL timeset(routinen, handle)
1638 subsection_name=
"DFT%PRINT%HYPERFINE_COUPLING_TENSOR")
1644 CALL timestop(handle)
1646 END SUBROUTINE qs_scf_post_epr
1655 SUBROUTINE write_available_results(qs_env, scf_env)
1659 CHARACTER(len=*),
PARAMETER :: routinen =
'write_available_results'
1663 CALL timeset(routinen, handle)
1671 CALL timestop(handle)
1673 END SUBROUTINE write_available_results
1686 CHARACTER(len=*),
PARAMETER :: routinen =
'write_mo_dependent_results'
1688 INTEGER :: handle, homo, ispin, nmo, output_unit
1689 LOGICAL :: all_equal, do_kpoints
1690 REAL(kind=
dp) :: maxocc, s_square, s_square_ideal, &
1692 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_eigenvalues, occupation_numbers
1698 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_rmpv, matrix_s
1699 TYPE(dbcsr_type),
POINTER :: mo_coeff_deriv
1712 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1718 CALL timeset(routinen, handle)
1720 NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, mo_coeff, &
1721 mo_coeff_deriv, mo_eigenvalues, mos, atomic_kind_set, qs_kind_set, &
1722 particle_set, rho, ks_rmpv, matrix_s, scf_control, dft_section, &
1723 molecule_set, input, particles, subsys, rho_r)
1728 cpassert(
ASSOCIATED(qs_env))
1730 dft_control=dft_control, &
1731 molecule_set=molecule_set, &
1732 atomic_kind_set=atomic_kind_set, &
1733 particle_set=particle_set, &
1734 qs_kind_set=qs_kind_set, &
1735 admm_env=admm_env, &
1736 scf_control=scf_control, &
1745 CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
1749 IF (.NOT. qs_env%run_rtp)
THEN
1751 IF (.NOT. do_kpoints)
THEN
1752 CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
1758 dft_section,
"PRINT%CHARGEMOL"), &
1767 IF (do_kpoints)
THEN
1768 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
1779 IF (do_kpoints)
THEN
1780 cpwarn(
"Projected density of states (pDOS) is not implemented for k points")
1785 DO ispin = 1, dft_control%nspins
1787 IF (dft_control%do_admm)
THEN
1790 IF (
PRESENT(scf_env))
THEN
1792 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
1793 eigenvalues=mo_eigenvalues)
1794 IF (
ASSOCIATED(qs_env%mo_derivs))
THEN
1795 mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
1797 mo_coeff_deriv => null()
1800 do_rotation=.true., &
1801 co_rotate_dbcsr=mo_coeff_deriv)
1805 IF (dft_control%nspins == 2)
THEN
1807 qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin)
1810 qs_kind_set, particle_set, qs_env, dft_section)
1813 IF (dft_control%do_admm)
THEN
1822 IF (dft_control%nspins == 2)
THEN
1824 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1825 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1827 CALL auxbas_pw_pool%create_pw(wf_r)
1829 CALL pw_axpy(rho_r(2), wf_r, alpha=-1._dp)
1831 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
'(/,(T3,A,T61,F20.10))') &
1832 "Integrated absolute spin density : ", total_abs_spin_dens
1833 CALL auxbas_pw_pool%give_back_pw(wf_r)
1839 IF (do_kpoints)
THEN
1840 cpwarn(
"Spin contamination estimate not implemented for k-points.")
1843 DO ispin = 1, dft_control%nspins
1845 occupation_numbers=occupation_numbers, &
1850 all_equal = all_equal .AND. &
1851 (all(occupation_numbers(1:homo) == maxocc) .AND. &
1852 all(occupation_numbers(homo + 1:nmo) == 0.0_dp))
1855 IF (.NOT. all_equal)
THEN
1856 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"(T3,A)") &
1857 "WARNING: S**2 computation does not yet treat fractional occupied orbitals"
1860 matrix_s=matrix_s, &
1863 s_square_ideal=s_square_ideal)
1864 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
'(T3,A,T51,2F15.6)') &
1865 "Ideal and single determinant S**2 : ", s_square_ideal, s_square
1866 energy%s_square = s_square
1871 CALL timestop(handle)
1883 CHARACTER(len=*),
PARAMETER :: routinen =
'write_mo_free_results'
1884 CHARACTER(len=1),
DIMENSION(3),
PARAMETER :: cdir = (/
"x",
"y",
"z"/)
1886 CHARACTER(LEN=2) :: element_symbol
1887 CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube, &
1889 CHARACTER(LEN=default_string_length) :: name, print_density
1890 INTEGER :: after, handle, i, iat, id, ikind, img, iso, ispin, iw, l, n_rep_hf, natom, nd(3), &
1891 ngto, niso, nkind, np, nr, output_unit, print_level, should_print_bqb, should_print_voro, &
1892 unit_nr, unit_nr_voro
1893 LOGICAL :: append_cube, append_voro, do_hfx, do_kpoints, mpi_io, omit_headers, print_it, &
1894 rho_r_valid, voro_print_txt, write_ks, write_xc, xrd_interface
1895 REAL(kind=
dp) :: norm_factor, q_max, rho_hard, rho_soft, &
1896 rho_total, rho_total_rspace, udvol, &
1898 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: bfun
1899 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: aedens, ccdens, ppdens
1900 REAL(kind=
dp),
DIMENSION(3) :: dr
1901 REAL(kind=
dp),
DIMENSION(:),
POINTER :: my_q0
1906 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_hr
1907 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_rmpv, matrix_vxc, rho_ao
1922 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1930 print_key, print_key_bqb, &
1931 print_key_voro, xc_section
1933 CALL timeset(routinen, handle)
1934 NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, hfx_section, &
1935 atomic_kind_set, qs_kind_set, particle_set, rho, ks_rmpv, rho_ao, rho_r, &
1936 dft_section, xc_section, input, particles, subsys, matrix_vxc, v_hartree_rspace, &
1942 cpassert(
ASSOCIATED(qs_env))
1944 atomic_kind_set=atomic_kind_set, &
1945 qs_kind_set=qs_kind_set, &
1946 particle_set=particle_set, &
1948 para_env=para_env, &
1949 dft_control=dft_control, &
1951 do_kpoints=do_kpoints, &
1961 "DFT%PRINT%TOT_DENSITY_CUBE"),
cp_p_file))
THEN
1962 NULLIFY (rho_core, rho0_s_gs)
1964 my_pos_cube =
"REWIND"
1965 IF (append_cube)
THEN
1966 my_pos_cube =
"APPEND"
1969 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho_core=rho_core, &
1970 rho0_s_gs=rho0_s_gs)
1971 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1973 CALL auxbas_pw_pool%create_pw(wf_r)
1974 IF (dft_control%qs_control%gapw)
THEN
1975 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw)
THEN
1976 CALL pw_axpy(rho_core, rho0_s_gs)
1978 CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
1985 DO ispin = 1, dft_control%nspins
1986 CALL pw_axpy(rho_r(ispin), wf_r)
1988 filename =
"TOTAL_DENSITY"
1991 extension=
".cube", middle_name=trim(filename), file_position=my_pos_cube, &
1992 log_filename=.false., mpi_io=mpi_io)
1994 particles=particles, &
1995 stride=
section_get_ivals(dft_section,
"PRINT%TOT_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
1997 "DFT%PRINT%TOT_DENSITY_CUBE", mpi_io=mpi_io)
1998 CALL auxbas_pw_pool%give_back_pw(wf_r)
2003 "DFT%PRINT%E_DENSITY_CUBE"),
cp_p_file))
THEN
2005 keyword_name=
"PRINT%E_DENSITY_CUBE%DENSITY_INCLUDE", &
2006 c_val=print_density)
2007 print_density = trim(print_density)
2009 my_pos_cube =
"REWIND"
2010 IF (append_cube)
THEN
2011 my_pos_cube =
"APPEND"
2015 xrd_interface =
section_get_lval(input,
"DFT%PRINT%E_DENSITY_CUBE%XRD_INTERFACE")
2016 IF (xrd_interface)
THEN
2018 IF (dft_control%qs_control%gapw) print_density =
"SOFT_DENSITY"
2020 filename =
"ELECTRON_DENSITY"
2022 extension=
".xrd", middle_name=trim(filename), &
2023 file_position=my_pos_cube, log_filename=.false.)
2025 IF (output_unit > 0)
THEN
2026 INQUIRE (unit=unit_nr, name=filename)
2027 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2028 "The electron density (atomic part) is written to the file:", &
2033 nkind =
SIZE(atomic_kind_set)
2034 IF (unit_nr > 0)
THEN
2035 WRITE (unit_nr, *)
"Atomic (core) densities"
2036 WRITE (unit_nr, *)
"Unit cell"
2037 WRITE (unit_nr, fmt=
"(3F20.12)") cell%hmat(1, 1), cell%hmat(1, 2), cell%hmat(1, 3)
2038 WRITE (unit_nr, fmt=
"(3F20.12)") cell%hmat(2, 1), cell%hmat(2, 2), cell%hmat(2, 3)
2039 WRITE (unit_nr, fmt=
"(3F20.12)") cell%hmat(3, 1), cell%hmat(3, 2), cell%hmat(3, 3)
2040 WRITE (unit_nr, *)
"Atomic types"
2041 WRITE (unit_nr, *) nkind
2044 ALLOCATE (ppdens(ngto, 2, nkind), aedens(ngto, 2, nkind), ccdens(ngto, 2, nkind))
2046 atomic_kind => atomic_kind_set(ikind)
2047 qs_kind => qs_kind_set(ikind)
2048 CALL get_atomic_kind(atomic_kind, name=name, element_symbol=element_symbol)
2050 iunit=output_unit, confine=.true.)
2052 iunit=output_unit, allelectron=.true., confine=.true.)
2053 ccdens(:, 1, ikind) = aedens(:, 1, ikind)
2054 ccdens(:, 2, ikind) = 0._dp
2055 CALL project_function_a(ccdens(1:ngto, 2, ikind), ccdens(1:ngto, 1, ikind), &
2056 ppdens(1:ngto, 2, ikind), ppdens(1:ngto, 1, ikind), 0)
2057 ccdens(:, 2, ikind) = aedens(:, 2, ikind) - ccdens(:, 2, ikind)
2058 IF (unit_nr > 0)
THEN
2059 WRITE (unit_nr, fmt=
"(I6,A10,A20)") ikind, trim(element_symbol), trim(name)
2060 WRITE (unit_nr, fmt=
"(I6)") ngto
2061 WRITE (unit_nr, *)
" Total density"
2062 WRITE (unit_nr, fmt=
"(2G24.12)") (aedens(i, 1, ikind), aedens(i, 2, ikind), i=1, ngto)
2063 WRITE (unit_nr, *)
" Core density"
2064 WRITE (unit_nr, fmt=
"(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
2066 NULLIFY (atomic_kind)
2069 IF (dft_control%qs_control%gapw)
THEN
2070 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
2072 IF (unit_nr > 0)
THEN
2073 WRITE (unit_nr, *)
"Coordinates and GAPW density"
2075 np = particles%n_els
2077 CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
2078 CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
2079 rho_atom => rho_atom_set(iat)
2080 IF (
ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef))
THEN
2081 nr =
SIZE(rho_atom%rho_rad_h(1)%r_coef, 1)
2082 niso =
SIZE(rho_atom%rho_rad_h(1)%r_coef, 2)
2087 CALL para_env%sum(nr)
2088 CALL para_env%sum(niso)
2090 ALLOCATE (bfun(nr, niso))
2092 DO ispin = 1, dft_control%nspins
2093 IF (
ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef))
THEN
2094 bfun(:, :) = bfun + rho_atom%rho_rad_h(ispin)%r_coef - rho_atom%rho_rad_s(ispin)%r_coef
2097 CALL para_env%sum(bfun)
2098 ccdens(:, 1, ikind) = ppdens(:, 1, ikind)
2099 ccdens(:, 2, ikind) = 0._dp
2100 IF (unit_nr > 0)
THEN
2101 WRITE (unit_nr,
'(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
2105 CALL project_function_b(ccdens(:, 2, ikind), ccdens(:, 1, ikind), bfun(:, iso), grid_atom, l)
2106 IF (unit_nr > 0)
THEN
2107 WRITE (unit_nr, fmt=
"(3I6)") iso, l, ngto
2108 WRITE (unit_nr, fmt=
"(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
2114 IF (unit_nr > 0)
THEN
2115 WRITE (unit_nr, *)
"Coordinates"
2116 np = particles%n_els
2118 CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
2119 WRITE (unit_nr,
'(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
2124 DEALLOCATE (ppdens, aedens, ccdens)
2127 "DFT%PRINT%E_DENSITY_CUBE")
2130 IF (dft_control%qs_control%gapw .AND. print_density ==
"TOTAL_DENSITY")
THEN
2132 cpassert(.NOT. do_kpoints)
2137 auxbas_pw_pool=auxbas_pw_pool, &
2139 CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
2141 CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
2146 q_max = sqrt(sum((
pi/dr(:))**2))
2148 auxbas_pw_pool=auxbas_pw_pool, &
2149 rhotot_elec_gspace=rho_elec_gspace, &
2151 rho_hard=rho_hard, &
2153 rho_total = rho_hard + rho_soft
2158 CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
2160 CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
2162 filename =
"TOTAL_ELECTRON_DENSITY"
2165 extension=
".cube", middle_name=trim(filename), &
2166 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2168 IF (output_unit > 0)
THEN
2169 IF (.NOT. mpi_io)
THEN
2170 INQUIRE (unit=unit_nr, name=filename)
2172 filename = mpi_filename
2174 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2175 "The total electron density is written in cube file format to the file:", &
2177 WRITE (unit=output_unit, fmt=
"(/,(T2,A,F20.10))") &
2178 "q(max) [1/Angstrom] :", q_max/
angstrom, &
2179 "Soft electronic charge (G-space) :", rho_soft, &
2180 "Hard electronic charge (G-space) :", rho_hard, &
2181 "Total electronic charge (G-space):", rho_total, &
2182 "Total electronic charge (R-space):", rho_total_rspace
2184 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"TOTAL ELECTRON DENSITY", &
2185 particles=particles, &
2186 stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
2188 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2190 IF (dft_control%nspins > 1)
THEN
2194 auxbas_pw_pool=auxbas_pw_pool, &
2195 rhotot_elec_gspace=rho_elec_gspace, &
2197 rho_hard=rho_hard, &
2198 rho_soft=rho_soft, &
2200 rho_total = rho_hard + rho_soft
2204 CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
2206 CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
2208 filename =
"TOTAL_SPIN_DENSITY"
2211 extension=
".cube", middle_name=trim(filename), &
2212 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2214 IF (output_unit > 0)
THEN
2215 IF (.NOT. mpi_io)
THEN
2216 INQUIRE (unit=unit_nr, name=filename)
2218 filename = mpi_filename
2220 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2221 "The total spin density is written in cube file format to the file:", &
2223 WRITE (unit=output_unit, fmt=
"(/,(T2,A,F20.10))") &
2224 "q(max) [1/Angstrom] :", q_max/
angstrom, &
2225 "Soft part of the spin density (G-space):", rho_soft, &
2226 "Hard part of the spin density (G-space):", rho_hard, &
2227 "Total spin density (G-space) :", rho_total, &
2228 "Total spin density (R-space) :", rho_total_rspace
2230 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"TOTAL SPIN DENSITY", &
2231 particles=particles, &
2232 stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
2234 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2236 CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
2237 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
2239 ELSE IF (print_density ==
"SOFT_DENSITY" .OR. .NOT. dft_control%qs_control%gapw)
THEN
2240 IF (dft_control%nspins > 1)
THEN
2244 auxbas_pw_pool=auxbas_pw_pool, &
2246 CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
2247 CALL pw_copy(rho_r(1), rho_elec_rspace)
2248 CALL pw_axpy(rho_r(2), rho_elec_rspace)
2249 filename =
"ELECTRON_DENSITY"
2252 extension=
".cube", middle_name=trim(filename), &
2253 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2255 IF (output_unit > 0)
THEN
2256 IF (.NOT. mpi_io)
THEN
2257 INQUIRE (unit=unit_nr, name=filename)
2259 filename = mpi_filename
2261 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2262 "The sum of alpha and beta density is written in cube file format to the file:", &
2265 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"SUM OF ALPHA AND BETA DENSITY", &
2266 particles=particles, stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), &
2269 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2270 CALL pw_copy(rho_r(1), rho_elec_rspace)
2271 CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
2272 filename =
"SPIN_DENSITY"
2275 extension=
".cube", middle_name=trim(filename), &
2276 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2278 IF (output_unit > 0)
THEN
2279 IF (.NOT. mpi_io)
THEN
2280 INQUIRE (unit=unit_nr, name=filename)
2282 filename = mpi_filename
2284 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2285 "The spin density is written in cube file format to the file:", &
2288 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"SPIN DENSITY", &
2289 particles=particles, &
2290 stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
2292 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2293 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
2295 filename =
"ELECTRON_DENSITY"
2298 extension=
".cube", middle_name=trim(filename), &
2299 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2301 IF (output_unit > 0)
THEN
2302 IF (.NOT. mpi_io)
THEN
2303 INQUIRE (unit=unit_nr, name=filename)
2305 filename = mpi_filename
2307 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2308 "The electron density is written in cube file format to the file:", &
2312 particles=particles, &
2313 stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
2315 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2318 ELSE IF (dft_control%qs_control%gapw .AND. print_density ==
"TOTAL_HARD_APPROX")
THEN
2319 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole, natom=natom)
2320 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
2321 CALL auxbas_pw_pool%create_pw(rho_elec_rspace)
2324 ALLOCATE (my_q0(natom))
2328 norm_factor = sqrt((rho0_mpole%zet0_h/
pi)**3)
2332 my_q0(iat) = sum(rho0_mpole%mp_rho(iat)%Q0(1:dft_control%nspins))*norm_factor
2336 CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
2340 DO ispin = 1, dft_control%nspins
2341 CALL pw_axpy(rho_r(ispin), rho_elec_rspace)
2345 rho_total_rspace = rho_soft + rho_hard
2347 filename =
"ELECTRON_DENSITY"
2350 extension=
".cube", middle_name=trim(filename), &
2351 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2353 IF (output_unit > 0)
THEN
2354 IF (.NOT. mpi_io)
THEN
2355 INQUIRE (unit=unit_nr, name=filename)
2357 filename = mpi_filename
2359 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2360 "The electron density is written in cube file format to the file:", &
2362 WRITE (unit=output_unit, fmt=
"(/,(T2,A,F20.10))") &
2363 "Soft electronic charge (R-space) :", rho_soft, &
2364 "Hard electronic charge (R-space) :", rho_hard, &
2365 "Total electronic charge (R-space):", rho_total_rspace
2367 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"ELECTRON DENSITY", &
2368 particles=particles, stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), &
2371 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2374 IF (dft_control%nspins > 1)
THEN
2376 my_q0(iat) = (rho0_mpole%mp_rho(iat)%Q0(1) - rho0_mpole%mp_rho(iat)%Q0(2))*norm_factor
2379 CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
2382 CALL pw_axpy(rho_r(1), rho_elec_rspace)
2383 CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
2387 rho_total_rspace = rho_soft + rho_hard
2389 filename =
"SPIN_DENSITY"
2392 extension=
".cube", middle_name=trim(filename), &
2393 file_position=my_pos_cube, log_filename=.false., mpi_io=mpi_io, &
2395 IF (output_unit > 0)
THEN
2396 IF (.NOT. mpi_io)
THEN
2397 INQUIRE (unit=unit_nr, name=filename)
2399 filename = mpi_filename
2401 WRITE (unit=output_unit, fmt=
"(/,T2,A,/,/,T2,A)") &
2402 "The spin density is written in cube file format to the file:", &
2404 WRITE (unit=output_unit, fmt=
"(/,(T2,A,F20.10))") &
2405 "Soft part of the spin density :", rho_soft, &
2406 "Hard part of the spin density :", rho_hard, &
2407 "Total spin density (R-space) :", rho_total_rspace
2409 CALL cp_pw_to_cube(rho_elec_rspace, unit_nr,
"SPIN DENSITY", &
2410 particles=particles, &
2411 stride=
section_get_ivals(dft_section,
"PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io)
2413 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io)
2415 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
2421 dft_section,
"PRINT%ENERGY_WINDOWS"),
cp_p_file) .AND. .NOT. do_kpoints)
THEN
2427 "DFT%PRINT%V_HARTREE_CUBE"),
cp_p_file))
THEN
2431 v_hartree_rspace=v_hartree_rspace)
2432 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2433 CALL auxbas_pw_pool%create_pw(aux_r)
2436 my_pos_cube =
"REWIND"
2437 IF (append_cube)
THEN
2438 my_pos_cube =
"APPEND"
2441 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
2444 extension=
".cube", middle_name=
"v_hartree", file_position=my_pos_cube, mpi_io=mpi_io)
2445 udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
2447 CALL pw_copy(v_hartree_rspace, aux_r)
2450 CALL cp_pw_to_cube(aux_r, unit_nr,
"HARTREE POTENTIAL", particles=particles, &
2452 "PRINT%V_HARTREE_CUBE%STRIDE"), mpi_io=mpi_io)
2454 "DFT%PRINT%V_HARTREE_CUBE", mpi_io=mpi_io)
2456 CALL auxbas_pw_pool%give_back_pw(aux_r)
2461 "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE"),
cp_p_file))
THEN
2462 IF (dft_control%apply_external_potential)
THEN
2463 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, vee=vee)
2464 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2465 CALL auxbas_pw_pool%create_pw(aux_r)
2467 append_cube =
section_get_lval(input,
"DFT%PRINT%EXTERNAL_POTENTIAL_CUBE%APPEND")
2468 my_pos_cube =
"REWIND"
2469 IF (append_cube)
THEN
2470 my_pos_cube =
"APPEND"
2475 extension=
".cube", middle_name=
"ext_pot", file_position=my_pos_cube, mpi_io=mpi_io)
2479 CALL cp_pw_to_cube(aux_r, unit_nr,
"EXTERNAL POTENTIAL", particles=particles, &
2481 "PRINT%EXTERNAL_POTENTIAL_CUBE%STRIDE"), mpi_io=mpi_io)
2483 "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", mpi_io=mpi_io)
2485 CALL auxbas_pw_pool%give_back_pw(aux_r)
2491 "DFT%PRINT%EFIELD_CUBE"),
cp_p_file))
THEN
2493 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
2494 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2495 CALL auxbas_pw_pool%create_pw(aux_r)
2496 CALL auxbas_pw_pool%create_pw(aux_g)
2499 my_pos_cube =
"REWIND"
2500 IF (append_cube)
THEN
2501 my_pos_cube =
"APPEND"
2503 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, &
2504 v_hartree_rspace=v_hartree_rspace)
2506 udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
2510 extension=
".cube", middle_name=
"efield_"//cdir(id), file_position=my_pos_cube, &
2521 unit_nr,
"ELECTRIC FIELD", particles=particles, &
2523 "PRINT%EFIELD_CUBE%STRIDE"), mpi_io=mpi_io)
2525 "DFT%PRINT%EFIELD_CUBE", mpi_io=mpi_io)
2528 CALL auxbas_pw_pool%give_back_pw(aux_r)
2529 CALL auxbas_pw_pool%give_back_pw(aux_g)
2533 CALL qs_scf_post_local_energy(input, logger, qs_env)
2536 CALL qs_scf_post_local_stress(input, logger, qs_env)
2539 CALL qs_scf_post_ps_implicit(input, logger, qs_env)
2547 "DFT%PRINT%AO_MATRICES/DENSITY"),
cp_p_file))
THEN
2552 after = min(max(after, 1), 16)
2553 DO ispin = 1, dft_control%nspins
2554 DO img = 1, dft_control%nimages
2556 para_env, output_unit=iw, omit_headers=omit_headers)
2560 "DFT%PRINT%AO_MATRICES/DENSITY")
2565 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"),
cp_p_file)
2567 "DFT%PRINT%AO_MATRICES/MATRIX_VXC"),
cp_p_file)
2569 IF (write_ks .OR. write_xc)
THEN
2570 IF (write_xc) qs_env%requires_matrix_vxc = .true.
2573 just_energy=.false.)
2574 IF (write_xc) qs_env%requires_matrix_vxc = .false.
2581 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=ks_rmpv)
2583 after = min(max(after, 1), 16)
2584 DO ispin = 1, dft_control%nspins
2585 DO img = 1, dft_control%nimages
2587 para_env, output_unit=iw, omit_headers=omit_headers)
2591 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
2596 IF (.NOT. dft_control%qs_control%pao)
THEN
2602 CALL write_adjacency_matrix(qs_env, input)
2606 CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc)
2607 cpassert(
ASSOCIATED(matrix_vxc))
2611 after = min(max(after, 1), 16)
2612 DO ispin = 1, dft_control%nspins
2613 DO img = 1, dft_control%nimages
2615 para_env, output_unit=iw, omit_headers=omit_headers)
2619 "DFT%PRINT%AO_MATRICES/MATRIX_VXC")
2624 "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR"),
cp_p_file))
THEN
2630 after = min(max(after, 1), 16)
2633 para_env, output_unit=iw, omit_headers=omit_headers)
2637 "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR")
2643 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%MULLIKEN", extension=
".mulliken", log_filename=.false.)
2646 IF (print_it) print_level = 2
2648 IF (print_it) print_level = 3
2659 CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
2660 IF (rho_r_valid)
THEN
2661 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%HIRSHFELD", extension=
".hirshfeld", log_filename=.false.)
2662 CALL hirshfeld_charges(qs_env, print_key, unit_nr)
2671 should_print_voro = 1
2673 should_print_voro = 0
2676 should_print_bqb = 1
2678 should_print_bqb = 0
2680 IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0))
THEN
2685 CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
2686 IF (rho_r_valid)
THEN
2688 IF (dft_control%nspins > 1)
THEN
2692 auxbas_pw_pool=auxbas_pw_pool, &
2696 CALL auxbas_pw_pool%create_pw(pw=mb_rho)
2697 CALL pw_copy(rho_r(1), mb_rho)
2698 CALL pw_axpy(rho_r(2), mb_rho)
2705 IF (should_print_voro /= 0)
THEN
2707 IF (voro_print_txt)
THEN
2709 my_pos_voro =
"REWIND"
2710 IF (append_voro)
THEN
2711 my_pos_voro =
"APPEND"
2713 unit_nr_voro =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%VORONOI", extension=
".voronoi", &
2714 file_position=my_pos_voro, log_filename=.false.)
2722 CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
2723 unit_nr_voro, qs_env, mb_rho)
2725 IF (dft_control%nspins > 1)
THEN
2726 CALL auxbas_pw_pool%give_back_pw(mb_rho)
2730 IF (unit_nr_voro > 0)
THEN
2740 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%MAO_ANALYSIS", extension=
".mao", log_filename=.false.)
2748 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%MINBAS_ANALYSIS", extension=
".mao", log_filename=.false.)
2756 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%IAO_ANALYSIS", extension=
".iao", log_filename=.false.)
2758 IF (iao_env%do_iao)
THEN
2768 extension=
".mao", log_filename=.false.)
2779 IF (qs_env%x_data(i, 1)%do_hfx_ri)
CALL print_ri_hfx(qs_env%x_data(i, 1)%ri_data, qs_env)
2783 CALL timestop(handle)
2793 SUBROUTINE hirshfeld_charges(qs_env, input_section, unit_nr)
2796 INTEGER,
INTENT(IN) :: unit_nr
2798 INTEGER :: i, iat, ikind, natom, nkind, nspin, &
2799 radius_type, refc, shapef
2800 INTEGER,
DIMENSION(:),
POINTER :: atom_list
2801 LOGICAL :: do_radius, do_sc, paw_atom
2802 REAL(kind=
dp) :: zeff
2803 REAL(kind=
dp),
DIMENSION(:),
POINTER :: radii
2804 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: charges
2807 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
2813 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2817 NULLIFY (hirshfeld_env)
2821 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
2822 ALLOCATE (hirshfeld_env%charges(natom))
2831 IF (.NOT.
SIZE(radii) == nkind) &
2832 CALL cp_abort(__location__, &
2833 "Length of keyword HIRSHFELD\ATOMIC_RADII does not "// &
2834 "match number of atomic kinds in the input coordinate file.")
2839 iterative=do_sc, ref_charge=refc, &
2840 radius_type=radius_type)
2842 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
2848 nspin =
SIZE(matrix_p, 1)
2849 ALLOCATE (charges(natom, nspin))
2854 atomic_kind => atomic_kind_set(ikind)
2856 DO iat = 1,
SIZE(atom_list)
2858 hirshfeld_env%charges(i) = zeff
2862 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
2865 hirshfeld_env%charges(iat) = sum(charges(iat, :))
2868 cpabort(
"Unknown type of reference charge for Hirshfeld partitioning.")
2872 IF (hirshfeld_env%iterative)
THEN
2879 CALL get_qs_env(qs_env, particle_set=particle_set, dft_control=dft_control)
2880 IF (dft_control%qs_control%gapw)
THEN
2882 CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
2885 atomic_kind => particle_set(iat)%atomic_kind
2887 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
2889 charges(iat, 1:nspin) = charges(iat, 1:nspin) + mp_rho(iat)%q0(1:nspin)
2894 IF (unit_nr > 0)
THEN
2896 qs_kind_set, unit_nr)
2902 DEALLOCATE (charges)
2904 END SUBROUTINE hirshfeld_charges
2914 SUBROUTINE project_function_a(ca, a, cb, b, l)
2916 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: ca
2917 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: a, cb, b
2918 INTEGER,
INTENT(IN) :: l
2921 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
2922 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: smat, tmat, v
2925 ALLOCATE (smat(n, n), tmat(n, n), v(n, 1), ipiv(n))
2929 v(:, 1) = matmul(tmat, cb)
2934 DEALLOCATE (smat, tmat, v, ipiv)
2936 END SUBROUTINE project_function_a
2946 SUBROUTINE project_function_b(ca, a, bfun, grid_atom, l)
2948 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: ca
2949 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: a, bfun
2951 INTEGER,
INTENT(IN) :: l
2953 INTEGER :: i, info, n, nr
2954 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
2955 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: afun
2956 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: smat, v
2960 ALLOCATE (smat(n, n), v(n, 1), ipiv(n), afun(nr))
2964 afun(:) = grid_atom%rad(:)**l*exp(-a(i)*grid_atom%rad2(:))
2965 v(i, 1) = sum(afun(:)*bfun(:)*grid_atom%wr(:))
2971 DEALLOCATE (smat, v, ipiv, afun)
2973 END SUBROUTINE project_function_b
2984 SUBROUTINE qs_scf_post_local_energy(input, logger, qs_env)
2989 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_local_energy'
2991 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
2992 INTEGER :: handle, io_unit, natom, unit_nr
2993 LOGICAL :: append_cube, gapw, gapw_xc, mpi_io
3002 CALL timeset(routinen, handle)
3005 "DFT%PRINT%LOCAL_ENERGY_CUBE"),
cp_p_file))
THEN
3007 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
3008 gapw = dft_control%qs_control%gapw
3009 gapw_xc = dft_control%qs_control%gapw_xc
3010 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3012 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3013 CALL auxbas_pw_pool%create_pw(eden)
3017 append_cube =
section_get_lval(input,
"DFT%PRINT%LOCAL_ENERGY_CUBE%APPEND")
3018 IF (append_cube)
THEN
3019 my_pos_cube =
"APPEND"
3021 my_pos_cube =
"REWIND"
3025 extension=
".cube", middle_name=
"local_energy", &
3026 file_position=my_pos_cube, mpi_io=mpi_io)
3028 unit_nr,
"LOCAL ENERGY", particles=particles, &
3030 "PRINT%LOCAL_ENERGY_CUBE%STRIDE"), mpi_io=mpi_io)
3031 IF (io_unit > 0)
THEN
3032 INQUIRE (unit=unit_nr, name=filename)
3033 IF (gapw .OR. gapw_xc)
THEN
3034 WRITE (unit=io_unit, fmt=
"(/,T3,A,A)") &
3035 "The soft part of the local energy is written to the file: ", trim(adjustl(filename))
3037 WRITE (unit=io_unit, fmt=
"(/,T3,A,A)") &
3038 "The local energy is written to the file: ", trim(adjustl(filename))
3042 "DFT%PRINT%LOCAL_ENERGY_CUBE", mpi_io=mpi_io)
3044 CALL auxbas_pw_pool%give_back_pw(eden)
3046 CALL timestop(handle)
3048 END SUBROUTINE qs_scf_post_local_energy
3059 SUBROUTINE qs_scf_post_local_stress(input, logger, qs_env)
3064 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_local_stress'
3066 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
3067 INTEGER :: handle, io_unit, natom, unit_nr
3068 LOGICAL :: append_cube, gapw, gapw_xc, mpi_io
3069 REAL(kind=
dp) :: beta
3078 CALL timeset(routinen, handle)
3081 "DFT%PRINT%LOCAL_STRESS_CUBE"),
cp_p_file))
THEN
3083 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
3084 gapw = dft_control%qs_control%gapw
3085 gapw_xc = dft_control%qs_control%gapw_xc
3086 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3088 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3089 CALL auxbas_pw_pool%create_pw(stress)
3095 append_cube =
section_get_lval(input,
"DFT%PRINT%LOCAL_STRESS_CUBE%APPEND")
3096 IF (append_cube)
THEN
3097 my_pos_cube =
"APPEND"
3099 my_pos_cube =
"REWIND"
3103 extension=
".cube", middle_name=
"local_stress", &
3104 file_position=my_pos_cube, mpi_io=mpi_io)
3106 unit_nr,
"LOCAL STRESS", particles=particles, &
3108 "PRINT%LOCAL_STRESS_CUBE%STRIDE"), mpi_io=mpi_io)
3109 IF (io_unit > 0)
THEN
3110 INQUIRE (unit=unit_nr, name=filename)
3111 WRITE (unit=io_unit, fmt=
"(/,T3,A)")
"Write 1/3*Tr(sigma) to cube file"
3112 IF (gapw .OR. gapw_xc)
THEN
3113 WRITE (unit=io_unit, fmt=
"(T3,A,A)") &
3114 "The soft part of the local stress is written to the file: ", trim(adjustl(filename))
3116 WRITE (unit=io_unit, fmt=
"(T3,A,A)") &
3117 "The local stress is written to the file: ", trim(adjustl(filename))
3121 "DFT%PRINT%LOCAL_STRESS_CUBE", mpi_io=mpi_io)
3123 CALL auxbas_pw_pool%give_back_pw(stress)
3126 CALL timestop(handle)
3128 END SUBROUTINE qs_scf_post_local_stress
3139 SUBROUTINE qs_scf_post_ps_implicit(input, logger, qs_env)
3144 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_scf_post_ps_implicit'
3146 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
3147 INTEGER :: boundary_condition, handle, i, j, &
3148 n_cstr, n_tiles, unit_nr
3149 LOGICAL :: append_cube, do_cstr_charge_cube, do_dielectric_cube, do_dirichlet_bc_cube, &
3150 has_dirichlet_bc, has_implicit_ps, mpi_io, tile_cubes
3160 CALL timeset(routinen, handle)
3162 NULLIFY (pw_env, auxbas_pw_pool, dft_section, particles)
3165 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3167 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3169 has_implicit_ps = .false.
3170 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
3171 IF (pw_env%poisson_env%parameters%solver .EQ.
pw_poisson_implicit) has_implicit_ps = .true.
3175 "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE"),
cp_p_file)
3176 IF (has_implicit_ps .AND. do_dielectric_cube)
THEN
3177 append_cube =
section_get_lval(input,
"DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%APPEND")
3178 my_pos_cube =
"REWIND"
3179 IF (append_cube)
THEN
3180 my_pos_cube =
"APPEND"
3184 extension=
".cube", middle_name=
"DIELECTRIC_CONSTANT", file_position=my_pos_cube, &
3186 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3187 CALL auxbas_pw_pool%create_pw(aux_r)
3189 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3190 SELECT CASE (boundary_condition)
3192 CALL pw_copy(poisson_env%implicit_env%dielectric%eps, aux_r)
3194 CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
3195 pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
3196 pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
3197 pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
3198 poisson_env%implicit_env%dielectric%eps, aux_r)
3201 CALL cp_pw_to_cube(aux_r, unit_nr,
"DIELECTRIC CONSTANT", particles=particles, &
3202 stride=
section_get_ivals(dft_section,
"PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%STRIDE"), &
3205 "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", mpi_io=mpi_io)
3207 CALL auxbas_pw_pool%give_back_pw(aux_r)
3212 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE"),
cp_p_file)
3214 has_dirichlet_bc = .false.
3215 IF (has_implicit_ps)
THEN
3216 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3218 has_dirichlet_bc = .true.
3222 IF (has_implicit_ps .AND. do_cstr_charge_cube .AND. has_dirichlet_bc)
THEN
3224 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%APPEND")
3225 my_pos_cube =
"REWIND"
3226 IF (append_cube)
THEN
3227 my_pos_cube =
"APPEND"
3231 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", &
3232 extension=
".cube", middle_name=
"dirichlet_cstr_charge", file_position=my_pos_cube, &
3234 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3235 CALL auxbas_pw_pool%create_pw(aux_r)
3237 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3238 SELECT CASE (boundary_condition)
3240 CALL pw_copy(poisson_env%implicit_env%cstr_charge, aux_r)
3242 CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
3243 pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
3244 pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
3245 pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
3246 poisson_env%implicit_env%cstr_charge, aux_r)
3249 CALL cp_pw_to_cube(aux_r, unit_nr,
"DIRICHLET CONSTRAINT CHARGE", particles=particles, &
3250 stride=
section_get_ivals(dft_section,
"PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%STRIDE"), &
3253 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", mpi_io=mpi_io)
3255 CALL auxbas_pw_pool%give_back_pw(aux_r)
3260 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE"),
cp_p_file)
3261 has_dirichlet_bc = .false.
3262 IF (has_implicit_ps)
THEN
3263 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3265 has_dirichlet_bc = .true.
3269 IF (has_implicit_ps .AND. has_dirichlet_bc .AND. do_dirichlet_bc_cube)
THEN
3270 append_cube =
section_get_lval(input,
"DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%APPEND")
3271 my_pos_cube =
"REWIND"
3272 IF (append_cube)
THEN
3273 my_pos_cube =
"APPEND"
3275 tile_cubes =
section_get_lval(input,
"DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%TILE_CUBES")
3277 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3278 CALL auxbas_pw_pool%create_pw(aux_r)
3281 IF (tile_cubes)
THEN
3283 n_cstr =
SIZE(poisson_env%implicit_env%contacts)
3285 n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
3287 filename =
"dirichlet_cstr_"//trim(adjustl(
cp_to_string(j)))// &
3290 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
3291 extension=
".cube", middle_name=filename, file_position=my_pos_cube, &
3294 CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, aux_r)
3296 CALL cp_pw_to_cube(aux_r, unit_nr,
"DIRICHLET TYPE CONSTRAINT", particles=particles, &
3297 stride=
section_get_ivals(dft_section,
"PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
3300 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
3305 NULLIFY (dirichlet_tile)
3306 ALLOCATE (dirichlet_tile)
3307 CALL auxbas_pw_pool%create_pw(dirichlet_tile)
3310 unit_nr =
cp_print_key_unit_nr(logger, input,
"DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
3311 extension=
".cube", middle_name=
"DIRICHLET_CSTR", file_position=my_pos_cube, &
3314 n_cstr =
SIZE(poisson_env%implicit_env%contacts)
3316 n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
3318 CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, dirichlet_tile)
3319 CALL pw_axpy(dirichlet_tile, aux_r)
3323 CALL cp_pw_to_cube(aux_r, unit_nr,
"DIRICHLET TYPE CONSTRAINT", particles=particles, &
3324 stride=
section_get_ivals(dft_section,
"PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
3327 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
3328 CALL auxbas_pw_pool%give_back_pw(dirichlet_tile)
3329 DEALLOCATE (dirichlet_tile)
3332 CALL auxbas_pw_pool%give_back_pw(aux_r)
3335 CALL timestop(handle)
3337 END SUBROUTINE qs_scf_post_ps_implicit
3345 SUBROUTINE write_adjacency_matrix(qs_env, input)
3349 CHARACTER(len=*),
PARAMETER :: routinen =
'write_adjacency_matrix'
3351 INTEGER :: adjm_size, colind, handle, iatom, ikind, &
3352 ind, jatom, jkind, k, natom, nkind, &
3353 output_unit, rowind, unit_nr
3354 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: interact_adjm
3355 LOGICAL :: do_adjm_write, do_symmetric
3361 DIMENSION(:),
POINTER :: nl_iterator
3364 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3367 CALL timeset(routinen, handle)
3369 NULLIFY (dft_section)
3378 IF (do_adjm_write)
THEN
3379 NULLIFY (qs_kind_set, nl_iterator)
3380 NULLIFY (basis_set_list_a, basis_set_list_b, basis_set_a, basis_set_b)
3382 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_orb=nl, natom=natom, para_env=para_env)
3384 nkind =
SIZE(qs_kind_set)
3385 cpassert(
SIZE(nl) .GT. 0)
3387 cpassert(do_symmetric)
3388 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
3392 adjm_size = ((natom + 1)*natom)/2
3393 ALLOCATE (interact_adjm(4*adjm_size))
3396 NULLIFY (nl_iterator)
3400 ikind=ikind, jkind=jkind, &
3401 iatom=iatom, jatom=jatom)
3403 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
3404 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
3405 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
3406 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
3409 IF (iatom .LE. jatom)
THEN
3416 ikind = ikind + jkind
3417 jkind = ikind - jkind
3418 ikind = ikind - jkind
3422 ind = adjm_size - (natom - rowind + 1)*((natom - rowind + 1) + 1)/2 + colind - rowind + 1
3425 interact_adjm((ind - 1)*4 + 1) = rowind
3426 interact_adjm((ind - 1)*4 + 2) = colind
3427 interact_adjm((ind - 1)*4 + 3) = ikind
3428 interact_adjm((ind - 1)*4 + 4) = jkind
3431 CALL para_env%sum(interact_adjm)
3434 extension=
".adjmat", file_form=
"FORMATTED", &
3435 file_status=
"REPLACE")
3436 IF (unit_nr .GT. 0)
THEN
3437 WRITE (unit_nr,
"(1A,2X,1A,5X,1A,4X,A5,3X,A5)")
"#",
"iatom",
"jatom",
"ikind",
"jkind"
3438 DO k = 1, 4*adjm_size, 4
3440 IF (interact_adjm(k) .GT. 0 .AND. interact_adjm(k + 1) .GT. 0)
THEN
3441 WRITE (unit_nr,
"(I8,2X,I8,3X,I6,2X,I6)") interact_adjm(k:k + 3)
3449 DEALLOCATE (basis_set_list_a, basis_set_list_b)
3452 CALL timestop(handle)
3454 END SUBROUTINE write_adjacency_matrix
3462 SUBROUTINE update_hartree_with_mp2(rho, qs_env)
3466 LOGICAL :: use_virial
3476 NULLIFY (auxbas_pw_pool, pw_env, poisson_env, energy, rho_core, v_hartree_rspace, virial)
3477 CALL get_qs_env(qs_env, pw_env=pw_env, energy=energy, &
3478 rho_core=rho_core, virial=virial, &
3479 v_hartree_rspace=v_hartree_rspace)
3481 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
3483 IF (.NOT. use_virial)
THEN
3485 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
3486 poisson_env=poisson_env)
3487 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
3488 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
3492 v_hartree_gspace, rho_core=rho_core)
3494 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
3495 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
3497 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
3498 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
3501 END SUBROUTINE update_hartree_with_mp2
Types and set/get functions for auxiliary density matrix methods.
Contains methods used in the context of density fitting.
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
subroutine, public sg_overlap(smat, l, pa, pb)
...
calculate the orbitals for a given atomic kind type
subroutine, public calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, allelectron, confine)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
recursive subroutine, public get_ddapc(qs_env, calc_force, density_fit_section, density_type, qout1, qout2, out_radii, dq_out, ext_rho_tot_g, itype_of_density, iwc)
Computes the Density Derived Atomic Point Charges.
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
the type I Discrete Cosine Transform (DCT-I)
subroutine, public pw_shrink(neumann_directions, dests_shrink, srcs_shrink, bounds_local_shftd, pw_in, pw_shrinked)
shrinks an evenly symmetric pw_r3d_rs_type data to a pw_r3d_rs_type data that is 8 times smaller (the...
Calculate Energy Decomposition analysis.
subroutine, public edmf_analysis(qs_env, input_section, unit_nr)
...
Definition and initialisation of the et_coupling data type.
subroutine, public set_et_coupling_type(et_coupling, et_mo_coeff, rest_mat)
...
subroutine, public print_ri_hfx(ri_data, qs_env)
Print RI-HFX quantities, as required by the PRINT subsection.
Calculate Hirshfeld charges and related functions.
subroutine, public comp_hirshfeld_charges(qs_env, hirshfeld_env, charges)
...
subroutine, public create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, radius, radii_list)
creates kind specific shape functions for Hirshfeld charges
subroutine, public write_hirshfeld_charges(charges, hirshfeld_env, particle_set, qs_kind_set, unit_nr)
...
subroutine, public comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, ounit)
...
subroutine, public save_hirshfeld_charges(charges, particle_set, qs_kind_set, qs_env)
saves the Hirshfeld charges to the results structure
The types needed for the calculation of Hirshfeld charges and related functions.
subroutine, public create_hirshfeld_type(hirshfeld_env)
...
subroutine, public set_hirshfeld_info(hirshfeld_env, shape_function_type, iterative, ref_charge, fnorm, radius_type, use_bohr)
Set values of a Hirshfeld env.
subroutine, public release_hirshfeld_type(hirshfeld_env)
...
Calculate intrinsic atomic orbitals and analyze wavefunctions.
subroutine, public iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
...
Calculate ntrinsic atomic orbitals and analyze wavefunctions.
subroutine, public iao_read_input(iao_env, iao_section, cell)
...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Types and basic routines needed for a kpoint calculation.
Interface to the LAPACK F77 library.
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_analysis(qs_env, input_section, unit_nr)
...
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Utility routines for the memory handling.
Interface to the message passing library MPI.
Calculate localized minimal basis and analyze wavefunctions.
subroutine, public minbas_analysis(qs_env, input_section, unit_nr)
...
Functions handling the MOLDEN format. Split from mode_selective.
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section)
Write out the MOs in molden format for visualisation.
Define the data structure for the molecule information.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:, :), allocatable, public indso
represent a simple array based list of the given type
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public evolt
real(kind=dp), parameter, public angstrom
Provide various population analyses and print the requested output information.
subroutine, public lowdin_population_analysis(qs_env, output_unit, print_level)
Perform a Lowdin population analysis based on a symmetric orthogonalisation of the density matrix usi...
subroutine, public mulliken_population_analysis(qs_env, output_unit, print_level)
Perform a Mulliken population analysis.
computes preconditioners, and implements methods to apply them currently used in qs_ot
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
This module defines the grid data type and some basic operations on it.
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_implicit
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Write wfx file, works as interface to chargemol and multiwfn.
subroutine, public write_wfx(qs_env, dft_section)
...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
Calculation of commutator [H,r] matrices.
subroutine, public build_com_hr_matrix(qs_env, matrix_hr)
Calculation of the [H,r] commutators matrices over Cartesian Gaussian functions.
Calculation of the energies concerning the core charge distribution.
Calculation and writing of density of states.
subroutine, public calculate_dos_kp(kpoints, qs_env, dft_section)
Compute and write density of states (kpoints)
subroutine, public calculate_dos(mos, dft_section)
Compute and write density of states.
Calculates electric field gradients H.M. Petrili, P.E. Blochl, P. Blaha, K. Schwarz,...
subroutine, public qs_efg_calc(qs_env)
...
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public qs_elf_calc(qs_env, elf_r, rho_cutoff)
...
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public energy_windows(qs_env)
...
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Calculates hyperfine values.
subroutine, public qs_epr_hyp_calc(qs_env)
...
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
subroutine, public loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
Computes and prints the Dipole (using localized charges)
subroutine, public get_localization_info(qs_env, qs_loc_env, loc_section, mo_local, wf_r, wf_g, particles, coeff, evals, marked_states)
Performs localization of the orbitals.
New version of the module for the localization of the molecular orbitals This should be able to use d...
subroutine, public qs_loc_env_release(qs_loc_env)
...
subroutine, public qs_loc_env_create(qs_loc_env)
...
Some utilities for the construction of the localization environment.
subroutine, public loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, do_homo, evals, do_mixed)
...
subroutine, public qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, loc_coeff, mo_loc_history)
allocates the data, and initializes the operators
subroutine, public qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, do_xas, nloc_xas, spin_xas)
initializes everything needed for localization of the HOMOs
subroutine, public retain_history(mo_loc_history, mo_loc)
copy old mos to new ones, allocating as necessary
subroutine, public qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, do_homo, do_mo_cubes, mo_loc_history, evals, tot_zeff_corr, do_mixed)
initializes everything needed for localization of the molecular orbitals
Routines for calculating local energy and stress tensor.
subroutine, public qs_local_stress(qs_env, stress_tensor, beta)
Routine to calculate the local stress.
subroutine, public qs_local_energy(qs_env, energy_density)
Routine to calculate the local energy.
Definition and initialisation of the mo data type.
subroutine, public write_dm_binary_restart(mo_array, dft_section, tmpl_matrix)
calculates density matrix from mo set and writes the density matrix into a binary restart file
collects routines that perform operations directly related to MOs
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env)
Calculate KS eigenvalues starting from OF MOS.
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
subroutine, public qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
...
subroutine, public qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Calculation and writing of projected density of states The DOS is computed per angular momentum and p...
subroutine, public calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, dft_section, ispin, xas_mittle, external_matrix_shalf)
Compute and write projected density of states.
provides a resp fit for gas phase systems
subroutine, public resp_fit(qs_env)
performs resp fit and generates RESP charges
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, qlm_gg, qlm_car, qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Functions to print the KS and S matrix in the CSR format to file.
subroutine, public write_s_matrix_csr(qs_env, input)
writing the overlap matrix in csr format into a file
subroutine, public write_ks_matrix_csr(qs_env, input)
writing the KS matrix in csr format into a file
subroutine, public qs_scf_write_mos(qs_env, scf_env, final_mos)
Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit.
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
Gets the lumos, and eigenvalues for the lumos.
subroutine, public write_mo_free_results(qs_env)
Write QS results always available (if switched on through the print_keys) Can be called from ls_scf.
subroutine, public qs_scf_post_moments(input, logger, qs_env, output_unit)
Computes and prints electric moments.
subroutine, public write_mo_dependent_results(qs_env, scf_env)
Write QS results available if MO's are present (if switched on through the print_keys) Writes only MO...
subroutine, public scf_post_calculation_gpw(qs_env, wf_type, do_mp2)
collects possible post - scf calculations and prints info / computes properties.
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, output_unit, unoccupied_orbs, scf_env, matrix_s, marked_states, for_rtp)
writes a new 'mixed' set of mos to restart file, without touching the current MOs
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Interface to Wannier90 code.
subroutine, public wannier90_interface(input, logger, qs_env)
...
Methods related to (\cal S)^2 (i.e. spin)
subroutine, public compute_s_square(mos, matrix_s, s_square, s_square_ideal, mo_derivs, strength)
Compute the expectation value <(\cal S)^2> of the single determinant defined by the spin up (alpha) a...
parameters that control an scf iteration
Calculation of STM image as post processing of an electronic structure calculation,...
subroutine, public th_stm_image(qs_env, stm_section, particles, unoccupied_orbs, unoccupied_evals)
Driver for the calculation of STM image, as post processing of a ground-state electronic structure ca...
routines for DFT+NEGF calculations (coupling with the quantum transport code OMEN)
subroutine, public qs_scf_post_transport(qs_env)
post scf calculations for transport
Interface for Voronoi Integration and output of BQB files.
subroutine, public entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
Does a Voronoi integration of density or stores the density to compressed BQB format.
subroutine, public xray_diffraction_spectrum(qs_env, unit_number, q_max)
Calculate the coherent X-ray diffraction spectrum using the total electronic density in reciprocal sp...
subroutine, public calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, rhotot_elec_gspace, q_max, rho_hard, rho_soft, fsign)
The total electronic density in reciprocal space (g-space) is calculated.
stores some data used in wavefunction fitting
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
represent a pointer to a 1d array
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
quantities needed for a Hirshfeld based partitioning of real space
Contains information about kpoints.
stores all the informations relevant to an mpi environment
represent a list of objects
contained for different pw related things
environment for the poisson solver
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
keeps the density in various representations, keeping track of which ones are valid.