diff --git a/data/EMSL_BASIS_SETS b/data/EMSL_BASIS_SETS index ce319143fc..075418a24b 100644 --- a/data/EMSL_BASIS_SETS +++ b/data/EMSL_BASIS_SETS @@ -2858,8 +2858,8 @@ S 6-311Gxx 6-311G** 902.45600000 0.11902800 0.00000000 297.15800000 0.36843200 0.00000000 108.70200000 0.57729900 0.14318600 - 43.15530000 0.00000000 0.62446500 - 18.10790000 0.00000000 0.28336600 + 43.155300000 0.00000000 0.62446500 + 18.107900000 0.00000000 0.28336600 1 0 0 1 1 5.56009000 1.00000000 1 0 0 1 1 diff --git a/src/input_cp2k_exstate.F b/src/input_cp2k_exstate.F index 70c3e1de8b..f4b2af4848 100644 --- a/src/input_cp2k_exstate.F +++ b/src/input_cp2k_exstate.F @@ -74,6 +74,14 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="OVERLAP_DELTAT", & + description="Keyword for the computation of the overlap matrix between two consecutive time steps.", & + usage="OVERLAP_DELTAT", & + default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_exstate_section END MODULE input_cp2k_exstate diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index 1e55ab52e8..c80aa01292 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -1600,6 +1600,16 @@ CONTAINS default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PRINT_PHASES", & + description="Print phases of occupied and virtuals MOs.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="SCALE_WITH_PHASES", & + description="Scale ES eigenvectors with phases of occupied and virtuals MOs.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) CALL section_add_subsection(subsection, print_key) CALL section_release(print_key) diff --git a/src/motion/input_cp2k_vib.F b/src/motion/input_cp2k_vib.F index ff4d60e449..3b8a6b942d 100644 --- a/src/motion/input_cp2k_vib.F +++ b/src/motion/input_cp2k_vib.F @@ -134,12 +134,6 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="NAMD_PRINT", & - description="Adjust cartesian eigenvalues / vectors to NewtonX format.", & - default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - CALL create_mode_selective_section(subsection) CALL section_add_subsection(section, subsection) CALL section_release(subsection) @@ -217,6 +211,19 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, name="NAMD_PRINT", & + description="Adjust cartesian eigenvalues / vectors to NewtonX format.", & + print_level=debug_print_level + 1, add_last=add_last_numeric, & + filename="FullNormalizedCartesian") + CALL keyword_create(keyword, __LOCATION__, name="BACKUP_COPIES", & + description="Specifies the maximum index of backup copies.", & + usage="BACKUP_COPIES {int}", & + default_i_val=3) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "HESSIAN", & description="Write the Hessian matrix from a vibrational analysis calculation "// & "into a binary file.", & diff --git a/src/motion/vibrational_analysis.F b/src/motion/vibrational_analysis.F index d72ac0bbbd..be64201a8a 100644 --- a/src/motion/vibrational_analysis.F +++ b/src/motion/vibrational_analysis.F @@ -110,10 +110,10 @@ CONTAINS CHARACTER(LEN=default_string_length) :: description_d, description_p INTEGER :: handle, i, icoord, icoordm, icoordp, ierr, imap, iounit, ip1, ip2, iparticle1, & iparticle2, iseq, iw, j, k, natoms, ncoord, nfrozen, nrep, nres, nRotTrM, nvib, & - output_unit, output_unit_eig, prep, print_grrm, print_scine, proc_dist_type + output_unit, output_unit_eig, prep, print_grrm, print_namd, print_scine, proc_dist_type INTEGER, DIMENSION(:), POINTER :: Clist, Mlist LOGICAL :: calc_intens, calc_thchdata, do_mode_tracking, intens_ir, intens_raman, & - keep_rotations, lprint_namd, row_force, something_frozen + keep_rotations, row_force, something_frozen REAL(KIND=dp) :: a1, a2, a3, conver, dummy, dx, & inertia(3), minimum_energy, norm, & tc_press, tc_temp, tmp @@ -167,7 +167,7 @@ CONTAINS CALL section_vals_val_get(vib_section, "THERMOCHEMISTRY", l_val=calc_thchdata) CALL section_vals_val_get(vib_section, "TC_TEMPERATURE", r_val=tc_temp) CALL section_vals_val_get(vib_section, "TC_PRESSURE", r_val=tc_press) - CALL section_vals_val_get(vib_section, "NAMD_PRINT", l_val=lprint_namd) + tc_temp = tc_temp*kelvin tc_press = tc_press*pascal @@ -520,6 +520,52 @@ CONTAINS END IF CALL cp_print_key_finished_output(print_scine, logger, force_env_section, "PRINT%SCINE") ! + ! Print NEWTONX interface file + print_namd = cp_print_key_unit_nr(logger, print_section, "NAMD_PRINT", & + extension=".eig", file_status="REPLACE", & + file_action="WRITE", do_backup=.TRUE., & + file_form="UNFORMATTED") + IF (print_namd > 0) THEN + ! NewtonX requires normalized Cartesian frequencies and eigenvectors + ! in full matrix format (ncoord x ncoord) + NULLIFY (Dfull) + ALLOCATE (Dfull(ncoord, ncoord)) + ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2))) + ALLOCATE (HeigvalDfull(SIZE(Dfull, 2))) + ALLOCATE (MatM(ncoord, ncoord)) + ALLOCATE (rmass(SIZE(Dfull, 2))) + Dfull = 0.0_dp + ! Dfull in dimension of degrees of freedom + CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms) + ! TEST MatM = MATMUL(TRANSPOSE(Dfull),Dfull)= 1 + ! Hessian in MWC -> Hessian in INT (Hint2Dfull) + Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull)) + ! Heig = L^T Hint2Dfull L + CALL diamat_all(Hint2Dfull, HeigvalDfull) + ! TEST MatM = MATMUL(TRANSPOSE(Hint2Dfull),Hint2Dfull) = 1 + ! TEST MatM=MATMUL(TRANSPOSE(MATMUL(Dfull,Hint2Dfull)),MATMUL(Dfull,Hint2Dfull)) = 1 + MatM = 0.0_dp + DO i = 1, natoms + DO j = 1, 3 + MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i) ! mass is sqrt(mass) + END DO + END DO + ! Dfull = Cartesian displacements of the normal modes + Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !Dfull=D L / sqrt(m) + DO i = 1, ncoord + ! Renormalize displacements + norm = 1.0_dp/SUM(Dfull(:, i)*Dfull(:, i)) + rmass(i) = norm/massunit + Dfull(:, i) = SQRT(norm)*(Dfull(:, i)) + END DO + CALL write_eigs_unformatted(print_namd, ncoord, HeigvalDfull, Dfull) + DEALLOCATE (HeigvalDfull) + DEALLOCATE (Hint2Dfull) + DEALLOCATE (Dfull) + DEALLOCATE (MatM) + DEALLOCATE (rmass) + END IF !print_namd + ! nvib = ncoord - nRotTrM ALLOCATE (H_eigval1(ncoord)) ALLOCATE (H_eigval2(SIZE(D, 2))) @@ -552,35 +598,7 @@ CONTAINS END IF ! write frequencies and eigenvectors to cartesian eig file IF (output_unit_eig > 0) THEN - IF (lprint_namd) THEN - ! NewtonX requires normalized Cartesian frequencies and eigenvectors - ! in full matrix format (ncoord x ncoord) - NULLIFY (Dfull) - ALLOCATE (Dfull(ncoord, ncoord)) - ALLOCATE (Hint2Dfull(SIZE(Dfull, 2), SIZE(Dfull, 2))) - ALLOCATE (HeigvalDfull(SIZE(Dfull, 2))) - ALLOCATE (MatM(ncoord, ncoord)) - Dfull = 0.0_dp - CALL build_D_matrix(RotTrM, nRotTrM, Dfull, full=.TRUE., natoms=natoms) - Hint2Dfull(:, :) = MATMUL(TRANSPOSE(Dfull), MATMUL(Hessian, Dfull)) ! MWC -> INT - CALL diamat_all(Hint2Dfull, HeigvalDfull) - MatM = 0.0_dp - DO i = 1, natoms - DO j = 1, 3 - MatM((i - 1)*3 + j, (i - 1)*3 + j) = 1.0_dp/mass(i) - END DO - END DO - ! Cartesian displacements of the normal modes - Dfull = MATMUL(MatM, MATMUL(Dfull, Hint2Dfull)) !!! this is !D L / sqrt(m) - ! CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Dfull) - CALL write_eigs_unformatted(output_unit_eig, ncoord, HeigvalDfull, Dfull) - DEALLOCATE (HeigvalDfull) - DEALLOCATE (Hint2Dfull) - DEALLOCATE (Dfull) - DEALLOCATE (MatM) - END IF - CALL write_eigs_unformatted(output_unit_eig, & - ncoord, H_eigval1, Hint1) + CALL write_eigs_unformatted(output_unit_eig, ncoord, H_eigval1, Hint1) END IF IF (nvib /= 0) THEN Hint2(:, :) = MATMUL(TRANSPOSE(D), MATMUL(Hessian, D)) diff --git a/src/qs_energy.F b/src/qs_energy.F index dc534c8327..601bf093b6 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -18,6 +18,10 @@ MODULE qs_energy USE dm_ls_scf, ONLY: ls_scf USE energy_corrections, ONLY: energy_correction USE excited_states, ONLY: excited_state_energy + USE input_section_types, ONLY: section_vals_get,& + section_vals_get_subs_vals,& + section_vals_type,& + section_vals_val_get USE lri_environment_methods, ONLY: lri_print_stat USE mp2, ONLY: mp2_main USE qs_energy_init, ONLY: qs_energies_init @@ -63,10 +67,12 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies' INTEGER :: handle - LOGICAL :: do_consistent_energies, my_calc_forces, & - run_rtp + LOGICAL :: do_consistent_energies, & + do_excited_state, loverlap_deltat, & + my_calc_forces, run_rtp TYPE(dft_control_type), POINTER :: dft_control TYPE(qs_energy_type), POINTER :: energy + TYPE(section_vals_type), POINTER :: excited_state_section CALL timeset(routineN, handle) @@ -85,13 +91,25 @@ CONTAINS CALL qs_energies_init(qs_env, my_calc_forces) CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, energy=energy) + ! *** check if only overlap matrix is needed for couplings + loverlap_deltat = .FALSE. + NULLIFY (excited_state_section) + excited_state_section => section_vals_get_subs_vals(qs_env%input, "DFT%EXCITED_STATES") + CALL section_vals_get(excited_state_section, explicit=do_excited_state) + IF (do_excited_state) THEN + CALL section_vals_val_get(excited_state_section, "OVERLAP_DELTAT", & + l_val=loverlap_deltat) + END IF + ! *** Perform a SCF run *** - IF (dft_control%qs_control%do_ls_scf) THEN - CALL ls_scf(qs_env=qs_env) - ELSE IF (dft_control%qs_control%do_almo_scf) THEN - CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces) - ELSE - CALL scf(qs_env=qs_env) + IF (.NOT. loverlap_deltat) THEN + IF (dft_control%qs_control%do_ls_scf) THEN + CALL ls_scf(qs_env=qs_env) + ELSE IF (dft_control%qs_control%do_almo_scf) THEN + CALL almo_entry_scf(qs_env=qs_env, calc_forces=my_calc_forces) + ELSE + CALL scf(qs_env=qs_env) + END IF END IF IF (do_consistent_energies) THEN @@ -114,9 +132,11 @@ CONTAINS CALL energy_correction(qs_env, ec_init=.TRUE., calculate_forces=.FALSE.) END IF - CALL qs_energies_properties(qs_env, calc_forces) + IF (.NOT. loverlap_deltat) THEN + CALL qs_energies_properties(qs_env, calc_forces) - CALL excited_state_energy(qs_env, calculate_forces=.FALSE.) + CALL excited_state_energy(qs_env, calculate_forces=.FALSE.) + END IF IF (dft_control%qs_control%lrigpw) THEN CALL lri_print_stat(qs_env) diff --git a/src/qs_tddfpt2_methods.F b/src/qs_tddfpt2_methods.F index 08b69ad83a..9cc2b2f688 100644 --- a/src/qs_tddfpt2_methods.F +++ b/src/qs_tddfpt2_methods.F @@ -143,7 +143,7 @@ CONTAINS old_state INTEGER, DIMENSION(maxspins) :: nactive LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, & - is_restarted, state_change + explicit, is_restarted, state_change REAL(kind=dp) :: conv REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength TYPE(admm_type), POINTER :: admm_env @@ -159,8 +159,9 @@ CONTAINS TYPE(kernel_env_type) :: kernel_env TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(qs_scf_env_type), POINTER :: scf_env - TYPE(section_vals_type), POINTER :: lri_section, tddfpt_print_section, & - tddfpt_section, xc_section + TYPE(section_vals_type), POINTER :: lri_section, namd_print_section, & + tddfpt_print_section, tddfpt_section, & + xc_section TYPE(stda_env_type), TARGET :: stda_kernel TYPE(tddfpt2_control_type), POINTER :: tddfpt_control TYPE(tddfpt_ground_state_mos), DIMENSION(:), & @@ -433,8 +434,6 @@ CONTAINS CALL cp_warn(__LOCATION__, "Skipping TDDFPT wavefunction optimization") END IF - CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_env) - IF (ASSOCIATED(matrix_ks_oep) .AND. tddfpt_control%dipole_form == tddfpt_dipole_velocity) THEN CALL cp_warn(__LOCATION__, & "Transition dipole moments and oscillator strengths are likely to be incorrect "// & @@ -444,6 +443,13 @@ CONTAINS ! *** print summary information *** log_unit = cp_logger_get_default_io_unit() + + namd_print_section => section_vals_get_subs_vals(tddfpt_print_section, "NAMD_PRINT") + CALL section_vals_get(namd_print_section, explicit=explicit) + IF (explicit) THEN + CALL tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, & + matrix_s(1)%matrix, S_evects, sub_env) !qs_env, sub_env) + END IF ALLOCATE (ostrength(nstates)) ostrength = 0.0_dp CALL tddfpt_print_summary(log_unit, evects, evals, ostrength, mult, & diff --git a/src/qs_tddfpt2_restart.F b/src/qs_tddfpt2_restart.F index 4b15dadb35..4534018ff1 100644 --- a/src/qs_tddfpt2_restart.F +++ b/src/qs_tddfpt2_restart.F @@ -8,17 +8,20 @@ MODULE qs_tddfpt2_restart USE cp_blacs_env, ONLY: cp_blacs_env_type,& get_blacs_info + USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply USE cp_files, ONLY: close_file,& open_file - USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale + USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& + cp_fm_scale_and_add,& + cp_fm_trace USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& fm_pool_create_fm - USE cp_fm_types, ONLY: cp_fm_get_info,& - cp_fm_p_type,& - cp_fm_read_unformatted,& - cp_fm_write_formatted,& - cp_fm_write_info,& - cp_fm_write_unformatted + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: & + cp_fm_create, cp_fm_get_info, cp_fm_p_type, cp_fm_read_unformatted, cp_fm_release, & + cp_fm_type, cp_fm_write_formatted, cp_fm_write_info, cp_fm_write_unformatted USE cp_log_handling, ONLY: cp_logger_type USE cp_output_handling, ONLY: cp_p_file,& cp_print_key_finished_output,& @@ -26,14 +29,17 @@ MODULE qs_tddfpt2_restart cp_print_key_should_output,& cp_print_key_unit_nr USE cp_para_types, ONLY: cp_para_env_type + USE dbcsr_api, ONLY: dbcsr_type USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get USE kinds, ONLY: default_path_length,& dp - USE message_passing, ONLY: mp_bcast - USE qs_environment_types, ONLY: get_qs_env,& - qs_environment_type + USE message_passing, ONLY: mp_bcast,& + mp_min,& + mp_sum + USE parallel_gemm_api, ONLY: parallel_gemm + USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos USE string_utilities, ONLY: integer_to_string #include "./base/base_uses.f90" @@ -300,30 +306,39 @@ CONTAINS !> molecular orbitals !> \param logger a logger object !> \param tddfpt_print_section TDDFPT%PRINT input section -!> \param qs_env ... +!> \param matrix_s ... +!> \param S_evects ... +!> \param sub_env ... ! ************************************************************************************************** - SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, qs_env) + SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, & + matrix_s, S_evects, sub_env) TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals TYPE(tddfpt_ground_state_mos), DIMENSION(:), & INTENT(in) :: gs_mos - TYPE(cp_logger_type), POINTER :: logger - TYPE(section_vals_type), POINTER :: tddfpt_print_section - TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env + TYPE(cp_logger_type), INTENT(in), POINTER :: logger + TYPE(section_vals_type), INTENT(in), POINTER :: tddfpt_print_section + TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s + TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects + TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_newtonx_output' INTEGER :: handle, iocc, ispin, istate, ivirt, nao, & nspins, nstates, ounit INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt - LOGICAL :: print_virtuals - TYPE(cp_para_env_type), POINTER :: para_env + LOGICAL :: print_phases, print_virtuals, & + scale_with_phases + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase_evects + TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: evects_mo + TYPE(cp_fm_struct_type), POINTER :: fmstruct IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "NAMD_PRINT"), cp_p_file)) THEN CALL timeset(routineN, handle) - NULLIFY (para_env) - CALL get_qs_env(qs_env, para_env=para_env) + CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals) + CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases) + CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%SCALE_WITH_PHASES", l_val=scale_with_phases) nspins = SIZE(evects, 1) nstates = SIZE(evects, 2) @@ -336,36 +351,101 @@ CONTAINS CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) + IF (sub_env%is_split) THEN + CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when states"// & + " are distributed to different CPU pools.") + END IF + ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "NAMD_PRINT", & extension=".inp", file_form="FORMATTED", file_action="WRITE", file_status="REPLACE") + IF (debug_this_module) CALL tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s) ! print eigenvectors + IF (print_virtuals) THEN + ALLOCATE (evects_mo(nspins, nstates)) + DO istate = 1, nstates + DO ispin = 1, nspins + + ! transform eigenvectors + NULLIFY (fmstruct) + nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) + nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt) + CALL cp_fm_struct_create(fmstruct, para_env=sub_env%para_env, & + context=sub_env%blacs_env, & + nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin)) + NULLIFY (evects_mo(ispin, istate)%matrix) + ALLOCATE (evects_mo(ispin, istate)%matrix) + CALL cp_fm_create(evects_mo(ispin, istate)%matrix, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, istate)%matrix, S_evects(ispin, istate)%matrix, & + ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp) + END DO + END DO + DO istate = 1, nstates + DO ispin = 1, nspins + CALL parallel_gemm("T", "N", & + nmo_virt(ispin), & + nmo_occ(ispin), & + nao, & + 1.0_dp, & + gs_mos(ispin)%mos_virt, & + S_evects(ispin, istate)%matrix, & !this also needs to be orthogonalized + 0.0_dp, & + evects_mo(ispin, istate)%matrix) + END DO + END DO + END IF + DO istate = 1, nstates DO ispin = 1, nspins - ! TDDFPT wave function is actually stored as a linear combination of virtual MOs - ! that replaces the corresponding deoccupied MO. Unfortunately, the phase - ! of the occupied MOs varies depending on the eigensolver used as well as - ! how eigenvectors are distributed across computational cores. The phase is important - ! because TDDFPT wave functions are used to compute a response electron density - ! \rho^{-} = 1/2 * [C_{0} * evect^T + evect * C_{0}^{-}], where C_{0} is the expansion - ! coefficients of the reference ground-state wave function. To make the restart file - ! transferable, TDDFPT wave functions are stored in assumption that all ground state - ! MOs have a positive phase. - CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) - IF (ounit > 0) THEN - WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE" - CALL cp_fm_write_info(evects(ispin, istate)%matrix, ounit) + IF (.NOT. print_virtuals) THEN + CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) + IF (ounit > 0) THEN + WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE" + CALL cp_fm_write_info(evects(ispin, istate)%matrix, ounit) + END IF + CALL cp_fm_write_formatted(evects(ispin, istate)%matrix, ounit, "ES EIGENVECTORS") + ELSE + CALL cp_fm_column_scale(evects_mo(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) + IF (ounit > 0) THEN + WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE" + CALL cp_fm_write_info(evects_mo(ispin, istate)%matrix, ounit) + END IF + CALL cp_fm_write_formatted(evects_mo(ispin, istate)%matrix, ounit, "ES EIGENVECTORS") END IF - CALL cp_fm_write_formatted(evects(ispin, istate)%matrix, ounit, "ES EIGENVECTORS") - CALL cp_fm_column_scale(evects(ispin, istate)%matrix, gs_mos(ispin)%phases_occ) + ! compute and print phase of eigenvectors + nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ) + ALLOCATE (phase_evects(nmo_occ(ispin))) + IF (print_virtuals) THEN + CALL compute_phase_eigenvectors(evects_mo(ispin, istate)%matrix, phase_evects, sub_env) + ELSE + CALL compute_phase_eigenvectors(evects(ispin, istate)%matrix, phase_evects, sub_env) + END IF + IF (ounit > 0) THEN + WRITE (ounit, "(/,A,/)") "PHASES ES EIGENVECTORS" + DO iocc = 1, nmo_occ(ispin) + WRITE (ounit, "(F20.14)") phase_evects(iocc) + END DO + END IF + DEALLOCATE (phase_evects) + END DO END DO - ! print molecular orbitals + IF (print_virtuals) THEN + IF (ASSOCIATED(evects_mo)) THEN + DO istate = 1, nstates + DO ispin = 1, nspins + CALL cp_fm_release(matrix=evects_mo(ispin, istate)%matrix) + DEALLOCATE (evects_mo(ispin, istate)%matrix) + END DO + END DO + END IF + DEALLOCATE (evects_mo) + END IF - CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals) DO ispin = 1, nspins IF (ounit > 0) THEN WRITE (ounit, "(/,A)") "OCCUPIED MOS SIZE" @@ -404,6 +484,27 @@ CONTAINS END IF END IF + ! print phases of molecular orbitals + + IF (print_phases) THEN + IF (ounit > 0) THEN + WRITE (ounit, "(A)") "PHASES OCCUPIED ORBITALS" + DO ispin = 1, nspins + DO iocc = 1, nmo_occ(ispin) + WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_occ(iocc) + END DO + END DO + IF (print_virtuals) THEN + WRITE (ounit, "(A)") "PHASES VIRTUAL ORBITALS" + DO ispin = 1, nspins + DO ivirt = 1, nmo_virt(ispin) + WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_virt(ivirt) + END DO + END DO + END IF + END IF + END IF + CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "NAMD_PRINT") CALL timestop(handle) @@ -411,5 +512,166 @@ CONTAINS END SUBROUTINE tddfpt_write_newtonx_output ! ************************************************************************************************** +!> \brief ... +!> \param evects ... +!> \param ounit ... +!> \param S_evects ... +!> \param matrix_s ... +! ************************************************************************************************** + SUBROUTINE tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s) + + TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: evects + INTEGER, INTENT(in) :: ounit + TYPE(cp_fm_p_type), DIMENSION(:, :), INTENT(in) :: S_evects + TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s + + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_check_orthonormality' + + INTEGER :: handle, ispin, ivect, jvect, nspins, & + nvects_total + INTEGER, DIMENSION(maxspins) :: nactive + REAL(kind=dp) :: norm + REAL(kind=dp), DIMENSION(maxspins) :: weights + + CALL timeset(routineN, handle) + + nspins = SIZE(evects, 1) + nvects_total = SIZE(evects, 2) + + IF (debug_this_module) THEN + CPASSERT(SIZE(S_evects, 1) == nspins) + CPASSERT(SIZE(S_evects, 2) == nvects_total) + END IF + + DO ispin = 1, nspins + CALL cp_fm_get_info(matrix=evects(ispin, 1)%matrix, ncol_global=nactive(ispin)) + END DO + + DO jvect = 1, nvects_total + ! + DO ivect = 1, jvect - 1 + CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.) + norm = SUM(weights(1:nspins)) + + DO ispin = 1, nspins + CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect)%matrix, -norm, evects(ispin, ivect)%matrix) + END DO + END DO + + ! + DO ispin = 1, nspins + CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect)%matrix, S_evects(ispin, jvect)%matrix, & + ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp) + END DO + + CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.) + + norm = SUM(weights(1:nspins)) + norm = 1.0_dp/SQRT(norm) + + IF ((ounit > 0) .AND. debug_this_module) WRITE (ounit, '(A,F10.8)') "norm", norm + + END DO + + CALL timestop(handle) + + END SUBROUTINE tddfpt_check_orthonormality +! ************************************************************************************************** +!> \brief ... +!> \param evects ... +!> \param phase_evects ... +!> \param sub_env ... +! ************************************************************************************************** + SUBROUTINE compute_phase_eigenvectors(evects, phase_evects, sub_env) + + ! copied from parts of tddgpt_init_ground_state_mos by S. Chulkov + + TYPE(cp_fm_type), INTENT(in) :: evects + REAL(kind=dp), DIMENSION(:), INTENT(out) :: phase_evects + TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env + + CHARACTER(len=*), PARAMETER :: routineN = 'compute_phase_eigenvectors' + REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp) + + INTEGER :: handle, icol_global, icol_local, irow_global, irow_local, ncol_global, & + ncol_local, nrow_global, nrow_local, sign_int + INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, & + sum_sign_array + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + REAL(kind=dp) :: element + REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), & + POINTER :: my_block + + CALL timeset(routineN, handle) + + ! compute and print the phase of excited-state eigenvectors: + CALL cp_fm_get_info(evects, nrow_global=nrow_global, ncol_global=ncol_global, & + nrow_local=nrow_local, ncol_local=ncol_local, local_data=my_block, & + row_indices=row_indices, col_indices=col_indices) ! nrow_global either nao or nocc + + ALLOCATE (minrow_neg_array(ncol_global), minrow_pos_array(ncol_global), sum_sign_array(ncol_global)) + minrow_neg_array(:) = nrow_global + minrow_pos_array(:) = nrow_global + sum_sign_array(:) = 0 + + DO icol_local = 1, ncol_local + icol_global = col_indices(icol_local) + + DO irow_local = 1, nrow_local + irow_global = row_indices(irow_local) + + element = my_block(irow_local, icol_local) + + sign_int = 0 + IF (element >= eps_dp) THEN + sign_int = 1 + ELSE IF (element <= -eps_dp) THEN + sign_int = -1 + END IF + + sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int + + IF (sign_int > 0) THEN + IF (minrow_pos_array(icol_global) > irow_global) & + minrow_pos_array(icol_global) = irow_global + ELSE IF (sign_int < 0) THEN + IF (minrow_neg_array(icol_global) > irow_global) & + minrow_neg_array(icol_global) = irow_global + END IF + + END DO + END DO + + CALL mp_sum(sum_sign_array, sub_env%para_env%group) + CALL mp_min(minrow_neg_array, sub_env%para_env%group) + CALL mp_min(minrow_pos_array, sub_env%para_env%group) + + DO icol_global = 1, ncol_global + + IF (sum_sign_array(icol_global) > 0) THEN + ! most of the expansion coefficients are positive => MO's phase = +1 + phase_evects(icol_global) = 1.0_dp + ELSE IF (sum_sign_array(icol_global) < 0) THEN + ! most of the expansion coefficients are negative => MO's phase = -1 + phase_evects(icol_global) = -1.0_dp + ELSE + ! equal number of positive and negative expansion coefficients + IF (minrow_pos_array(icol_global) <= minrow_neg_array(icol_global)) THEN + ! the first positive expansion coefficient has a lower index then + ! the first negative expansion coefficient; MO's phase = +1 + phase_evects(icol_global) = 1.0_dp + ELSE + ! MO's phase = -1 + phase_evects(icol_global) = -1.0_dp + END IF + END IF + + END DO + + DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array) + + CALL timestop(handle) + + END SUBROUTINE compute_phase_eigenvectors END MODULE qs_tddfpt2_restart diff --git a/src/qs_tddfpt2_types.F b/src/qs_tddfpt2_types.F index 5418ba97f1..53e7a51a69 100644 --- a/src/qs_tddfpt2_types.F +++ b/src/qs_tddfpt2_types.F @@ -106,6 +106,7 @@ MODULE qs_tddfpt2_types !> phase of occupied MOs; +1.0 -- positive, -1.0 -- negative; !> it is mainly needed to make the restart file transferable REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phases_occ + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phases_virt END TYPE tddfpt_ground_state_mos ! ************************************************************************************************** diff --git a/src/qs_tddfpt2_utils.F b/src/qs_tddfpt2_utils.F index 5a2e5a3631..dee6d64fab 100644 --- a/src/qs_tddfpt2_utils.F +++ b/src/qs_tddfpt2_utils.F @@ -182,7 +182,8 @@ CONTAINS nlumo=tddfpt_control%nlumo, & blacs_env=blacs_env, cholesky_method=cholesky_restore, & matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, & - mos_virt=mos_virt_spin, evals_virt=evals_virt_spin) + mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, & + qs_env=qs_env) END DO moc = 0 @@ -236,13 +237,14 @@ CONTAINS !> (in addition to the ADDED_MOS, if present). NULL when no OT is in use. !> \param evals_virt orbital energies of precomputed (OT) virtual molecular orbitals. !> NULL when no OT is in use. +!> \param qs_env ... !> \par History !> * 05.2016 created as tddfpt_lumos() [Sergey Chulkov] !> * 06.2016 renamed, altered prototype [Sergey Chulkov] !> * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov] ! ************************************************************************************************** SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, & - mos_virt, evals_virt) + mos_virt, evals_virt, qs_env) TYPE(tddfpt_ground_state_mos) :: gs_mos TYPE(mo_set_type), INTENT(IN) :: mo_set INTEGER, INTENT(in) :: nlumo @@ -251,16 +253,17 @@ CONTAINS TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s TYPE(cp_fm_type), INTENT(IN), POINTER :: mos_virt REAL(kind=dp), DIMENSION(:), POINTER :: evals_virt + TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_ground_state_mos' REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp) - INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, irow_global, & + INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, & irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, & sum_sign_array INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices - LOGICAL :: do_eigen + LOGICAL :: do_eigen, print_phases REAL(kind=dp) :: element, maxocc REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), & POINTER :: my_block @@ -268,18 +271,28 @@ CONTAINS mo_occ_scf TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct, ao_mo_occ_fm_struct, & ao_mo_virt_fm_struct, wfn_fm_struct - TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm + TYPE(cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, & + work_fm_virt TYPE(cp_fm_type), POINTER :: mo_coeff_extended + TYPE(cp_logger_type), POINTER :: logger TYPE(cp_para_env_type), POINTER :: para_env TYPE(mo_set_type), POINTER :: mos_extended + TYPE(section_vals_type), POINTER :: print_section CALL timeset(routineN, handle) + NULLIFY (logger) + logger => cp_get_default_logger() + iounit = cp_logger_get_default_io_unit(logger) + CALL get_blacs_info(blacs_env, para_env=para_env) CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, & nelectron=nelectrons, occupation_numbers=mo_occ_scf) + print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT") + CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases) + nmo_virt = nao - nmo_occ IF (nlumo >= 0) & nmo_virt = MIN(nmo_virt, nlumo) @@ -291,6 +304,7 @@ CONTAINS do_eigen = .FALSE. ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small IF (ASSOCIATED(evals_virt)) THEN + CPASSERT(ASSOCIATED(mos_virt)) IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE. ELSE IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE. @@ -309,8 +323,7 @@ CONTAINS ALLOCATE (gs_mos%evals_occ(nmo_occ)) ALLOCATE (gs_mos%evals_virt(nmo_virt)) ALLOCATE (gs_mos%phases_occ(nmo_occ)) - - CALL cp_fm_struct_release(ao_mo_virt_fm_struct) + ALLOCATE (gs_mos%phases_virt(nmo_virt)) ! ++ nullify pointers NULLIFY (ao_ao_fm_struct, wfn_fm_struct) @@ -457,6 +470,73 @@ CONTAINS gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt) END IF + IF (print_phases) THEN + ! compute the phase of molecular orbitals; + ! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors + !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env) + CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct) + + CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1) + CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, & + row_indices=row_indices, col_indices=col_indices, local_data=my_block) + + ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt)) + minrow_neg_array(:) = nao + minrow_pos_array(:) = nao + sum_sign_array(:) = 0 + DO icol_local = 1, ncol_local + icol_global = col_indices(icol_local) + + DO irow_local = 1, nrow_local + element = my_block(irow_local, icol_local) + + sign_int = 0 + IF (element >= eps_dp) THEN + sign_int = 1 + ELSE IF (element <= -eps_dp) THEN + sign_int = -1 + END IF + + sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int + + irow_global = row_indices(irow_local) + IF (sign_int > 0) THEN + IF (minrow_pos_array(icol_global) > irow_global) & + minrow_pos_array(icol_global) = irow_global + ELSE IF (sign_int < 0) THEN + IF (minrow_neg_array(icol_global) > irow_global) & + minrow_neg_array(icol_global) = irow_global + END IF + END DO + END DO + + CALL mp_sum(sum_sign_array, para_env%group) + CALL mp_min(minrow_neg_array, para_env%group) + CALL mp_min(minrow_pos_array, para_env%group) + DO icol_local = 1, nmo_virt + IF (sum_sign_array(icol_local) > 0) THEN + ! most of the expansion coefficients are positive => MO's phase = +1 + gs_mos%phases_virt(icol_local) = 1.0_dp + ELSE IF (sum_sign_array(icol_local) < 0) THEN + ! most of the expansion coefficients are negative => MO's phase = -1 + gs_mos%phases_virt(icol_local) = -1.0_dp + ELSE + ! equal number of positive and negative expansion coefficients + IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN + ! the first positive expansion coefficient has a lower index then + ! the first negative expansion coefficient; MO's phase = +1 + gs_mos%phases_virt(icol_local) = 1.0_dp + ELSE + ! MO's phase = -1 + gs_mos%phases_virt(icol_local) = -1.0_dp + END IF + END IF + END DO + DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array) + CALL cp_fm_release(work_fm_virt) + END IF !print_phases + CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases + CALL cp_fm_release(work_fm) IF (do_eigen) THEN @@ -492,6 +572,9 @@ CONTAINS IF (ALLOCATED(gs_mos%evals_occ)) & DEALLOCATE (gs_mos%evals_occ) + IF (ALLOCATED(gs_mos%phases_virt)) & + DEALLOCATE (gs_mos%phases_virt) + IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN CALL cp_fm_release(gs_mos%evals_occ_matrix) DEALLOCATE (gs_mos%evals_occ_matrix) diff --git a/tests/QS/regtest-gpw-9/h2o_NEWTONX.inp b/tests/QS/regtest-gpw-9/h2o_NEWTONX.inp index 8a2c5a8c62..bf9f8f8f97 100644 --- a/tests/QS/regtest-gpw-9/h2o_NEWTONX.inp +++ b/tests/QS/regtest-gpw-9/h2o_NEWTONX.inp @@ -43,9 +43,8 @@ &END FORCE_EVAL &VIBRATIONAL_ANALYSIS &PRINT - &CARTESIAN_EIGS - &END CARTESIAN_EIGS + &NAMD_PRINT + &END NAMD_PRINT &END PRINT - NAMD_PRINT DX 0.001 &END VIBRATIONAL_ANALYSIS