mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-27 05:35:28 -04:00
Allow virtual orbital output in molden files with OT (#4924)
Get virtual orbital energies and coefficients from existing OT eigensolver and output it to .molden, which can meet a part of the feature request in #4353. Users should generate this function through FORCE_EVAL%DFT%PRINT%MO_MOLDEN%OT_NLMO keyword (e.g. OT_NLUMO 30).
This commit is contained in:
SY Wang
GitHub
parent
0ca8b53e62
commit
ff0bbe84df
3 changed files with 173 additions and 14 deletions
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue