Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Simplified Tamm Dancoff approach (sTDA).
9 : ! **************************************************************************************************
10 : MODULE qs_tddfpt2_stda_utils
11 :
12 : USE atomic_kind_types, ONLY: atomic_kind_type,&
13 : get_atomic_kind_set
14 : USE cell_types, ONLY: cell_type,&
15 : pbc
16 : USE cp_control_types, ONLY: stda_control_type,&
17 : tddfpt2_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, &
20 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
21 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_release, dbcsr_set, &
22 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
23 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag
24 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : copy_fm_to_dbcsr,&
27 : cp_dbcsr_plus_fm_fm_t,&
28 : cp_dbcsr_sm_fm_multiply,&
29 : dbcsr_allocate_matrix_set
30 : USE cp_fm_basic_linalg, ONLY: cp_fm_row_scale,&
31 : cp_fm_schur_product
32 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
33 : cp_fm_power
34 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 : cp_fm_struct_release,&
36 : cp_fm_struct_type
37 : USE cp_fm_types, ONLY: cp_fm_create,&
38 : cp_fm_get_info,&
39 : cp_fm_release,&
40 : cp_fm_set_all,&
41 : cp_fm_set_submatrix,&
42 : cp_fm_to_fm,&
43 : cp_fm_type,&
44 : cp_fm_vectorssum
45 : USE cp_log_handling, ONLY: cp_get_default_logger,&
46 : cp_logger_get_default_io_unit,&
47 : cp_logger_type
48 : USE ewald_environment_types, ONLY: ewald_env_create,&
49 : ewald_env_get,&
50 : ewald_env_set,&
51 : ewald_environment_type,&
52 : read_ewald_section_tb
53 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
54 : tb_spme_evaluate
55 : USE ewald_pw_types, ONLY: ewald_pw_create,&
56 : ewald_pw_type
57 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
58 : section_vals_type
59 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
60 : USE kinds, ONLY: dp
61 : USE mathconstants, ONLY: oorootpi
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE particle_methods, ONLY: get_particle_set
64 : USE particle_types, ONLY: particle_type
65 : USE qs_environment_types, ONLY: get_qs_env,&
66 : qs_environment_type
67 : USE qs_kind_types, ONLY: get_qs_kind_set,&
68 : qs_kind_type
69 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
70 : neighbor_list_iterate,&
71 : neighbor_list_iterator_create,&
72 : neighbor_list_iterator_p_type,&
73 : neighbor_list_iterator_release,&
74 : neighbor_list_set_p_type
75 : USE qs_tddfpt2_stda_types, ONLY: stda_env_type
76 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
77 : USE qs_tddfpt2_types, ONLY: tddfpt_work_matrices
78 : USE scf_control_types, ONLY: scf_control_type
79 : USE util, ONLY: get_limit
80 : USE virial_types, ONLY: virial_type
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_stda_utils'
88 :
89 : PUBLIC:: stda_init_matrices, stda_calculate_kernel, get_lowdin_x, get_lowdin_mo_coefficients, &
90 : setup_gamma
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief Calculate sTDA matrices
96 : !> \param qs_env ...
97 : !> \param stda_kernel ...
98 : !> \param sub_env ...
99 : !> \param work ...
100 : !> \param tddfpt_control ...
101 : ! **************************************************************************************************
102 408 : SUBROUTINE stda_init_matrices(qs_env, stda_kernel, sub_env, work, tddfpt_control)
103 :
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 : TYPE(stda_env_type) :: stda_kernel
106 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
107 : TYPE(tddfpt_work_matrices) :: work
108 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
109 :
110 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_init_matrices'
111 :
112 : INTEGER :: handle
113 : LOGICAL :: do_coulomb
114 : TYPE(cell_type), POINTER :: cell, cell_ref
115 : TYPE(ewald_environment_type), POINTER :: ewald_env
116 : TYPE(ewald_pw_type), POINTER :: ewald_pw
117 : TYPE(section_vals_type), POINTER :: ewald_section, poisson_section, &
118 : print_section
119 :
120 408 : CALL timeset(routineN, handle)
121 :
122 408 : do_coulomb = .NOT. tddfpt_control%rks_triplets
123 408 : IF (do_coulomb) THEN
124 : ! calculate exchange gamma matrix
125 314 : CALL setup_gamma(qs_env, stda_kernel, sub_env, work%gamma_exchange)
126 : END IF
127 :
128 : ! calculate S_half and Lowdin MO coefficients
129 408 : CALL get_lowdin_mo_coefficients(qs_env, sub_env, work)
130 :
131 : ! initialize Ewald for sTDA
132 408 : IF (tddfpt_control%stda_control%do_ewald) THEN
133 94 : NULLIFY (ewald_env, ewald_pw)
134 1692 : ALLOCATE (ewald_env)
135 94 : CALL ewald_env_create(ewald_env, sub_env%para_env)
136 94 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
137 94 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
138 94 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
139 94 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
140 94 : CALL get_qs_env(qs_env, cell=cell, cell_ref=cell_ref)
141 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat, &
142 94 : cell_periodic=cell%perd)
143 94 : ALLOCATE (ewald_pw)
144 94 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
145 94 : work%ewald_env => ewald_env
146 94 : work%ewald_pw => ewald_pw
147 : END IF
148 :
149 408 : CALL timestop(handle)
150 :
151 408 : END SUBROUTINE stda_init_matrices
152 : ! **************************************************************************************************
153 : !> \brief Calculate sTDA exchange-type contributions
154 : !> \param qs_env ...
155 : !> \param stda_env ...
156 : !> \param sub_env ...
157 : !> \param gamma_matrix sTDA exchange-type contributions
158 : !> \param ndim ...
159 : !> \note Note the specific sTDA notation exchange-type integrals (ia|jb) refer to Coulomb interaction
160 : ! **************************************************************************************************
161 484 : SUBROUTINE setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim)
162 :
163 : TYPE(qs_environment_type), POINTER :: qs_env
164 : TYPE(stda_env_type) :: stda_env
165 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
166 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: gamma_matrix
167 : INTEGER, INTENT(IN), OPTIONAL :: ndim
168 :
169 : CHARACTER(len=*), PARAMETER :: routineN = 'setup_gamma'
170 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
171 :
172 : INTEGER :: handle, i, iatom, icol, ikind, imat, &
173 : irow, jatom, jkind, natom, nmat
174 484 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
175 : LOGICAL :: found
176 : REAL(KIND=dp) :: dfcut, dgb, dr, eta, fcut, r, rcut, &
177 : rcuta, rcutb, x
178 : REAL(KIND=dp), DIMENSION(3) :: rij
179 484 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dgblock, gblock
180 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
181 : TYPE(neighbor_list_iterator_p_type), &
182 484 : DIMENSION(:), POINTER :: nl_iterator
183 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
184 484 : POINTER :: n_list
185 :
186 484 : CALL timeset(routineN, handle)
187 :
188 484 : CALL get_qs_env(qs_env=qs_env, natom=natom)
189 484 : dbcsr_dist => sub_env%dbcsr_dist
190 : ! Using the overlap list here can have a considerable effect on the number of
191 : ! terms calculated. This makes gamma also dependent on EPS_DEFAULT -> Overlap
192 484 : n_list => sub_env%sab_orb
193 :
194 484 : IF (PRESENT(ndim)) THEN
195 170 : nmat = ndim
196 : ELSE
197 314 : nmat = 1
198 : END IF
199 484 : CPASSERT(nmat == 1 .OR. nmat == 4)
200 484 : CPASSERT(.NOT. ASSOCIATED(gamma_matrix))
201 484 : CALL dbcsr_allocate_matrix_set(gamma_matrix, nmat)
202 :
203 1452 : ALLOCATE (row_blk_sizes(natom))
204 3120 : row_blk_sizes(1:natom) = 1
205 1478 : DO imat = 1, nmat
206 1478 : ALLOCATE (gamma_matrix(imat)%matrix)
207 : END DO
208 :
209 : CALL dbcsr_create(gamma_matrix(1)%matrix, name="gamma", dist=dbcsr_dist, &
210 : matrix_type=dbcsr_type_symmetric, row_blk_size=row_blk_sizes, &
211 484 : col_blk_size=row_blk_sizes)
212 994 : DO imat = 2, nmat
213 : CALL dbcsr_create(gamma_matrix(imat)%matrix, name="dgamma", dist=dbcsr_dist, &
214 : matrix_type=dbcsr_type_antisymmetric, row_blk_size=row_blk_sizes, &
215 994 : col_blk_size=row_blk_sizes)
216 : END DO
217 :
218 484 : DEALLOCATE (row_blk_sizes)
219 :
220 : ! setup the matrices using the neighbor list
221 1478 : DO imat = 1, nmat
222 994 : CALL cp_dbcsr_alloc_block_from_nbl(gamma_matrix(imat)%matrix, n_list)
223 1478 : CALL dbcsr_set(gamma_matrix(imat)%matrix, 0.0_dp)
224 : END DO
225 :
226 484 : NULLIFY (nl_iterator)
227 484 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
228 85654 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
230 85170 : iatom=iatom, jatom=jatom, r=rij)
231 :
232 340680 : dr = SQRT(SUM(rij(:)**2)) ! interatomic distance
233 :
234 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
235 85170 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
236 :
237 85170 : icol = MAX(iatom, jatom)
238 85170 : irow = MIN(iatom, jatom)
239 :
240 85170 : NULLIFY (gblock)
241 : CALL dbcsr_get_block_p(matrix=gamma_matrix(1)%matrix, &
242 85170 : row=irow, col=icol, BLOCK=gblock, found=found)
243 85170 : CPASSERT(found)
244 :
245 : ! get rcuta and rcutb
246 85170 : rcuta = stda_env%kind_param_set(ikind)%kind_param%rcut
247 85170 : rcutb = stda_env%kind_param_set(jkind)%kind_param%rcut
248 85170 : rcut = rcuta + rcutb
249 :
250 : !> Computes the short-range gamma parameter from
251 : !> Nataga-Mishimoto-Ohno-Klopman formula equivalently as it is done for xTB
252 85170 : IF (dr < 1.e-6) THEN
253 : ! on site terms
254 3954 : gblock(:, :) = gblock(:, :) + eta
255 83852 : ELSEIF (dr > rcut) THEN
256 : ! do nothing
257 : ELSE
258 40087 : IF (dr < rcut - rsmooth) THEN
259 : fcut = 1.0_dp
260 : ELSE
261 8759 : r = dr - (rcut - rsmooth)
262 8759 : x = r/rsmooth
263 8759 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
264 : END IF
265 : gblock(:, :) = gblock(:, :) + &
266 : fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
267 120261 : **(1._dp/stda_env%alpha_param) - fcut/dr
268 : END IF
269 :
270 85654 : IF (nmat > 1) THEN
271 : !> Computes the short-range gamma parameter from
272 : !> Nataga-Mishimoto-Ohno-Klopman formula equivalently as it is done for xTB
273 : !> Derivatives
274 16218 : IF (dr < 1.e-6 .OR. dr > rcut) THEN
275 : ! on site terms or beyond cutoff
276 : dgb = 0.0_dp
277 : ELSE
278 8257 : IF (dr < rcut - rsmooth) THEN
279 : fcut = 1.0_dp
280 : dfcut = 0.0_dp
281 : ELSE
282 1824 : r = dr - (rcut - rsmooth)
283 1824 : x = r/rsmooth
284 1824 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
285 1824 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
286 1824 : dfcut = dfcut/rsmooth
287 : END IF
288 : dgb = dfcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
289 8257 : **(1._dp/stda_env%alpha_param)
290 8257 : dgb = dgb - dfcut/dr + fcut/dr**2
291 : dgb = dgb - fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
292 8257 : **(1._dp/stda_env%alpha_param + 1._dp)*dr**(stda_env%alpha_param - 1._dp)
293 : END IF
294 64872 : DO imat = 2, nmat
295 48654 : NULLIFY (dgblock)
296 : CALL dbcsr_get_block_p(matrix=gamma_matrix(imat)%matrix, &
297 48654 : row=irow, col=icol, BLOCK=dgblock, found=found)
298 64872 : IF (found) THEN
299 48654 : IF (dr > 1.e-6) THEN
300 47592 : i = imat - 1
301 47592 : IF (irow == iatom) THEN
302 75357 : dgblock(:, :) = dgblock(:, :) + dgb*rij(i)/dr
303 : ELSE
304 67419 : dgblock(:, :) = dgblock(:, :) - dgb*rij(i)/dr
305 : END IF
306 : END IF
307 : END IF
308 : END DO
309 : END IF
310 :
311 : END DO
312 :
313 484 : CALL neighbor_list_iterator_release(nl_iterator)
314 :
315 1478 : DO imat = 1, nmat
316 1478 : CALL dbcsr_finalize(gamma_matrix(imat)%matrix)
317 : END DO
318 :
319 484 : CALL timestop(handle)
320 :
321 484 : END SUBROUTINE setup_gamma
322 :
323 : ! **************************************************************************************************
324 : !> \brief Calculate Lowdin MO coefficients
325 : !> \param qs_env ...
326 : !> \param sub_env ...
327 : !> \param work ...
328 : ! **************************************************************************************************
329 414 : SUBROUTINE get_lowdin_mo_coefficients(qs_env, sub_env, work)
330 :
331 : TYPE(qs_environment_type), POINTER :: qs_env
332 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
333 : TYPE(tddfpt_work_matrices) :: work
334 :
335 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_mo_coefficients'
336 :
337 : INTEGER :: handle, i, iounit, ispin, j, &
338 : max_iter_lanczos, nactive, ndep, nsgf, &
339 : nspins, order_lanczos
340 : LOGICAL :: converged
341 : REAL(KIND=dp) :: eps_lanczos, sij, threshold
342 414 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: slam
343 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
344 414 : POINTER :: local_data
345 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
346 : TYPE(cp_fm_type) :: fm_s_half, fm_work1
347 : TYPE(cp_logger_type), POINTER :: logger
348 414 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
349 : TYPE(dbcsr_type) :: sm_hinv
350 : TYPE(dbcsr_type), POINTER :: sm_h, sm_s
351 414 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
352 : TYPE(scf_control_type), POINTER :: scf_control
353 :
354 414 : CALL timeset(routineN, handle)
355 :
356 414 : NULLIFY (logger) !get output_unit
357 414 : logger => cp_get_default_logger()
358 414 : iounit = cp_logger_get_default_io_unit(logger)
359 :
360 : ! Calculate S^1/2 matrix
361 414 : IF (iounit > 0) THEN
362 207 : WRITE (iounit, "(1X,A)") "", &
363 207 : "-------------------------------------------------------------------------------", &
364 207 : "- Create Matrix SQRT(S) -", &
365 414 : "-------------------------------------------------------------------------------"
366 : END IF
367 :
368 414 : IF (sub_env%is_split) THEN
369 0 : CPABORT('SPLIT')
370 : ELSE
371 414 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrixkp_s)
372 414 : CPASSERT(ASSOCIATED(matrixkp_s))
373 414 : CPWARN_IF(SIZE(matrixkp_s, 2) > 1, "not implemented for k-points.")
374 414 : sm_s => matrixkp_s(1, 1)%matrix
375 : END IF
376 414 : sm_h => work%shalf
377 :
378 414 : CALL dbcsr_create(sm_hinv, template=sm_s)
379 414 : CALL dbcsr_add_on_diag(sm_h, 1.0_dp)
380 414 : threshold = 1.0e-8_dp
381 414 : order_lanczos = 3
382 414 : eps_lanczos = 1.0e-4_dp
383 414 : max_iter_lanczos = 40
384 : CALL matrix_sqrt_Newton_Schulz(sm_h, sm_hinv, sm_s, &
385 : threshold, order_lanczos, eps_lanczos, max_iter_lanczos, &
386 414 : converged=converged)
387 414 : CALL dbcsr_release(sm_hinv)
388 : !
389 414 : NULLIFY (qs_kind_set)
390 414 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
391 : ! Get the total number of contracted spherical Gaussian basis functions
392 414 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
393 : !
394 414 : IF (.NOT. converged) THEN
395 0 : IF (iounit > 0) THEN
396 0 : WRITE (iounit, "(T3,A)") "STDA| Newton-Schulz iteration did not converge"
397 0 : WRITE (iounit, "(T3,A)") "STDA| Calculate SQRT(S) from diagonalization"
398 : END IF
399 0 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control)
400 : ! Provide full size work matrices
401 : CALL cp_fm_struct_create(fmstruct=fmstruct, &
402 : para_env=sub_env%para_env, &
403 : context=sub_env%blacs_env, &
404 : nrow_global=nsgf, &
405 0 : ncol_global=nsgf)
406 0 : CALL cp_fm_create(matrix=fm_s_half, matrix_struct=fmstruct, name="S^(1/2) MATRIX")
407 0 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, name="TMP MATRIX")
408 0 : CALL cp_fm_struct_release(fmstruct=fmstruct)
409 0 : CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
410 0 : CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
411 0 : IF (ndep /= 0) &
412 : CALL cp_warn(__LOCATION__, &
413 : "Overlap matrix exhibits linear dependencies. At least some "// &
414 0 : "eigenvalues have been quenched.")
415 0 : CALL copy_fm_to_dbcsr(fm_s_half, sm_h)
416 0 : CALL cp_fm_release(fm_s_half)
417 0 : CALL cp_fm_release(fm_work1)
418 0 : IF (iounit > 0) WRITE (iounit, *)
419 : END IF
420 :
421 414 : nspins = SIZE(sub_env%mos_occ)
422 :
423 860 : DO ispin = 1, nspins
424 446 : CALL cp_fm_get_info(work%ctransformed(ispin), ncol_global=nactive)
425 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, sub_env%mos_active(ispin), &
426 860 : work%ctransformed(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
427 : END DO
428 :
429 : ! for Lowdin forces
430 414 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=work%S_eigenvectors%matrix_struct, name="TMP MATRIX")
431 414 : CALL copy_dbcsr_to_fm(sm_s, fm_work1)
432 414 : CALL choose_eigv_solver(fm_work1, work%S_eigenvectors, work%S_eigenvalues)
433 414 : CALL cp_fm_release(fm_work1)
434 : !
435 1242 : ALLOCATE (slam(nsgf, 1))
436 9690 : DO i = 1, nsgf
437 9690 : IF (work%S_eigenvalues(i) > 0._dp) THEN
438 9276 : slam(i, 1) = SQRT(work%S_eigenvalues(i))
439 : ELSE
440 0 : CPABORT("S matrix not positive definit")
441 : END IF
442 : END DO
443 9690 : DO i = 1, nsgf
444 9690 : CALL cp_fm_set_submatrix(work%slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
445 : END DO
446 9690 : DO i = 1, nsgf
447 9690 : CALL cp_fm_set_submatrix(work%slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .TRUE.)
448 : END DO
449 414 : CALL cp_fm_get_info(work%slambda, local_data=local_data)
450 9690 : DO i = 1, SIZE(local_data, 2)
451 422512 : DO j = 1, SIZE(local_data, 1)
452 412822 : sij = local_data(j, i)
453 412822 : IF (sij > 0.0_dp) sij = 1.0_dp/sij
454 422098 : local_data(j, i) = sij
455 : END DO
456 : END DO
457 414 : DEALLOCATE (slam)
458 :
459 414 : CALL timestop(handle)
460 :
461 1242 : END SUBROUTINE get_lowdin_mo_coefficients
462 :
463 : ! **************************************************************************************************
464 : !> \brief Calculate Lowdin transformed Davidson trial vector X
465 : !> shalf (dbcsr), xvec, xt (fm) are defined in the same sub_env
466 : !> \param shalf ...
467 : !> \param xvec ...
468 : !> \param xt ...
469 : ! **************************************************************************************************
470 6710 : SUBROUTINE get_lowdin_x(shalf, xvec, xt)
471 :
472 : TYPE(dbcsr_type), INTENT(IN) :: shalf
473 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: xvec
474 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: xt
475 :
476 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_x'
477 :
478 : INTEGER :: handle, ispin, nactive, nspins
479 :
480 6710 : CALL timeset(routineN, handle)
481 :
482 6710 : nspins = SIZE(xvec)
483 :
484 : ! Build Lowdin transformed tilde(X)= S^1/2 X for each spin
485 14116 : DO ispin = 1, nspins
486 7406 : CALL cp_fm_get_info(xt(ispin), ncol_global=nactive)
487 : CALL cp_dbcsr_sm_fm_multiply(shalf, xvec(ispin), &
488 14116 : xt(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
489 : END DO
490 :
491 6710 : CALL timestop(handle)
492 :
493 6710 : END SUBROUTINE get_lowdin_x
494 :
495 : ! **************************************************************************************************
496 : !> \brief ...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients --
497 : !> transition charges with the Coulomb-type or exchange-type integrals
498 : !> \param qs_env ...
499 : !> \param stda_control ...
500 : !> \param stda_env ...
501 : !> \param sub_env ...
502 : !> \param work ...
503 : !> \param is_rks_triplets ...
504 : !> \param X ...
505 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
506 : !> eigenvalue problem A*X = omega*X
507 : ! **************************************************************************************************
508 6478 : SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
509 6478 : work, is_rks_triplets, X, res)
510 :
511 : TYPE(qs_environment_type), POINTER :: qs_env
512 : TYPE(stda_control_type) :: stda_control
513 : TYPE(stda_env_type) :: stda_env
514 : TYPE(tddfpt_subgroup_env_type) :: sub_env
515 : TYPE(tddfpt_work_matrices) :: work
516 : LOGICAL, INTENT(IN) :: is_rks_triplets
517 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
518 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
519 :
520 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_calculate_kernel'
521 :
522 : INTEGER :: ewald_type, handle, ia, iatom, ikind, &
523 : is, ispin, jatom, jkind, jspin, natom, &
524 : nsgf, nspins
525 6478 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
526 : INTEGER, DIMENSION(2) :: nactive, nlim
527 : LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
528 : do_exchange, use_virial
529 : REAL(KIND=dp) :: alpha, bp, dr, eta, gabr, hfx, rbeta, &
530 : spinfac
531 6478 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
532 6478 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
533 : REAL(KIND=dp), DIMENSION(3) :: rij
534 6478 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gab, pblock
535 6478 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
536 : TYPE(cell_type), POINTER :: cell
537 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstructjspin
538 : TYPE(cp_fm_type) :: cvec, cvecjspin
539 6478 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
540 : TYPE(cp_fm_type), POINTER :: ct, ctjspin
541 : TYPE(dbcsr_iterator_type) :: iter
542 : TYPE(dbcsr_type) :: pdens
543 : TYPE(dbcsr_type), POINTER :: tempmat
544 : TYPE(ewald_environment_type), POINTER :: ewald_env
545 : TYPE(ewald_pw_type), POINTER :: ewald_pw
546 : TYPE(mp_para_env_type), POINTER :: para_env
547 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
548 6478 : POINTER :: n_list
549 6478 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
550 6478 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
551 : TYPE(virial_type), POINTER :: virial
552 :
553 6478 : CALL timeset(routineN, handle)
554 :
555 19434 : nactive(:) = stda_env%nactive(:)
556 6478 : nspins = SIZE(X)
557 6478 : spinfac = 2.0_dp
558 6478 : IF (nspins == 2) spinfac = 1.0_dp
559 :
560 5808 : IF (nspins == 1 .AND. is_rks_triplets) THEN
561 : do_coulomb = .FALSE.
562 : ELSE
563 : do_coulomb = .TRUE.
564 : END IF
565 6478 : do_ewald = stda_control%do_ewald
566 6478 : do_exchange = stda_control%do_exchange
567 :
568 6478 : para_env => sub_env%para_env
569 :
570 : CALL get_qs_env(qs_env, natom=natom, cell=cell, &
571 6478 : particle_set=particle_set, qs_kind_set=qs_kind_set)
572 19434 : ALLOCATE (first_sgf(natom))
573 12956 : ALLOCATE (last_sgf(natom))
574 6478 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
575 :
576 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
577 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
578 26582 : ALLOCATE (xtransformed(nspins))
579 13626 : DO ispin = 1, nspins
580 7148 : NULLIFY (fmstruct)
581 7148 : ct => work%ctransformed(ispin)
582 7148 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
583 13626 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
584 : END DO
585 6478 : CALL get_lowdin_x(work%shalf, X, xtransformed)
586 :
587 25912 : ALLOCATE (tcharge(natom), gtcharge(natom, 1))
588 :
589 13626 : DO ispin = 1, nspins
590 13626 : CALL cp_fm_set_all(res(ispin), 0.0_dp)
591 : END DO
592 :
593 13626 : DO ispin = 1, nspins
594 7148 : ct => work%ctransformed(ispin)
595 7148 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
596 21444 : ALLOCATE (tv(nsgf))
597 7148 : CALL cp_fm_create(cvec, fmstruct)
598 : !
599 : ! *** Coulomb contribution
600 : !
601 7148 : IF (do_coulomb) THEN
602 57080 : tcharge(:) = 0.0_dp
603 12520 : DO jspin = 1, nspins
604 6930 : ctjspin => work%ctransformed(jspin)
605 6930 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
606 6930 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
607 6930 : CALL cp_fm_create(cvecjspin, fmstructjspin)
608 : ! CV(mu,j) = CT(mu,j)*XT(mu,j)
609 6930 : CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
610 : ! TV(mu) = SUM_j CV(mu,j)
611 6930 : CALL cp_fm_vectorssum(cvecjspin, tv, "R")
612 : ! contract charges
613 : ! TC(a) = SUM_(mu of a) TV(mu)
614 62488 : DO ia = 1, natom
615 271930 : DO is = first_sgf(ia), last_sgf(ia)
616 265000 : tcharge(ia) = tcharge(ia) + tv(is)
617 : END DO
618 : END DO
619 19450 : CALL cp_fm_release(cvecjspin)
620 : END DO !jspin
621 : ! Apply tcharge*gab -> gtcharge
622 : ! gT(b) = SUM_a g(a,b)*TC(a)
623 : ! gab = work%gamma_exchange(1)%matrix
624 62670 : gtcharge = 0.0_dp
625 : ! short range contribution
626 5590 : tempmat => work%gamma_exchange(1)%matrix
627 5590 : CALL dbcsr_iterator_start(iter, tempmat)
628 878993 : DO WHILE (dbcsr_iterator_blocks_left(iter))
629 873403 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab)
630 873403 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
631 878993 : IF (iatom /= jatom) THEN
632 847658 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
633 : END IF
634 : END DO
635 5590 : CALL dbcsr_iterator_stop(iter)
636 : ! Ewald long range contribution
637 5590 : IF (do_ewald) THEN
638 740 : ewald_env => work%ewald_env
639 740 : ewald_pw => work%ewald_pw
640 740 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
641 740 : CALL get_qs_env(qs_env=qs_env, virial=virial)
642 740 : use_virial = .FALSE.
643 740 : calculate_forces = .FALSE.
644 740 : n_list => sub_env%sab_orb
645 740 : CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
646 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
647 740 : gtcharge, tcharge, calculate_forces, virial, use_virial)
648 : ! add self charge interaction contribution
649 740 : IF (para_env%is_source()) THEN
650 18724 : gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
651 : END IF
652 : ELSE
653 4850 : nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
654 12241 : DO iatom = nlim(1), nlim(2)
655 19864 : DO jatom = 1, iatom - 1
656 30492 : rij = particle_set(iatom)%r - particle_set(jatom)%r
657 30492 : rij = pbc(rij, cell)
658 30492 : dr = SQRT(SUM(rij(:)**2))
659 15014 : IF (dr > 1.e-6_dp) THEN
660 7623 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
661 7623 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
662 : END IF
663 : END DO
664 : END DO
665 : END IF
666 5590 : CALL para_env%sum(gtcharge)
667 : ! expand charges
668 : ! TV(mu) = TC(a of mu)
669 186372 : tv(1:nsgf) = 0.0_dp
670 57080 : DO ia = 1, natom
671 237862 : DO is = first_sgf(ia), last_sgf(ia)
672 232272 : tv(is) = gtcharge(ia, 1)
673 : END DO
674 : END DO
675 : ! CV(mu,i) = TV(mu)*CV(mu,i)
676 5590 : ct => work%ctransformed(ispin)
677 5590 : CALL cp_fm_to_fm(ct, cvec)
678 5590 : CALL cp_fm_row_scale(cvec, tv)
679 : ! rho(nu,i) = rho(nu,i) + Shalf(nu,mu)*CV(mu,i)
680 5590 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), spinfac, 1.0_dp)
681 : END IF
682 : !
683 : ! *** Exchange contribution
684 : !
685 7148 : IF (do_exchange) THEN ! option to explicitly switch off exchange
686 : ! (exchange contributes also if hfx_fraction=0)
687 6704 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
688 6704 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
689 : !
690 6704 : tempmat => work%shalf
691 6704 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
692 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
693 6704 : ct => work%ctransformed(ispin)
694 6704 : CALL dbcsr_set(pdens, 0.0_dp)
695 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
696 6704 : 1.0_dp, keep_sparsity=.FALSE.)
697 6704 : CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
698 : ! Apply PP*gab -> PP; gab = gamma_coulomb
699 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
700 6704 : bp = stda_env%beta_param
701 6704 : hfx = stda_env%hfx_fraction
702 6704 : CALL dbcsr_iterator_start(iter, pdens)
703 1738129 : DO WHILE (dbcsr_iterator_blocks_left(iter))
704 1731425 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
705 6925700 : rij = particle_set(iatom)%r - particle_set(jatom)%r
706 6925700 : rij = pbc(rij, cell)
707 6925700 : dr = SQRT(SUM(rij(:)**2))
708 1731425 : ikind = kind_of(iatom)
709 1731425 : jkind = kind_of(jatom)
710 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
711 1731425 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
712 1731425 : rbeta = dr**bp
713 1731425 : IF (hfx > 0.0_dp) THEN
714 1729225 : gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
715 : ELSE
716 2200 : IF (dr < 1.e-6) THEN
717 : gabr = 0.0_dp
718 : ELSE
719 1524 : gabr = 1._dp/dr
720 : END IF
721 : END IF
722 19425165 : pblock = gabr*pblock
723 : END DO
724 6704 : CALL dbcsr_iterator_stop(iter)
725 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
726 6704 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
727 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
728 6704 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), -1.0_dp, 1.0_dp)
729 : !
730 6704 : CALL dbcsr_release(pdens)
731 6704 : DEALLOCATE (kind_of)
732 : END IF
733 : !
734 7148 : CALL cp_fm_release(cvec)
735 27922 : DEALLOCATE (tv)
736 : END DO
737 :
738 6478 : CALL cp_fm_release(xtransformed)
739 6478 : DEALLOCATE (tcharge, gtcharge)
740 6478 : DEALLOCATE (first_sgf, last_sgf)
741 :
742 6478 : CALL timestop(handle)
743 :
744 12956 : END SUBROUTINE stda_calculate_kernel
745 :
746 : ! **************************************************************************************************
747 :
748 : END MODULE qs_tddfpt2_stda_utils
|
