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 LRIGPW method
10 : !> lri : local resolution of the identity
11 : !> \par History
12 : !> created JGH [08.2012]
13 : !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
14 : !> (2) heavily debugged
15 : !> \authors JGH
16 : !> Dorothea Golze
17 : ! **************************************************************************************************
18 : MODULE lri_environment_methods
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind,&
21 : get_atomic_kind_set
22 : USE basis_set_types, ONLY: gto_basis_set_type
23 : USE cell_types, ONLY: cell_type
24 : USE core_ppl, ONLY: build_core_ppl_ri
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
27 : dbcsr_get_block_p,&
28 : dbcsr_p_type,&
29 : dbcsr_release,&
30 : dbcsr_replicate_all,&
31 : dbcsr_type
32 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
33 : dbcsr_get_block_diag
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_type
36 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
37 : cp_print_key_unit_nr
38 : USE generic_os_integrals, ONLY: int_overlap_aabb_os
39 : USE input_constants, ONLY: do_lri_inv,&
40 : do_lri_inv_auto,&
41 : do_lri_pseudoinv_diag,&
42 : do_lri_pseudoinv_svd
43 : USE input_section_types, ONLY: section_vals_type
44 : USE kinds, ONLY: dp
45 : USE lri_compression, ONLY: lri_comp,&
46 : lri_cont_mem,&
47 : lri_decomp_i
48 : USE lri_environment_types, ONLY: &
49 : allocate_lri_coefs, allocate_lri_ints, allocate_lri_ints_rho, allocate_lri_ppl_ints, &
50 : allocate_lri_rhos, deallocate_lri_ints, deallocate_lri_ints_rho, deallocate_lri_ppl_ints, &
51 : lri_density_create, lri_density_release, lri_density_type, lri_environment_type, &
52 : lri_int_rho_type, lri_int_type, lri_kind_type, lri_list_type, lri_rhoab_type
53 : USE lri_integrals, ONLY: allocate_int_type,&
54 : deallocate_int_type,&
55 : int_type,&
56 : lri_int2
57 : USE mathlib, ONLY: get_pseudo_inverse_diag,&
58 : get_pseudo_inverse_svd,&
59 : invmat_symm
60 : USE message_passing, ONLY: mp_para_env_type
61 : USE particle_types, ONLY: particle_type
62 : USE pw_types, ONLY: pw_c1d_gs_type,&
63 : pw_r3d_rs_type
64 : USE qs_collocate_density, ONLY: calculate_lri_rho_elec
65 : USE qs_environment_types, ONLY: get_qs_env,&
66 : qs_environment_type
67 : USE qs_force_types, ONLY: qs_force_type
68 : USE qs_kind_types, ONLY: qs_kind_type
69 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
70 : neighbor_list_iterate,&
71 : neighbor_list_iterator_create,&
72 : neighbor_list_iterator_p_type,&
73 : neighbor_list_iterator_release,&
74 : neighbor_list_set_p_type
75 : USE qs_rho_types, ONLY: qs_rho_get,&
76 : qs_rho_set,&
77 : qs_rho_type
78 : USE virial_types, ONLY: virial_type
79 :
80 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : ! **************************************************************************************************
88 :
89 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_environment_methods'
90 :
91 : PUBLIC :: build_lri_matrices, calculate_lri_densities, &
92 : calculate_lri_integrals, calculate_avec_lri, &
93 : v_int_ppl_update, v_int_ppl_energy, lri_kg_rho_update, lri_print_stat
94 :
95 : ! **************************************************************************************************
96 :
97 : CONTAINS
98 :
99 : ! **************************************************************************************************
100 : !> \brief creates and initializes an lri_env
101 : !> \param lri_env the lri_environment you want to create
102 : !> \param qs_env ...
103 : ! **************************************************************************************************
104 68 : SUBROUTINE build_lri_matrices(lri_env, qs_env)
105 :
106 : TYPE(lri_environment_type), POINTER :: lri_env
107 : TYPE(qs_environment_type), POINTER :: qs_env
108 :
109 : ! calculate the integrals needed to do the local (2-center) expansion
110 : ! of the (pair) densities
111 68 : CALL calculate_lri_integrals(lri_env, qs_env)
112 :
113 : ! calculate integrals for local pp (if used in LRI)
114 68 : IF (lri_env%ppl_ri) THEN
115 2 : CALL calculate_lri_ppl_integrals(lri_env, qs_env, .FALSE.)
116 : END IF
117 :
118 68 : END SUBROUTINE build_lri_matrices
119 :
120 : ! **************************************************************************************************
121 : !> \brief calculates integrals needed for the LRI density fitting,
122 : !> integrals are calculated once, before the SCF starts
123 : !> \param lri_env ...
124 : !> \param qs_env ...
125 : ! **************************************************************************************************
126 88 : SUBROUTINE calculate_lri_integrals(lri_env, qs_env)
127 :
128 : TYPE(lri_environment_type), POINTER :: lri_env
129 : TYPE(qs_environment_type), POINTER :: qs_env
130 :
131 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_integrals'
132 :
133 : INTEGER :: handle, iac, iatom, ikind, ilist, jatom, &
134 : jkind, jneighbor, mepos, nba, nbb, &
135 : nfa, nfb, nkind, nlist, nn, nneighbor, &
136 : nthread
137 : LOGICAL :: e1c, use_virial
138 : REAL(KIND=dp) :: cmem, cpff, cpsr, cptt, dab
139 : REAL(KIND=dp), DIMENSION(3) :: rab
140 : TYPE(cell_type), POINTER :: cell
141 : TYPE(dft_control_type), POINTER :: dft_control
142 : TYPE(gto_basis_set_type), POINTER :: fbasa, fbasb, obasa, obasb
143 88 : TYPE(int_type) :: lriint
144 : TYPE(lri_int_type), POINTER :: lrii
145 : TYPE(lri_list_type), POINTER :: lri_ints
146 : TYPE(neighbor_list_iterator_p_type), &
147 88 : DIMENSION(:), POINTER :: nl_iterator
148 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 88 : POINTER :: soo_list
150 88 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 : TYPE(virial_type), POINTER :: virial
152 :
153 88 : CALL timeset(routineN, handle)
154 88 : NULLIFY (cell, dft_control, fbasa, fbasb, lrii, lri_ints, nl_iterator, &
155 88 : obasa, obasb, particle_set, soo_list, virial)
156 :
157 88 : lri_env%stat%pairs_tt = 0.0_dp
158 88 : lri_env%stat%pairs_sr = 0.0_dp
159 88 : lri_env%stat%pairs_ff = 0.0_dp
160 88 : lri_env%stat%overlap_error = 0.0_dp
161 88 : lri_env%stat%abai_mem = 0.0_dp
162 :
163 88 : IF (ASSOCIATED(lri_env%soo_list)) THEN
164 88 : soo_list => lri_env%soo_list
165 :
166 : CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, &
167 88 : nkind=nkind, particle_set=particle_set, virial=virial)
168 88 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
169 : nthread = 1
170 88 : !$ nthread = omp_get_max_threads()
171 :
172 88 : IF (ASSOCIATED(lri_env%lri_ints)) THEN
173 30 : CALL deallocate_lri_ints(lri_env%lri_ints)
174 : END IF
175 :
176 : ! allocate matrices storing the LRI integrals
177 88 : CALL allocate_lri_ints(lri_env, lri_env%lri_ints, nkind)
178 88 : lri_ints => lri_env%lri_ints
179 :
180 88 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
181 : !$OMP PARALLEL DEFAULT(NONE)&
182 : !$OMP SHARED (nthread,nl_iterator,lri_env,lri_ints,nkind,use_virial)&
183 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,rab,dab,&
184 88 : !$OMP e1c,iac,obasa,obasb,fbasa,fbasb,lrii,lriint,nba,nbb,nfa,nfb,nn,cptt,cpsr,cpff,cmem)
185 :
186 : mepos = 0
187 : !$ mepos = omp_get_thread_num()
188 :
189 : cptt = 0.0_dp
190 : cpsr = 0.0_dp
191 : cpff = 0.0_dp
192 : cmem = 0.0_dp
193 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
194 :
195 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
196 : nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
197 : iatom=iatom, jatom=jatom, r=rab)
198 : iac = ikind + nkind*(jkind - 1)
199 : dab = SQRT(SUM(rab*rab))
200 :
201 : obasa => lri_env%orb_basis(ikind)%gto_basis_set
202 : obasb => lri_env%orb_basis(jkind)%gto_basis_set
203 : fbasa => lri_env%ri_basis(ikind)%gto_basis_set
204 : fbasb => lri_env%ri_basis(jkind)%gto_basis_set
205 :
206 : IF (.NOT. ASSOCIATED(obasa)) CYCLE
207 : IF (.NOT. ASSOCIATED(obasb)) CYCLE
208 :
209 : lrii => lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
210 :
211 : ! exact 1 center approximation
212 : e1c = .FALSE.
213 : IF (iatom == jatom .AND. dab < lri_env%delta) e1c = .TRUE.
214 : ! forces: not every pair is giving contributions
215 : ! no forces for self-pair aa
216 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
217 : lrii%calc_force_pair = .FALSE.
218 : ELSE
219 : ! forces for periodic self-pair aa' required for virial
220 : IF (iatom == jatom .AND. .NOT. use_virial) THEN
221 : lrii%calc_force_pair = .FALSE.
222 : ELSE
223 : lrii%calc_force_pair = .TRUE.
224 : END IF
225 : END IF
226 :
227 : IF (e1c .AND. lri_env%exact_1c_terms) THEN
228 : ! nothing to do
229 : ELSE
230 : cptt = cptt + 1.0_dp
231 :
232 : nba = obasa%nsgf
233 : nbb = obasb%nsgf
234 : nfa = fbasa%nsgf
235 : nfb = fbasb%nsgf
236 :
237 : lrii%nba = nba
238 : lrii%nbb = nbb
239 : lrii%nfa = nfa
240 : lrii%nfb = nfb
241 :
242 : ! full method is used
243 : ! calculate integrals (fa,fb), (a,b), (a,b,fa) and (a,b,fb)
244 : CALL allocate_int_type(lriint=lriint, &
245 : nba=nba, nbb=nbb, nfa=nfa, nfb=nfb, skip_sab=(.NOT. lrii%lrisr))
246 : CALL lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, &
247 : iatom, jatom, ikind, jkind)
248 : ! store abaint/abbint in compressed form
249 : IF (e1c) THEN
250 : CALL lri_comp(lriint%abaint, lrii%abascr, lrii%cabai)
251 : cmem = cmem + lri_cont_mem(lrii%cabai)
252 : ELSE
253 : CALL lri_comp(lriint%abaint, lrii%abascr, lrii%cabai)
254 : cmem = cmem + lri_cont_mem(lrii%cabai)
255 : CALL lri_comp(lriint%abbint, lrii%abbscr, lrii%cabbi)
256 : cmem = cmem + lri_cont_mem(lrii%cabbi)
257 : END IF
258 : ! store overlap
259 : lrii%soo(1:nba, 1:nbb) = lriint%sooint(1:nba, 1:nbb)
260 :
261 : ! Full LRI method
262 : IF (lrii%lrisr) THEN
263 : cpsr = cpsr + 1.0_dp
264 : ! construct and invert S matrix
265 : ! calculate Sinv*n and n*Sinv*n
266 : IF (e1c) THEN
267 : lrii%sinv(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp_inv(1:nfa, 1:nfa)
268 : lrii%n(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
269 : CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%sinv(1, 1), nfa, &
270 : lrii%n(1), 1, 0.0_dp, lrii%sn, 1)
271 : lrii%nsn = SUM(lrii%sn(1:nfa)*lrii%n(1:nfa))
272 : ELSE
273 : nn = nfa + nfb
274 : CALL inverse_lri_overlap(lri_env, lrii%sinv, lri_env%bas_prop(ikind)%ri_ovlp, &
275 : lri_env%bas_prop(jkind)%ri_ovlp, lriint%sabint)
276 : lrii%n(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
277 : lrii%n(nfa + 1:nn) = lri_env%bas_prop(jkind)%int_fbas(1:nfb)
278 : CALL dgemv("N", nn, nn, 1.0_dp, lrii%sinv(1, 1), nn, &
279 : lrii%n(1), 1, 0.0_dp, lrii%sn, 1)
280 : lrii%nsn = SUM(lrii%sn(1:nn)*lrii%n(1:nn))
281 : END IF
282 : IF (lri_env%store_integrals) THEN
283 : IF (ALLOCATED(lrii%sab)) DEALLOCATE (lrii%sab)
284 : ALLOCATE (lrii%sab(nfa, nfb))
285 : lrii%sab(1:nfa, 1:nfb) = lriint%sabint(1:nfa, 1:nfb)
286 : END IF
287 : END IF
288 :
289 : ! Distant Pair methods
290 : IF (lrii%lriff) THEN
291 : cpff = cpff + 1.0_dp
292 : CPASSERT(.NOT. e1c)
293 : CPASSERT(.NOT. lri_env%store_integrals)
294 : ! calculate Sinv*n and n*Sinv*n for A and B centers
295 : lrii%na(1:nfa) = lri_env%bas_prop(ikind)%int_fbas(1:nfa)
296 : lrii%nb(1:nfb) = lri_env%bas_prop(jkind)%int_fbas(1:nfb)
297 : CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%asinv(1, 1), nfa, &
298 : lrii%na(1), 1, 0.0_dp, lrii%sna, 1)
299 : lrii%nsna = SUM(lrii%sna(1:nfa)*lrii%na(1:nfa))
300 : CALL dgemv("N", nfb, nfb, 1.0_dp, lrii%bsinv(1, 1), nfb, &
301 : lrii%nb(1), 1, 0.0_dp, lrii%snb, 1)
302 : lrii%nsnb = SUM(lrii%snb(1:nfb)*lrii%nb(1:nfb))
303 : END IF
304 :
305 : CALL deallocate_int_type(lriint=lriint)
306 :
307 : END IF
308 :
309 : END DO
310 : !$OMP CRITICAL(UPDATE)
311 : lri_env%stat%pairs_tt = lri_env%stat%pairs_tt + cptt
312 : lri_env%stat%pairs_sr = lri_env%stat%pairs_sr + cpsr
313 : lri_env%stat%pairs_ff = lri_env%stat%pairs_ff + cpff
314 : lri_env%stat%abai_mem = lri_env%stat%abai_mem + cmem
315 : !$OMP END CRITICAL(UPDATE)
316 :
317 : !$OMP END PARALLEL
318 :
319 88 : CALL neighbor_list_iterator_release(nl_iterator)
320 :
321 88 : IF (lri_env%debug) THEN
322 2 : CALL output_debug_info(lri_env, qs_env, lri_ints, soo_list)
323 : END IF
324 :
325 : END IF
326 :
327 88 : CALL timestop(handle)
328 :
329 176 : END SUBROUTINE calculate_lri_integrals
330 :
331 : ! **************************************************************************************************
332 : !> \brief ...
333 : !> \param rho_struct ...
334 : !> \param qs_env ...
335 : !> \param lri_env ...
336 : !> \param lri_density ...
337 : !> \param atomlist ...
338 : ! **************************************************************************************************
339 100 : SUBROUTINE lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
340 :
341 : TYPE(qs_rho_type), POINTER :: rho_struct
342 : TYPE(qs_environment_type), POINTER :: qs_env
343 : TYPE(lri_environment_type), POINTER :: lri_env
344 : TYPE(lri_density_type), POINTER :: lri_density
345 : INTEGER, DIMENSION(:), INTENT(IN) :: atomlist
346 :
347 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_kg_rho_update'
348 :
349 : INTEGER :: handle, ispin, nspins
350 100 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
351 : TYPE(dbcsr_type) :: pmat_diag
352 : TYPE(dft_control_type), POINTER :: dft_control
353 : TYPE(mp_para_env_type), POINTER :: para_env
354 100 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
355 100 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
356 :
357 100 : CALL timeset(routineN, handle)
358 :
359 100 : CPASSERT(ASSOCIATED(rho_struct))
360 :
361 100 : CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
362 :
363 100 : CALL qs_rho_get(rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
364 :
365 100 : nspins = dft_control%nspins
366 :
367 100 : CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
368 100 : CPASSERT(.NOT. lri_env%exact_1c_terms)
369 :
370 200 : DO ispin = 1, nspins
371 : CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, &
372 : lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
373 200 : "LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag, atomlist=atomlist)
374 : END DO
375 :
376 100 : CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
377 :
378 100 : CALL timestop(handle)
379 :
380 100 : END SUBROUTINE lri_kg_rho_update
381 :
382 : ! **************************************************************************************************
383 : !> \brief calculates integrals needed for the LRI ppl method,
384 : !> integrals are calculated once, before the SCF starts
385 : !> For forces we assume integrals are available and will not be updated
386 : !> \param lri_env ...
387 : !> \param qs_env ...
388 : !> \param calculate_forces ...
389 : ! **************************************************************************************************
390 4 : SUBROUTINE calculate_lri_ppl_integrals(lri_env, qs_env, calculate_forces)
391 :
392 : TYPE(lri_environment_type), POINTER :: lri_env
393 : TYPE(qs_environment_type), POINTER :: qs_env
394 : LOGICAL, INTENT(IN) :: calculate_forces
395 :
396 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_ppl_integrals'
397 :
398 : INTEGER :: handle, ikind, ispin, na, nb, nkind, &
399 : nspin
400 : LOGICAL :: use_virial
401 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
402 : TYPE(lri_density_type), POINTER :: lri_density
403 4 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_ppl_coef
404 : TYPE(mp_para_env_type), POINTER :: para_env
405 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
406 4 : POINTER :: sac_lri
407 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
408 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
409 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
410 : TYPE(virial_type), POINTER :: virial
411 :
412 4 : CALL timeset(routineN, handle)
413 :
414 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, &
415 4 : particle_set=particle_set, atomic_kind_set=atomic_kind_set)
416 4 : IF (calculate_forces) THEN
417 2 : CPASSERT(ASSOCIATED(lri_env%lri_ppl_ints))
418 2 : CALL get_qs_env(qs_env, lri_density=lri_density)
419 2 : nspin = SIZE(lri_density%lri_coefs, 1)
420 : ELSE
421 2 : IF (ASSOCIATED(lri_env%lri_ppl_ints)) THEN
422 0 : CALL deallocate_lri_ppl_ints(lri_env%lri_ppl_ints)
423 : END IF
424 2 : CALL allocate_lri_ppl_ints(lri_env, lri_env%lri_ppl_ints, atomic_kind_set)
425 : END IF
426 : !
427 4 : CALL get_qs_env(qs_env, sac_lri=sac_lri, force=force, virial=virial, para_env=para_env)
428 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
429 4 : use_virial = use_virial .AND. calculate_forces
430 : !
431 4 : CALL get_qs_env(qs_env, nkind=nkind)
432 20 : ALLOCATE (lri_ppl_coef(nkind))
433 12 : DO ikind = 1, nkind
434 8 : na = SIZE(lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int, 1)
435 8 : nb = SIZE(lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int, 2)
436 8 : NULLIFY (lri_ppl_coef(ikind)%acoef)
437 8 : NULLIFY (lri_ppl_coef(ikind)%v_int)
438 8 : NULLIFY (lri_ppl_coef(ikind)%v_dadr)
439 8 : NULLIFY (lri_ppl_coef(ikind)%v_dfdr)
440 32 : ALLOCATE (lri_ppl_coef(ikind)%v_int(na, nb))
441 2140 : lri_ppl_coef(ikind)%v_int = 0.0_dp
442 12 : IF (calculate_forces) THEN
443 12 : ALLOCATE (lri_ppl_coef(ikind)%acoef(na, nb))
444 1070 : lri_ppl_coef(ikind)%acoef = 0.0_dp
445 8 : DO ispin = 1, nspin
446 : lri_ppl_coef(ikind)%acoef(1:na, 1:nb) = lri_ppl_coef(ikind)%acoef(1:na, 1:nb) + &
447 2144 : lri_density%lri_coefs(ispin)%lri_kinds(ikind)%acoef(1:na, 1:nb)
448 : END DO
449 : END IF
450 : END DO
451 : !
452 : CALL build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, &
453 : qs_kind_set, atomic_kind_set, particle_set, sac_lri, &
454 4 : "LRI_AUX")
455 : !
456 4 : IF (.NOT. calculate_forces) THEN
457 6 : DO ikind = 1, nkind
458 2136 : CALL para_env%sum(lri_ppl_coef(ikind)%v_int)
459 2142 : lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int = lri_ppl_coef(ikind)%v_int
460 : END DO
461 : END IF
462 : !
463 12 : DO ikind = 1, nkind
464 8 : IF (ASSOCIATED(lri_ppl_coef(ikind)%acoef)) DEALLOCATE (lri_ppl_coef(ikind)%acoef)
465 8 : IF (ASSOCIATED(lri_ppl_coef(ikind)%v_int)) DEALLOCATE (lri_ppl_coef(ikind)%v_int)
466 8 : IF (ASSOCIATED(lri_ppl_coef(ikind)%v_dadr)) DEALLOCATE (lri_ppl_coef(ikind)%v_dadr)
467 12 : IF (ASSOCIATED(lri_ppl_coef(ikind)%v_dfdr)) DEALLOCATE (lri_ppl_coef(ikind)%v_dfdr)
468 : END DO
469 4 : DEALLOCATE (lri_ppl_coef)
470 :
471 4 : CALL timestop(handle)
472 :
473 4 : END SUBROUTINE calculate_lri_ppl_integrals
474 :
475 : ! **************************************************************************************************
476 : !> \brief calculates overlap integrals (aabb) of the orbital basis set,
477 : !> required for LRI basis set optimization
478 : !> \param lri_env ...
479 : !> \param qs_env ...
480 : ! **************************************************************************************************
481 0 : SUBROUTINE calculate_lri_overlap_aabb(lri_env, qs_env)
482 :
483 : TYPE(lri_environment_type), POINTER :: lri_env
484 : TYPE(qs_environment_type), POINTER :: qs_env
485 :
486 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_overlap_aabb'
487 :
488 : INTEGER :: handle, iac, iatom, ikind, ilist, jatom, &
489 : jkind, jneighbor, nba, nbb, nkind, &
490 : nlist, nneighbor
491 : REAL(KIND=dp) :: dab
492 : REAL(KIND=dp), DIMENSION(3) :: rab
493 : TYPE(cell_type), POINTER :: cell
494 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb
495 : TYPE(lri_int_rho_type), POINTER :: lriir
496 : TYPE(lri_list_type), POINTER :: lri_ints_rho
497 : TYPE(neighbor_list_iterator_p_type), &
498 0 : DIMENSION(:), POINTER :: nl_iterator
499 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
500 0 : POINTER :: soo_list
501 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
502 :
503 0 : CALL timeset(routineN, handle)
504 0 : NULLIFY (cell, lriir, lri_ints_rho, nl_iterator, obasa, obasb, &
505 0 : particle_set, soo_list)
506 :
507 0 : IF (ASSOCIATED(lri_env%soo_list)) THEN
508 0 : soo_list => lri_env%soo_list
509 :
510 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, particle_set=particle_set, &
511 0 : cell=cell)
512 :
513 0 : IF (ASSOCIATED(lri_env%lri_ints_rho)) THEN
514 0 : CALL deallocate_lri_ints_rho(lri_env%lri_ints_rho)
515 : END IF
516 :
517 0 : CALL allocate_lri_ints_rho(lri_env, lri_env%lri_ints_rho, nkind)
518 0 : lri_ints_rho => lri_env%lri_ints_rho
519 :
520 0 : CALL neighbor_list_iterator_create(nl_iterator, soo_list)
521 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
522 :
523 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
524 : nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
525 0 : iatom=iatom, jatom=jatom, r=rab)
526 :
527 0 : iac = ikind + nkind*(jkind - 1)
528 0 : dab = SQRT(SUM(rab*rab))
529 :
530 0 : obasa => lri_env%orb_basis(ikind)%gto_basis_set
531 0 : obasb => lri_env%orb_basis(jkind)%gto_basis_set
532 0 : IF (.NOT. ASSOCIATED(obasa)) CYCLE
533 0 : IF (.NOT. ASSOCIATED(obasb)) CYCLE
534 :
535 0 : lriir => lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
536 :
537 0 : nba = obasa%nsgf
538 0 : nbb = obasb%nsgf
539 :
540 : ! calculate integrals (aa,bb)
541 : CALL int_overlap_aabb_os(lriir%soaabb, obasa, obasb, rab, lri_env%debug, &
542 0 : lriir%dmax_aabb)
543 :
544 : END DO
545 :
546 0 : CALL neighbor_list_iterator_release(nl_iterator)
547 :
548 : END IF
549 :
550 0 : CALL timestop(handle)
551 :
552 0 : END SUBROUTINE calculate_lri_overlap_aabb
553 :
554 : ! **************************************************************************************************
555 : !> \brief performs the fitting of the density and distributes the fitted
556 : !> density on the grid
557 : !> \param lri_env the lri environment
558 : !> \param lri_density ...
559 : !> \param qs_env ...
560 : !> \param pmatrix ...
561 : !> \param cell_to_index ...
562 : !> \param lri_rho_struct ...
563 : !> \param atomic_kind_set ...
564 : !> \param para_env ...
565 : !> \param response_density ...
566 : ! **************************************************************************************************
567 792 : SUBROUTINE calculate_lri_densities(lri_env, lri_density, qs_env, pmatrix, cell_to_index, &
568 : lri_rho_struct, atomic_kind_set, para_env, response_density)
569 :
570 : TYPE(lri_environment_type), POINTER :: lri_env
571 : TYPE(lri_density_type), POINTER :: lri_density
572 : TYPE(qs_environment_type), POINTER :: qs_env
573 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
574 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
575 : TYPE(qs_rho_type), INTENT(INOUT) :: lri_rho_struct
576 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
577 : TYPE(mp_para_env_type), POINTER :: para_env
578 : LOGICAL, INTENT(IN) :: response_density
579 :
580 792 : CALL calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index, response_density)
581 : CALL distribute_lri_density_on_the_grid(lri_env, lri_density, qs_env, &
582 : lri_rho_struct, atomic_kind_set, para_env, &
583 792 : response_density)
584 :
585 792 : END SUBROUTINE calculate_lri_densities
586 :
587 : ! **************************************************************************************************
588 : !> \brief performs the fitting of the density; solves the linear system of
589 : !> equations; yield the expansion coefficients avec
590 : !> \param lri_env the lri environment
591 : !> lri_density the environment for the fitting
592 : !> pmatrix density matrix
593 : !> \param lri_density ...
594 : !> \param pmatrix ...
595 : !> \param cell_to_index ...
596 : !> \param response_density ...
597 : ! **************************************************************************************************
598 810 : SUBROUTINE calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index, response_density)
599 :
600 : TYPE(lri_environment_type), POINTER :: lri_env
601 : TYPE(lri_density_type), POINTER :: lri_density
602 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
603 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
604 : LOGICAL, INTENT(IN), OPTIONAL :: response_density
605 :
606 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_avec_lri'
607 :
608 : INTEGER :: handle, i, iac, iatom, ic, ikind, ilist, ispin, jatom, jkind, jneighbor, mepos, &
609 : nba, nbb, nfa, nfb, nkind, nlist, nn, nneighbor, nspin, nthread
610 : INTEGER, DIMENSION(3) :: cell
611 : LOGICAL :: found, lresponse, trans, use_cell_mapping
612 : REAL(KIND=dp) :: dab, rab(3), threshold
613 810 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: m
614 810 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: int3, paa, pab, pbb
615 810 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pbij
616 : TYPE(dbcsr_type), POINTER :: pmat
617 : TYPE(lri_int_type), POINTER :: lrii
618 : TYPE(lri_list_type), POINTER :: lri_rho
619 : TYPE(lri_rhoab_type), POINTER :: lrho
620 : TYPE(neighbor_list_iterator_p_type), &
621 810 : DIMENSION(:), POINTER :: nl_iterator
622 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
623 810 : POINTER :: soo_list
624 :
625 810 : CALL timeset(routineN, handle)
626 810 : NULLIFY (lrii, lri_rho, nl_iterator, pbij, pmat, soo_list)
627 :
628 810 : lresponse = .FALSE.
629 810 : IF (PRESENT(response_density)) lresponse = response_density
630 :
631 810 : IF (ASSOCIATED(lri_env%soo_list)) THEN
632 810 : soo_list => lri_env%soo_list
633 :
634 810 : nspin = SIZE(pmatrix, 1)
635 810 : use_cell_mapping = (SIZE(pmatrix, 2) > 1)
636 810 : nkind = lri_env%lri_ints%nkind
637 : nthread = 1
638 810 : !$ nthread = omp_get_max_threads()
639 :
640 810 : IF (ASSOCIATED(lri_density)) THEN
641 752 : CALL lri_density_release(lri_density)
642 : ELSE
643 58 : ALLOCATE (lri_density)
644 : END IF
645 810 : CALL lri_density_create(lri_density)
646 810 : lri_density%nspin = nspin
647 :
648 : ! allocate structure lri_rhos and vectors tvec and avec
649 810 : CALL allocate_lri_rhos(lri_env, lri_density%lri_rhos, nspin, nkind)
650 :
651 1654 : DO ispin = 1, nspin
652 844 : lri_rho => lri_density%lri_rhos(ispin)%lri_list
653 :
654 844 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
655 : !$OMP PARALLEL DEFAULT(NONE)&
656 : !$OMP SHARED (nthread,nl_iterator,lri_env,lri_rho,pmatrix,nkind,cell_to_index,use_cell_mapping,ispin,&
657 : !$OMP lresponse)&
658 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,rab,dab,iac,&
659 844 : !$OMP trans,found,pmat,pbij,pab,paa,pbb,int3,lrho,lrii,nba,nbb,nfa,nfb,nn,threshold,i,m,cell,ic)
660 :
661 : mepos = 0
662 : !$ mepos = omp_get_thread_num()
663 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
664 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
665 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
666 : r=rab, cell=cell)
667 :
668 : iac = ikind + nkind*(jkind - 1)
669 : dab = SQRT(SUM(rab*rab))
670 :
671 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
672 : IF (iatom == jatom .AND. dab < lri_env%delta .AND. lri_env%exact_1c_terms) CYCLE
673 :
674 : IF (use_cell_mapping) THEN
675 : ic = cell_to_index(cell(1), cell(2), cell(3))
676 : CPASSERT(ic > 0)
677 : ELSE
678 : ic = 1
679 : END IF
680 : pmat => pmatrix(ispin, ic)%matrix
681 : ! get the density matrix Pab
682 : ! this needs to be response density for LRI-TDDFT
683 : NULLIFY (pbij)
684 : IF (iatom <= jatom) THEN
685 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pbij, found=found)
686 : trans = .FALSE.
687 : ELSE
688 : CALL dbcsr_get_block_p(matrix=pmat, row=jatom, col=iatom, block=pbij, found=found)
689 : trans = .TRUE.
690 : END IF
691 : CPASSERT(found)
692 :
693 : lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
694 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
695 :
696 : nba = lrii%nba
697 : nbb = lrii%nbb
698 : nfa = lrii%nfa
699 : nfb = lrii%nfb
700 :
701 : nn = nfa + nfb
702 :
703 : ! compute tvec = SUM_ab Pab *(a,b,x) and charge constraint
704 : ALLOCATE (pab(nba, nbb))
705 : IF (trans) THEN
706 : pab(1:nba, 1:nbb) = TRANSPOSE(pbij(1:nbb, 1:nba))
707 : ELSE
708 : pab(1:nba, 1:nbb) = pbij(1:nba, 1:nbb)
709 : END IF
710 :
711 : ALLOCATE (int3(nba, nbb))
712 : IF (lrii%lrisr) THEN
713 : ! full LRI method
714 : lrho%charge = SUM(pab(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
715 : DO i = 1, nfa
716 : CALL lri_decomp_i(int3, lrii%cabai, i)
717 : lrho%tvec(i) = SUM(pab(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
718 : END DO
719 : IF (dab > lri_env%delta) THEN
720 : DO i = 1, nfb
721 : CALL lri_decomp_i(int3, lrii%cabbi, i)
722 : lrho%tvec(nfa + i) = SUM(pab(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
723 : END DO
724 : END IF
725 : !
726 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
727 : lrho%nst = SUM(lrho%tvec(1:nfa)*lrii%sn(1:nfa))
728 : ELSE
729 : lrho%nst = SUM(lrho%tvec(1:nn)*lrii%sn(1:nn))
730 : END IF
731 : lrho%lambda = (lrho%charge - lrho%nst)/lrii%nsn
732 : !
733 : ! solve the linear system of equations
734 : ALLOCATE (m(nn))
735 : m = 0._dp
736 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
737 : m(1:nfa) = lrho%tvec(1:nfa) + lrho%lambda*lrii%n(1:nfa)
738 : CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%sinv(1, 1), nfa, &
739 : m(1), 1, 0.0_dp, lrho%avec, 1)
740 : ELSE
741 : m(1:nn) = lrho%tvec(1:nn) + lrho%lambda*lrii%n(1:nn)
742 : CALL dgemv("N", nn, nn, 1.0_dp, lrii%sinv(1, 1), nn, &
743 : m(1), 1, 0.0_dp, lrho%avec, 1)
744 : END IF
745 : DEALLOCATE (m)
746 : END IF
747 :
748 : IF (lrii%lriff) THEN
749 : ! distant pair approximations
750 : ALLOCATE (paa(nba, nbb), pbb(nba, nbb))
751 : paa(1:nba, 1:nbb) = pab(1:nba, 1:nbb)*lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
752 : pbb(1:nba, 1:nbb) = pab(1:nba, 1:nbb)*(1._dp - lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb))
753 : !
754 : threshold = lri_env%eps_o3_int/MAX(SUM(ABS(paa(1:nba, 1:nbb))), 1.0e-14_dp)
755 : lrho%chargea = SUM(paa(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
756 : DO i = 1, nfa
757 : IF (lrii%abascr(i) > threshold) THEN
758 : CALL lri_decomp_i(int3, lrii%cabai, i)
759 : lrho%tveca(i) = SUM(paa(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
760 : ELSE
761 : lrho%tveca(i) = 0.0_dp
762 : END IF
763 : END DO
764 : threshold = lri_env%eps_o3_int/MAX(SUM(ABS(pbb(1:nba, 1:nbb))), 1.0e-14_dp)
765 : lrho%chargeb = SUM(pbb(1:nba, 1:nbb)*lrii%soo(1:nba, 1:nbb))
766 : DO i = 1, nfb
767 : IF (lrii%abbscr(i) > threshold) THEN
768 : CALL lri_decomp_i(int3, lrii%cabbi, i)
769 : lrho%tvecb(i) = SUM(pbb(1:nba, 1:nbb)*int3(1:nba, 1:nbb))
770 : ELSE
771 : lrho%tvecb(i) = 0.0_dp
772 : END IF
773 : END DO
774 : !
775 : lrho%nsta = SUM(lrho%tveca(1:nfa)*lrii%sna(1:nfa))
776 : lrho%nstb = SUM(lrho%tvecb(1:nfb)*lrii%snb(1:nfb))
777 : lrho%lambdaa = (lrho%chargea - lrho%nsta)/lrii%nsna
778 : lrho%lambdab = (lrho%chargeb - lrho%nstb)/lrii%nsnb
779 : ! solve the linear system of equations
780 : ALLOCATE (m(nfa))
781 : m(1:nfa) = lrho%tveca(1:nfa) + lrho%lambdaa*lrii%na(1:nfa)
782 : CALL dgemv("N", nfa, nfa, 1.0_dp, lrii%asinv(1, 1), nfa, &
783 : m(1), 1, 0.0_dp, lrho%aveca, 1)
784 : DEALLOCATE (m)
785 : ALLOCATE (m(nfb))
786 : m(1:nfb) = lrho%tvecb(1:nfb) + lrho%lambdab*lrii%nb(1:nfb)
787 : CALL dgemv("N", nfb, nfb, 1.0_dp, lrii%bsinv(1, 1), nfb, &
788 : m(1), 1, 0.0_dp, lrho%avecb, 1)
789 : DEALLOCATE (m)
790 : !
791 : DEALLOCATE (paa, pbb)
792 : END IF
793 :
794 : DEALLOCATE (pab, int3)
795 :
796 : END DO
797 : !$OMP END PARALLEL
798 1654 : CALL neighbor_list_iterator_release(nl_iterator)
799 :
800 : END DO
801 :
802 : END IF
803 :
804 810 : CALL timestop(handle)
805 :
806 1620 : END SUBROUTINE calculate_avec_lri
807 :
808 : ! **************************************************************************************************
809 : !> \brief sums up avec and distributes the fitted density on the grid
810 : !> \param lri_env the lri environment
811 : !> \param lri_density ...
812 : !> \param qs_env ...
813 : !> \param lri_rho_struct ...
814 : !> \param atomic_kind_set ...
815 : !> \param para_env ...
816 : !> \param response_density ...
817 : ! **************************************************************************************************
818 792 : SUBROUTINE distribute_lri_density_on_the_grid(lri_env, lri_density, qs_env, &
819 : lri_rho_struct, atomic_kind_set, para_env, &
820 : response_density)
821 :
822 : TYPE(lri_environment_type), POINTER :: lri_env
823 : TYPE(lri_density_type), POINTER :: lri_density
824 : TYPE(qs_environment_type), POINTER :: qs_env
825 : TYPE(qs_rho_type), INTENT(INOUT) :: lri_rho_struct
826 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
827 : TYPE(mp_para_env_type), POINTER :: para_env
828 : LOGICAL, INTENT(IN) :: response_density
829 :
830 : CHARACTER(LEN=*), PARAMETER :: routineN = 'distribute_lri_density_on_the_grid'
831 :
832 : INTEGER :: atom_a, atom_b, handle, iac, iatom, &
833 : ikind, ilist, ispin, jatom, jkind, &
834 : jneighbor, natom, nfa, nfb, nkind, &
835 : nspin
836 792 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
837 : REAL(KIND=dp) :: dab, fw, rab(3), str
838 792 : REAL(KIND=dp), DIMENSION(:), POINTER :: aci, acj, tot_rho_r
839 : TYPE(atomic_kind_type), POINTER :: atomic_kind
840 792 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
841 : TYPE(dbcsr_type) :: pmat_diag
842 : TYPE(lri_int_type), POINTER :: lrii
843 792 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
844 : TYPE(lri_list_type), POINTER :: lri_rho
845 : TYPE(lri_rhoab_type), POINTER :: lrho
846 : TYPE(neighbor_list_iterator_p_type), &
847 792 : DIMENSION(:), POINTER :: nl_iterator
848 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
849 792 : POINTER :: soo_list
850 792 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
851 792 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
852 : TYPE(qs_rho_type), POINTER :: rho
853 :
854 792 : CALL timeset(routineN, handle)
855 792 : NULLIFY (aci, acj, atomic_kind, lri_coef, lri_rho, &
856 792 : nl_iterator, soo_list, rho_r, rho_g, tot_rho_r)
857 :
858 792 : IF (ASSOCIATED(lri_env%soo_list)) THEN
859 792 : soo_list => lri_env%soo_list
860 :
861 792 : lri_env%stat%rho_tt = 0.0_dp
862 792 : lri_env%stat%rho_sr = 0.0_dp
863 792 : lri_env%stat%rho_ff = 0.0_dp
864 792 : lri_env%stat%rho_1c = 0.0_dp
865 :
866 792 : nspin = lri_density%nspin
867 792 : nkind = lri_env%lri_ints%nkind
868 :
869 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
870 792 : atom_of_kind=atom_of_kind)
871 :
872 : ! allocate the arrays to hold RI expansion coefficients lri_coefs
873 792 : CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set)
874 1600 : DO ispin = 1, nspin
875 :
876 808 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
877 808 : lri_rho => lri_density%lri_rhos(ispin)%lri_list
878 :
879 : ! sum up expansion coefficients
880 808 : CALL neighbor_list_iterator_create(nl_iterator, soo_list)
881 40328 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
882 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
883 39520 : iatom=iatom, jatom=jatom, ilist=ilist, inode=jneighbor, r=rab)
884 158080 : dab = SQRT(SUM(rab*rab))
885 39520 : atom_a = atom_of_kind(iatom)
886 39520 : atom_b = atom_of_kind(jatom)
887 39520 : aci => lri_coef(ikind)%acoef(atom_a, :)
888 39520 : acj => lri_coef(jkind)%acoef(atom_b, :)
889 39520 : iac = ikind + nkind*(jkind - 1)
890 39520 : lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
891 39520 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
892 39520 : nfa = lrho%nfa
893 39520 : nfb = lrho%nfb
894 :
895 39520 : IF (lrii%lrisr) THEN
896 14848 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
897 : !self pair aa
898 1461 : IF (.NOT. lri_env%exact_1c_terms) THEN
899 272245 : aci(1:nfa) = aci(1:nfa) + lrho%avec(1:nfa)
900 136781 : lri_env%stat%rho_sr = lri_env%stat%rho_sr + SUM(lrho%avec(:)*lrii%n(:))
901 : END IF
902 : ELSE
903 13387 : IF (iatom == jatom) THEN
904 : !periodic self pair aa'
905 5036 : fw = lrii%wsr
906 : ELSE
907 : !pairs ab
908 8351 : fw = 2.0_dp*lrii%wsr
909 : END IF
910 1726289 : aci(1:nfa) = aci(1:nfa) + fw*lrho%avec(1:nfa)
911 1726289 : acj(1:nfb) = acj(1:nfb) + fw*lrho%avec(nfa + 1:nfa + nfb)
912 1726289 : lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*SUM(lrho%avec(:)*lrii%n(:))
913 : END IF
914 : END IF
915 : !
916 40328 : IF (lrii%lriff) THEN
917 26412 : IF (iatom == jatom) THEN
918 5952 : fw = lrii%wff
919 : ELSE
920 20460 : fw = 2.0_dp*lrii%wff
921 : END IF
922 4984428 : aci(1:nfa) = aci(1:nfa) + fw*lrho%aveca(1:nfa)
923 4984428 : acj(1:nfb) = acj(1:nfb) + fw*lrho%avecb(1:nfb)
924 2505420 : lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*SUM(lrho%aveca(:)*lrii%na(:))
925 2505420 : lri_env%stat%rho_sr = lri_env%stat%rho_sr + fw*SUM(lrho%avecb(:)*lrii%nb(:))
926 : END IF
927 :
928 : END DO
929 808 : CALL neighbor_list_iterator_release(nl_iterator)
930 :
931 : ! replicate the acoef information
932 3210 : DO ikind = 1, nkind
933 1610 : atomic_kind => atomic_kind_set(ikind)
934 1610 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
935 5340 : DO iatom = 1, natom
936 2922 : aci => lri_coef(ikind)%acoef(iatom, :)
937 644308 : CALL para_env%sum(aci)
938 : END DO
939 : END DO
940 :
941 : END DO
942 :
943 : !distribute fitted density on the grid
944 : CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, &
945 792 : tot_rho_r=tot_rho_r)
946 1600 : DO ispin = 1, nspin
947 808 : IF (lri_env%exact_1c_terms) THEN
948 : ! 1 center terms (if requested)
949 48 : CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s)
950 48 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
951 48 : CALL dbcsr_create(pmat_diag, name="Block diagonal density", template=matrix_p(1, 1)%matrix)
952 48 : CALL dbcsr_get_block_diag(matrix_p(ispin, 1)%matrix, pmat_diag)
953 : str = 0.0_dp
954 48 : CALL dbcsr_dot(matrix_s(1, 1)%matrix, pmat_diag, str)
955 48 : lri_env%stat%rho_1c = lri_env%stat%rho_1c + str
956 48 : CALL dbcsr_replicate_all(pmat_diag)
957 : END IF
958 : !
959 808 : IF (.NOT. response_density) THEN
960 : CALL calculate_lri_rho_elec(rho_g(ispin), &
961 : rho_r(ispin), qs_env, &
962 : lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
963 636 : "LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag)
964 : ELSE
965 : CALL calculate_lri_rho_elec(rho_g(ispin), &
966 : rho_r(ispin), qs_env, &
967 : lri_density%lri_coefs(ispin)%lri_kinds, tot_rho_r(ispin), &
968 172 : "P_LRI_AUX", lri_env%exact_1c_terms, pmat=pmat_diag)
969 : END IF
970 808 : lri_env%stat%rho_tt = lri_env%stat%rho_tt + tot_rho_r(ispin)
971 : !
972 1600 : IF (lri_env%exact_1c_terms) CALL dbcsr_release(pmat_diag)
973 : END DO
974 :
975 : END IF
976 :
977 792 : CALL timestop(handle)
978 :
979 1584 : END SUBROUTINE distribute_lri_density_on_the_grid
980 :
981 : ! **************************************************************************************************
982 : !> \brief get inverse or pseudoinverse of lri overlap matrix for aux basis set
983 : !> \param lri_env lri environment
984 : !> \param sinv on exit its inverse
985 : !> \param sa ...
986 : !> \param sb ...
987 : !> \param sab ...
988 : ! **************************************************************************************************
989 1991 : SUBROUTINE inverse_lri_overlap(lri_env, sinv, sa, sb, sab)
990 :
991 : TYPE(lri_environment_type) :: lri_env
992 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sinv
993 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sa, sb, sab
994 :
995 : CHARACTER(LEN=*), PARAMETER :: routineN = 'inverse_lri_overlap'
996 :
997 : INTEGER :: handle, info, n, nfa, nfb, nn
998 1991 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
999 : REAL(KIND=dp) :: anorm, delta, rcond, rskip
1000 1991 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
1001 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: s
1002 : REAL(KIND=dp), EXTERNAL :: dlange
1003 :
1004 1991 : CALL timeset(routineN, handle)
1005 :
1006 1991 : nfa = SIZE(sa, 1)
1007 1991 : nfb = SIZE(sb, 1)
1008 1991 : nn = nfa + nfb
1009 1991 : n = SIZE(sinv, 1)
1010 1991 : CPASSERT(n == nn)
1011 7964 : ALLOCATE (s(n, n))
1012 10979115 : s(1:nfa, 1:nfa) = sa(1:nfa, 1:nfa)
1013 10190395 : s(1:nfa, nfa + 1:nn) = sab(1:nfa, 1:nfb)
1014 10190395 : s(nfa + 1:nn, 1:nfa) = TRANSPOSE(sab(1:nfa, 1:nfb))
1015 10979115 : s(nfa + 1:nn, nfa + 1:nn) = sb(1:nfb, 1:nfb)
1016 :
1017 1991 : rskip = 1.E-8_dp ! parameter for pseudo inverse
1018 :
1019 42076897 : sinv(:, :) = s
1020 3761 : SELECT CASE (lri_env%lri_overlap_inv)
1021 : CASE (do_lri_inv)
1022 1770 : CALL invmat_symm(sinv)
1023 : CASE (do_lri_pseudoinv_svd)
1024 0 : CALL get_pseudo_inverse_svd(s, sinv, rskip)
1025 : CASE (do_lri_pseudoinv_diag)
1026 1 : CALL get_pseudo_inverse_diag(s, sinv, rskip)
1027 : CASE (do_lri_inv_auto)
1028 : ! decide whether to calculate inverse or pseudoinverse
1029 1100 : ALLOCATE (work(3*n), iwork(n))
1030 : ! norm of matrix
1031 220 : anorm = dlange('1', n, n, sinv, n, work)
1032 : ! Cholesky factorization (fails if matrix not positive definite)
1033 220 : CALL dpotrf('U', n, sinv, n, info)
1034 220 : IF (info == 0) THEN
1035 : ! condition number
1036 220 : CALL dpocon('U', n, sinv, n, anorm, rcond, work, iwork, info)
1037 220 : IF (info /= 0) THEN
1038 0 : CPABORT("DPOCON failed")
1039 : END IF
1040 220 : IF (LOG(1._dp/rcond) > lri_env%cond_max) THEN
1041 3 : CALL get_pseudo_inverse_svd(s, sinv, rskip)
1042 : ELSE
1043 217 : CALL invmat_symm(sinv, potrf=.FALSE.)
1044 : END IF
1045 : ELSE
1046 0 : CALL get_pseudo_inverse_svd(s, sinv, rskip)
1047 : END IF
1048 220 : DEALLOCATE (work, iwork)
1049 : CASE DEFAULT
1050 1991 : CPABORT("No initialization for LRI overlap available?")
1051 : END SELECT
1052 :
1053 1991 : delta = inv_test(s, sinv)
1054 1991 : !$OMP CRITICAL(sum_critical)
1055 1991 : lri_env%stat%overlap_error = MAX(delta, lri_env%stat%overlap_error)
1056 : !$OMP END CRITICAL(sum_critical)
1057 :
1058 1991 : DEALLOCATE (s)
1059 :
1060 1991 : CALL timestop(handle)
1061 :
1062 1991 : END SUBROUTINE inverse_lri_overlap
1063 :
1064 : ! **************************************************************************************************
1065 : !> \brief ...
1066 : !> \param amat ...
1067 : !> \param ainv ...
1068 : !> \return ...
1069 : ! **************************************************************************************************
1070 1991 : FUNCTION inv_test(amat, ainv) RESULT(delta)
1071 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: amat, ainv
1072 : REAL(KIND=dp) :: delta
1073 :
1074 : INTEGER :: i, n
1075 1991 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1076 :
1077 1991 : n = SIZE(amat, 1)
1078 7964 : ALLOCATE (work(n, n))
1079 42122429 : work(1:n, 1:n) = MATMUL(amat(1:n, 1:n), ainv(1:n, 1:n))
1080 258141 : DO i = 1, n
1081 258141 : work(i, i) = work(i, i) - 1.0_dp
1082 : END DO
1083 42076897 : delta = MAXVAL(ABS(work))
1084 1991 : DEALLOCATE (work)
1085 1991 : END FUNCTION inv_test
1086 :
1087 : ! **************************************************************************************************
1088 : !> \brief debug output
1089 : !> \param lri_env ...
1090 : !> \param qs_env ...
1091 : !> \param lri_ints ...
1092 : !> \param soo_list ...
1093 : ! **************************************************************************************************
1094 2 : SUBROUTINE output_debug_info(lri_env, qs_env, lri_ints, soo_list)
1095 :
1096 : TYPE(lri_environment_type), POINTER :: lri_env
1097 : TYPE(qs_environment_type), POINTER :: qs_env
1098 : TYPE(lri_list_type), POINTER :: lri_ints
1099 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1100 : POINTER :: soo_list
1101 :
1102 : CHARACTER(LEN=*), PARAMETER :: routineN = 'output_debug_info'
1103 :
1104 : INTEGER :: handle, iac, ikind, ilist, iunit, jkind, &
1105 : jneighbor, nkind
1106 : REAL(KIND=dp) :: dmax_aabb, dmax_ab, dmax_aba, dmax_abb, &
1107 : dmax_oo
1108 : TYPE(cp_logger_type), POINTER :: logger
1109 : TYPE(dft_control_type), POINTER :: dft_control
1110 : TYPE(lri_int_rho_type), POINTER :: lriir
1111 : TYPE(lri_int_type), POINTER :: lrii
1112 : TYPE(mp_para_env_type), POINTER :: para_env
1113 : TYPE(neighbor_list_iterator_p_type), &
1114 2 : DIMENSION(:), POINTER :: nl_iterator
1115 : TYPE(section_vals_type), POINTER :: input
1116 :
1117 2 : CALL timeset(routineN, handle)
1118 2 : NULLIFY (input, logger, lrii, lriir, nl_iterator, para_env)
1119 : CALL get_qs_env(qs_env, dft_control=dft_control, input=input, nkind=nkind, &
1120 2 : para_env=para_env)
1121 2 : dmax_ab = 0._dp
1122 2 : dmax_oo = 0._dp
1123 2 : dmax_aba = 0._dp
1124 2 : dmax_abb = 0._dp
1125 2 : dmax_aabb = 0._dp
1126 :
1127 2 : CALL neighbor_list_iterator_create(nl_iterator, soo_list)
1128 5 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1129 :
1130 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
1131 3 : ilist=ilist, inode=jneighbor)
1132 :
1133 3 : iac = ikind + nkind*(jkind - 1)
1134 3 : lrii => lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
1135 :
1136 3 : dmax_ab = MAX(dmax_ab, lrii%dmax_ab)
1137 3 : dmax_oo = MAX(dmax_oo, lrii%dmax_oo)
1138 3 : dmax_aba = MAX(dmax_aba, lrii%dmax_aba)
1139 3 : dmax_abb = MAX(dmax_abb, lrii%dmax_abb)
1140 :
1141 5 : IF (dft_control%qs_control%lri_optbas) THEN
1142 3 : lriir => lri_env%lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
1143 3 : dmax_aabb = MAX(dmax_aabb, lriir%dmax_aabb)
1144 : END IF
1145 :
1146 : END DO
1147 :
1148 2 : CALL neighbor_list_iterator_release(nl_iterator)
1149 2 : CALL para_env%max(dmax_ab)
1150 2 : CALL para_env%max(dmax_oo)
1151 2 : CALL para_env%max(dmax_aba)
1152 2 : CALL para_env%max(dmax_abb)
1153 2 : CALL para_env%max(dmax_aabb)
1154 :
1155 2 : logger => cp_get_default_logger()
1156 : iunit = cp_print_key_unit_nr(logger, input, "PRINT%PROGRAM_RUN_INFO", &
1157 2 : extension=".lridebug")
1158 :
1159 2 : IF (iunit > 0) THEN
1160 1 : WRITE (iunit, FMT="(/,T2,A)") "DEBUG INFO FOR LRI INTEGRALS"
1161 : WRITE (iunit, FMT="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1162 1 : "[ai|bi]; fit basis", dmax_ab
1163 : WRITE (iunit, FMT="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1164 1 : "[a|b]; orbital basis", dmax_oo
1165 : WRITE (iunit, FMT="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1166 1 : "[a|b|ai]", dmax_aba
1167 : WRITE (iunit, FMT="(T2,A,T69,ES12.5)") "Maximal deviation of integrals "// &
1168 1 : "[a|b|bi]", dmax_abb
1169 1 : IF (dft_control%qs_control%lri_optbas) THEN
1170 : WRITE (iunit, FMT="(T2,A,T69,ES12.5,/)") "Maximal deviation of integrals "// &
1171 1 : "[aa|bb]; orbital basis", &
1172 2 : dmax_aabb
1173 : END IF
1174 : END IF
1175 :
1176 : CALL cp_print_key_finished_output(iunit, logger, input, &
1177 2 : "PRINT%PROGRAM_RUN_INFO")
1178 2 : CALL timestop(handle)
1179 :
1180 2 : END SUBROUTINE output_debug_info
1181 :
1182 : ! **************************************************************************************************
1183 : !> \brief ...
1184 : !> \param qs_env ...
1185 : !> \param lri_v_int ...
1186 : !> \param calculate_forces ...
1187 : ! **************************************************************************************************
1188 8 : SUBROUTINE v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1189 : TYPE(qs_environment_type), POINTER :: qs_env
1190 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1191 : LOGICAL, INTENT(IN) :: calculate_forces
1192 :
1193 : INTEGER :: ikind, natom, nfa, nkind
1194 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: v_int
1195 : TYPE(lri_environment_type), POINTER :: lri_env
1196 :
1197 8 : CALL get_qs_env(qs_env, lri_env=lri_env, nkind=nkind)
1198 :
1199 24 : DO ikind = 1, nkind
1200 16 : natom = SIZE(lri_v_int(ikind)%v_int, 1)
1201 16 : nfa = SIZE(lri_v_int(ikind)%v_int, 2)
1202 16 : v_int => lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int
1203 16 : CPASSERT(SIZE(v_int, 1) == natom)
1204 16 : CPASSERT(SIZE(v_int, 2) == nfa)
1205 8568 : lri_v_int(ikind)%v_int(:, :) = lri_v_int(ikind)%v_int(:, :) + v_int(:, :)
1206 : END DO
1207 :
1208 8 : IF (calculate_forces) THEN
1209 2 : CALL calculate_lri_ppl_integrals(lri_env, qs_env, calculate_forces)
1210 : END IF
1211 :
1212 8 : END SUBROUTINE v_int_ppl_update
1213 :
1214 : ! **************************************************************************************************
1215 : !> \brief ...
1216 : !> \param qs_env ...
1217 : !> \param lri_v_int ...
1218 : !> \param ecore_ppl_ri ...
1219 : ! **************************************************************************************************
1220 8 : SUBROUTINE v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl_ri)
1221 : TYPE(qs_environment_type), POINTER :: qs_env
1222 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1223 : REAL(KIND=dp), INTENT(INOUT) :: ecore_ppl_ri
1224 :
1225 : INTEGER :: ikind, natom, nfa, nkind
1226 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: v_int
1227 : TYPE(lri_environment_type), POINTER :: lri_env
1228 :
1229 8 : CALL get_qs_env(qs_env, lri_env=lri_env, nkind=nkind)
1230 24 : DO ikind = 1, nkind
1231 16 : natom = SIZE(lri_v_int(ikind)%v_int, 1)
1232 16 : nfa = SIZE(lri_v_int(ikind)%v_int, 2)
1233 16 : v_int => lri_env%lri_ppl_ints%lri_ppl(ikind)%v_int
1234 16 : CPASSERT(SIZE(v_int, 1) == natom)
1235 16 : CPASSERT(SIZE(v_int, 2) == nfa)
1236 4288 : ecore_ppl_ri = ecore_ppl_ri + SUM(v_int(:, :)*lri_v_int(ikind)%acoef(:, :))
1237 : END DO
1238 :
1239 8 : END SUBROUTINE v_int_ppl_energy
1240 :
1241 : ! **************************************************************************************************
1242 : !> \brief ...
1243 : !> \param qs_env ...
1244 : !> \param ltddfpt ...
1245 : !> \param tddfpt_lri_env ...
1246 : ! **************************************************************************************************
1247 68 : SUBROUTINE lri_print_stat(qs_env, ltddfpt, tddfpt_lri_env)
1248 :
1249 : TYPE(qs_environment_type), POINTER :: qs_env
1250 : LOGICAL, OPTIONAL :: ltddfpt
1251 : TYPE(lri_environment_type), OPTIONAL, POINTER :: tddfpt_lri_env
1252 :
1253 : INTEGER :: iunit
1254 : LOGICAL :: my_ltddfpt
1255 : REAL(KIND=dp) :: abai_mem, coef_mem, oint_mem, overlap_error, pairs_ff, pairs_sr, pairs_tt, &
1256 : ppli_mem, ppx, rho_1c, rho_ff, rho_sr, rho_tt, rhos_mem
1257 : TYPE(cp_logger_type), POINTER :: logger
1258 : TYPE(lri_environment_type), POINTER :: lri_env
1259 : TYPE(mp_para_env_type), POINTER :: para_env
1260 : TYPE(section_vals_type), POINTER :: input
1261 :
1262 68 : my_ltddfpt = .FALSE.
1263 68 : IF (PRESENT(ltddfpt)) my_ltddfpt = ltddfpt
1264 :
1265 10 : IF (.NOT. my_ltddfpt) THEN
1266 58 : CALL get_qs_env(qs_env, lri_env=lri_env, input=input, para_env=para_env)
1267 : ELSE
1268 10 : CALL get_qs_env(qs_env, input=input, para_env=para_env)
1269 : NULLIFY (lri_env)
1270 10 : lri_env => tddfpt_lri_env
1271 : END IF
1272 :
1273 68 : IF (lri_env%statistics) THEN
1274 12 : logger => cp_get_default_logger()
1275 12 : iunit = cp_print_key_unit_nr(logger, input, "PRINT%PROGRAM_RUN_INFO", extension=".Log")
1276 12 : pairs_tt = lri_env%stat%pairs_tt
1277 12 : CALL para_env%sum(pairs_tt)
1278 12 : pairs_sr = lri_env%stat%pairs_sr
1279 12 : CALL para_env%sum(pairs_sr)
1280 12 : pairs_ff = lri_env%stat%pairs_ff
1281 12 : CALL para_env%sum(pairs_ff)
1282 12 : overlap_error = lri_env%stat%overlap_error
1283 12 : CALL para_env%sum(overlap_error)
1284 12 : rho_tt = -lri_env%stat%rho_tt
1285 12 : rho_sr = lri_env%stat%rho_sr
1286 12 : CALL para_env%sum(rho_sr)
1287 12 : rho_ff = lri_env%stat%rho_ff
1288 12 : CALL para_env%sum(rho_ff)
1289 12 : IF (lri_env%exact_1c_terms) THEN
1290 0 : rho_1c = lri_env%stat%rho_1c
1291 : ELSE
1292 12 : rho_1c = 0.0_dp
1293 : END IF
1294 12 : ppx = 8.e-6_dp
1295 12 : coef_mem = lri_env%stat%coef_mem*ppx
1296 12 : CALL para_env%sum(coef_mem)
1297 12 : oint_mem = lri_env%stat%oint_mem*ppx
1298 12 : CALL para_env%sum(oint_mem)
1299 12 : rhos_mem = lri_env%stat%rhos_mem*ppx
1300 12 : CALL para_env%sum(rhos_mem)
1301 12 : abai_mem = lri_env%stat%abai_mem*ppx
1302 12 : CALL para_env%sum(abai_mem)
1303 12 : IF (lri_env%ppl_ri) THEN
1304 2 : ppli_mem = lri_env%stat%ppli_mem*ppx
1305 2 : CALL para_env%sum(ppli_mem)
1306 : ELSE
1307 10 : ppli_mem = 0.0_dp
1308 : END IF
1309 12 : IF (iunit > 0) THEN
1310 : !
1311 6 : WRITE (iunit, FMT="(/,T2,A,A,A)") "********************************", &
1312 12 : " LRI STATISTICS ", "*******************************"
1313 : !
1314 6 : WRITE (iunit, FMT="(T4,A,T63,F16.0)") "Total number of atom pairs", pairs_tt
1315 6 : ppx = pairs_sr/pairs_tt*100._dp
1316 6 : WRITE (iunit, FMT="(T4,A,T52,A,F5.1,A,T63,F16.0)") "Near field atom pairs", &
1317 12 : "[", ppx, "%]", pairs_sr
1318 6 : ppx = pairs_ff/pairs_tt*100._dp
1319 6 : WRITE (iunit, FMT="(T4,A,T52,A,F5.1,A,T63,F16.0)") "Far field atom pairs", &
1320 12 : "[", ppx, "%]", pairs_ff
1321 : !
1322 6 : WRITE (iunit, FMT="(T4,A,T63,G16.8)") "Max. error in pair overlap inversion", overlap_error
1323 : !
1324 6 : WRITE (iunit, FMT="(T4,A,T63,F16.2)") "Total charge approximated", rho_tt
1325 6 : ppx = rho_sr/rho_tt*100._dp
1326 6 : WRITE (iunit, FMT="(T4,A,T52,A,F5.1,A,T63,F16.2)") "Near field charge", &
1327 12 : "[", ppx, "%]", rho_sr
1328 6 : ppx = rho_ff/rho_tt*100._dp
1329 6 : WRITE (iunit, FMT="(T4,A,T52,A,F5.1,A,T63,F16.2)") "Far field charge", &
1330 12 : "[", ppx, "%]", rho_ff
1331 6 : IF (lri_env%exact_1c_terms) THEN
1332 0 : ppx = rho_1c/rho_tt*100._dp
1333 0 : WRITE (iunit, FMT="(T4,A,T52,A,F5.1,A,T63,F16.2)") "One center charge", &
1334 0 : "[", ppx, "%]", rho_1c
1335 : END IF
1336 : !
1337 6 : WRITE (iunit, FMT="(T4,A,T63,F9.0,A)") "Max. memory/task for expansion coeficients", coef_mem, " Mbytes"
1338 6 : WRITE (iunit, FMT="(T4,A,T63,F9.0,A)") "Max. memory/task for overlap matrices", oint_mem, " Mbytes"
1339 6 : WRITE (iunit, FMT="(T4,A,T63,F9.0,A)") "Max. memory/task for density expansions", rhos_mem, " Mbytes"
1340 6 : WRITE (iunit, FMT="(T4,A,T63,F9.0,A)") "Max. memory/task for aba/abb integrals", abai_mem, " Mbytes"
1341 6 : IF (lri_env%ppl_ri) THEN
1342 1 : WRITE (iunit, FMT="(T4,A,T63,F9.0,A)") "Max. memory/task for ppl integrals", ppli_mem, " Mbytes"
1343 : END IF
1344 : !
1345 6 : WRITE (iunit, FMT="(T2,A,A)") "********************************", &
1346 12 : "***********************************************"
1347 : !
1348 : END IF
1349 :
1350 12 : CALL cp_print_key_finished_output(iunit, logger, input, "PRINT%PROGRAM_RUN_INFO")
1351 : END IF
1352 :
1353 68 : END SUBROUTINE lri_print_stat
1354 :
1355 : END MODULE lri_environment_methods
|
