Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates integral matrices for RIGPW method
10 : !> \par History
11 : !> created JGH [08.2012]
12 : !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
13 : !> (2) heavily debugged
14 : !> updated for RI JGH [08.2017]
15 : !> \authors JGH
16 : !> Dorothea Golze
17 : ! **************************************************************************************************
18 : MODULE ri_environment_methods
19 : USE arnoldi_api, ONLY: arnoldi_conjugate_gradient
20 : USE atomic_kind_types, ONLY: atomic_kind_type,&
21 : get_atomic_kind_set
22 : USE cp_blacs_env, ONLY: cp_blacs_env_type
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_create, dbcsr_dot, 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_scale_by_vector, dbcsr_set, dbcsr_type, &
28 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
29 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
30 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
31 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
32 : USE iterate_matrix, ONLY: invert_Hotelling
33 : USE kinds, ONLY: dp
34 : USE lri_environment_types, ONLY: allocate_lri_coefs,&
35 : lri_density_create,&
36 : lri_density_release,&
37 : lri_density_type,&
38 : lri_environment_type,&
39 : lri_kind_type
40 : USE message_passing, ONLY: mp_comm_type,&
41 : mp_para_env_type
42 : USE pw_types, ONLY: pw_c1d_gs_type,&
43 : pw_r3d_rs_type
44 : USE qs_collocate_density, ONLY: calculate_lri_rho_elec
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type,&
47 : set_qs_env
48 : USE qs_ks_types, ONLY: qs_ks_env_type
49 : USE qs_o3c_methods, ONLY: calculate_o3c_integrals,&
50 : contract12_o3c
51 : USE qs_o3c_types, ONLY: get_o3c_iterator_info,&
52 : init_o3c_container,&
53 : o3c_container_type,&
54 : o3c_iterate,&
55 : o3c_iterator_create,&
56 : o3c_iterator_release,&
57 : o3c_iterator_type,&
58 : release_o3c_container
59 : USE qs_overlap, ONLY: build_overlap_matrix
60 : USE qs_rho_types, ONLY: qs_rho_get,&
61 : qs_rho_type
62 : USE util, ONLY: get_limit
63 :
64 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 :
69 : PRIVATE
70 :
71 : ! **************************************************************************************************
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_environment_methods'
74 :
75 : PUBLIC :: build_ri_matrices, calculate_ri_densities, ri_metric_solver
76 :
77 : ! **************************************************************************************************
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief creates and initializes an lri_env
83 : !> \param lri_env the lri_environment you want to create
84 : !> \param qs_env ...
85 : !> \param calculate_forces ...
86 : ! **************************************************************************************************
87 0 : 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 0 : CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces)
94 :
95 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
116 : REAL(KIND=dp) :: fpre
117 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval
118 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout
119 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
120 : TYPE(dbcsr_iterator_type) :: iter
121 0 : 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 0 : CALL timeset(routineN, handle)
130 :
131 0 : CPASSERT(ASSOCIATED(lri_env))
132 0 : CPASSERT(ASSOCIATED(qs_env))
133 :
134 : ! overlap matrices
135 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
136 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
137 0 : 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 0 : IF (calculate_forces) THEN
141 : ! charge constraint forces are calculated here
142 0 : CALL get_qs_env(qs_env=qs_env, rho=rho)
143 0 : CALL qs_rho_get(rho, rho_ao=matrix_p)
144 0 : ALLOCATE (fmat)
145 0 : CALL dbcsr_create(fmat, template=matrix_p(1)%matrix)
146 0 : DO ispin = 1, dft_control%nspins
147 0 : fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n
148 0 : 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 0 : calculate_forces=.TRUE., matrix_p=fmat)
153 0 : CALL dbcsr_release(fmat)
154 0 : DEALLOCATE (fmat)
155 : ELSE
156 : CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, &
157 0 : 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 0 : basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list)
162 0 : IF (calculate_forces) THEN
163 : ! calculate f*a(T) pseudo density matrix for forces
164 0 : bas_ptr => lri_env%ri_fit%bas_ptr
165 0 : avec => lri_env%ri_fit%avec
166 0 : fout => lri_env%ri_fit%fout
167 0 : ALLOCATE (fmat)
168 0 : CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix)
169 0 : CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list)
170 0 : CALL dbcsr_set(fmat, 0.0_dp)
171 0 : CALL dbcsr_iterator_start(iter, fmat)
172 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
173 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk)
174 0 : i1 = bas_ptr(1, iatom)
175 0 : i2 = bas_ptr(2, iatom)
176 0 : j1 = bas_ptr(1, jatom)
177 0 : j2 = bas_ptr(2, jatom)
178 0 : IF (iatom <= jatom) THEN
179 0 : DO ispin = 1, dft_control%nspins
180 0 : DO j = j1, j2
181 0 : m = j - j1 + 1
182 0 : DO i = i1, i2
183 0 : n = i - i1 + 1
184 0 : 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 0 : DO ispin = 1, dft_control%nspins
190 0 : DO i = i1, i2
191 0 : n = i - i1 + 1
192 0 : DO j = j1, j2
193 0 : m = j - j1 + 1
194 0 : 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 0 : IF (iatom == jatom) THEN
200 0 : fblk(:, :) = 0.25_dp*fblk(:, :)
201 : ELSE
202 0 : fblk(:, :) = 0.5_dp*fblk(:, :)
203 : END IF
204 : END DO
205 0 : 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 0 : calculate_forces=.TRUE., matrix_p=fmat)
210 0 : CALL dbcsr_release(fmat)
211 0 : DEALLOCATE (fmat)
212 : END IF
213 : ! approximation (preconditioner) or exact inverse of RI overlap
214 0 : CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1)
215 0 : ALLOCATE (lri_env%ri_sinv(1)%matrix)
216 0 : SELECT CASE (lri_env%ri_sinv_app)
217 : CASE ("NONE")
218 : !
219 : CASE ("INVS")
220 0 : 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 0 : norm_convergence=1.e-10_dp, filter_eps=1.e-12_dp, silent=.FALSE.)
224 : CASE ("INVF")
225 0 : CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix)
226 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
227 0 : CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas)
228 0 : ALLOCATE (eval(nbas))
229 0 : CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env)
230 0 : izero = 0
231 0 : DO i = 1, nbas
232 0 : IF (eval(i) < 1.0e-10_dp) THEN
233 0 : eval(i) = 0.0_dp
234 0 : izero = izero + 1
235 : ELSE
236 0 : eval(i) = SQRT(1.0_dp/eval(i))
237 : END IF
238 : END DO
239 0 : CALL dbcsr_scale_by_vector(emat, eval, side='right')
240 0 : CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix)
241 0 : CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix)
242 0 : DEALLOCATE (eval)
243 0 : CALL dbcsr_release(emat)
244 : CASE ("AINV")
245 0 : 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 0 : norm_convergence=1.e-4_dp, filter_eps=1.e-4_dp, silent=.FALSE.)
249 : CASE DEFAULT
250 0 : 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 0 : ptr=lri_env%ri_fit%bas_ptr)
260 :
261 : ! calculate n(t)x
262 0 : lri_env%ri_fit%ntrm1n = SUM(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:))
263 :
264 : ! calculate 3c-overlap integrals
265 0 : IF (ASSOCIATED(lri_env%o3c)) THEN
266 0 : CALL release_o3c_container(lri_env%o3c)
267 : ELSE
268 0 : 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 0 : lri_env%soo_list, lri_env%soa_list)
273 :
274 0 : CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p)
275 :
276 0 : CALL timestop(handle)
277 :
278 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vect
303 :
304 0 : CALL timeset(routineN, handle)
305 :
306 0 : threshold = 1.e-10_dp
307 0 : max_iter = 100
308 :
309 0 : vecx(:) = vecr(:)
310 :
311 0 : SELECT CASE (solver)
312 : CASE ("NONE")
313 : CALL arnoldi_conjugate_gradient(mat, vecx, &
314 0 : converged=converged, threshold=threshold, max_iter=max_iter)
315 : CASE ("INVS", "INVF")
316 0 : converged = .FALSE.
317 0 : CALL ri_matvec(matp, vecr, vecx, ptr)
318 : CASE ("AINV")
319 : CALL arnoldi_conjugate_gradient(mat, vecx, matrix_p=matp, &
320 0 : converged=converged, threshold=threshold, max_iter=max_iter)
321 : CASE DEFAULT
322 0 : CPABORT("Unknown RI solver")
323 : END SELECT
324 :
325 0 : IF (.NOT. converged) THEN
326 : ! get error
327 0 : rerror = 0.0_dp
328 0 : n = SIZE(vecr)
329 0 : ALLOCATE (vect(n))
330 0 : CALL ri_matvec(mat, vecx, vect, ptr)
331 0 : vect(:) = vect(:) - vecr(:)
332 0 : rerror = MAXVAL(ABS(vect(:)))
333 0 : DEALLOCATE (vect)
334 0 : CPWARN_IF(rerror > threshold, "RI solver: CG did not converge properly")
335 : END IF
336 :
337 0 : CALL timestop(handle)
338 :
339 0 : END SUBROUTINE ri_metric_solver
340 :
341 : ! **************************************************************************************************
342 : !> \brief performs the fitting of the density and distributes the fitted
343 : !> density on the grid
344 : !> \param lri_env the lri environment
345 : !> lri_density the environment for the fitting
346 : !> pmatrix density matrix
347 : !> lri_rho_struct where the fitted density is stored
348 : !> \param qs_env ...
349 : !> \param pmatrix ...
350 : !> \param lri_rho_struct ...
351 : !> \param atomic_kind_set ...
352 : !> \param para_env ...
353 : ! **************************************************************************************************
354 0 : SUBROUTINE calculate_ri_densities(lri_env, qs_env, pmatrix, &
355 : lri_rho_struct, atomic_kind_set, para_env)
356 :
357 : TYPE(lri_environment_type), POINTER :: lri_env
358 : TYPE(qs_environment_type), POINTER :: qs_env
359 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
360 : TYPE(qs_rho_type), INTENT(IN) :: lri_rho_struct
361 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
362 : TYPE(mp_para_env_type), POINTER :: para_env
363 :
364 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_ri_densities'
365 :
366 : INTEGER :: atom_a, handle, i1, i2, iatom, ikind, &
367 : ispin, n, natom, nspin
368 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
369 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
370 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
371 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec
372 : TYPE(lri_density_type), POINTER :: lri_density
373 0 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
374 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
375 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
376 :
377 0 : CALL timeset(routineN, handle)
378 :
379 0 : nspin = SIZE(pmatrix, 1)
380 0 : CALL contract12_o3c(lri_env%o3c, pmatrix)
381 0 : CALL calculate_tvec_ri(lri_env, lri_env%o3c, para_env)
382 0 : CALL calculate_avec_ri(lri_env, pmatrix)
383 : !
384 0 : CALL get_qs_env(qs_env, lri_density=lri_density, natom=natom)
385 0 : IF (ASSOCIATED(lri_density)) THEN
386 0 : CALL lri_density_release(lri_density)
387 : ELSE
388 0 : ALLOCATE (lri_density)
389 : END IF
390 0 : CALL lri_density_create(lri_density)
391 0 : lri_density%nspin = nspin
392 : ! allocate the arrays to hold RI expansion coefficients lri_coefs
393 0 : CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
394 : ! reassign avec
395 0 : avec => lri_env%ri_fit%avec
396 0 : bas_ptr => lri_env%ri_fit%bas_ptr
397 0 : ALLOCATE (atom_of_kind(natom), kind_of(natom))
398 0 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
399 0 : DO ispin = 1, nspin
400 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
401 0 : DO iatom = 1, natom
402 0 : ikind = kind_of(iatom)
403 0 : atom_a = atom_of_kind(iatom)
404 0 : i1 = bas_ptr(1, iatom)
405 0 : i2 = bas_ptr(2, iatom)
406 0 : n = i2 - i1 + 1
407 0 : lri_coef(ikind)%acoef(atom_a, 1:n) = avec(i1:i2, ispin)
408 : END DO
409 : END DO
410 0 : CALL set_qs_env(qs_env, lri_density=lri_density)
411 0 : DEALLOCATE (atom_of_kind, kind_of)
412 : !
413 0 : CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
414 0 : DO ispin = 1, nspin
415 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
416 : CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
417 0 : lri_coef, tot_rho_r(ispin), "RI_HXC", .FALSE.)
418 : END DO
419 :
420 0 : CALL timestop(handle)
421 :
422 0 : END SUBROUTINE calculate_ri_densities
423 :
424 : ! **************************************************************************************************
425 : !> \brief assembles the global vector t=<P.T>
426 : !> \param lri_env the lri environment
427 : !> \param o3c ...
428 : !> \param para_env ...
429 : ! **************************************************************************************************
430 0 : SUBROUTINE calculate_tvec_ri(lri_env, o3c, para_env)
431 :
432 : TYPE(lri_environment_type), POINTER :: lri_env
433 : TYPE(o3c_container_type), POINTER :: o3c
434 : TYPE(mp_para_env_type), POINTER :: para_env
435 :
436 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_tvec_ri'
437 : INTEGER, PARAMETER :: msweep = 32
438 :
439 : INTEGER :: handle, i1, i2, ibl, ibu, il, ispin, &
440 : isweep, it, iu, katom, m, ma, mba, &
441 : mepos, natom, nspin, nsweep, nthread
442 : INTEGER, DIMENSION(2, msweep) :: nlimit
443 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
444 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ta
445 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rm1t, tvec, tvl
446 : TYPE(o3c_iterator_type) :: o3c_iterator
447 :
448 0 : CALL timeset(routineN, handle)
449 :
450 0 : nspin = SIZE(lri_env%ri_fit%tvec, 2)
451 0 : bas_ptr => lri_env%ri_fit%bas_ptr
452 0 : tvec => lri_env%ri_fit%tvec
453 0 : tvec = 0.0_dp
454 0 : natom = SIZE(bas_ptr, 2)
455 : nthread = 1
456 0 : !$ nthread = omp_get_max_threads()
457 :
458 0 : IF (natom < 1000) THEN
459 0 : nsweep = 1
460 : ELSE
461 0 : nsweep = MIN(nthread, msweep)
462 : END IF
463 :
464 0 : nlimit = 0
465 0 : DO isweep = 1, nsweep
466 0 : nlimit(1:2, isweep) = get_limit(natom, nsweep, isweep - 1)
467 : END DO
468 :
469 0 : DO ispin = 1, nspin
470 0 : DO isweep = 1, nsweep
471 0 : il = nlimit(1, isweep)
472 0 : iu = nlimit(2, isweep)
473 0 : ma = iu - il + 1
474 0 : IF (ma < 1) CYCLE
475 0 : ibl = bas_ptr(1, il)
476 0 : ibu = bas_ptr(2, iu)
477 0 : mba = ibu - ibl + 1
478 0 : ALLOCATE (ta(mba, nthread))
479 0 : ta = 0.0_dp
480 :
481 0 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
482 :
483 : !$OMP PARALLEL DEFAULT(NONE)&
484 : !$OMP SHARED (nthread,o3c_iterator,ispin,il,iu,ibl,ta,bas_ptr)&
485 0 : !$OMP PRIVATE (mepos,katom,tvl,i1,i2,m)
486 :
487 : mepos = 0
488 : !$ mepos = omp_get_thread_num()
489 :
490 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
491 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, katom=katom, tvec=tvl)
492 : IF (katom < il .OR. katom > iu) CYCLE
493 : i1 = bas_ptr(1, katom) - ibl + 1
494 : i2 = bas_ptr(2, katom) - ibl + 1
495 : m = i2 - i1 + 1
496 : ta(i1:i2, mepos + 1) = ta(i1:i2, mepos + 1) + tvl(1:m, ispin)
497 : END DO
498 : !$OMP END PARALLEL
499 0 : CALL o3c_iterator_release(o3c_iterator)
500 :
501 : ! sum over threads
502 0 : DO it = 1, nthread
503 0 : tvec(ibl:ibu, ispin) = tvec(ibl:ibu, ispin) + ta(1:mba, it)
504 : END DO
505 0 : DEALLOCATE (ta)
506 : END DO
507 : END DO
508 :
509 : ! global summation
510 0 : CALL para_env%sum(tvec)
511 :
512 0 : rm1t => lri_env%ri_fit%rm1t
513 0 : DO ispin = 1, nspin
514 : ! solve Rx=t
515 : CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
516 : vecr=tvec(:, ispin), &
517 : vecx=rm1t(:, ispin), &
518 : matp=lri_env%ri_sinv(1)%matrix, &
519 : solver=lri_env%ri_sinv_app, &
520 0 : ptr=bas_ptr)
521 : END DO
522 :
523 0 : CALL timestop(handle)
524 :
525 0 : END SUBROUTINE calculate_tvec_ri
526 : ! **************************************************************************************************
527 : !> \brief performs the fitting of the density in the RI method
528 : !> \param lri_env the lri environment
529 : !> lri_density the environment for the fitting
530 : !> pmatrix density matrix
531 : !> \param pmatrix ...
532 : ! **************************************************************************************************
533 0 : SUBROUTINE calculate_avec_ri(lri_env, pmatrix)
534 :
535 : TYPE(lri_environment_type), POINTER :: lri_env
536 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
537 :
538 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_avec_ri'
539 :
540 : INTEGER :: handle, ispin, nspin
541 : REAL(KIND=dp) :: etr, nelec, nrm1t
542 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: nvec, rm1n
543 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, rm1t, tvec
544 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: smatrix
545 :
546 0 : CALL timeset(routineN, handle)
547 :
548 0 : nspin = SIZE(pmatrix)
549 : ! number of electrons
550 0 : smatrix => lri_env%ob_smat
551 0 : DO ispin = 1, nspin
552 0 : etr = 0.0_dp
553 0 : CALL dbcsr_dot(smatrix(1)%matrix, pmatrix(ispin)%matrix, etr)
554 0 : lri_env%ri_fit%echarge(ispin) = etr
555 : END DO
556 :
557 0 : tvec => lri_env%ri_fit%tvec
558 0 : rm1t => lri_env%ri_fit%rm1t
559 0 : nvec => lri_env%ri_fit%nvec
560 0 : rm1n => lri_env%ri_fit%rm1n
561 :
562 : ! calculate lambda
563 0 : DO ispin = 1, nspin
564 0 : nelec = lri_env%ri_fit%echarge(ispin)
565 0 : nrm1t = SUM(nvec(:)*rm1t(:, ispin))
566 0 : lri_env%ri_fit%lambda(ispin) = 2.0_dp*(nrm1t - nelec)/lri_env%ri_fit%ntrm1n
567 : END DO
568 :
569 : ! calculate avec = rm1t - lambda/2 * rm1n
570 0 : avec => lri_env%ri_fit%avec
571 0 : DO ispin = 1, nspin
572 0 : avec(:, ispin) = rm1t(:, ispin) - 0.5_dp*lri_env%ri_fit%lambda(ispin)*rm1n(:)
573 : END DO
574 :
575 0 : CALL timestop(handle)
576 :
577 0 : END SUBROUTINE calculate_avec_ri
578 :
579 : ! **************************************************************************************************
580 : !> \brief Mutiplies a replicated vector with a DBCSR matrix: vo = mat*vi
581 : !> \param mat ...
582 : !> \param vi ...
583 : !> \param vo ...
584 : !> \param ptr ...
585 : ! **************************************************************************************************
586 0 : SUBROUTINE ri_matvec(mat, vi, vo, ptr)
587 :
588 : TYPE(dbcsr_type) :: mat
589 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vi
590 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vo
591 : INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr
592 :
593 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_matvec'
594 :
595 : CHARACTER :: matrix_type
596 : INTEGER :: group_handle, handle, iatom, jatom, m1, &
597 : m2, mb, n1, n2, nb
598 : LOGICAL :: symm
599 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
600 : TYPE(dbcsr_iterator_type) :: iter
601 : TYPE(mp_comm_type) :: group
602 :
603 0 : CALL timeset(routineN, handle)
604 :
605 0 : CALL dbcsr_get_info(mat, matrix_type=matrix_type, group=group_handle)
606 0 : CALL group%set_handle(group_handle)
607 :
608 : SELECT CASE (matrix_type)
609 : CASE (dbcsr_type_no_symmetry)
610 0 : symm = .FALSE.
611 : CASE (dbcsr_type_symmetric)
612 0 : symm = .TRUE.
613 : CASE (dbcsr_type_antisymmetric)
614 0 : CPABORT("NYI, antisymmetric matrix not permitted")
615 : CASE DEFAULT
616 0 : CPABORT("Unknown matrix type, ...")
617 : END SELECT
618 :
619 0 : vo(:) = 0.0_dp
620 0 : CALL dbcsr_iterator_start(iter, mat)
621 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
622 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
623 0 : n1 = ptr(1, iatom)
624 0 : n2 = ptr(2, iatom)
625 0 : nb = n2 - n1 + 1
626 0 : m1 = ptr(1, jatom)
627 0 : m2 = ptr(2, jatom)
628 0 : mb = m2 - m1 + 1
629 0 : CPASSERT(nb == SIZE(block, 1))
630 0 : CPASSERT(mb == SIZE(block, 2))
631 0 : vo(n1:n2) = vo(n1:n2) + MATMUL(block, vi(m1:m2))
632 0 : IF (symm .AND. (iatom /= jatom)) THEN
633 0 : vo(m1:m2) = vo(m1:m2) + MATMUL(TRANSPOSE(block), vi(n1:n2))
634 : END IF
635 : END DO
636 0 : CALL dbcsr_iterator_stop(iter)
637 :
638 0 : CALL group%sum(vo)
639 :
640 0 : CALL timestop(handle)
641 :
642 0 : END SUBROUTINE ri_matvec
643 :
644 0 : END MODULE ri_environment_methods
|
