diff --git a/src/input_cp2k_print_dft.F b/src/input_cp2k_print_dft.F index eb52faa598..59f9b21145 100644 --- a/src/input_cp2k_print_dft.F +++ b/src/input_cp2k_print_dft.F @@ -616,6 +616,17 @@ CONTAINS enum_i_vals=[gto_cartesian, gto_spherical]) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="OT_NLUMO", & + description="Number of unoccupied orbitals to include in the Molden file. "// & + "Triggers a post-SCF eigensolver to obtain virtual orbital energies and "// & + "coefficients. Only available with OT method, since diagonalization "// & + "can already compute unoccupied orbitals via `ADDED_MOS` during "// & + "SCF and output them to the Molden file after convergence. "// & + "0 means no virtual orbitals, -1 means all available.", & + usage="OT_NLUMO {int}", & + default_i_val=0) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) CALL section_add_subsection(section, print_key) CALL section_release(print_key) diff --git a/src/molden_utils.F b/src/molden_utils.F index bacd7ea169..04fcdc47ee 100644 --- a/src/molden_utils.F +++ b/src/molden_utils.F @@ -13,9 +13,11 @@ MODULE molden_utils USE atomic_kind_types, ONLY: get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& gto_basis_set_type + USE cp_array_utils, ONLY: cp_1d_r_p_type USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm USE cp_fm_types, ONLY: cp_fm_get_info,& - cp_fm_get_submatrix + cp_fm_get_submatrix,& + cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_type USE cp_output_handling, ONLY: cp_p_file,& @@ -59,13 +61,20 @@ CONTAINS !> \param qs_kind_set for basis set info !> \param particle_set particles data structure, for positions and kinds !> \param print_section input section containing relevant print key +!> \param unoccupied_orbs optional: unoccupied orbital coefficients from make_lumo_gpw +!> \param unoccupied_evals optional: unoccupied orbital eigenvalues !> \author MattW, IainB ! ************************************************************************************************** - SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section) + SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section, & + unoccupied_orbs, unoccupied_evals) TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(section_vals_type), POINTER :: print_section + TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), & + OPTIONAL :: unoccupied_orbs + TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(IN), & + OPTIONAL :: unoccupied_evals CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mos_molden' CHARACTER(LEN=molden_lmax+1), PARAMETER :: angmom = "spdfg" @@ -73,7 +82,7 @@ CONTAINS CHARACTER(LEN=15) :: fmtstr1, fmtstr2 CHARACTER(LEN=2) :: element_symbol INTEGER :: gto_kind, handle, i, iatom, icgf, icol, ikind, ipgf, irow, irow_in, iset, isgf, & - ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, z + ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, numos, z INTEGER, DIMENSION(:), POINTER :: npgf, nshell INTEGER, DIMENSION(:, :), POINTER :: l INTEGER, DIMENSION(molden_ncomax, 0:molden_lmax) :: orbmap @@ -345,6 +354,103 @@ CONTAINS IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix) END DO + ! Write unoccupied (virtual) orbitals if provided; only used with OT + IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN + DO ispin = 1, SIZE(unoccupied_orbs) + CALL cp_fm_get_info(unoccupied_orbs(ispin), & + nrow_global=nrow_global, & + ncol_global=numos) + ALLOCATE (smatrix(nrow_global, numos)) + CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix) + + IF (iw > 0) THEN + IF (gto_kind == gto_cartesian) THEN + CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf) + ALLOCATE (cmatrix(ncgf, numos)) + cmatrix = 0.0_dp + + icgf = 1 + isgf = 1 + DO iatom = 1, SIZE(particle_set) + NULLIFY (orb_basis_set) + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) + CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) + IF (ASSOCIATED(orb_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + nset=nset, nshell=nshell, l=l) + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + CALL dgemm("T", "N", nco(lshell), numos, nso(lshell), 1.0_dp, & + orbtramat(lshell)%c2s, nso(lshell), & + smatrix(isgf, 1), nsgf, 0.0_dp, & + cmatrix(icgf, 1), ncgf) + icgf = icgf + nco(lshell) + isgf = isgf + nso(lshell) + END DO + END DO + END IF + END DO + END IF + + DO icol = 1, numos + irow = 1 + irow_in = 1 + + WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol) + IF (ispin < 2) THEN + WRITE (iw, '(A)') 'Spin= Alpha' + ELSE + WRITE (iw, '(A)') 'Spin= Beta' + END IF + WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp + + DO iatom = 1, SIZE(particle_set) + NULLIFY (orb_basis_set) + CALL get_atomic_kind(particle_set(iatom)%atomic_kind, & + element_symbol=element_symbol, kind_number=ikind) + CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) + IF (ASSOCIATED(orb_basis_set)) THEN + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + nset=nset, nshell=nshell, l=l) + + IF (gto_kind == gto_cartesian) THEN + icgf = 1 + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + IF (lshell <= molden_lmax) THEN + CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), & + cmatrix(irow:irow + nco(lshell) - 1, icol)) + irow_in = irow_in + nco(lshell) + END IF + irow = irow + nco(lshell) + END DO + END DO + ELSE IF (gto_kind == gto_spherical) THEN + DO iset = 1, nset + DO ishell = 1, nshell(iset) + lshell = l(ishell, iset) + IF (lshell <= molden_lmax) THEN + CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), & + smatrix(irow:irow + nso(lshell) - 1, icol)) + irow_in = irow_in + nso(lshell) + END IF + irow = irow + nso(lshell) + END DO + END DO + END IF + + END IF + END DO ! iatom + END DO ! icol + END IF + + IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix) + IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix) + END DO ! ispin + END IF + CALL cp_print_key_finished_output(iw, logger, print_section, "") END IF diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 6bf26d8133..d7be9fbd46 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -548,10 +548,8 @@ CONTAINS warning_cube_kpoint = "Print MO cubes not implemented for k-point calculations", & warning_openpmd_kpoint = "Writing to openPMD not implemented for k-point calculations" - INTEGER :: handle, homo, ispin, min_lumos, n_rep, & - nchk_nmoloc, nhomo, nlumo, nlumo_stm, & - nlumos, nmo, nspins, output_unit, & - unit_nr + INTEGER :: handle, homo, ispin, min_lumos, n_rep, nchk_nmoloc, nhomo, nlumo, nlumo_molden, & + nlumo_stm, nlumos, nmo, nspins, output_unit, unit_nr INTEGER, DIMENSION(:, :, :), POINTER :: marked_states LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mixed, do_stm, & do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_do_mp2, & @@ -595,7 +593,7 @@ CONTAINS TYPE(scf_control_type), POINTER :: scf_control TYPE(section_vals_type), POINTER :: dft_section, input, loc_print_section, & localize_section, print_key, & - stm_section + sprint_section, stm_section CALL timeset(routineN, handle) @@ -1099,6 +1097,33 @@ CONTAINS END IF END IF + ! Write molden file including unoccupied orbitals for OT calculations + IF (.NOT. do_kpoints) THEN + sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN") + CALL section_vals_val_get(sprint_section, "OT_NLUMO", i_val=nlumo_molden) + IF (nlumo_molden /= 0 .AND. scf_env%method == ot_method_nr) THEN + CALL get_qs_env(qs_env, mos=mos, qs_kind_set=qs_kind_set, & + particle_set=particle_set) + ALLOCATE (unoccupied_orbs(dft_control%nspins)) + ALLOCATE (unoccupied_evals(dft_control%nspins)) + CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, & + nlumo_molden, nlumos) + IF (output_unit > 0) THEN + WRITE (output_unit, '(/,T2,A,I6,A)') & + "MO_MOLDEN| Writing ", nlumos, " unoccupied orbitals to molden file" + END IF + CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section, & + unoccupied_orbs=unoccupied_orbs, & + unoccupied_evals=unoccupied_evals) + DO ispin = 1, dft_control%nspins + DEALLOCATE (unoccupied_evals(ispin)%array) + CALL cp_fm_release(unoccupied_orbs(ispin)) + END DO + DEALLOCATE (unoccupied_evals) + DEALLOCATE (unoccupied_orbs) + END IF + END IF + ! Print coherent X-ray diffraction spectrum CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit) @@ -2025,8 +2050,10 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_dependent_results' - INTEGER :: handle, homo, ispin, nmo, output_unit - LOGICAL :: all_equal, do_kpoints, explicit + INTEGER :: handle, homo, ispin, nlumo_molden, nmo, & + output_unit + LOGICAL :: all_equal, defer_molden, do_kpoints, & + explicit REAL(KIND=dp) :: maxocc, s_square, s_square_ideal, & total_abs_spin_dens, total_spin_dens REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation_numbers @@ -2095,17 +2122,32 @@ CONTAINS IF (explicit) THEN CALL write_trexio(qs_env, trexio_section) END IF + sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN") IF (.NOT. do_kpoints) THEN CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv) CALL write_dm_binary_restart(mos, dft_section, ks_rmpv) - sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN") - CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section) + ! Check if molden write should be deferred for OT unoccupied orbitals + defer_molden = .FALSE. + CALL section_vals_val_get(sprint_section, "OT_NLUMO", i_val=nlumo_molden) + IF (nlumo_molden /= 0 .AND. PRESENT(scf_env)) THEN + IF (scf_env%method == ot_method_nr) defer_molden = .TRUE. + END IF + IF (.NOT. defer_molden) THEN + CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section) + END IF ! Write Chargemol .wfx - IF (BTEST(cp_print_key_should_output(logger%iter_info, & - dft_section, "PRINT%CHARGEMOL"), & + IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%CHARGEMOL"), & cp_p_file)) THEN CALL write_wfx(qs_env, dft_section) END IF + ELSE + IF (BTEST(cp_print_key_should_output(logger%iter_info, sprint_section, ""), cp_p_file)) THEN + CPWARN("Molden format output is not possible for k-point calculations.") + END IF + IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%CHARGEMOL"), & + cp_p_file)) THEN + CPWARN("Chargemol .wfx format output is not possible for k-point calculations.") + END IF END IF ! K-point MO wavefunction dump