Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE qs_tddfpt2_lri_utils
9 :
10 : USE atomic_kind_types, ONLY: atomic_kind_type
11 : USE auto_basis, ONLY: create_lri_aux_basis_set
12 : USE basis_set_container_types, ONLY: add_basis_set_to_container
13 : USE basis_set_output, ONLY: print_basis_set_file
14 : USE basis_set_types, ONLY: get_gto_basis_set,&
15 : gto_basis_set_type,&
16 : init_orb_basis_set,&
17 : sort_gto_basis_set
18 : USE cell_types, ONLY: cell_type
19 : USE cp_control_types, ONLY: dft_control_type,&
20 : tddfpt2_control_type
21 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
22 : USE dbcsr_api, ONLY: dbcsr_copy,&
23 : dbcsr_deallocate_matrix,&
24 : dbcsr_p_type
25 : USE distribution_1d_types, ONLY: distribution_1d_type
26 : USE distribution_2d_types, ONLY: distribution_2d_type
27 : USE input_section_types, ONLY: section_vals_get,&
28 : section_vals_get_subs_vals,&
29 : section_vals_type,&
30 : section_vals_val_get
31 : USE kinds, ONLY: dp
32 : USE lri_environment_init, ONLY: lri_env_basis,&
33 : lri_env_init
34 : USE lri_environment_methods, ONLY: build_lri_matrices
35 : USE lri_environment_types, ONLY: lri_environment_type,&
36 : lri_kind_type
37 : USE lri_ks_methods, ONLY: calculate_lri_ks_matrix
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE molecule_types, ONLY: molecule_type
40 : USE orbital_pointers, ONLY: init_orbital_pointers
41 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
42 : USE particle_types, ONLY: particle_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
46 : USE qs_kernel_types, ONLY: kernel_env_type
47 : USE qs_kind_types, ONLY: get_qs_kind,&
48 : get_qs_kind_set,&
49 : qs_kind_type
50 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
51 : USE qs_neighbor_lists, ONLY: atom2d_build,&
52 : atom2d_cleanup,&
53 : build_neighbor_lists,&
54 : local_atoms_type,&
55 : pair_radius_setup,&
56 : write_neighbor_lists
57 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_dbcsr_create_by_dist,&
58 : tddfpt_subgroup_env_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_lri_utils'
66 :
67 : PUBLIC:: tddfpt2_lri_init, tddfpt2_lri_Amat
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Initialize LRI environment, basis, neighborlists and matrices
73 : !> \param qs_env ...
74 : !> \param kernel_env ...
75 : !> \param lri_section ...
76 : !> \param tddfpt_print_section ...
77 : ! **************************************************************************************************
78 10 : SUBROUTINE tddfpt2_lri_init(qs_env, kernel_env, lri_section, tddfpt_print_section)
79 :
80 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
81 : TYPE(kernel_env_type) :: kernel_env
82 : TYPE(section_vals_type), INTENT(IN), POINTER :: lri_section, tddfpt_print_section
83 :
84 : CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt2_lri_init'
85 :
86 : INTEGER :: basis_sort, handle, ikind, lmax_sphere, &
87 : maxl, maxlgto, maxlgto_lri, nkind
88 : LOGICAL :: explicit, mic, molecule_only, &
89 : redefine_interaction_radii
90 10 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_present
91 : REAL(dp) :: subcells
92 10 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: orb_radius
93 10 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
94 10 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
95 : TYPE(cell_type), POINTER :: cell
96 : TYPE(dft_control_type), POINTER :: dft_control
97 : TYPE(distribution_1d_type), POINTER :: distribution_1d
98 : TYPE(distribution_2d_type), POINTER :: distribution_2d
99 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set, p_lri_aux_basis
100 10 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
101 : TYPE(lri_environment_type), POINTER :: lri_env
102 10 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
103 : TYPE(mp_para_env_type), POINTER :: para_env
104 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
105 10 : POINTER :: soo_list
106 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
108 : TYPE(qs_kind_type), POINTER :: qs_kind
109 : TYPE(section_vals_type), POINTER :: neighbor_list_section, print_section
110 : TYPE(tddfpt2_control_type), POINTER :: tddfpt2_control
111 :
112 10 : CALL timeset(routineN, handle)
113 :
114 10 : CALL get_qs_env(qs_env, dft_control=dft_control)
115 10 : tddfpt2_control => dft_control%tddfpt2_control
116 :
117 10 : NULLIFY (kernel_env%full_kernel%lri_env)
118 : ! initialize lri_env
119 10 : CALL lri_env_init(kernel_env%full_kernel%lri_env, lri_section)
120 : NULLIFY (lri_env)
121 10 : lri_env => kernel_env%full_kernel%lri_env
122 10 : redefine_interaction_radii = .FALSE.
123 :
124 : ! exact_1c_terms not implemented
125 10 : IF (lri_env%exact_1c_terms) THEN
126 0 : CPABORT("TDDFT with LRI and exact one-center terms not implemented")
127 : END IF
128 :
129 : ! check if LRI_AUX basis is available
130 10 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
131 10 : nkind = SIZE(qs_kind_set)
132 30 : DO ikind = 1, nkind
133 20 : NULLIFY (p_lri_aux_basis)
134 20 : qs_kind => qs_kind_set(ikind)
135 20 : CALL get_qs_kind(qs_kind, basis_set=p_lri_aux_basis, basis_type="P_LRI_AUX")
136 30 : IF (.NOT. (ASSOCIATED(p_lri_aux_basis))) THEN
137 : ! Generate a default basis
138 12 : redefine_interaction_radii = .TRUE.
139 : CALL create_lri_aux_basis_set(p_lri_aux_basis, qs_kind, &
140 : dft_control%auto_basis_p_lri_aux, lri_env%exact_1c_terms, &
141 12 : tda_kernel=.TRUE.)
142 12 : CALL add_basis_set_to_container(qs_kind%basis_sets, p_lri_aux_basis, "P_LRI_AUX")
143 : END IF
144 : END DO !nkind
145 : ! needs to be done here if p_lri_aux_basis is not specified explicitly
146 10 : IF (redefine_interaction_radii) THEN
147 18 : DO ikind = 1, nkind
148 12 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=p_lri_aux_basis, basis_type="P_LRI_AUX")
149 18 : IF (ASSOCIATED(p_lri_aux_basis)) THEN
150 12 : CALL init_orb_basis_set(p_lri_aux_basis) ! standardly done in init_qs_kind
151 12 : basis_sort = 0
152 12 : CALL sort_gto_basis_set(p_lri_aux_basis, basis_sort)
153 12 : CALL init_interaction_radii_orb_basis(p_lri_aux_basis, dft_control%qs_control%eps_pgf_orb)
154 : END IF
155 : END DO
156 : END IF
157 :
158 10 : print_section => section_vals_get_subs_vals(tddfpt_print_section, "BASIS_SET_FILE")
159 10 : CALL section_vals_get(print_section, explicit=explicit)
160 10 : IF (explicit) THEN
161 0 : CALL print_basis_set_file(qs_env, print_section)
162 : END IF
163 :
164 : !set maxl as in qs_environment for gs lri
165 10 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
166 10 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="P_LRI_AUX")
167 : !take maxlgto from lri basis if larger (usually)
168 10 : maxlgto = MAX(maxlgto, maxlgto_lri)
169 10 : lmax_sphere = 2*maxlgto
170 10 : maxl = MAX(2*maxlgto, lmax_sphere) + 1
171 :
172 10 : CALL init_orbital_pointers(maxl)
173 10 : CALL init_spherical_harmonics(maxl, 0)
174 :
175 : ! initialize LRI basis sets
176 10 : CALL lri_env_basis("P_LRI", qs_env, lri_env, qs_kind_set)
177 :
178 : ! check for debugging that automatically generated basis is available
179 : ! DO ikind= 1, nkind
180 : ! qs_kind => qs_kind_set(ikind)
181 : ! CALL get_qs_kind(qs_kind, basis_set=p_lri_aux_basis, basis_type="P_LRI_AUX")
182 : ! CALL get_gto_basis_set(gto_basis_set=p_lri_aux_basis, set_radius=set_radius,&
183 : ! pgf_radius=pgf_radius)
184 : ! END DO !nkind
185 :
186 : ! set up LRI neighbor list soo_list => same as in qs_neighbor_lists for ground-state LRI
187 10 : NULLIFY (cell, para_env, particle_set)
188 10 : CALL get_qs_env(qs_env, para_env=para_env, cell=cell, particle_set=particle_set)
189 :
190 10 : NULLIFY (distribution_1d, distribution_2d, atomic_kind_set, molecule_set)
191 : CALL get_qs_env(qs_env, local_particles=distribution_1d, distribution_2d=distribution_2d, &
192 10 : atomic_kind_set=atomic_kind_set, molecule_set=molecule_set)
193 :
194 30 : ALLOCATE (atom2d(nkind))
195 10 : molecule_only = .FALSE. !this still needs to be checked
196 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
197 10 : molecule_set, molecule_only, particle_set=particle_set)
198 :
199 80 : ALLOCATE (orb_present(nkind), orb_radius(nkind), pair_radius(nkind, nkind))
200 30 : orb_radius(:) = 0.0_dp
201 70 : pair_radius(:, :) = 0.0_dp
202 30 : orb_present(:) = .FALSE.
203 30 : DO ikind = 1, nkind
204 20 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
205 30 : IF (ASSOCIATED(orb_basis_set)) THEN
206 20 : orb_present(ikind) = .TRUE.
207 20 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
208 : ELSE
209 0 : orb_present(ikind) = .FALSE.
210 : END IF
211 : END DO ! ikind
212 :
213 10 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
214 :
215 : NULLIFY (soo_list)
216 10 : soo_list => lri_env%soo_list
217 10 : mic = .TRUE. !enforcing minimum image convention
218 10 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
219 : CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
220 10 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
221 :
222 : ! make this a TDDFT neighborlist
223 10 : neighbor_list_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%NEIGHBOR_LISTS")
224 : CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, &
225 10 : "/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)")
226 10 : lri_env%soo_list => soo_list
227 :
228 10 : CALL atom2d_cleanup(atom2d)
229 :
230 10 : DEALLOCATE (orb_present, orb_radius, pair_radius)
231 :
232 : ! calculate LRI integrals
233 10 : lri_env%ppl_ri = .FALSE. ! make sure that option is not available for ES
234 10 : CALL build_lri_matrices(lri_env, qs_env)
235 10 : kernel_env%full_kernel%lri_env => lri_env
236 :
237 : ! CALL get_condition_number_of_overlap(lri_env)
238 :
239 10 : CALL timestop(handle)
240 :
241 50 : END SUBROUTINE tddfpt2_lri_init
242 : ! **************************************************************************************************
243 : !> \brief Calculate contribution to response vector for LRI
244 : !> \param qs_env ...
245 : !> \param sub_env ...
246 : !> \param lri_env ...
247 : !> \param lri_v_int ...
248 : !> \param A_ia_munu_sub ...
249 : ! **************************************************************************************************
250 206 : SUBROUTINE tddfpt2_lri_Amat(qs_env, sub_env, lri_env, lri_v_int, A_ia_munu_sub)
251 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
252 : TYPE(tddfpt_subgroup_env_type), INTENT(IN) :: sub_env
253 : TYPE(lri_environment_type), INTENT(IN), POINTER :: lri_env
254 : TYPE(lri_kind_type), DIMENSION(:), INTENT(IN), &
255 : POINTER :: lri_v_int
256 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: A_ia_munu_sub
257 :
258 : CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt2_lri_Amat'
259 :
260 : INTEGER :: handle, ispin, nspins
261 206 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
262 206 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dummymat, matrix_s
263 :
264 206 : CALL timeset(routineN, handle)
265 206 : NULLIFY (atomic_kind_set)
266 206 : NULLIFY (matrix_s)
267 206 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, matrix_s=matrix_s)
268 206 : nspins = SIZE(A_ia_munu_sub)
269 412 : DO ispin = 1, nspins!
270 : !no kpoints for excited states => using dummy matrix with no cell index
271 206 : NULLIFY (dummymat)
272 206 : CALL dbcsr_allocate_matrix_set(dummymat, 1)
273 : CALL tddfpt_dbcsr_create_by_dist(dummymat(1)%matrix, template=matrix_s(1)%matrix, &
274 206 : dbcsr_dist=sub_env%dbcsr_dist, sab=sub_env%sab_orb)
275 :
276 206 : CALL dbcsr_copy(A_ia_munu_sub(ispin)%matrix, dummymat(1)%matrix)
277 :
278 206 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, dummymat, atomic_kind_set)
279 :
280 206 : CALL dbcsr_copy(A_ia_munu_sub(ispin)%matrix, dummymat(1)%matrix)
281 :
282 206 : CALL dbcsr_deallocate_matrix(dummymat(1)%matrix)
283 412 : DEALLOCATE (dummymat)
284 : END DO
285 206 : CALL timestop(handle)
286 :
287 206 : END SUBROUTINE tddfpt2_lri_Amat
288 : ! **************************************************************************************************
289 : END MODULE qs_tddfpt2_lri_utils
|
