(git:e7e05ae)
ri_environment_methods.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 ! **************************************************************************************************
9 !> \brief Calculates integral matrices for RIGPW method
10 !> \par History
11 !> created JGH [08.2012]
12 !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
13 !> (2) heavily debugged
14 !> updated for RI JGH [08.2017]
15 !> \authors JGH
16 !> Dorothea Golze
17 ! **************************************************************************************************
18 MODULE ri_environment_methods
19  USE arnoldi_api, ONLY: arnoldi_conjugate_gradient
20  USE atomic_kind_types, ONLY: atomic_kind_type,&
21  get_atomic_kind_set
22  USE cp_blacs_env, ONLY: cp_blacs_env_type
23  USE cp_control_types, ONLY: dft_control_type
24  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25  USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
26  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
27  USE dbcsr_api, ONLY: &
28  dbcsr_add, dbcsr_create, dbcsr_dot, dbcsr_get_info, dbcsr_iterator_blocks_left, &
29  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
30  dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale_by_vector, dbcsr_set, dbcsr_type, &
31  dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
32  USE iterate_matrix, ONLY: invert_hotelling
33  USE kinds, ONLY: dp
34  USE lri_environment_types, ONLY: allocate_lri_coefs,&
35  lri_density_create,&
36  lri_density_release,&
37  lri_density_type,&
38  lri_environment_type,&
39  lri_kind_type
40  USE message_passing, ONLY: mp_comm_type,&
41  mp_para_env_type
42  USE pw_types, ONLY: pw_c1d_gs_type,&
43  pw_r3d_rs_type
44  USE qs_collocate_density, ONLY: calculate_lri_rho_elec
45  USE qs_environment_types, ONLY: get_qs_env,&
46  qs_environment_type,&
47  set_qs_env
48  USE qs_ks_types, ONLY: qs_ks_env_type
49  USE qs_o3c_methods, ONLY: calculate_o3c_integrals,&
50  contract12_o3c
51  USE qs_o3c_types, ONLY: get_o3c_iterator_info,&
52  init_o3c_container,&
53  o3c_container_type,&
54  o3c_iterate,&
55  o3c_iterator_create,&
56  o3c_iterator_release,&
57  o3c_iterator_type,&
58  release_o3c_container
59  USE qs_overlap, ONLY: build_overlap_matrix
60  USE qs_rho_types, ONLY: qs_rho_get,&
61  qs_rho_type
62  USE util, ONLY: get_limit
63 
64 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
65 #include "./base/base_uses.f90"
66 
67  IMPLICIT NONE
68 
69  PRIVATE
70 
71 ! **************************************************************************************************
72 
73  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_environment_methods'
74 
75  PUBLIC :: build_ri_matrices, calculate_ri_densities, ri_metric_solver
76 
77 ! **************************************************************************************************
78 
79 CONTAINS
80 
81 ! **************************************************************************************************
82 !> \brief creates and initializes an lri_env
83 !> \param lri_env the lri_environment you want to create
84 !> \param qs_env ...
85 !> \param calculate_forces ...
86 ! **************************************************************************************************
87  SUBROUTINE build_ri_matrices(lri_env, qs_env, calculate_forces)
88 
89  TYPE(lri_environment_type), POINTER :: lri_env
90  TYPE(qs_environment_type), POINTER :: qs_env
91  LOGICAL, INTENT(IN) :: calculate_forces
92 
93  CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces)
94 
95  END SUBROUTINE build_ri_matrices
96 
97 ! **************************************************************************************************
98 !> \brief calculates integrals needed for the RI density fitting,
99 !> integrals are calculated once, before the SCF starts
100 !> \param lri_env ...
101 !> \param qs_env ...
102 !> \param calculate_forces ...
103 ! **************************************************************************************************
104  SUBROUTINE calculate_ri_integrals(lri_env, qs_env, calculate_forces)
105 
106  TYPE(lri_environment_type), POINTER :: lri_env
107  TYPE(qs_environment_type), POINTER :: qs_env
108  LOGICAL, INTENT(IN) :: calculate_forces
109 
110  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_ri_integrals'
111  REAL(kind=dp), DIMENSION(2), PARAMETER :: fx = (/0.0_dp, 1.0_dp/)
112 
113  INTEGER :: handle, i, i1, i2, iatom, ispin, izero, &
114  j, j1, j2, jatom, m, n, nbas
115  INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
116  REAL(kind=dp) :: fpre
117  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eval
118  REAL(kind=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout
119  TYPE(cp_blacs_env_type), POINTER :: blacs_env
120  TYPE(dbcsr_iterator_type) :: iter
121  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
122  TYPE(dbcsr_type) :: emat
123  TYPE(dbcsr_type), POINTER :: fmat
124  TYPE(dft_control_type), POINTER :: dft_control
125  TYPE(mp_para_env_type), POINTER :: para_env
126  TYPE(qs_ks_env_type), POINTER :: ks_env
127  TYPE(qs_rho_type), POINTER :: rho
128 
129  CALL timeset(routinen, handle)
130 
131  cpassert(ASSOCIATED(lri_env))
132  cpassert(ASSOCIATED(qs_env))
133 
134  ! overlap matrices
135  CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
136  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
137  NULLIFY (rho, matrix_p)
138  ! orbital overlap matrix (needed for N constraints and forces)
139  ! we replicate this in order to not directly interact qith qs_env
140  IF (calculate_forces) THEN
141  ! charge constraint forces are calculated here
142  CALL get_qs_env(qs_env=qs_env, rho=rho)
143  CALL qs_rho_get(rho, rho_ao=matrix_p)
144  ALLOCATE (fmat)
145  CALL dbcsr_create(fmat, template=matrix_p(1)%matrix)
146  DO ispin = 1, dft_control%nspins
147  fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n
148  CALL dbcsr_add(fmat, matrix_p(ispin)%matrix, fx(ispin), -fpre)
149  END DO
150  CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
151  basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list, &
152  calculate_forces=.true., matrix_p=fmat)
153  CALL dbcsr_release(fmat)
154  DEALLOCATE (fmat)
155  ELSE
156  CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
157  basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list)
158  END IF
159  ! RI overlap matrix
160  CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, &
161  basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list)
162  IF (calculate_forces) THEN
163  ! calculate f*a(T) pseudo density matrix for forces
164  bas_ptr => lri_env%ri_fit%bas_ptr
165  avec => lri_env%ri_fit%avec
166  fout => lri_env%ri_fit%fout
167  ALLOCATE (fmat)
168  CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix)
169  CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list)
170  CALL dbcsr_set(fmat, 0.0_dp)
171  CALL dbcsr_iterator_start(iter, fmat)
172  DO WHILE (dbcsr_iterator_blocks_left(iter))
173  CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk)
174  i1 = bas_ptr(1, iatom)
175  i2 = bas_ptr(2, iatom)
176  j1 = bas_ptr(1, jatom)
177  j2 = bas_ptr(2, jatom)
178  IF (iatom <= jatom) THEN
179  DO ispin = 1, dft_control%nspins
180  DO j = j1, j2
181  m = j - j1 + 1
182  DO i = i1, i2
183  n = i - i1 + 1
184  fblk(n, m) = fblk(n, m) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin)
185  END DO
186  END DO
187  END DO
188  ELSE
189  DO ispin = 1, dft_control%nspins
190  DO i = i1, i2
191  n = i - i1 + 1
192  DO j = j1, j2
193  m = j - j1 + 1
194  fblk(m, n) = fblk(m, n) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin)
195  END DO
196  END DO
197  END DO
198  END IF
199  IF (iatom == jatom) THEN
200  fblk(:, :) = 0.25_dp*fblk(:, :)
201  ELSE
202  fblk(:, :) = 0.5_dp*fblk(:, :)
203  END IF
204  END DO
205  CALL dbcsr_iterator_stop(iter)
206  !
207  CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, &
208  basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list, &
209  calculate_forces=.true., matrix_p=fmat)
210  CALL dbcsr_release(fmat)
211  DEALLOCATE (fmat)
212  END IF
213  ! approximation (preconditioner) or exact inverse of RI overlap
214  CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1)
215  ALLOCATE (lri_env%ri_sinv(1)%matrix)
216  SELECT CASE (lri_env%ri_sinv_app)
217  CASE ("NONE")
218  !
219  CASE ("INVS")
220  CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
221  CALL invert_hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, &
222  threshold=1.e-10_dp, use_inv_as_guess=.false., &
223  norm_convergence=1.e-10_dp, filter_eps=1.e-12_dp, silent=.false.)
224  CASE ("INVF")
225  CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix)
226  CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
227  CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas)
228  ALLOCATE (eval(nbas))
229  CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env)
230  izero = 0
231  DO i = 1, nbas
232  IF (eval(i) < 1.0e-10_dp) THEN
233  eval(i) = 0.0_dp
234  izero = izero + 1
235  ELSE
236  eval(i) = sqrt(1.0_dp/eval(i))
237  END IF
238  END DO
239  CALL dbcsr_scale_by_vector(emat, eval, side='right')
240  CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
241  CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix)
242  DEALLOCATE (eval)
243  CALL dbcsr_release(emat)
244  CASE ("AINV")
245  CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
246  CALL invert_hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, &
247  threshold=1.e-5_dp, use_inv_as_guess=.false., &
248  norm_convergence=1.e-4_dp, filter_eps=1.e-4_dp, silent=.false.)
249  CASE DEFAULT
250  cpabort("Unknown RI_SINV type")
251  END SELECT
252 
253  ! solve Rx=n
254  CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
255  vecr=lri_env%ri_fit%nvec, &
256  vecx=lri_env%ri_fit%rm1n, &
257  matp=lri_env%ri_sinv(1)%matrix, &
258  solver=lri_env%ri_sinv_app, &
259  ptr=lri_env%ri_fit%bas_ptr)
260 
261  ! calculate n(t)x
262  lri_env%ri_fit%ntrm1n = sum(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:))
263 
264  ! calculate 3c-overlap integrals
265  IF (ASSOCIATED(lri_env%o3c)) THEN
266  CALL release_o3c_container(lri_env%o3c)
267  ELSE
268  ALLOCATE (lri_env%o3c)
269  END IF
270  CALL init_o3c_container(lri_env%o3c, dft_control%nspins, &
271  lri_env%orb_basis, lri_env%orb_basis, lri_env%ri_basis, &
272  lri_env%soo_list, lri_env%soa_list)
273 
274  CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p)
275 
276  CALL timestop(handle)
277 
278  END SUBROUTINE calculate_ri_integrals
279 ! **************************************************************************************************
280 !> \brief solver for RI systems (R*x=n)
281 !> \param mat ...
282 !> \param vecr ...
283 !> \param vecx ...
284 !> \param matp ...
285 !> \param solver ...
286 !> \param ptr ...
287 ! **************************************************************************************************
288  SUBROUTINE ri_metric_solver(mat, vecr, vecx, matp, solver, ptr)
289 
290  TYPE(dbcsr_type) :: mat
291  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vecr
292  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: vecx
293  TYPE(dbcsr_type) :: matp
294  CHARACTER(LEN=*), INTENT(IN) :: solver
295  INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr
296 
297  CHARACTER(LEN=*), PARAMETER :: routinen = 'ri_metric_solver'
298 
299  INTEGER :: handle, max_iter, n
300  LOGICAL :: converged
301  REAL(kind=dp) :: rerror, threshold
302  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vect
303 
304  CALL timeset(routinen, handle)
305 
306  threshold = 1.e-10_dp
307  max_iter = 100
308 
309  vecx(:) = vecr(:)
310 
311  SELECT CASE (solver)
312  CASE ("NONE")
313  CALL arnoldi_conjugate_gradient(mat, vecx, &
314  converged=converged, threshold=threshold, max_iter=max_iter)
315  CASE ("INVS", "INVF")
316  converged = .false.
317  CALL ri_matvec(matp, vecr, vecx, ptr)
318  CASE ("AINV")
319  CALL arnoldi_conjugate_gradient(mat, vecx, matrix_p=matp, &
320  converged=converged, threshold=threshold, max_iter=max_iter)
321  CASE DEFAULT
322  cpabort("Unknown RI solver")
323  END SELECT
324 
325  IF (.NOT. converged) THEN
326  ! get error
327  rerror = 0.0_dp
328  n = SIZE(vecr)
329  ALLOCATE (vect(n))
330  CALL ri_matvec(mat, vecx, vect, ptr)
331  vect(:) = vect(:) - vecr(:)
332  rerror = maxval(abs(vect(:)))
333  DEALLOCATE (vect)
334  IF (rerror > threshold) THEN
335  cpwarn("RI solver: CG did not converge properly")
336  END IF
337  END IF
338 
339  CALL timestop(handle)
340 
341  END SUBROUTINE ri_metric_solver
342 
343 ! **************************************************************************************************
344 !> \brief performs the fitting of the density and distributes the fitted
345 !> density on the grid
346 !> \param lri_env the lri environment
347 !> lri_density the environment for the fitting
348 !> pmatrix density matrix
349 !> lri_rho_struct where the fitted density is stored
350 !> \param qs_env ...
351 !> \param pmatrix ...
352 !> \param lri_rho_struct ...
353 !> \param atomic_kind_set ...
354 !> \param para_env ...
355 ! **************************************************************************************************
356  SUBROUTINE calculate_ri_densities(lri_env, qs_env, pmatrix, &
357  lri_rho_struct, atomic_kind_set, para_env)
358 
359  TYPE(lri_environment_type), POINTER :: lri_env
360  TYPE(qs_environment_type), POINTER :: qs_env
361  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
362  TYPE(qs_rho_type), INTENT(IN) :: lri_rho_struct
363  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
364  TYPE(mp_para_env_type), POINTER :: para_env
365 
366  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_ri_densities'
367 
368  INTEGER :: atom_a, handle, i1, i2, iatom, ikind, &
369  ispin, n, natom, nspin
370  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
371  INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
372  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
373  REAL(kind=dp), DIMENSION(:, :), POINTER :: avec
374  TYPE(lri_density_type), POINTER :: lri_density
375  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
376  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
377  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
378 
379  CALL timeset(routinen, handle)
380 
381  nspin = SIZE(pmatrix, 1)
382  CALL contract12_o3c(lri_env%o3c, pmatrix)
383  CALL calculate_tvec_ri(lri_env, lri_env%o3c, para_env)
384  CALL calculate_avec_ri(lri_env, pmatrix)
385  !
386  CALL get_qs_env(qs_env, lri_density=lri_density, natom=natom)
387  IF (ASSOCIATED(lri_density)) THEN
388  CALL lri_density_release(lri_density)
389  ELSE
390  ALLOCATE (lri_density)
391  END IF
392  CALL lri_density_create(lri_density)
393  lri_density%nspin = nspin
394  ! allocate the arrays to hold RI expansion coefficients lri_coefs
395  CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
396  ! reassign avec
397  avec => lri_env%ri_fit%avec
398  bas_ptr => lri_env%ri_fit%bas_ptr
399  ALLOCATE (atom_of_kind(natom), kind_of(natom))
400  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
401  DO ispin = 1, nspin
402  lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
403  DO iatom = 1, natom
404  ikind = kind_of(iatom)
405  atom_a = atom_of_kind(iatom)
406  i1 = bas_ptr(1, iatom)
407  i2 = bas_ptr(2, iatom)
408  n = i2 - i1 + 1
409  lri_coef(ikind)%acoef(atom_a, 1:n) = avec(i1:i2, ispin)
410  END DO
411  END DO
412  CALL set_qs_env(qs_env, lri_density=lri_density)
413  DEALLOCATE (atom_of_kind, kind_of)
414  !
415  CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
416  DO ispin = 1, nspin
417  lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
418  CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
419  lri_coef, tot_rho_r(ispin), "RI_HXC", .false.)
420  END DO
421 
422  CALL timestop(handle)
423 
424  END SUBROUTINE calculate_ri_densities
425 
426 ! **************************************************************************************************
427 !> \brief assembles the global vector t=<P.T>
428 !> \param lri_env the lri environment
429 !> \param o3c ...
430 !> \param para_env ...
431 ! **************************************************************************************************
432  SUBROUTINE calculate_tvec_ri(lri_env, o3c, para_env)
433 
434  TYPE(lri_environment_type), POINTER :: lri_env
435  TYPE(o3c_container_type), POINTER :: o3c
436  TYPE(mp_para_env_type), POINTER :: para_env
437 
438  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_tvec_ri'
439  INTEGER, PARAMETER :: msweep = 32
440 
441  INTEGER :: handle, i1, i2, ibl, ibu, il, ispin, &
442  isweep, it, iu, katom, m, ma, mba, &
443  mepos, natom, nspin, nsweep, nthread
444  INTEGER, DIMENSION(2, msweep) :: nlimit
445  INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
446  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ta
447  REAL(kind=dp), DIMENSION(:, :), POINTER :: rm1t, tvec, tvl
448  TYPE(o3c_iterator_type) :: o3c_iterator
449 
450  CALL timeset(routinen, handle)
451 
452  nspin = SIZE(lri_env%ri_fit%tvec, 2)
453  bas_ptr => lri_env%ri_fit%bas_ptr
454  tvec => lri_env%ri_fit%tvec
455  tvec = 0.0_dp
456  natom = SIZE(bas_ptr, 2)
457  nthread = 1
458 !$ nthread = omp_get_max_threads()
459 
460  IF (natom < 1000) THEN
461  nsweep = 1
462  ELSE
463  nsweep = min(nthread, msweep)
464  END IF
465 
466  nlimit = 0
467  DO isweep = 1, nsweep
468  nlimit(1:2, isweep) = get_limit(natom, nsweep, isweep - 1)
469  END DO
470 
471  DO ispin = 1, nspin
472  DO isweep = 1, nsweep
473  il = nlimit(1, isweep)
474  iu = nlimit(2, isweep)
475  ma = iu - il + 1
476  IF (ma < 1) cycle
477  ibl = bas_ptr(1, il)
478  ibu = bas_ptr(2, iu)
479  mba = ibu - ibl + 1
480  ALLOCATE (ta(mba, nthread))
481  ta = 0.0_dp
482 
483  CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
484 
485 !$OMP PARALLEL DEFAULT(NONE)&
486 !$OMP SHARED (nthread,o3c_iterator,ispin,il,iu,ibl,ta,bas_ptr)&
487 !$OMP PRIVATE (mepos,katom,tvl,i1,i2,m)
488 
489  mepos = 0
490 !$ mepos = omp_get_thread_num()
491 
492  DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
493  CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, katom=katom, tvec=tvl)
494  IF (katom < il .OR. katom > iu) cycle
495  i1 = bas_ptr(1, katom) - ibl + 1
496  i2 = bas_ptr(2, katom) - ibl + 1
497  m = i2 - i1 + 1
498  ta(i1:i2, mepos + 1) = ta(i1:i2, mepos + 1) + tvl(1:m, ispin)
499  END DO
500 !$OMP END PARALLEL
501  CALL o3c_iterator_release(o3c_iterator)
502 
503  ! sum over threads
504  DO it = 1, nthread
505  tvec(ibl:ibu, ispin) = tvec(ibl:ibu, ispin) + ta(1:mba, it)
506  END DO
507  DEALLOCATE (ta)
508  END DO
509  END DO
510 
511  ! global summation
512  CALL para_env%sum(tvec)
513 
514  rm1t => lri_env%ri_fit%rm1t
515  DO ispin = 1, nspin
516  ! solve Rx=t
517  CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
518  vecr=tvec(:, ispin), &
519  vecx=rm1t(:, ispin), &
520  matp=lri_env%ri_sinv(1)%matrix, &
521  solver=lri_env%ri_sinv_app, &
522  ptr=bas_ptr)
523  END DO
524 
525  CALL timestop(handle)
526 
527  END SUBROUTINE calculate_tvec_ri
528 ! **************************************************************************************************
529 !> \brief performs the fitting of the density in the RI method
530 !> \param lri_env the lri environment
531 !> lri_density the environment for the fitting
532 !> pmatrix density matrix
533 !> \param pmatrix ...
534 ! **************************************************************************************************
535  SUBROUTINE calculate_avec_ri(lri_env, pmatrix)
536 
537  TYPE(lri_environment_type), POINTER :: lri_env
538  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
539 
540  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_avec_ri'
541 
542  INTEGER :: handle, ispin, nspin
543  REAL(kind=dp) :: etr, nelec, nrm1t
544  REAL(kind=dp), DIMENSION(:), POINTER :: nvec, rm1n
545  REAL(kind=dp), DIMENSION(:, :), POINTER :: avec, rm1t, tvec
546  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: smatrix
547 
548  CALL timeset(routinen, handle)
549 
550  nspin = SIZE(pmatrix)
551  ! number of electrons
552  smatrix => lri_env%ob_smat
553  DO ispin = 1, nspin
554  etr = 0.0_dp
555  CALL dbcsr_dot(smatrix(1)%matrix, pmatrix(ispin)%matrix, etr)
556  lri_env%ri_fit%echarge(ispin) = etr
557  END DO
558 
559  tvec => lri_env%ri_fit%tvec
560  rm1t => lri_env%ri_fit%rm1t
561  nvec => lri_env%ri_fit%nvec
562  rm1n => lri_env%ri_fit%rm1n
563 
564  ! calculate lambda
565  DO ispin = 1, nspin
566  nelec = lri_env%ri_fit%echarge(ispin)
567  nrm1t = sum(nvec(:)*rm1t(:, ispin))
568  lri_env%ri_fit%lambda(ispin) = 2.0_dp*(nrm1t - nelec)/lri_env%ri_fit%ntrm1n
569  END DO
570 
571  ! calculate avec = rm1t - lambda/2 * rm1n
572  avec => lri_env%ri_fit%avec
573  DO ispin = 1, nspin
574  avec(:, ispin) = rm1t(:, ispin) - 0.5_dp*lri_env%ri_fit%lambda(ispin)*rm1n(:)
575  END DO
576 
577  CALL timestop(handle)
578 
579  END SUBROUTINE calculate_avec_ri
580 
581 ! **************************************************************************************************
582 !> \brief Mutiplies a replicated vector with a DBCSR matrix: vo = mat*vi
583 !> \param mat ...
584 !> \param vi ...
585 !> \param vo ...
586 !> \param ptr ...
587 ! **************************************************************************************************
588  SUBROUTINE ri_matvec(mat, vi, vo, ptr)
589 
590  TYPE(dbcsr_type) :: mat
591  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vi
592  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: vo
593  INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr
594 
595  CHARACTER(LEN=*), PARAMETER :: routinen = 'ri_matvec'
596 
597  CHARACTER :: matrix_type
598  INTEGER :: group_handle, handle, iatom, jatom, m1, &
599  m2, mb, n1, n2, nb
600  LOGICAL :: symm
601  REAL(kind=dp), DIMENSION(:, :), POINTER :: block
602  TYPE(dbcsr_iterator_type) :: iter
603  TYPE(mp_comm_type) :: group
604 
605  CALL timeset(routinen, handle)
606 
607  CALL dbcsr_get_info(mat, matrix_type=matrix_type, group=group_handle)
608  CALL group%set_handle(group_handle)
609 
610  SELECT CASE (matrix_type)
611  CASE (dbcsr_type_no_symmetry)
612  symm = .false.
613  CASE (dbcsr_type_symmetric)
614  symm = .true.
615  CASE (dbcsr_type_antisymmetric)
616  cpabort("NYI, antisymmetric matrix not permitted")
617  CASE DEFAULT
618  cpabort("Unknown matrix type, ...")
619  END SELECT
620 
621  vo(:) = 0.0_dp
622  CALL dbcsr_iterator_start(iter, mat)
623  DO WHILE (dbcsr_iterator_blocks_left(iter))
624  CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
625  n1 = ptr(1, iatom)
626  n2 = ptr(2, iatom)
627  nb = n2 - n1 + 1
628  m1 = ptr(1, jatom)
629  m2 = ptr(2, jatom)
630  mb = m2 - m1 + 1
631  cpassert(nb == SIZE(block, 1))
632  cpassert(mb == SIZE(block, 2))
633  vo(n1:n2) = vo(n1:n2) + matmul(block, vi(m1:m2))
634  IF (symm .AND. (iatom /= jatom)) THEN
635  vo(m1:m2) = vo(m1:m2) + matmul(transpose(block), vi(n1:n2))
636  END IF
637  END DO
638  CALL dbcsr_iterator_stop(iter)
639 
640  CALL group%sum(vo)
641 
642  CALL timestop(handle)
643 
644  END SUBROUTINE ri_matvec
645 
646 END MODULE ri_environment_methods
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
arnoldi_api
arnoldi iteration using dbcsr
Definition: arnoldi_api.F:15
arnoldi_api::arnoldi_conjugate_gradient
subroutine, public arnoldi_conjugate_gradient(matrix_a, vec_x, matrix_p, converged, threshold, max_iter)
Wrapper for conjugated gradient algorithm for Ax=b.
Definition: arnoldi_api.F:303
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
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.
Definition: atomic_kind_types.F:223
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition: cp_dbcsr_diag.F:67
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
iterate_matrix
Routines useful for iterative matrix calculations.
Definition: iterate_matrix.F:13
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition: iterate_matrix.F:472
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
lri_environment_types::allocate_lri_coefs
subroutine, public allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
creates and initializes lri_coefs
Definition: lri_environment_types.F:974
lri_environment_types::lri_density_release
subroutine, public lri_density_release(lri_density)
releases the given lri_density
Definition: lri_environment_types.F:552
lri_environment_types::lri_density_create
subroutine, public lri_density_create(lri_density)
creates and initializes an lri_density environment
Definition: lri_environment_types.F:537
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_lri_rho_elec
subroutine, public calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
Collocates the fitted lri density on a grid.
Definition: qs_collocate_density.F:517
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
qs_ks_types
Definition: qs_ks_types.F:15
qs_o3c_methods
Methods used with 3-center overlap type integrals containers.
Definition: qs_o3c_methods.F:13
qs_o3c_methods::calculate_o3c_integrals
subroutine, public calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
...
Definition: qs_o3c_methods.F:48
qs_o3c_methods::contract12_o3c
subroutine, public contract12_o3c(o3c, matrix_p)
Contraction of 3-tensor over indices 1 and 2 (assuming symmetry) t(k) = sum_ij (ijk)*p(ij)
Definition: qs_o3c_methods.F:315
qs_o3c_types
3-center overlap type integrals containers
Definition: qs_o3c_types.F:12
qs_o3c_types::get_o3c_iterator_info
subroutine, public get_o3c_iterator_info(o3c_iterator, mepos, iatom, jatom, katom, ikind, jkind, kkind, rij, rik, cellj, cellk, integral, tvec, force_i, force_j, force_k)
...
Definition: qs_o3c_types.F:418
qs_o3c_types::o3c_iterator_create
subroutine, public o3c_iterator_create(o3c, o3c_iterator, nthread)
...
Definition: qs_o3c_types.F:357
qs_o3c_types::release_o3c_container
subroutine, public release_o3c_container(o3c_container)
...
Definition: qs_o3c_types.F:225
qs_o3c_types::o3c_iterator_release
subroutine, public o3c_iterator_release(o3c_iterator)
...
Definition: qs_o3c_types.F:384
qs_o3c_types::o3c_iterate
integer function, public o3c_iterate(o3c_iterator, mepos)
...
Definition: qs_o3c_types.F:473
qs_o3c_types::init_o3c_container
subroutine, public init_o3c_container(o3c, nspin, basis_set_list_a, basis_set_list_b, basis_set_list_c, sab_nl, sac_nl, only_bc_same_center)
...
Definition: qs_o3c_types.F:117
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:120
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
ri_environment_methods
Calculates integral matrices for RIGPW method.
Definition: ri_environment_methods.F:18
ri_environment_methods::build_ri_matrices
subroutine, public build_ri_matrices(lri_env, qs_env, calculate_forces)
creates and initializes an lri_env
Definition: ri_environment_methods.F:88
ri_environment_methods::ri_metric_solver
subroutine, public ri_metric_solver(mat, vecr, vecx, matp, solver, ptr)
solver for RI systems (R*x=n)
Definition: ri_environment_methods.F:289
ri_environment_methods::calculate_ri_densities
subroutine, public calculate_ri_densities(lri_env, qs_env, pmatrix, lri_rho_struct, atomic_kind_set, para_env)
performs the fitting of the density and distributes the fitted density on the grid
Definition: ri_environment_methods.F:358
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
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