(git:0de0cc2)
qs_tddfpt2_stda_utils.F
Go to the documentation of this file.
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 ! **************************************************************************************************
11 
12  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  USE cell_types, ONLY: cell_type,&
15  pbc
16  USE cp_control_types, ONLY: stda_control_type,&
17  tddfpt2_control_type
26  USE cp_fm_diag, ONLY: choose_eigv_solver,&
30  cp_fm_struct_type
31  USE cp_fm_types, ONLY: cp_fm_create,&
33  cp_fm_release,&
36  cp_fm_to_fm,&
37  cp_fm_type,&
41  cp_logger_type
42  USE dbcsr_api, ONLY: &
43  dbcsr_add_on_diag, dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
44  dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
45  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
46  dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
47  dbcsr_type_symmetric
51  ewald_environment_type,&
55  USE ewald_pw_types, ONLY: ewald_pw_create,&
56  ewald_pw_type
58  section_vals_type
60  USE kinds, ONLY: dp
61  USE mathconstants, ONLY: oorootpi
62  USE message_passing, ONLY: mp_para_env_type
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
72  neighbor_list_iterator_p_type,&
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 
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  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  CALL timeset(routinen, handle)
121 
122  do_coulomb = .NOT. tddfpt_control%rks_triplets
123  IF (do_coulomb) THEN
124  ! calculate exchange gamma matrix
125  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  CALL get_lowdin_mo_coefficients(qs_env, sub_env, work)
130 
131  ! initialize Ewald for sTDA
132  IF (tddfpt_control%stda_control%do_ewald) THEN
133  NULLIFY (ewald_env, ewald_pw)
134  ALLOCATE (ewald_env)
135  CALL ewald_env_create(ewald_env, sub_env%para_env)
136  poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
137  CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
138  ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
139  print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
140  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  ALLOCATE (ewald_pw)
143  CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
144  work%ewald_env => ewald_env
145  work%ewald_pw => ewald_pw
146  END IF
147 
148  CALL timestop(handle)
149 
150  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  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  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  REAL(kind=dp), DIMENSION(:, :), POINTER :: dgblock, gblock
179  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
180  TYPE(neighbor_list_iterator_p_type), &
181  DIMENSION(:), POINTER :: nl_iterator
182  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
183  POINTER :: n_list
184 
185  CALL timeset(routinen, handle)
186 
187  CALL get_qs_env(qs_env=qs_env, natom=natom)
188  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  n_list => sub_env%sab_orb
192 
193  IF (PRESENT(ndim)) THEN
194  nmat = ndim
195  ELSE
196  nmat = 1
197  END IF
198  cpassert(nmat == 1 .OR. nmat == 4)
199  cpassert(.NOT. ASSOCIATED(gamma_matrix))
200  CALL dbcsr_allocate_matrix_set(gamma_matrix, nmat)
201 
202  ALLOCATE (row_blk_sizes(natom))
203  row_blk_sizes(1:natom) = 1
204  DO imat = 1, nmat
205  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  col_blk_size=row_blk_sizes, nze=0)
211  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  col_blk_size=row_blk_sizes, nze=0)
215  END DO
216 
217  DEALLOCATE (row_blk_sizes)
218 
219  ! setup the matrices using the neighbor list
220  DO imat = 1, nmat
221  CALL cp_dbcsr_alloc_block_from_nbl(gamma_matrix(imat)%matrix, n_list)
222  CALL dbcsr_set(gamma_matrix(imat)%matrix, 0.0_dp)
223  END DO
224 
225  NULLIFY (nl_iterator)
226  CALL neighbor_list_iterator_create(nl_iterator, n_list)
227  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
228  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
229  iatom=iatom, jatom=jatom, r=rij)
230 
231  dr = sqrt(sum(rij(:)**2)) ! interatomic distance
232 
233  eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
234  stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
235 
236  icol = max(iatom, jatom)
237  irow = min(iatom, jatom)
238 
239  NULLIFY (gblock)
240  CALL dbcsr_get_block_p(matrix=gamma_matrix(1)%matrix, &
241  row=irow, col=icol, block=gblock, found=found)
242  cpassert(found)
243 
244  ! get rcuta and rcutb
245  rcuta = stda_env%kind_param_set(ikind)%kind_param%rcut
246  rcutb = stda_env%kind_param_set(jkind)%kind_param%rcut
247  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  IF (dr < 1.e-6) THEN
252  ! on site terms
253  gblock(:, :) = gblock(:, :) + eta
254  ELSEIF (dr > rcut) THEN
255  ! do nothing
256  ELSE
257  IF (dr < rcut - rsmooth) THEN
258  fcut = 1.0_dp
259  ELSE
260  r = dr - (rcut - rsmooth)
261  x = r/rsmooth
262  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  **(1._dp/stda_env%alpha_param) - fcut/dr
267  END IF
268 
269  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  IF (dr < 1.e-6 .OR. dr > rcut) THEN
274  ! on site terms or beyond cutoff
275  dgb = 0.0_dp
276  ELSE
277  IF (dr < rcut - rsmooth) THEN
278  fcut = 1.0_dp
279  dfcut = 0.0_dp
280  ELSE
281  r = dr - (rcut - rsmooth)
282  x = r/rsmooth
283  fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
284  dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
285  dfcut = dfcut/rsmooth
286  END IF
287  dgb = dfcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
288  **(1._dp/stda_env%alpha_param)
289  dgb = dgb - dfcut/dr + fcut/dr**2
290  dgb = dgb - fcut*(1._dp/(dr**(stda_env%alpha_param) + eta**(-stda_env%alpha_param))) &
291  **(1._dp/stda_env%alpha_param + 1._dp)*dr**(stda_env%alpha_param - 1._dp)
292  END IF
293  DO imat = 2, nmat
294  NULLIFY (dgblock)
295  CALL dbcsr_get_block_p(matrix=gamma_matrix(imat)%matrix, &
296  row=irow, col=icol, block=dgblock, found=found)
297  IF (found) THEN
298  IF (dr > 1.e-6) THEN
299  i = imat - 1
300  IF (irow == iatom) THEN
301  dgblock(:, :) = dgblock(:, :) + dgb*rij(i)/dr
302  ELSE
303  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  CALL neighbor_list_iterator_release(nl_iterator)
313 
314  DO imat = 1, nmat
315  CALL dbcsr_finalize(gamma_matrix(imat)%matrix)
316  END DO
317 
318  CALL timestop(handle)
319 
320  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: slam
342  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
343  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  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
348  TYPE(dbcsr_type) :: sm_hinv
349  TYPE(dbcsr_type), POINTER :: sm_h, sm_s
350  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
351  TYPE(scf_control_type), POINTER :: scf_control
352 
353  CALL timeset(routinen, handle)
354 
355  NULLIFY (logger) !get output_unit
356  logger => cp_get_default_logger()
357  iounit = cp_logger_get_default_io_unit(logger)
358 
359  ! Calculate S^1/2 matrix
360  IF (iounit > 0) THEN
361  WRITE (iounit, "(1X,A)") "", &
362  "-------------------------------------------------------------------------------", &
363  "- Create Matrix SQRT(S) -", &
364  "-------------------------------------------------------------------------------"
365  END IF
366 
367  IF (sub_env%is_split) THEN
368  cpabort('SPLIT')
369  ELSE
370  CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrixkp_s)
371  cpassert(ASSOCIATED(matrixkp_s))
372  IF (SIZE(matrixkp_s, 2) > 1) cpwarn("not implemented for k-points.")
373  sm_s => matrixkp_s(1, 1)%matrix
374  END IF
375  sm_h => work%shalf
376 
377  CALL dbcsr_create(sm_hinv, template=sm_s)
378  CALL dbcsr_add_on_diag(sm_h, 1.0_dp)
379  threshold = 1.0e-8_dp
380  order_lanczos = 3
381  eps_lanczos = 1.0e-4_dp
382  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  converged=converged)
386  CALL dbcsr_release(sm_hinv)
387  !
388  NULLIFY (qs_kind_set)
389  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  CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
392  !
393  IF (.NOT. converged) THEN
394  IF (iounit > 0) THEN
395  WRITE (iounit, "(T3,A)") "STDA| Newton-Schulz iteration did not converge"
396  WRITE (iounit, "(T3,A)") "STDA| Calculate SQRT(S) from diagonalization"
397  END IF
398  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  ncol_global=nsgf)
405  CALL cp_fm_create(matrix=fm_s_half, matrix_struct=fmstruct, name="S^(1/2) MATRIX")
406  CALL cp_fm_create(matrix=fm_work1, matrix_struct=fmstruct, name="TMP MATRIX")
407  CALL cp_fm_struct_release(fmstruct=fmstruct)
408  CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
409  CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
410  IF (ndep /= 0) &
411  CALL cp_warn(__location__, &
412  "Overlap matrix exhibits linear dependencies. At least some "// &
413  "eigenvalues have been quenched.")
414  CALL copy_fm_to_dbcsr(fm_s_half, sm_h)
415  CALL cp_fm_release(fm_s_half)
416  CALL cp_fm_release(fm_work1)
417  IF (iounit > 0) WRITE (iounit, *)
418  END IF
419 
420  nspins = SIZE(sub_env%mos_occ)
421 
422  DO ispin = 1, nspins
423  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  work%ctransformed(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
426  END DO
427 
428  ! for Lowdin forces
429  CALL cp_fm_create(matrix=fm_work1, matrix_struct=work%S_eigenvectors%matrix_struct, name="TMP MATRIX")
430  CALL copy_dbcsr_to_fm(sm_s, fm_work1)
431  CALL choose_eigv_solver(fm_work1, work%S_eigenvectors, work%S_eigenvalues)
432  CALL cp_fm_release(fm_work1)
433  !
434  ALLOCATE (slam(nsgf, 1))
435  DO i = 1, nsgf
436  IF (work%S_eigenvalues(i) > 0._dp) THEN
437  slam(i, 1) = sqrt(work%S_eigenvalues(i))
438  ELSE
439  cpabort("S matrix not positive definit")
440  END IF
441  END DO
442  DO i = 1, nsgf
443  CALL cp_fm_set_submatrix(work%slambda, slam, 1, i, nsgf, 1, 1.0_dp, 0.0_dp)
444  END DO
445  DO i = 1, nsgf
446  CALL cp_fm_set_submatrix(work%slambda, slam, i, 1, 1, nsgf, 1.0_dp, 1.0_dp, .true.)
447  END DO
448  CALL cp_fm_get_info(work%slambda, local_data=local_data)
449  DO i = 1, SIZE(local_data, 2)
450  DO j = 1, SIZE(local_data, 1)
451  sij = local_data(j, i)
452  IF (sij > 0.0_dp) sij = 1.0_dp/sij
453  local_data(j, i) = sij
454  END DO
455  END DO
456  DEALLOCATE (slam)
457 
458  CALL timestop(handle)
459 
460  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  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  CALL timeset(routinen, handle)
479 
480  nspins = SIZE(xvec)
481 
482  ! Build Lowdin transformed tilde(X)= S^1/2 X for each spin
483  DO ispin = 1, nspins
484  CALL cp_fm_get_info(xt(ispin), ncol_global=nactive)
485  CALL cp_dbcsr_sm_fm_multiply(shalf, xvec(ispin), &
486  xt(ispin), nactive, alpha=1.0_dp, beta=0.0_dp)
487  END DO
488 
489  CALL timestop(handle)
490 
491  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  SUBROUTINE stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, &
507  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
529  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
530  REAL(kind=dp), DIMENSION(3) :: rij
531  REAL(kind=dp), DIMENSION(:, :), POINTER :: gab, pblock
532  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  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  POINTER :: n_list
546  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
547  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
548  TYPE(virial_type), POINTER :: virial
549 
550  CALL timeset(routinen, handle)
551 
552  nactive(:) = stda_env%nactive(:)
553  nspins = SIZE(x)
554  spinfac = 2.0_dp
555  IF (nspins == 2) spinfac = 1.0_dp
556 
557  IF (nspins == 1 .AND. is_rks_triplets) THEN
558  do_coulomb = .false.
559  ELSE
560  do_coulomb = .true.
561  END IF
562  do_ewald = stda_control%do_ewald
563  do_exchange = stda_control%do_exchange
564 
565  para_env => sub_env%para_env
566 
567  CALL get_qs_env(qs_env, natom=natom, cell=cell, &
568  particle_set=particle_set, qs_kind_set=qs_kind_set)
569  ALLOCATE (first_sgf(natom))
570  ALLOCATE (last_sgf(natom))
571  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  ALLOCATE (xtransformed(nspins))
576  DO ispin = 1, nspins
577  NULLIFY (fmstruct)
578  ct => work%ctransformed(ispin)
579  CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
580  CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
581  END DO
582  CALL get_lowdin_x(work%shalf, x, xtransformed)
583 
584  ALLOCATE (tcharge(natom), gtcharge(natom, 1))
585 
586  DO ispin = 1, nspins
587  CALL cp_fm_set_all(res(ispin), 0.0_dp)
588  END DO
589 
590  DO ispin = 1, nspins
591  ct => work%ctransformed(ispin)
592  CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
593  ALLOCATE (tv(nsgf))
594  CALL cp_fm_create(cvec, fmstruct)
595  !
596  ! *** Coulomb contribution
597  !
598  IF (do_coulomb) THEN
599  tcharge(:) = 0.0_dp
600  DO jspin = 1, nspins
601  ctjspin => work%ctransformed(jspin)
602  CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
603  CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
604  CALL cp_fm_create(cvecjspin, fmstructjspin)
605  ! CV(mu,j) = CT(mu,j)*XT(mu,j)
606  CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
607  ! TV(mu) = SUM_j CV(mu,j)
608  CALL cp_fm_vectorssum(cvecjspin, tv, "R")
609  ! contract charges
610  ! TC(a) = SUM_(mu of a) TV(mu)
611  DO ia = 1, natom
612  DO is = first_sgf(ia), last_sgf(ia)
613  tcharge(ia) = tcharge(ia) + tv(is)
614  END DO
615  END DO
616  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  gtcharge = 0.0_dp
622  ! short range contribution
623  tempmat => work%gamma_exchange(1)%matrix
624  CALL dbcsr_iterator_start(iter, tempmat)
625  DO WHILE (dbcsr_iterator_blocks_left(iter))
626  CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab, blk)
627  gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
628  IF (iatom /= jatom) THEN
629  gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
630  END IF
631  END DO
632  CALL dbcsr_iterator_stop(iter)
633  ! Ewald long range contribution
634  IF (do_ewald) THEN
635  ewald_env => work%ewald_env
636  ewald_pw => work%ewald_pw
637  CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
638  CALL get_qs_env(qs_env=qs_env, virial=virial)
639  use_virial = .false.
640  calculate_forces = .false.
641  n_list => sub_env%sab_orb
642  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  gtcharge, tcharge, calculate_forces, virial, use_virial)
645  ! add self charge interaction contribution
646  IF (para_env%is_source()) THEN
647  gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
648  END IF
649  ELSE
650  nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
651  DO iatom = nlim(1), nlim(2)
652  DO jatom = 1, iatom - 1
653  rij = particle_set(iatom)%r - particle_set(jatom)%r
654  rij = pbc(rij, cell)
655  dr = sqrt(sum(rij(:)**2))
656  IF (dr > 1.e-6_dp) THEN
657  gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
658  gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
659  END IF
660  END DO
661  END DO
662  END IF
663  CALL para_env%sum(gtcharge)
664  ! expand charges
665  ! TV(mu) = TC(a of mu)
666  tv(1:nsgf) = 0.0_dp
667  DO ia = 1, natom
668  DO is = first_sgf(ia), last_sgf(ia)
669  tv(is) = gtcharge(ia, 1)
670  END DO
671  END DO
672  ! CV(mu,i) = TV(mu)*CV(mu,i)
673  ct => work%ctransformed(ispin)
674  CALL cp_fm_to_fm(ct, cvec)
675  CALL cp_fm_row_scale(cvec, tv)
676  ! rho(nu,i) = rho(nu,i) + Shalf(nu,mu)*CV(mu,i)
677  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  IF (do_exchange) THEN ! option to explicitly switch off exchange
683  ! (exchange contributes also if hfx_fraction=0)
684  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
685  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
686  !
687  tempmat => work%shalf
688  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  ct => work%ctransformed(ispin)
691  CALL dbcsr_set(pdens, 0.0_dp)
692  CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
693  1.0_dp, keep_sparsity=.false.)
694  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  bp = stda_env%beta_param
698  hfx = stda_env%hfx_fraction
699  CALL dbcsr_iterator_start(iter, pdens)
700  DO WHILE (dbcsr_iterator_blocks_left(iter))
701  CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock, blk)
702  rij = particle_set(iatom)%r - particle_set(jatom)%r
703  rij = pbc(rij, cell)
704  dr = sqrt(sum(rij(:)**2))
705  ikind = kind_of(iatom)
706  jkind = kind_of(jatom)
707  eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
708  stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
709  rbeta = dr**bp
710  IF (hfx > 0.0_dp) THEN
711  gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
712  ELSE
713  IF (dr < 1.e-6) THEN
714  gabr = 0.0_dp
715  ELSE
716  gabr = 1._dp/dr
717  END IF
718  END IF
719  pblock = gabr*pblock
720  END DO
721  CALL dbcsr_iterator_stop(iter)
722  ! CV(mu,i) = P(nu,mu)*CT(mu,i)
723  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  CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), -1.0_dp, 1.0_dp)
726  !
727  CALL dbcsr_release(pdens)
728  DEALLOCATE (kind_of)
729  END IF
730  !
731  CALL cp_fm_release(cvec)
732  DEALLOCATE (tv)
733  END DO
734 
735  CALL cp_fm_release(xtransformed)
736  DEALLOCATE (tcharge, gtcharge)
737  DEALLOCATE (first_sgf, last_sgf)
738 
739  CALL timestop(handle)
740 
741  END SUBROUTINE stda_calculate_kernel
742 
743 ! **************************************************************************************************
744 
745 END MODULE qs_tddfpt2_stda_utils
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_row_scale(matrixa, scaling)
scales row i of matrix a with scaling(i)
subroutine, public cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
computes the schur product of two matrices c_ij = a_ij * b_ij
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition: cp_fm_diag.F:896
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_vectorssum(matrix, sum_array, dir)
summing up all the elements along the matrix's i-th index or
Definition: cp_fm_types.F:1252
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition: cp_fm_types.F:768
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
subroutine, public ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, gmax, ns_max, precs, o_spline, para_env, poisson_section, interaction_cutoffs, cell_hmat)
Purpose: Set the EWALD environment.
subroutine, public ewald_env_create(ewald_env, para_env)
allocates and intitializes a ewald_env
subroutine, public read_ewald_section_tb(ewald_env, ewald_section, hmat)
Purpose: read the EWALD section for TB methods.
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Calculation of Ewald contributions in DFTB.
subroutine, public tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial, atprop)
...
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
subroutine, public ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
creates the structure ewald_pw_type
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Routines useful for iterative matrix calculations.
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public oorootpi
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Simplified Tamm Dancoff approach (sTDA).
Simplified Tamm Dancoff approach (sTDA).
subroutine, public get_lowdin_mo_coefficients(qs_env, sub_env, work)
Calculate Lowdin MO coefficients.
subroutine, public get_lowdin_x(shalf, xvec, xt)
Calculate Lowdin transformed Davidson trial vector X shalf (dbcsr), xvec, xt (fm) are defined in the ...
subroutine, public setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim)
Calculate sTDA exchange-type contributions.
subroutine, public stda_init_matrices(qs_env, stda_kernel, sub_env, work, tddfpt_control)
Calculate sTDA matrices.
subroutine, public stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work, is_rks_triplets, X, res)
...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients – transition char...
parameters that control an scf iteration
All kind of helpful little routines.
Definition: util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333