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 440 : 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 440 : CALL timeset(routineN, handle)
121 :
122 440 : do_coulomb = .NOT. tddfpt_control%rks_triplets
123 440 : IF (do_coulomb) THEN
124 : ! calculate exchange gamma matrix
125 346 : 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 440 : CALL get_lowdin_mo_coefficients(qs_env, sub_env, work)
130 :
131 : ! initialize Ewald for sTDA
132 440 : IF (tddfpt_control%stda_control%do_ewald) THEN
133 106 : NULLIFY (ewald_env, ewald_pw)
134 1908 : ALLOCATE (ewald_env)
135 106 : CALL ewald_env_create(ewald_env, sub_env%para_env)
136 106 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
137 106 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
138 106 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
139 106 : print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
140 106 : 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 106 : cell_periodic=cell%perd)
143 106 : ALLOCATE (ewald_pw)
144 106 : CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
145 106 : work%ewald_env => ewald_env
146 106 : work%ewald_pw => ewald_pw
147 : END IF
148 :
149 440 : CALL timestop(handle)
150 :
151 440 : 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 536 : 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 536 : 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 536 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dgblock, gblock
180 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
181 : TYPE(neighbor_list_iterator_p_type), &
182 536 : DIMENSION(:), POINTER :: nl_iterator
183 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
184 536 : POINTER :: n_list
185 :
186 536 : CALL timeset(routineN, handle)
187 :
188 536 : CALL get_qs_env(qs_env=qs_env, natom=natom)
189 536 : 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 536 : n_list => sub_env%sab_orb
193 :
194 536 : IF (PRESENT(ndim)) THEN
195 190 : nmat = ndim
196 : ELSE
197 346 : nmat = 1
198 : END IF
199 536 : CPASSERT(nmat == 1 .OR. nmat == 4)
200 536 : CPASSERT(.NOT. ASSOCIATED(gamma_matrix))
201 536 : CALL dbcsr_allocate_matrix_set(gamma_matrix, nmat)
202 :
203 1608 : ALLOCATE (row_blk_sizes(natom))
204 3330 : row_blk_sizes(1:natom) = 1
205 1642 : DO imat = 1, nmat
206 1642 : 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 536 : col_blk_size=row_blk_sizes)
212 1106 : 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 1106 : col_blk_size=row_blk_sizes)
216 : END DO
217 :
218 536 : DEALLOCATE (row_blk_sizes)
219 :
220 : ! setup the matrices using the neighbor list
221 1642 : DO imat = 1, nmat
222 1106 : CALL cp_dbcsr_alloc_block_from_nbl(gamma_matrix(imat)%matrix, n_list)
223 1642 : CALL dbcsr_set(gamma_matrix(imat)%matrix, 0.0_dp)
224 : END DO
225 :
226 536 : NULLIFY (nl_iterator)
227 536 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
228 86606 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
230 86070 : iatom=iatom, jatom=jatom, r=rij)
231 :
232 344280 : dr = SQRT(SUM(rij(:)**2)) ! interatomic distance
233 :
234 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
235 86070 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
236 :
237 86070 : icol = MAX(iatom, jatom)
238 86070 : irow = MIN(iatom, jatom)
239 :
240 86070 : NULLIFY (gblock)
241 : CALL dbcsr_get_block_p(matrix=gamma_matrix(1)%matrix, &
242 86070 : row=irow, col=icol, BLOCK=gblock, found=found)
243 86070 : CPASSERT(found)
244 :
245 : ! get rcuta and rcutb
246 86070 : rcuta = stda_env%kind_param_set(ikind)%kind_param%rcut
247 86070 : rcutb = stda_env%kind_param_set(jkind)%kind_param%rcut
248 86070 : 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 86070 : IF (dr < 1.e-6) THEN
253 : ! on site terms
254 4191 : gblock(:, :) = gblock(:, :) + eta
255 84673 : ELSE IF (dr > rcut) THEN
256 : ! do nothing
257 : ELSE
258 40418 : IF (dr < rcut - rsmooth) THEN
259 : fcut = 1.0_dp
260 : ELSE
261 8799 : r = dr - (rcut - rsmooth)
262 8799 : x = r/rsmooth
263 8799 : 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 121254 : **(1._dp/stda_env%alpha_param) - fcut/dr
268 : END IF
269 :
270 86606 : 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 16574 : IF (dr < 1.e-6 .OR. dr > rcut) THEN
275 : ! on site terms or beyond cutoff
276 : dgb = 0.0_dp
277 : ELSE
278 8387 : IF (dr < rcut - rsmooth) THEN
279 : fcut = 1.0_dp
280 : dfcut = 0.0_dp
281 : ELSE
282 1840 : r = dr - (rcut - rsmooth)
283 1840 : x = r/rsmooth
284 1840 : fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
285 1840 : dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
286 1840 : dfcut = dfcut/rsmooth
287 : END IF
288 : dgb = dfcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
289 8387 : **(1._dp/stda_env%alpha_param)
290 8387 : dgb = dgb - dfcut/dr + fcut/dr**2
291 : dgb = dgb - fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
292 8387 : **(1._dp/stda_env%alpha_param + 1._dp)*dr**(stda_env%alpha_param - 1._dp)
293 : END IF
294 66296 : DO imat = 2, nmat
295 49722 : NULLIFY (dgblock)
296 : CALL dbcsr_get_block_p(matrix=gamma_matrix(imat)%matrix, &
297 49722 : row=irow, col=icol, BLOCK=dgblock, found=found)
298 66296 : IF (found) THEN
299 49722 : IF (dr > 1.e-6) THEN
300 48570 : i = imat - 1
301 48570 : IF (irow == iatom) THEN
302 77283 : dgblock(:, :) = dgblock(:, :) + dgb*rij(i)/dr
303 : ELSE
304 68427 : 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 536 : CALL neighbor_list_iterator_release(nl_iterator)
314 :
315 1642 : DO imat = 1, nmat
316 1642 : CALL dbcsr_finalize(gamma_matrix(imat)%matrix)
317 : END DO
318 :
319 536 : CALL timestop(handle)
320 :
321 536 : 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 446 : 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 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: slam
343 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
344 446 : 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 446 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
349 : TYPE(dbcsr_type) :: sm_hinv
350 : TYPE(dbcsr_type), POINTER :: sm_h, sm_s
351 446 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
352 : TYPE(scf_control_type), POINTER :: scf_control
353 :
354 446 : CALL timeset(routineN, handle)
355 :
356 446 : NULLIFY (logger) !get output_unit
357 446 : logger => cp_get_default_logger()
358 446 : iounit = cp_logger_get_default_io_unit(logger)
359 :
360 : ! Calculate S^1/2 matrix
361 446 : IF (iounit > 0) THEN
362 223 : WRITE (iounit, "(1X,A)") "", &
363 223 : "-------------------------------------------------------------------------------", &
364 223 : "- Create Matrix SQRT(S) -", &
365 446 : "-------------------------------------------------------------------------------"
366 : END IF
367 :
368 446 : IF (sub_env%is_split) THEN
369 0 : CPABORT('SPLIT')
370 : ELSE
371 446 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrixkp_s)
372 446 : CPASSERT(ASSOCIATED(matrixkp_s))
373 446 : CPWARN_IF(SIZE(matrixkp_s, 2) > 1, "not implemented for k-points.")
374 446 : sm_s => matrixkp_s(1, 1)%matrix
375 : END IF
376 446 : sm_h => work%shalf
377 :
378 446 : CALL dbcsr_create(sm_hinv, template=sm_s)
379 446 : CALL dbcsr_add_on_diag(sm_h, 1.0_dp)
380 446 : threshold = 1.0e-8_dp
381 446 : order_lanczos = 3
382 446 : eps_lanczos = 1.0e-4_dp
383 446 : 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 446 : converged=converged)
387 446 : CALL dbcsr_release(sm_hinv)
388 : !
389 446 : NULLIFY (qs_kind_set)
390 446 : 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 446 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
393 : !
394 446 : 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) THEN
412 : CALL cp_warn(__LOCATION__, &
413 : "Overlap matrix exhibits linear dependencies. At least some "// &
414 0 : "eigenvalues have been quenched.")
415 : END IF
416 0 : CALL copy_fm_to_dbcsr(fm_s_half, sm_h)
417 0 : CALL cp_fm_release(fm_s_half)
418 0 : CALL cp_fm_release(fm_work1)
419 0 : IF (iounit > 0) WRITE (iounit, *)
420 : END IF
421 :
422 446 : nspins = SIZE(sub_env%mos_occ)
423 :
424 954 : DO ispin = 1, nspins
425 508 : CALL cp_fm_get_info(work%ctransformed(ispin), ncol_global=nactive)
426 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, sub_env%mos_active(ispin), &
427 954 : work%ctransformed(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
428 : END DO
429 :
430 : ! for Lowdin forces
431 446 : CALL cp_fm_create(matrix=fm_work1, matrix_struct=work%S_eigenvectors%matrix_struct, name="TMP MATRIX")
432 446 : CALL copy_dbcsr_to_fm(sm_s, fm_work1)
433 446 : CALL choose_eigv_solver(fm_work1, work%S_eigenvectors, work%S_eigenvalues)
434 446 : CALL cp_fm_release(fm_work1)
435 : !
436 1338 : ALLOCATE (slam(nsgf, 1))
437 9986 : DO i = 1, nsgf
438 9986 : IF (work%S_eigenvalues(i) > 0._dp) THEN
439 9540 : slam(i, 1) = SQRT(work%S_eigenvalues(i))
440 : ELSE
441 0 : CPABORT("S matrix not positive definit")
442 : END IF
443 : END DO
444 9986 : DO i = 1, nsgf
445 9986 : CALL cp_fm_set_submatrix(work%slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
446 : END DO
447 9986 : DO i = 1, nsgf
448 9986 : CALL cp_fm_set_submatrix(work%slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .TRUE.)
449 : END DO
450 446 : CALL cp_fm_get_info(work%slambda, local_data=local_data)
451 9986 : DO i = 1, SIZE(local_data, 2)
452 423912 : DO j = 1, SIZE(local_data, 1)
453 413926 : sij = local_data(j, i)
454 413926 : IF (sij > 0.0_dp) sij = 1.0_dp/sij
455 423466 : local_data(j, i) = sij
456 : END DO
457 : END DO
458 446 : DEALLOCATE (slam)
459 :
460 446 : CALL timestop(handle)
461 :
462 1338 : END SUBROUTINE get_lowdin_mo_coefficients
463 :
464 : ! **************************************************************************************************
465 : !> \brief Calculate Lowdin transformed Davidson trial vector X
466 : !> shalf (dbcsr), xvec, xt (fm) are defined in the same sub_env
467 : !> \param shalf ...
468 : !> \param xvec ...
469 : !> \param xt ...
470 : ! **************************************************************************************************
471 7786 : SUBROUTINE get_lowdin_x(shalf, xvec, xt)
472 :
473 : TYPE(dbcsr_type), INTENT(IN) :: shalf
474 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: xvec
475 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: xt
476 :
477 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_x'
478 :
479 : INTEGER :: handle, ispin, nactive, nspins
480 :
481 7786 : CALL timeset(routineN, handle)
482 :
483 7786 : nspins = SIZE(xvec)
484 :
485 : ! Build Lowdin transformed tilde(X)= S^1/2 X for each spin
486 17318 : DO ispin = 1, nspins
487 9532 : CALL cp_fm_get_info(xt(ispin), ncol_global=nactive)
488 : CALL cp_dbcsr_sm_fm_multiply(shalf, xvec(ispin), &
489 17318 : xt(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
490 : END DO
491 :
492 7786 : CALL timestop(handle)
493 :
494 7786 : END SUBROUTINE get_lowdin_x
495 :
496 : ! **************************************************************************************************
497 : !> \brief ...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients --
498 : !> transition charges with the Coulomb-type or exchange-type integrals
499 : !> \param qs_env ...
500 : !> \param stda_control ...
501 : !> \param stda_env ...
502 : !> \param sub_env ...
503 : !> \param work ...
504 : !> \param is_rks_triplets ...
505 : !> \param X ...
506 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
507 : !> eigenvalue problem A*X = omega*X
508 : ! **************************************************************************************************
509 7544 : SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
510 7544 : work, is_rks_triplets, X, res)
511 :
512 : TYPE(qs_environment_type), POINTER :: qs_env
513 : TYPE(stda_control_type) :: stda_control
514 : TYPE(stda_env_type) :: stda_env
515 : TYPE(tddfpt_subgroup_env_type) :: sub_env
516 : TYPE(tddfpt_work_matrices) :: work
517 : LOGICAL, INTENT(IN) :: is_rks_triplets
518 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
519 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
520 :
521 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_calculate_kernel'
522 :
523 : INTEGER :: ewald_type, handle, ia, iatom, ikind, &
524 : is, ispin, jatom, jkind, jspin, natom, &
525 : nsgf, nspins
526 7544 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
527 : INTEGER, DIMENSION(2) :: nactive, nlim
528 : LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
529 : do_exchange, use_virial
530 : REAL(KIND=dp) :: alpha, bp, dr, eta, gabr, hfx, rbeta, &
531 : spinfac
532 7544 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
533 7544 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
534 : REAL(KIND=dp), DIMENSION(3) :: rij
535 7544 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gab, pblock
536 7544 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
537 : TYPE(cell_type), POINTER :: cell
538 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstructjspin
539 : TYPE(cp_fm_type) :: cvec, cvecjspin
540 7544 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
541 : TYPE(cp_fm_type), POINTER :: ct, ctjspin
542 : TYPE(dbcsr_iterator_type) :: iter
543 : TYPE(dbcsr_type) :: pdens
544 : TYPE(dbcsr_type), POINTER :: tempmat
545 : TYPE(ewald_environment_type), POINTER :: ewald_env
546 : TYPE(ewald_pw_type), POINTER :: ewald_pw
547 : TYPE(mp_para_env_type), POINTER :: para_env
548 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
549 7544 : POINTER :: n_list
550 7544 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
551 7544 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
552 : TYPE(virial_type), POINTER :: virial
553 :
554 7544 : CALL timeset(routineN, handle)
555 :
556 22632 : nactive(:) = stda_env%nactive(:)
557 7544 : nspins = SIZE(X)
558 7544 : spinfac = 2.0_dp
559 7544 : IF (nspins == 2) spinfac = 1.0_dp
560 :
561 5834 : IF (nspins == 1 .AND. is_rks_triplets) THEN
562 : do_coulomb = .FALSE.
563 : ELSE
564 : do_coulomb = .TRUE.
565 : END IF
566 7544 : do_ewald = stda_control%do_ewald
567 7544 : do_exchange = stda_control%do_exchange
568 :
569 7544 : para_env => sub_env%para_env
570 :
571 : CALL get_qs_env(qs_env, natom=natom, cell=cell, &
572 7544 : particle_set=particle_set, qs_kind_set=qs_kind_set)
573 22632 : ALLOCATE (first_sgf(natom))
574 15088 : ALLOCATE (last_sgf(natom))
575 7544 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
576 :
577 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
578 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
579 31886 : ALLOCATE (xtransformed(nspins))
580 16798 : DO ispin = 1, nspins
581 9254 : NULLIFY (fmstruct)
582 9254 : ct => work%ctransformed(ispin)
583 9254 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
584 16798 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
585 : END DO
586 7544 : CALL get_lowdin_x(work%shalf, X, xtransformed)
587 :
588 30176 : ALLOCATE (tcharge(natom), gtcharge(natom, 1))
589 :
590 16798 : DO ispin = 1, nspins
591 16798 : CALL cp_fm_set_all(res(ispin), 0.0_dp)
592 : END DO
593 :
594 16798 : DO ispin = 1, nspins
595 9254 : ct => work%ctransformed(ispin)
596 9254 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
597 27762 : ALLOCATE (tv(nsgf))
598 9254 : CALL cp_fm_create(cvec, fmstruct)
599 : !
600 : ! *** Coulomb contribution
601 : !
602 9254 : IF (do_coulomb) THEN
603 7716 : tcharge(:) = 0.0_dp
604 18852 : DO jspin = 1, nspins
605 11136 : ctjspin => work%ctransformed(jspin)
606 11136 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
607 11136 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
608 11136 : CALL cp_fm_create(cvecjspin, fmstructjspin)
609 : ! CV(mu,j) = CT(mu,j)*XT(mu,j)
610 11136 : CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
611 : ! TV(mu) = SUM_j CV(mu,j)
612 11136 : CALL cp_fm_vectorssum(cvecjspin, tv, "R")
613 : ! contract charges
614 : ! TC(a) = SUM_(mu of a) TV(mu)
615 79318 : DO ia = 1, natom
616 322432 : DO is = first_sgf(ia), last_sgf(ia)
617 311296 : tcharge(ia) = tcharge(ia) + tv(is)
618 : END DO
619 : END DO
620 29988 : CALL cp_fm_release(cvecjspin)
621 : END DO !jspin
622 : ! Apply tcharge*gab -> gtcharge
623 : ! gT(b) = SUM_a g(a,b)*TC(a)
624 : ! gab = work%gamma_exchange(1)%matrix
625 7716 : gtcharge = 0.0_dp
626 : ! short range contribution
627 7716 : tempmat => work%gamma_exchange(1)%matrix
628 7716 : CALL dbcsr_iterator_start(iter, tempmat)
629 887509 : DO WHILE (dbcsr_iterator_blocks_left(iter))
630 879793 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab)
631 879793 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
632 887509 : IF (iatom /= jatom) THEN
633 850856 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
634 : END IF
635 : END DO
636 7716 : CALL dbcsr_iterator_stop(iter)
637 : ! Ewald long range contribution
638 7716 : IF (do_ewald) THEN
639 1464 : ewald_env => work%ewald_env
640 1464 : ewald_pw => work%ewald_pw
641 1464 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
642 1464 : CALL get_qs_env(qs_env=qs_env, virial=virial)
643 1464 : use_virial = .FALSE.
644 1464 : calculate_forces = .FALSE.
645 1464 : n_list => sub_env%sab_orb
646 1464 : CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
647 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
648 1464 : gtcharge, tcharge, calculate_forces, virial, use_virial)
649 : ! add self charge interaction contribution
650 1464 : IF (para_env%is_source()) THEN
651 20172 : gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
652 : END IF
653 : ELSE
654 6252 : nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
655 15749 : DO iatom = nlim(1), nlim(2)
656 25484 : DO jatom = 1, iatom - 1
657 38940 : rij = particle_set(iatom)%r - particle_set(jatom)%r
658 38940 : rij = pbc(rij, cell)
659 38940 : dr = SQRT(SUM(rij(:)**2))
660 19232 : IF (dr > 1.e-6_dp) THEN
661 9735 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
662 9735 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
663 : END IF
664 : END DO
665 : END DO
666 : END IF
667 7716 : CALL para_env%sum(gtcharge)
668 : ! expand charges
669 : ! TV(mu) = TC(a of mu)
670 205530 : tv(1:nsgf) = 0.0_dp
671 65590 : DO ia = 1, natom
672 263404 : DO is = first_sgf(ia), last_sgf(ia)
673 255688 : tv(is) = gtcharge(ia, 1)
674 : END DO
675 : END DO
676 : ! CV(mu,i) = TV(mu)*CV(mu,i)
677 7716 : ct => work%ctransformed(ispin)
678 7716 : CALL cp_fm_to_fm(ct, cvec)
679 7716 : CALL cp_fm_row_scale(cvec, tv)
680 : ! rho(nu,i) = rho(nu,i) + Shalf(nu,mu)*CV(mu,i)
681 7716 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), spinfac, 1.0_dp)
682 : END IF
683 : !
684 : ! *** Exchange contribution
685 : !
686 9254 : IF (do_exchange) THEN ! option to explicitly switch off exchange
687 : ! (exchange contributes also if hfx_fraction=0)
688 8810 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
689 8810 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
690 : !
691 8810 : tempmat => work%shalf
692 8810 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
693 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
694 8810 : ct => work%ctransformed(ispin)
695 8810 : CALL dbcsr_set(pdens, 0.0_dp)
696 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
697 8810 : 1.0_dp, keep_sparsity=.FALSE.)
698 8810 : CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
699 : ! Apply PP*gab -> PP; gab = gamma_coulomb
700 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
701 8810 : bp = stda_env%beta_param
702 8810 : hfx = stda_env%hfx_fraction
703 8810 : CALL dbcsr_iterator_start(iter, pdens)
704 1748515 : DO WHILE (dbcsr_iterator_blocks_left(iter))
705 1739705 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
706 6958820 : rij = particle_set(iatom)%r - particle_set(jatom)%r
707 6958820 : rij = pbc(rij, cell)
708 6958820 : dr = SQRT(SUM(rij(:)**2))
709 1739705 : ikind = kind_of(iatom)
710 1739705 : jkind = kind_of(jatom)
711 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
712 1739705 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
713 1739705 : rbeta = dr**bp
714 1739705 : IF (hfx > 0.0_dp) THEN
715 1737583 : gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
716 : ELSE
717 2122 : IF (dr < 1.e-6) THEN
718 : gabr = 0.0_dp
719 : ELSE
720 1472 : gabr = 1._dp/dr
721 : END IF
722 : END IF
723 19517747 : pblock = gabr*pblock
724 : END DO
725 8810 : CALL dbcsr_iterator_stop(iter)
726 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
727 8810 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
728 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
729 8810 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), -1.0_dp, 1.0_dp)
730 : !
731 8810 : CALL dbcsr_release(pdens)
732 8810 : DEALLOCATE (kind_of)
733 : END IF
734 : !
735 9254 : CALL cp_fm_release(cvec)
736 35306 : DEALLOCATE (tv)
737 : END DO
738 :
739 7544 : CALL cp_fm_release(xtransformed)
740 7544 : DEALLOCATE (tcharge, gtcharge)
741 7544 : DEALLOCATE (first_sgf, last_sgf)
742 :
743 7544 : CALL timestop(handle)
744 :
745 15088 : END SUBROUTINE stda_calculate_kernel
746 :
747 : ! **************************************************************************************************
748 :
749 : END MODULE qs_tddfpt2_stda_utils
|
