(git:d18deda)
Loading...
Searching...
No Matches
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-2025 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! **************************************************************************************************
18MODULE 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_api, ONLY: &
25 dbcsr_add, dbcsr_create, dbcsr_get_info, dbcsr_iterator_blocks_left, &
26 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
27 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, &
28 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
29 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
30 dbcsr_scale_by_vector
31 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
32 USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
33 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
34 USE iterate_matrix, ONLY: invert_hotelling
35 USE kinds, ONLY: dp
36 USE lri_environment_types, ONLY: allocate_lri_coefs,&
37 lri_density_create,&
38 lri_density_release,&
39 lri_density_type,&
40 lri_environment_type,&
41 lri_kind_type
42 USE message_passing, ONLY: mp_comm_type,&
43 mp_para_env_type
44 USE pw_types, ONLY: pw_c1d_gs_type,&
45 pw_r3d_rs_type
46 USE qs_collocate_density, ONLY: calculate_lri_rho_elec
47 USE qs_environment_types, ONLY: get_qs_env,&
48 qs_environment_type,&
49 set_qs_env
50 USE qs_ks_types, ONLY: qs_ks_env_type
51 USE qs_o3c_methods, ONLY: calculate_o3c_integrals,&
52 contract12_o3c
53 USE qs_o3c_types, ONLY: get_o3c_iterator_info,&
54 init_o3c_container,&
55 o3c_container_type,&
56 o3c_iterate,&
57 o3c_iterator_create,&
58 o3c_iterator_release,&
59 o3c_iterator_type,&
60 release_o3c_container
61 USE qs_overlap, ONLY: build_overlap_matrix
62 USE qs_rho_types, ONLY: qs_rho_get,&
63 qs_rho_type
64 USE util, ONLY: get_limit
65
66!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
67#include "./base/base_uses.f90"
68
69 IMPLICIT NONE
70
71 PRIVATE
72
73! **************************************************************************************************
74
75 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_environment_methods'
76
77 PUBLIC :: build_ri_matrices, calculate_ri_densities, ri_metric_solver
78
79! **************************************************************************************************
80
81CONTAINS
82
83! **************************************************************************************************
84!> \brief creates and initializes an lri_env
85!> \param lri_env the lri_environment you want to create
86!> \param qs_env ...
87!> \param calculate_forces ...
88! **************************************************************************************************
89 SUBROUTINE build_ri_matrices(lri_env, qs_env, calculate_forces)
90
91 TYPE(lri_environment_type), POINTER :: lri_env
92 TYPE(qs_environment_type), POINTER :: qs_env
93 LOGICAL, INTENT(IN) :: calculate_forces
94
95 CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces)
96
97 END SUBROUTINE build_ri_matrices
98
99! **************************************************************************************************
100!> \brief calculates integrals needed for the RI density fitting,
101!> integrals are calculated once, before the SCF starts
102!> \param lri_env ...
103!> \param qs_env ...
104!> \param calculate_forces ...
105! **************************************************************************************************
106 SUBROUTINE calculate_ri_integrals(lri_env, qs_env, calculate_forces)
107
108 TYPE(lri_environment_type), POINTER :: lri_env
109 TYPE(qs_environment_type), POINTER :: qs_env
110 LOGICAL, INTENT(IN) :: calculate_forces
111
112 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_ri_integrals'
113 REAL(kind=dp), DIMENSION(2), PARAMETER :: fx = (/0.0_dp, 1.0_dp/)
114
115 INTEGER :: handle, i, i1, i2, iatom, ispin, izero, &
116 j, j1, j2, jatom, m, n, nbas
117 INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
118 REAL(kind=dp) :: fpre
119 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eval
120 REAL(kind=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout
121 TYPE(cp_blacs_env_type), POINTER :: blacs_env
122 TYPE(dbcsr_iterator_type) :: iter
123 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
124 TYPE(dbcsr_type) :: emat
125 TYPE(dbcsr_type), POINTER :: fmat
126 TYPE(dft_control_type), POINTER :: dft_control
127 TYPE(mp_para_env_type), POINTER :: para_env
128 TYPE(qs_ks_env_type), POINTER :: ks_env
129 TYPE(qs_rho_type), POINTER :: rho
130
131 CALL timeset(routinen, handle)
132
133 cpassert(ASSOCIATED(lri_env))
134 cpassert(ASSOCIATED(qs_env))
135
136 ! overlap matrices
137 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
138 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
139 NULLIFY (rho, matrix_p)
140 ! orbital overlap matrix (needed for N constraints and forces)
141 ! we replicate this in order to not directly interact qith qs_env
142 IF (calculate_forces) THEN
143 ! charge constraint forces are calculated here
144 CALL get_qs_env(qs_env=qs_env, rho=rho)
145 CALL qs_rho_get(rho, rho_ao=matrix_p)
146 ALLOCATE (fmat)
147 CALL dbcsr_create(fmat, template=matrix_p(1)%matrix)
148 DO ispin = 1, dft_control%nspins
149 fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n
150 CALL dbcsr_add(fmat, matrix_p(ispin)%matrix, fx(ispin), -fpre)
151 END DO
152 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
153 basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list, &
154 calculate_forces=.true., matrix_p=fmat)
155 CALL dbcsr_release(fmat)
156 DEALLOCATE (fmat)
157 ELSE
158 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
159 basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list)
160 END IF
161 ! RI overlap matrix
162 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, &
163 basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list)
164 IF (calculate_forces) THEN
165 ! calculate f*a(T) pseudo density matrix for forces
166 bas_ptr => lri_env%ri_fit%bas_ptr
167 avec => lri_env%ri_fit%avec
168 fout => lri_env%ri_fit%fout
169 ALLOCATE (fmat)
170 CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix)
171 CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list)
172 CALL dbcsr_set(fmat, 0.0_dp)
173 CALL dbcsr_iterator_start(iter, fmat)
174 DO WHILE (dbcsr_iterator_blocks_left(iter))
175 CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk)
176 i1 = bas_ptr(1, iatom)
177 i2 = bas_ptr(2, iatom)
178 j1 = bas_ptr(1, jatom)
179 j2 = bas_ptr(2, jatom)
180 IF (iatom <= jatom) THEN
181 DO ispin = 1, dft_control%nspins
182 DO j = j1, j2
183 m = j - j1 + 1
184 DO i = i1, i2
185 n = i - i1 + 1
186 fblk(n, m) = fblk(n, m) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin)
187 END DO
188 END DO
189 END DO
190 ELSE
191 DO ispin = 1, dft_control%nspins
192 DO i = i1, i2
193 n = i - i1 + 1
194 DO j = j1, j2
195 m = j - j1 + 1
196 fblk(m, n) = fblk(m, n) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin)
197 END DO
198 END DO
199 END DO
200 END IF
201 IF (iatom == jatom) THEN
202 fblk(:, :) = 0.25_dp*fblk(:, :)
203 ELSE
204 fblk(:, :) = 0.5_dp*fblk(:, :)
205 END IF
206 END DO
207 CALL dbcsr_iterator_stop(iter)
208 !
209 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, &
210 basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list, &
211 calculate_forces=.true., matrix_p=fmat)
212 CALL dbcsr_release(fmat)
213 DEALLOCATE (fmat)
214 END IF
215 ! approximation (preconditioner) or exact inverse of RI overlap
216 CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1)
217 ALLOCATE (lri_env%ri_sinv(1)%matrix)
218 SELECT CASE (lri_env%ri_sinv_app)
219 CASE ("NONE")
220 !
221 CASE ("INVS")
222 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
223 CALL invert_hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, &
224 threshold=1.e-10_dp, use_inv_as_guess=.false., &
225 norm_convergence=1.e-10_dp, filter_eps=1.e-12_dp, silent=.false.)
226 CASE ("INVF")
227 CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix)
228 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
229 CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas)
230 ALLOCATE (eval(nbas))
231 CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env)
232 izero = 0
233 DO i = 1, nbas
234 IF (eval(i) < 1.0e-10_dp) THEN
235 eval(i) = 0.0_dp
236 izero = izero + 1
237 ELSE
238 eval(i) = sqrt(1.0_dp/eval(i))
239 END IF
240 END DO
241 CALL dbcsr_scale_by_vector(emat, eval, side='right')
242 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
243 CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix)
244 DEALLOCATE (eval)
245 CALL dbcsr_release(emat)
246 CASE ("AINV")
247 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
248 CALL invert_hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, &
249 threshold=1.e-5_dp, use_inv_as_guess=.false., &
250 norm_convergence=1.e-4_dp, filter_eps=1.e-4_dp, silent=.false.)
251 CASE DEFAULT
252 cpabort("Unknown RI_SINV type")
253 END SELECT
254
255 ! solve Rx=n
256 CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
257 vecr=lri_env%ri_fit%nvec, &
258 vecx=lri_env%ri_fit%rm1n, &
259 matp=lri_env%ri_sinv(1)%matrix, &
260 solver=lri_env%ri_sinv_app, &
261 ptr=lri_env%ri_fit%bas_ptr)
262
263 ! calculate n(t)x
264 lri_env%ri_fit%ntrm1n = sum(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:))
265
266 ! calculate 3c-overlap integrals
267 IF (ASSOCIATED(lri_env%o3c)) THEN
268 CALL release_o3c_container(lri_env%o3c)
269 ELSE
270 ALLOCATE (lri_env%o3c)
271 END IF
272 CALL init_o3c_container(lri_env%o3c, dft_control%nspins, &
273 lri_env%orb_basis, lri_env%orb_basis, lri_env%ri_basis, &
274 lri_env%soo_list, lri_env%soa_list)
275
276 CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p)
277
278 CALL timestop(handle)
279
280 END SUBROUTINE calculate_ri_integrals
281! **************************************************************************************************
282!> \brief solver for RI systems (R*x=n)
283!> \param mat ...
284!> \param vecr ...
285!> \param vecx ...
286!> \param matp ...
287!> \param solver ...
288!> \param ptr ...
289! **************************************************************************************************
290 SUBROUTINE ri_metric_solver(mat, vecr, vecx, matp, solver, ptr)
291
292 TYPE(dbcsr_type) :: mat
293 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vecr
294 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: vecx
295 TYPE(dbcsr_type) :: matp
296 CHARACTER(LEN=*), INTENT(IN) :: solver
297 INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr
298
299 CHARACTER(LEN=*), PARAMETER :: routinen = 'ri_metric_solver'
300
301 INTEGER :: handle, max_iter, n
302 LOGICAL :: converged
303 REAL(kind=dp) :: rerror, threshold
304 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vect
305
306 CALL timeset(routinen, handle)
307
308 threshold = 1.e-10_dp
309 max_iter = 100
310
311 vecx(:) = vecr(:)
312
313 SELECT CASE (solver)
314 CASE ("NONE")
315 CALL arnoldi_conjugate_gradient(mat, vecx, &
316 converged=converged, threshold=threshold, max_iter=max_iter)
317 CASE ("INVS", "INVF")
318 converged = .false.
319 CALL ri_matvec(matp, vecr, vecx, ptr)
320 CASE ("AINV")
321 CALL arnoldi_conjugate_gradient(mat, vecx, matrix_p=matp, &
322 converged=converged, threshold=threshold, max_iter=max_iter)
323 CASE DEFAULT
324 cpabort("Unknown RI solver")
325 END SELECT
326
327 IF (.NOT. converged) THEN
328 ! get error
329 rerror = 0.0_dp
330 n = SIZE(vecr)
331 ALLOCATE (vect(n))
332 CALL ri_matvec(mat, vecx, vect, ptr)
333 vect(:) = vect(:) - vecr(:)
334 rerror = maxval(abs(vect(:)))
335 DEALLOCATE (vect)
336 cpwarn_if(rerror > threshold, "RI solver: CG did not converge properly")
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 :: handle, iatom, jatom, m1, m2, mb, n1, &
599 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)
608
609 SELECT CASE (matrix_type)
610 CASE (dbcsr_type_no_symmetry)
611 symm = .false.
612 CASE (dbcsr_type_symmetric)
613 symm = .true.
614 CASE (dbcsr_type_antisymmetric)
615 cpabort("NYI, antisymmetric matrix not permitted")
616 CASE DEFAULT
617 cpabort("Unknown matrix type, ...")
618 END SELECT
619
620 vo(:) = 0.0_dp
621 CALL dbcsr_iterator_start(iter, mat)
622 DO WHILE (dbcsr_iterator_blocks_left(iter))
623 CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
624 n1 = ptr(1, iatom)
625 n2 = ptr(2, iatom)
626 nb = n2 - n1 + 1
627 m1 = ptr(1, jatom)
628 m2 = ptr(2, jatom)
629 mb = m2 - m1 + 1
630 cpassert(nb == SIZE(block, 1))
631 cpassert(mb == SIZE(block, 2))
632 vo(n1:n2) = vo(n1:n2) + matmul(block, vi(m1:m2))
633 IF (symm .AND. (iatom /= jatom)) THEN
634 vo(m1:m2) = vo(m1:m2) + matmul(transpose(block), vi(n1:n2))
635 END IF
636 END DO
637 CALL dbcsr_iterator_stop(iter)
638
639 CALL group%sum(vo)
640
641 CALL timestop(handle)
642
643 END SUBROUTINE ri_matvec
644
645END MODULE ri_environment_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1073
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:794
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_scale_by_vector
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
Definition cp_dbcsr_contrib.F:486
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:473
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:976
lri_environment_types::lri_density_release
subroutine, public lri_density_release(lri_density)
releases the given lri_density
Definition lri_environment_types.F:554
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:539
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:518
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, harris_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, eeq, rhs)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1055
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:90
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:291
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
lri_environment_types::lri_density_type
Definition lri_environment_types.F:359
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_o3c_types::o3c_container_type
Definition qs_o3c_types.F:56
qs_o3c_types::o3c_iterator_type
Definition qs_o3c_types.F:75
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62