Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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_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 :
81 : CONTAINS
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 4 : 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 4 : CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces)
96 :
97 4 : 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 4 : 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 4 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
118 : REAL(KIND=dp) :: fpre
119 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval
120 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout
121 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
122 : TYPE(dbcsr_iterator_type) :: iter
123 4 : 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 4 : CALL timeset(routineN, handle)
132 :
133 4 : CPASSERT(ASSOCIATED(lri_env))
134 4 : CPASSERT(ASSOCIATED(qs_env))
135 :
136 : ! overlap matrices
137 4 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
138 4 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
139 4 : 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 4 : IF (calculate_forces) THEN
143 : ! charge constraint forces are calculated here
144 2 : CALL get_qs_env(qs_env=qs_env, rho=rho)
145 2 : CALL qs_rho_get(rho, rho_ao=matrix_p)
146 2 : ALLOCATE (fmat)
147 2 : CALL dbcsr_create(fmat, template=matrix_p(1)%matrix)
148 4 : DO ispin = 1, dft_control%nspins
149 2 : fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n
150 4 : 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 2 : calculate_forces=.TRUE., matrix_p=fmat)
155 2 : CALL dbcsr_release(fmat)
156 2 : DEALLOCATE (fmat)
157 : ELSE
158 : CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
159 2 : 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 4 : basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list)
164 4 : IF (calculate_forces) THEN
165 : ! calculate f*a(T) pseudo density matrix for forces
166 2 : bas_ptr => lri_env%ri_fit%bas_ptr
167 2 : avec => lri_env%ri_fit%avec
168 2 : fout => lri_env%ri_fit%fout
169 2 : ALLOCATE (fmat)
170 2 : CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix)
171 2 : CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list)
172 2 : CALL dbcsr_set(fmat, 0.0_dp)
173 2 : CALL dbcsr_iterator_start(iter, fmat)
174 5 : DO WHILE (dbcsr_iterator_blocks_left(iter))
175 3 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk)
176 3 : i1 = bas_ptr(1, iatom)
177 3 : i2 = bas_ptr(2, iatom)
178 3 : j1 = bas_ptr(1, jatom)
179 3 : j2 = bas_ptr(2, jatom)
180 3 : IF (iatom <= jatom) THEN
181 6 : DO ispin = 1, dft_control%nspins
182 48 : DO j = j1, j2
183 42 : m = j - j1 + 1
184 633 : DO i = i1, i2
185 588 : n = i - i1 + 1
186 630 : 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 0 : DO ispin = 1, dft_control%nspins
192 0 : DO i = i1, i2
193 0 : n = i - i1 + 1
194 0 : DO j = j1, j2
195 0 : m = j - j1 + 1
196 0 : 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 5 : IF (iatom == jatom) THEN
202 422 : fblk(:, :) = 0.25_dp*fblk(:, :)
203 : ELSE
204 211 : fblk(:, :) = 0.5_dp*fblk(:, :)
205 : END IF
206 : END DO
207 2 : 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 2 : calculate_forces=.TRUE., matrix_p=fmat)
212 2 : CALL dbcsr_release(fmat)
213 2 : DEALLOCATE (fmat)
214 : END IF
215 : ! approximation (preconditioner) or exact inverse of RI overlap
216 4 : CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1)
217 4 : ALLOCATE (lri_env%ri_sinv(1)%matrix)
218 4 : SELECT CASE (lri_env%ri_sinv_app)
219 : CASE ("NONE")
220 : !
221 : CASE ("INVS")
222 0 : 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 0 : norm_convergence=1.e-10_dp, filter_eps=1.e-12_dp, silent=.FALSE.)
226 : CASE ("INVF")
227 4 : CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix)
228 4 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
229 4 : CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas)
230 12 : ALLOCATE (eval(nbas))
231 4 : CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env)
232 4 : izero = 0
233 116 : DO i = 1, nbas
234 116 : IF (eval(i) < 1.0e-10_dp) THEN
235 0 : eval(i) = 0.0_dp
236 0 : izero = izero + 1
237 : ELSE
238 112 : eval(i) = SQRT(1.0_dp/eval(i))
239 : END IF
240 : END DO
241 4 : CALL dbcsr_scale_by_vector(emat, eval, side='right')
242 4 : CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
243 4 : CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix)
244 4 : DEALLOCATE (eval)
245 8 : CALL dbcsr_release(emat)
246 : CASE ("AINV")
247 0 : 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 0 : norm_convergence=1.e-4_dp, filter_eps=1.e-4_dp, silent=.FALSE.)
251 : CASE DEFAULT
252 4 : 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 4 : ptr=lri_env%ri_fit%bas_ptr)
262 :
263 : ! calculate n(t)x
264 116 : lri_env%ri_fit%ntrm1n = SUM(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:))
265 :
266 : ! calculate 3c-overlap integrals
267 4 : IF (ASSOCIATED(lri_env%o3c)) THEN
268 2 : CALL release_o3c_container(lri_env%o3c)
269 : ELSE
270 2 : 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 4 : lri_env%soo_list, lri_env%soa_list)
275 :
276 4 : CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p)
277 :
278 4 : CALL timestop(handle)
279 :
280 8 : 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 56 : 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 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vect
305 :
306 56 : CALL timeset(routineN, handle)
307 :
308 56 : threshold = 1.e-10_dp
309 56 : max_iter = 100
310 :
311 1624 : vecx(:) = vecr(:)
312 :
313 0 : SELECT CASE (solver)
314 : CASE ("NONE")
315 : CALL arnoldi_conjugate_gradient(mat, vecx, &
316 0 : converged=converged, threshold=threshold, max_iter=max_iter)
317 : CASE ("INVS", "INVF")
318 56 : converged = .FALSE.
319 56 : CALL ri_matvec(matp, vecr, vecx, ptr)
320 : CASE ("AINV")
321 : CALL arnoldi_conjugate_gradient(mat, vecx, matrix_p=matp, &
322 0 : converged=converged, threshold=threshold, max_iter=max_iter)
323 : CASE DEFAULT
324 56 : CPABORT("Unknown RI solver")
325 : END SELECT
326 :
327 56 : IF (.NOT. converged) THEN
328 : ! get error
329 56 : rerror = 0.0_dp
330 56 : n = SIZE(vecr)
331 168 : ALLOCATE (vect(n))
332 56 : CALL ri_matvec(mat, vecx, vect, ptr)
333 1624 : vect(:) = vect(:) - vecr(:)
334 1680 : rerror = MAXVAL(ABS(vect(:)))
335 56 : DEALLOCATE (vect)
336 56 : CPWARN_IF(rerror > threshold, "RI solver: CG did not converge properly")
337 : END IF
338 :
339 56 : CALL timestop(handle)
340 :
341 56 : 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 26 : 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 26 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
371 26 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
372 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
373 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec
374 : TYPE(lri_density_type), POINTER :: lri_density
375 26 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
376 26 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
377 26 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
378 :
379 26 : CALL timeset(routineN, handle)
380 :
381 26 : nspin = SIZE(pmatrix, 1)
382 26 : CALL contract12_o3c(lri_env%o3c, pmatrix)
383 26 : CALL calculate_tvec_ri(lri_env, lri_env%o3c, para_env)
384 26 : CALL calculate_avec_ri(lri_env, pmatrix)
385 : !
386 26 : CALL get_qs_env(qs_env, lri_density=lri_density, natom=natom)
387 26 : IF (ASSOCIATED(lri_density)) THEN
388 24 : CALL lri_density_release(lri_density)
389 : ELSE
390 2 : ALLOCATE (lri_density)
391 : END IF
392 26 : CALL lri_density_create(lri_density)
393 26 : lri_density%nspin = nspin
394 : ! allocate the arrays to hold RI expansion coefficients lri_coefs
395 26 : CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
396 : ! reassign avec
397 26 : avec => lri_env%ri_fit%avec
398 26 : bas_ptr => lri_env%ri_fit%bas_ptr
399 26 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
400 52 : DO ispin = 1, nspin
401 26 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
402 104 : DO iatom = 1, natom
403 52 : ikind = kind_of(iatom)
404 52 : atom_a = atom_of_kind(iatom)
405 52 : i1 = bas_ptr(1, iatom)
406 52 : i2 = bas_ptr(2, iatom)
407 52 : n = i2 - i1 + 1
408 1586 : lri_coef(ikind)%acoef(atom_a, 1:n) = avec(i1:i2, ispin)
409 : END DO
410 : END DO
411 26 : CALL set_qs_env(qs_env, lri_density=lri_density)
412 26 : DEALLOCATE (atom_of_kind, kind_of)
413 : !
414 26 : CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
415 52 : DO ispin = 1, nspin
416 26 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
417 : CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
418 52 : lri_coef, tot_rho_r(ispin), "RI_HXC", .FALSE.)
419 : END DO
420 :
421 26 : CALL timestop(handle)
422 :
423 26 : END SUBROUTINE calculate_ri_densities
424 :
425 : ! **************************************************************************************************
426 : !> \brief assembles the global vector t=<P.T>
427 : !> \param lri_env the lri environment
428 : !> \param o3c ...
429 : !> \param para_env ...
430 : ! **************************************************************************************************
431 26 : SUBROUTINE calculate_tvec_ri(lri_env, o3c, para_env)
432 :
433 : TYPE(lri_environment_type), POINTER :: lri_env
434 : TYPE(o3c_container_type), POINTER :: o3c
435 : TYPE(mp_para_env_type), POINTER :: para_env
436 :
437 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_tvec_ri'
438 : INTEGER, PARAMETER :: msweep = 32
439 :
440 : INTEGER :: handle, i1, i2, ibl, ibu, il, ispin, &
441 : isweep, it, iu, katom, m, ma, mba, &
442 : mepos, natom, nspin, nsweep, nthread
443 : INTEGER, DIMENSION(2, msweep) :: nlimit
444 26 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
445 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ta
446 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rm1t, tvec, tvl
447 : TYPE(o3c_iterator_type) :: o3c_iterator
448 :
449 26 : CALL timeset(routineN, handle)
450 :
451 26 : nspin = SIZE(lri_env%ri_fit%tvec, 2)
452 26 : bas_ptr => lri_env%ri_fit%bas_ptr
453 26 : tvec => lri_env%ri_fit%tvec
454 780 : tvec = 0.0_dp
455 26 : natom = SIZE(bas_ptr, 2)
456 : nthread = 1
457 26 : !$ nthread = omp_get_max_threads()
458 :
459 26 : IF (natom < 1000) THEN
460 26 : nsweep = 1
461 : ELSE
462 0 : nsweep = MIN(nthread, msweep)
463 : END IF
464 :
465 26 : nlimit = 0
466 52 : DO isweep = 1, nsweep
467 52 : nlimit(1:2, isweep) = get_limit(natom, nsweep, isweep - 1)
468 : END DO
469 :
470 52 : DO ispin = 1, nspin
471 78 : DO isweep = 1, nsweep
472 26 : il = nlimit(1, isweep)
473 26 : iu = nlimit(2, isweep)
474 26 : ma = iu - il + 1
475 26 : IF (ma < 1) CYCLE
476 26 : ibl = bas_ptr(1, il)
477 26 : ibu = bas_ptr(2, iu)
478 26 : mba = ibu - ibl + 1
479 104 : ALLOCATE (ta(mba, nthread))
480 780 : ta = 0.0_dp
481 :
482 26 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
483 :
484 : !$OMP PARALLEL DEFAULT(NONE)&
485 : !$OMP SHARED (nthread,o3c_iterator,ispin,il,iu,ibl,ta,bas_ptr)&
486 26 : !$OMP PRIVATE (mepos,katom,tvl,i1,i2,m)
487 :
488 : mepos = 0
489 : !$ mepos = omp_get_thread_num()
490 :
491 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
492 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, katom=katom, tvec=tvl)
493 : IF (katom < il .OR. katom > iu) CYCLE
494 : i1 = bas_ptr(1, katom) - ibl + 1
495 : i2 = bas_ptr(2, katom) - ibl + 1
496 : m = i2 - i1 + 1
497 : ta(i1:i2, mepos + 1) = ta(i1:i2, mepos + 1) + tvl(1:m, ispin)
498 : END DO
499 : !$OMP END PARALLEL
500 26 : CALL o3c_iterator_release(o3c_iterator)
501 :
502 : ! sum over threads
503 52 : DO it = 1, nthread
504 780 : tvec(ibl:ibu, ispin) = tvec(ibl:ibu, ispin) + ta(1:mba, it)
505 : END DO
506 52 : DEALLOCATE (ta)
507 : END DO
508 : END DO
509 :
510 : ! global summation
511 1534 : CALL para_env%sum(tvec)
512 :
513 26 : rm1t => lri_env%ri_fit%rm1t
514 52 : DO ispin = 1, nspin
515 : ! solve Rx=t
516 : CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
517 : vecr=tvec(:, ispin), &
518 : vecx=rm1t(:, ispin), &
519 : matp=lri_env%ri_sinv(1)%matrix, &
520 : solver=lri_env%ri_sinv_app, &
521 52 : ptr=bas_ptr)
522 : END DO
523 :
524 26 : CALL timestop(handle)
525 :
526 52 : END SUBROUTINE calculate_tvec_ri
527 : ! **************************************************************************************************
528 : !> \brief performs the fitting of the density in the RI method
529 : !> \param lri_env the lri environment
530 : !> lri_density the environment for the fitting
531 : !> pmatrix density matrix
532 : !> \param pmatrix ...
533 : ! **************************************************************************************************
534 26 : SUBROUTINE calculate_avec_ri(lri_env, pmatrix)
535 :
536 : TYPE(lri_environment_type), POINTER :: lri_env
537 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
538 :
539 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_avec_ri'
540 :
541 : INTEGER :: handle, ispin, nspin
542 : REAL(KIND=dp) :: etr, nelec, nrm1t
543 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: nvec, rm1n
544 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, rm1t, tvec
545 26 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: smatrix
546 :
547 26 : CALL timeset(routineN, handle)
548 :
549 26 : nspin = SIZE(pmatrix)
550 : ! number of electrons
551 26 : smatrix => lri_env%ob_smat
552 52 : DO ispin = 1, nspin
553 : etr = 0.0_dp
554 26 : CALL dbcsr_dot(smatrix(1)%matrix, pmatrix(ispin)%matrix, etr)
555 52 : lri_env%ri_fit%echarge(ispin) = etr
556 : END DO
557 :
558 26 : tvec => lri_env%ri_fit%tvec
559 26 : rm1t => lri_env%ri_fit%rm1t
560 26 : nvec => lri_env%ri_fit%nvec
561 26 : rm1n => lri_env%ri_fit%rm1n
562 :
563 : ! calculate lambda
564 52 : DO ispin = 1, nspin
565 26 : nelec = lri_env%ri_fit%echarge(ispin)
566 754 : nrm1t = SUM(nvec(:)*rm1t(:, ispin))
567 52 : lri_env%ri_fit%lambda(ispin) = 2.0_dp*(nrm1t - nelec)/lri_env%ri_fit%ntrm1n
568 : END DO
569 :
570 : ! calculate avec = rm1t - lambda/2 * rm1n
571 26 : avec => lri_env%ri_fit%avec
572 52 : DO ispin = 1, nspin
573 1534 : avec(:, ispin) = rm1t(:, ispin) - 0.5_dp*lri_env%ri_fit%lambda(ispin)*rm1n(:)
574 : END DO
575 :
576 26 : CALL timestop(handle)
577 :
578 26 : END SUBROUTINE calculate_avec_ri
579 :
580 : ! **************************************************************************************************
581 : !> \brief Mutiplies a replicated vector with a DBCSR matrix: vo = mat*vi
582 : !> \param mat ...
583 : !> \param vi ...
584 : !> \param vo ...
585 : !> \param ptr ...
586 : ! **************************************************************************************************
587 112 : SUBROUTINE ri_matvec(mat, vi, vo, ptr)
588 :
589 : TYPE(dbcsr_type) :: mat
590 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vi
591 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vo
592 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr
593 :
594 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_matvec'
595 :
596 : CHARACTER :: matrix_type
597 : INTEGER :: handle, iatom, jatom, m1, m2, mb, n1, &
598 : n2, nb
599 : LOGICAL :: symm
600 112 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
601 : TYPE(dbcsr_iterator_type) :: iter
602 : TYPE(mp_comm_type) :: group
603 :
604 112 : CALL timeset(routineN, handle)
605 :
606 112 : CALL dbcsr_get_info(mat, matrix_type=matrix_type, group=group)
607 :
608 : SELECT CASE (matrix_type)
609 : CASE (dbcsr_type_no_symmetry)
610 112 : symm = .FALSE.
611 : CASE (dbcsr_type_symmetric)
612 112 : symm = .TRUE.
613 : CASE (dbcsr_type_antisymmetric)
614 0 : CPABORT("NYI, antisymmetric matrix not permitted")
615 : CASE DEFAULT
616 112 : CPABORT("Unknown matrix type, ...")
617 : END SELECT
618 :
619 3248 : vo(:) = 0.0_dp
620 112 : CALL dbcsr_iterator_start(iter, mat)
621 280 : DO WHILE (dbcsr_iterator_blocks_left(iter))
622 168 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
623 168 : n1 = ptr(1, iatom)
624 168 : n2 = ptr(2, iatom)
625 168 : nb = n2 - n1 + 1
626 168 : m1 = ptr(1, jatom)
627 168 : m2 = ptr(2, jatom)
628 168 : mb = m2 - m1 + 1
629 168 : CPASSERT(nb == SIZE(block, 1))
630 168 : CPASSERT(mb == SIZE(block, 2))
631 73584 : vo(n1:n2) = vo(n1:n2) + MATMUL(block, vi(m1:m2))
632 336 : IF (symm .AND. (iatom /= jatom)) THEN
633 840 : vo(m1:m2) = vo(m1:m2) + MATMUL(TRANSPOSE(block), vi(n1:n2))
634 : END IF
635 : END DO
636 112 : CALL dbcsr_iterator_stop(iter)
637 :
638 6384 : CALL group%sum(vo)
639 :
640 112 : CALL timestop(handle)
641 :
642 112 : END SUBROUTINE ri_matvec
643 :
644 224 : END MODULE ri_environment_methods
|
