Extension to the XAS code based on the Slater transition potential me…

…thod.

The extension allows multiple initial core excited states, of any type
(s, p, d, ..). Code by Axel Erbung (Dep. Physics Stockholm)



git-svn-id: svn://svn.code.sf.net/p/cp2k/code/trunk@18514 bac0169d-1eea-4375-b5f6-f89d2ced9994

cp2k/src/input_cp2k_dft.F

@@ -89,7 +89,8 @@ MODULE input_cp2k_dft
         weight_type_unit, wfi_aspc_nr, wfi_frozen_method_nr, wfi_linear_p_method_nr, &
         wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, wfi_ps_method_nr, &
         wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, &
-        wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, &
+        wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, xas_2s_type, xas_3d_type, &
+        xas_3p_type, xas_3s_type, xas_4d_type, xas_4f_type, xas_4p_type, xas_4s_type, xas_dip_len, &
         xas_dip_vel, xas_dscf, xas_none, xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, &
         xas_tp_xhh, xes_tp_val
    USE input_cp2k_almo,                 ONLY: create_almo_scf_section
@@ -7585,9 +7586,11 @@ SUBROUTINE create_xas_section(section)
                           description="Type of the orbitals that are excited for the xas spectra calculation", &
                           usage="STATE_TYPE 1S", &
                           default_i_val=xas_1s_type, &
-                          enum_c_vals=s2a("1S", "2S", "2P"), &
-                          enum_desc=s2a("1s orbitals", "2s orbitals", "2p orbitals"), &
-                          enum_i_vals=(/xas_1s_type, xas_2s_type, xas_2p_type/))
+                          enum_c_vals=s2a("1S", "2S", "2P", "3S", "3P", "3D", "4S", "4P", "4D", "4F"), &
+                          enum_desc=s2a("1s orbitals", "2s orbitals", "2p orbitals", "3s orbitals", "3p orbitals", &
+                                        "3d orbitals", "4s orbitals", "4p orbitals", "4d orbitals", "4f orbitals"), &
+                          enum_i_vals=(/xas_1s_type, xas_2s_type, xas_2p_type, xas_3s_type, xas_3p_type, xas_3d_type, &
+                                        xas_4s_type, xas_4p_type, xas_4d_type, xas_4f_type/))
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
@@ -7615,6 +7618,23 @@ SUBROUTINE create_xas_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, name="OVERLAP_THRESHOLD", &
+                          description="Threshold for including more than one initial core excited state "// &
+                          "per atom. The threshold is taken relative to the maximum overlap.", &
+                          usage="OVERLAP_THRESHOLD 8.e-1", default_r_val=1.0_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="ORBITAL_LIST", &
+                          variants=(/"ORBITAL_LIST"/), &
+                          description="Indices of the localized orbitals to be excited. "// &
+                          "This keyword can be repeated several times"// &
+                          "(useful if you have to specify many indexes).", &
+                          usage="ORBITAL_LIST {integer}  {integer} ..  {integer} ", &
+                          n_var=-1, type_of_var=integer_t, repeats=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, name="ADDED_MOS", &
                           description="Number of additional MOS added spin up only", &
                           usage="ADDED_MOS {integer}", default_i_val=-1)

cp2k/src/xas_control.F

@@ -48,11 +48,12 @@ MODULE xas_control
       INTEGER                             :: ngauss
       INTEGER                             :: stride
       INTEGER, DIMENSION(:), POINTER      :: exc_atoms
+      INTEGER, DIMENSION(:), POINTER      :: orbital_list
       LOGICAL                             :: cubes, do_centers
       LOGICAL                             :: xas_restart
       INTEGER, DIMENSION(:), POINTER      :: list_cubes
 !
-      REAL(dp)                            :: eps_added
+      REAL(dp)                            :: eps_added, overlap_threshold
       REAL(dp)                            :: xes_core_occupation
       REAL(dp)                            :: xes_homo_occupation
       REAL(dp)                            :: nel_tot, xas_core_occupation
@@ -87,7 +88,7 @@ SUBROUTINE read_xas_control(xas_control, xas_section)
       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_xas_control', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: i, ir, n_rep, nex_at
+      INTEGER                                            :: i, ir, n_rep, nex_at, nex_st
       INTEGER, DIMENSION(:), POINTER                     :: list
       LOGICAL                                            :: hempty, was_present
 
@@ -172,6 +173,31 @@ SUBROUTINE read_xas_control(xas_control, xas_section)
       CALL section_vals_val_get(xas_section, "NGAUSS", &
                                 i_val=xas_control%ngauss)
 
+      CALL section_vals_val_get(xas_section, "OVERLAP_THRESHOLD", &
+                                r_val=xas_control%overlap_threshold)
+
+      CALL section_vals_val_get(xas_section, "ORBITAL_LIST", &
+                                n_rep_val=n_rep)
+      IF (n_rep > 0) THEN
+         nex_st = 0
+         DO ir = 1, n_rep
+            NULLIFY (list)
+            CALL section_vals_val_get(xas_section, "ORBITAL_LIST", &
+                                      i_rep_val=ir, i_vals=list)
+
+            IF (ASSOCIATED(list)) THEN
+               CALL reallocate(xas_control%orbital_list, 1, nex_st+SIZE(list))
+               DO i = 1, SIZE(list)
+                  xas_control%orbital_list(i+nex_st) = list(i)
+               END DO
+               nex_st = nex_st+SIZE(list)
+            END IF
+         END DO ! ir
+      ELSE
+         ALLOCATE (xas_control%orbital_list(1))
+         xas_control%orbital_list(1) = -1
+      END IF
+
    END SUBROUTINE read_xas_control
 
 ! **************************************************************************************************
@@ -267,6 +293,7 @@ SUBROUTINE xas_control_create(xas_control)
       xas_control%xes_core_occupation = 1.0_dp
       xas_control%xes_homo_occupation = 1.0_dp
       NULLIFY (xas_control%exc_atoms)
+      NULLIFY (xas_control%orbital_list)
       xas_control%cubes = .FALSE.
       xas_control%do_centers = .FALSE.
       NULLIFY (xas_control%list_cubes)
@@ -291,6 +318,9 @@ SUBROUTINE xas_control_release(xas_control)
             IF (ASSOCIATED(xas_control%exc_atoms)) THEN
                DEALLOCATE (xas_control%exc_atoms)
             END IF
+            IF (ASSOCIATED(xas_control%orbital_list)) THEN
+               DEALLOCATE (xas_control%orbital_list)
+            END IF
             IF (ASSOCIATED(xas_control%list_cubes)) THEN
                DEALLOCATE (xas_control%list_cubes)
             END IF
@@ -299,7 +329,6 @@ SUBROUTINE xas_control_release(xas_control)
       END IF
 
    END SUBROUTINE xas_control_release
-
 ! **************************************************************************************************
 !> \brief ...
 !> \param xas_control ...

cp2k/src/xas_env_types.F

@@ -53,10 +53,12 @@ MODULE xas_env_types
 !> \param nao number of atomic orbitals in the basis
 !> \param exc_state state that is now excited (this change atom by atom)
 !> \param nvirtual number of empy states to take into account for the spectrum
-!> \param state_of_atom for each atom the state that have to be excited (global index)
-!>                      dimension is the number of atoms to be excited
+!> \param state_of_atom for each atom the states that have to be excited (global index)
+!>        dimension is the number of atoms to be excited by the largest number of included states
 !> \param atom_of_state atom to which each state is assigned,
-!>                     dimension is the number of states occupied that might be excited
+!>        dimension is the number of states occupied that might be excited
+!> \param nexc_states number of states to be excited per atom
+!>        dimension is the number of atoms to be excited
 !> \param type_of_state character of the state (1s,2s,2p...)
 !> \param spectrum for each excitation the energy and the oscillator strength
 !> \param centers_wfn for each wfn the center of charge (optimized by localization)
@@ -76,12 +78,15 @@ MODULE xas_env_types
       INTEGER :: spin_channel
       INTEGER :: nvirtual, nvirtual2
       INTEGER :: unoccupied_max_iter
-      INTEGER, DIMENSION(:), POINTER :: state_of_atom
+
       INTEGER, DIMENSION(:), POINTER :: atom_of_state
       INTEGER, DIMENSION(:), POINTER :: type_of_state
       INTEGER, DIMENSION(:), POINTER :: mykind_of_atom
       INTEGER, DIMENSION(:), POINTER :: mykind_of_kind
       INTEGER, DIMENSION(:), POINTER :: exc_atoms
+      INTEGER, DIMENSION(:), POINTER :: nexc_states
+      INTEGER, DIMENSION(:, :), POINTER :: state_of_atom
+
       REAL(dp) :: ip_energy, occ_estate, unoccupied_eps, xas_nelectron, homo_occ
       REAL(dp), DIMENSION(:), POINTER :: all_evals
       REAL(dp), DIMENSION(:), POINTER :: unoccupied_evals
@@ -116,6 +121,7 @@ MODULE xas_env_types
 !> \param centers_wfn ...
 !> \param atom_of_state ...
 !> \param exc_atoms ...
+!> \param nexc_states ...
 !> \param type_of_state ...
 !> \param mykind_of_atom ...
 !> \param mykind_of_kind ...
@@ -145,8 +151,8 @@ MODULE xas_env_types
 !> \param scf_control ...
 ! **************************************************************************************************
    SUBROUTINE get_xas_env(xas_env, iter_count, exc_state, nao, nvirtual, nvirtual2, &
-                          centers_wfn, atom_of_state, exc_atoms, type_of_state, mykind_of_atom, mykind_of_kind, &
-                          state_of_atom, spectrum, groundstate_coeff, ostrength_sm, &
+                          centers_wfn, atom_of_state, exc_atoms, nexc_states, type_of_state, mykind_of_atom, &
+                          mykind_of_kind, state_of_atom, spectrum, groundstate_coeff, ostrength_sm, &
                           dip_fm_set, excvec_coeff, excvec_overlap, &
                           unoccupied_orbs, unoccupied_evals, unoccupied_max_iter, unoccupied_eps, &
                           all_vectors, all_evals, my_gto_basis, qs_loc_env, &
@@ -157,9 +163,10 @@ SUBROUTINE get_xas_env(xas_env, iter_count, exc_state, nao, nvirtual, nvirtual2,
       INTEGER, INTENT(OUT), OPTIONAL                     :: iter_count, exc_state, nao, nvirtual, &
                                                             nvirtual2
       REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER       :: centers_wfn
-      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: atom_of_state, exc_atoms, type_of_state, &
-                                                            mykind_of_atom, mykind_of_kind, &
-                                                            state_of_atom
+      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: atom_of_state, exc_atoms, nexc_states, &
+                                                            type_of_state, mykind_of_atom, &
+                                                            mykind_of_kind
+      INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: state_of_atom
       REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER       :: spectrum
       TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, &
          POINTER                                         :: groundstate_coeff
@@ -198,6 +205,7 @@ SUBROUTINE get_xas_env(xas_env, iter_count, exc_state, nao, nvirtual, nvirtual2,
       IF (PRESENT(occ_estate)) occ_estate = xas_env%occ_estate
       IF (PRESENT(xas_estate)) xas_estate = xas_env%xas_estate
       IF (PRESENT(nexc_search)) nexc_search = xas_env%nexc_search
+      IF (PRESENT(nexc_states)) nexc_states => xas_env%nexc_states
       IF (PRESENT(spin_channel)) spin_channel = xas_env%spin_channel
       IF (PRESENT(nexc_atoms)) nexc_atoms = xas_env%nexc_atoms
       IF (PRESENT(unoccupied_eps)) unoccupied_eps = xas_env%unoccupied_eps
@@ -311,7 +319,9 @@ SUBROUTINE xas_env_create(xas_env)
 
       NULLIFY (xas_env%ao_mo_fm_pools)
       NULLIFY (xas_env%my_gto_basis)
-      NULLIFY (xas_env%state_of_atom, xas_env%atom_of_state)
+      NULLIFY (xas_env%atom_of_state)
+      NULLIFY (xas_env%nexc_states)
+      NULLIFY (xas_env%state_of_atom)
       NULLIFY (xas_env%exc_atoms)
       NULLIFY (xas_env%excvec_coeff, xas_env%excvec_overlap)
       NULLIFY (xas_env%type_of_state, xas_env%mykind_of_atom)
@@ -345,6 +355,7 @@ SUBROUTINE xas_env_release(xas_env)
          IF (xas_env%ref_count == 0) THEN
 
             DEALLOCATE (xas_env%state_of_atom, xas_env%atom_of_state)
+            DEALLOCATE (xas_env%nexc_states)
             DEALLOCATE (xas_env%type_of_state)
             DEALLOCATE (xas_env%mykind_of_atom)
             DEALLOCATE (xas_env%mykind_of_kind)