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