(git:0de0cc2)
qs_p_env_methods.F
Go to the documentation of this file.
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 ! **************************************************************************************************
18  USE admm_types, ONLY: admm_gapw_r3d_rs_type,&
19  admm_type,&
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
28  USE cp_fm_basic_linalg, ONLY: cp_fm_symm,&
31  USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
32  cp_fm_pool_type,&
36  fm_pools_create_fm_vect
40  cp_fm_struct_type
41  USE cp_fm_types, ONLY: cp_fm_create,&
43  cp_fm_release,&
45  cp_fm_to_fm,&
46  cp_fm_type
48  cp_logger_type,&
49  cp_to_string
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
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
71  USE pw_env_types, ONLY: pw_env_type
72  USE pw_types, ONLY: pw_c1d_gs_type,&
73  pw_r3d_rs_type
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
81  kpp1_create,&
84  USE qs_ks_types, ONLY: qs_ks_did_change,&
85  qs_ks_env_type
86  USE qs_linres_types, ONLY: linres_control_type
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
94  USE qs_rho0_methods, ONLY: init_rho0
97  USE qs_rho_methods, ONLY: qs_rho_rebuild,&
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
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  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  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
147  TYPE(cp_blacs_env_type), POINTER :: blacs_env
148  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  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  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
155 
156  CALL timeset(routinen, handle)
157  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  blacs_env=blacs_env)
163 
164  n_spins = dft_control%nspins
165 
166  p_env%new_preconditioner = .true.
167 
168  ALLOCATE (p_env%rho1)
169  CALL qs_rho_create(p_env%rho1)
170  ALLOCATE (p_env%rho1_xc)
171  CALL qs_rho_create(p_env%rho1_xc)
172 
173  ALLOCATE (p_env%kpp1_env)
174  CALL kpp1_create(p_env%kpp1_env)
175 
176  IF (PRESENT(p1_option)) THEN
177  p_env%p1 => p1_option
178  ELSE
179  CALL dbcsr_allocate_matrix_set(p_env%p1, n_spins)
180  DO spin = 1, n_spins
181  ALLOCATE (p_env%p1(spin)%matrix)
182  CALL dbcsr_copy(p_env%p1(spin)%matrix, matrix_s(1)%matrix, &
183  name="p_env%p1-"//trim(adjustl(cp_to_string(spin))))
184  CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp)
185  END DO
186  END IF
187 
188  IF (dft_control%do_admm) THEN
189  CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
190  IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
191  ALLOCATE (p_env%rho1_admm)
192  CALL qs_rho_create(p_env%rho1_admm)
193  END IF
194 
195  IF (PRESENT(p1_admm_option)) THEN
196  p_env%p1_admm => p1_admm_option
197  ELSE
198  CALL dbcsr_allocate_matrix_set(p_env%p1_admm, n_spins)
199  DO spin = 1, n_spins
200  ALLOCATE (p_env%p1_admm(spin)%matrix)
201  CALL dbcsr_copy(p_env%p1_admm(spin)%matrix, matrix_s_aux_fit(1)%matrix, &
202  name="p_env%p1_admm-"//trim(adjustl(cp_to_string(spin))))
203  CALL dbcsr_set(p_env%p1_admm(spin)%matrix, 0.0_dp)
204  END DO
205  END IF
206  CALL get_qs_env(qs_env, admm_env=admm_env)
207  IF (admm_env%do_gapw) THEN
208  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
209  admm_gapw_env => admm_env%admm_gapw_env
210  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  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  mo_mo_fm_pools=mo_mo_fm_pools)
218 
219  p_env%n_mo = 0
220  p_env%n_ao = 0
221  DO spin = 1, n_spins
222  CALL get_mo_set(qs_env%mos(spin), mo_coeff=qs_env_c)
223  CALL cp_fm_get_info(qs_env_c, &
224  ncol_global=n_mo, nrow_global=n_ao)
225  p_env%n_mo(spin) = n_mo
226  p_env%n_ao(spin) = n_ao
227  END DO
228 
229  p_env%orthogonal_orbitals = .false.
230  IF (PRESENT(orthogonal_orbitals)) &
231  p_env%orthogonal_orbitals = orthogonal_orbitals
232 
233  CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%S_psi0, &
234  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  name="p_env%m_epsilon")
239 
240  ! alloc Smo_inv
241  IF (.NOT. p_env%orthogonal_orbitals) THEN
242  CALL fm_pools_create_fm_vect(mo_mo_fm_pools, elements=p_env%Smo_inv, &
243  name="p_env%Smo_inv")
244  END IF
245 
246  IF (.NOT. p_env%orthogonal_orbitals) THEN
247  CALL fm_pools_create_fm_vect(ao_mo_fm_pools, &
248  elements=p_env%psi0d, &
249  name="p_env%psi0d")
250  END IF
251 
252  !------------------------------!
253  ! GAPW/GAPW_XC initializations !
254  !------------------------------!
255  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  qs_kind_set=qs_kind_set)
261 
262  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  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  zcore=0.0_dp)
268  CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole)
269  CALL hartree_local_create(p_env%hartree_local)
270  CALL init_coulomb_local(p_env%hartree_local, natom)
271  ELSEIF (dft_control%qs_control%gapw_xc) THEN
272  CALL get_qs_env(qs_env, &
273  atomic_kind_set=atomic_kind_set, &
274  qs_kind_set=qs_kind_set)
275  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  qs_kind_set, dft_control, para_env)
278  END IF
279 
280  !------------------------!
281  ! LINRES initializations !
282  !------------------------!
283  IF (PRESENT(linres_control)) THEN
284 
285  IF (linres_control%preconditioner_type /= ot_precond_none) THEN
286  ! Initialize the preconditioner matrix
287  IF (.NOT. ASSOCIATED(p_env%preconditioner)) THEN
288 
289  ALLOCATE (p_env%preconditioner(n_spins))
290  DO spin = 1, n_spins
291  CALL init_preconditioner(p_env%preconditioner(spin), &
292  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  name="p_env%PS_psi0")
297  END IF
298  END IF
299 
300  END IF
301 
302  CALL timestop(handle)
303 
304  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
326  TYPE(dft_control_type), POINTER :: dft_control
327 
328  CALL timeset(routinen, handle)
329 
330  NULLIFY (dft_control, matrix_s)
331 
332  CALL get_qs_env(qs_env, dft_control=dft_control)
333  gapw_xc = dft_control%qs_control%gapw_xc
334  IF (.NOT. ASSOCIATED(p_env%kpp1)) THEN
335  CALL get_qs_env(qs_env, matrix_s=matrix_s)
336  nspins = dft_control%nspins
337 
338  CALL dbcsr_allocate_matrix_set(p_env%kpp1, nspins)
339  name = "p_env%kpp1-"
340  CALL compress(name, full=.true.)
341  DO ispin = 1, nspins
342  ALLOCATE (p_env%kpp1(ispin)%matrix)
343  CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, matrix_s(1)%matrix, &
344  name=trim(name)//adjustl(cp_to_string(ispin)))
345  CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
346  END DO
347 
348  CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env)
349  IF (gapw_xc) THEN
350  CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env)
351  END IF
352 
353  END IF
354 
355  IF (dft_control%do_admm .AND. .NOT. ASSOCIATED(p_env%kpp1_admm)) THEN
356  CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
357  nspins = dft_control%nspins
358 
359  CALL dbcsr_allocate_matrix_set(p_env%kpp1_admm, nspins)
360  name = "p_env%kpp1_admm-"
361  CALL compress(name, full=.true.)
362  DO ispin = 1, nspins
363  ALLOCATE (p_env%kpp1_admm(ispin)%matrix)
364  CALL dbcsr_copy(p_env%kpp1_admm(ispin)%matrix, matrix_s(1)%matrix, &
365  name=trim(name)//adjustl(cp_to_string(ispin)))
366  CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
367  END DO
368 
369  IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
370  CALL qs_rho_rebuild(p_env%rho1_admm, qs_env=qs_env, admm=.true.)
371  END IF
372 
373  END IF
374 
375  IF (.NOT. ASSOCIATED(p_env%rho1)) THEN
376  CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env)
377  IF (gapw_xc) THEN
378  CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env)
379  END IF
380 
381  IF (dft_control%do_admm) THEN
382  IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
383  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  CALL timestop(handle)
390  END SUBROUTINE p_env_check_i_alloc
391 
392 ! **************************************************************************************************
393 !> \brief ...
394 !> \param p_env ...
395 !> \param qs_env ...
396 ! **************************************************************************************************
397  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  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  POINTER :: sab_aux_fit
411  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux
412  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  CALL timeset(routinen, handle)
417 
418  CALL get_qs_env(qs_env, dft_control=dft_control)
419 
420  IF (dft_control%do_admm) CALL admm_aux_response_density(qs_env, p_env%p1, p_env%p1_admm)
421 
422  CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
423  DO ispin = 1, SIZE(rho1_ao)
424  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  qs_env=qs_env)
431 
432  IF (dft_control%do_admm) THEN
433  IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
434  NULLIFY (ks_env, rho1_ao, rho_g_aux, rho_r_aux, task_list)
435 
436  CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
437  basis_type = "AUX_FIT"
438  CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
439  IF (admm_env%do_gapw) THEN
440  basis_type = "AUX_FIT_SOFT"
441  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  rho_r=rho_r_aux)
447  DO ispin = 1, SIZE(rho1_ao)
448  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  task_list_external=task_list)
456  END DO
457  IF (admm_env%do_gapw) THEN
458  CALL get_qs_env(qs_env, para_env=para_env)
459  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  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  CALL timestop(handle)
469 
470  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  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  TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
492  TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0
493  TYPE(cp_fm_type), POINTER :: mo_coeff
494  TYPE(cp_logger_type), POINTER :: logger
495  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
496  TYPE(dft_control_type), POINTER :: dft_control
497  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  CALL timeset(routinen, handle)
505 
506  NULLIFY (ao_mo_fm_pools, mos, psi0, matrix_s, mos, para_env, ks_env, rho, &
507  logger, input, lr_section, energy, matrix_ks, dft_control, rho_ao)
508  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  dft_control=dft_control)
520 
521  CALL qs_rho_get(rho, rho_ao=rho_ao)
522 
523  n_spins = dft_control%nspins
524  CALL mpools_get(qs_env%mpools, &
525  ao_mo_fm_pools=ao_mo_fm_pools)
526  ALLOCATE (psi0(n_spins))
527  DO spin = 1, n_spins
528  CALL get_mo_set(mos(spin), mo_coeff=mo_coeff)
529  CALL cp_fm_create(psi0(spin), mo_coeff%matrix_struct)
530  CALL cp_fm_to_fm(mo_coeff, psi0(spin))
531  END DO
532 
533  lr_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES")
534  ! def psi0d
535  IF (p_env%orthogonal_orbitals) THEN
536  IF (ASSOCIATED(p_env%psi0d)) THEN
537  CALL cp_fm_release(p_env%psi0d)
538  END IF
539  p_env%psi0d => psi0
540  ELSE
541 
542  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  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  beta=0.0_dp, matrix_c=p_env%m_epsilon(spin))
554  CALL cp_fm_cholesky_decompose(p_env%m_epsilon(spin), &
555  n=p_env%n_mo(spin))
556 
557  ! Smo_inv= (psi0^T S psi0)^-1
558  CALL cp_fm_set_all(p_env%Smo_inv(spin), 0.0_dp, 1.0_dp)
559  ! faster using cp_fm_cholesky_invert ?
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  n_cols=p_env%n_mo(spin))
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  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  p_env%psi0d(spin))
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  n_cols=p_env%n_mo(spin))
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  n_cols=p_env%n_mo(spin))
587 
588  ! updates P
589  CALL get_mo_set(mos(spin), lfomo=lfomo, &
590  nmo=nmo, maxocc=maxocc)
591  IF (lfomo > nmo) THEN
592  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  ncol=p_env%n_mo(spin))
597  CALL dbcsr_scale(rho_ao(spin)%matrix, alpha_scalar=maxocc)
598  ELSE
599  cpabort("symmetrized onesided smearing to do")
600  END IF
601  END DO
602 
603  ! updates rho
604  CALL qs_rho_update_rho(rho_struct=rho, qs_env=qs_env)
605 
606  ! tells ks_env that p changed
607  CALL qs_ks_did_change(ks_env=ks_env, rho_changed=.true.)
608 
609  END IF
610 
611  ! updates K (if necessary)
612  CALL qs_ks_update_qs_env(qs_env)
613  iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
614  extension=".linresLog")
615  IF (iounit > 0) THEN
616  CALL section_vals_get(lr_section, explicit=was_present)
617  IF (was_present) THEN
618  WRITE (unit=iounit, fmt="(/,(T3,A,T55,F25.14))") &
619  "Total energy ground state: ", energy%total
620  END IF
621  END IF
622  CALL cp_print_key_finished_output(iounit, logger, lr_section, &
623  "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  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  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  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  CALL get_qs_env(qs_env, matrix_s=matrix_s)
648  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  p_env%n_mo(spin))
653  END DO
654 
655  ! releases psi0
656  IF (p_env%orthogonal_orbitals) THEN
657  NULLIFY (psi0)
658  ELSE
659  CALL cp_fm_release(psi0)
660  END IF
661 
662  ! tells kpp1_env about the change of psi0
663  CALL kpp1_did_change(p_env%kpp1_env)
664 
665  CALL timestop(handle)
666 
667  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  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  NULLIFY (dft_control)
691  CALL get_qs_env(qs_env, dft_control=dft_control)
692 
693  n_spins = dft_control%nspins
694  DO spin = 1, n_spins
695  CALL p_op_l1_spin(p_env, qs_env, spin, v(spin), &
696  res(spin))
697  END DO
698 
699  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  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  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  CALL timeset(routinen, handle)
733  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  dft_control=dft_control)
740 
741  cpassert(0 < spin)
742  cpassert(spin <= dft_control%nspins)
743 
744  CALL cp_fm_create(tmp, v%matrix_struct)
745 
746  k_p => matrix_ks(spin)%matrix
747  CALL cp_fm_get_info(v, ncol_global=ncol)
748 
749  IF (p_env%orthogonal_orbitals) THEN
750  CALL cp_dbcsr_sm_fm_multiply(k_p, v, res, ncol)
751  ELSE
752  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  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  p_env%m_epsilon(spin), v, 0.0_dp, tmp)
759  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, tmp, &
760  res, p_env%n_mo(spin), alpha=1.0_dp, beta=1.0_dp)
761 
762  CALL cp_fm_release(tmp)
763  CALL timestop(handle)
764 
765  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  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  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  CALL timeset(routinen, handle)
803  NULLIFY (dft_control, rho, rho1_ao)
804 
805  CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
806 
807  gapw = dft_control%qs_control%gapw
808  gapw_xc = dft_control%qs_control%gapw_xc
809  my_alpha = 1.0_dp
810  IF (PRESENT(alpha)) my_alpha = alpha
811  my_beta = 0.0_dp
812  IF (PRESENT(beta)) my_beta = beta
813 
814  CALL p_env_check_i_alloc(p_env, qs_env=qs_env)
815  n_spins = dft_control%nspins
816 
817  CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
818  DO ispin = 1, SIZE(p1)
819  ! hack to avoid crashes in ep regs
820  IF (.NOT. ASSOCIATED(rho1_ao(ispin)%matrix, p1(ispin)%matrix)) THEN
821  CALL dbcsr_copy(rho1_ao(ispin)%matrix, p1(ispin)%matrix)
822  END IF
823  END DO
824 
825  IF (gapw) THEN
826  CALL qs_rho_update_rho(rho_struct=p_env%rho1, qs_env=qs_env, local_rho_set=p_env%local_rho_set)
827  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  local_rho_set=p_env%local_rho_set)
830  ELSE
831  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  rho1=p_env%rho1, rho1_xc=p_env%rho1_xc)
841  END IF
842 
843  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  alpha=my_alpha, beta=my_beta)
848  END DO
849 
850  CALL timestop(handle)
851 
852  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  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  CALL timeset(routinen, handle)
882 
883  NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
884  dft_control)
885 
886  CALL get_qs_env(qs_env, dft_control=dft_control)
887  CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
888  n_spins = dft_control%nspins
889  maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
890  CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
891  cpassert(SIZE(v) >= n_spins)
892  ! alloc tmp storage
893  IF (PRESENT(n_cols)) THEN
894  max_cols = maxval(n_cols(1:n_spins))
895  ELSE
896  max_cols = 0
897  DO spin = 1, n_spins
898  CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
899  max_cols = max(max_cols, v_cols)
900  END DO
901  END IF
902  IF (max_cols <= nmo2) THEN
903  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  ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
907  CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
908  CALL cp_fm_struct_release(tmp_fmstruct)
909  END IF
910 
911  DO spin = 1, n_spins
912 
913  CALL cp_fm_get_info(v(spin), &
914  nrow_global=v_rows, ncol_global=v_cols)
915  cpassert(v_rows >= p_env%n_ao(spin))
916  cols = v_cols
917  IF (PRESENT(n_cols)) THEN
918  cpassert(n_cols(spin) <= cols)
919  cols = n_cols(spin)
920  END IF
921  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  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  beta=1.0_dp, matrix_c=v(spin))
933 
934  END DO
935 
936  IF (max_cols <= nmo2) THEN
937  CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
938  ELSE
939  CALL cp_fm_release(tmp_matrix)
940  END IF
941 
942  CALL timestop(handle)
943 
944  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  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  CALL timeset(routinen, handle)
974 
975  NULLIFY (maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, &
976  dft_control)
977 
978  CALL get_qs_env(qs_env, dft_control=dft_control)
979  CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool)
980  n_spins = dft_control%nspins
981  maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool)
982  CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo)
983  cpassert(SIZE(v) >= n_spins)
984  ! alloc tmp storage
985  IF (PRESENT(n_cols)) THEN
986  max_cols = maxval(n_cols(1:n_spins))
987  ELSE
988  max_cols = 0
989  DO spin = 1, n_spins
990  CALL cp_fm_get_info(v(spin), ncol_global=v_cols)
991  max_cols = max(max_cols, v_cols)
992  END DO
993  END IF
994  IF (max_cols <= nmo2) THEN
995  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  ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct)
999  CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct)
1000  CALL cp_fm_struct_release(tmp_fmstruct)
1001  END IF
1002 
1003  DO spin = 1, n_spins
1004 
1005  CALL cp_fm_get_info(v(spin), &
1006  nrow_global=v_rows, ncol_global=v_cols)
1007  cpassert(v_rows >= p_env%n_ao(spin))
1008  cols = v_cols
1009  IF (PRESENT(n_cols)) THEN
1010  cpassert(n_cols(spin) <= cols)
1011  cols = n_cols(spin)
1012  END IF
1013  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  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  beta=1.0_dp, matrix_c=v(spin))
1025 
1026  END DO
1027 
1028  IF (max_cols <= nmo2) THEN
1029  CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix)
1030  ELSE
1031  CALL cp_fm_release(tmp_matrix)
1032  END IF
1033 
1034  CALL timestop(handle)
1035 
1036  END SUBROUTINE p_postortho
1037 
1038 ! **************************************************************************************************
1039 !> \brief ...
1040 !> \param qs_env ...
1041 !> \param p_env ...
1042 ! **************************************************************************************************
1043  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  CALL timeset(routinen, handle)
1055 
1056  CALL get_qs_env(qs_env, dft_control=dft_control, admm_env=admm_env)
1057 
1058  IF (dft_control%do_admm) THEN
1059  CALL dbcsr_copy(work_hmat, p_env%kpp1(1)%matrix)
1060 
1061  CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
1062  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  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  admm_env%work_aux_orb, 0.0_dp, admm_env%work_orb_orb)
1067  CALL dbcsr_set(work_hmat, 0.0_dp)
1068  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, work_hmat, keep_sparsity=.true.)
1069  CALL dbcsr_add(p_env%kpp1(ispin)%matrix, work_hmat, 1.0_dp, 1.0_dp)
1070  END DO
1071 
1072  CALL dbcsr_release(work_hmat)
1073  END IF
1074 
1075  CALL timestop(handle)
1076 
1077  END SUBROUTINE p_env_finish_kpp1
1078 
1079 END MODULE qs_p_env_methods
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
subroutine, public admm_aux_response_density(qs_env, dm, dm_admm)
Calculate ADMM auxiliary response density.
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
Define the atomic kind types and their sub types.
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
type(cp_fm_struct_type) function, pointer, public fm_pool_get_el_struct(pool)
returns the structure of the elements in this pool
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
Definition: cp_fm_struct.F:409
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
subroutine, public init_coulomb_local(hartree_local, natom)
...
subroutine, public hartree_local_create(hartree_local)
...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public ot_precond_none
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
types of preconditioners
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
container for various plainwaves related things
Definition: pw_env_types.F:14
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
module that builds the second order perturbation kernel kpp1 = delta_rho|_P delta_rho|_P E drho(P1) d...
subroutine, public kpp1_calc_k_p_p1(p_env, qs_env, rho1, rho1_xc)
calculates the k_p_p1 kernel of the perturbation theory
subroutine, public kpp1_create(kpp1_env)
allocates and initializes a kpp1_env
subroutine, public kpp1_calc_k_p_p1_fdiff(qs_env, k_p_p1, rho, rho1, diff)
calcualtes the k_p_p1 kernel of the perturbation theory with finite differences
subroutine, public kpp1_did_change(kpp1_env)
function to advise of changes either in the grids
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition: qs_ks_types.F:919
Type definitiona for linear response calculations.
subroutine, public local_rho_set_create(local_rho_set)
...
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
Define the neighbor list data types and the corresponding functionality.
Utility functions for the perturbation calculations.
subroutine, public p_op_l1(p_env, qs_env, v, res)
Evaluates Fv (S_mo)^-1 - Sv(epsilon) and stores it in res.
subroutine, public p_postortho(p_env, qs_env, v, n_cols)
does a postorthogonalization on the given matrix vector: v = (I-SP) v
subroutine, public p_env_psi0_changed(p_env, qs_env)
To be called after the value of psi0 has changed. Recalculates the quantities S_psi0 and m_epsilon.
subroutine, public p_env_finish_kpp1(qs_env, p_env)
...
subroutine, public p_op_l2(p_env, qs_env, p1, res, alpha, beta)
evaluates res = alpha kpp1(v)*psi0 + beta res with kpp1 evaluated with p=qs_envrhorho_ao,...
subroutine, public p_env_create(p_env, qs_env, p1_option, p1_admm_option, orthogonal_orbitals, linres_control)
allocates and initializes the perturbation environment (no setup)
subroutine, public p_env_update_rho(p_env, qs_env)
...
subroutine, public p_preortho(p_env, qs_env, v, n_cols)
does a preorthogonalization of the given matrix: v = (I-PS)v
subroutine, public p_env_check_i_alloc(p_env, qs_env)
checks that the intenal storage is allocated, and allocs it if needed
basis types for the calculation of the perturbation of density theory.
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
subroutine, public init_rho0(local_rho_set, qs_env, gapw_control, zcore)
...
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
subroutine, public qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
rebuilds rho (if necessary allocating and initializing it)
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Definition: qs_rho_types.F:99
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
types for task lists