From 9e9dfdb71266d8fd1d4091853def9925083ed0f4 Mon Sep 17 00:00:00 2001 From: Sergey Chulkov Date: Mon, 13 Aug 2018 11:30:11 +0000 Subject: [PATCH] Fix TDDFPT oscillator strengths and print them in scientific format svn-origin-rev: 18584 --- src/fm/cp_fm_basic_linalg.F | 75 ++++ src/qs_tddfpt2_methods.F | 869 +++++++++++++++++++++--------------- 2 files changed, 589 insertions(+), 355 deletions(-) diff --git a/src/fm/cp_fm_basic_linalg.F b/src/fm/cp_fm_basic_linalg.F index 487a62da52..660e061251 100644 --- a/src/fm/cp_fm_basic_linalg.F +++ b/src/fm/cp_fm_basic_linalg.F @@ -64,6 +64,7 @@ MODULE cp_fm_basic_linalg INTERFACE cp_fm_trace MODULE PROCEDURE cp_fm_trace_a0b0t0 + MODULE PROCEDURE cp_fm_trace_a1b0t1 MODULE PROCEDURE cp_fm_trace_a1b1t1 END INTERFACE cp_fm_trace @@ -755,6 +756,80 @@ CONTAINS END SUBROUTINE cp_fm_trace_a0b0t0 +! ************************************************************************************************** +!> \brief Compute trace(k) = Tr (matrix_a(k)^T matrix_b) for each pair of matrices A_k and B. +!> \param matrix_a list of A matrices +!> \param matrix_b B matrix +!> \param trace computed traces +!> \par History +!> * 08.2018 forked from cp_fm_trace() [Sergey Chulkov] +!> \note \parblock +!> Computing the trace requires collective communication between involved MPI processes +!> that implies a synchronisation point between them. The aim of this subroutine is to reduce +!> the amount of time wasted in such synchronisation by performing one large collective +!> operation which involves all the matrices in question. +!> +!> The subroutine's suffix reflects dimensionality of dummy arrays; 'a1b0t1' means that +!> the dummy variables 'matrix_a' and 'trace' are 1-dimensional arrays, while the variable +!> 'matrix_b' is a single matrix. +!> \endparblock +! ************************************************************************************************** + SUBROUTINE cp_fm_trace_a1b0t1(matrix_a, matrix_b, trace) + TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: matrix_a + TYPE(cp_fm_type), POINTER :: matrix_b + REAL(kind=dp), DIMENSION(:), INTENT(out) :: trace + + CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_trace_a1b0t1', & + routineP = moduleN//':'//routineN + + INTEGER :: group, handle, imatrix, n_matrices, & + ncols_local, nrows_local + LOGICAL :: use_sp_a, use_sp_b + REAL(kind=dp), DIMENSION(:, :), POINTER :: ldata_a, ldata_b + REAL(kind=sp), DIMENSION(:, :), POINTER :: ldata_a_sp, ldata_b_sp + + CALL timeset(routineN, handle) + + n_matrices = SIZE(trace) + CPASSERT(SIZE(matrix_a) == n_matrices) + + CALL cp_fm_get_info(matrix_b, nrow_local=nrows_local, ncol_local=ncols_local) + use_sp_b = matrix_b%use_sp + + IF (use_sp_b) THEN + ldata_b_sp => matrix_b%local_data_sp(1:nrows_local, 1:ncols_local) + ELSE + ldata_b => matrix_b%local_data(1:nrows_local, 1:ncols_local) + END IF + +!$OMP PARALLEL DO DEFAULT(NONE), & +!$OMP PRIVATE(imatrix, ldata_a, ldata_a_sp, use_sp_a), & +!$OMP SHARED(ldata_b, ldata_b_sp, matrix_a, matrix_b), & +!$OMP SHARED(ncols_local, nrows_local, n_matrices, trace, use_sp_b) + + DO imatrix = 1, n_matrices + + use_sp_a = matrix_a(imatrix)%matrix%use_sp + + ! assume that the matrices A(i) and B have identical shapes and distribution schemes + IF (use_sp_a .AND. use_sp_b) THEN + ldata_a_sp => matrix_a(imatrix)%matrix%local_data_sp(1:nrows_local, 1:ncols_local) + trace(imatrix) = accurate_dot_product(ldata_a_sp, ldata_b_sp) + ELSE IF (.NOT. use_sp_a .AND. .NOT. use_sp_b) THEN + ldata_a => matrix_a(imatrix)%matrix%local_data(1:nrows_local, 1:ncols_local) + trace(imatrix) = accurate_dot_product(ldata_a, ldata_b) + ELSE + CPABORT("Matrices A and B are of different types") + END IF + END DO +!$OMP END PARALLEL DO + + group = matrix_b%matrix_struct%para_env%group + CALL mp_sum(trace, group) + + CALL timestop(handle) + END SUBROUTINE cp_fm_trace_a1b0t1 + ! ************************************************************************************************** !> \brief Compute trace(k) = Tr (matrix_a(k)^T matrix_b(k)) for each pair of matrices A_k and B_k. !> \param matrix_a list of A matrices diff --git a/src/qs_tddfpt2_methods.F b/src/qs_tddfpt2_methods.F index f4de560cd3..975357a5dc 100644 --- a/src/qs_tddfpt2_methods.F +++ b/src/qs_tddfpt2_methods.F @@ -30,7 +30,6 @@ MODULE qs_tddfpt2_methods cp_fm_contracted_trace,& cp_fm_scale,& cp_fm_scale_and_add,& - cp_fm_schur_product,& cp_fm_trace,& cp_fm_triangular_invert USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose @@ -86,8 +85,7 @@ MODULE qs_tddfpt2_methods int_8 USE machine, ONLY: m_flush,& m_walltime - USE mathconstants, ONLY: pi,& - twopi,& + USE mathconstants, ONLY: twopi,& z_one,& z_zero USE message_passing, ONLY: mp_bcast,& @@ -314,23 +312,19 @@ CONTAINS CHARACTER(len=20) :: nstates_str INTEGER :: energy_unit, handle, ispin, istate, & - iter, log_unit, mult, ndim_periodic, & - niters, nspins, nstates, nstates_read + iter, log_unit, mult, niters, nspins, & + nstates, nstates_read LOGICAL :: do_admm, do_hfx, is_restarted REAL(kind=dp) :: C_hf, conv REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals - REAL(kind=dp), DIMENSION(3) :: kvec, reference_point TYPE(admm_type), POINTER :: admm_env TYPE(cell_type), POINTER :: cell TYPE(cp_blacs_env_type), POINTER :: blacs_env - TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: evects, S_evects + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ, evects, S_evects TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: berry_cossin_xyz, dipole_op_xyz, & - matrix_ks, matrix_s, rRc_xyz + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos - TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_poisson_type), POINTER :: poisson_env TYPE(section_vals_type), POINTER :: hfx_section, input, & tddfpt_print_section, tddfpt_section, & xc_section @@ -346,9 +340,9 @@ CONTAINS CALL cite_reference(Iannuzzi2005) - NULLIFY (blacs_env, cell, dft_control, input, matrix_ks, matrix_s, mos, pw_env) + NULLIFY (blacs_env, cell, dft_control, input, matrix_ks, matrix_s, mos) CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, & - input=input, matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, pw_env=pw_env) + input=input, matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos) tddfpt_control => dft_control%tddfpt2_control nspins = dft_control%nspins @@ -383,76 +377,6 @@ CONTAINS CALL cp_warn(__LOCATION__, "Keyword RKS_TRIPLETS has been ignored for spin-polarised calculations") END IF - ! check that the chosen dipole operator is consistent with the periodic boundary conditions used - CALL pw_env_get(pw_env, poisson_env=poisson_env) - ndim_periodic = COUNT(poisson_env%parameters%periodic == 1) - - ! choose appropriate components of the dipole operator - NULLIFY (dipole_op_xyz) - SELECT CASE (tddfpt_control%dipole_form) - CASE (tddfpt_dipole_berry) - IF (ndim_periodic /= 3) THEN - CALL cp_warn(__LOCATION__, & - "Fully periodic Poisson solver (PERIODIC xyz) is needed "// & - "for oscillator strengths based on the Berry phase formula") - END IF - - NULLIFY (berry_cossin_xyz) - ! index: 1 = Re[exp(-i * G_x * x)], 2 = Im[exp(-i * G_x * x)], - ! 3 = Re[exp(-i * G_y * y)], 4 = Im[exp(-i * G_y * y)], - ! 5 = Re[exp(-i * G_z * z)], 4 = Im[exp(-i * G_z * z)]. - CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2*nderivs) - - DO ispin = 1, nderivs - ! cos component - CALL dbcsr_init_p(berry_cossin_xyz(2*ispin-1)%matrix) - CALL dbcsr_copy(berry_cossin_xyz(2*ispin-1)%matrix, matrix_s(1)%matrix) - CALL dbcsr_set(berry_cossin_xyz(2*ispin-1)%matrix, 0.0_dp) - - ! sin component - CALL dbcsr_init_p(berry_cossin_xyz(2*ispin)%matrix) - CALL dbcsr_copy(berry_cossin_xyz(2*ispin)%matrix, matrix_s(1)%matrix) - CALL dbcsr_set(berry_cossin_xyz(2*ispin)%matrix, 0.0_dp) - - kvec(:) = twopi*cell%h_inv(ispin, :) - CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(2*ispin-1)%matrix, berry_cossin_xyz(2*ispin)%matrix, kvec) - END DO - - ! assign computed components to the generic variable - dipole_op_xyz => berry_cossin_xyz - - CASE (tddfpt_dipole_length) - IF (ndim_periodic /= 0) THEN - CALL cp_warn(__LOCATION__, & - "Non-periodic Poisson solver (PERIODIC none) is needed "// & - "for oscillator strengths based on the length operator") - END IF - - ! compute components of the dipole operator in the length form - NULLIFY (rRc_xyz) - CALL dbcsr_allocate_matrix_set(rRc_xyz, nderivs) - - DO ispin = 1, nderivs - CALL dbcsr_init_p(rRc_xyz(ispin)%matrix) - CALL dbcsr_copy(rRc_xyz(ispin)%matrix, matrix_s(1)%matrix) - END DO - - CALL get_reference_point(reference_point, qs_env=qs_env, & - reference=tddfpt_control%dipole_reference, & - ref_point=tddfpt_control%dipole_ref_point) - - CALL rRc_xyz_ao(op=rRc_xyz, qs_env=qs_env, rc=reference_point, order=1, minimum_image=.FALSE., soft=.FALSE.) - - ! assign computed components to the generic variable - dipole_op_xyz => rRc_xyz - - CASE (tddfpt_dipole_velocity) - dipole_op_xyz => matrix_s(2:nderivs+1) - - CASE DEFAULT - CPABORT("Unimplemented form of the dipole operator") - END SELECT - ! recompute ground-state Kohn-Sham orbitals ALLOCATE (gs_mos(nspins)) DO ispin = 1, nspins @@ -461,6 +385,9 @@ CONTAINS matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix) END DO + ! components of the dipole operator + CALL tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env) + ! multiplicity of molecular system IF (nspins > 1) THEN mult = ABS(SIZE(gs_mos(1)%evals_occ)-SIZE(gs_mos(2)%evals_occ))+1 @@ -611,10 +538,9 @@ CONTAINS ! *** print summary information *** log_unit = cp_logger_get_default_io_unit() - CALL tddfpt_print_summary(log_unit, evects, evals, mult, gs_mos, matrix_s=matrix_s(1)%matrix, & - dipole_op_xyz=dipole_op_xyz, dipole_form=tddfpt_control%dipole_form, & - cell=cell, min_amplitude=tddfpt_control%min_excitation_amplitude, & - fm_pool_ao_mo_occ=work_matrices%fm_pool_ao_mo_occ) + CALL tddfpt_print_summary(log_unit, evects, evals, mult, dipole_op_mos_occ) + CALL tddfpt_print_excitation_analysis(log_unit, evects, gs_mos, matrix_s(1)%matrix, & + min_amplitude=tddfpt_control%min_excitation_amplitude) ! -- clean up all useless stuff DO istate = SIZE(evects, 2), 1, -1 @@ -632,16 +558,20 @@ CONTAINS CALL tddfpt_release_work_matrices(work_matrices, sub_env) CALL tddfpt_sub_env_release(sub_env) + IF (ALLOCATED(dipole_op_mos_occ)) THEN + DO ispin = nspins, 1, -1 + DO istate = SIZE(dipole_op_mos_occ, 1), 1, -1 + CALL cp_fm_release(dipole_op_mos_occ(istate, ispin)%matrix) + END DO + END DO + DEALLOCATE (dipole_op_mos_occ) + END IF + DO ispin = nspins, 1, -1 CALL tddfpt_release_ground_state_mos(gs_mos(ispin)) END DO DEALLOCATE (gs_mos) - IF (tddfpt_control%dipole_form == tddfpt_dipole_length) & - CALL dbcsr_deallocate_matrix_set(rRc_xyz) - IF (tddfpt_control%dipole_form == tddfpt_dipole_berry) & - CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz) - CALL timestop(handle) END SUBROUTINE tddfpt @@ -2781,23 +2711,17 @@ CONTAINS END FUNCTION tddfpt_davidson_solver ! ************************************************************************************************** -!> \brief Print final TDDFPT excitation energies and analysis. -!> \param log_unit output unit -!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins; -!> SIZE(evects,2) -- number of excited states to print) -!> \param evals TDDFPT eigenvalues -!> \param mult multiplicity +!> \brief Compute the action of the dipole operator on the ground state wave function. +!> \param dipole_op_mos_occ 2-D array [x,y,z ; spin] of matrices where to put the computed quantity +!> (allocated and initialised on exit) +!> \param tddfpt_control TDDFPT control parameters !> \param gs_mos molecular orbitals optimised for the ground state -!> \param matrix_s overlap matrix -!> \param dipole_op_xyz components of the dipole operator in terms of contracted basis functions -!> \param dipole_form form of the dipole operator (berry-phase / length / vector) -!> \param cell unit cell -!> \param min_amplitude the smallest excitation amplitude to print -!> \param fm_pool_ao_mo_occ pools of dense matrices with shape [nao x nmo_occ(spin)] +!> \param qs_env Quickstep environment !> \par History -!> * 05.2016 created [Sergey Chulkov] -!> * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov] -!> * 07.2016 spin-unpolarised electron density [Sergey Chulkov] +!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov] +!> * 06.2018 dipole operator based on the Berry-phase formula [Sergey Chulkov] +!> * 08.2018 splited of from 'tddfpt_print_summary' and merged with code from 'tddfpt' +!> [Sergey Chulkov] !> \note \parblock !> Adapted version of the subroutine find_contributions() which was originally created !> by Thomas Chassaing on 02.2005. @@ -2809,29 +2733,481 @@ CONTAINS !> \f[\vec{r}\hat{H} - \hat{H}\vec{r} = [\vec{r},\hat{H}] = [\vec{r},-1/2 \nabla^2] = \nabla\f] . !> \endparblock ! ************************************************************************************************** - SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, mult, gs_mos, matrix_s, dipole_op_xyz, & - dipole_form, cell, min_amplitude, fm_pool_ao_mo_occ) + SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env) + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :), & + INTENT(inout) :: dipole_op_mos_occ + TYPE(tddfpt2_control_type), POINTER :: tddfpt_control + TYPE(tddfpt_ground_state_mos), DIMENSION(:), & + INTENT(in) :: gs_mos + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_dipole_operator', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, i_cos_sin, icol, ideriv, irow, & + ispin, jderiv, nao, ncols_local, & + ndim_periodic, nrows_local, nspins + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt + REAL(kind=dp) :: eval_occ + REAL(kind=dp), DIMENSION(3) :: kvec, reference_point + REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data_ediff, local_data_wfm + TYPE(cell_type), POINTER :: cell + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: gamma_00, gamma_inv_00 + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: dBerry_mos_occ, gamma_real_imag, opvec + TYPE(cp_fm_struct_type), POINTER :: fm_struct + TYPE(cp_fm_type), POINTER :: ediff_inv, rRc_mos_occ, wfm_ao_ao, & + wfm_mo_virt_mo_occ + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: berry_cossin_xyz, matrix_s, rRc_xyz + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_poisson_type), POINTER :: poisson_env + + CALL timeset(routineN, handle) + + NULLIFY (blacs_env, cell, matrix_s, pw_env) + CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, matrix_s=matrix_s, pw_env=pw_env) + + nspins = SIZE(gs_mos) + CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao) + DO ispin = 1, nspins + nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) + nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) + END DO + + ! +++ allocate dipole operator matrices (must be deallocated elsewhere) + ALLOCATE (dipole_op_mos_occ(nderivs, nspins)) + DO ispin = 1, nspins + CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) + + DO ideriv = 1, nderivs + NULLIFY (dipole_op_mos_occ(ideriv, ispin)%matrix) + CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin)%matrix, fm_struct) + END DO + END DO + + ! +++ allocate work matrices + ALLOCATE (S_mos_virt(nspins)) + DO ispin = 1, nspins + NULLIFY (S_mos_virt(ispin)%matrix) + CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct) + CALL cp_fm_create(S_mos_virt(ispin)%matrix, fm_struct) + CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, & + gs_mos(ispin)%mos_virt, & + S_mos_virt(ispin)%matrix, & + ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp) + END DO + + ! check that the chosen dipole operator is consistent with the periodic boundary conditions used + CALL pw_env_get(pw_env, poisson_env=poisson_env) + ndim_periodic = COUNT(poisson_env%parameters%periodic == 1) + + SELECT CASE (tddfpt_control%dipole_form) + CASE (tddfpt_dipole_berry) + IF (ndim_periodic /= 3) THEN + CALL cp_warn(__LOCATION__, & + "Fully periodic Poisson solver (PERIODIC xyz) is needed "// & + "for oscillator strengths based on the Berry phase formula") + END IF + + NULLIFY (berry_cossin_xyz) + ! index: 1 = Re[exp(-i * G_t * t)], + ! 2 = Im[exp(-i * G_t * t)]; + ! t = x,y,z + CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2) + + DO i_cos_sin = 1, 2 + CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix) + CALL dbcsr_copy(berry_cossin_xyz(i_cos_sin)%matrix, matrix_s(1)%matrix) + END DO + + ! +++ allocate berry-phase-related work matrices + ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins)) + ALLOCATE (dBerry_mos_occ(nderivs, nspins)) + DO ispin = 1, nspins + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) + + NULLIFY (gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) + CALL cp_cfm_create(gamma_00(ispin)%matrix, fm_struct) + CALL cp_cfm_create(gamma_inv_00(ispin)%matrix, fm_struct) + + DO i_cos_sin = 1, 2 + NULLIFY (gamma_real_imag(i_cos_sin, ispin)%matrix) + CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin)%matrix, fm_struct) + END DO + CALL cp_fm_struct_release(fm_struct) + + ! G_real C_0, G_imag C_0 + CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) + DO i_cos_sin = 1, 2 + NULLIFY (opvec(i_cos_sin, ispin)%matrix) + CALL cp_fm_create(opvec(i_cos_sin, ispin)%matrix, fm_struct) + END DO + + ! dBerry * C_0 + DO ideriv = 1, nderivs + NULLIFY (dBerry_mos_occ(ideriv, ispin)%matrix) + CALL cp_fm_create(dBerry_mos_occ(ideriv, ispin)%matrix, fm_struct) + CALL cp_fm_set_all(dBerry_mos_occ(ideriv, ispin)%matrix, 0.0_dp) + END DO + END DO + + DO ideriv = 1, nderivs + kvec(:) = twopi*cell%h_inv(ideriv, :) + DO i_cos_sin = 1, 2 + CALL dbcsr_set(berry_cossin_xyz(i_cos_sin)%matrix, 0.0_dp) + END DO + CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(1)%matrix, berry_cossin_xyz(2)%matrix, kvec) + + DO ispin = 1, nspins + ! i_cos_sin = 1: cos (real) component; opvec(1) = gamma_real C_0 + ! i_cos_sin = 2: sin (imaginary) component; opvec(2) = gamma_imag C_0 + DO i_cos_sin = 1, 2 + CALL cp_dbcsr_sm_fm_multiply(berry_cossin_xyz(i_cos_sin)%matrix, & + gs_mos(ispin)%mos_occ, & + opvec(i_cos_sin, ispin)%matrix, & + ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) + END DO + + CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, & + 1.0_dp, gs_mos(ispin)%mos_occ, opvec(1, ispin)%matrix, & + 0.0_dp, gamma_real_imag(1, ispin)%matrix) + + CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, & + -1.0_dp, gs_mos(ispin)%mos_occ, opvec(2, ispin)%matrix, & + 0.0_dp, gamma_real_imag(2, ispin)%matrix) + + CALL cp_fm_to_cfm(msourcer=gamma_real_imag(1, ispin)%matrix, & + msourcei=gamma_real_imag(2, ispin)%matrix, & + mtarget=gamma_00(ispin)%matrix) + + ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1} + CALL cp_cfm_set_all(gamma_inv_00(ispin)%matrix, z_zero, z_one) + CALL cp_cfm_solve(gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) + + CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin)%matrix, & + mtargetr=gamma_real_imag(1, ispin)%matrix, & + mtargeti=gamma_real_imag(2, ispin)%matrix) + + ! dBerry_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators() + CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), & + 1.0_dp, opvec(1, ispin)%matrix, gamma_real_imag(2, ispin)%matrix, & + 0.0_dp, dipole_op_mos_occ(1, ispin)%matrix) + CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), & + -1.0_dp, opvec(2, ispin)%matrix, gamma_real_imag(1, ispin)%matrix, & + 1.0_dp, dipole_op_mos_occ(1, ispin)%matrix) + + DO jderiv = 1, nderivs + CALL cp_fm_scale_and_add(1.0_dp, dBerry_mos_occ(jderiv, ispin)%matrix, & + cell%hmat(jderiv, ideriv), dipole_op_mos_occ(1, ispin)%matrix) + END DO + END DO + END DO + + ! --- release berry-phase-related work matrices + DO ispin = nspins, 1, -1 + DO i_cos_sin = SIZE(opvec, 1), 1, -1 + CALL cp_fm_release(opvec(i_cos_sin, ispin)%matrix) + END DO + + DO i_cos_sin = SIZE(gamma_real_imag, 1), 1, -1 + CALL cp_fm_release(gamma_real_imag(i_cos_sin, ispin)%matrix) + END DO + + CALL cp_cfm_release(gamma_inv_00(ispin)%matrix) + CALL cp_cfm_release(gamma_00(ispin)%matrix) + END DO + DEALLOCATE (gamma_00, gamma_inv_00, gamma_real_imag, opvec) + CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz) + + NULLIFY (fm_struct, wfm_ao_ao) + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) + CALL cp_fm_create(wfm_ao_ao, fm_struct) + CALL cp_fm_struct_release(fm_struct) + + ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) + + ! 2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) , + ! + ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact + ! that the response wave-function is a real-valued function, the above expression can be simplified as + ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) + ! + ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) . + DO ispin = 1, nspins + ! wfm_ao_ao = S * mos_virt * mos_virt^T + CALL cp_gemm('N', 'T', nao, nao, nmo_virt(ispin), & + 1.0_dp/twopi, S_mos_virt(ispin)%matrix, gs_mos(ispin)%mos_virt, & + 0.0_dp, wfm_ao_ao) + + DO ideriv = 1, nderivs + CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nao, & + 1.0_dp, wfm_ao_ao, dBerry_mos_occ(ideriv, ispin)%matrix, & + 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) + END DO + END DO + + CALL cp_fm_release(wfm_ao_ao) + + DO ispin = nspins, 1, -1 + DO ideriv = SIZE(dBerry_mos_occ, 1), 1, -1 + CALL cp_fm_release(dBerry_mos_occ(ideriv, ispin)%matrix) + END DO + END DO + DEALLOCATE (dBerry_mos_occ) + + CASE (tddfpt_dipole_length) + IF (ndim_periodic /= 0) THEN + CALL cp_warn(__LOCATION__, & + "Non-periodic Poisson solver (PERIODIC none) is needed "// & + "for oscillator strengths based on the length operator") + END IF + + ! compute components of the dipole operator in the length form + NULLIFY (rRc_xyz) + CALL dbcsr_allocate_matrix_set(rRc_xyz, nderivs) + + DO ideriv = 1, nderivs + CALL dbcsr_init_p(rRc_xyz(ideriv)%matrix) + CALL dbcsr_copy(rRc_xyz(ideriv)%matrix, matrix_s(1)%matrix) + END DO + + CALL get_reference_point(reference_point, qs_env=qs_env, & + reference=tddfpt_control%dipole_reference, & + ref_point=tddfpt_control%dipole_ref_point) + + CALL rRc_xyz_ao(op=rRc_xyz, qs_env=qs_env, rc=reference_point, order=1, minimum_image=.FALSE., soft=.FALSE.) + + NULLIFY (fm_struct, wfm_ao_ao) + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) + CALL cp_fm_create(wfm_ao_ao, fm_struct) + CALL cp_fm_struct_release(fm_struct) + + DO ispin = 1, nspins + NULLIFY (rRc_mos_occ) + CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct) + CALL cp_fm_create(rRc_mos_occ, fm_struct) + + ! wfm_ao_ao = S * mos_virt * mos_virt^T + CALL cp_gemm('N', 'T', nao, nao, nmo_virt(ispin), & + 1.0_dp, S_mos_virt(ispin)%matrix, gs_mos(ispin)%mos_virt, & + 0.0_dp, wfm_ao_ao) + + DO ideriv = 1, nderivs + CALL cp_dbcsr_sm_fm_multiply(rRc_xyz(ideriv)%matrix, & + gs_mos(ispin)%mos_occ, & + rRc_mos_occ, & + ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) + + CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nao, & + 1.0_dp, wfm_ao_ao, rRc_mos_occ, & + 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) + END DO + + CALL cp_fm_release(rRc_mos_occ) + END DO + + CALL cp_fm_release(wfm_ao_ao) + CALL dbcsr_deallocate_matrix_set(rRc_xyz) + + CASE (tddfpt_dipole_velocity) + IF (SIZE(matrix_s) < nderivs+1) THEN + CPABORT("Where is the derivative?") + END IF + + DO ispin = 1, nspins + NULLIFY (fm_struct, ediff_inv, wfm_mo_virt_mo_occ) + CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) + CALL cp_fm_create(ediff_inv, fm_struct) + CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct) + CALL cp_fm_struct_release(fm_struct) + + CALL cp_fm_get_info(ediff_inv, nrow_local=nrows_local, ncol_local=ncols_local, & + row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff) + CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm) + +!$OMP PARALLEL DO DEFAULT(NONE), & +!$OMP PRIVATE(eval_occ, icol, irow), & +!$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices) + DO icol = 1, ncols_local + ! E_occ_i ; imo_occ = col_indices(icol) + eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol)) + + DO irow = 1, nrows_local + ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i) + ! imo_virt = row_indices(irow) + local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow))-eval_occ) + END DO + END DO +!$OMP END PARALLEL DO + + DO ideriv = 1, nderivs + CALL cp_dbcsr_sm_fm_multiply(matrix_s(ideriv+1)%matrix, & + gs_mos(ispin)%mos_occ, & + dipole_op_mos_occ(ideriv, ispin)%matrix, & + ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) + + CALL cp_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, & + 1.0_dp, gs_mos(ispin)%mos_virt, dipole_op_mos_occ(ideriv, ispin)%matrix, & + 0.0_dp, wfm_mo_virt_mo_occ) + + ! in-place element-wise (Schur) product; + ! avoid allocation of a temporary [nmo_virt x nmo_occ] matrix which is needed for cp_fm_schur_product() subroutine call + +!$OMP PARALLEL DO DEFAULT(NONE), & +!$OMP PRIVATE(icol, irow), & +!$OMP SHARED(ispin, local_data_ediff, local_data_wfm, ncols_local, nrows_local) + DO icol = 1, ncols_local + DO irow = 1, nrows_local + local_data_wfm(irow, icol) = local_data_wfm(irow, icol)*local_data_ediff(irow, icol) + END DO + END DO +!$OMP END PARALLEL DO + + CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), & + 1.0_dp, S_mos_virt(ispin)%matrix, wfm_mo_virt_mo_occ, & + 0.0_dp, dipole_op_mos_occ(ideriv, ispin)%matrix) + END DO + + CALL cp_fm_release(wfm_mo_virt_mo_occ) + CALL cp_fm_release(ediff_inv) + END DO + + CASE DEFAULT + CPABORT("Unimplemented form of the dipole operator") + END SELECT + + ! --- release work matrices + DO ispin = nspins, 1, -1 + CALL cp_fm_release(S_mos_virt(ispin)%matrix) + END DO + DEALLOCATE (S_mos_virt) + + CALL timestop(handle) + END SUBROUTINE tddfpt_dipole_operator + +! ************************************************************************************************** +!> \brief Print final TDDFPT excitation energies and oscillator strengths. +!> \param log_unit output unit +!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins; +!> SIZE(evects,2) -- number of excited states to print) +!> \param evals TDDFPT eigenvalues +!> \param mult multiplicity +!> \param dipole_op_mos_occ action of the dipole operator on the ground state wave function +!> [x,y,z ; spin] +!> \par History +!> * 05.2016 created [Sergey Chulkov] +!> * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov] +!> * 07.2016 spin-unpolarised electron density [Sergey Chulkov] +!> * 08.2018 compute 'dipole_op_mos_occ' in a separate subroutine [Sergey Chulkov] +!> \note \parblock +!> Adapted version of the subroutine find_contributions() which was originally created +!> by Thomas Chassaing on 02.2005. +!> +!> Transition dipole moment along direction 'd' is computed as following: +!> \f[ t_d(spin) = Tr[evects^T dipole\_op\_mos\_occ(d, spin)] .\f] +!> \endparblock +! ************************************************************************************************** + SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, mult, dipole_op_mos_occ) INTEGER, INTENT(in) :: log_unit TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals INTEGER, INTENT(in) :: mult - TYPE(tddfpt_ground_state_mos), DIMENSION(:), & - INTENT(in) :: gs_mos - TYPE(dbcsr_type), POINTER :: matrix_s - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipole_op_xyz - INTEGER, INTENT(in) :: dipole_form - TYPE(cell_type), POINTER :: cell - REAL(kind=dp), INTENT(in) :: min_amplitude - TYPE(cp_fm_pool_p_type), DIMENSION(:), INTENT(in) :: fm_pool_ao_mo_occ + TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: dipole_op_mos_occ CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_summary', & routineP = moduleN//':'//routineN CHARACTER(len=1) :: lsd_str CHARACTER(len=20) :: mult_str + INTEGER :: handle, ideriv, ispin, istate, nspins, & + nstates + REAL(kind=dp) :: osc_strength + REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_dipoles + + CALL timeset(routineN, handle) + + nspins = SIZE(evects, 1) + nstates = SIZE(evects, 2) + + IF (nspins > 1) THEN + lsd_str = 'U' + ELSE + lsd_str = 'R' + END IF + + ! *** summary header *** + IF (log_unit > 0) THEN + CALL integer_to_string(mult, mult_str) + WRITE (log_unit, '(/,1X,A1,A,1X,A,/)') lsd_str, "-TDDFPT states of multiplicity", TRIM(mult_str) + + WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", "Transition dipole (a.u.)", "Oscillator" + WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", "x", "y", "z", "strength (a.u.)" + WRITE (log_unit, '(T10,72("-"))') + END IF + + ! transition dipole moment + ALLOCATE (trans_dipoles(nstates, nderivs, nspins)) + trans_dipoles(:, :, :) = 0.0_dp + + ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons + IF (nspins > 1 .OR. mult == 1) THEN + DO ispin = 1, nspins + DO ideriv = 1, nderivs + CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin)%matrix, trans_dipoles(:, ideriv, ispin)) + END DO + END DO + + IF (nspins == 1) THEN + trans_dipoles(:, :, 1) = SQRT(2.0_dp)*trans_dipoles(:, :, 1) + ELSE + trans_dipoles(:, :, 1) = SQRT(trans_dipoles(:, :, 1)**2+trans_dipoles(:, :, 2)**2) + END IF + END IF + + ! *** summary information *** + DO istate = 1, nstates + IF (log_unit > 0) THEN + osc_strength = 2.0_dp/3.0_dp*evals(istate)* & + accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1)) + + WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') & + "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength + END IF + END DO + + DEALLOCATE (trans_dipoles) + + CALL timestop(handle) + END SUBROUTINE tddfpt_print_summary + +! ************************************************************************************************** +!> \brief Print excitation analysis. +!> \param log_unit output unit +!> \param evects TDDFPT trial vectors (SIZE(evects,1) -- number of spins; +!> SIZE(evects,2) -- number of excited states to print) +!> \param gs_mos molecular orbitals optimised for the ground state +!> \param matrix_s overlap matrix +!> \param min_amplitude the smallest excitation amplitude to print +!> \par History +!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov] +!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov] +! ************************************************************************************************** + SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, gs_mos, matrix_s, min_amplitude) + INTEGER, INTENT(in) :: log_unit + TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects + TYPE(tddfpt_ground_state_mos), DIMENSION(:), & + INTENT(in) :: gs_mos + TYPE(dbcsr_type), POINTER :: matrix_s + REAL(kind=dp), INTENT(in) :: min_amplitude + + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_excitation_analysis', & + routineP = moduleN//':'//routineN + CHARACTER(len=5) :: spin_label - INTEGER :: handle, icol, ideriv, iproc, irow, ispin, istate, jderiv, nao, ncols_local, & - nrows_local, nspins, nstates, send_handler, send_handler2, state_spin + INTEGER :: handle, icol, iproc, irow, ispin, istate, nao, ncols_local, nrows_local, nspins, & + nstates, send_handler, send_handler2, state_spin INTEGER(kind=int_8) :: iexc, imo_occ, imo_virt, ind, nexcs, & nexcs_local, nexcs_max_local, & nmo_virt_occ, nmo_virt_occ_alpha @@ -2842,17 +3218,11 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt LOGICAL :: do_exc_analysis - REAL(kind=dp) :: eval_occ, osc_strength, trans_dipole REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights_local, weights_neg_abs_recv, & weights_recv REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data - REAL(kind=dp), DIMENSION(nderivs) :: trans_dipoles TYPE(cp_blacs_env_type), POINTER :: blacs_env - TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: gamma_00, gamma_inv_00 - TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: ediff_inv, ediff_inv_weights, & - gamma_imag, gamma_real, S_mos_virt, & - scaled_evect, weights_fm - TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: dS_mos_occ, opvec + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt, weights_fm TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_para_env_type), POINTER :: para_env @@ -2861,16 +3231,9 @@ CONTAINS nspins = SIZE(evects, 1) nstates = SIZE(evects, 2) do_exc_analysis = min_amplitude < 1.0_dp - osc_strength = 0.0_dp + CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env) CALL dbcsr_get_info(matrix_s, nfullrows_total=nao) - CALL cp_fm_get_info(gs_mos(1)%mos_virt, context=blacs_env, para_env=para_env) - - IF (nspins > 1) THEN - lsd_str = 'U' - ELSE - lsd_str = 'R' - END IF DO ispin = 1, nspins nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) @@ -2879,226 +3242,6 @@ CONTAINS nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8) END DO - ALLOCATE (dS_mos_occ(nderivs, nspins), S_mos_virt(nspins), weights_fm(nspins)) - DO ispin = 1, nspins - DO ideriv = 1, nderivs - NULLIFY (dS_mos_occ(ideriv, ispin)%matrix) - CALL fm_pool_create_fm(fm_pool_ao_mo_occ(ispin)%pool, dS_mos_occ(ideriv, ispin)%matrix) - END DO - - NULLIFY (S_mos_virt(ispin)%matrix) - CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct) - CALL cp_fm_create(S_mos_virt(ispin)%matrix, fm_struct) - - NULLIFY (fm_struct, weights_fm(ispin)%matrix) - CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) - CALL cp_fm_create(weights_fm(ispin)%matrix, fm_struct) - CALL cp_fm_struct_release(fm_struct) - END DO - - IF (dipole_form == tddfpt_dipole_berry) THEN - ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_imag(nspins), gamma_real(nspins), opvec(2, nspins)) - NULLIFY (fm_struct) - - DO ispin = 1, nspins - CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) - - NULLIFY (gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) - NULLIFY (gamma_imag(ispin)%matrix, gamma_real(ispin)%matrix) - - CALL cp_cfm_create(gamma_00(ispin)%matrix, fm_struct) - CALL cp_cfm_create(gamma_inv_00(ispin)%matrix, fm_struct) - CALL cp_fm_create(gamma_imag(ispin)%matrix, fm_struct) - CALL cp_fm_create(gamma_real(ispin)%matrix, fm_struct) - - CALL cp_fm_struct_release(fm_struct) - - DO ideriv = 1, 2 ! G_real C_0, G_imag C_0 - NULLIFY (opvec(ideriv, ispin)%matrix) - CALL fm_pool_create_fm(fm_pool_ao_mo_occ(ispin)%pool, opvec(ideriv, ispin)%matrix) - END DO - END DO - - DO ispin = 1, nspins - DO ideriv = 1, nderivs - ! cos (real) component; opvec(1) = gamma_real C_0 - CALL cp_dbcsr_sm_fm_multiply(dipole_op_xyz(2*ideriv-1)%matrix, gs_mos(ispin)%mos_occ, & - opvec(1, ispin)%matrix, ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) - CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_occ, & - opvec(1, ispin)%matrix, 0.0_dp, gamma_real(ispin)%matrix) - - ! sin (imaginary) component; opvec(2) = gamma_imag C_0 - CALL cp_dbcsr_sm_fm_multiply(dipole_op_xyz(2*ideriv)%matrix, gs_mos(ispin)%mos_occ, & - opvec(2, ispin)%matrix, ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) - CALL cp_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, -1.0_dp, gs_mos(ispin)%mos_occ, & - opvec(2, ispin)%matrix, 0.0_dp, gamma_imag(ispin)%matrix) - - CALL cp_fm_to_cfm(msourcer=gamma_real(ispin)%matrix, & - msourcei=gamma_imag(ispin)%matrix, & - mtarget=gamma_00(ispin)%matrix) - - ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1} - CALL cp_cfm_set_all(gamma_inv_00(ispin)%matrix, z_zero, z_one) - CALL cp_cfm_solve(gamma_00(ispin)%matrix, gamma_inv_00(ispin)%matrix) - - CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin)%matrix, & - mtargetr=gamma_real(ispin)%matrix, & - mtargeti=gamma_imag(ispin)%matrix) - - ! dS_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators() - CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), 1.0_dp, opvec(1, ispin)%matrix, & - gamma_imag(ispin)%matrix, 0.0_dp, dS_mos_occ(ideriv, ispin)%matrix) - CALL cp_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), -1.0_dp, opvec(2, ispin)%matrix, & - gamma_real(ispin)%matrix, 1.0_dp, dS_mos_occ(ideriv, ispin)%matrix) - END DO - END DO - - DO ispin = nspins, 1, -1 - DO ideriv = 2, 1, -1 - CALL fm_pool_give_back_fm(fm_pool_ao_mo_occ(ispin)%pool, opvec(ideriv, ispin)%matrix) - END DO - - CALL cp_fm_release(gamma_real(ispin)%matrix) - CALL cp_fm_release(gamma_imag(ispin)%matrix) - CALL cp_cfm_release(gamma_inv_00(ispin)%matrix) - CALL cp_cfm_release(gamma_00(ispin)%matrix) - END DO - DEALLOCATE (gamma_00, gamma_inv_00, gamma_imag, gamma_real, opvec) - END IF - - IF (dipole_form == tddfpt_dipole_velocity) THEN - ALLOCATE (ediff_inv(nspins), ediff_inv_weights(nspins), scaled_evect(nspins)) - DO ispin = 1, nspins - NULLIFY (ediff_inv(ispin)%matrix, ediff_inv_weights(ispin)%matrix, scaled_evect(ispin)%matrix) - CALL cp_fm_get_info(weights_fm(ispin)%matrix, matrix_struct=fm_struct) - - CALL cp_fm_create(ediff_inv(ispin)%matrix, fm_struct) - CALL cp_fm_create(ediff_inv_weights(ispin)%matrix, fm_struct) - - CALL fm_pool_create_fm(fm_pool_ao_mo_occ(ispin)%pool, scaled_evect(ispin)%matrix) - END DO - - ! compute (E_virt_a - E_occ_i)^{-1} - DO ispin = 1, nspins - CALL cp_fm_get_info(ediff_inv(ispin)%matrix, nrow_local=nrows_local, ncol_local=ncols_local, & - row_indices=row_indices, col_indices=col_indices, local_data=local_data) - -!$OMP PARALLEL DO DEFAULT(NONE), & -!$OMP PRIVATE(eval_occ, icol, irow), & -!$OMP SHARED(col_indices, gs_mos, ispin, local_data, ncols_local, nrows_local, row_indices) - DO icol = 1, ncols_local - ! E_occ_i ; imo_occ = col_indices(icol) - eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol)) - - DO irow = 1, nrows_local - ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i) - ! imo_virt = row_indices(irow) - local_data(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow))-eval_occ) - END DO - END DO -!$OMP END PARALLEL DO - END DO - END IF - - DO ispin = 1, nspins - CALL cp_dbcsr_sm_fm_multiply(matrix_s, gs_mos(ispin)%mos_virt, S_mos_virt(ispin)%matrix, & - ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp) - - IF (dipole_form == tddfpt_dipole_length .OR. dipole_form == tddfpt_dipole_velocity) THEN - DO ideriv = 1, nderivs - CALL cp_dbcsr_sm_fm_multiply(dipole_op_xyz(ideriv)%matrix, gs_mos(ispin)%mos_occ, & - dS_mos_occ(ideriv, ispin)%matrix, ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) - END DO - END IF - END DO - - ! *** summary header *** - IF (log_unit > 0) THEN - CALL integer_to_string(mult, mult_str) - WRITE (log_unit, '(/,1X,A1,A,1X,A,/)') lsd_str, "-TDDFPT states of multiplicity", TRIM(mult_str) - - WRITE (log_unit, '(T11,A,T21,A,T37,A,T69,A)') "State", "Excitation", "Transition dipole (a.u.)", "Oscillator" - WRITE (log_unit, '(T11,A,T21,A,T39,A,T49,A,T59,A,T67,A)') "number", "energy (eV)", "x", "y", "z", "strength (a.u.)" - WRITE (log_unit, '(T10,72("-"))') - END IF - - DO istate = 1, nstates - trans_dipoles(:) = 0.0_dp - - ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons - IF (nspins > 1 .OR. mult == 1) THEN - DO ispin = 1, nspins - CALL cp_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, 1.0_dp, S_mos_virt(ispin)%matrix, & - evects(ispin, istate)%matrix, 0.0_dp, weights_fm(ispin)%matrix) - - ! compute electric dipole moments - IF (dipole_form == tddfpt_dipole_velocity) THEN - CALL cp_fm_schur_product(weights_fm(ispin)%matrix, ediff_inv(ispin)%matrix, ediff_inv_weights(ispin)%matrix) - CALL cp_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, & - ediff_inv_weights(ispin)%matrix, 0.0_dp, scaled_evect(ispin)%matrix) - END IF - - DO ideriv = 1, nderivs - SELECT CASE (dipole_form) - CASE (tddfpt_dipole_berry) - ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) + - ! 2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) , - ! - ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact - ! that the response wave-function is a real-valued function, the above expression can be simplified as - ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) - CALL cp_fm_trace(dS_mos_occ(ideriv, ispin)%matrix, evects(ispin, istate)%matrix, trans_dipole) - - ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) . - ! In case of an isolated molecule in a large unit cell, the scale factor 1.0_dp/(3.0_dp*pi) - ! gives the same oscillator strengths as the ones obtained using the length operator. - trans_dipole = 1.0_dp/(3.0_dp*pi)*trans_dipole - DO jderiv = 1, nderivs - trans_dipoles(jderiv) = trans_dipoles(jderiv)+cell%hmat(jderiv, ideriv)*trans_dipole - END DO - - CASE (tddfpt_dipole_length) - CALL cp_fm_trace(dS_mos_occ(ideriv, ispin)%matrix, evects(ispin, istate)%matrix, trans_dipole) - trans_dipoles(ideriv) = trans_dipoles(ideriv)+trans_dipole - - CASE (tddfpt_dipole_velocity) - CALL cp_fm_trace(dS_mos_occ(ideriv, ispin)%matrix, scaled_evect(ispin)%matrix, trans_dipole) - trans_dipoles(ideriv) = trans_dipoles(ideriv)+trans_dipole - - CASE DEFAULT - CPABORT("Unimplemented form of the dipole operator") - END SELECT - END DO - END DO - END IF - - ! *** summary information *** - IF (log_unit > 0) THEN - IF (nspins == 1 .AND. dipole_form /= tddfpt_dipole_velocity) & - trans_dipoles(:) = SQRT(2.0_dp)*trans_dipoles(:) - osc_strength = 2.0_dp/3.0_dp*evals(istate)*accurate_dot_product(trans_dipoles, trans_dipoles) - - WRITE (log_unit, '(1X,A,T9,I8,T21,F11.5,T34,3(1X,F9.4),T69,F10.5)') & - "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(1:nderivs), osc_strength - END IF - END DO - - IF (dipole_form == tddfpt_dipole_velocity) THEN - DO ispin = nspins, 1, -1 - CALL fm_pool_give_back_fm(fm_pool_ao_mo_occ(ispin)%pool, scaled_evect(ispin)%matrix) - CALL cp_fm_release(ediff_inv_weights(ispin)%matrix) - CALL cp_fm_release(ediff_inv(ispin)%matrix) - END DO - DEALLOCATE (ediff_inv, ediff_inv_weights, scaled_evect) - END IF - - DO ispin = nspins, 1, -1 - DO ideriv = nderivs, 1, -1 - CALL fm_pool_give_back_fm(fm_pool_ao_mo_occ(ispin)%pool, dS_mos_occ(ideriv, ispin)%matrix) - END DO - END DO - DEALLOCATE (dS_mos_occ) - ! *** excitation analysis *** IF (do_exc_analysis) THEN CPASSERT(log_unit <= 0 .OR. para_env%mepos == para_env%source) @@ -3117,6 +3260,22 @@ CONTAINS END IF END IF + ALLOCATE (S_mos_virt(nspins), weights_fm(nspins)) + DO ispin = 1, nspins + NULLIFY (S_mos_virt(ispin)%matrix) + CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct) + CALL cp_fm_create(S_mos_virt(ispin)%matrix, fm_struct) + CALL cp_dbcsr_sm_fm_multiply(matrix_s, & + gs_mos(ispin)%mos_virt, & + S_mos_virt(ispin)%matrix, & + ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp) + + NULLIFY (fm_struct, weights_fm(ispin)%matrix) + CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), context=blacs_env) + CALL cp_fm_create(weights_fm(ispin)%matrix, fm_struct) + CALL cp_fm_struct_release(fm_struct) + END DO + nexcs_max_local = 0 DO ispin = 1, nspins CALL cp_fm_get_info(weights_fm(ispin)%matrix, nrow_local=nrows_local, ncol_local=ncols_local) @@ -3274,7 +3433,7 @@ CONTAINS DEALLOCATE (S_mos_virt, weights_fm) CALL timestop(handle) - END SUBROUTINE tddfpt_print_summary + END SUBROUTINE tddfpt_print_excitation_analysis ! ************************************************************************************************** !> \brief Write Ritz vectors to a binary restart file.