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 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 1636 : 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 1636 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
144 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
145 1636 : 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 1636 : 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 1636 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
152 :
153 1636 : CALL timeset(routineN, handle)
154 1636 : 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 1636 : blacs_env=blacs_env)
160 :
161 1636 : n_spins = dft_control%nspins
162 :
163 1636 : p_env%new_preconditioner = .TRUE.
164 :
165 1636 : ALLOCATE (p_env%rho1)
166 1636 : CALL qs_rho_create(p_env%rho1)
167 1636 : ALLOCATE (p_env%rho1_xc)
168 1636 : CALL qs_rho_create(p_env%rho1_xc)
169 :
170 1636 : ALLOCATE (p_env%kpp1_env)
171 1636 : CALL kpp1_create(p_env%kpp1_env)
172 :
173 1636 : IF (PRESENT(p1_option)) THEN
174 264 : p_env%p1 => p1_option
175 : ELSE
176 1372 : CALL dbcsr_allocate_matrix_set(p_env%p1, n_spins)
177 2940 : DO spin = 1, n_spins
178 1568 : ALLOCATE (p_env%p1(spin)%matrix)
179 : CALL dbcsr_copy(p_env%p1(spin)%matrix, matrix_s(1)%matrix, &
180 1568 : name="p_env%p1-"//TRIM(ADJUSTL(cp_to_string(spin))))
181 2940 : CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp)
182 : END DO
183 : END IF
184 :
185 1636 : IF (dft_control%do_admm) THEN
186 314 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
187 314 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
188 194 : ALLOCATE (p_env%rho1_admm)
189 194 : CALL qs_rho_create(p_env%rho1_admm)
190 : END IF
191 :
192 314 : IF (PRESENT(p1_admm_option)) THEN
193 0 : p_env%p1_admm => p1_admm_option
194 : ELSE
195 314 : CALL dbcsr_allocate_matrix_set(p_env%p1_admm, n_spins)
196 674 : DO spin = 1, n_spins
197 360 : ALLOCATE (p_env%p1_admm(spin)%matrix)
198 : CALL dbcsr_copy(p_env%p1_admm(spin)%matrix, matrix_s_aux_fit(1)%matrix, &
199 360 : name="p_env%p1_admm-"//TRIM(ADJUSTL(cp_to_string(spin))))
200 674 : CALL dbcsr_set(p_env%p1_admm(spin)%matrix, 0.0_dp)
201 : END DO
202 : END IF
203 314 : CALL get_qs_env(qs_env, admm_env=admm_env)
204 314 : IF (admm_env%do_gapw) THEN
205 28 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
206 28 : admm_gapw_env => admm_env%admm_gapw_env
207 28 : 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 28 : 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 1636 : mo_mo_fm_pools=mo_mo_fm_pools)
215 :
216 4908 : p_env%n_mo = 0
217 4908 : p_env%n_ao = 0
218 3554 : DO spin = 1, n_spins
219 1918 : CALL get_mo_set(qs_env%mos(spin), mo_coeff=qs_env_c)
220 : CALL cp_fm_get_info(qs_env_c, &
221 1918 : ncol_global=n_mo, nrow_global=n_ao)
222 1918 : p_env%n_mo(spin) = n_mo
223 3554 : p_env%n_ao(spin) = n_ao
224 : END DO
225 :
226 1636 : p_env%orthogonal_orbitals = .FALSE.
227 1636 : IF (PRESENT(orthogonal_orbitals)) &
228 1636 : p_env%orthogonal_orbitals = orthogonal_orbitals
229 :
230 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%S_psi0, &
231 1636 : 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 1636 : name="p_env%m_epsilon")
236 :
237 : ! alloc Smo_inv
238 1636 : 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 1636 : 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 1636 : 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 192 : qs_kind_set=qs_kind_set)
258 :
259 192 : 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 192 : 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 192 : zcore=0.0_dp)
265 192 : CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole)
266 192 : CALL hartree_local_create(p_env%hartree_local)
267 192 : CALL init_coulomb_local(p_env%hartree_local, natom)
268 1444 : ELSEIF (dft_control%qs_control%gapw_xc) THEN
269 : CALL get_qs_env(qs_env, &
270 : atomic_kind_set=atomic_kind_set, &
271 28 : qs_kind_set=qs_kind_set)
272 28 : 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 28 : qs_kind_set, dft_control, para_env)
275 : END IF
276 :
277 : !------------------------!
278 : ! LINRES initializations !
279 : !------------------------!
280 1636 : IF (PRESENT(linres_control)) THEN
281 :
282 1636 : IF (linres_control%preconditioner_type /= ot_precond_none) THEN
283 : ! Initialize the preconditioner matrix
284 1632 : IF (.NOT. ASSOCIATED(p_env%preconditioner)) THEN
285 :
286 6810 : ALLOCATE (p_env%preconditioner(n_spins))
287 3546 : DO spin = 1, n_spins
288 : CALL init_preconditioner(p_env%preconditioner(spin), &
289 3546 : 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 1632 : name="p_env%PS_psi0")
294 : END IF
295 : END IF
296 :
297 : END IF
298 :
299 1636 : CALL timestop(handle)
300 :
301 1636 : 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 8570 : 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 8570 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
323 : TYPE(dft_control_type), POINTER :: dft_control
324 :
325 8570 : CALL timeset(routineN, handle)
326 :
327 8570 : NULLIFY (dft_control, matrix_s)
328 :
329 8570 : CALL get_qs_env(qs_env, dft_control=dft_control)
330 8570 : gapw_xc = dft_control%qs_control%gapw_xc
331 8570 : IF (.NOT. ASSOCIATED(p_env%kpp1)) THEN
332 1420 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
333 1420 : nspins = dft_control%nspins
334 :
335 1420 : CALL dbcsr_allocate_matrix_set(p_env%kpp1, nspins)
336 1420 : name = "p_env%kpp1-"
337 1420 : CALL compress(name, full=.TRUE.)
338 3044 : DO ispin = 1, nspins
339 1624 : ALLOCATE (p_env%kpp1(ispin)%matrix)
340 : CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, matrix_s(1)%matrix, &
341 1624 : name=TRIM(name)//ADJUSTL(cp_to_string(ispin)))
342 3044 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
343 : END DO
344 :
345 1420 : CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env)
346 1420 : IF (gapw_xc) THEN
347 28 : CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env)
348 : END IF
349 :
350 : END IF
351 :
352 8570 : IF (dft_control%do_admm .AND. .NOT. ASSOCIATED(p_env%kpp1_admm)) THEN
353 314 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
354 314 : nspins = dft_control%nspins
355 :
356 314 : CALL dbcsr_allocate_matrix_set(p_env%kpp1_admm, nspins)
357 314 : name = "p_env%kpp1_admm-"
358 314 : CALL compress(name, full=.TRUE.)
359 674 : DO ispin = 1, nspins
360 360 : ALLOCATE (p_env%kpp1_admm(ispin)%matrix)
361 : CALL dbcsr_copy(p_env%kpp1_admm(ispin)%matrix, matrix_s(1)%matrix, &
362 360 : name=TRIM(name)//ADJUSTL(cp_to_string(ispin)))
363 674 : CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
364 : END DO
365 :
366 314 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
367 194 : CALL qs_rho_rebuild(p_env%rho1_admm, qs_env=qs_env, admm=.TRUE.)
368 : END IF
369 :
370 : END IF
371 :
372 8570 : 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 8570 : CALL timestop(handle)
387 8570 : END SUBROUTINE p_env_check_i_alloc
388 :
389 : ! **************************************************************************************************
390 : !> \brief ...
391 : !> \param p_env ...
392 : !> \param qs_env ...
393 : ! **************************************************************************************************
394 10040 : 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 10040 : 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 10040 : POINTER :: sab_aux_fit
408 10040 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
409 10040 : 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 10040 : CALL timeset(routineN, handle)
414 :
415 10040 : CALL get_qs_env(qs_env, dft_control=dft_control)
416 :
417 10040 : IF (dft_control%do_admm) CALL admm_aux_response_density(qs_env, p_env%p1, p_env%p1_admm)
418 :
419 10040 : CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
420 21486 : DO ispin = 1, SIZE(rho1_ao)
421 21486 : 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 10040 : qs_env=qs_env)
428 :
429 10040 : IF (dft_control%do_admm) THEN
430 2072 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
431 1160 : NULLIFY (ks_env, rho1_ao, rho_g_aux, rho_r_aux, task_list)
432 :
433 1160 : CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
434 1160 : basis_type = "AUX_FIT"
435 1160 : CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
436 1160 : IF (admm_env%do_gapw) THEN
437 160 : basis_type = "AUX_FIT_SOFT"
438 160 : 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 1160 : rho_r=rho_r_aux)
444 2454 : DO ispin = 1, SIZE(rho1_ao)
445 1294 : 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 2454 : task_list_external=task_list)
453 : END DO
454 1160 : IF (admm_env%do_gapw) THEN
455 160 : CALL get_qs_env(qs_env, para_env=para_env)
456 160 : 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 160 : 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 10040 : CALL timestop(handle)
466 :
467 10040 : 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 1636 : 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 1636 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
489 1636 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0
490 : TYPE(cp_fm_type), POINTER :: mo_coeff
491 : TYPE(cp_logger_type), POINTER :: logger
492 1636 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
493 : TYPE(dft_control_type), POINTER :: dft_control
494 1636 : 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 1636 : CALL timeset(routineN, handle)
502 :
503 1636 : NULLIFY (ao_mo_fm_pools, mos, psi0, matrix_s, mos, para_env, ks_env, rho, &
504 1636 : logger, input, lr_section, energy, matrix_ks, dft_control, rho_ao)
505 1636 : 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 1636 : dft_control=dft_control)
517 :
518 1636 : CALL qs_rho_get(rho, rho_ao=rho_ao)
519 :
520 1636 : n_spins = dft_control%nspins
521 : CALL mpools_get(qs_env%mpools, &
522 1636 : ao_mo_fm_pools=ao_mo_fm_pools)
523 6826 : ALLOCATE (psi0(n_spins))
524 3554 : DO spin = 1, n_spins
525 1918 : CALL get_mo_set(mos(spin), mo_coeff=mo_coeff)
526 1918 : CALL cp_fm_create(psi0(spin), mo_coeff%matrix_struct)
527 3554 : CALL cp_fm_to_fm(mo_coeff, psi0(spin))
528 : END DO
529 :
530 1636 : lr_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES")
531 : ! def psi0d
532 1636 : IF (p_env%orthogonal_orbitals) THEN
533 1636 : IF (ASSOCIATED(p_env%psi0d)) THEN
534 0 : CALL cp_fm_release(p_env%psi0d)
535 : END IF
536 1636 : 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 1636 : CALL qs_ks_update_qs_env(qs_env)
610 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
611 1636 : extension=".linresLog")
612 1636 : IF (iounit > 0) THEN
613 793 : CALL section_vals_get(lr_section, explicit=was_present)
614 793 : IF (was_present) THEN
615 : WRITE (UNIT=iounit, FMT="(/,(T3,A,T55,F25.14))") &
616 159 : "Total energy ground state: ", energy%total
617 : END IF
618 : END IF
619 : CALL cp_print_key_finished_output(iounit, logger, lr_section, &
620 1636 : "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 3554 : 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 1918 : 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 3554 : 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 1636 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
645 3554 : 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 3554 : p_env%n_mo(spin))
650 : END DO
651 :
652 : ! releases psi0
653 1636 : IF (p_env%orthogonal_orbitals) THEN
654 1636 : NULLIFY (psi0)
655 : ELSE
656 0 : CALL cp_fm_release(psi0)
657 : END IF
658 :
659 : ! tells kpp1_env about the change of psi0
660 1636 : CALL kpp1_did_change(p_env%kpp1_env)
661 :
662 1636 : CALL timestop(handle)
663 :
664 1636 : END SUBROUTINE p_env_psi0_changed
665 :
666 : ! **************************************************************************************************
667 : !> \brief does a preorthogonalization of the given matrix:
668 : !> v = (I-PS)v
669 : !> \param p_env the perturbation environment
670 : !> \param qs_env the qs_env that is perturbed by this p_env
671 : !> \param v matrix to orthogonalize
672 : !> \param n_cols the number of columns of C to multiply (defaults to size(v,2))
673 : !> \par History
674 : !> 02.09.2002 adapted for new qs_p_env_type (TC)
675 : !> \author Fawzi Mohamed
676 : ! **************************************************************************************************
677 0 : SUBROUTINE p_preortho(p_env, qs_env, v, n_cols)
678 :
679 : TYPE(qs_p_env_type) :: p_env
680 : TYPE(qs_environment_type), POINTER :: qs_env
681 : TYPE(cp_fm_type), DIMENSION(:), INTENT(inout) :: v
682 : INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols
683 :
684 : CHARACTER(len=*), PARAMETER :: routineN = 'p_preortho'
685 :
686 : INTEGER :: cols, handle, max_cols, maxnmo, n_spins, &
687 : nmo2, spin, v_cols, v_rows
688 : TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool
689 : TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct
690 : TYPE(cp_fm_type) :: tmp_matrix
691 : TYPE(dft_control_type), POINTER :: dft_control
692 :
693 0 : CALL timeset(routineN, handle)
694 :
695 0 : NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
696 0 : dft_control)
697 :
698 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
699 0 : CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
700 0 : n_spins = dft_control%nspins
701 0 : maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
702 0 : CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
703 0 : CPASSERT(SIZE(v) >= n_spins)
704 : ! alloc tmp storage
705 0 : IF (PRESENT(n_cols)) THEN
706 0 : max_cols = MAXVAL(n_cols(1:n_spins))
707 : ELSE
708 0 : max_cols = 0
709 0 : DO spin = 1, n_spins
710 0 : CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
711 0 : max_cols = MAX(max_cols, v_cols)
712 : END DO
713 : END IF
714 0 : IF (max_cols <= nmo2) THEN
715 0 : CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix)
716 : ELSE
717 : CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, &
718 0 : ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
719 0 : CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
720 0 : CALL cp_fm_struct_release(tmp_fmstruct)
721 : END IF
722 :
723 0 : DO spin = 1, n_spins
724 :
725 : CALL cp_fm_get_info(v(spin), &
726 0 : nrow_global=v_rows, ncol_global=v_cols)
727 0 : CPASSERT(v_rows >= p_env%n_ao(spin))
728 0 : cols = v_cols
729 0 : IF (PRESENT(n_cols)) THEN
730 0 : CPASSERT(n_cols(spin) <= cols)
731 0 : cols = n_cols(spin)
732 : END IF
733 0 : CPASSERT(cols <= max_cols)
734 :
735 : ! tmp_matrix = v^T (S psi0)
736 : CALL parallel_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), &
737 : k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin), &
738 : matrix_b=p_env%S_psi0(spin), beta=0.0_dp, &
739 0 : matrix_c=tmp_matrix)
740 : ! v = v - psi0d tmp_matrix^T = v - psi0d psi0^T S v
741 : CALL parallel_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, &
742 : k=p_env%n_mo(spin), alpha=-1.0_dp, &
743 : matrix_a=p_env%psi0d(spin), matrix_b=tmp_matrix, &
744 0 : beta=1.0_dp, matrix_c=v(spin))
745 :
746 : END DO
747 :
748 0 : IF (max_cols <= nmo2) THEN
749 0 : CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
750 : ELSE
751 0 : CALL cp_fm_release(tmp_matrix)
752 : END IF
753 :
754 0 : CALL timestop(handle)
755 :
756 0 : END SUBROUTINE p_preortho
757 :
758 : ! **************************************************************************************************
759 : !> \brief does a postorthogonalization on the given matrix vector:
760 : !> v = (I-SP) v
761 : !> \param p_env the perturbation environment
762 : !> \param qs_env the qs_env that is perturbed by this p_env
763 : !> \param v matrix to orthogonalize
764 : !> \param n_cols the number of columns of C to multiply (defaults to size(v,2))
765 : !> \par History
766 : !> 07.2002 created [fawzi]
767 : !> \author Fawzi Mohamed
768 : ! **************************************************************************************************
769 0 : SUBROUTINE p_postortho(p_env, qs_env, v, n_cols)
770 :
771 : TYPE(qs_p_env_type) :: p_env
772 : TYPE(qs_environment_type), POINTER :: qs_env
773 : TYPE(cp_fm_type), DIMENSION(:), INTENT(inout) :: v
774 : INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols
775 :
776 : CHARACTER(len=*), PARAMETER :: routineN = 'p_postortho'
777 :
778 : INTEGER :: cols, handle, max_cols, maxnmo, n_spins, &
779 : nmo2, spin, v_cols, v_rows
780 : TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool
781 : TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct
782 : TYPE(cp_fm_type) :: tmp_matrix
783 : TYPE(dft_control_type), POINTER :: dft_control
784 :
785 0 : CALL timeset(routineN, handle)
786 :
787 0 : NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
788 0 : dft_control)
789 :
790 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
791 0 : CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
792 0 : n_spins = dft_control%nspins
793 0 : maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
794 0 : CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
795 0 : CPASSERT(SIZE(v) >= n_spins)
796 : ! alloc tmp storage
797 0 : IF (PRESENT(n_cols)) THEN
798 0 : max_cols = MAXVAL(n_cols(1:n_spins))
799 : ELSE
800 0 : max_cols = 0
801 0 : DO spin = 1, n_spins
802 0 : CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
803 0 : max_cols = MAX(max_cols, v_cols)
804 : END DO
805 : END IF
806 0 : IF (max_cols <= nmo2) THEN
807 0 : CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix)
808 : ELSE
809 : CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, &
810 0 : ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
811 0 : CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
812 0 : CALL cp_fm_struct_release(tmp_fmstruct)
813 : END IF
814 :
815 0 : DO spin = 1, n_spins
816 :
817 : CALL cp_fm_get_info(v(spin), &
818 0 : nrow_global=v_rows, ncol_global=v_cols)
819 0 : CPASSERT(v_rows >= p_env%n_ao(spin))
820 0 : cols = v_cols
821 0 : IF (PRESENT(n_cols)) THEN
822 0 : CPASSERT(n_cols(spin) <= cols)
823 0 : cols = n_cols(spin)
824 : END IF
825 0 : CPASSERT(cols <= max_cols)
826 :
827 : ! tmp_matrix = v^T psi0d
828 : CALL parallel_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), &
829 : k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin), &
830 : matrix_b=p_env%psi0d(spin), beta=0.0_dp, &
831 0 : matrix_c=tmp_matrix)
832 : ! v = v - (S psi0) tmp_matrix^T = v - S psi0 psi0d^T v
833 : CALL parallel_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, &
834 : k=p_env%n_mo(spin), alpha=-1.0_dp, &
835 : matrix_a=p_env%S_psi0(spin), matrix_b=tmp_matrix, &
836 0 : beta=1.0_dp, matrix_c=v(spin))
837 :
838 : END DO
839 :
840 0 : IF (max_cols <= nmo2) THEN
841 0 : CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
842 : ELSE
843 0 : CALL cp_fm_release(tmp_matrix)
844 : END IF
845 :
846 0 : CALL timestop(handle)
847 :
848 0 : END SUBROUTINE p_postortho
849 :
850 : ! **************************************************************************************************
851 : !> \brief ...
852 : !> \param qs_env ...
853 : !> \param p_env ...
854 : ! **************************************************************************************************
855 2492 : SUBROUTINE p_env_finish_kpp1(qs_env, p_env)
856 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
857 : TYPE(qs_p_env_type), INTENT(IN) :: p_env
858 :
859 : CHARACTER(len=*), PARAMETER :: routineN = 'p_env_finish_kpp1'
860 :
861 : INTEGER :: handle, ispin, nao, nao_aux
862 : TYPE(admm_type), POINTER :: admm_env
863 : TYPE(dbcsr_type) :: work_hmat
864 : TYPE(dft_control_type), POINTER :: dft_control
865 :
866 2492 : CALL timeset(routineN, handle)
867 :
868 2492 : CALL get_qs_env(qs_env, dft_control=dft_control, admm_env=admm_env)
869 :
870 2492 : IF (dft_control%do_admm) THEN
871 2048 : CALL dbcsr_copy(work_hmat, p_env%kpp1(1)%matrix)
872 :
873 2048 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
874 4368 : DO ispin = 1, SIZE(p_env%kpp1)
875 : CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1_admm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
876 2320 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
877 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
878 2320 : admm_env%work_aux_orb, 0.0_dp, admm_env%work_orb_orb)
879 2320 : CALL dbcsr_set(work_hmat, 0.0_dp)
880 2320 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, work_hmat, keep_sparsity=.TRUE.)
881 4368 : CALL dbcsr_add(p_env%kpp1(ispin)%matrix, work_hmat, 1.0_dp, 1.0_dp)
882 : END DO
883 :
884 2048 : CALL dbcsr_release(work_hmat)
885 : END IF
886 :
887 2492 : CALL timestop(handle)
888 :
889 2492 : END SUBROUTINE p_env_finish_kpp1
890 :
891 : END MODULE qs_p_env_methods
|
