Skip to content

Commit

Permalink
print real and imaginary parts of the MO coefficients during RTP
Browse files Browse the repository at this point in the history
  • Loading branch information
@mattiatj @fstein93
mattiatj authored and fstein93 committed Nov 28, 2022
1 parent 6e16aa4 commit 9aeaa09
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 15 deletions.
84 changes: 83 additions & 1 deletion src/emd/rt_propagation_output.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ MODULE rt_propagation_output
cp_fm_p_type,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
Expand All @@ -50,6 +51,7 @@ MODULE rt_propagation_output
section_vals_type
USE kahan_sum, ONLY: accurate_sum
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE machine, ONLY: m_flush
USE message_passing, ONLY: mp_max
Expand All @@ -73,7 +75,13 @@ MODULE rt_propagation_output
qs_kind_type
USE qs_linres_current, ONLY: calculate_jrho_resp
USE qs_linres_types, ONLY: current_env_type
USE qs_mo_io, ONLY: write_rt_mos_to_restart
USE qs_mo_io, ONLY: write_mo_set_to_output_unit,&
write_rt_mos_to_restart
USE qs_mo_types, ONLY: allocate_mo_set,&
deallocate_mo_set,&
get_mo_set,&
init_mo_set,&
mo_set_type
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_post_gpw, ONLY: qs_scf_post_moments,&
Expand Down Expand Up @@ -240,6 +248,7 @@ SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)
END IF
CALL write_rt_mos_to_restart(qs_env%mos, mos_new, particle_set, &
dft_section, qs_kind_set)
CALL write_rtp_mos_to_output_unit(qs_env, mos_new)
END IF
END IF

Expand All @@ -251,6 +260,79 @@ SUBROUTINE rt_prop_output(qs_env, run_type, delta_iter, used_time)

END SUBROUTINE rt_prop_output

! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param mos_new ...
! **************************************************************************************************
SUBROUTINE write_rtp_mos_to_output_unit(qs_env, mos_new)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new

CHARACTER(len=*), PARAMETER :: routineN = 'write_rtp_mos_to_output_unit'

CHARACTER(LEN=10) :: spin
CHARACTER(LEN=2*default_string_length) :: name
INTEGER :: i, ispin, nao, nelectron, nmo, nspins
REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(mo_set_type) :: mo_set_rtp
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: dft_section, input

NULLIFY (atomic_kind_set, particle_set, qs_kind_set, input, mos)

CALL get_qs_env(qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
input=input, &
mos=mos)

NULLIFY (dft_section)
dft_section => section_vals_get_subs_vals(input, "DFT")
nspins = SIZE(mos_new)/2

DO ispin = 1, nspins
! initiate mo_set
CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo, nelectron=nelectron, &
n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)

CALL allocate_mo_set(mo_set_rtp, &
nao=nao, &
nmo=nmo, &
nelectron=nelectron, &
n_el_f=n_el_f, &
maxocc=maxocc, &
flexible_electron_count=flexible_electron_count)

WRITE (name, fmt="(a,I6)") "RTP MO SET, SPIN ", ispin
CALL init_mo_set(mo_set_rtp, fm_ref=mos_new(2*ispin - 1)%matrix, name=name)

IF (nspins > 1) THEN
IF (ispin == 1) THEN
spin = "ALPHA SPIN"
ELSE
spin = "BETA SPIN"
END IF
ELSE
spin = ""
END IF

!loop for real (odd) and imaginary (even) parts
DO i = 1, 0, -1
CALL cp_fm_to_fm(mos_new(2*ispin - i)%matrix, mo_set_rtp%mo_coeff)
CALL write_mo_set_to_output_unit(mo_set_rtp, atomic_kind_set, qs_kind_set, particle_set, &
dft_section, 4, 0, rtp=.TRUE., spin=TRIM(spin), cpart=MOD(i, 2), sim_step=qs_env%sim_step)
END DO

CALL deallocate_mo_set(mo_set_rtp)
END DO

END SUBROUTINE write_rtp_mos_to_output_unit

! **************************************************************************************************
!> \brief computes the effective orthonormality of a set of mos given an s-matrix
!> orthonormality is the max deviation from unity of the C^T S C
Expand Down
69 changes: 55 additions & 14 deletions src/qs_mo_io.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -968,6 +968,9 @@ END SUBROUTINE read_mos_restart_low
!> \param final_mos ...
!> \param spin ...
!> \param solver_method ...
!> \param rtp ...
!> \param cpart ...
!> \param sim_step ...
!> \date 15.05.2001
!> \par History:
!> - Optionally print Cartesian MOs (20.04.2005, MK)
Expand All @@ -980,7 +983,7 @@ END SUBROUTINE read_mos_restart_low
! **************************************************************************************************
SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
dft_section, before, kpoint, final_mos, spin, &
solver_method)
solver_method, rtp, cpart, sim_step)
TYPE(mo_set_type), INTENT(IN) :: mo_set
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
Expand All @@ -991,27 +994,32 @@ SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, par
LOGICAL, INTENT(IN), OPTIONAL :: final_mos
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: spin
CHARACTER(LEN=2), INTENT(IN), OPTIONAL :: solver_method
LOGICAL, INTENT(IN), OPTIONAL :: rtp
INTEGER, INTENT(IN), OPTIONAL :: cpart, sim_step
CHARACTER(LEN=12) :: symbol
CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
CHARACTER(LEN=14) :: fmtstr5
CHARACTER(LEN=15) :: energy_str, orbital_str, vector_str
CHARACTER(LEN=15) :: energy_str, orbital_str, step_string
CHARACTER(LEN=2) :: element_symbol, my_solver_method
CHARACTER(LEN=2*default_string_length) :: name
CHARACTER(LEN=21) :: vector_str
CHARACTER(LEN=22) :: fmtstr4
CHARACTER(LEN=24) :: fmtstr2
CHARACTER(LEN=25) :: fmtstr1
CHARACTER(LEN=29) :: fmtstr6
CHARACTER(LEN=4) :: reim
CHARACTER(LEN=40) :: fmtstr3
CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
INTEGER :: after, first_mo, from, iatom, icgf, ico, icol, ikind, imo, irow, iset, isgf, &
ishell, iso, iw, jcol, last_mo, left, lmax, lshell, natom, ncgf, ncol, ncol_global, &
nrow_global, nset, nsgf, right, scf_step, to, width
INTEGER, DIMENSION(:), POINTER :: mo_index_range, nshell
INTEGER, DIMENSION(:, :), POINTER :: l
LOGICAL :: ionode, my_final, print_cartesian, &
print_eigvals, print_eigvecs, &
print_occup, should_output
LOGICAL :: ionode, my_final, my_rtp, &
print_cartesian, print_eigvals, &
print_eigvecs, print_occup, &
should_output
REAL(KIND=dp) :: gap
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
TYPE(cp_logger_type), POINTER :: logger
Expand Down Expand Up @@ -1041,12 +1049,38 @@ SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, par
my_final = .FALSE.
END IF
! complex MOS for RTP, no eigenvalues
my_rtp = .FALSE.
IF (PRESENT(rtp)) THEN
my_rtp = rtp
END IF
should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
"PRINT%MO"), cp_p_file) .OR. my_final
IF ((.NOT. should_output) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup))) RETURN
scf_step = MAX(0, logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1)
IF (my_rtp) THEN
CPASSERT(PRESENT(sim_step))
CPASSERT(PRESENT(cpart))
scf_step = sim_step
IF (cpart == 0) THEN
reim = "IMAG"
ELSE
reim = "REAL"
END IF
print_eigvals = .FALSE.
ELSE
scf_step = MAX(0, logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1)
END IF
IF (.NOT. my_final) THEN
IF (.NOT. my_rtp) THEN
step_string = " AFTER SCF STEP"
ELSE
step_string = " AFTER RTP STEP"
END IF
END IF
IF (PRESENT(solver_method)) THEN
my_solver_method = solver_method
Expand Down Expand Up @@ -1107,6 +1141,11 @@ SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, par
vector_str = "COEFFICIENTS"
END IF
IF (my_rtp) THEN
energy_str = "ZEROS"
vector_str = TRIM(reim)//" RTP COEFFICIENTS"
END IF
IF (print_eigvecs) THEN
IF (print_cartesian) THEN
Expand Down Expand Up @@ -1168,13 +1207,13 @@ SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, par
TRIM(orbital_str)//" "//TRIM(vector_str)
IF (.NOT. my_final) &
WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//" AFTER SCF STEP", scf_step
WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
ELSE IF (print_occup .OR. print_eigvals) THEN
name = TRIM(energy_str)//" AND OCCUPATION NUMBERS"
IF (.NOT. my_final) &
WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//" AFTER SCF STEP", scf_step
WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
END IF ! print_eigvecs
! Print headline
Expand Down Expand Up @@ -1358,12 +1397,14 @@ SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, par
END IF ! print_eigvecs
fmtstr6 = "(A,T18,F17. ,A,T41,F17. ,A)"
WRITE (UNIT=fmtstr6(12:13), FMT="(I2)") after
WRITE (UNIT=fmtstr6(25:26), FMT="(I2)") after
WRITE (UNIT=iw, FMT=fmtstr6) &
" MO| E(Fermi):", mo_set%mu, " a.u.", mo_set%mu*evolt, " eV"
IF ((mo_set%homo > 0)) THEN
IF (.NOT. my_rtp) THEN
fmtstr6 = "(A,T18,F17. ,A,T41,F17. ,A)"
WRITE (UNIT=fmtstr6(12:13), FMT="(I2)") after
WRITE (UNIT=fmtstr6(25:26), FMT="(I2)") after
WRITE (UNIT=iw, FMT=fmtstr6) &
" MO| E(Fermi):", mo_set%mu, " a.u.", mo_set%mu*evolt, " eV"
END IF
IF ((mo_set%homo > 0) .AND. .NOT. my_rtp) THEN
IF ((mo_set%occupation_numbers(mo_set%homo) == mo_set%maxocc) .AND. (last_mo > mo_set%homo)) THEN
gap = mo_set%eigenvalues(mo_set%homo + 1) - &
mo_set%eigenvalues(mo_set%homo)
Expand Down

0 comments on commit 9aeaa09

Please sign in to comment.