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