projection on reference in RTP (cp2k#3298)

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

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

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.)

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)

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. "// &

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

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

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

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

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

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

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