Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 402 : 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 402 : CALL timeset(routineN, handle)
121 :
122 402 : do_coulomb = .NOT. tddfpt_control%rks_triplets
123 402 : IF (do_coulomb) THEN
124 : ! calculate exchange gamma matrix
125 308 : 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 402 : CALL get_lowdin_mo_coefficients(qs_env, sub_env, work)
130 :
131 : ! initialize Ewald for sTDA
132 402 : 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 94 : CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat)
142 94 : ALLOCATE (ewald_pw)
143 94 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
144 94 : work%ewald_env => ewald_env
145 94 : work%ewald_pw => ewald_pw
146 : END IF
147 :
148 402 : CALL timestop(handle)
149 :
150 402 : END SUBROUTINE stda_init_matrices
151 : ! **************************************************************************************************
152 : !> \brief Calculate sTDA exchange-type contributions
153 : !> \param qs_env ...
154 : !> \param stda_env ...
155 : !> \param sub_env ...
156 : !> \param gamma_matrix sTDA exchange-type contributions
157 : !> \param ndim ...
158 : !> \note Note the specific sTDA notation exchange-type integrals (ia|jb) refer to Coulomb interaction
159 : ! **************************************************************************************************
160 476 : SUBROUTINE setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim)
161 :
162 : TYPE(qs_environment_type), POINTER :: qs_env
163 : TYPE(stda_env_type) :: stda_env
164 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
165 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: gamma_matrix
166 : INTEGER, INTENT(IN), OPTIONAL :: ndim
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'setup_gamma'
169 : REAL(KIND=dp), PARAMETER :: rsmooth = 1.0_dp
170 :
171 : INTEGER :: handle, i, iatom, icol, ikind, imat, &
172 : irow, jatom, jkind, natom, nmat
173 476 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
174 : LOGICAL :: found
175 : REAL(KIND=dp) :: dfcut, dgb, dr, eta, fcut, r, rcut, &
176 : rcuta, rcutb, x
177 : REAL(KIND=dp), DIMENSION(3) :: rij
178 476 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dgblock, gblock
179 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
180 : TYPE(neighbor_list_iterator_p_type), &
181 476 : DIMENSION(:), POINTER :: nl_iterator
182 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
183 476 : POINTER :: n_list
184 :
185 476 : CALL timeset(routineN, handle)
186 :
187 476 : CALL get_qs_env(qs_env=qs_env, natom=natom)
188 476 : dbcsr_dist => sub_env%dbcsr_dist
189 : ! Using the overlap list here can have a considerable effect on the number of
190 : ! terms calculated. This makes gamma also dependent on EPS_DEFAULT -> Overlap
191 476 : n_list => sub_env%sab_orb
192 :
193 476 : IF (PRESENT(ndim)) THEN
194 168 : nmat = ndim
195 : ELSE
196 308 : nmat = 1
197 : END IF
198 476 : CPASSERT(nmat == 1 .OR. nmat == 4)
199 476 : CPASSERT(.NOT. ASSOCIATED(gamma_matrix))
200 476 : CALL dbcsr_allocate_matrix_set(gamma_matrix, nmat)
201 :
202 1428 : ALLOCATE (row_blk_sizes(natom))
203 3088 : row_blk_sizes(1:natom) = 1
204 1456 : DO imat = 1, nmat
205 1456 : ALLOCATE (gamma_matrix(imat)%matrix)
206 : END DO
207 :
208 : CALL dbcsr_create(gamma_matrix(1)%matrix, name="gamma", dist=dbcsr_dist, &
209 : matrix_type=dbcsr_type_symmetric, row_blk_size=row_blk_sizes, &
210 476 : col_blk_size=row_blk_sizes)
211 980 : DO imat = 2, nmat
212 : CALL dbcsr_create(gamma_matrix(imat)%matrix, name="dgamma", dist=dbcsr_dist, &
213 : matrix_type=dbcsr_type_antisymmetric, row_blk_size=row_blk_sizes, &
214 980 : col_blk_size=row_blk_sizes)
215 : END DO
216 :
217 476 : DEALLOCATE (row_blk_sizes)
218 :
219 : ! setup the matrices using the neighbor list
220 1456 : DO imat = 1, nmat
221 980 : CALL cp_dbcsr_alloc_block_from_nbl(gamma_matrix(imat)%matrix, n_list)
222 1456 : CALL dbcsr_set(gamma_matrix(imat)%matrix, 0.0_dp)
223 : END DO
224 :
225 476 : NULLIFY (nl_iterator)
226 476 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
227 85622 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
228 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
229 85146 : iatom=iatom, jatom=jatom, r=rij)
230 :
231 340584 : dr = SQRT(SUM(rij(:)**2)) ! interatomic distance
232 :
233 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
234 85146 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
235 :
236 85146 : icol = MAX(iatom, jatom)
237 85146 : irow = MIN(iatom, jatom)
238 :
239 85146 : NULLIFY (gblock)
240 : CALL dbcsr_get_block_p(matrix=gamma_matrix(1)%matrix, &
241 85146 : row=irow, col=icol, BLOCK=gblock, found=found)
242 85146 : CPASSERT(found)
243 :
244 : ! get rcuta and rcutb
245 85146 : rcuta = stda_env%kind_param_set(ikind)%kind_param%rcut
246 85146 : rcutb = stda_env%kind_param_set(jkind)%kind_param%rcut
247 85146 : rcut = rcuta + rcutb
248 :
249 : !> Computes the short-range gamma parameter from
250 : !> Nataga-Mishimoto-Ohno-Klopman formula equivalently as it is done for xTB
251 85146 : IF (dr < 1.e-6) THEN
252 : ! on site terms
253 3918 : gblock(:, :) = gblock(:, :) + eta
254 83840 : ELSEIF (dr > rcut) THEN
255 : ! do nothing
256 : ELSE
257 40075 : IF (dr < rcut - rsmooth) THEN
258 : fcut = 1.0_dp
259 : ELSE
260 8759 : r = dr - (rcut - rsmooth)
261 8759 : x = r/rsmooth
262 8759 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
263 : END IF
264 : gblock(:, :) = gblock(:, :) + &
265 : fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
266 120225 : **(1._dp/stda_env%alpha_param) - fcut/dr
267 : END IF
268 :
269 85622 : IF (nmat > 1) THEN
270 : !> Computes the short-range gamma parameter from
271 : !> Nataga-Mishimoto-Ohno-Klopman formula equivalently as it is done for xTB
272 : !> Derivatives
273 16212 : IF (dr < 1.e-6 .OR. dr > rcut) THEN
274 : ! on site terms or beyond cutoff
275 : dgb = 0.0_dp
276 : ELSE
277 8254 : IF (dr < rcut - rsmooth) THEN
278 : fcut = 1.0_dp
279 : dfcut = 0.0_dp
280 : ELSE
281 1824 : r = dr - (rcut - rsmooth)
282 1824 : x = r/rsmooth
283 1824 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
284 1824 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
285 1824 : dfcut = dfcut/rsmooth
286 : END IF
287 : dgb = dfcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
288 8254 : **(1._dp/stda_env%alpha_param)
289 8254 : dgb = dgb - dfcut/dr + fcut/dr**2
290 : dgb = dgb - fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
291 8254 : **(1._dp/stda_env%alpha_param + 1._dp)*dr**(stda_env%alpha_param - 1._dp)
292 : END IF
293 64848 : DO imat = 2, nmat
294 48636 : NULLIFY (dgblock)
295 : CALL dbcsr_get_block_p(matrix=gamma_matrix(imat)%matrix, &
296 48636 : row=irow, col=icol, BLOCK=dgblock, found=found)
297 64848 : IF (found) THEN
298 48636 : IF (dr > 1.e-6) THEN
299 47583 : i = imat - 1
300 47583 : IF (irow == iatom) THEN
301 75348 : dgblock(:, :) = dgblock(:, :) + dgb*rij(i)/dr
302 : ELSE
303 67401 : dgblock(:, :) = dgblock(:, :) - dgb*rij(i)/dr
304 : END IF
305 : END IF
306 : END IF
307 : END DO
308 : END IF
309 :
310 : END DO
311 :
312 476 : CALL neighbor_list_iterator_release(nl_iterator)
313 :
314 1456 : DO imat = 1, nmat
315 1456 : CALL dbcsr_finalize(gamma_matrix(imat)%matrix)
316 : END DO
317 :
318 476 : CALL timestop(handle)
319 :
320 476 : END SUBROUTINE setup_gamma
321 :
322 : ! **************************************************************************************************
323 : !> \brief Calculate Lowdin MO coefficients
324 : !> \param qs_env ...
325 : !> \param sub_env ...
326 : !> \param work ...
327 : ! **************************************************************************************************
328 408 : SUBROUTINE get_lowdin_mo_coefficients(qs_env, sub_env, work)
329 :
330 : TYPE(qs_environment_type), POINTER :: qs_env
331 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
332 : TYPE(tddfpt_work_matrices) :: work
333 :
334 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_mo_coefficients'
335 :
336 : INTEGER :: handle, i, iounit, ispin, j, &
337 : max_iter_lanczos, nactive, ndep, nsgf, &
338 : nspins, order_lanczos
339 : LOGICAL :: converged
340 : REAL(KIND=dp) :: eps_lanczos, sij, threshold
341 408 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: slam
342 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
343 408 : POINTER :: local_data
344 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
345 : TYPE(cp_fm_type) :: fm_s_half, fm_work1
346 : TYPE(cp_logger_type), POINTER :: logger
347 408 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
348 : TYPE(dbcsr_type) :: sm_hinv
349 : TYPE(dbcsr_type), POINTER :: sm_h, sm_s
350 408 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
351 : TYPE(scf_control_type), POINTER :: scf_control
352 :
353 408 : CALL timeset(routineN, handle)
354 :
355 408 : NULLIFY (logger) !get output_unit
356 408 : logger => cp_get_default_logger()
357 408 : iounit = cp_logger_get_default_io_unit(logger)
358 :
359 : ! Calculate S^1/2 matrix
360 408 : IF (iounit > 0) THEN
361 204 : WRITE (iounit, "(1X,A)") "", &
362 204 : "-------------------------------------------------------------------------------", &
363 204 : "- Create Matrix SQRT(S) -", &
364 408 : "-------------------------------------------------------------------------------"
365 : END IF
366 :
367 408 : IF (sub_env%is_split) THEN
368 0 : CPABORT('SPLIT')
369 : ELSE
370 408 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrixkp_s)
371 408 : CPASSERT(ASSOCIATED(matrixkp_s))
372 408 : CPWARN_IF(SIZE(matrixkp_s, 2) > 1, "not implemented for k-points.")
373 408 : sm_s => matrixkp_s(1, 1)%matrix
374 : END IF
375 408 : sm_h => work%shalf
376 :
377 408 : CALL dbcsr_create(sm_hinv, template=sm_s)
378 408 : CALL dbcsr_add_on_diag(sm_h, 1.0_dp)
379 408 : threshold = 1.0e-8_dp
380 408 : order_lanczos = 3
381 408 : eps_lanczos = 1.0e-4_dp
382 408 : max_iter_lanczos = 40
383 : CALL matrix_sqrt_Newton_Schulz(sm_h, sm_hinv, sm_s, &
384 : threshold, order_lanczos, eps_lanczos, max_iter_lanczos, &
385 408 : converged=converged)
386 408 : CALL dbcsr_release(sm_hinv)
387 : !
388 408 : NULLIFY (qs_kind_set)
389 408 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
390 : ! Get the total number of contracted spherical Gaussian basis functions
391 408 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
392 : !
393 408 : IF (.NOT. converged) THEN
394 0 : IF (iounit > 0) THEN
395 0 : WRITE (iounit, "(T3,A)") "STDA| Newton-Schulz iteration did not converge"
396 0 : WRITE (iounit, "(T3,A)") "STDA| Calculate SQRT(S) from diagonalization"
397 : END IF
398 0 : CALL get_qs_env(qs_env=qs_env, scf_control=scf_control)
399 : ! Provide full size work matrices
400 : CALL cp_fm_struct_create(fmstruct=fmstruct, &
401 : para_env=sub_env%para_env, &
402 : context=sub_env%blacs_env, &
403 : nrow_global=nsgf, &
404 0 : ncol_global=nsgf)
405 0 : CALL cp_fm_create(matrix=fm_s_half, matrix_struct=fmstruct, name="S^(1/2) MATRIX")
406 0 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, name="TMP MATRIX")
407 0 : CALL cp_fm_struct_release(fmstruct=fmstruct)
408 0 : CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
409 0 : CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
410 0 : IF (ndep /= 0) &
411 : CALL cp_warn(__LOCATION__, &
412 : "Overlap matrix exhibits linear dependencies. At least some "// &
413 0 : "eigenvalues have been quenched.")
414 0 : CALL copy_fm_to_dbcsr(fm_s_half, sm_h)
415 0 : CALL cp_fm_release(fm_s_half)
416 0 : CALL cp_fm_release(fm_work1)
417 0 : IF (iounit > 0) WRITE (iounit, *)
418 : END IF
419 :
420 408 : nspins = SIZE(sub_env%mos_occ)
421 :
422 848 : DO ispin = 1, nspins
423 440 : CALL cp_fm_get_info(work%ctransformed(ispin), ncol_global=nactive)
424 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, sub_env%mos_occ(ispin), &
425 848 : work%ctransformed(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
426 : END DO
427 :
428 : ! for Lowdin forces
429 408 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=work%S_eigenvectors%matrix_struct, name="TMP MATRIX")
430 408 : CALL copy_dbcsr_to_fm(sm_s, fm_work1)
431 408 : CALL choose_eigv_solver(fm_work1, work%S_eigenvectors, work%S_eigenvalues)
432 408 : CALL cp_fm_release(fm_work1)
433 : !
434 1224 : ALLOCATE (slam(nsgf, 1))
435 9546 : DO i = 1, nsgf
436 9546 : IF (work%S_eigenvalues(i) > 0._dp) THEN
437 9138 : slam(i, 1) = SQRT(work%S_eigenvalues(i))
438 : ELSE
439 0 : CPABORT("S matrix not positive definit")
440 : END IF
441 : END DO
442 9546 : DO i = 1, nsgf
443 9546 : CALL cp_fm_set_submatrix(work%slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
444 : END DO
445 9546 : DO i = 1, nsgf
446 9546 : CALL cp_fm_set_submatrix(work%slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .TRUE.)
447 : END DO
448 408 : CALL cp_fm_get_info(work%slambda, local_data=local_data)
449 9546 : DO i = 1, SIZE(local_data, 2)
450 420781 : DO j = 1, SIZE(local_data, 1)
451 411235 : sij = local_data(j, i)
452 411235 : IF (sij > 0.0_dp) sij = 1.0_dp/sij
453 420373 : local_data(j, i) = sij
454 : END DO
455 : END DO
456 408 : DEALLOCATE (slam)
457 :
458 408 : CALL timestop(handle)
459 :
460 1224 : END SUBROUTINE get_lowdin_mo_coefficients
461 :
462 : ! **************************************************************************************************
463 : !> \brief Calculate Lowdin transformed Davidson trial vector X
464 : !> shalf (dbcsr), xvec, xt (fm) are defined in the same sub_env
465 : !> \param shalf ...
466 : !> \param xvec ...
467 : !> \param xt ...
468 : ! **************************************************************************************************
469 6692 : SUBROUTINE get_lowdin_x(shalf, xvec, xt)
470 :
471 : TYPE(dbcsr_type), INTENT(IN) :: shalf
472 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: xvec
473 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: xt
474 :
475 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_x'
476 :
477 : INTEGER :: handle, ispin, nactive, nspins
478 :
479 6692 : CALL timeset(routineN, handle)
480 :
481 6692 : nspins = SIZE(xvec)
482 :
483 : ! Build Lowdin transformed tilde(X)= S^1/2 X for each spin
484 14080 : DO ispin = 1, nspins
485 7388 : CALL cp_fm_get_info(xt(ispin), ncol_global=nactive)
486 : CALL cp_dbcsr_sm_fm_multiply(shalf, xvec(ispin), &
487 14080 : xt(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
488 : END DO
489 :
490 6692 : CALL timestop(handle)
491 :
492 6692 : END SUBROUTINE get_lowdin_x
493 :
494 : ! **************************************************************************************************
495 : !> \brief ...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients --
496 : !> transition charges with the Coulomb-type or exchange-type integrals
497 : !> \param qs_env ...
498 : !> \param stda_control ...
499 : !> \param stda_env ...
500 : !> \param sub_env ...
501 : !> \param work ...
502 : !> \param is_rks_triplets ...
503 : !> \param X ...
504 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
505 : !> eigenvalue problem A*X = omega*X
506 : ! **************************************************************************************************
507 6462 : SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
508 6462 : work, is_rks_triplets, X, res)
509 :
510 : TYPE(qs_environment_type), POINTER :: qs_env
511 : TYPE(stda_control_type) :: stda_control
512 : TYPE(stda_env_type) :: stda_env
513 : TYPE(tddfpt_subgroup_env_type) :: sub_env
514 : TYPE(tddfpt_work_matrices) :: work
515 : LOGICAL, INTENT(IN) :: is_rks_triplets
516 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
517 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
518 :
519 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_calculate_kernel'
520 :
521 : INTEGER :: ewald_type, handle, ia, iatom, ikind, &
522 : is, ispin, jatom, jkind, jspin, natom, &
523 : nsgf, nspins
524 6462 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
525 : INTEGER, DIMENSION(2) :: nactive, nlim
526 : LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
527 : do_exchange, use_virial
528 : REAL(KIND=dp) :: alpha, bp, dr, eta, gabr, hfx, rbeta, &
529 : spinfac
530 6462 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
531 6462 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
532 : REAL(KIND=dp), DIMENSION(3) :: rij
533 6462 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gab, pblock
534 6462 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
535 : TYPE(cell_type), POINTER :: cell
536 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstructjspin
537 : TYPE(cp_fm_type) :: cvec, cvecjspin
538 6462 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
539 : TYPE(cp_fm_type), POINTER :: ct, ctjspin
540 : TYPE(dbcsr_iterator_type) :: iter
541 : TYPE(dbcsr_type) :: pdens
542 : TYPE(dbcsr_type), POINTER :: tempmat
543 : TYPE(ewald_environment_type), POINTER :: ewald_env
544 : TYPE(ewald_pw_type), POINTER :: ewald_pw
545 : TYPE(mp_para_env_type), POINTER :: para_env
546 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
547 6462 : POINTER :: n_list
548 6462 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
549 6462 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
550 : TYPE(virial_type), POINTER :: virial
551 :
552 6462 : CALL timeset(routineN, handle)
553 :
554 19386 : nactive(:) = stda_env%nactive(:)
555 6462 : nspins = SIZE(X)
556 6462 : spinfac = 2.0_dp
557 6462 : IF (nspins == 2) spinfac = 1.0_dp
558 :
559 5792 : IF (nspins == 1 .AND. is_rks_triplets) THEN
560 : do_coulomb = .FALSE.
561 : ELSE
562 : do_coulomb = .TRUE.
563 : END IF
564 6462 : do_ewald = stda_control%do_ewald
565 6462 : do_exchange = stda_control%do_exchange
566 :
567 6462 : para_env => sub_env%para_env
568 :
569 : CALL get_qs_env(qs_env, natom=natom, cell=cell, &
570 6462 : particle_set=particle_set, qs_kind_set=qs_kind_set)
571 19386 : ALLOCATE (first_sgf(natom))
572 12924 : ALLOCATE (last_sgf(natom))
573 6462 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
574 :
575 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
576 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
577 26518 : ALLOCATE (xtransformed(nspins))
578 13594 : DO ispin = 1, nspins
579 7132 : NULLIFY (fmstruct)
580 7132 : ct => work%ctransformed(ispin)
581 7132 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
582 13594 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
583 : END DO
584 6462 : CALL get_lowdin_x(work%shalf, X, xtransformed)
585 :
586 25848 : ALLOCATE (tcharge(natom), gtcharge(natom, 1))
587 :
588 13594 : DO ispin = 1, nspins
589 13594 : CALL cp_fm_set_all(res(ispin), 0.0_dp)
590 : END DO
591 :
592 13594 : DO ispin = 1, nspins
593 7132 : ct => work%ctransformed(ispin)
594 7132 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
595 21396 : ALLOCATE (tv(nsgf))
596 7132 : CALL cp_fm_create(cvec, fmstruct)
597 : !
598 : ! *** Coulomb contribution
599 : !
600 7132 : IF (do_coulomb) THEN
601 57016 : tcharge(:) = 0.0_dp
602 12488 : DO jspin = 1, nspins
603 6914 : ctjspin => work%ctransformed(jspin)
604 6914 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
605 6914 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
606 6914 : CALL cp_fm_create(cvecjspin, fmstructjspin)
607 : ! CV(mu,j) = CT(mu,j)*XT(mu,j)
608 6914 : CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
609 : ! TV(mu) = SUM_j CV(mu,j)
610 6914 : CALL cp_fm_vectorssum(cvecjspin, tv, "R")
611 : ! contract charges
612 : ! TC(a) = SUM_(mu of a) TV(mu)
613 62424 : DO ia = 1, natom
614 271198 : DO is = first_sgf(ia), last_sgf(ia)
615 264284 : tcharge(ia) = tcharge(ia) + tv(is)
616 : END DO
617 : END DO
618 19402 : CALL cp_fm_release(cvecjspin)
619 : END DO !jspin
620 : ! Apply tcharge*gab -> gtcharge
621 : ! gT(b) = SUM_a g(a,b)*TC(a)
622 : ! gab = work%gamma_exchange(1)%matrix
623 62590 : gtcharge = 0.0_dp
624 : ! short range contribution
625 5574 : tempmat => work%gamma_exchange(1)%matrix
626 5574 : CALL dbcsr_iterator_start(iter, tempmat)
627 878929 : DO WHILE (dbcsr_iterator_blocks_left(iter))
628 873355 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab)
629 873355 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
630 878929 : IF (iatom /= jatom) THEN
631 847634 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
632 : END IF
633 : END DO
634 5574 : CALL dbcsr_iterator_stop(iter)
635 : ! Ewald long range contribution
636 5574 : IF (do_ewald) THEN
637 740 : ewald_env => work%ewald_env
638 740 : ewald_pw => work%ewald_pw
639 740 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
640 740 : CALL get_qs_env(qs_env=qs_env, virial=virial)
641 740 : use_virial = .FALSE.
642 740 : calculate_forces = .FALSE.
643 740 : n_list => sub_env%sab_orb
644 740 : CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
645 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
646 740 : gtcharge, tcharge, calculate_forces, virial, use_virial)
647 : ! add self charge interaction contribution
648 740 : IF (para_env%is_source()) THEN
649 18724 : gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
650 : END IF
651 : ELSE
652 4834 : nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
653 12201 : DO iatom = nlim(1), nlim(2)
654 19800 : DO jatom = 1, iatom - 1
655 30396 : rij = particle_set(iatom)%r - particle_set(jatom)%r
656 30396 : rij = pbc(rij, cell)
657 30396 : dr = SQRT(SUM(rij(:)**2))
658 14966 : IF (dr > 1.e-6_dp) THEN
659 7599 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
660 7599 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
661 : END IF
662 : END DO
663 : END DO
664 : END IF
665 5574 : CALL para_env%sum(gtcharge)
666 : ! expand charges
667 : ! TV(mu) = TC(a of mu)
668 185688 : tv(1:nsgf) = 0.0_dp
669 57016 : DO ia = 1, natom
670 237130 : DO is = first_sgf(ia), last_sgf(ia)
671 231556 : tv(is) = gtcharge(ia, 1)
672 : END DO
673 : END DO
674 : ! CV(mu,i) = TV(mu)*CV(mu,i)
675 5574 : ct => work%ctransformed(ispin)
676 5574 : CALL cp_fm_to_fm(ct, cvec)
677 5574 : CALL cp_fm_row_scale(cvec, tv)
678 : ! rho(nu,i) = rho(nu,i) + Shalf(nu,mu)*CV(mu,i)
679 5574 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), spinfac, 1.0_dp)
680 : END IF
681 : !
682 : ! *** Exchange contribution
683 : !
684 7132 : IF (do_exchange) THEN ! option to explicitly switch off exchange
685 : ! (exchange contributes also if hfx_fraction=0)
686 6688 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
687 6688 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
688 : !
689 6688 : tempmat => work%shalf
690 6688 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
691 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
692 6688 : ct => work%ctransformed(ispin)
693 6688 : CALL dbcsr_set(pdens, 0.0_dp)
694 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
695 6688 : 1.0_dp, keep_sparsity=.FALSE.)
696 6688 : CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
697 : ! Apply PP*gab -> PP; gab = gamma_coulomb
698 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
699 6688 : bp = stda_env%beta_param
700 6688 : hfx = stda_env%hfx_fraction
701 6688 : CALL dbcsr_iterator_start(iter, pdens)
702 1738029 : DO WHILE (dbcsr_iterator_blocks_left(iter))
703 1731341 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
704 6925364 : rij = particle_set(iatom)%r - particle_set(jatom)%r
705 6925364 : rij = pbc(rij, cell)
706 6925364 : dr = SQRT(SUM(rij(:)**2))
707 1731341 : ikind = kind_of(iatom)
708 1731341 : jkind = kind_of(jatom)
709 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
710 1731341 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
711 1731341 : rbeta = dr**bp
712 1731341 : IF (hfx > 0.0_dp) THEN
713 1729303 : gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
714 : ELSE
715 2038 : IF (dr < 1.e-6) THEN
716 : gabr = 0.0_dp
717 : ELSE
718 1416 : gabr = 1._dp/dr
719 : END IF
720 : END IF
721 19415093 : pblock = gabr*pblock
722 : END DO
723 6688 : CALL dbcsr_iterator_stop(iter)
724 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
725 6688 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
726 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
727 6688 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), -1.0_dp, 1.0_dp)
728 : !
729 6688 : CALL dbcsr_release(pdens)
730 6688 : DEALLOCATE (kind_of)
731 : END IF
732 : !
733 7132 : CALL cp_fm_release(cvec)
734 27858 : DEALLOCATE (tv)
735 : END DO
736 :
737 6462 : CALL cp_fm_release(xtransformed)
738 6462 : DEALLOCATE (tcharge, gtcharge)
739 6462 : DEALLOCATE (first_sgf, last_sgf)
740 :
741 6462 : CALL timestop(handle)
742 :
743 12924 : END SUBROUTINE stda_calculate_kernel
744 :
745 : ! **************************************************************************************************
746 :
747 : END MODULE qs_tddfpt2_stda_utils
|
