diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 7d492b1853..6a23fa3fc3 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -44,8 +44,10 @@ MODULE cp_control_types TYPE proj_mo_type INTEGER, DIMENSION(:), ALLOCATABLE :: ref_mo_index INTEGER :: ref_mo_spin = 1 + INTEGER :: ref_nlumo = 0 LOGICAL :: sum_on_all_ref = .FALSE. INTEGER, DIMENSION(:), ALLOCATABLE :: td_mo_index + REAL(dp), DIMENSION(:), ALLOCATABLE :: td_mo_occ INTEGER :: td_mo_spin = 1 LOGICAL :: sum_on_all_td = .FALSE. CHARACTER(LEN=default_path_length) :: ref_mo_file_name = "" @@ -980,6 +982,8 @@ SUBROUTINE proj_mo_list_release(proj_mo_list) END IF IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_index)) & DEALLOCATE (proj_mo_list(i)%proj_mo%td_mo_index) + IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_occ)) & + DEALLOCATE (proj_mo_list(i)%proj_mo%td_mo_occ) DEALLOCATE (proj_mo_list(i)%proj_mo) END IF END DO diff --git a/src/emd/rt_projection_mo_utils.F b/src/emd/rt_projection_mo_utils.F index 6acc5914d8..92bf89a3ed 100644 --- a/src/emd/rt_projection_mo_utils.F +++ b/src/emd/rt_projection_mo_utils.F @@ -50,7 +50,6 @@ MODULE rt_projection_mo_utils USE qs_kind_types, ONLY: qs_kind_type USE qs_mo_io, ONLY: read_mos_restart_low USE qs_mo_types, ONLY: deallocate_mo_set,& - duplicate_mo_set,& mo_set_type USE rt_propagation_types, ONLY: get_rtp,& rt_prop_type @@ -105,6 +104,9 @@ SUBROUTINE init_mo_projection(qs_env, rtp_control) CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, & c_val=proj_mo%ref_mo_file_name) + CALL section_vals_val_get(proj_mo_section, "REF_ADD_LUMO", i_rep_section=i_rep, & + i_val=proj_mo%ref_nlumo) + ! Relevent only in EMD IF (.NOT. rtp_control%fixed_ions) & CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, & @@ -165,20 +167,26 @@ SUBROUTINE init_mo_projection(qs_env, rtp_control) CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, & i_vals=tmp_ints) - ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:)) nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global + + ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:)) IF (proj_mo%td_mo_index(1) == -1) THEN DEALLOCATE (proj_mo%td_mo_index) ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max)) + ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max)) DO j_td = 1, nbr_mo_td_max proj_mo%td_mo_index(j_td) = j_td + proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td)) END DO ELSE + ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index))) + proj_mo%td_mo_occ(:) = 0.0_dp DO j_td = 1, SIZE(proj_mo%td_mo_index) IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) & CALL cp_abort(__LOCATION__, & "The MO number available in the Time Dependent run "// & "is smaller than the MO number you have required in TD_MO_INDEX.") + proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td)) END DO END IF @@ -210,11 +218,12 @@ SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref) TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_qs, mo_ref_temp + TYPE(mo_set_type), POINTER :: mo_set TYPE(mp_para_env_type), POINTER :: para_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - NULLIFY (mo_qs, mo_ref_temp, qs_kind_set, particle_set, para_env, dft_control, & + NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, & mo_ref_fmstruct, matrix_s) my_xasref = .FALSE. @@ -230,20 +239,44 @@ SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref) natom = SIZE(particle_set, 1) + nspins = SIZE(mo_qs) ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved - IF (my_xasref) THEN + IF (my_xasref .AND. nspins < 2) THEN nspins = 2 - ALLOCATE (mo_ref_temp(nspins)) - DO ispin = 1, nspins - CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1)) - END DO - ELSE - nspins = SIZE(mo_qs) - ALLOCATE (mo_ref_temp(nspins)) - DO ispin = 1, nspins - CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin)) - END DO END IF + ALLOCATE (mo_ref_temp(nspins)) + + DO ispin = 1, nspins + IF (my_xasref) THEN + mo_set => mo_qs(1) + ELSE + mo_set => mo_qs(ispin) + END IF + mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo + NULLIFY (mo_ref_fmstruct) + CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, & + ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context) + NULLIFY (mo_ref_temp(ispin)%mo_coeff) + ALLOCATE (mo_ref_temp(ispin)%mo_coeff) + CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct) + CALL cp_fm_struct_release(mo_ref_fmstruct) + + mo_ref_temp(ispin)%nao = mo_set%nao + mo_ref_temp(ispin)%homo = mo_set%homo + mo_ref_temp(ispin)%nelectron = mo_set%nelectron + ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo)) + ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo)) + NULLIFY (mo_set) + END DO + +! DO ispin = 1, nspins +! CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1)) +! END DO +! ELSE +! DO ispin = 1, nspins +! CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin)) +! END DO +! END IF IF (para_env%is_source()) THEN INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file) @@ -490,7 +523,6 @@ SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref) ncol=1, & source_start=proj_mo%td_mo_index(j_td), & target_start=1) - IF (is_emd) THEN ! The reference MO have to be propagated in the new basis, so the projection CALL cp_fm_trace(mo_coeff_temp, S_mo_ref, real_proj) @@ -514,7 +546,7 @@ SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref) ! Store the result phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians - popu(j_td) = real_proj**2 + imag_proj**2 + popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2) END DO CALL cp_fm_struct_release(mo_ref_fmstruct) @@ -572,9 +604,9 @@ SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, IF (print_unit /= output_unit) THEN INQUIRE (UNIT=print_unit, NAME=filename) - WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") & - "PROJECTION MO", "The projection of the TD MOs is done in the file:", & - TRIM(filename) +! WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") & +! "PROJECTION MO", "The projection of the TD MOs is done in the file:", & +! TRIM(filename) WRITE (UNIT=print_unit, FMT="(/,(T2,A,T40,I6))") & "Real time propagation step:", qs_env%sim_step ELSE diff --git a/src/emd/rt_propagation_methods.F b/src/emd/rt_propagation_methods.F index 88d40653db..649e132429 100644 --- a/src/emd/rt_propagation_methods.F +++ b/src/emd/rt_propagation_methods.F @@ -13,6 +13,7 @@ MODULE rt_propagation_methods USE bibliography, ONLY: Kolafa2004,& Kuhne2007,& cite_reference + USE cell_types, ONLY: cell_type USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,& cp_cfm_triangular_multiply USE cp_cfm_types, ONLY: cp_cfm_create,& @@ -109,6 +110,7 @@ SUBROUTINE propagation_step(qs_env, rtp, rtp_control) CHARACTER(len=*), PARAMETER :: routineN = 'propagation_step' INTEGER :: aspc_order, handle, i, im, re, unit_nr + TYPE(cell_type), POINTER :: cell TYPE(cp_fm_type), DIMENSION(:), POINTER :: delta_mos, mos_new TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P, H_last_iter, ks_mix, ks_mix_im, & @@ -125,10 +127,10 @@ SUBROUTINE propagation_step(qs_env, rtp, rtp_control) unit_nr = -1 END IF - NULLIFY (delta_P, rho_new, delta_mos, mos_new) + NULLIFY (cell, delta_P, rho_new, delta_mos, mos_new) NULLIFY (ks_mix, ks_mix_im) ! get everything needed and set some values - CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control) + CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, dft_control=dft_control) IF (rtp%iter == 1) THEN CALL qs_energies_init(qs_env, .FALSE.) @@ -137,10 +139,12 @@ SUBROUTINE propagation_step(qs_env, rtp, rtp_control) ! add additional terms to matrix_h and matrix_h_im in the case of applied electric field, ! either in the lengh or velocity gauge. - ! should be called imediately after qs_energies_init and before qs_ks_update_qs_env - IF (dft_control%apply_efield_field) & + ! should be called after qs_energies_init and before qs_ks_update_qs_env + IF (dft_control%apply_efield_field) THEN + IF (ANY(cell%perd(1:3) /= 0)) & + CPABORT("Length gauge (efield) and periodicity are not compatible") CALL efield_potential_lengh_gauge(qs_env) - IF (rtp_control%velocity_gauge) THEN + ELSE IF (rtp_control%velocity_gauge) THEN IF (dft_control%apply_vector_potential) & CALL update_vector_potential(qs_env, dft_control) CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.) diff --git a/src/emd/rt_propagation_utils.F b/src/emd/rt_propagation_utils.F index 3fe6d5f792..19fac8484c 100644 --- a/src/emd/rt_propagation_utils.F +++ b/src/emd/rt_propagation_utils.F @@ -263,7 +263,7 @@ SUBROUTINE get_restart_wfn(qs_env) REAL(KIND=dp) :: alpha, cs_pos TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_fm_type) :: mos_occ - TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old + TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_distribution_type) :: dist TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv, rho_new, rho_old @@ -379,7 +379,7 @@ SUBROUTINE get_restart_wfn(qs_env) END DO CALL calc_update_rho_sparse(qs_env) ELSE - CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new) + CALL get_rtp(rtp=rtp, mos_old=mos_old) CALL read_rt_mos_from_restart(mo_array, mos_old, atomic_kind_set, qs_kind_set, particle_set, para_env, & id_nr, dft_control%multiplicity, dft_section) CALL set_uniform_occupation_mo_array(mo_array, nspin) diff --git a/src/input_cp2k_projection_rtp.F b/src/input_cp2k_projection_rtp.F index 653d0d23de..25c9ff406b 100644 --- a/src/input_cp2k_projection_rtp.F +++ b/src/input_cp2k_projection_rtp.F @@ -117,6 +117,16 @@ SUBROUTINE create_projection_rtp_section(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="REF_ADD_LUMO", & + description="If the reference MOs include more empty states that are not propagated, "// & + "using this keyword it is possible to read them as well and thus compute the corresponding projection. ", & + usage="REF_ADD_LUMO 10", & + default_i_val=0, & + n_var=1, type_of_var=integer_t, repeats=.FALSE.) + + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="SUM_ON_ALL_REF", & description="Set to .TRUE. in order to sum all the projections done "// & "over the required MO_REF_INDEX for each TD MOs. "// & diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp index 409bd4f3f6..8b4091f9e0 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp @@ -65,6 +65,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD H 6.000000 4.5000000 6.000000 diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp b/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp index 634587bbd2..d127105f40 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp @@ -69,6 +69,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD H 6.000000 4.5000000 6.000000 diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield.inp b/tests/QS/regtest-rtp-2/H2-emd-efield.inp index 20422462f0..f449c31e3b 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield.inp @@ -67,6 +67,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD H 6.000000 4.5000000 6.000000 diff --git a/tests/QS/regtest-rtp-2/H2-emd_ETRS_ARNOLDI.inp b/tests/QS/regtest-rtp-2/H2-emd_ETRS_ARNOLDI.inp index 76245ea1db..895859d6fe 100644 --- a/tests/QS/regtest-rtp-2/H2-emd_ETRS_ARNOLDI.inp +++ b/tests/QS/regtest-rtp-2/H2-emd_ETRS_ARNOLDI.inp @@ -69,6 +69,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD O 0.000000 0.000000 -0.065587 H2O diff --git a/tests/QS/regtest-rtp-2/H2-rtp-efield.inp b/tests/QS/regtest-rtp-2/H2-rtp-efield.inp index 2597815999..cfdf4614c2 100644 --- a/tests/QS/regtest-rtp-2/H2-rtp-efield.inp +++ b/tests/QS/regtest-rtp-2/H2-rtp-efield.inp @@ -108,6 +108,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD H 6.000000 4.5000000 6.000000 diff --git a/tests/QS/regtest-rtp-2/H2-rtp_ETRS_ARNOLDI.inp b/tests/QS/regtest-rtp-2/H2-rtp_ETRS_ARNOLDI.inp index f90611f260..4523e4453f 100644 --- a/tests/QS/regtest-rtp-2/H2-rtp_ETRS_ARNOLDI.inp +++ b/tests/QS/regtest-rtp-2/H2-rtp_ETRS_ARNOLDI.inp @@ -59,6 +59,7 @@ &SUBSYS &CELL ABC 6.0 6.0 6.0 + PERIODIC NONE &END CELL &COORD H 6.000000 4.5000000 6.000000 diff --git a/tests/QS/regtest-rtp-2/TEST_FILES b/tests/QS/regtest-rtp-2/TEST_FILES index 36d0419fd3..55301cee9f 100644 --- a/tests/QS/regtest-rtp-2/TEST_FILES +++ b/tests/QS/regtest-rtp-2/TEST_FILES @@ -7,13 +7,13 @@ H2-emd_restart.inp 2 1.0E-14 - H2-emd_restart-1.inp 2 1.0E-14 -0.902241027707E+00 H2-rtp_restart.inp 1 3e-13 -0.90223968349591 H2-rtp_restart-1.inp 1 3e-13 -0.90223968361437 -H2-rtp-efield.inp 1 2e-10 -0.82613324527108 +H2-rtp-efield.inp 1 2e-10 -0.8312583089 H2-rtp-efield-vg.inp 1 2e-10 -0.83550227784861 H2-rtp-efield-vg-restart.inp 1 2e-10 -0.79589429986106 -H2-emd-efield.inp 2 4e-11 -0.894818188949 -H2-emd-efield-ramp.inp 2 5e-12 -0.885822019185 -H2-emd-efield-custom.inp 2 1e-14 -0.902238043295 -H2-rtp_ETRS_ARNOLDI.inp 1 3e-13 -0.90223968349550 -H2-emd_ETRS_ARNOLDI.inp 1 1e-10 -17.08556435992424 +H2-emd-efield.inp 2 4e-11 -0.8760125223 +H2-emd-efield-ramp.inp 2 5e-12 -0.8704554332 +H2-emd-efield-custom.inp 2 1e-10 -0.8805787785 +H2-rtp_ETRS_ARNOLDI.inp 1 3e-10 -0.8805676355 +H2-emd_ETRS_ARNOLDI.inp 1 1e-10 -17.08579286 H2-emd_VG_ETRS_ARNOLDI.inp 1 1e-10 -17.08551342644315 #EOF