XAS_TDP| Full TDDFT reformulated as an Hermitian eigenvalue problem

src/cp_control_types.F

@@ -478,7 +478,8 @@ MODULE cp_control_types
       INTEGER                              :: auto_basis_ri_aux = 1, &
                                               auto_basis_aux_fit = 1, &
                                               auto_basis_lri_aux = 1, &
-                                              auto_basis_ri_hxc = 1
+                                              auto_basis_ri_hxc = 1, &
+                                              auto_basis_ri_xas = 1
       REAL(KIND=dp)                        :: relax_multiplicity, &
                                               sic_scaling_a, &
                                               sic_scaling_b

src/cp_control_utils.F

@@ -206,6 +206,8 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
                dft_control%auto_basis_lri_aux = isize
             CASE ("RI_HXC")
                dft_control%auto_basis_ri_hxc = isize
+            CASE ("RI_XAS")
+               dft_control%auto_basis_ri_xas = isize
             CASE DEFAULT
                CPWARN("Unknown basis type in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
             END SELECT

src/input_constants.F

@@ -626,6 +626,12 @@ MODULE input_constants
    INTEGER, PARAMETER, PUBLIC               :: xas_scf_default = 0, &
                                                xas_scf_general = 1
 
+   ! Time-dependent XAS
+   INTEGER, PARAMETER, PUBLIC               :: xas_tdp_by_index = 1, &
+                                               xas_tdp_by_kind = 2, &
+                                               xas_tdp_hfx = 1, &
+                                               xas_tdp_pade = 2 
+
    ! Form of dipole operator for TDDFPT oscillator strength calculation
    INTEGER, PARAMETER, PUBLIC               :: tddfpt_dipole_berry = 1, &
                                                tddfpt_dipole_length = 2, &

src/input_cp2k_dft.F

@@ -96,7 +96,7 @@ MODULE input_cp2k_dft
         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, xas_not_excited, xas_tdp_by_kind, &
-        xas_tdp_by_index, xas_tdp_exop_coulomb
+        xas_tdp_by_index, xas_tdp_hfx, xas_tdp_pade
    USE input_cp2k_almo,                 ONLY: create_almo_scf_section
    USE input_cp2k_distribution,         ONLY: create_distribution_section
    USE input_cp2k_field,                ONLY: create_efield_section,&
@@ -7854,7 +7854,7 @@ SUBROUTINE create_xas_section(section)
       CALL section_add_subsection(subsection, print_key)
       CALL section_release(print_key)
 
-      CALL section_add_subsection(section, subsection)
+      CALL section_add_subsection(section, subsection) 
       CALL section_release(subsection)
 
    END SUBROUTINE create_xas_section
@@ -7880,7 +7880,7 @@ SUBROUTINE create_xas_tdp_section(section)
                           description="XAS simulations using TDDFPT. Excitation from specified "// &
                           "core orbitals are considered one at a time. In case of high symmetry "// &
                           "structures, donor core orbitals should be localized.", &
-                          n_keywords=11, n_subsections=3, repeats=.FALSE.)
+                          n_keywords=11, n_subsections=4, repeats=.FALSE.)
 
       NULLIFY (keyword, subsection, print_key)
 
@@ -8006,6 +8006,39 @@ SUBROUTINE create_xas_tdp_section(section)
                           default_l_val=.FALSE.)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
+
+      !*** XC_KERNEL subsection ***!
+      CALL section_create(subsection, "XC_KERNEL", &
+                          description="This subsection sets parameter for the xc kernel "//&
+                                      "evaluation such as the functional to use and the atomic "//&
+                                      "grids specification", &
+                                      n_keywords=2, n_subsections=0, repeats=.FALSE.)
+
+      CALL keyword_create(keyword, name="FUNCTIONAL", &
+                          description="This keyword is used to chose the functional with which "//&
+                                      "the XC kernel should be evaluated. Currently only exact "//&
+                                      "exchange and the PADE(LDA) are available.", &
+                          usage="FUNCTIONAL string", &
+                          default_i_val=xas_tdp_hfx, &
+                          enum_c_vals=s2a("HFX", "PADE"), &
+                          enum_i_vals=(/xas_tdp_hfx, xas_tdp_pade/))
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="GRID", &
+                          description="Specification of the atomic grids for the excited kinds."//&
+                                      "This keyword can/should be repeated for each excited kind. "//&
+                                      "The default grid dimensions are those set for the GAPW "// &
+                                      "ground state calculation.", &
+                          usage="GRID kind na nr", &
+                          n_var=3, type_of_var=char_t, repeats=.TRUE.)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
+      !*** END of XC_KERNEL subsection ***!
 
       CALL create_localize_section(subsection)
       CALL section_add_subsection(section, subsection)

src/particle_methods.F

@@ -87,6 +87,7 @@ MODULE particle_methods
 !> \par History
 !>      - particle type cleaned (13.10.2003,MK)
 !>      - refactoring and add basis set option (17.08.2010,jhu)
+!>      - if basis option present but not associated, return nsgf=1 for this basis (06.08.2018 AB)
 !> \author  MK
 !> \version 1.0
 ! **************************************************************************************************
@@ -125,7 +126,11 @@ SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf
          DO iparticle = 1, nparticle
             CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
             IF (PRESENT(basis)) THEN
-               CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
+               IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
+                  CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
+               ELSE 
+                  ns = 1
+               END IF
             ELSE
                CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns)
             END IF

src/pw_env_methods.F

@@ -757,7 +757,7 @@ SUBROUTINE compute_max_radius(radius, pw_env, qs_env)
       CHARACTER(LEN=8), DIMENSION(4), PARAMETER :: &
          pbas = (/"ORB     ", "AUX_FIT ", "MAO     ", "HARRIS  "/)
       CHARACTER(LEN=8), DIMENSION(8), PARAMETER :: sbas = (/"ORB     ", "AUX     ", "RI_AUX  ", &
-         "MAO     ", "HARRIS  ", "RI_HXC  ", "RI_K    ", "LRI_AUX "/)
+                                    "MAO     ", "HARRIS  ", "RI_HXC  ", "RI_K    ", "LRI_AUX "/)
       REAL(KIND=dp), PARAMETER                           :: safety_factor = 1.015_dp
 
       INTEGER :: handle, ibasis_set_type, igrid_level, igrid_zet0_s, ikind, ipgf, iset, ishell, &

src/qs_energy_utils.F

@@ -291,6 +291,10 @@ SUBROUTINE qs_energies_properties(qs_env)
          CALL xas(qs_env, dft_control)
       END IF
 
+      IF (dft_control%do_xas_tdp_calculation) THEN
+         CALL xas_tdp(qs_env)
+      END IF  
+
       ! Compute Linear Response properties as post-scf
       IF (.NOT. qs_env%linres_run) THEN
          CALL linres_calculation_low(qs_env)

src/qs_environment.F

@@ -506,7 +506,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
       TYPE(fist_nonbond_env_type), POINTER               :: se_nonbond_env
       TYPE(gapw_control_type), POINTER                   :: gapw_control
       TYPE(gto_basis_set_type), POINTER                  :: aux_fit_basis, lri_aux_basis, &
-                                                            ri_hfx_basis, tmp_basis_set
+                                                            ri_hfx_basis, tmp_basis_set, ri_xas_basis
       TYPE(lri_environment_type), POINTER                :: lri_env
       TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos, mos_aux_fit
       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
@@ -792,6 +792,22 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          END DO
       END IF
 
+      IF (dft_control%do_xas_tdp_calculation) THEN
+         ! Check if RI_XAS basis is given, auto-generate if not
+         CALL get_qs_env(qs_env, nkind=nkind)
+         DO ikind = 1, nkind
+            NULLIFY(ri_xas_basis)
+            qs_kind => qs_kind_set(ikind)
+            CALL get_qs_kind(qs_kind, basis_Set=ri_xas_basis, basis_type="RI_XAS")
+            IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN
+               CALL cp_warn(__LOCATION__, "Automatic Generation of RI_XAS basis. "// & 
+                            "This is experimental code.")
+               CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas)
+               CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS")
+            END IF
+         END DO
+      END IF
+
       !   *** Initialize the spherical harmonics and ***
       !   *** the orbital transformation matrices    ***
       CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl)
@@ -809,6 +825,11 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC")
          maxlgto = MAX(maxlgto, maxlgto_lri)
       END IF
+      IF (dft_control%do_xas_tdp_calculation) THEN
+         !done as a precaution
+         CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS")
+         maxlgto = MAX(maxlgto, maxlgto_lri)
+      END IF
       maxl = MAX(2*maxlgto, maxlppl, maxlppnl, lmax_sphere)+1
 
       CALL init_orbital_pointers(maxl)

src/qs_neighbor_lists.F

@@ -959,15 +959,17 @@ END SUBROUTINE build_qs_neighbor_lists
 !> \param current_subset ...
 !> \param operator_type ...
 !> \param nlname ...
+!> \param atomb_to_keep the list of atom indices to keep for pairs from the atom2d%b_list
 !> \date    20.03.2002
 !> \par History
 !>          - Major refactoring (25.07.2010,jhu)
+!>          - Added option to filter out atoms from list_b (08.2018, A.  Bussy)
 !> \author  MK
 !> \version 2.0
 ! **************************************************************************************************
    SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, &
                                    mic, symmetric, molecular, subset_of_mol, current_subset, &
-                                   operator_type, nlname)
+                                   operator_type, nlname, atomb_to_keep)
 
       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
          POINTER                                         :: ab_list
@@ -981,6 +983,7 @@ SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius,
       INTEGER, OPTIONAL                                  :: current_subset
       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: operator_type
       CHARACTER(LEN=*), INTENT(IN)                       :: nlname
+      INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: atomb_to_keep
 
       CHARACTER(len=*), PARAMETER :: routineN = 'build_neighbor_lists', &
          routineP = moduleN//':'//routineN
@@ -995,7 +998,7 @@ SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius,
                                                   nsubcell, periodic
       INTEGER, DIMENSION(:), POINTER           :: index_list
       LOGICAL                                  :: include_ab, my_mic, &
-                                                  my_molecular, my_symmetric
+                                                  my_molecular, my_symmetric, my_sort_atomb
       LOGICAL, ALLOCATABLE, DIMENSION(:)       :: pres_a, pres_b
       REAL(dp)                                 :: rab2, rab2_max, rab_max, rabm, deth, subcell_scale
       REAL(dp), DIMENSION(3)                   :: r, rab, ra, rb, sab_max, sb, &
@@ -1040,6 +1043,10 @@ SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius,
          ! default is a simple AB neighbor list
          otype = 1
       END IF
+      my_sort_atomb = .FALSE.
+      IF (PRESENT(atomb_to_keep)) THEN
+         my_sort_atomb = .TRUE.
+      END IF
 
       ! Deallocate the old neighbor list structure
       IF (ASSOCIATED(ab_list)) THEN
@@ -1204,6 +1211,9 @@ SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius,
 
             DO jatom_local = 1, nlistb(jkind)
                jatom = listb(jkind)%list(jatom_local)
+               IF (my_sort_atomb) THEN
+                  IF (.NOT. ANY(atomb_to_keep == jatom)) CYCLE
+               END IF
                atom_b = atom(jkind)%list(jatom)
                IF (my_molecular) THEN
                   mol_b = atom(jkind)%list_b_mol(jatom_local)

src/qs_o3c_methods.F

@@ -236,8 +236,8 @@ SUBROUTINE calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
                      END IF
 
                      DEALLOCATE (sabc_contr)
+                     DEALLOCATE (sabc)
                   END IF
-                  DEALLOCATE (sabc)
                   IF (do_force) THEN
                      DEALLOCATE (sdabc, sabdc)
                   END IF