Line data Source code
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 : ! **************************************************************************************************
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 0 : 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 0 : CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces)
96 :
97 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
118 : REAL(KIND=dp) :: fpre
119 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval
120 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout
121 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
122 : TYPE(dbcsr_iterator_type) :: iter
123 0 : 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 0 : CALL timeset(routineN, handle)
132 :
133 0 : CPASSERT(ASSOCIATED(lri_env))
134 0 : CPASSERT(ASSOCIATED(qs_env))
135 :
136 : ! overlap matrices
137 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
138 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
139 0 : 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 0 : IF (calculate_forces) THEN
143 : ! charge constraint forces are calculated here
144 0 : CALL get_qs_env(qs_env=qs_env, rho=rho)
145 0 : CALL qs_rho_get(rho, rho_ao=matrix_p)
146 0 : ALLOCATE (fmat)
147 0 : CALL dbcsr_create(fmat, template=matrix_p(1)%matrix)
148 0 : DO ispin = 1, dft_control%nspins
149 0 : fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n
150 0 : 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 0 : calculate_forces=.TRUE., matrix_p=fmat)
155 0 : CALL dbcsr_release(fmat)
156 0 : DEALLOCATE (fmat)
157 : ELSE
158 : CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
159 0 : 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 0 : basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list)
164 0 : IF (calculate_forces) THEN
165 : ! calculate f*a(T) pseudo density matrix for forces
166 0 : bas_ptr => lri_env%ri_fit%bas_ptr
167 0 : avec => lri_env%ri_fit%avec
168 0 : fout => lri_env%ri_fit%fout
169 0 : ALLOCATE (fmat)
170 0 : CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix)
171 0 : CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list)
172 0 : CALL dbcsr_set(fmat, 0.0_dp)
173 0 : CALL dbcsr_iterator_start(iter, fmat)
174 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
175 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk)
176 0 : i1 = bas_ptr(1, iatom)
177 0 : i2 = bas_ptr(2, iatom)
178 0 : j1 = bas_ptr(1, jatom)
179 0 : j2 = bas_ptr(2, jatom)
180 0 : IF (iatom <= jatom) THEN
181 0 : DO ispin = 1, dft_control%nspins
182 0 : DO j = j1, j2
183 0 : m = j - j1 + 1
184 0 : DO i = i1, i2
185 0 : n = i - i1 + 1
186 0 : 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 0 : IF (iatom == jatom) THEN
202 0 : fblk(:, :) = 0.25_dp*fblk(:, :)
203 : ELSE
204 0 : fblk(:, :) = 0.5_dp*fblk(:, :)
205 : END IF
206 : END DO
207 0 : 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 0 : calculate_forces=.TRUE., matrix_p=fmat)
212 0 : CALL dbcsr_release(fmat)
213 0 : DEALLOCATE (fmat)
214 : END IF
215 : ! approximation (preconditioner) or exact inverse of RI overlap
216 0 : CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1)
217 0 : ALLOCATE (lri_env%ri_sinv(1)%matrix)
218 0 : 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 0 : CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix)
228 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
229 0 : CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas)
230 0 : ALLOCATE (eval(nbas))
231 0 : CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env)
232 0 : izero = 0
233 0 : DO i = 1, nbas
234 0 : IF (eval(i) < 1.0e-10_dp) THEN
235 0 : eval(i) = 0.0_dp
236 0 : izero = izero + 1
237 : ELSE
238 0 : eval(i) = SQRT(1.0_dp/eval(i))
239 : END IF
240 : END DO
241 0 : CALL dbcsr_scale_by_vector(emat, eval, side='right')
242 0 : CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
243 0 : CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix)
244 0 : DEALLOCATE (eval)
245 0 : 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 0 : 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 0 : ptr=lri_env%ri_fit%bas_ptr)
262 :
263 : ! calculate n(t)x
264 0 : lri_env%ri_fit%ntrm1n = SUM(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:))
265 :
266 : ! calculate 3c-overlap integrals
267 0 : IF (ASSOCIATED(lri_env%o3c)) THEN
268 0 : CALL release_o3c_container(lri_env%o3c)
269 : ELSE
270 0 : 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 0 : lri_env%soo_list, lri_env%soa_list)
275 :
276 0 : CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p)
277 :
278 0 : CALL timestop(handle)
279 :
280 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vect
305 :
306 0 : CALL timeset(routineN, handle)
307 :
308 0 : threshold = 1.e-10_dp
309 0 : max_iter = 100
310 :
311 0 : 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 0 : converged = .FALSE.
319 0 : 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 0 : CPABORT("Unknown RI solver")
325 : END SELECT
326 :
327 0 : IF (.NOT. converged) THEN
328 : ! get error
329 0 : rerror = 0.0_dp
330 0 : n = SIZE(vecr)
331 0 : ALLOCATE (vect(n))
332 0 : CALL ri_matvec(mat, vecx, vect, ptr)
333 0 : vect(:) = vect(:) - vecr(:)
334 0 : rerror = MAXVAL(ABS(vect(:)))
335 0 : DEALLOCATE (vect)
336 0 : CPWARN_IF(rerror > threshold, "RI solver: CG did not converge properly")
337 : END IF
338 :
339 0 : CALL timestop(handle)
340 :
341 0 : 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 0 : 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 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
371 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
372 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
373 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec
374 : TYPE(lri_density_type), POINTER :: lri_density
375 0 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
376 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
377 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
378 :
379 0 : CALL timeset(routineN, handle)
380 :
381 0 : nspin = SIZE(pmatrix, 1)
382 0 : CALL contract12_o3c(lri_env%o3c, pmatrix)
383 0 : CALL calculate_tvec_ri(lri_env, lri_env%o3c, para_env)
384 0 : CALL calculate_avec_ri(lri_env, pmatrix)
385 : !
386 0 : CALL get_qs_env(qs_env, lri_density=lri_density, natom=natom)
387 0 : IF (ASSOCIATED(lri_density)) THEN
388 0 : CALL lri_density_release(lri_density)
389 : ELSE
390 0 : ALLOCATE (lri_density)
391 : END IF
392 0 : CALL lri_density_create(lri_density)
393 0 : lri_density%nspin = nspin
394 : ! allocate the arrays to hold RI expansion coefficients lri_coefs
395 0 : CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
396 : ! reassign avec
397 0 : avec => lri_env%ri_fit%avec
398 0 : bas_ptr => lri_env%ri_fit%bas_ptr
399 0 : ALLOCATE (atom_of_kind(natom), kind_of(natom))
400 0 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
401 0 : DO ispin = 1, nspin
402 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
403 0 : DO iatom = 1, natom
404 0 : ikind = kind_of(iatom)
405 0 : atom_a = atom_of_kind(iatom)
406 0 : i1 = bas_ptr(1, iatom)
407 0 : i2 = bas_ptr(2, iatom)
408 0 : n = i2 - i1 + 1
409 0 : lri_coef(ikind)%acoef(atom_a, 1:n) = avec(i1:i2, ispin)
410 : END DO
411 : END DO
412 0 : CALL set_qs_env(qs_env, lri_density=lri_density)
413 0 : DEALLOCATE (atom_of_kind, kind_of)
414 : !
415 0 : CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
416 0 : DO ispin = 1, nspin
417 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
418 : CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
419 0 : lri_coef, tot_rho_r(ispin), "RI_HXC", .FALSE.)
420 : END DO
421 :
422 0 : CALL timestop(handle)
423 :
424 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
446 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ta
447 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rm1t, tvec, tvl
448 : TYPE(o3c_iterator_type) :: o3c_iterator
449 :
450 0 : CALL timeset(routineN, handle)
451 :
452 0 : nspin = SIZE(lri_env%ri_fit%tvec, 2)
453 0 : bas_ptr => lri_env%ri_fit%bas_ptr
454 0 : tvec => lri_env%ri_fit%tvec
455 0 : tvec = 0.0_dp
456 0 : natom = SIZE(bas_ptr, 2)
457 : nthread = 1
458 0 : !$ nthread = omp_get_max_threads()
459 :
460 0 : IF (natom < 1000) THEN
461 0 : nsweep = 1
462 : ELSE
463 0 : nsweep = MIN(nthread, msweep)
464 : END IF
465 :
466 0 : nlimit = 0
467 0 : DO isweep = 1, nsweep
468 0 : nlimit(1:2, isweep) = get_limit(natom, nsweep, isweep - 1)
469 : END DO
470 :
471 0 : DO ispin = 1, nspin
472 0 : DO isweep = 1, nsweep
473 0 : il = nlimit(1, isweep)
474 0 : iu = nlimit(2, isweep)
475 0 : ma = iu - il + 1
476 0 : IF (ma < 1) CYCLE
477 0 : ibl = bas_ptr(1, il)
478 0 : ibu = bas_ptr(2, iu)
479 0 : mba = ibu - ibl + 1
480 0 : ALLOCATE (ta(mba, nthread))
481 0 : ta = 0.0_dp
482 :
483 0 : 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 0 : !$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 0 : CALL o3c_iterator_release(o3c_iterator)
502 :
503 : ! sum over threads
504 0 : DO it = 1, nthread
505 0 : tvec(ibl:ibu, ispin) = tvec(ibl:ibu, ispin) + ta(1:mba, it)
506 : END DO
507 0 : DEALLOCATE (ta)
508 : END DO
509 : END DO
510 :
511 : ! global summation
512 0 : CALL para_env%sum(tvec)
513 :
514 0 : rm1t => lri_env%ri_fit%rm1t
515 0 : 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 0 : ptr=bas_ptr)
523 : END DO
524 :
525 0 : CALL timestop(handle)
526 :
527 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: nvec, rm1n
545 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, rm1t, tvec
546 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: smatrix
547 :
548 0 : CALL timeset(routineN, handle)
549 :
550 0 : nspin = SIZE(pmatrix)
551 : ! number of electrons
552 0 : smatrix => lri_env%ob_smat
553 0 : DO ispin = 1, nspin
554 : etr = 0.0_dp
555 0 : CALL dbcsr_dot(smatrix(1)%matrix, pmatrix(ispin)%matrix, etr)
556 0 : lri_env%ri_fit%echarge(ispin) = etr
557 : END DO
558 :
559 0 : tvec => lri_env%ri_fit%tvec
560 0 : rm1t => lri_env%ri_fit%rm1t
561 0 : nvec => lri_env%ri_fit%nvec
562 0 : rm1n => lri_env%ri_fit%rm1n
563 :
564 : ! calculate lambda
565 0 : DO ispin = 1, nspin
566 0 : nelec = lri_env%ri_fit%echarge(ispin)
567 0 : nrm1t = SUM(nvec(:)*rm1t(:, ispin))
568 0 : 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 0 : avec => lri_env%ri_fit%avec
573 0 : DO ispin = 1, nspin
574 0 : avec(:, ispin) = rm1t(:, ispin) - 0.5_dp*lri_env%ri_fit%lambda(ispin)*rm1n(:)
575 : END DO
576 :
577 0 : CALL timestop(handle)
578 :
579 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
602 : TYPE(dbcsr_iterator_type) :: iter
603 : TYPE(mp_comm_type) :: group
604 :
605 0 : CALL timeset(routineN, handle)
606 :
607 0 : CALL dbcsr_get_info(mat, matrix_type=matrix_type, group=group)
608 :
609 : SELECT CASE (matrix_type)
610 : CASE (dbcsr_type_no_symmetry)
611 0 : symm = .FALSE.
612 : CASE (dbcsr_type_symmetric)
613 0 : symm = .TRUE.
614 : CASE (dbcsr_type_antisymmetric)
615 0 : CPABORT("NYI, antisymmetric matrix not permitted")
616 : CASE DEFAULT
617 0 : CPABORT("Unknown matrix type, ...")
618 : END SELECT
619 :
620 0 : vo(:) = 0.0_dp
621 0 : CALL dbcsr_iterator_start(iter, mat)
622 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
623 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
624 0 : n1 = ptr(1, iatom)
625 0 : n2 = ptr(2, iatom)
626 0 : nb = n2 - n1 + 1
627 0 : m1 = ptr(1, jatom)
628 0 : m2 = ptr(2, jatom)
629 0 : mb = m2 - m1 + 1
630 0 : CPASSERT(nb == SIZE(block, 1))
631 0 : CPASSERT(mb == SIZE(block, 2))
632 0 : vo(n1:n2) = vo(n1:n2) + MATMUL(block, vi(m1:m2))
633 0 : IF (symm .AND. (iatom /= jatom)) THEN
634 0 : vo(m1:m2) = vo(m1:m2) + MATMUL(TRANSPOSE(block), vi(n1:n2))
635 : END IF
636 : END DO
637 0 : CALL dbcsr_iterator_stop(iter)
638 :
639 0 : CALL group%sum(vo)
640 :
641 0 : CALL timestop(handle)
642 :
643 0 : END SUBROUTINE ri_matvec
644 :
645 0 : END MODULE ri_environment_methods
|
