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 : !> \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_add_on_diag, dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
20 : dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
21 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
22 : dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
23 : dbcsr_type_symmetric
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, nze=0)
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, nze=0)
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 6712 : SUBROUTINE get_lowdin_x(shalf, xvec, xt)
470 :
471 : TYPE(dbcsr_type) :: shalf
472 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: xvec, xt
473 :
474 : CHARACTER(len=*), PARAMETER :: routineN = 'get_lowdin_x'
475 :
476 : INTEGER :: handle, ispin, nactive, nspins
477 :
478 6712 : CALL timeset(routineN, handle)
479 :
480 6712 : nspins = SIZE(xvec)
481 :
482 : ! Build Lowdin transformed tilde(X)= S^1/2 X for each spin
483 14120 : DO ispin = 1, nspins
484 7408 : CALL cp_fm_get_info(xt(ispin), ncol_global=nactive)
485 : CALL cp_dbcsr_sm_fm_multiply(shalf, xvec(ispin), &
486 14120 : xt(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
487 : END DO
488 :
489 6712 : CALL timestop(handle)
490 :
491 6712 : END SUBROUTINE get_lowdin_x
492 :
493 : ! **************************************************************************************************
494 : !> \brief ...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients --
495 : !> transition charges with the Coulomb-type or exchange-type integrals
496 : !> \param qs_env ...
497 : !> \param stda_control ...
498 : !> \param stda_env ...
499 : !> \param sub_env ...
500 : !> \param work ...
501 : !> \param is_rks_triplets ...
502 : !> \param X ...
503 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
504 : !> eigenvalue problem A*X = omega*X
505 : ! **************************************************************************************************
506 6482 : SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
507 6482 : work, is_rks_triplets, X, res)
508 :
509 : TYPE(qs_environment_type), POINTER :: qs_env
510 : TYPE(stda_control_type) :: stda_control
511 : TYPE(stda_env_type) :: stda_env
512 : TYPE(tddfpt_subgroup_env_type) :: sub_env
513 : TYPE(tddfpt_work_matrices) :: work
514 : LOGICAL, INTENT(IN) :: is_rks_triplets
515 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X, res
516 :
517 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_calculate_kernel'
518 :
519 : INTEGER :: blk, ewald_type, handle, ia, iatom, &
520 : ikind, is, ispin, jatom, jkind, jspin, &
521 : natom, nsgf, nspins
522 6482 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
523 : INTEGER, DIMENSION(2) :: nactive, nlim
524 : LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
525 : do_exchange, use_virial
526 : REAL(KIND=dp) :: alpha, bp, dr, eta, gabr, hfx, rbeta, &
527 : spinfac
528 6482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
529 6482 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
530 : REAL(KIND=dp), DIMENSION(3) :: rij
531 6482 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gab, pblock
532 6482 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
533 : TYPE(cell_type), POINTER :: cell
534 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstructjspin
535 : TYPE(cp_fm_type) :: cvec, cvecjspin
536 6482 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
537 : TYPE(cp_fm_type), POINTER :: ct, ctjspin
538 : TYPE(dbcsr_iterator_type) :: iter
539 : TYPE(dbcsr_type) :: pdens
540 : TYPE(dbcsr_type), POINTER :: tempmat
541 : TYPE(ewald_environment_type), POINTER :: ewald_env
542 : TYPE(ewald_pw_type), POINTER :: ewald_pw
543 : TYPE(mp_para_env_type), POINTER :: para_env
544 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
545 6482 : POINTER :: n_list
546 6482 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
547 6482 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
548 : TYPE(virial_type), POINTER :: virial
549 :
550 6482 : CALL timeset(routineN, handle)
551 :
552 19446 : nactive(:) = stda_env%nactive(:)
553 6482 : nspins = SIZE(X)
554 6482 : spinfac = 2.0_dp
555 6482 : IF (nspins == 2) spinfac = 1.0_dp
556 :
557 5812 : IF (nspins == 1 .AND. is_rks_triplets) THEN
558 : do_coulomb = .FALSE.
559 : ELSE
560 : do_coulomb = .TRUE.
561 : END IF
562 6482 : do_ewald = stda_control%do_ewald
563 6482 : do_exchange = stda_control%do_exchange
564 :
565 6482 : para_env => sub_env%para_env
566 :
567 : CALL get_qs_env(qs_env, natom=natom, cell=cell, &
568 6482 : particle_set=particle_set, qs_kind_set=qs_kind_set)
569 19446 : ALLOCATE (first_sgf(natom))
570 12964 : ALLOCATE (last_sgf(natom))
571 6482 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
572 :
573 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
574 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
575 26598 : ALLOCATE (xtransformed(nspins))
576 13634 : DO ispin = 1, nspins
577 7152 : NULLIFY (fmstruct)
578 7152 : ct => work%ctransformed(ispin)
579 7152 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
580 13634 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
581 : END DO
582 6482 : CALL get_lowdin_x(work%shalf, X, xtransformed)
583 :
584 25928 : ALLOCATE (tcharge(natom), gtcharge(natom, 1))
585 :
586 13634 : DO ispin = 1, nspins
587 13634 : CALL cp_fm_set_all(res(ispin), 0.0_dp)
588 : END DO
589 :
590 13634 : DO ispin = 1, nspins
591 7152 : ct => work%ctransformed(ispin)
592 7152 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
593 21456 : ALLOCATE (tv(nsgf))
594 7152 : CALL cp_fm_create(cvec, fmstruct)
595 : !
596 : ! *** Coulomb contribution
597 : !
598 7152 : IF (do_coulomb) THEN
599 57096 : tcharge(:) = 0.0_dp
600 12528 : DO jspin = 1, nspins
601 6934 : ctjspin => work%ctransformed(jspin)
602 6934 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
603 6934 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
604 6934 : CALL cp_fm_create(cvecjspin, fmstructjspin)
605 : ! CV(mu,j) = CT(mu,j)*XT(mu,j)
606 6934 : CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
607 : ! TV(mu) = SUM_j CV(mu,j)
608 6934 : CALL cp_fm_vectorssum(cvecjspin, tv, "R")
609 : ! contract charges
610 : ! TC(a) = SUM_(mu of a) TV(mu)
611 62504 : DO ia = 1, natom
612 271438 : DO is = first_sgf(ia), last_sgf(ia)
613 264504 : tcharge(ia) = tcharge(ia) + tv(is)
614 : END DO
615 : END DO
616 19462 : CALL cp_fm_release(cvecjspin)
617 : END DO !jspin
618 : ! Apply tcharge*gab -> gtcharge
619 : ! gT(b) = SUM_a g(a,b)*TC(a)
620 : ! gab = work%gamma_exchange(1)%matrix
621 62690 : gtcharge = 0.0_dp
622 : ! short range contribution
623 5594 : tempmat => work%gamma_exchange(1)%matrix
624 5594 : CALL dbcsr_iterator_start(iter, tempmat)
625 879009 : DO WHILE (dbcsr_iterator_blocks_left(iter))
626 873415 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab, blk)
627 873415 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
628 879009 : IF (iatom /= jatom) THEN
629 847664 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
630 : END IF
631 : END DO
632 5594 : CALL dbcsr_iterator_stop(iter)
633 : ! Ewald long range contribution
634 5594 : IF (do_ewald) THEN
635 740 : ewald_env => work%ewald_env
636 740 : ewald_pw => work%ewald_pw
637 740 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
638 740 : CALL get_qs_env(qs_env=qs_env, virial=virial)
639 740 : use_virial = .FALSE.
640 740 : calculate_forces = .FALSE.
641 740 : n_list => sub_env%sab_orb
642 740 : CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
643 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
644 740 : gtcharge, tcharge, calculate_forces, virial, use_virial)
645 : ! add self charge interaction contribution
646 740 : IF (para_env%is_source()) THEN
647 18724 : gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
648 : END IF
649 : ELSE
650 4854 : nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
651 12251 : DO iatom = nlim(1), nlim(2)
652 19880 : DO jatom = 1, iatom - 1
653 30516 : rij = particle_set(iatom)%r - particle_set(jatom)%r
654 30516 : rij = pbc(rij, cell)
655 30516 : dr = SQRT(SUM(rij(:)**2))
656 15026 : IF (dr > 1.e-6_dp) THEN
657 7629 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
658 7629 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
659 : END IF
660 : END DO
661 : END DO
662 : END IF
663 5594 : CALL para_env%sum(gtcharge)
664 : ! expand charges
665 : ! TV(mu) = TC(a of mu)
666 185868 : tv(1:nsgf) = 0.0_dp
667 57096 : DO ia = 1, natom
668 237370 : DO is = first_sgf(ia), last_sgf(ia)
669 231776 : tv(is) = gtcharge(ia, 1)
670 : END DO
671 : END DO
672 : ! CV(mu,i) = TV(mu)*CV(mu,i)
673 5594 : ct => work%ctransformed(ispin)
674 5594 : CALL cp_fm_to_fm(ct, cvec)
675 5594 : CALL cp_fm_row_scale(cvec, tv)
676 : ! rho(nu,i) = rho(nu,i) + Shalf(nu,mu)*CV(mu,i)
677 5594 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), spinfac, 1.0_dp)
678 : END IF
679 : !
680 : ! *** Exchange contribution
681 : !
682 7152 : IF (do_exchange) THEN ! option to explicitly switch off exchange
683 : ! (exchange contributes also if hfx_fraction=0)
684 6708 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
685 6708 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
686 : !
687 6708 : tempmat => work%shalf
688 6708 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
689 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
690 6708 : ct => work%ctransformed(ispin)
691 6708 : CALL dbcsr_set(pdens, 0.0_dp)
692 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
693 6708 : 1.0_dp, keep_sparsity=.FALSE.)
694 6708 : CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
695 : ! Apply PP*gab -> PP; gab = gamma_coulomb
696 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
697 6708 : bp = stda_env%beta_param
698 6708 : hfx = stda_env%hfx_fraction
699 6708 : CALL dbcsr_iterator_start(iter, pdens)
700 1738125 : DO WHILE (dbcsr_iterator_blocks_left(iter))
701 1731417 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock, blk)
702 6925668 : rij = particle_set(iatom)%r - particle_set(jatom)%r
703 6925668 : rij = pbc(rij, cell)
704 6925668 : dr = SQRT(SUM(rij(:)**2))
705 1731417 : ikind = kind_of(iatom)
706 1731417 : jkind = kind_of(jatom)
707 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
708 1731417 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
709 1731417 : rbeta = dr**bp
710 1731417 : IF (hfx > 0.0_dp) THEN
711 1729379 : gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
712 : ELSE
713 2038 : IF (dr < 1.e-6) THEN
714 : gabr = 0.0_dp
715 : ELSE
716 1416 : gabr = 1._dp/dr
717 : END IF
718 : END IF
719 19415969 : pblock = gabr*pblock
720 : END DO
721 6708 : CALL dbcsr_iterator_stop(iter)
722 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
723 6708 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
724 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
725 6708 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), -1.0_dp, 1.0_dp)
726 : !
727 6708 : CALL dbcsr_release(pdens)
728 6708 : DEALLOCATE (kind_of)
729 : END IF
730 : !
731 7152 : CALL cp_fm_release(cvec)
732 27938 : DEALLOCATE (tv)
733 : END DO
734 :
735 6482 : CALL cp_fm_release(xtransformed)
736 6482 : DEALLOCATE (tcharge, gtcharge)
737 6482 : DEALLOCATE (first_sgf, last_sgf)
738 :
739 6482 : CALL timestop(handle)
740 :
741 12964 : END SUBROUTINE stda_calculate_kernel
742 :
743 : ! **************************************************************************************************
744 :
745 : END MODULE qs_tddfpt2_stda_utils
|
