Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility functions for the perturbation calculations.
10 : !> \note
11 : !> - routines are programmed with spins in mind
12 : !> but are as of now not tested with them
13 : !> \par History
14 : !> 22-08-2002, TCH, started development
15 : ! **************************************************************************************************
16 : MODULE qs_p_env_methods
17 : USE admm_methods, ONLY: admm_aux_response_density
18 : USE admm_types, ONLY: admm_gapw_r3d_rs_type,&
19 : admm_type,&
20 : get_admm_env
21 : USE atomic_kind_types, ONLY: atomic_kind_type
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: dbcsr_add,&
25 : dbcsr_copy,&
26 : dbcsr_p_type,&
27 : dbcsr_release,&
28 : dbcsr_scale,&
29 : dbcsr_set,&
30 : dbcsr_type
31 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
32 : cp_dbcsr_plus_fm_fm_t,&
33 : cp_dbcsr_sm_fm_multiply,&
34 : dbcsr_allocate_matrix_set
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_multiply
36 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
37 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
38 : cp_fm_pool_type,&
39 : fm_pool_create_fm,&
40 : fm_pool_get_el_struct,&
41 : fm_pool_give_back_fm,&
42 : fm_pools_create_fm_vect
43 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
44 : cp_fm_struct_get,&
45 : cp_fm_struct_release,&
46 : cp_fm_struct_type
47 : USE cp_fm_types, ONLY: cp_fm_create,&
48 : cp_fm_get_info,&
49 : cp_fm_release,&
50 : cp_fm_set_all,&
51 : cp_fm_to_fm,&
52 : cp_fm_type
53 : USE cp_log_handling, ONLY: cp_get_default_logger,&
54 : cp_logger_type,&
55 : cp_to_string
56 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
57 : cp_print_key_unit_nr
58 : USE hartree_local_methods, ONLY: init_coulomb_local
59 : USE hartree_local_types, ONLY: hartree_local_create
60 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
61 : ot_precond_none
62 : USE input_section_types, ONLY: section_vals_get,&
63 : section_vals_get_subs_vals,&
64 : section_vals_type
65 : USE kinds, ONLY: default_string_length,&
66 : dp
67 : USE message_passing, ONLY: mp_para_env_type
68 : USE parallel_gemm_api, ONLY: parallel_gemm
69 : USE preconditioner_types, ONLY: init_preconditioner
70 : USE pw_env_types, ONLY: pw_env_type
71 : USE pw_types, ONLY: pw_c1d_gs_type,&
72 : pw_r3d_rs_type
73 : USE qs_collocate_density, ONLY: calculate_rho_elec
74 : USE qs_energy_types, ONLY: qs_energy_type
75 : USE qs_environment_types, ONLY: get_qs_env,&
76 : qs_environment_type
77 : USE qs_kind_types, ONLY: qs_kind_type
78 : USE qs_kpp1_env_methods, ONLY: kpp1_create,&
79 : kpp1_did_change
80 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
81 : USE qs_ks_types, ONLY: qs_ks_did_change,&
82 : qs_ks_env_type
83 : USE qs_linres_types, ONLY: linres_control_type
84 : USE qs_local_rho_types, ONLY: local_rho_set_create
85 : USE qs_matrix_pools, ONLY: mpools_get
86 : USE qs_mo_types, ONLY: get_mo_set,&
87 : mo_set_type
88 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
89 : USE qs_p_env_types, ONLY: qs_p_env_type
90 : USE qs_rho0_ggrid, ONLY: rho0_s_grid_create
91 : USE qs_rho0_methods, ONLY: init_rho0
92 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
93 : calculate_rho_atom_coeff
94 : USE qs_rho_methods, ONLY: qs_rho_rebuild,&
95 : qs_rho_update_rho
96 : USE qs_rho_types, ONLY: qs_rho_create,&
97 : qs_rho_get,&
98 : qs_rho_type
99 : USE string_utilities, ONLY: compress
100 : USE task_list_types, ONLY: task_list_type
101 : #include "./base/base_uses.f90"
102 :
103 : IMPLICIT NONE
104 :
105 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_p_env_methods'
106 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
107 :
108 : PRIVATE
109 : PUBLIC :: p_env_create, p_env_psi0_changed
110 : PUBLIC :: p_preortho, p_postortho
111 : PUBLIC :: p_env_check_i_alloc, p_env_update_rho
112 : PUBLIC :: p_env_finish_kpp1
113 :
114 : CONTAINS
115 :
116 : ! **************************************************************************************************
117 : !> \brief allocates and initializes the perturbation environment (no setup)
118 : !> \param p_env the environment to initialize
119 : !> \param qs_env the qs_environment for the system
120 : !> \param p1_option ...
121 : !> \param p1_admm_option ...
122 : !> \param orthogonal_orbitals if the orbitals are orthogonal
123 : !> \param linres_control ...
124 : !> \par History
125 : !> 07.2002 created [fawzi]
126 : !> \author Fawzi Mohamed
127 : ! **************************************************************************************************
128 1830 : SUBROUTINE p_env_create(p_env, qs_env, p1_option, p1_admm_option, &
129 : orthogonal_orbitals, linres_control)
130 :
131 : TYPE(qs_p_env_type) :: p_env
132 : TYPE(qs_environment_type), POINTER :: qs_env
133 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
134 : POINTER :: p1_option, p1_admm_option
135 : LOGICAL, INTENT(in), OPTIONAL :: orthogonal_orbitals
136 : TYPE(linres_control_type), OPTIONAL, POINTER :: linres_control
137 :
138 : CHARACTER(len=*), PARAMETER :: routineN = 'p_env_create'
139 :
140 : INTEGER :: handle, n_ao, n_mo, n_spins, natom, spin
141 : TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
142 : TYPE(admm_type), POINTER :: admm_env
143 1830 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
144 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
145 1830 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools, mo_mo_fm_pools
146 : TYPE(cp_fm_type), POINTER :: qs_env_c
147 1830 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux_fit
148 : TYPE(dft_control_type), POINTER :: dft_control
149 : TYPE(mp_para_env_type), POINTER :: para_env
150 : TYPE(pw_env_type), POINTER :: pw_env
151 1830 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
152 :
153 1830 : CALL timeset(routineN, handle)
154 1830 : NULLIFY (ao_mo_fm_pools, mo_mo_fm_pools, matrix_s, dft_control, para_env, blacs_env)
155 : CALL get_qs_env(qs_env, &
156 : matrix_s=matrix_s, &
157 : dft_control=dft_control, &
158 : para_env=para_env, &
159 1830 : blacs_env=blacs_env)
160 :
161 1830 : n_spins = dft_control%nspins
162 :
163 1830 : p_env%new_preconditioner = .TRUE.
164 :
165 1830 : ALLOCATE (p_env%rho1)
166 1830 : CALL qs_rho_create(p_env%rho1)
167 1830 : ALLOCATE (p_env%rho1_xc)
168 1830 : CALL qs_rho_create(p_env%rho1_xc)
169 :
170 1830 : ALLOCATE (p_env%kpp1_env)
171 1830 : CALL kpp1_create(p_env%kpp1_env)
172 :
173 1830 : IF (PRESENT(p1_option)) THEN
174 272 : p_env%p1 => p1_option
175 : ELSE
176 1558 : CALL dbcsr_allocate_matrix_set(p_env%p1, n_spins)
177 3342 : DO spin = 1, n_spins
178 1784 : ALLOCATE (p_env%p1(spin)%matrix)
179 : CALL dbcsr_copy(p_env%p1(spin)%matrix, matrix_s(1)%matrix, &
180 1784 : name="p_env%p1-"//TRIM(ADJUSTL(cp_to_string(spin))))
181 3342 : CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp)
182 : END DO
183 : END IF
184 :
185 1830 : IF (dft_control%do_admm) THEN
186 340 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
187 340 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
188 206 : ALLOCATE (p_env%rho1_admm)
189 206 : CALL qs_rho_create(p_env%rho1_admm)
190 : END IF
191 :
192 340 : IF (PRESENT(p1_admm_option)) THEN
193 0 : p_env%p1_admm => p1_admm_option
194 : ELSE
195 340 : CALL dbcsr_allocate_matrix_set(p_env%p1_admm, n_spins)
196 726 : DO spin = 1, n_spins
197 386 : ALLOCATE (p_env%p1_admm(spin)%matrix)
198 : CALL dbcsr_copy(p_env%p1_admm(spin)%matrix, matrix_s_aux_fit(1)%matrix, &
199 386 : name="p_env%p1_admm-"//TRIM(ADJUSTL(cp_to_string(spin))))
200 726 : CALL dbcsr_set(p_env%p1_admm(spin)%matrix, 0.0_dp)
201 : END DO
202 : END IF
203 340 : CALL get_qs_env(qs_env, admm_env=admm_env)
204 340 : IF (admm_env%do_gapw) THEN
205 54 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
206 54 : admm_gapw_env => admm_env%admm_gapw_env
207 54 : CALL local_rho_set_create(p_env%local_rho_set_admm)
208 : CALL allocate_rho_atom_internals(p_env%local_rho_set_admm%rho_atom_set, atomic_kind_set, &
209 54 : admm_gapw_env%admm_kind_set, dft_control, para_env)
210 : END IF
211 : END IF
212 :
213 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools, &
214 1830 : mo_mo_fm_pools=mo_mo_fm_pools)
215 :
216 5490 : p_env%n_mo = 0
217 5490 : p_env%n_ao = 0
218 3972 : DO spin = 1, n_spins
219 2142 : CALL get_mo_set(qs_env%mos(spin), mo_coeff=qs_env_c)
220 : CALL cp_fm_get_info(qs_env_c, &
221 2142 : ncol_global=n_mo, nrow_global=n_ao)
222 2142 : p_env%n_mo(spin) = n_mo
223 3972 : p_env%n_ao(spin) = n_ao
224 : END DO
225 :
226 1830 : p_env%orthogonal_orbitals = .FALSE.
227 1830 : IF (PRESENT(orthogonal_orbitals)) THEN
228 1830 : p_env%orthogonal_orbitals = orthogonal_orbitals
229 : END IF
230 :
231 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%S_psi0, &
232 1830 : name="p_env%S_psi0")
233 :
234 : ! alloc m_epsilon
235 : CALL fm_pools_create_fm_vect(mo_mo_fm_pools, elements=p_env%m_epsilon, &
236 1830 : name="p_env%m_epsilon")
237 :
238 : ! alloc Smo_inv
239 1830 : IF (.NOT. p_env%orthogonal_orbitals) THEN
240 : CALL fm_pools_create_fm_vect(mo_mo_fm_pools, elements=p_env%Smo_inv, &
241 0 : name="p_env%Smo_inv")
242 : END IF
243 :
244 1830 : IF (.NOT. p_env%orthogonal_orbitals) THEN
245 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, &
246 : elements=p_env%psi0d, &
247 0 : name="p_env%psi0d")
248 : END IF
249 :
250 : !------------------------------!
251 : ! GAPW/GAPW_XC initializations !
252 : !------------------------------!
253 1830 : IF (dft_control%qs_control%gapw) THEN
254 : CALL get_qs_env(qs_env, &
255 : atomic_kind_set=atomic_kind_set, &
256 : natom=natom, &
257 : pw_env=pw_env, &
258 314 : qs_kind_set=qs_kind_set)
259 :
260 314 : CALL local_rho_set_create(p_env%local_rho_set)
261 : CALL allocate_rho_atom_internals(p_env%local_rho_set%rho_atom_set, atomic_kind_set, &
262 314 : qs_kind_set, dft_control, para_env)
263 :
264 : CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
265 314 : zcore=0.0_dp)
266 314 : CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole)
267 314 : CALL hartree_local_create(p_env%hartree_local)
268 314 : CALL init_coulomb_local(p_env%hartree_local, natom)
269 1516 : ELSE IF (dft_control%qs_control%gapw_xc) THEN
270 : CALL get_qs_env(qs_env, &
271 : atomic_kind_set=atomic_kind_set, &
272 54 : qs_kind_set=qs_kind_set)
273 54 : CALL local_rho_set_create(p_env%local_rho_set)
274 : CALL allocate_rho_atom_internals(p_env%local_rho_set%rho_atom_set, atomic_kind_set, &
275 54 : qs_kind_set, dft_control, para_env)
276 : END IF
277 :
278 : !------------------------!
279 : ! LINRES initializations !
280 : !------------------------!
281 1830 : IF (PRESENT(linres_control)) THEN
282 :
283 1830 : IF (linres_control%preconditioner_type /= ot_precond_none) THEN
284 : ! Initialize the preconditioner matrix
285 1826 : IF (.NOT. ASSOCIATED(p_env%preconditioner)) THEN
286 :
287 7616 : ALLOCATE (p_env%preconditioner(n_spins))
288 3964 : DO spin = 1, n_spins
289 : CALL init_preconditioner(p_env%preconditioner(spin), &
290 3964 : para_env=para_env, blacs_env=blacs_env)
291 : END DO
292 :
293 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%PS_psi0, &
294 1826 : name="p_env%PS_psi0")
295 : END IF
296 : END IF
297 :
298 : END IF
299 :
300 1830 : CALL timestop(handle)
301 :
302 1830 : END SUBROUTINE p_env_create
303 :
304 : ! **************************************************************************************************
305 : !> \brief checks that the intenal storage is allocated, and allocs it if needed
306 : !> \param p_env the environment to check
307 : !> \param qs_env the qs environment this p_env lives in
308 : !> \par History
309 : !> 12.2002 created [fawzi]
310 : !> \author Fawzi Mohamed
311 : !> \note
312 : !> private routine
313 : ! **************************************************************************************************
314 10142 : SUBROUTINE p_env_check_i_alloc(p_env, qs_env)
315 : TYPE(qs_p_env_type) :: p_env
316 : TYPE(qs_environment_type), POINTER :: qs_env
317 :
318 : CHARACTER(len=*), PARAMETER :: routineN = 'p_env_check_i_alloc'
319 :
320 : CHARACTER(len=25) :: name
321 : INTEGER :: handle, ispin, nspins
322 : LOGICAL :: gapw_xc
323 10142 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
324 : TYPE(dft_control_type), POINTER :: dft_control
325 :
326 10142 : CALL timeset(routineN, handle)
327 :
328 10142 : NULLIFY (dft_control, matrix_s)
329 :
330 10142 : CALL get_qs_env(qs_env, dft_control=dft_control)
331 10142 : gapw_xc = dft_control%qs_control%gapw_xc
332 10142 : IF (.NOT. ASSOCIATED(p_env%kpp1)) THEN
333 1610 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
334 1610 : nspins = dft_control%nspins
335 :
336 1610 : CALL dbcsr_allocate_matrix_set(p_env%kpp1, nspins)
337 1610 : name = "p_env%kpp1-"
338 1610 : CALL compress(name, full=.TRUE.)
339 3454 : DO ispin = 1, nspins
340 1844 : ALLOCATE (p_env%kpp1(ispin)%matrix)
341 : CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, matrix_s(1)%matrix, &
342 1844 : name=TRIM(name)//ADJUSTL(cp_to_string(ispin)))
343 3454 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
344 : END DO
345 :
346 1610 : CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env)
347 1610 : IF (gapw_xc) THEN
348 52 : CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env)
349 : END IF
350 :
351 : END IF
352 :
353 10142 : IF (dft_control%do_admm .AND. .NOT. ASSOCIATED(p_env%kpp1_admm)) THEN
354 340 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
355 340 : nspins = dft_control%nspins
356 :
357 340 : CALL dbcsr_allocate_matrix_set(p_env%kpp1_admm, nspins)
358 340 : name = "p_env%kpp1_admm-"
359 340 : CALL compress(name, full=.TRUE.)
360 726 : DO ispin = 1, nspins
361 386 : ALLOCATE (p_env%kpp1_admm(ispin)%matrix)
362 : CALL dbcsr_copy(p_env%kpp1_admm(ispin)%matrix, matrix_s(1)%matrix, &
363 386 : name=TRIM(name)//ADJUSTL(cp_to_string(ispin)))
364 726 : CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
365 : END DO
366 :
367 340 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
368 206 : CALL qs_rho_rebuild(p_env%rho1_admm, qs_env=qs_env, admm=.TRUE.)
369 : END IF
370 :
371 : END IF
372 :
373 10142 : IF (.NOT. ASSOCIATED(p_env%rho1)) THEN
374 0 : CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env)
375 0 : IF (gapw_xc) THEN
376 0 : CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env)
377 : END IF
378 :
379 0 : IF (dft_control%do_admm) THEN
380 0 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
381 0 : CALL qs_rho_rebuild(p_env%rho1_admm, qs_env=qs_env, admm=.TRUE.)
382 : END IF
383 : END IF
384 :
385 : END IF
386 :
387 10142 : CALL timestop(handle)
388 10142 : END SUBROUTINE p_env_check_i_alloc
389 :
390 : ! **************************************************************************************************
391 : !> \brief ...
392 : !> \param p_env ...
393 : !> \param qs_env ...
394 : ! **************************************************************************************************
395 11654 : SUBROUTINE p_env_update_rho(p_env, qs_env)
396 : TYPE(qs_p_env_type), INTENT(IN) :: p_env
397 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
398 :
399 : CHARACTER(LEN=*), PARAMETER :: routineN = 'p_env_update_rho'
400 :
401 : CHARACTER(LEN=default_string_length) :: basis_type
402 : INTEGER :: handle, ispin
403 : TYPE(admm_type), POINTER :: admm_env
404 11654 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao
405 : TYPE(dft_control_type), POINTER :: dft_control
406 : TYPE(mp_para_env_type), POINTER :: para_env
407 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
408 11654 : POINTER :: sab_aux_fit
409 11654 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
410 11654 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux
411 : TYPE(qs_ks_env_type), POINTER :: ks_env
412 : TYPE(task_list_type), POINTER :: task_list
413 :
414 11654 : CALL timeset(routineN, handle)
415 :
416 11654 : CALL get_qs_env(qs_env, dft_control=dft_control)
417 :
418 11654 : IF (dft_control%do_admm) CALL admm_aux_response_density(qs_env, p_env%p1, p_env%p1_admm)
419 :
420 11654 : CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
421 24972 : DO ispin = 1, SIZE(rho1_ao)
422 24972 : CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
423 : END DO
424 :
425 : CALL qs_rho_update_rho(rho_struct=p_env%rho1, &
426 : rho_xc_external=p_env%rho1_xc, &
427 : local_rho_set=p_env%local_rho_set, &
428 11654 : qs_env=qs_env)
429 :
430 11654 : IF (dft_control%do_admm) THEN
431 2426 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
432 1398 : NULLIFY (ks_env, rho1_ao, rho_g_aux, rho_r_aux, task_list)
433 :
434 1398 : CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
435 1398 : basis_type = "AUX_FIT"
436 1398 : CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
437 1398 : IF (admm_env%do_gapw) THEN
438 396 : basis_type = "AUX_FIT_SOFT"
439 396 : task_list => admm_env%admm_gapw_env%task_list
440 : END IF
441 : CALL qs_rho_get(p_env%rho1_admm, &
442 : rho_ao=rho1_ao, &
443 : rho_g=rho_g_aux, &
444 1398 : rho_r=rho_r_aux)
445 2958 : DO ispin = 1, SIZE(rho1_ao)
446 1560 : CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
447 : CALL calculate_rho_elec(ks_env=ks_env, &
448 : matrix_p=rho1_ao(ispin)%matrix, &
449 : rho=rho_r_aux(ispin), &
450 : rho_gspace=rho_g_aux(ispin), &
451 : soft_valid=.FALSE., &
452 : basis_type=basis_type, &
453 2958 : task_list_external=task_list)
454 : END DO
455 1398 : IF (admm_env%do_gapw) THEN
456 396 : CALL get_qs_env(qs_env, para_env=para_env)
457 396 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
458 : CALL calculate_rho_atom_coeff(qs_env, rho1_ao, &
459 : rho_atom_set=p_env%local_rho_set_admm%rho_atom_set, &
460 : qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
461 396 : oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=para_env)
462 : END IF
463 : END IF
464 : END IF
465 :
466 11654 : CALL timestop(handle)
467 :
468 11654 : END SUBROUTINE p_env_update_rho
469 :
470 : ! **************************************************************************************************
471 : !> \brief To be called after the value of psi0 has changed.
472 : !> Recalculates the quantities S_psi0 and m_epsilon.
473 : !> \param p_env the perturbation environment to set
474 : !> \param qs_env ...
475 : !> \par History
476 : !> 07.2002 created [fawzi]
477 : !> \author Fawzi Mohamed
478 : ! **************************************************************************************************
479 1830 : SUBROUTINE p_env_psi0_changed(p_env, qs_env)
480 :
481 : TYPE(qs_p_env_type) :: p_env
482 : TYPE(qs_environment_type), POINTER :: qs_env
483 :
484 : CHARACTER(len=*), PARAMETER :: routineN = 'p_env_psi0_changed'
485 :
486 : INTEGER :: handle, iounit, lfomo, n_spins, nmo, spin
487 : LOGICAL :: was_present
488 : REAL(KIND=dp) :: maxocc
489 1830 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
490 1830 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0
491 : TYPE(cp_fm_type), POINTER :: mo_coeff
492 : TYPE(cp_logger_type), POINTER :: logger
493 1830 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
494 : TYPE(dft_control_type), POINTER :: dft_control
495 1830 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
496 : TYPE(mp_para_env_type), POINTER :: para_env
497 : TYPE(qs_energy_type), POINTER :: energy
498 : TYPE(qs_ks_env_type), POINTER :: ks_env
499 : TYPE(qs_rho_type), POINTER :: rho
500 : TYPE(section_vals_type), POINTER :: input, lr_section
501 :
502 1830 : CALL timeset(routineN, handle)
503 :
504 1830 : NULLIFY (ao_mo_fm_pools, mos, psi0, matrix_s, mos, para_env, ks_env, rho, &
505 1830 : logger, input, lr_section, energy, matrix_ks, dft_control, rho_ao)
506 1830 : logger => cp_get_default_logger()
507 :
508 : CALL get_qs_env(qs_env, &
509 : ks_env=ks_env, &
510 : mos=mos, &
511 : matrix_s=matrix_s, &
512 : matrix_ks=matrix_ks, &
513 : para_env=para_env, &
514 : rho=rho, &
515 : input=input, &
516 : energy=energy, &
517 1830 : dft_control=dft_control)
518 :
519 1830 : CALL qs_rho_get(rho, rho_ao=rho_ao)
520 :
521 1830 : n_spins = dft_control%nspins
522 : CALL mpools_get(qs_env%mpools, &
523 1830 : ao_mo_fm_pools=ao_mo_fm_pools)
524 7632 : ALLOCATE (psi0(n_spins))
525 3972 : DO spin = 1, n_spins
526 2142 : CALL get_mo_set(mos(spin), mo_coeff=mo_coeff)
527 2142 : CALL cp_fm_create(psi0(spin), mo_coeff%matrix_struct)
528 3972 : CALL cp_fm_to_fm(mo_coeff, psi0(spin))
529 : END DO
530 :
531 1830 : lr_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES")
532 : ! def psi0d
533 1830 : IF (p_env%orthogonal_orbitals) THEN
534 1830 : IF (ASSOCIATED(p_env%psi0d)) THEN
535 0 : CALL cp_fm_release(p_env%psi0d)
536 : END IF
537 1830 : p_env%psi0d => psi0
538 : ELSE
539 :
540 0 : DO spin = 1, n_spins
541 : ! m_epsilon=cholesky_decomposition(psi0^T S psi0)^-1
542 : ! could be optimized by combining next two calls
543 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
544 : psi0(spin), &
545 : p_env%S_psi0(spin), &
546 0 : ncol=p_env%n_mo(spin), alpha=1.0_dp)
547 : CALL parallel_gemm(transa='T', transb='N', n=p_env%n_mo(spin), &
548 : m=p_env%n_mo(spin), k=p_env%n_ao(spin), alpha=1.0_dp, &
549 : matrix_a=psi0(spin), &
550 : matrix_b=p_env%S_psi0(spin), &
551 0 : beta=0.0_dp, matrix_c=p_env%m_epsilon(spin))
552 : CALL cp_fm_cholesky_decompose(p_env%m_epsilon(spin), &
553 0 : n=p_env%n_mo(spin))
554 :
555 : ! Smo_inv= (psi0^T S psi0)^-1
556 0 : CALL cp_fm_set_all(p_env%Smo_inv(spin), 0.0_dp, 1.0_dp)
557 : ! faster using cp_fm_cholesky_invert ?
558 : CALL cp_fm_triangular_multiply( &
559 : triangular_matrix=p_env%m_epsilon(spin), &
560 : matrix_b=p_env%Smo_inv(spin), side='R', &
561 : invert_tr=.TRUE., n_rows=p_env%n_mo(spin), &
562 0 : n_cols=p_env%n_mo(spin))
563 : CALL cp_fm_triangular_multiply( &
564 : triangular_matrix=p_env%m_epsilon(spin), &
565 : matrix_b=p_env%Smo_inv(spin), side='R', &
566 : transpose_tr=.TRUE., &
567 : invert_tr=.TRUE., n_rows=p_env%n_mo(spin), &
568 0 : n_cols=p_env%n_mo(spin))
569 :
570 : ! psi0d=psi0 (psi0^T S psi0)^-1
571 : ! faster using cp_fm_cholesky_invert ?
572 : CALL cp_fm_to_fm(psi0(spin), &
573 0 : p_env%psi0d(spin))
574 : CALL cp_fm_triangular_multiply( &
575 : triangular_matrix=p_env%m_epsilon(spin), &
576 : matrix_b=p_env%psi0d(spin), side='R', &
577 : invert_tr=.TRUE., n_rows=p_env%n_ao(spin), &
578 0 : n_cols=p_env%n_mo(spin))
579 : CALL cp_fm_triangular_multiply( &
580 : triangular_matrix=p_env%m_epsilon(spin), &
581 : matrix_b=p_env%psi0d(spin), side='R', &
582 : transpose_tr=.TRUE., &
583 : invert_tr=.TRUE., n_rows=p_env%n_ao(spin), &
584 0 : n_cols=p_env%n_mo(spin))
585 :
586 : ! updates P
587 : CALL get_mo_set(mos(spin), lfomo=lfomo, &
588 0 : nmo=nmo, maxocc=maxocc)
589 0 : IF (lfomo > nmo) THEN
590 0 : CALL dbcsr_set(rho_ao(spin)%matrix, 0.0_dp)
591 : CALL cp_dbcsr_plus_fm_fm_t(rho_ao(spin)%matrix, &
592 : matrix_v=psi0(spin), &
593 : matrix_g=p_env%psi0d(spin), &
594 0 : ncol=p_env%n_mo(spin))
595 0 : CALL dbcsr_scale(rho_ao(spin)%matrix, alpha_scalar=maxocc)
596 : ELSE
597 0 : CPABORT("symmetrized onesided smearing to do")
598 : END IF
599 : END DO
600 :
601 : ! updates rho
602 0 : CALL qs_rho_update_rho(rho_struct=rho, qs_env=qs_env)
603 :
604 : ! tells ks_env that p changed
605 0 : CALL qs_ks_did_change(ks_env=ks_env, rho_changed=.TRUE.)
606 :
607 : END IF
608 :
609 : ! updates K (if necessary)
610 1830 : CALL qs_ks_update_qs_env(qs_env)
611 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
612 1830 : extension=".linresLog")
613 1830 : IF (iounit > 0) THEN
614 890 : CALL section_vals_get(lr_section, explicit=was_present)
615 890 : IF (was_present) THEN
616 : WRITE (UNIT=iounit, FMT="(/,(T3,A,T55,F25.14))") &
617 168 : "Total energy ground state: ", energy%total
618 : END IF
619 : END IF
620 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
621 1830 : "PRINT%PROGRAM_RUN_INFO")
622 : !-----------------------------------------------------------------------|
623 : ! calculates |
624 : ! m_epsilon = - psi0d^T times K times psi0d |
625 : ! = - [K times psi0d]^T times psi0d (because K is symmetric) |
626 : !-----------------------------------------------------------------------|
627 3972 : DO spin = 1, n_spins
628 : ! S_psi0 = k times psi0d
629 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(spin)%matrix, &
630 : p_env%psi0d(spin), &
631 2142 : p_env%S_psi0(spin), p_env%n_mo(spin))
632 : ! m_epsilon = -1 times S_psi0^T times psi0d
633 : CALL parallel_gemm('T', 'N', &
634 : p_env%n_mo(spin), p_env%n_mo(spin), p_env%n_ao(spin), &
635 : -1.0_dp, p_env%S_psi0(spin), p_env%psi0d(spin), &
636 3972 : 0.0_dp, p_env%m_epsilon(spin))
637 : END DO
638 :
639 : !----------------------------------|
640 : ! calculates S_psi0 = S * psi0 |
641 : !----------------------------------|
642 : ! calculating this reduces the mat mult without storing a full aoxao
643 : ! matrix (for P). If nspin>1 you might consider calculating it on the
644 : ! fly to spare some memory
645 1830 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
646 3972 : DO spin = 1, n_spins
647 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
648 : psi0(spin), &
649 : p_env%S_psi0(spin), &
650 3972 : p_env%n_mo(spin))
651 : END DO
652 :
653 : ! releases psi0
654 1830 : IF (p_env%orthogonal_orbitals) THEN
655 1830 : NULLIFY (psi0)
656 : ELSE
657 0 : CALL cp_fm_release(psi0)
658 : END IF
659 :
660 : ! tells kpp1_env about the change of psi0
661 1830 : CALL kpp1_did_change(p_env%kpp1_env)
662 :
663 1830 : CALL timestop(handle)
664 :
665 1830 : END SUBROUTINE p_env_psi0_changed
666 :
667 : ! **************************************************************************************************
668 : !> \brief does a preorthogonalization of the given matrix:
669 : !> v = (I-PS)v
670 : !> \param p_env the perturbation environment
671 : !> \param qs_env the qs_env that is perturbed by this p_env
672 : !> \param v matrix to orthogonalize
673 : !> \param n_cols the number of columns of C to multiply (defaults to size(v,2))
674 : !> \par History
675 : !> 02.09.2002 adapted for new qs_p_env_type (TC)
676 : !> \author Fawzi Mohamed
677 : ! **************************************************************************************************
678 0 : SUBROUTINE p_preortho(p_env, qs_env, v, n_cols)
679 :
680 : TYPE(qs_p_env_type) :: p_env
681 : TYPE(qs_environment_type), POINTER :: qs_env
682 : TYPE(cp_fm_type), DIMENSION(:), INTENT(inout) :: v
683 : INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols
684 :
685 : CHARACTER(len=*), PARAMETER :: routineN = 'p_preortho'
686 :
687 : INTEGER :: cols, handle, max_cols, maxnmo, n_spins, &
688 : nmo2, spin, v_cols, v_rows
689 : TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool
690 : TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct
691 : TYPE(cp_fm_type) :: tmp_matrix
692 : TYPE(dft_control_type), POINTER :: dft_control
693 :
694 0 : CALL timeset(routineN, handle)
695 :
696 0 : NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
697 0 : dft_control)
698 :
699 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
700 0 : CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
701 0 : n_spins = dft_control%nspins
702 0 : maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
703 0 : CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
704 0 : CPASSERT(SIZE(v) >= n_spins)
705 : ! alloc tmp storage
706 0 : IF (PRESENT(n_cols)) THEN
707 0 : max_cols = MAXVAL(n_cols(1:n_spins))
708 : ELSE
709 0 : max_cols = 0
710 0 : DO spin = 1, n_spins
711 0 : CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
712 0 : max_cols = MAX(max_cols, v_cols)
713 : END DO
714 : END IF
715 0 : IF (max_cols <= nmo2) THEN
716 0 : CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix)
717 : ELSE
718 : CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, &
719 0 : ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
720 0 : CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
721 0 : CALL cp_fm_struct_release(tmp_fmstruct)
722 : END IF
723 :
724 0 : DO spin = 1, n_spins
725 :
726 : CALL cp_fm_get_info(v(spin), &
727 0 : nrow_global=v_rows, ncol_global=v_cols)
728 0 : CPASSERT(v_rows >= p_env%n_ao(spin))
729 0 : cols = v_cols
730 0 : IF (PRESENT(n_cols)) THEN
731 0 : CPASSERT(n_cols(spin) <= cols)
732 0 : cols = n_cols(spin)
733 : END IF
734 0 : CPASSERT(cols <= max_cols)
735 :
736 : ! tmp_matrix = v^T (S psi0)
737 : CALL parallel_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), &
738 : k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin), &
739 : matrix_b=p_env%S_psi0(spin), beta=0.0_dp, &
740 0 : matrix_c=tmp_matrix)
741 : ! v = v - psi0d tmp_matrix^T = v - psi0d psi0^T S v
742 : CALL parallel_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, &
743 : k=p_env%n_mo(spin), alpha=-1.0_dp, &
744 : matrix_a=p_env%psi0d(spin), matrix_b=tmp_matrix, &
745 0 : beta=1.0_dp, matrix_c=v(spin))
746 :
747 : END DO
748 :
749 0 : IF (max_cols <= nmo2) THEN
750 0 : CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
751 : ELSE
752 0 : CALL cp_fm_release(tmp_matrix)
753 : END IF
754 :
755 0 : CALL timestop(handle)
756 :
757 0 : END SUBROUTINE p_preortho
758 :
759 : ! **************************************************************************************************
760 : !> \brief does a postorthogonalization on the given matrix vector:
761 : !> v = (I-SP) v
762 : !> \param p_env the perturbation environment
763 : !> \param qs_env the qs_env that is perturbed by this p_env
764 : !> \param v matrix to orthogonalize
765 : !> \param n_cols the number of columns of C to multiply (defaults to size(v,2))
766 : !> \par History
767 : !> 07.2002 created [fawzi]
768 : !> \author Fawzi Mohamed
769 : ! **************************************************************************************************
770 0 : SUBROUTINE p_postortho(p_env, qs_env, v, n_cols)
771 :
772 : TYPE(qs_p_env_type) :: p_env
773 : TYPE(qs_environment_type), POINTER :: qs_env
774 : TYPE(cp_fm_type), DIMENSION(:), INTENT(inout) :: v
775 : INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols
776 :
777 : CHARACTER(len=*), PARAMETER :: routineN = 'p_postortho'
778 :
779 : INTEGER :: cols, handle, max_cols, maxnmo, n_spins, &
780 : nmo2, spin, v_cols, v_rows
781 : TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool
782 : TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct
783 : TYPE(cp_fm_type) :: tmp_matrix
784 : TYPE(dft_control_type), POINTER :: dft_control
785 :
786 0 : CALL timeset(routineN, handle)
787 :
788 0 : NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
789 0 : dft_control)
790 :
791 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
792 0 : CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
793 0 : n_spins = dft_control%nspins
794 0 : maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
795 0 : CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
796 0 : CPASSERT(SIZE(v) >= n_spins)
797 : ! alloc tmp storage
798 0 : IF (PRESENT(n_cols)) THEN
799 0 : max_cols = MAXVAL(n_cols(1:n_spins))
800 : ELSE
801 0 : max_cols = 0
802 0 : DO spin = 1, n_spins
803 0 : CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
804 0 : max_cols = MAX(max_cols, v_cols)
805 : END DO
806 : END IF
807 0 : IF (max_cols <= nmo2) THEN
808 0 : CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix)
809 : ELSE
810 : CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, &
811 0 : ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
812 0 : CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
813 0 : CALL cp_fm_struct_release(tmp_fmstruct)
814 : END IF
815 :
816 0 : DO spin = 1, n_spins
817 :
818 : CALL cp_fm_get_info(v(spin), &
819 0 : nrow_global=v_rows, ncol_global=v_cols)
820 0 : CPASSERT(v_rows >= p_env%n_ao(spin))
821 0 : cols = v_cols
822 0 : IF (PRESENT(n_cols)) THEN
823 0 : CPASSERT(n_cols(spin) <= cols)
824 0 : cols = n_cols(spin)
825 : END IF
826 0 : CPASSERT(cols <= max_cols)
827 :
828 : ! tmp_matrix = v^T psi0d
829 : CALL parallel_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), &
830 : k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin), &
831 : matrix_b=p_env%psi0d(spin), beta=0.0_dp, &
832 0 : matrix_c=tmp_matrix)
833 : ! v = v - (S psi0) tmp_matrix^T = v - S psi0 psi0d^T v
834 : CALL parallel_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, &
835 : k=p_env%n_mo(spin), alpha=-1.0_dp, &
836 : matrix_a=p_env%S_psi0(spin), matrix_b=tmp_matrix, &
837 0 : beta=1.0_dp, matrix_c=v(spin))
838 :
839 : END DO
840 :
841 0 : IF (max_cols <= nmo2) THEN
842 0 : CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
843 : ELSE
844 0 : CALL cp_fm_release(tmp_matrix)
845 : END IF
846 :
847 0 : CALL timestop(handle)
848 :
849 0 : END SUBROUTINE p_postortho
850 :
851 : ! **************************************************************************************************
852 : !> \brief ...
853 : !> \param qs_env ...
854 : !> \param p_env ...
855 : ! **************************************************************************************************
856 2846 : SUBROUTINE p_env_finish_kpp1(qs_env, p_env)
857 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
858 : TYPE(qs_p_env_type), INTENT(IN) :: p_env
859 :
860 : CHARACTER(len=*), PARAMETER :: routineN = 'p_env_finish_kpp1'
861 :
862 : INTEGER :: handle, ispin, nao, nao_aux
863 : TYPE(admm_type), POINTER :: admm_env
864 : TYPE(dbcsr_type) :: work_hmat
865 : TYPE(dft_control_type), POINTER :: dft_control
866 :
867 2846 : CALL timeset(routineN, handle)
868 :
869 2846 : CALL get_qs_env(qs_env, dft_control=dft_control, admm_env=admm_env)
870 :
871 2846 : IF (dft_control%do_admm) THEN
872 2402 : CALL dbcsr_copy(work_hmat, p_env%kpp1(1)%matrix)
873 :
874 2402 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
875 5104 : DO ispin = 1, SIZE(p_env%kpp1)
876 : CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1_admm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
877 2702 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
878 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
879 2702 : admm_env%work_aux_orb, 0.0_dp, admm_env%work_orb_orb)
880 2702 : CALL dbcsr_set(work_hmat, 0.0_dp)
881 2702 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, work_hmat, keep_sparsity=.TRUE.)
882 5104 : CALL dbcsr_add(p_env%kpp1(ispin)%matrix, work_hmat, 1.0_dp, 1.0_dp)
883 : END DO
884 :
885 2402 : CALL dbcsr_release(work_hmat)
886 : END IF
887 :
888 2846 : CALL timestop(handle)
889 :
890 2846 : END SUBROUTINE p_env_finish_kpp1
891 :
892 : END MODULE qs_p_env_methods
|
