(git:374b731)
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-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! **************************************************************************************************
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_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
79CONTAINS
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
646END 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::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_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_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
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:559
lri_environment_types::lri_density_type
Definition lri_environment_types.F:357
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:269
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:330
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:215
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:135
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