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_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 :
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 : IF (rerror > threshold) THEN
335 0 : CPWARN("RI solver: CG did not converge properly")
336 : END IF
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 0 : 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 :: group_handle, handle, iatom, jatom, m1, &
599 : m2, mb, n1, 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_handle)
608 0 : CALL group%set_handle(group_handle)
609 :
610 : SELECT CASE (matrix_type)
611 : CASE (dbcsr_type_no_symmetry)
612 0 : symm = .FALSE.
613 : CASE (dbcsr_type_symmetric)
614 0 : symm = .TRUE.
615 : CASE (dbcsr_type_antisymmetric)
616 0 : CPABORT("NYI, antisymmetric matrix not permitted")
617 : CASE DEFAULT
618 0 : CPABORT("Unknown matrix type, ...")
619 : END SELECT
620 :
621 0 : vo(:) = 0.0_dp
622 0 : CALL dbcsr_iterator_start(iter, mat)
623 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
624 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, block)
625 0 : n1 = ptr(1, iatom)
626 0 : n2 = ptr(2, iatom)
627 0 : nb = n2 - n1 + 1
628 0 : m1 = ptr(1, jatom)
629 0 : m2 = ptr(2, jatom)
630 0 : mb = m2 - m1 + 1
631 0 : CPASSERT(nb == SIZE(block, 1))
632 0 : CPASSERT(mb == SIZE(block, 2))
633 0 : vo(n1:n2) = vo(n1:n2) + MATMUL(block, vi(m1:m2))
634 0 : IF (symm .AND. (iatom /= jatom)) THEN
635 0 : vo(m1:m2) = vo(m1:m2) + MATMUL(TRANSPOSE(block), vi(n1:n2))
636 : END IF
637 : END DO
638 0 : CALL dbcsr_iterator_stop(iter)
639 :
640 0 : CALL group%sum(vo)
641 :
642 0 : CALL timestop(handle)
643 :
644 0 : END SUBROUTINE ri_matvec
645 :
646 0 : END MODULE ri_environment_methods
|
