(git:0de0cc2)
energy_corrections.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 Routines for an energy correction on top of a Kohn-Sham calculation
10 !> \par History
11 !> 03.2014 created
12 !> 09.2019 Moved from KG to Kohn-Sham
13 !> 08.2022 Add Density-Corrected DFT methods
14 !> 04.2023 Add hybrid functionals for DC-DFT
15 !> \author JGH
16 ! **************************************************************************************************
20  USE atomic_kind_types, ONLY: atomic_kind_type,&
23  gto_basis_set_type
24  USE bibliography, ONLY: belleflamme2023,&
25  cite_reference
26  USE cell_types, ONLY: cell_type,&
27  pbc
28  USE core_ae, ONLY: build_core_ae
29  USE core_ppl, ONLY: build_core_ppl
30  USE core_ppnl, ONLY: build_core_ppnl
31  USE cp_blacs_env, ONLY: cp_blacs_env_type
32  USE cp_control_types, ONLY: dft_control_type
45  cp_fm_struct_type
46  USE cp_fm_types, ONLY: cp_fm_create,&
48  cp_fm_release,&
51  cp_fm_type
54  cp_logger_type
55  USE cp_output_handling, ONLY: cp_p_file,&
60  put_results
61  USE cp_result_types, ONLY: cp_result_type
62  USE dbcsr_api, ONLY: &
63  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, &
64  dbcsr_dot, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
65  dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
66  dbcsr_type_symmetric
67  USE distribution_1d_types, ONLY: distribution_1d_type
68  USE distribution_2d_types, ONLY: distribution_2d_type
69  USE dm_ls_scf, ONLY: calculate_w_matrix_ls,&
76  USE dm_ls_scf_qs, ONLY: matrix_ls_create,&
79  USE dm_ls_scf_types, ONLY: ls_scf_env_type
82  USE ec_env_types, ONLY: energy_correction_type
83  USE ec_methods, ONLY: create_kernel,&
85  USE external_potential_types, ONLY: get_potential,&
86  gth_potential_type,&
87  sgp_potential_type
88  USE hfx_exx, ONLY: add_exx_to_rhs,&
90  USE input_constants, ONLY: &
98  section_vals_type,&
100  USE kinds, ONLY: default_path_length,&
102  dp
103  USE mao_basis, ONLY: mao_generate_basis
104  USE mathlib, ONLY: det_3x3
105  USE message_passing, ONLY: mp_para_env_type
106  USE molecule_types, ONLY: molecule_type
108  USE particle_types, ONLY: particle_type
109  USE periodic_table, ONLY: ptable
110  USE physcon, ONLY: bohr,&
111  debye,&
112  pascal
116  preconditioner_type
117  USE pw_env_types, ONLY: pw_env_get,&
118  pw_env_type
119  USE pw_grid_types, ONLY: pw_grid_type
120  USE pw_methods, ONLY: pw_axpy,&
121  pw_copy,&
122  pw_integral_ab,&
123  pw_scale,&
124  pw_transfer,&
125  pw_zero
126  USE pw_poisson_methods, ONLY: pw_poisson_solve
127  USE pw_poisson_types, ONLY: pw_poisson_type
128  USE pw_pool_types, ONLY: pw_pool_p_type,&
129  pw_pool_type
130  USE pw_types, ONLY: pw_c1d_gs_type,&
131  pw_r3d_rs_type
135  calculate_ptrace
136  USE qs_density_matrices, ONLY: calculate_density_matrix,&
137  calculate_w_matrix
139  USE qs_dispersion_types, ONLY: qs_dispersion_type
140  USE qs_energy_types, ONLY: qs_energy_type
141  USE qs_environment_types, ONLY: get_qs_env,&
142  qs_environment_type
143  USE qs_force_types, ONLY: qs_force_type,&
146  USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
147  integrate_v_rspace
148  USE qs_kind_types, ONLY: get_qs_kind,&
150  qs_kind_type
154  USE qs_ks_types, ONLY: qs_ks_env_type
155  USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,&
157  USE qs_mo_types, ONLY: deallocate_mo_set,&
158  get_mo_set,&
159  mo_set_type
161  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
162  USE qs_neighbor_lists, ONLY: atom2d_build,&
165  local_atoms_type,&
169  USE qs_rho_types, ONLY: qs_rho_get,&
170  qs_rho_type
171  USE qs_vxc, ONLY: qs_vxc_create
175  USE string_utilities, ONLY: uppercase
182  USE virial_types, ONLY: symmetrize_virial,&
183  virial_type
185 #include "./base/base_uses.f90"
186 
187  IMPLICIT NONE
188 
189  PRIVATE
190 
191 ! *** Global parameters ***
192 
193  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'energy_corrections'
194 
195  PUBLIC :: energy_correction
196 
197 CONTAINS
198 
199 ! **************************************************************************************************
200 !> \brief Energy Correction to a Kohn-Sham simulation
201 !> Available energy corrections: (1) Harris energy functional
202 !> (2) Density-corrected DFT
203 !>
204 !> \param qs_env ...
205 !> \param ec_init ...
206 !> \param calculate_forces ...
207 !> \par History
208 !> 03.2014 created
209 !> \author JGH
210 ! **************************************************************************************************
211  SUBROUTINE energy_correction(qs_env, ec_init, calculate_forces)
212  TYPE(qs_environment_type), POINTER :: qs_env
213  LOGICAL, INTENT(IN), OPTIONAL :: ec_init, calculate_forces
214 
215  CHARACTER(len=*), PARAMETER :: routinen = 'energy_correction'
216 
217  INTEGER :: handle, unit_nr
218  LOGICAL :: my_calc_forces
219  TYPE(cp_logger_type), POINTER :: logger
220  TYPE(energy_correction_type), POINTER :: ec_env
221  TYPE(qs_energy_type), POINTER :: energy
222  TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force
223 
224  CALL timeset(routinen, handle)
225 
226  logger => cp_get_default_logger()
227  IF (logger%para_env%is_source()) THEN
228  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
229  ELSE
230  unit_nr = -1
231  END IF
232 
233  CALL cite_reference(belleflamme2023)
234 
235  NULLIFY (ec_env)
236  CALL get_qs_env(qs_env, ec_env=ec_env)
237 
238  ! Skip energy correction if ground-state is NOT converged
239  IF (.NOT. ec_env%do_skip) THEN
240 
241  ec_env%should_update = .true.
242  IF (PRESENT(ec_init)) ec_env%should_update = ec_init
243 
244  my_calc_forces = .false.
245  IF (PRESENT(calculate_forces)) my_calc_forces = calculate_forces
246 
247  IF (ec_env%should_update) THEN
248  ec_env%old_etotal = 0.0_dp
249  ec_env%etotal = 0.0_dp
250  ec_env%eband = 0.0_dp
251  ec_env%ehartree = 0.0_dp
252  ec_env%ex = 0.0_dp
253  ec_env%exc = 0.0_dp
254  ec_env%vhxc = 0.0_dp
255  ec_env%edispersion = 0.0_dp
256  ec_env%exc_aux_fit = 0.0_dp
257 
258  ! Save total energy of reference calculation
259  CALL get_qs_env(qs_env, energy=energy)
260  ec_env%old_etotal = energy%total
261 
262  END IF
263 
264  IF (my_calc_forces) THEN
265  IF (unit_nr > 0) THEN
266  WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 25), &
267  " Energy Correction Forces ", repeat("-", 26), "!"
268  END IF
269  CALL get_qs_env(qs_env, force=ks_force)
270  CALL zero_qs_force(ks_force)
271  ELSE
272  IF (unit_nr > 0) THEN
273  WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 29), &
274  " Energy Correction ", repeat("-", 29), "!"
275  END IF
276  END IF
277 
278  ! Perform the energy correction
279  CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr)
280 
281  ! Update total energy in qs environment and amount fo correction
282  IF (ec_env%should_update) THEN
283  energy%nonscf_correction = ec_env%etotal - ec_env%old_etotal
284  energy%total = ec_env%etotal
285  END IF
286 
287  IF (.NOT. my_calc_forces .AND. unit_nr > 0) THEN
288  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Energy Correction ", energy%nonscf_correction
289  END IF
290  IF (unit_nr > 0) THEN
291  WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
292  END IF
293 
294  ELSE
295 
296  ! Ground-state energy calculation did not converge,
297  ! do not calculate energy correction
298  IF (unit_nr > 0) THEN
299  WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
300  WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 26), &
301  " Skip Energy Correction ", repeat("-", 27), "!"
302  WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
303  END IF
304 
305  END IF
306 
307  CALL timestop(handle)
308 
309  END SUBROUTINE energy_correction
310 
311 ! **************************************************************************************************
312 !> \brief Energy Correction to a Kohn-Sham simulation
313 !>
314 !> \param qs_env ...
315 !> \param ec_env ...
316 !> \param calculate_forces ...
317 !> \param unit_nr ...
318 !> \par History
319 !> 03.2014 created
320 !> \author JGH
321 ! **************************************************************************************************
322  SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr)
323  TYPE(qs_environment_type), POINTER :: qs_env
324  TYPE(energy_correction_type), POINTER :: ec_env
325  LOGICAL, INTENT(IN) :: calculate_forces
326  INTEGER, INTENT(IN) :: unit_nr
327 
328  INTEGER :: ispin, nspins
329  LOGICAL :: debug_f
330  REAL(kind=dp) :: exc
331  TYPE(dft_control_type), POINTER :: dft_control
332  TYPE(pw_env_type), POINTER :: pw_env
333  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
334 
335  IF (ec_env%should_update) THEN
336  CALL ec_build_neighborlist(qs_env, ec_env)
337  CALL ks_ref_potential(qs_env, &
338  ec_env%vh_rspace, &
339  ec_env%vxc_rspace, &
340  ec_env%vtau_rspace, &
341  ec_env%vadmm_rspace, &
342  ec_env%ehartree, exc)
343 
344  SELECT CASE (ec_env%energy_functional)
345  CASE (ec_functional_harris)
346 
347  CALL ec_build_core_hamiltonian(qs_env, ec_env)
348  CALL ec_build_ks_matrix(qs_env, ec_env)
349 
350  IF (ec_env%mao) THEN
351  ! MAO basis
352  IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef)
353  NULLIFY (ec_env%mao_coef)
354  CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.true., &
355  max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, &
356  eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr)
357  END IF
358 
359  CALL ec_ks_solver(qs_env, ec_env)
360 
361  CALL evaluate_ec_core_matrix_traces(qs_env, ec_env)
362 
363  CASE (ec_functional_dc)
364 
365  ! Prepare Density-corrected DFT (DC-DFT) calculation
366  CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.false.)
367 
368  ! Rebuild KS matrix with DC-DFT XC functional evaluated in ground-state density.
369  ! KS matrix might contain unwanted contributions
370  ! Calculate Hartree and XC related energies here
371  CALL ec_build_ks_matrix(qs_env, ec_env)
372 
373  CASE DEFAULT
374  cpabort("unknown energy correction")
375  END SELECT
376 
377  ! dispersion through pairpotentials
378  CALL ec_disp(qs_env, ec_env, calculate_forces=.false.)
379 
380  ! Calculate total energy
381  CALL ec_energy(ec_env, unit_nr)
382 
383  END IF
384 
385  IF (calculate_forces) THEN
386 
387  debug_f = ec_env%debug_forces .OR. ec_env%debug_stress
388 
389  CALL ec_disp(qs_env, ec_env, calculate_forces=.true.)
390 
391  SELECT CASE (ec_env%energy_functional)
392  CASE (ec_functional_harris)
393 
394  CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
395  ec_env%matrix_p, &
396  ec_env%matrix_s, &
397  ec_env%matrix_w)
398  CALL ec_build_ks_matrix_force(qs_env, ec_env)
399  CASE (ec_functional_dc)
400 
401  ! Prepare Density-corrected DFT (DC-DFT) calculation
402  ! by getting ground-state matrices
403  CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.true.)
404 
405  CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
406  ec_env%matrix_p, &
407  ec_env%matrix_s, &
408  ec_env%matrix_w)
409  CALL ec_dc_build_ks_matrix_force(qs_env, ec_env)
410 
411  CASE DEFAULT
412  cpabort("unknown energy correction")
413  END SELECT
414 
415  CALL response_calculation(qs_env, ec_env)
416 
417  ! Allocate response density on real space grid for use in properties
418  ! Calculated in response_force
419  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
420  nspins = dft_control%nspins
421 
422  cpassert(ASSOCIATED(pw_env))
423  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
424  ALLOCATE (ec_env%rhoz_r(nspins))
425  DO ispin = 1, nspins
426  CALL auxbas_pw_pool%create_pw(ec_env%rhoz_r(ispin))
427  END DO
428 
429  CALL response_force(qs_env, &
430  vh_rspace=ec_env%vh_rspace, &
431  vxc_rspace=ec_env%vxc_rspace, &
432  vtau_rspace=ec_env%vtau_rspace, &
433  vadmm_rspace=ec_env%vadmm_rspace, &
434  matrix_hz=ec_env%matrix_hz, &
435  matrix_pz=ec_env%matrix_z, &
436  matrix_pz_admm=ec_env%z_admm, &
437  matrix_wz=ec_env%matrix_wz, &
438  rhopz_r=ec_env%rhoz_r, &
439  zehartree=ec_env%ehartree, &
440  zexc=ec_env%exc, &
441  zexc_aux_fit=ec_env%exc_aux_fit, &
442  p_env=ec_env%p_env, &
443  debug=debug_f)
444 
445  CALL ec_properties(qs_env, ec_env)
446 
447  ! Deallocate Harris density and response density on grid
448  DO ispin = 1, nspins
449  CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin))
450  CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin))
451  END DO
452  DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r)
453 
454  ! Deallocate matrices
455  IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
456  IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
457  IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
458  IF (ASSOCIATED(ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_t)
459  IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
460  IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
461  IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
462  IF (ASSOCIATED(ec_env%matrix_wz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_wz)
463  IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z)
464 
465  END IF
466 
467  END SUBROUTINE energy_correction_low
468 
469 ! **************************************************************************************************
470 !> \brief Calculates the traces of the core matrices and the density matrix.
471 !> \param qs_env ...
472 !> \param ec_env ...
473 !> \author Ole Schuett
474 !> adapted for energy correction fbelle
475 ! **************************************************************************************************
476  SUBROUTINE evaluate_ec_core_matrix_traces(qs_env, ec_env)
477  TYPE(qs_environment_type), POINTER :: qs_env
478  TYPE(energy_correction_type), POINTER :: ec_env
479 
480  CHARACTER(LEN=*), PARAMETER :: routinen = 'evaluate_ec_core_matrix_traces'
481 
482  INTEGER :: handle
483  TYPE(dft_control_type), POINTER :: dft_control
484  TYPE(qs_energy_type), POINTER :: energy
485 
486  CALL timeset(routinen, handle)
487  NULLIFY (energy)
488 
489  CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
490 
491  ! Core hamiltonian energy
492  CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, energy%core, dft_control%nspins)
493 
494  ! kinetic energy
495  CALL calculate_ptrace(ec_env%matrix_t, ec_env%matrix_p, energy%kinetic, dft_control%nspins)
496 
497  CALL timestop(handle)
498 
499  END SUBROUTINE evaluate_ec_core_matrix_traces
500 
501 ! **************************************************************************************************
502 !> \brief Prepare DC-DFT calculation by copying unaffected ground-state matrices (core Hamiltonian,
503 !> density matrix) into energy correction environment and rebuild the overlap matrix
504 !>
505 !> \param qs_env ...
506 !> \param ec_env ...
507 !> \param calculate_forces ...
508 !> \par History
509 !> 07.2022 created
510 !> \author fbelle
511 ! **************************************************************************************************
512  SUBROUTINE ec_dc_energy(qs_env, ec_env, calculate_forces)
513  TYPE(qs_environment_type), POINTER :: qs_env
514  TYPE(energy_correction_type), POINTER :: ec_env
515  LOGICAL, INTENT(IN) :: calculate_forces
516 
517  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_dc_energy'
518 
519  CHARACTER(LEN=default_string_length) :: headline
520  INTEGER :: handle, ispin, nspins
521  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
522  TYPE(dft_control_type), POINTER :: dft_control
523  TYPE(qs_ks_env_type), POINTER :: ks_env
524  TYPE(qs_rho_type), POINTER :: rho
525 
526  CALL timeset(routinen, handle)
527 
528  NULLIFY (dft_control, ks_env, matrix_h, matrix_p, matrix_s, matrix_w, rho)
529  CALL get_qs_env(qs_env=qs_env, &
530  dft_control=dft_control, &
531  ks_env=ks_env, &
532  matrix_h_kp=matrix_h, &
533  matrix_s_kp=matrix_s, &
534  matrix_w_kp=matrix_w, &
535  rho=rho)
536  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
537  nspins = dft_control%nspins
538 
539  ! For density-corrected DFT only the ground-state matrices are required
540  ! Comply with ec_env environment for property calculations later
541  CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
542  matrix_name="OVERLAP MATRIX", &
543  basis_type_a="HARRIS", &
544  basis_type_b="HARRIS", &
545  sab_nl=ec_env%sab_orb)
546 
547  ! Core Hamiltonian matrix
548  IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
549  CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
550  headline = "CORE HAMILTONIAN MATRIX"
551  ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
552  CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=trim(headline), &
553  template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
554  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, ec_env%sab_orb)
555  CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
556 
557  ! Density matrix
558  IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
559  CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
560  headline = "DENSITY MATRIX"
561  DO ispin = 1, nspins
562  ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
563  CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
564  template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
565  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
566  CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
567  END DO
568 
569  IF (calculate_forces) THEN
570 
571  ! Energy-weighted density matrix
572  IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
573  CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
574  headline = "ENERGY-WEIGHTED DENSITY MATRIX"
575  DO ispin = 1, nspins
576  ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
577  CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
578  template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
579  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
580  CALL dbcsr_copy(ec_env%matrix_w(ispin, 1)%matrix, matrix_w(ispin, 1)%matrix)
581  END DO
582 
583  END IF
584 
585  ! External field (nonperiodic case)
586  ec_env%efield_nuclear = 0.0_dp
587  ec_env%efield_elec = 0.0_dp
588  CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces=.false.)
589 
590  CALL timestop(handle)
591 
592  END SUBROUTINE ec_dc_energy
593 
594 ! **************************************************************************************************
595 !> \brief Kohn-Sham matrix contributions to force in DC-DFT
596 !> also calculate right-hand-side matrix B for response equations AX=B
597 !> \param qs_env ...
598 !> \param ec_env ...
599 !> \par History
600 !> 08.2022 adapted from qs_ks_build_kohn_sham_matrix
601 !> \author fbelle
602 ! **************************************************************************************************
603  SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
604  TYPE(qs_environment_type), POINTER :: qs_env
605  TYPE(energy_correction_type), POINTER :: ec_env
606 
607  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_dc_build_ks_matrix_force'
608 
609  INTEGER :: handle, i, iounit, ispin, natom, nspins
610  LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
611  use_virial
612  REAL(dp) :: dummy_real, dummy_real2(2), ehartree, &
613  eovrl, exc, fconv
614  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
615  REAL(dp), DIMENSION(3) :: fodeb, fodeb2
616  REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
617  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
618  TYPE(cell_type), POINTER :: cell
619  TYPE(cp_logger_type), POINTER :: logger
620  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, scrm
621  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
622  TYPE(dft_control_type), POINTER :: dft_control
623  TYPE(mp_para_env_type), POINTER :: para_env
624  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
625  POINTER :: sab_orb
626  TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
627  TYPE(pw_env_type), POINTER :: pw_env
628  TYPE(pw_grid_type), POINTER :: pw_grid
629  TYPE(pw_poisson_type), POINTER :: poisson_env
630  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
631  TYPE(pw_r3d_rs_type) :: v_hartree_rspace
632  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace, v_rspace_in, &
633  v_tau_rspace
634  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
635  TYPE(qs_ks_env_type), POINTER :: ks_env
636  TYPE(qs_rho_type), POINTER :: rho
637  TYPE(section_vals_type), POINTER :: ec_hfx_sections
638  TYPE(virial_type), POINTER :: virial
639 
640  CALL timeset(routinen, handle)
641 
642  debug_forces = ec_env%debug_forces
643  debug_stress = ec_env%debug_stress
644 
645  logger => cp_get_default_logger()
646  IF (logger%para_env%is_source()) THEN
647  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
648  ELSE
649  iounit = -1
650  END IF
651 
652  NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, matrix_ks, &
653  matrix_p, matrix_s, para_env, pw_env, rho, sab_orb, virial)
654  CALL get_qs_env(qs_env=qs_env, &
655  cell=cell, &
656  dft_control=dft_control, &
657  force=force, &
658  ks_env=ks_env, &
659  matrix_ks=matrix_ks, &
660  matrix_s=matrix_s, &
661  para_env=para_env, &
662  pw_env=pw_env, &
663  rho=rho, &
664  sab_orb=sab_orb, &
665  virial=virial)
666  cpassert(ASSOCIATED(pw_env))
667 
668  nspins = dft_control%nspins
669  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
670 
671  fconv = 1.0e-9_dp*pascal/cell%deth
672  IF (debug_stress .AND. use_virial) THEN
673  sttot = virial%pv_virial
674  END IF
675 
676  ! Get density matrix of reference calculation
677  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
678 
679  NULLIFY (auxbas_pw_pool, poisson_env)
680  ! gets the tmp grids
681  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
682  poisson_env=poisson_env)
683 
684  ! Calculate the Hartree potential
685  CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
686  CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
687  CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
688 
689  ! Get the total input density in g-space [ions + electrons]
690  CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
691 
692  ! v_H[n_in]
693  IF (use_virial) THEN
694 
695  ! Stress tensor - Volume and Green function contribution
696  h_stress(:, :) = 0.0_dp
697  CALL pw_poisson_solve(poisson_env, &
698  density=rho_tot_gspace, &
699  ehartree=ehartree, &
700  vhartree=v_hartree_gspace, &
701  h_stress=h_stress)
702 
703  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
704  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
705 
706  IF (debug_stress) THEN
707  stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
708  CALL para_env%sum(stdeb)
709  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
710  'STRESS| GREEN 1st V_H[n_in]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
711  END IF
712 
713  ELSE
714  CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
715  v_hartree_gspace)
716  END IF
717 
718  CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
719  CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
720 
721  ! Save density on real space grid for use in properties
722  CALL qs_rho_get(rho, rho_r=rho_r)
723  ALLOCATE (ec_env%rhoout_r(nspins))
724  DO ispin = 1, nspins
725  CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
726  CALL pw_copy(rho_r(ispin), ec_env%rhoout_r(ispin))
727  END DO
728 
729  ! Getting nuclear force contribution from the core charge density
730  ! Vh(rho_c + rho_in)
731  IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
732  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
733  CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
734  IF (debug_forces) THEN
735  fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
736  CALL para_env%sum(fodeb)
737  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb
738  END IF
739  IF (debug_stress .AND. use_virial) THEN
740  stdeb = fconv*(virial%pv_ehartree - stdeb)
741  CALL para_env%sum(stdeb)
742  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
743  'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb)
744  END IF
745 
746  ! v_XC[n_in]_DC
747  ! v_rspace and v_tau_rspace are generated from the auxbas pool
748  NULLIFY (v_rspace, v_tau_rspace)
749 
750  ! only activate stress calculation if
751  IF (use_virial) virial%pv_calculate = .true.
752 
753  ! Exchange-correlation potential
754  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
755  vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
756 
757  IF (.NOT. ASSOCIATED(v_rspace)) THEN
758  ALLOCATE (v_rspace(nspins))
759  DO ispin = 1, nspins
760  CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
761  CALL pw_zero(v_rspace(ispin))
762  END DO
763  END IF
764 
765  IF (use_virial) THEN
766  virial%pv_exc = virial%pv_exc - virial%pv_xc
767  virial%pv_virial = virial%pv_virial - virial%pv_xc
768  ! virial%pv_xc will be zeroed in the xc routines
769  END IF
770 
771  ! initialize srcm matrix
772  NULLIFY (scrm)
773  CALL dbcsr_allocate_matrix_set(scrm, nspins)
774  DO ispin = 1, nspins
775  ALLOCATE (scrm(ispin)%matrix)
776  CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
777  CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
778  CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
779  END DO
780 
781  pw_grid => v_hartree_rspace%pw_grid
782  ALLOCATE (v_rspace_in(nspins))
783  DO ispin = 1, nspins
784  CALL v_rspace_in(ispin)%create(pw_grid)
785  END DO
786 
787  ! v_rspace_in = v_H[n_in] + v_xc[n_in] calculated in ks_ref_potential
788  DO ispin = 1, nspins
789  ! v_xc[n_in]_GS
790  CALL pw_transfer(ec_env%vxc_rspace(ispin), v_rspace_in(ispin))
791  ! add v_H[n_in]
792  CALL pw_axpy(ec_env%vh_rspace, v_rspace_in(ispin))
793  END DO
794 
795 !------------------------------------------------
796 
797  ! If hybrid functional in DC-DFT
798  ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF")
799  CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
800 
801  IF (do_ec_hfx) THEN
802 
803  IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
804  IF (debug_forces) fodeb2(1:3) = force(1)%overlap_admm(1:3, 1)
805 
806  ! Calculate direct HFX forces here
807  ! Virial contribution (fock_4c) done inside calculate_exx
808  dummy_real = 0.0_dp
809  CALL calculate_exx(qs_env=qs_env, &
810  unit_nr=iounit, &
811  hfx_sections=ec_hfx_sections, &
812  x_data=ec_env%x_data, &
813  do_gw=.false., &
814  do_admm=ec_env%do_ec_admm, &
815  calc_forces=.true., &
816  reuse_hfx=ec_env%reuse_hfx, &
817  do_im_time=.false., &
818  e_ex_from_gw=dummy_real, &
819  e_admm_from_gw=dummy_real2, &
820  t3=dummy_real)
821 
822  IF (debug_forces) THEN
823  fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
824  CALL para_env%sum(fodeb)
825  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC ", fodeb
826 
827  fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb2(1:3)
828  CALL para_env%sum(fodeb2)
829  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC*S ", fodeb2
830  END IF
831  IF (debug_stress .AND. use_virial) THEN
832  stdeb = -1.0_dp*fconv*virial%pv_fock_4c
833  CALL para_env%sum(stdeb)
834  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
835  'STRESS| P*hfx_DC ', one_third_sum_diag(stdeb), det_3x3(stdeb)
836  END IF
837 
838  END IF
839 
840 !------------------------------------------------
841 
842  ! Stress-tensor contribution derivative of integrand
843  ! int v_Hxc[n^în]*n^out
844  IF (use_virial) THEN
845  pv_loc = virial%pv_virial
846  END IF
847 
848  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
849  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
850 
851  DO ispin = 1, nspins
852  ! Add v_H[n_in] + v_xc[n_in] = v_rspace
853  CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
854  CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
855  ! integrate over potential <a|V|b>
856  CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
857  hmat=scrm(ispin), &
858  pmat=matrix_p(ispin, 1), &
859  qs_env=qs_env, &
860  calculate_forces=.true., &
861  basis_type="HARRIS", &
862  task_list_external=ec_env%task_list)
863  END DO
864 
865  IF (debug_forces) THEN
866  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
867  CALL para_env%sum(fodeb)
868  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb
869  END IF
870  IF (debug_stress .AND. use_virial) THEN
871  stdeb = fconv*(virial%pv_virial - stdeb)
872  CALL para_env%sum(stdeb)
873  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
874  'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
875  END IF
876 
877  IF (ASSOCIATED(v_tau_rspace)) THEN
878  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
879  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
880  DO ispin = 1, nspins
881  CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
882  ! integrate over Tau-potential <nabla.a|V|nabla.b>
883  CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
884  hmat=scrm(ispin), &
885  pmat=matrix_p(ispin, 1), &
886  qs_env=qs_env, &
887  calculate_forces=.true., &
888  compute_tau=.true., &
889  basis_type="HARRIS", &
890  task_list_external=ec_env%task_list)
891  END DO
892 
893  IF (debug_forces) THEN
894  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
895  CALL para_env%sum(fodeb)
896  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb
897  END IF
898  IF (debug_stress .AND. use_virial) THEN
899  stdeb = fconv*(virial%pv_virial - stdeb)
900  CALL para_env%sum(stdeb)
901  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
902  'STRESS| INT Pout*dVhxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb)
903  END IF
904  END IF
905 
906  ! Stress-tensor
907  IF (use_virial) THEN
908  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
909  END IF
910 
911  ! delete scrm matrix
912  CALL dbcsr_deallocate_matrix_set(scrm)
913 
914  !----------------------------------------------------
915  ! Right-hand-side matrix B for linear response equations AX = B
916  !----------------------------------------------------
917 
918  ! RHS = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC * E_X[n] - alpha_gs * E_X[n]
919  ! = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC / alpha_GS * E_X[n]_GS - E_X[n]_GS
920  !
921  ! with v_Hxc[n] = v_H[n] + v_xc[n]
922  !
923  ! Actually v_H[n_in] same for DC and GS, just there for convenience
924  ! v_xc[n_in]_GS = 0 if GS is HF BUT =/0 if hybrid
925  ! so, we keep this general form
926 
927  NULLIFY (ec_env%matrix_hz)
928  CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
929  DO ispin = 1, nspins
930  ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
931  CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
932  CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
933  CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
934  END DO
935 
936  DO ispin = 1, nspins
937  ! v_rspace = v_rspace - v_rspace_in
938  ! = v_Hxc[n_in]_DC - v_Hxc[n_in]_GS
939  CALL pw_axpy(v_rspace_in(ispin), v_rspace(ispin), -1.0_dp)
940  END DO
941 
942  DO ispin = 1, nspins
943  CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
944  hmat=ec_env%matrix_hz(ispin), &
945  pmat=matrix_p(ispin, 1), &
946  qs_env=qs_env, &
947  calculate_forces=.false., &
948  basis_type="HARRIS", &
949  task_list_external=ec_env%task_list)
950  END DO
951 
952  ! Check if mGGA functionals are used
953  IF (dft_control%use_kinetic_energy_density) THEN
954 
955  ! If DC-DFT without mGGA functional, this needs to be allocated now.
956  IF (.NOT. ASSOCIATED(v_tau_rspace)) THEN
957  ALLOCATE (v_tau_rspace(nspins))
958  DO ispin = 1, nspins
959  CALL auxbas_pw_pool%create_pw(v_tau_rspace(ispin))
960  CALL pw_zero(v_tau_rspace(ispin))
961  END DO
962  END IF
963 
964  DO ispin = 1, nspins
965  ! v_tau_rspace = v_Hxc_tau[n_in]_DC - v_Hxc_tau[n_in]_GS
966  IF (ASSOCIATED(ec_env%vtau_rspace)) THEN
967  CALL pw_axpy(ec_env%vtau_rspace(ispin), v_tau_rspace(ispin), -1.0_dp)
968  END IF
969  ! integrate over Tau-potential <nabla.a|V|nabla.b>
970  CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
971  hmat=ec_env%matrix_hz(ispin), &
972  pmat=matrix_p(ispin, 1), &
973  qs_env=qs_env, &
974  calculate_forces=.false., compute_tau=.true., &
975  basis_type="HARRIS", &
976  task_list_external=ec_env%task_list)
977  END DO
978  END IF
979 
980  ! Need to also subtract HFX contribution of reference calculation from ec_env%matrix_hz
981  ! and/or add HFX contribution if DC-DFT ueses hybrid XC-functional
982  CALL add_exx_to_rhs(rhs=ec_env%matrix_hz, &
983  qs_env=qs_env, &
984  ext_hfx_section=ec_hfx_sections, &
985  x_data=ec_env%x_data, &
986  recalc_integrals=.false., &
987  do_admm=ec_env%do_ec_admm, &
988  do_ec=.true., &
989  do_exx=.false., &
990  reuse_hfx=ec_env%reuse_hfx)
991 
992  ! Core overlap
993  IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
994  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
995  CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
996  IF (debug_forces) THEN
997  fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
998  CALL para_env%sum(fodeb)
999  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb
1000  END IF
1001  IF (debug_stress .AND. use_virial) THEN
1002  stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
1003  CALL para_env%sum(stdeb)
1004  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1005  'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1006  END IF
1007 
1008  IF (debug_forces) THEN
1009  CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
1010  ALLOCATE (ftot(3, natom))
1011  CALL total_qs_force(ftot, force, atomic_kind_set)
1012  fodeb(1:3) = ftot(1:3, 1)
1013  DEALLOCATE (ftot)
1014  CALL para_env%sum(fodeb)
1015  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Force Explicit", fodeb
1016  END IF
1017 
1018  ! return pw grids
1019  DO ispin = 1, nspins
1020  CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1021  CALL auxbas_pw_pool%give_back_pw(v_rspace_in(ispin))
1022  IF (ASSOCIATED(v_tau_rspace)) THEN
1023  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1024  END IF
1025  END DO
1026 
1027  DEALLOCATE (v_rspace, v_rspace_in)
1028  IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
1029  !
1030  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1031  CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1032  CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1033 
1034  ! Stress tensor - volume terms need to be stored,
1035  ! for a sign correction in QS at the end of qs_force
1036  IF (use_virial) THEN
1037  IF (qs_env%energy_correction) THEN
1038  ec_env%ehartree = ehartree
1039  ec_env%exc = exc
1040  END IF
1041  END IF
1042 
1043  IF (debug_stress .AND. use_virial) THEN
1044  ! In total: -1.0*E_H
1045  stdeb = -1.0_dp*fconv*ehartree
1046  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1047  'STRESS| VOL 1st v_H[n_in]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb)
1048 
1049  stdeb = -1.0_dp*fconv*exc
1050  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1051  'STRESS| VOL 1st E_XC_DC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb)
1052 
1053  ! For debugging, create a second virial environment,
1054  ! apply volume terms immediately
1055  block
1056  TYPE(virial_type) :: virdeb
1057  virdeb = virial
1058 
1059  CALL para_env%sum(virdeb%pv_overlap)
1060  CALL para_env%sum(virdeb%pv_ekinetic)
1061  CALL para_env%sum(virdeb%pv_ppl)
1062  CALL para_env%sum(virdeb%pv_ppnl)
1063  CALL para_env%sum(virdeb%pv_ecore_overlap)
1064  CALL para_env%sum(virdeb%pv_ehartree)
1065  CALL para_env%sum(virdeb%pv_exc)
1066  CALL para_env%sum(virdeb%pv_exx)
1067  CALL para_env%sum(virdeb%pv_vdw)
1068  CALL para_env%sum(virdeb%pv_mp2)
1069  CALL para_env%sum(virdeb%pv_nlcc)
1070  CALL para_env%sum(virdeb%pv_gapw)
1071  CALL para_env%sum(virdeb%pv_lrigpw)
1072  CALL para_env%sum(virdeb%pv_virial)
1073  CALL symmetrize_virial(virdeb)
1074 
1075  ! apply stress-tensor 1st terms
1076  DO i = 1, 3
1077  virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*ehartree
1078  virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc &
1079  - 2.0_dp*ehartree
1080  virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc
1081  ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
1082  ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
1083  ! There should be a more elegant solution to that ...
1084  END DO
1085 
1086  CALL para_env%sum(sttot)
1087  stdeb = fconv*(virdeb%pv_virial - sttot)
1088  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1089  'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1090 
1091  stdeb = fconv*(virdeb%pv_virial)
1092  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1093  'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1094 
1095  CALL write_stress_tensor_components(virdeb, iounit, cell)
1096  CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .false.)
1097 
1098  END block
1099  END IF
1100 
1101  CALL timestop(handle)
1102 
1103  END SUBROUTINE ec_dc_build_ks_matrix_force
1104 
1105 ! **************************************************************************************************
1106 !> \brief ...
1107 !> \param qs_env ...
1108 !> \param ec_env ...
1109 !> \param calculate_forces ...
1110 ! **************************************************************************************************
1111  SUBROUTINE ec_disp(qs_env, ec_env, calculate_forces)
1112  TYPE(qs_environment_type), POINTER :: qs_env
1113  TYPE(energy_correction_type), POINTER :: ec_env
1114  LOGICAL, INTENT(IN) :: calculate_forces
1115 
1116  REAL(kind=dp) :: edisp
1117 
1118  CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, edisp, calculate_forces)
1119  IF (.NOT. calculate_forces) ec_env%edispersion = ec_env%edispersion + edisp
1120 
1121  END SUBROUTINE ec_disp
1122 
1123 ! **************************************************************************************************
1124 !> \brief Construction of the Core Hamiltonian Matrix
1125 !> Short version of qs_core_hamiltonian
1126 !> \param qs_env ...
1127 !> \param ec_env ...
1128 !> \author Creation (03.2014,JGH)
1129 ! **************************************************************************************************
1130  SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env)
1131  TYPE(qs_environment_type), POINTER :: qs_env
1132  TYPE(energy_correction_type), POINTER :: ec_env
1133 
1134  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_core_hamiltonian'
1135 
1136  INTEGER :: handle, nder, nimages
1137  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1138  LOGICAL :: calculate_forces, use_virial
1139  REAL(kind=dp) :: eps_ppnl
1140  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1141  TYPE(dft_control_type), POINTER :: dft_control
1142  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1143  POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1144  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1145  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1146  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1147  TYPE(qs_ks_env_type), POINTER :: ks_env
1148  TYPE(virial_type), POINTER :: virial
1149 
1150  CALL timeset(routinen, handle)
1151 
1152  NULLIFY (atomic_kind_set, cell_to_index, dft_control, ks_env, particle_set, qs_kind_set, virial)
1153 
1154  CALL get_qs_env(qs_env=qs_env, &
1155  atomic_kind_set=atomic_kind_set, &
1156  dft_control=dft_control, &
1157  particle_set=particle_set, &
1158  qs_kind_set=qs_kind_set, &
1159  ks_env=ks_env)
1160 
1161  ! no k-points possible
1162  nimages = dft_control%nimages
1163  IF (nimages /= 1) THEN
1164  cpabort("K-points for Harris functional not implemented")
1165  END IF
1166 
1167  ! check for GAPW/GAPW_XC
1168  IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1169  cpabort("Harris functional for GAPW not implemented")
1170  END IF
1171 
1172  ! Do not calculate forces or stress tensor here
1173  use_virial = .false.
1174  calculate_forces = .false.
1175 
1176  ! get neighbor lists, we need the full sab_orb list from the ec_env
1177  NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
1178  sab_orb => ec_env%sab_orb
1179  sac_ppl => ec_env%sac_ppl
1180  sap_ppnl => ec_env%sap_ppnl
1181 
1182  nder = 0
1183  ! Overlap and kinetic energy matrices
1184  CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
1185  matrix_name="OVERLAP MATRIX", &
1186  basis_type_a="HARRIS", &
1187  basis_type_b="HARRIS", &
1188  sab_nl=sab_orb)
1189  CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
1190  matrix_name="KINETIC ENERGY MATRIX", &
1191  basis_type="HARRIS", &
1192  sab_nl=sab_orb)
1193 
1194  ! initialize H matrix
1195  CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
1196  ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
1197  CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
1198  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
1199 
1200  ! add kinetic energy
1201  CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
1202  keep_sparsity=.true., name="CORE HAMILTONIAN MATRIX")
1203 
1204  ! Possible AE contributions (also ECP)
1205  IF (ASSOCIATED(sac_ae)) THEN
1206  cpabort("ECP/AE not available for energy corrections")
1207  ! missig code: sac_ae has to bee calculated and stored in ec_env
1208  ! build_core_ae: needs functionality to set the basis type at input
1209  CALL build_core_ae(ec_env%matrix_h, ec_env%matrix_p, force, &
1210  virial, calculate_forces, use_virial, nder, &
1211  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
1212  nimages, cell_to_index)
1213  END IF
1214  ! compute the ppl contribution to the core hamiltonian
1215  IF (ASSOCIATED(sac_ppl)) THEN
1216  CALL build_core_ppl(ec_env%matrix_h, ec_env%matrix_p, force, &
1217  virial, calculate_forces, use_virial, nder, &
1218  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
1219  nimages, cell_to_index, "HARRIS")
1220  END IF
1221 
1222  ! compute the ppnl contribution to the core hamiltonian ***
1223  eps_ppnl = dft_control%qs_control%eps_ppnl
1224  IF (ASSOCIATED(sap_ppnl)) THEN
1225  CALL build_core_ppnl(ec_env%matrix_h, ec_env%matrix_p, force, &
1226  virial, calculate_forces, use_virial, nder, &
1227  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
1228  eps_ppnl, nimages, cell_to_index, "HARRIS")
1229  END IF
1230 
1231  ! External field (nonperiodic case)
1232  ec_env%efield_nuclear = 0.0_dp
1233  CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1234 
1235  CALL timestop(handle)
1236 
1237  END SUBROUTINE ec_build_core_hamiltonian
1238 
1239 ! **************************************************************************************************
1240 !> \brief Solve KS equation for a given matrix
1241 !> calculate the complete KS matrix
1242 !> \param qs_env ...
1243 !> \param ec_env ...
1244 !> \par History
1245 !> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH]
1246 !> \author JGH
1247 ! **************************************************************************************************
1248  SUBROUTINE ec_build_ks_matrix(qs_env, ec_env)
1249  TYPE(qs_environment_type), POINTER :: qs_env
1250  TYPE(energy_correction_type), POINTER :: ec_env
1251 
1252  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_ks_matrix'
1253 
1254  CHARACTER(LEN=default_string_length) :: headline
1255  INTEGER :: handle, iounit, ispin, nspins
1256  LOGICAL :: calculate_forces, &
1257  do_adiabatic_rescaling, do_ec_hfx, &
1258  hfx_treat_lsd_in_core, use_virial
1259  REAL(dp) :: dummy_real, dummy_real2(2), eexc, evhxc, &
1260  t3
1261  TYPE(cp_logger_type), POINTER :: logger
1262  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_mat
1263  TYPE(dft_control_type), POINTER :: dft_control
1264  TYPE(pw_env_type), POINTER :: pw_env
1265  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1266  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace
1267  TYPE(qs_energy_type), POINTER :: energy
1268  TYPE(qs_ks_env_type), POINTER :: ks_env
1269  TYPE(qs_rho_type), POINTER :: rho
1270  TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
1271  ec_hfx_sections, ec_section
1272 
1273  CALL timeset(routinen, handle)
1274 
1275  logger => cp_get_default_logger()
1276  IF (logger%para_env%is_source()) THEN
1277  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1278  ELSE
1279  iounit = -1
1280  END IF
1281 
1282  ! get all information on the electronic density
1283  NULLIFY (auxbas_pw_pool, dft_control, energy, ks_env, rho, rho_r, tau_r)
1284  CALL get_qs_env(qs_env=qs_env, &
1285  dft_control=dft_control, &
1286  ks_env=ks_env, &
1287  rho=rho)
1288  nspins = dft_control%nspins
1289  calculate_forces = .false.
1290  use_virial = .false.
1291 
1292  ! Kohn-Sham matrix
1293  IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
1294  CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
1295  DO ispin = 1, nspins
1296  headline = "KOHN-SHAM MATRIX"
1297  ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
1298  CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=trim(headline), &
1299  template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1300  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb)
1301  CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp)
1302  END DO
1303 
1304  NULLIFY (pw_env)
1305  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1306  cpassert(ASSOCIATED(pw_env))
1307 
1308  ! Exact exchange contribution (hybrid functionals)
1309  ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
1310  ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF")
1311  CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
1312 
1313  IF (do_ec_hfx) THEN
1314 
1315  ! Check what works
1316  adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING")
1317  CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1318  IF (do_adiabatic_rescaling) THEN
1319  CALL cp_abort(__location__, "Adiabatic rescaling NYI for energy correction")
1320  END IF
1321  CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
1322  IF (hfx_treat_lsd_in_core) THEN
1323  CALL cp_abort(__location__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction")
1324  END IF
1325 
1326  ! calculate the density matrix for the fitted mo_coeffs
1327  IF (dft_control%do_admm) THEN
1328 
1329  IF (dft_control%do_admm_mo) THEN
1330  IF (qs_env%run_rtp) THEN
1331  CALL rtp_admm_calc_rho_aux(qs_env)
1332  ELSE
1333  CALL admm_mo_calc_rho_aux(qs_env)
1334  END IF
1335  ELSEIF (dft_control%do_admm_dm) THEN
1336  CALL admm_dm_calc_rho_aux(qs_env)
1337  END IF
1338  END IF
1339 
1340  ! Get exact exchange energy
1341  dummy_real = 0.0_dp
1342  t3 = 0.0_dp
1343  CALL get_qs_env(qs_env, energy=energy)
1344  CALL calculate_exx(qs_env=qs_env, &
1345  unit_nr=iounit, &
1346  hfx_sections=ec_hfx_sections, &
1347  x_data=ec_env%x_data, &
1348  do_gw=.false., &
1349  do_admm=ec_env%do_ec_admm, &
1350  calc_forces=.false., &
1351  reuse_hfx=ec_env%reuse_hfx, &
1352  do_im_time=.false., &
1353  e_ex_from_gw=dummy_real, &
1354  e_admm_from_gw=dummy_real2, &
1355  t3=dummy_real)
1356 
1357  ! Save exchange energy
1358  ec_env%ex = energy%ex
1359  ! Save EXX ADMM XC correction
1360  IF (ec_env%do_ec_admm) THEN
1361  ec_env%exc_aux_fit = energy%exc_aux_fit + energy%exc
1362  END IF
1363 
1364  ! Add exact echange contribution of EC to EC Hamiltonian
1365  ! do_ec = .FALSE prevents subtraction of HFX contribution of reference calculation
1366  ! do_exx = .FALSE. prevents subtraction of reference XC contribution
1367  ks_mat => ec_env%matrix_ks(:, 1)
1368  CALL add_exx_to_rhs(rhs=ks_mat, &
1369  qs_env=qs_env, &
1370  ext_hfx_section=ec_hfx_sections, &
1371  x_data=ec_env%x_data, &
1372  recalc_integrals=.false., &
1373  do_admm=ec_env%do_ec_admm, &
1374  do_ec=.false., &
1375  do_exx=.false., &
1376  reuse_hfx=ec_env%reuse_hfx)
1377 
1378  END IF
1379 
1380  ! v_rspace and v_tau_rspace are generated from the auxbas pool
1381  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1382  NULLIFY (v_rspace, v_tau_rspace)
1383  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1384  vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
1385 
1386  IF (.NOT. ASSOCIATED(v_rspace)) THEN
1387  ALLOCATE (v_rspace(nspins))
1388  DO ispin = 1, nspins
1389  CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
1390  CALL pw_zero(v_rspace(ispin))
1391  END DO
1392  END IF
1393 
1394  evhxc = 0.0_dp
1395  CALL qs_rho_get(rho, rho_r=rho_r)
1396  IF (ASSOCIATED(v_tau_rspace)) THEN
1397  CALL qs_rho_get(rho, tau_r=tau_r)
1398  END IF
1399  DO ispin = 1, nspins
1400  ! Add v_hartree + v_xc = v_rspace
1401  CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1402  CALL pw_axpy(ec_env%vh_rspace, v_rspace(ispin))
1403  ! integrate over potential <a|V|b>
1404  CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1405  hmat=ec_env%matrix_ks(ispin, 1), &
1406  qs_env=qs_env, &
1407  calculate_forces=.false., &
1408  basis_type="HARRIS", &
1409  task_list_external=ec_env%task_list)
1410 
1411  IF (ASSOCIATED(v_tau_rspace)) THEN
1412  ! integrate over Tau-potential <nabla.a|V|nabla.b>
1413  CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1414  CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1415  hmat=ec_env%matrix_ks(ispin, 1), &
1416  qs_env=qs_env, &
1417  calculate_forces=.false., &
1418  compute_tau=.true., &
1419  basis_type="HARRIS", &
1420  task_list_external=ec_env%task_list)
1421  END IF
1422 
1423  ! calclulate Int(vhxc*rho)dr and Int(vtau*tau)dr
1424  evhxc = evhxc + pw_integral_ab(rho_r(ispin), v_rspace(ispin))/ &
1425  v_rspace(1)%pw_grid%dvol
1426  IF (ASSOCIATED(v_tau_rspace)) THEN
1427  evhxc = evhxc + pw_integral_ab(tau_r(ispin), v_tau_rspace(ispin))/ &
1428  v_tau_rspace(ispin)%pw_grid%dvol
1429  END IF
1430 
1431  END DO
1432 
1433  ! return pw grids
1434  DO ispin = 1, nspins
1435  CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1436  IF (ASSOCIATED(v_tau_rspace)) THEN
1437  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1438  END IF
1439  END DO
1440  DEALLOCATE (v_rspace)
1441  IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
1442 
1443  ! energies
1444  ec_env%exc = eexc
1445  ec_env%vhxc = evhxc
1446 
1447  ! add the core matrix
1448  DO ispin = 1, nspins
1449  CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, &
1450  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1451  CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, &
1452  dft_control%qs_control%eps_filter_matrix)
1453  END DO
1454 
1455  CALL timestop(handle)
1456 
1457  END SUBROUTINE ec_build_ks_matrix
1458 
1459 ! **************************************************************************************************
1460 !> \brief Construction of the Core Hamiltonian Matrix
1461 !> Short version of qs_core_hamiltonian
1462 !> \param qs_env ...
1463 !> \param ec_env ...
1464 !> \param matrix_p ...
1465 !> \param matrix_s ...
1466 !> \param matrix_w ...
1467 !> \author Creation (03.2014,JGH)
1468 ! **************************************************************************************************
1469  SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env, matrix_p, matrix_s, matrix_w)
1470  TYPE(qs_environment_type), POINTER :: qs_env
1471  TYPE(energy_correction_type), POINTER :: ec_env
1472  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s, matrix_w
1473 
1474  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_core_hamiltonian_force'
1475 
1476  INTEGER :: handle, iounit, nder, nimages
1477  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1478  LOGICAL :: calculate_forces, debug_forces, &
1479  debug_stress, use_virial
1480  REAL(kind=dp) :: eps_ppnl, fconv
1481  REAL(kind=dp), DIMENSION(3) :: fodeb
1482  REAL(kind=dp), DIMENSION(3, 3) :: stdeb, sttot
1483  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1484  TYPE(cell_type), POINTER :: cell
1485  TYPE(cp_logger_type), POINTER :: logger
1486  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
1487  TYPE(dft_control_type), POINTER :: dft_control
1488  TYPE(mp_para_env_type), POINTER :: para_env
1489  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1490  POINTER :: sab_orb, sac_ppl, sap_ppnl
1491  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1492  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1493  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1494  TYPE(qs_ks_env_type), POINTER :: ks_env
1495  TYPE(virial_type), POINTER :: virial
1496 
1497  CALL timeset(routinen, handle)
1498 
1499  debug_forces = ec_env%debug_forces
1500  debug_stress = ec_env%debug_stress
1501 
1502  logger => cp_get_default_logger()
1503  IF (logger%para_env%is_source()) THEN
1504  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1505  ELSE
1506  iounit = -1
1507  END IF
1508 
1509  calculate_forces = .true.
1510 
1511  ! no k-points possible
1512  NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
1513  CALL get_qs_env(qs_env=qs_env, &
1514  cell=cell, &
1515  dft_control=dft_control, &
1516  force=force, &
1517  ks_env=ks_env, &
1518  para_env=para_env, &
1519  virial=virial)
1520  nimages = dft_control%nimages
1521  IF (nimages /= 1) THEN
1522  cpabort("K-points for Harris functional not implemented")
1523  END IF
1524 
1525  ! check for virial
1526  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1527 
1528  fconv = 1.0e-9_dp*pascal/cell%deth
1529  IF (debug_stress .AND. use_virial) THEN
1530  sttot = virial%pv_virial
1531  END IF
1532 
1533  ! check for GAPW/GAPW_XC
1534  IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1535  cpabort("Harris functional for GAPW not implemented")
1536  END IF
1537 
1538  ! get neighbor lists, we need the full sab_orb list from the ec_env
1539  NULLIFY (sab_orb, sac_ppl, sap_ppnl)
1540  sab_orb => ec_env%sab_orb
1541  sac_ppl => ec_env%sac_ppl
1542  sap_ppnl => ec_env%sap_ppnl
1543 
1544  ! initialize src matrix
1545  NULLIFY (scrm)
1546  CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
1547  ALLOCATE (scrm(1, 1)%matrix)
1548  CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_s(1, 1)%matrix)
1549  CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
1550 
1551  nder = 1
1552  IF (SIZE(matrix_p, 1) == 2) THEN
1553  CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
1554  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1555  CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
1556  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1557  END IF
1558 
1559  ! Overlap and kinetic energy matrices
1560  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
1561  IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
1562  CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
1563  matrix_name="OVERLAP MATRIX", &
1564  basis_type_a="HARRIS", &
1565  basis_type_b="HARRIS", &
1566  sab_nl=sab_orb, calculate_forces=.true., &
1567  matrixkp_p=matrix_w)
1568 
1569  IF (debug_forces) THEN
1570  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
1571  CALL para_env%sum(fodeb)
1572  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
1573  fodeb(1:3) = force(1)%kinetic(1:3, 1)
1574  END IF
1575  IF (debug_stress .AND. use_virial) THEN
1576  stdeb = fconv*(virial%pv_overlap - stdeb)
1577  CALL para_env%sum(stdeb)
1578  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1579  'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
1580  stdeb = virial%pv_ekinetic
1581  END IF
1582  CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
1583  matrix_name="KINETIC ENERGY MATRIX", &
1584  basis_type="HARRIS", &
1585  sab_nl=sab_orb, calculate_forces=.true., &
1586  matrixkp_p=matrix_p)
1587  IF (debug_forces) THEN
1588  fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
1589  CALL para_env%sum(fodeb)
1590  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb
1591  END IF
1592  IF (debug_stress .AND. use_virial) THEN
1593  stdeb = fconv*(virial%pv_ekinetic - stdeb)
1594  CALL para_env%sum(stdeb)
1595  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1596  'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
1597  END IF
1598  IF (SIZE(matrix_p, 1) == 2) THEN
1599  CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
1600  alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1601  END IF
1602 
1603  ! compute the ppl contribution to the core hamiltonian
1604  NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
1605  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1606  atomic_kind_set=atomic_kind_set)
1607 
1608  IF (ASSOCIATED(sac_ppl)) THEN
1609  IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
1610  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
1611  CALL build_core_ppl(scrm, matrix_p, force, &
1612  virial, calculate_forces, use_virial, nder, &
1613  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
1614  nimages, cell_to_index, "HARRIS")
1615  IF (calculate_forces .AND. debug_forces) THEN
1616  fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
1617  CALL para_env%sum(fodeb)
1618  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb
1619  END IF
1620  IF (debug_stress .AND. use_virial) THEN
1621  stdeb = fconv*(virial%pv_ppl - stdeb)
1622  CALL para_env%sum(stdeb)
1623  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1624  'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb)
1625  END IF
1626  END IF
1627 
1628  ! compute the ppnl contribution to the core hamiltonian ***
1629  eps_ppnl = dft_control%qs_control%eps_ppnl
1630  IF (ASSOCIATED(sap_ppnl)) THEN
1631  IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
1632  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
1633  CALL build_core_ppnl(scrm, matrix_p, force, &
1634  virial, calculate_forces, use_virial, nder, &
1635  qs_kind_set, atomic_kind_set, particle_set, &
1636  sab_orb, sap_ppnl, eps_ppnl, &
1637  nimages, cell_to_index, "HARRIS")
1638  IF (calculate_forces .AND. debug_forces) THEN
1639  fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
1640  CALL para_env%sum(fodeb)
1641  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb
1642  END IF
1643  IF (debug_stress .AND. use_virial) THEN
1644  stdeb = fconv*(virial%pv_ppnl - stdeb)
1645  CALL para_env%sum(stdeb)
1646  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1647  'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb)
1648  END IF
1649  END IF
1650 
1651  ! External field (nonperiodic case)
1652  ec_env%efield_nuclear = 0.0_dp
1653  IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%efield(1:3, 1)
1654  CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1655  IF (calculate_forces .AND. debug_forces) THEN
1656  fodeb(1:3) = force(1)%efield(1:3, 1) - fodeb(1:3)
1657  CALL para_env%sum(fodeb)
1658  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dEfield", fodeb
1659  END IF
1660  IF (debug_stress .AND. use_virial) THEN
1661  stdeb = fconv*(virial%pv_virial - sttot)
1662  CALL para_env%sum(stdeb)
1663  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1664  'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1665  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") ' '
1666  END IF
1667 
1668  ! delete scr matrix
1669  CALL dbcsr_deallocate_matrix_set(scrm)
1670 
1671  CALL timestop(handle)
1672 
1673  END SUBROUTINE ec_build_core_hamiltonian_force
1674 
1675 ! **************************************************************************************************
1676 !> \brief Solve KS equation for a given matrix
1677 !> \brief calculate the complete KS matrix
1678 !> \param qs_env ...
1679 !> \param ec_env ...
1680 !> \par History
1681 !> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH]
1682 !> \author JGH
1683 ! **************************************************************************************************
1684  SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
1685  TYPE(qs_environment_type), POINTER :: qs_env
1686  TYPE(energy_correction_type), POINTER :: ec_env
1687 
1688  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_ks_matrix_force'
1689 
1690  INTEGER :: handle, i, iounit, ispin, natom, nspins
1691  LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
1692  use_virial
1693  REAL(dp) :: dehartree, dummy_real, dummy_real2(2), &
1694  eexc, ehartree, eovrl, exc, fconv
1695  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
1696  REAL(dp), DIMENSION(3) :: fodeb
1697  REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
1698  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1699  TYPE(cell_type), POINTER :: cell
1700  TYPE(cp_logger_type), POINTER :: logger
1701  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, scrm
1702  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1703  TYPE(dft_control_type), POINTER :: dft_control
1704  TYPE(mp_para_env_type), POINTER :: para_env
1705  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1706  POINTER :: sab_orb
1707  TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhodn_tot_gspace, &
1708  v_hartree_gspace
1709  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rhoout_g
1710  TYPE(pw_c1d_gs_type), POINTER :: rho_core
1711  TYPE(pw_env_type), POINTER :: pw_env
1712  TYPE(pw_poisson_type), POINTER :: poisson_env
1713  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1714  TYPE(pw_r3d_rs_type) :: dv_hartree_rspace, v_hartree_rspace, &
1715  vtot_rspace
1716  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rhoout_r, tau_r, tauout_r, &
1717  v_rspace, v_tau_rspace, v_xc, v_xc_tau
1718  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1719  TYPE(qs_ks_env_type), POINTER :: ks_env
1720  TYPE(qs_rho_type), POINTER :: rho
1721  TYPE(section_vals_type), POINTER :: ec_hfx_sections, xc_section
1722  TYPE(virial_type), POINTER :: virial
1723 
1724  CALL timeset(routinen, handle)
1725 
1726  debug_forces = ec_env%debug_forces
1727  debug_stress = ec_env%debug_stress
1728 
1729  logger => cp_get_default_logger()
1730  IF (logger%para_env%is_source()) THEN
1731  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1732  ELSE
1733  iounit = -1
1734  END IF
1735 
1736  ! get all information on the electronic density
1737  NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, &
1738  matrix_ks, matrix_p, matrix_s, para_env, rho, rho_core, &
1739  rho_g, rho_r, sab_orb, tau_r, virial)
1740  CALL get_qs_env(qs_env=qs_env, &
1741  cell=cell, &
1742  dft_control=dft_control, &
1743  force=force, &
1744  ks_env=ks_env, &
1745  matrix_ks=matrix_ks, &
1746  para_env=para_env, &
1747  rho=rho, &
1748  sab_orb=sab_orb, &
1749  virial=virial)
1750 
1751  nspins = dft_control%nspins
1752  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1753 
1754  fconv = 1.0e-9_dp*pascal/cell%deth
1755  IF (debug_stress .AND. use_virial) THEN
1756  sttot = virial%pv_virial
1757  END IF
1758 
1759  NULLIFY (pw_env)
1760  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1761  cpassert(ASSOCIATED(pw_env))
1762 
1763  NULLIFY (auxbas_pw_pool, poisson_env)
1764  ! gets the tmp grids
1765  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1766  poisson_env=poisson_env)
1767 
1768  ! Calculate the Hartree potential
1769  CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
1770  CALL auxbas_pw_pool%create_pw(rhodn_tot_gspace)
1771  CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
1772 
1773  CALL pw_transfer(ec_env%vh_rspace, v_hartree_rspace)
1774 
1775  ! calculate output density on grid
1776  ! rho_in(R): CALL qs_rho_get(rho, rho_r=rho_r)
1777  ! rho_in(G): CALL qs_rho_get(rho, rho_g=rho_g)
1778  CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
1779  NULLIFY (rhoout_r, rhoout_g)
1780  ALLOCATE (rhoout_r(nspins), rhoout_g(nspins))
1781  DO ispin = 1, nspins
1782  CALL auxbas_pw_pool%create_pw(rhoout_r(ispin))
1783  CALL auxbas_pw_pool%create_pw(rhoout_g(ispin))
1784  END DO
1785  CALL auxbas_pw_pool%create_pw(dv_hartree_rspace)
1786  CALL auxbas_pw_pool%create_pw(vtot_rspace)
1787 
1788  CALL pw_zero(rhodn_tot_gspace)
1789  DO ispin = 1, nspins
1790  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1791  rho=rhoout_r(ispin), &
1792  rho_gspace=rhoout_g(ispin), &
1793  basis_type="HARRIS", &
1794  task_list_external=ec_env%task_list)
1795  END DO
1796 
1797  ! Save Harris on real space grid for use in properties
1798  ALLOCATE (ec_env%rhoout_r(nspins))
1799  DO ispin = 1, nspins
1800  CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
1801  CALL pw_copy(rhoout_r(ispin), ec_env%rhoout_r(ispin))
1802  END DO
1803 
1804  NULLIFY (tauout_r)
1805  IF (dft_control%use_kinetic_energy_density) THEN
1806  block
1807  TYPE(pw_c1d_gs_type) :: tauout_g
1808  ALLOCATE (tauout_r(nspins))
1809  DO ispin = 1, nspins
1810  CALL auxbas_pw_pool%create_pw(tauout_r(ispin))
1811  END DO
1812  CALL auxbas_pw_pool%create_pw(tauout_g)
1813 
1814  DO ispin = 1, nspins
1815  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1816  rho=tauout_r(ispin), &
1817  rho_gspace=tauout_g, &
1818  compute_tau=.true., &
1819  basis_type="HARRIS", &
1820  task_list_external=ec_env%task_list)
1821  END DO
1822 
1823  CALL auxbas_pw_pool%give_back_pw(tauout_g)
1824  END block
1825  END IF
1826 
1827  IF (use_virial) THEN
1828 
1829  ! Calculate the Hartree potential
1830  CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1831 
1832  ! Get the total input density in g-space [ions + electrons]
1833  CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1834 
1835  ! make rho_tot_gspace with output density
1836  CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1837  CALL pw_copy(rho_core, rhodn_tot_gspace)
1838  DO ispin = 1, dft_control%nspins
1839  CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1840  END DO
1841 
1842  ! Volume and Green function terms
1843  h_stress(:, :) = 0.0_dp
1844  CALL pw_poisson_solve(poisson_env, &
1845  density=rho_tot_gspace, & ! n_in
1846  ehartree=ehartree, &
1847  vhartree=v_hartree_gspace, & ! v_H[n_in]
1848  h_stress=h_stress, &
1849  aux_density=rhodn_tot_gspace) ! n_out
1850 
1851  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1852  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1853 
1854  IF (debug_stress) THEN
1855  stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1856  CALL para_env%sum(stdeb)
1857  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1858  'STRESS| GREEN 1st v_H[n_in]*n_out ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1859  END IF
1860 
1861  ! activate stress calculation
1862  virial%pv_calculate = .true.
1863 
1864  NULLIFY (v_rspace, v_tau_rspace)
1865  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1866  vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
1867 
1868  ! Stress tensor XC-functional GGA contribution
1869  virial%pv_exc = virial%pv_exc - virial%pv_xc
1870  virial%pv_virial = virial%pv_virial - virial%pv_xc
1871 
1872  IF (debug_stress) THEN
1873  stdeb = -1.0_dp*fconv*virial%pv_xc
1874  CALL para_env%sum(stdeb)
1875  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1876  'STRESS| GGA 1st E_xc[Pin] ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1877  END IF
1878 
1879  IF (ASSOCIATED(v_rspace)) THEN
1880  DO ispin = 1, nspins
1881  CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1882  END DO
1883  DEALLOCATE (v_rspace)
1884  END IF
1885  IF (ASSOCIATED(v_tau_rspace)) THEN
1886  DO ispin = 1, nspins
1887  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1888  END DO
1889  DEALLOCATE (v_tau_rspace)
1890  END IF
1891  CALL pw_zero(rhodn_tot_gspace)
1892 
1893  END IF
1894 
1895  ! rho_out - rho_in
1896  DO ispin = 1, nspins
1897  CALL pw_axpy(rho_r(ispin), rhoout_r(ispin), -1.0_dp)
1898  CALL pw_axpy(rho_g(ispin), rhoout_g(ispin), -1.0_dp)
1899  CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1900  IF (dft_control%use_kinetic_energy_density) CALL pw_axpy(tau_r(ispin), tauout_r(ispin), -1.0_dp)
1901  END DO
1902 
1903  ! calculate associated hartree potential
1904  IF (use_virial) THEN
1905 
1906  ! Stress tensor - 2nd derivative Volume and Green function contribution
1907  h_stress(:, :) = 0.0_dp
1908  CALL pw_poisson_solve(poisson_env, &
1909  density=rhodn_tot_gspace, & ! delta_n
1910  ehartree=dehartree, &
1911  vhartree=v_hartree_gspace, & ! v_H[delta_n]
1912  h_stress=h_stress, &
1913  aux_density=rho_tot_gspace) ! n_in
1914 
1915  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1916 
1917  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1918  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1919 
1920  IF (debug_stress) THEN
1921  stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1922  CALL para_env%sum(stdeb)
1923  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1924  'STRESS| GREEN 2nd V_H[dP]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1925  END IF
1926 
1927  ELSE
1928  ! v_H[dn]
1929  CALL pw_poisson_solve(poisson_env, rhodn_tot_gspace, dehartree, &
1930  v_hartree_gspace)
1931  END IF
1932 
1933  CALL pw_transfer(v_hartree_gspace, dv_hartree_rspace)
1934  CALL pw_scale(dv_hartree_rspace, dv_hartree_rspace%pw_grid%dvol)
1935  ! Getting nuclear force contribution from the core charge density
1936  ! Vh(rho_in + rho_c) + Vh(rho_out - rho_in)
1937  CALL pw_transfer(v_hartree_rspace, vtot_rspace)
1938  CALL pw_axpy(dv_hartree_rspace, vtot_rspace)
1939  IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
1940  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
1941  CALL integrate_v_core_rspace(vtot_rspace, qs_env)
1942  IF (debug_forces) THEN
1943  fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
1944  CALL para_env%sum(fodeb)
1945  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb
1946  END IF
1947  IF (debug_stress .AND. use_virial) THEN
1948  stdeb = fconv*(virial%pv_ehartree - stdeb)
1949  CALL para_env%sum(stdeb)
1950  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1951  'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb)
1952  END IF
1953  !
1954  ! Pulay force from Tr P_in (V_H(drho)+ Fxc(rho_in)*drho)
1955  ! RHS of CPKS equations: (V_H(drho)+ Fxc(rho_in)*drho)*C0
1956  ! Fxc*drho term
1957  xc_section => ec_env%xc_section
1958 
1959  IF (use_virial) virial%pv_xc = 0.0_dp
1960  NULLIFY (v_xc, v_xc_tau)
1961  CALL create_kernel(qs_env, &
1962  vxc=v_xc, &
1963  vxc_tau=v_xc_tau, &
1964  rho=rho, &
1965  rho1_r=rhoout_r, &
1966  rho1_g=rhoout_g, &
1967  tau1_r=tauout_r, &
1968  xc_section=xc_section, &
1969  compute_virial=use_virial, &
1970  virial_xc=virial%pv_xc)
1971 
1972  IF (use_virial) THEN
1973  ! Stress-tensor XC-functional 2nd GGA terms
1974  virial%pv_exc = virial%pv_exc + virial%pv_xc
1975  virial%pv_virial = virial%pv_virial + virial%pv_xc
1976  END IF
1977  IF (debug_stress .AND. use_virial) THEN
1978  stdeb = 1.0_dp*fconv*virial%pv_xc
1979  CALL para_env%sum(stdeb)
1980  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1981  'STRESS| GGA 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1982  END IF
1983  !
1984  CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s)
1985  NULLIFY (ec_env%matrix_hz)
1986  CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
1987  DO ispin = 1, nspins
1988  ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
1989  CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1, 1)%matrix)
1990  CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1, 1)%matrix)
1991  CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
1992  END DO
1993  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1994  ! vtot = v_xc(ispin) + dv_hartree
1995  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1996  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1997 
1998  ! Stress-tensor 2nd derivative integral contribution
1999  IF (use_virial) THEN
2000  pv_loc = virial%pv_virial
2001  END IF
2002 
2003  DO ispin = 1, nspins
2004  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2005  CALL pw_axpy(dv_hartree_rspace, v_xc(ispin))
2006  CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
2007  hmat=ec_env%matrix_hz(ispin), &
2008  pmat=matrix_p(ispin, 1), &
2009  qs_env=qs_env, &
2010  calculate_forces=.true.)
2011  END DO
2012 
2013  IF (debug_forces) THEN
2014  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2015  CALL para_env%sum(fodeb)
2016  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKdrho", fodeb
2017  END IF
2018  IF (debug_stress .AND. use_virial) THEN
2019  stdeb = fconv*(virial%pv_virial - stdeb)
2020  CALL para_env%sum(stdeb)
2021  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2022  'STRESS| INT 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2023  END IF
2024 
2025  IF (ASSOCIATED(v_xc_tau)) THEN
2026  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2027  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2028 
2029  DO ispin = 1, nspins
2030  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2031  CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
2032  hmat=ec_env%matrix_hz(ispin), &
2033  pmat=matrix_p(ispin, 1), &
2034  qs_env=qs_env, &
2035  compute_tau=.true., &
2036  calculate_forces=.true.)
2037  END DO
2038 
2039  IF (debug_forces) THEN
2040  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2041  CALL para_env%sum(fodeb)
2042  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKtaudtau", fodeb
2043  END IF
2044  IF (debug_stress .AND. use_virial) THEN
2045  stdeb = fconv*(virial%pv_virial - stdeb)
2046  CALL para_env%sum(stdeb)
2047  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2048  'STRESS| INT 2nd f_xctau[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2049  END IF
2050  END IF
2051  ! Stress-tensor 2nd derivative integral contribution
2052  IF (use_virial) THEN
2053  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2054  END IF
2055 
2056  ! initialize srcm matrix
2057  NULLIFY (scrm)
2058  CALL dbcsr_allocate_matrix_set(scrm, nspins)
2059  DO ispin = 1, nspins
2060  ALLOCATE (scrm(ispin)%matrix)
2061  CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
2062  CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
2063  CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
2064  END DO
2065 
2066  ! v_rspace and v_tau_rspace are generated from the auxbas pool
2067  NULLIFY (v_rspace, v_tau_rspace)
2068 
2069  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
2070  vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
2071 
2072  IF (use_virial) THEN
2073  eexc = 0.0_dp
2074  IF (ASSOCIATED(v_rspace)) THEN
2075  DO ispin = 1, nspins
2076  ! 2nd deriv xc-volume term
2077  eexc = eexc + pw_integral_ab(rhoout_r(ispin), v_rspace(ispin))
2078  END DO
2079  END IF
2080  IF (ASSOCIATED(v_tau_rspace)) THEN
2081  DO ispin = 1, nspins
2082  ! 2nd deriv xc-volume term
2083  eexc = eexc + pw_integral_ab(tauout_r(ispin), v_tau_rspace(ispin))
2084  END DO
2085  END IF
2086  END IF
2087 
2088  IF (.NOT. ASSOCIATED(v_rspace)) THEN
2089  ALLOCATE (v_rspace(nspins))
2090  DO ispin = 1, nspins
2091  CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
2092  CALL pw_zero(v_rspace(ispin))
2093  END DO
2094  END IF
2095 
2096  ! Stress-tensor contribution derivative of integrand
2097  ! int v_Hxc[n^în]*n^out
2098  IF (use_virial) THEN
2099  pv_loc = virial%pv_virial
2100  END IF
2101 
2102  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2103  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2104  DO ispin = 1, nspins
2105  ! Add v_hartree + v_xc = v_rspace
2106  CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
2107  CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
2108  ! integrate over potential <a|V|b>
2109  CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
2110  hmat=scrm(ispin), &
2111  pmat=ec_env%matrix_p(ispin, 1), &
2112  qs_env=qs_env, &
2113  calculate_forces=.true., &
2114  basis_type="HARRIS", &
2115  task_list_external=ec_env%task_list)
2116  END DO
2117 
2118  IF (debug_forces) THEN
2119  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2120  CALL para_env%sum(fodeb)
2121  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb
2122  END IF
2123  IF (debug_stress .AND. use_virial) THEN
2124  stdeb = fconv*(virial%pv_virial - stdeb)
2125  CALL para_env%sum(stdeb)
2126  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2127  'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2128  END IF
2129 
2130  ! Stress-tensor
2131  IF (use_virial) THEN
2132  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2133  END IF
2134 
2135  IF (ASSOCIATED(v_tau_rspace)) THEN
2136  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2137  DO ispin = 1, nspins
2138  ! integrate over Tau-potential <nabla.a|V|nabla.b>
2139  CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
2140  CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
2141  hmat=scrm(ispin), &
2142  pmat=ec_env%matrix_p(ispin, 1), &
2143  qs_env=qs_env, &
2144  calculate_forces=.true., &
2145  compute_tau=.true., &
2146  basis_type="HARRIS", &
2147  task_list_external=ec_env%task_list)
2148  END DO
2149  IF (debug_forces) THEN
2150  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2151  CALL para_env%sum(fodeb)
2152  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb
2153  END IF
2154  END IF
2155 
2156 !------------------------------------------------------------------------------
2157 ! HFX direct force
2158 !------------------------------------------------------------------------------
2159 
2160  ! If hybrid functional
2161  ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF")
2162  CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
2163 
2164  IF (do_ec_hfx) THEN
2165 
2166  IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2167  IF (use_virial) virial%pv_fock_4c = 0.0_dp
2168 
2169  CALL calculate_exx(qs_env=qs_env, &
2170  unit_nr=iounit, &
2171  hfx_sections=ec_hfx_sections, &
2172  x_data=ec_env%x_data, &
2173  do_gw=.false., &
2174  do_admm=ec_env%do_ec_admm, &
2175  calc_forces=.true., &
2176  reuse_hfx=ec_env%reuse_hfx, &
2177  do_im_time=.false., &
2178  e_ex_from_gw=dummy_real, &
2179  e_admm_from_gw=dummy_real2, &
2180  t3=dummy_real)
2181 
2182  IF (use_virial) THEN
2183  virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2184  virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2185  virial%pv_calculate = .false.
2186  END IF
2187  IF (debug_forces) THEN
2188  fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2189  CALL para_env%sum(fodeb)
2190  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*hfx ", fodeb
2191  END IF
2192  IF (debug_stress .AND. use_virial) THEN
2193  stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2194  CALL para_env%sum(stdeb)
2195  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2196  'STRESS| Pout*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2197  END IF
2198 
2199  END IF
2200 
2201 !------------------------------------------------------------------------------
2202 
2203  ! delete scrm matrix
2204  CALL dbcsr_deallocate_matrix_set(scrm)
2205 
2206  ! return pw grids
2207  CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
2208  DO ispin = 1, nspins
2209  CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
2210  IF (ASSOCIATED(v_tau_rspace)) THEN
2211  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
2212  END IF
2213  END DO
2214  IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
2215 
2216  ! Core overlap
2217  IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
2218  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
2219  CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
2220  IF (debug_forces) THEN
2221  fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
2222  CALL para_env%sum(fodeb)
2223  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb
2224  END IF
2225  IF (debug_stress .AND. use_virial) THEN
2226  stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
2227  CALL para_env%sum(stdeb)
2228  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2229  'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2230  END IF
2231 
2232  IF (debug_forces) THEN
2233  CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2234  ALLOCATE (ftot(3, natom))
2235  CALL total_qs_force(ftot, force, atomic_kind_set)
2236  fodeb(1:3) = ftot(1:3, 1)
2237  DEALLOCATE (ftot)
2238  CALL para_env%sum(fodeb)
2239  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Force Explicit", fodeb
2240  END IF
2241 
2242  DEALLOCATE (v_rspace)
2243  !
2244  CALL auxbas_pw_pool%give_back_pw(dv_hartree_rspace)
2245  CALL auxbas_pw_pool%give_back_pw(vtot_rspace)
2246  DO ispin = 1, nspins
2247  CALL auxbas_pw_pool%give_back_pw(rhoout_r(ispin))
2248  CALL auxbas_pw_pool%give_back_pw(rhoout_g(ispin))
2249  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2250  END DO
2251  DEALLOCATE (rhoout_r, rhoout_g, v_xc)
2252  IF (ASSOCIATED(tauout_r)) THEN
2253  DO ispin = 1, nspins
2254  CALL auxbas_pw_pool%give_back_pw(tauout_r(ispin))
2255  END DO
2256  DEALLOCATE (tauout_r)
2257  END IF
2258  IF (ASSOCIATED(v_xc_tau)) THEN
2259  DO ispin = 1, nspins
2260  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2261  END DO
2262  DEALLOCATE (v_xc_tau)
2263  END IF
2264  CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
2265  CALL auxbas_pw_pool%give_back_pw(rhodn_tot_gspace)
2266 
2267  ! Stress tensor - volume terms need to be stored,
2268  ! for a sign correction in QS at the end of qs_force
2269  IF (use_virial) THEN
2270  IF (qs_env%energy_correction) THEN
2271  ec_env%ehartree = ehartree + dehartree
2272  ec_env%exc = exc + eexc
2273  END IF
2274  END IF
2275 
2276  IF (debug_stress .AND. use_virial) THEN
2277  ! In total: -1.0*E_H
2278  stdeb = -1.0_dp*fconv*ehartree
2279  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2280  'STRESS| VOL 1st v_H[n_in]*n_out', one_third_sum_diag(stdeb), det_3x3(stdeb)
2281 
2282  stdeb = -1.0_dp*fconv*exc
2283  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2284  'STRESS| VOL 1st E_XC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb)
2285 
2286  stdeb = -1.0_dp*fconv*dehartree
2287  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2288  'STRESS| VOL 2nd v_H[dP]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb)
2289 
2290  stdeb = -1.0_dp*fconv*eexc
2291  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2292  'STRESS| VOL 2nd v_XC[n_in]*dP', one_third_sum_diag(stdeb), det_3x3(stdeb)
2293 
2294  ! For debugging, create a second virial environment,
2295  ! apply volume terms immediately
2296  block
2297  TYPE(virial_type) :: virdeb
2298  virdeb = virial
2299 
2300  CALL para_env%sum(virdeb%pv_overlap)
2301  CALL para_env%sum(virdeb%pv_ekinetic)
2302  CALL para_env%sum(virdeb%pv_ppl)
2303  CALL para_env%sum(virdeb%pv_ppnl)
2304  CALL para_env%sum(virdeb%pv_ecore_overlap)
2305  CALL para_env%sum(virdeb%pv_ehartree)
2306  CALL para_env%sum(virdeb%pv_exc)
2307  CALL para_env%sum(virdeb%pv_exx)
2308  CALL para_env%sum(virdeb%pv_vdw)
2309  CALL para_env%sum(virdeb%pv_mp2)
2310  CALL para_env%sum(virdeb%pv_nlcc)
2311  CALL para_env%sum(virdeb%pv_gapw)
2312  CALL para_env%sum(virdeb%pv_lrigpw)
2313  CALL para_env%sum(virdeb%pv_virial)
2314  CALL symmetrize_virial(virdeb)
2315 
2316  ! apply stress-tensor 1st and 2nd volume terms
2317  DO i = 1, 3
2318  virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*(ehartree + dehartree)
2319  virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc - eexc &
2320  - 2.0_dp*(ehartree + dehartree)
2321  virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc - eexc
2322  ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
2323  ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
2324  ! There should be a more elegant solution to that ...
2325  END DO
2326 
2327  CALL para_env%sum(sttot)
2328  stdeb = fconv*(virdeb%pv_virial - sttot)
2329  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2330  'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2331 
2332  stdeb = fconv*(virdeb%pv_virial)
2333  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2334  'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2335 
2336  CALL write_stress_tensor_components(virdeb, iounit, cell)
2337  CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .false.)
2338 
2339  END block
2340  END IF
2341 
2342  CALL timestop(handle)
2343 
2344  END SUBROUTINE ec_build_ks_matrix_force
2345 
2346 ! **************************************************************************************************
2347 !> \brief Solve KS equation for a given matrix
2348 !> \param qs_env ...
2349 !> \param ec_env ...
2350 !> \par History
2351 !> 03.2014 created [JGH]
2352 !> \author JGH
2353 ! **************************************************************************************************
2354  SUBROUTINE ec_ks_solver(qs_env, ec_env)
2355 
2356  TYPE(qs_environment_type), POINTER :: qs_env
2357  TYPE(energy_correction_type), POINTER :: ec_env
2358 
2359  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_ks_solver'
2360 
2361  CHARACTER(LEN=default_string_length) :: headline
2362  INTEGER :: handle, ispin, nspins
2363  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, pmat, smat, wmat
2364  TYPE(dft_control_type), POINTER :: dft_control
2365 
2366  CALL timeset(routinen, handle)
2367 
2368  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2369  nspins = dft_control%nspins
2370 
2371  ! create density matrix
2372  IF (.NOT. ASSOCIATED(ec_env%matrix_p)) THEN
2373  headline = "DENSITY MATRIX"
2374  CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
2375  DO ispin = 1, nspins
2376  ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
2377  CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
2378  template=ec_env%matrix_s(1, 1)%matrix)
2379  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
2380  END DO
2381  END IF
2382  ! create energy weighted density matrix
2383  IF (.NOT. ASSOCIATED(ec_env%matrix_w)) THEN
2384  headline = "ENERGY WEIGHTED DENSITY MATRIX"
2385  CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
2386  DO ispin = 1, nspins
2387  ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
2388  CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
2389  template=ec_env%matrix_s(1, 1)%matrix)
2390  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
2391  END DO
2392  END IF
2393 
2394  IF (ec_env%mao) THEN
2395  CALL mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2396  ELSE
2397  ksmat => ec_env%matrix_ks
2398  smat => ec_env%matrix_s
2399  pmat => ec_env%matrix_p
2400  wmat => ec_env%matrix_w
2401  END IF
2402 
2403  SELECT CASE (ec_env%ks_solver)
2404  CASE (ec_diagonalization)
2405  CALL ec_diag_solver(qs_env, ksmat, smat, pmat, wmat)
2406  CASE (ec_ot_diag)
2407  CALL ec_ot_diag_solver(qs_env, ec_env, ksmat, smat, pmat, wmat)
2408  CASE (ec_matrix_sign, ec_matrix_trs4, ec_matrix_tc2)
2409  CALL ec_ls_init(qs_env, ksmat, smat)
2410  CALL ec_ls_solver(qs_env, pmat, wmat, ec_ls_method=ec_env%ks_solver)
2411  CASE DEFAULT
2412  cpassert(.false.)
2413  END SELECT
2414 
2415  ! OUtput density available now
2416 
2417 ! HFX contribution to Harris functional and energy
2418 ! Can't calculate this earlier, cause ec_env%matrix_p doesnt exist yet
2419 !------------------------------------------------------------------------------
2420 
2421 ! ! Exact exchange contribution (hybrid functionals)
2422 ! ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
2423 ! ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF")
2424 ! CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
2425 !
2426 ! IF (do_ec_hfx) THEN
2427 !
2428 ! ! Check what works
2429 ! IF (dft_control%do_admm) THEN
2430 ! CALL cp_warn(__LOCATION__, "Energy correction with hybrid functional does not use ADMM.")
2431 ! END IF
2432 !
2433 ! adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING")
2434 ! CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
2435 ! IF (do_adiabatic_rescaling) THEN
2436 ! CALL cp_abort(__LOCATION__, "Adiabatic rescaling NYI for energy correction")
2437 ! END IF
2438 ! CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
2439 ! IF (hfx_treat_lsd_in_core) THEN
2440 ! CALL cp_abort(__LOCATION__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction")
2441 ! END IF
2442 !
2443 ! ! Exchange matrix
2444 ! IF (ASSOCIATED(ec_env%matrix_x)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_x)
2445 ! CALL dbcsr_allocate_matrix_set(ec_env%matrix_x, nspins)
2446 ! DO ispin = 1, nspins
2447 ! headline = "EXCHANGE MATRIX"
2448 ! ALLOCATE (ec_env%matrix_x(ispin)%matrix)
2449 ! CALL dbcsr_create(ec_env%matrix_x(ispin)%matrix, name=TRIM(headline), &
2450 ! template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
2451 ! CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_x(ispin)%matrix, ec_env%sab_orb)
2452 ! CALL dbcsr_set(ec_env%matrix_x(ispin)%matrix, 0.0_dp)
2453 ! END DO
2454 !
2455 ! ! Get exact exchange energy (fraction) and its contribution to the EC hamiltonian
2456 ! should_update=.TRUE.
2457 ! ks_mat => ec_env%matrix_ks(:,1)
2458 ! CALL ec_hfx_contributions(qs_env, ks_mat, matrix_p, &
2459 ! ec_hfx_sections, ec_env%x_data, use_virial, &
2460 ! should_update, calculate_forces, matrix_x = ec_env%matrix_x, ex = ec_env%ex)
2461 !
2462 ! END IF
2463 
2464 !------------------------------------------------------------------------------
2465 
2466  IF (ec_env%mao) THEN
2467  CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2468  END IF
2469 
2470  CALL timestop(handle)
2471 
2472  END SUBROUTINE ec_ks_solver
2473 
2474 ! **************************************************************************************************
2475 !> \brief Create matrices with MAO sizes
2476 !> \param ec_env ...
2477 !> \param ksmat ...
2478 !> \param smat ...
2479 !> \param pmat ...
2480 !> \param wmat ...
2481 !> \par History
2482 !> 08.2016 created [JGH]
2483 !> \author JGH
2484 ! **************************************************************************************************
2485  SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2486 
2487  TYPE(energy_correction_type), POINTER :: ec_env
2488  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat
2489 
2490  CHARACTER(LEN=*), PARAMETER :: routinen = 'mao_create_matrices'
2491 
2492  INTEGER :: handle, ispin, nspins
2493  INTEGER, DIMENSION(:), POINTER :: col_blk_sizes
2494  TYPE(dbcsr_distribution_type) :: dbcsr_dist
2495  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
2496  TYPE(dbcsr_type) :: cgmat
2497 
2498  CALL timeset(routinen, handle)
2499 
2500  mao_coef => ec_env%mao_coef
2501 
2502  NULLIFY (ksmat, smat, pmat, wmat)
2503  nspins = SIZE(ec_env%matrix_ks, 1)
2504  CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
2505  CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1)
2506  CALL dbcsr_allocate_matrix_set(smat, nspins, 1)
2507  DO ispin = 1, nspins
2508  ALLOCATE (ksmat(ispin, 1)%matrix)
2509  CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO KS mat", &
2510  matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2511  col_blk_size=col_blk_sizes, nze=0)
2512  ALLOCATE (smat(ispin, 1)%matrix)
2513  CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO S mat", &
2514  matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2515  col_blk_size=col_blk_sizes, nze=0)
2516  END DO
2517  !
2518  CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
2519  DO ispin = 1, nspins
2520  CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, &
2521  0.0_dp, cgmat)
2522  CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix)
2523  CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, &
2524  0.0_dp, cgmat)
2525  CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix)
2526  END DO
2527  CALL dbcsr_release(cgmat)
2528 
2529  CALL dbcsr_allocate_matrix_set(pmat, nspins, 1)
2530  DO ispin = 1, nspins
2531  ALLOCATE (pmat(ispin, 1)%matrix)
2532  CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name="MAO P mat")
2533  CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb)
2534  END DO
2535 
2536  CALL dbcsr_allocate_matrix_set(wmat, nspins, 1)
2537  DO ispin = 1, nspins
2538  ALLOCATE (wmat(ispin, 1)%matrix)
2539  CALL dbcsr_create(wmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name="MAO W mat")
2540  CALL cp_dbcsr_alloc_block_from_nbl(wmat(ispin, 1)%matrix, ec_env%sab_orb)
2541  END DO
2542 
2543  CALL timestop(handle)
2544 
2545  END SUBROUTINE mao_create_matrices
2546 
2547 ! **************************************************************************************************
2548 !> \brief Release matrices with MAO sizes
2549 !> \param ec_env ...
2550 !> \param ksmat ...
2551 !> \param smat ...
2552 !> \param pmat ...
2553 !> \param wmat ...
2554 !> \par History
2555 !> 08.2016 created [JGH]
2556 !> \author JGH
2557 ! **************************************************************************************************
2558  SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2559 
2560  TYPE(energy_correction_type), POINTER :: ec_env
2561  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat
2562 
2563  CHARACTER(LEN=*), PARAMETER :: routinen = 'mao_release_matrices'
2564 
2565  INTEGER :: handle, ispin, nspins
2566  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
2567  TYPE(dbcsr_type) :: cgmat
2568 
2569  CALL timeset(routinen, handle)
2570 
2571  mao_coef => ec_env%mao_coef
2572  nspins = SIZE(mao_coef, 1)
2573 
2574  ! save pmat/wmat in full basis format
2575  CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
2576  DO ispin = 1, nspins
2577  CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2578  CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2579  ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.true.)
2580  CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, wmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2581  CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2582  ec_env%matrix_w(ispin, 1)%matrix, retain_sparsity=.true.)
2583  END DO
2584  CALL dbcsr_release(cgmat)
2585 
2586  CALL dbcsr_deallocate_matrix_set(ksmat)
2587  CALL dbcsr_deallocate_matrix_set(smat)
2588  CALL dbcsr_deallocate_matrix_set(pmat)
2589  CALL dbcsr_deallocate_matrix_set(wmat)
2590 
2591  CALL timestop(handle)
2592 
2593  END SUBROUTINE mao_release_matrices
2594 
2595 ! **************************************************************************************************
2596 !> \brief Solve KS equation using diagonalization
2597 !> \param qs_env ...
2598 !> \param matrix_ks ...
2599 !> \param matrix_s ...
2600 !> \param matrix_p ...
2601 !> \param matrix_w ...
2602 !> \par History
2603 !> 03.2014 created [JGH]
2604 !> \author JGH
2605 ! **************************************************************************************************
2606  SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p, matrix_w)
2607 
2608  TYPE(qs_environment_type), POINTER :: qs_env
2609  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
2610 
2611  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_diag_solver'
2612 
2613  INTEGER :: handle, ispin, nmo(2), nsize, nspins
2614  REAL(kind=dp) :: eps_filter, focc(2)
2615  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
2616  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2617  TYPE(cp_fm_struct_type), POINTER :: fm_struct
2618  TYPE(cp_fm_type) :: fm_ks, fm_mo, fm_ortho
2619  TYPE(dbcsr_type), POINTER :: buf1_dbcsr, buf2_dbcsr, buf3_dbcsr, &
2620  ortho_dbcsr, ref_matrix
2621  TYPE(dft_control_type), POINTER :: dft_control
2622  TYPE(mp_para_env_type), POINTER :: para_env
2623 
2624  CALL timeset(routinen, handle)
2625 
2626  NULLIFY (blacs_env, para_env)
2627  CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
2628 
2629  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2630  eps_filter = dft_control%qs_control%eps_filter_matrix
2631  nspins = dft_control%nspins
2632 
2633  nmo = 0
2634  CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
2635  focc = 1._dp
2636  IF (nspins == 1) THEN
2637  focc = 2._dp
2638  nmo(1) = nmo(1)/2
2639  END IF
2640 
2641  CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
2642  ALLOCATE (eigenvalues(nsize))
2643 
2644  NULLIFY (fm_struct, ref_matrix)
2645  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
2646  ncol_global=nsize, para_env=para_env)
2647  CALL cp_fm_create(fm_ortho, fm_struct)
2648  CALL cp_fm_create(fm_ks, fm_struct)
2649  CALL cp_fm_create(fm_mo, fm_struct)
2650  CALL cp_fm_struct_release(fm_struct)
2651 
2652  ! factorization
2653  ref_matrix => matrix_s(1, 1)%matrix
2654  NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2655  CALL dbcsr_init_p(ortho_dbcsr)
2656  CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, &
2657  matrix_type=dbcsr_type_no_symmetry)
2658  CALL dbcsr_init_p(buf1_dbcsr)
2659  CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, &
2660  matrix_type=dbcsr_type_no_symmetry)
2661  CALL dbcsr_init_p(buf2_dbcsr)
2662  CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, &
2663  matrix_type=dbcsr_type_no_symmetry)
2664  CALL dbcsr_init_p(buf3_dbcsr)
2665  CALL dbcsr_create(buf3_dbcsr, template=ref_matrix, &
2666  matrix_type=dbcsr_type_no_symmetry)
2667 
2668  ref_matrix => matrix_s(1, 1)%matrix
2669  CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho)
2670  CALL cp_fm_cholesky_decompose(fm_ortho)
2671  CALL cp_fm_triangular_invert(fm_ortho)
2672  CALL cp_fm_set_all(fm_ks, 0.0_dp)
2673  CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks, "U")
2674  CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr)
2675  DO ispin = 1, nspins
2676  ! calculate ZHZ(T)
2677  CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr)
2678  CALL dbcsr_multiply("N", "N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, &
2679  0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2680  CALL dbcsr_multiply("T", "N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, &
2681  0.0_dp, buf1_dbcsr, filter_eps=eps_filter)
2682  ! copy to fm format
2683  CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks)
2684  CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues)
2685  ! back transform of mos c = Z(T)*c
2686  CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2687  CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2688  0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2689  ! density matrix
2690  CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp)
2691  CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf2_dbcsr, &
2692  1.0_dp, matrix_p(ispin, 1)%matrix, &
2693  retain_sparsity=.true., last_k=nmo(ispin))
2694 
2695  ! energy weighted density matrix
2696  CALL dbcsr_set(matrix_w(ispin, 1)%matrix, 0.0_dp)
2697  CALL cp_fm_column_scale(fm_mo, eigenvalues)
2698  CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2699  CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2700  0.0_dp, buf3_dbcsr, filter_eps=eps_filter)
2701  CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf3_dbcsr, &
2702  1.0_dp, matrix_w(ispin, 1)%matrix, &
2703  retain_sparsity=.true., last_k=nmo(ispin))
2704  END DO
2705 
2706  CALL cp_fm_release(fm_ks)
2707  CALL cp_fm_release(fm_mo)
2708  CALL cp_fm_release(fm_ortho)
2709  CALL dbcsr_release(ortho_dbcsr)
2710  CALL dbcsr_release(buf1_dbcsr)
2711  CALL dbcsr_release(buf2_dbcsr)
2712  CALL dbcsr_release(buf3_dbcsr)
2713  DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2714  DEALLOCATE (eigenvalues)
2715 
2716  CALL timestop(handle)
2717 
2718  END SUBROUTINE ec_diag_solver
2719 
2720 ! **************************************************************************************************
2721 !> \brief Calculate the energy correction
2722 !> \param ec_env ...
2723 !> \param unit_nr ...
2724 !> \author Creation (03.2014,JGH)
2725 ! **************************************************************************************************
2726  SUBROUTINE ec_energy(ec_env, unit_nr)
2727  TYPE(energy_correction_type) :: ec_env
2728  INTEGER, INTENT(IN) :: unit_nr
2729 
2730  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_energy'
2731 
2732  INTEGER :: handle, ispin, nspins
2733  REAL(kind=dp) :: eband, trace
2734 
2735  CALL timeset(routinen, handle)
2736 
2737  nspins = SIZE(ec_env%matrix_ks, 1)
2738  DO ispin = 1, nspins
2739  CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace)
2740  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T65,F16.10)') 'Tr[PS] ', trace
2741  END DO
2742 
2743  ! Total energy depends on energy correction method
2744  SELECT CASE (ec_env%energy_functional)
2745  CASE (ec_functional_harris)
2746 
2747  ! Get energy of "band structure" term
2748  eband = 0.0_dp
2749  DO ispin = 1, nspins
2750  CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace)
2751  eband = eband + trace
2752  END DO
2753  ec_env%eband = eband + ec_env%efield_nuclear
2754 
2755  ! Add Harris functional "correction" terms
2756  ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion - ec_env%ex
2757  IF (unit_nr > 0) THEN
2758  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Eband ", ec_env%eband
2759  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree
2760  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc
2761  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex
2762  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Evhxc ", ec_env%vhxc
2763  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion
2764  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Harris Functional ", ec_env%etotal
2765  END IF
2766 
2767  CASE (ec_functional_dc)
2768 
2769  ! Core hamiltonian energy
2770  CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore, SIZE(ec_env%matrix_p, 1))
2771 
2772  ec_env%ecore = ec_env%ecore + ec_env%efield_nuclear
2773  ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion &
2774  + ec_env%ex + ec_env%exc_aux_fit
2775 
2776  IF (unit_nr > 0) THEN
2777  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ecore ", ec_env%ecore
2778  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree
2779  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc
2780  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex
2781  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc_aux_fit", ec_env%exc_aux_fit
2782  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion
2783  WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal
2784  END IF
2785 
2786  CASE DEFAULT
2787 
2788  cpassert(.false.)
2789 
2790  END SELECT
2791 
2792  CALL timestop(handle)
2793 
2794  END SUBROUTINE ec_energy
2795 
2796 ! **************************************************************************************************
2797 !> \brief builds either the full neighborlist or neighborlists of molecular
2798 !> \brief subsets, depending on parameter values
2799 !> \param qs_env ...
2800 !> \param ec_env ...
2801 !> \par History
2802 !> 2012.07 created [Martin Haeufel]
2803 !> 2016.07 Adapted for Harris functional [JGH]
2804 !> \author Martin Haeufel
2805 ! **************************************************************************************************
2806  SUBROUTINE ec_build_neighborlist(qs_env, ec_env)
2807  TYPE(qs_environment_type), POINTER :: qs_env
2808  TYPE(energy_correction_type), POINTER :: ec_env
2809 
2810  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_neighborlist'
2811 
2812  INTEGER :: handle, ikind, nkind, zat
2813  LOGICAL :: gth_potential_present, &
2814  sgp_potential_present, &
2815  skip_load_balance_distributed
2816  LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present, orb_present, &
2817  ppl_present, ppnl_present
2818  REAL(dp) :: subcells
2819  REAL(dp), ALLOCATABLE, DIMENSION(:) :: c_radius, orb_radius, ppl_radius, &
2820  ppnl_radius
2821  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
2822  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2823  TYPE(cell_type), POINTER :: cell
2824  TYPE(dft_control_type), POINTER :: dft_control
2825  TYPE(distribution_1d_type), POINTER :: distribution_1d
2826  TYPE(distribution_2d_type), POINTER :: distribution_2d
2827  TYPE(gth_potential_type), POINTER :: gth_potential
2828  TYPE(gto_basis_set_type), POINTER :: basis_set
2829  TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
2830  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
2831  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2832  POINTER :: sab_cn, sab_vdw
2833  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2834  TYPE(qs_dispersion_type), POINTER :: dispersion_env
2835  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2836  TYPE(qs_kind_type), POINTER :: qs_kind
2837  TYPE(qs_ks_env_type), POINTER :: ks_env
2838  TYPE(sgp_potential_type), POINTER :: sgp_potential
2839 
2840  CALL timeset(routinen, handle)
2841 
2842  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
2843  CALL get_qs_kind_set(qs_kind_set, gth_potential_present=gth_potential_present, &
2844  sgp_potential_present=sgp_potential_present)
2845  nkind = SIZE(qs_kind_set)
2846  ALLOCATE (c_radius(nkind), default_present(nkind))
2847  ALLOCATE (orb_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind))
2848  ALLOCATE (orb_present(nkind), ppl_present(nkind), ppnl_present(nkind))
2849  ALLOCATE (pair_radius(nkind, nkind))
2850  ALLOCATE (atom2d(nkind))
2851 
2852  CALL get_qs_env(qs_env, &
2853  atomic_kind_set=atomic_kind_set, &
2854  cell=cell, &
2855  distribution_2d=distribution_2d, &
2856  local_particles=distribution_1d, &
2857  particle_set=particle_set, &
2858  molecule_set=molecule_set)
2859 
2860  CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
2861  molecule_set, .false., particle_set)
2862 
2863  DO ikind = 1, nkind
2864  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
2865  qs_kind => qs_kind_set(ikind)
2866  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="HARRIS")
2867  IF (ASSOCIATED(basis_set)) THEN
2868  orb_present(ikind) = .true.
2869  CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind))
2870  ELSE
2871  orb_present(ikind) = .false.
2872  orb_radius(ikind) = 0.0_dp
2873  END IF
2874  CALL get_qs_kind(qs_kind, gth_potential=gth_potential, sgp_potential=sgp_potential)
2875  IF (gth_potential_present .OR. sgp_potential_present) THEN
2876  IF (ASSOCIATED(gth_potential)) THEN
2877  CALL get_potential(potential=gth_potential, &
2878  ppl_present=ppl_present(ikind), &
2879  ppl_radius=ppl_radius(ikind), &
2880  ppnl_present=ppnl_present(ikind), &
2881  ppnl_radius=ppnl_radius(ikind))
2882  ELSE IF (ASSOCIATED(sgp_potential)) THEN
2883  CALL get_potential(potential=sgp_potential, &
2884  ppl_present=ppl_present(ikind), &
2885  ppl_radius=ppl_radius(ikind), &
2886  ppnl_present=ppnl_present(ikind), &
2887  ppnl_radius=ppnl_radius(ikind))
2888  ELSE
2889  ppl_present(ikind) = .false.
2890  ppl_radius(ikind) = 0.0_dp
2891  ppnl_present(ikind) = .false.
2892  ppnl_radius(ikind) = 0.0_dp
2893  END IF
2894  END IF
2895  END DO
2896 
2897  CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
2898 
2899  ! overlap
2900  CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
2901  CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
2902  subcells=subcells, nlname="sab_orb")
2903  ! pseudopotential
2904  IF (gth_potential_present .OR. sgp_potential_present) THEN
2905  IF (any(ppl_present)) THEN
2906  CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
2907  CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
2908  subcells=subcells, operator_type="ABC", nlname="sac_ppl")
2909  END IF
2910 
2911  IF (any(ppnl_present)) THEN
2912  CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
2913  CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
2914  subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
2915  END IF
2916  END IF
2917 
2918  ! Build the neighbor lists for the vdW pair potential
2919  c_radius(:) = 0.0_dp
2920  dispersion_env => ec_env%dispersion_env
2921  sab_vdw => dispersion_env%sab_vdw
2922  sab_cn => dispersion_env%sab_cn
2923  IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
2924  c_radius(:) = dispersion_env%rc_disp
2925  default_present = .true. !include all atoms in vdW (even without basis)
2926  CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
2927  CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
2928  subcells=subcells, operator_type="PP", nlname="sab_vdw")
2929  dispersion_env%sab_vdw => sab_vdw
2930  IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
2931  dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
2932  ! Build the neighbor lists for coordination numbers as needed by the DFT-D3 method
2933  DO ikind = 1, nkind
2934  CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
2935  c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
2936  END DO
2937  CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
2938  CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
2939  subcells=subcells, operator_type="PP", nlname="sab_cn")
2940  dispersion_env%sab_cn => sab_cn
2941  END IF
2942  END IF
2943 
2944  ! Release work storage
2945  CALL atom2d_cleanup(atom2d)
2946  DEALLOCATE (atom2d)
2947  DEALLOCATE (orb_present, default_present, ppl_present, ppnl_present)
2948  DEALLOCATE (orb_radius, ppl_radius, ppnl_radius, c_radius)
2949  DEALLOCATE (pair_radius)
2950 
2951  ! Task list
2952  CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
2953  skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
2954  IF (ASSOCIATED(ec_env%task_list)) CALL deallocate_task_list(ec_env%task_list)
2955  CALL allocate_task_list(ec_env%task_list)
2956  CALL generate_qs_task_list(ks_env, ec_env%task_list, &
2957  reorder_rs_grid_ranks=.false., soft_valid=.false., &
2958  skip_load_balance_distributed=skip_load_balance_distributed, &
2959  basis_type="HARRIS", sab_orb_external=ec_env%sab_orb)
2960 
2961  CALL timestop(handle)
2962 
2963  END SUBROUTINE ec_build_neighborlist
2964 
2965 ! **************************************************************************************************
2966 !> \brief ...
2967 !> \param qs_env ...
2968 !> \param ec_env ...
2969 ! **************************************************************************************************
2970  SUBROUTINE ec_properties(qs_env, ec_env)
2971  TYPE(qs_environment_type), POINTER :: qs_env
2972  TYPE(energy_correction_type), POINTER :: ec_env
2973 
2974  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_properties'
2975 
2976  CHARACTER(LEN=8), DIMENSION(3) :: rlab
2977  CHARACTER(LEN=default_path_length) :: filename, my_pos_voro
2978  CHARACTER(LEN=default_string_length) :: description
2979  INTEGER :: akind, handle, i, ia, iatom, idir, ikind, iounit, ispin, maxmom, nspins, &
2980  reference, should_print_bqb, should_print_voro, unit_nr, unit_nr_voro
2981  LOGICAL :: append_voro, magnetic, periodic, &
2982  voro_print_txt
2983  REAL(kind=dp) :: charge, dd, focc, tmp
2984  REAL(kind=dp), DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip
2985  REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
2986  TYPE(atomic_kind_type), POINTER :: atomic_kind
2987  TYPE(cell_type), POINTER :: cell
2988  TYPE(cp_logger_type), POINTER :: logger
2989  TYPE(cp_result_type), POINTER :: results
2990  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
2991  TYPE(dft_control_type), POINTER :: dft_control
2992  TYPE(distribution_1d_type), POINTER :: local_particles
2993  TYPE(mp_para_env_type), POINTER :: para_env
2994  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2995  TYPE(pw_env_type), POINTER :: pw_env
2996  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2997  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2998  TYPE(pw_r3d_rs_type) :: rho_elec_rspace
2999  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3000  TYPE(section_vals_type), POINTER :: ec_section, print_key, print_key_bqb, &
3001  print_key_voro
3002 
3003  CALL timeset(routinen, handle)
3004 
3005  rlab(1) = "X"
3006  rlab(2) = "Y"
3007  rlab(3) = "Z"
3008 
3009  logger => cp_get_default_logger()
3010  IF (logger%para_env%is_source()) THEN
3011  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
3012  ELSE
3013  iounit = -1
3014  END IF
3015 
3016  NULLIFY (dft_control)
3017  CALL get_qs_env(qs_env, dft_control=dft_control)
3018  nspins = dft_control%nspins
3019 
3020  ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
3021  print_key => section_vals_get_subs_vals(section_vals=ec_section, &
3022  subsection_name="PRINT%MOMENTS")
3023 
3024  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3025 
3026  maxmom = section_get_ival(section_vals=ec_section, &
3027  keyword_name="PRINT%MOMENTS%MAX_MOMENT")
3028  periodic = section_get_lval(section_vals=ec_section, &
3029  keyword_name="PRINT%MOMENTS%PERIODIC")
3030  reference = section_get_ival(section_vals=ec_section, &
3031  keyword_name="PRINT%MOMENTS%REFERENCE")
3032  magnetic = section_get_lval(section_vals=ec_section, &
3033  keyword_name="PRINT%MOMENTS%MAGNETIC")
3034  NULLIFY (ref_point)
3035  CALL section_vals_val_get(ec_section, "PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
3036  unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=ec_section, &
3037  print_key_path="PRINT%MOMENTS", extension=".dat", &
3038  middle_name="moments", log_filename=.false.)
3039 
3040  IF (iounit > 0) THEN
3041  IF (unit_nr /= iounit .AND. unit_nr > 0) THEN
3042  INQUIRE (unit=unit_nr, name=filename)
3043  WRITE (unit=iounit, fmt="(/,T2,A,2(/,T3,A),/)") &
3044  "MOMENTS", "The electric/magnetic moments are written to file:", &
3045  trim(filename)
3046  ELSE
3047  WRITE (unit=iounit, fmt="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
3048  END IF
3049  END IF
3050 
3051  IF (periodic) THEN
3052  cpabort("Periodic moments not implemented with EC")
3053  ELSE
3054  cpassert(maxmom < 2)
3055  cpassert(.NOT. magnetic)
3056  IF (maxmom == 1) THEN
3057  CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
3058  ! reference point
3059  CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3060  ! nuclear contribution
3061  cdip = 0.0_dp
3062  CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
3063  qs_kind_set=qs_kind_set, local_particles=local_particles)
3064  DO ikind = 1, SIZE(local_particles%n_el)
3065  DO ia = 1, local_particles%n_el(ikind)
3066  iatom = local_particles%list(ikind)%array(ia)
3067  ! fold atomic positions back into unit cell
3068  ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
3069  ria = ria - rcc
3070  atomic_kind => particle_set(iatom)%atomic_kind
3071  CALL get_atomic_kind(atomic_kind, kind_number=akind)
3072  CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
3073  cdip(1:3) = cdip(1:3) - charge*ria(1:3)
3074  END DO
3075  END DO
3076  CALL para_env%sum(cdip)
3077  !
3078  ! direct density contribution
3079  CALL ec_efield_integrals(qs_env, ec_env, rcc)
3080  !
3081  pdip = 0.0_dp
3082  DO ispin = 1, nspins
3083  DO idir = 1, 3
3084  CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, &
3085  ec_env%efield%dipmat(idir)%matrix, tmp)
3086  pdip(idir) = pdip(idir) + tmp
3087  END DO
3088  END DO
3089  !
3090  ! response contribution
3091  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
3092  NULLIFY (moments)
3093  CALL dbcsr_allocate_matrix_set(moments, 4)
3094  DO i = 1, 4
3095  ALLOCATE (moments(i)%matrix)
3096  CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
3097  CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
3098  END DO
3099  CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
3100  !
3101  focc = 2.0_dp
3102  IF (nspins == 2) focc = 1.0_dp
3103  rdip = 0.0_dp
3104  DO ispin = 1, nspins
3105  DO idir = 1, 3
3106  CALL dbcsr_dot(ec_env%matrix_z(ispin)%matrix, moments(idir)%matrix, tmp)
3107  rdip(idir) = rdip(idir) + tmp
3108  END DO
3109  END DO
3110  CALL dbcsr_deallocate_matrix_set(moments)
3111  !
3112  tdip = -(rdip + pdip + cdip)
3113  IF (unit_nr > 0) THEN
3114  WRITE (unit_nr, "(T3,A)") "Dipoles are based on the traditional operator."
3115  dd = sqrt(sum(tdip(1:3)**2))*debye
3116  WRITE (unit_nr, "(T3,A)") "Dipole moment [Debye]"
3117  WRITE (unit_nr, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
3118  (trim(rlab(i)), "=", tdip(i)*debye, i=1, 3), "Total=", dd
3119  END IF
3120  END IF
3121  END IF
3122 
3123  CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
3124  basis_section=ec_section, print_key_path="PRINT%MOMENTS")
3125  CALL get_qs_env(qs_env=qs_env, results=results)
3126  description = "[DIPOLE]"
3127  CALL cp_results_erase(results=results, description=description)
3128  CALL put_results(results=results, description=description, values=tdip(1:3))
3129  END IF
3130 
3131  ! Do a Voronoi Integration or write a compressed BQB File
3132  print_key_voro => section_vals_get_subs_vals(ec_section, "PRINT%VORONOI")
3133  print_key_bqb => section_vals_get_subs_vals(ec_section, "PRINT%E_DENSITY_BQB")
3134  IF (btest(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file)) THEN
3135  should_print_voro = 1
3136  ELSE
3137  should_print_voro = 0
3138  END IF
3139  IF (btest(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file)) THEN
3140  should_print_bqb = 1
3141  ELSE
3142  should_print_bqb = 0
3143  END IF
3144  IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0)) THEN
3145 
3146  CALL get_qs_env(qs_env=qs_env, &
3147  pw_env=pw_env)
3148  CALL pw_env_get(pw_env=pw_env, &
3149  auxbas_pw_pool=auxbas_pw_pool, &
3150  pw_pools=pw_pools)
3151  CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
3152 
3153  IF (dft_control%nspins > 1) THEN
3154 
3155  ! add Pout and Pz
3156  CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3157  CALL pw_axpy(ec_env%rhoout_r(2), rho_elec_rspace)
3158 
3159  CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3160  CALL pw_axpy(ec_env%rhoz_r(2), rho_elec_rspace)
3161  ELSE
3162 
3163  ! add Pout and Pz
3164  CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3165  CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3166  END IF ! nspins
3167 
3168  IF (should_print_voro /= 0) THEN
3169  CALL section_vals_val_get(print_key_voro, "OUTPUT_TEXT", l_val=voro_print_txt)
3170  IF (voro_print_txt) THEN
3171  append_voro = section_get_lval(ec_section, "PRINT%VORONOI%APPEND")
3172  my_pos_voro = "REWIND"
3173  IF (append_voro) THEN
3174  my_pos_voro = "APPEND"
3175  END IF
3176  unit_nr_voro = cp_print_key_unit_nr(logger, ec_section, "PRINT%VORONOI", extension=".voronoi", &
3177  file_position=my_pos_voro, log_filename=.false.)
3178  ELSE
3179  unit_nr_voro = 0
3180  END IF
3181  ELSE
3182  unit_nr_voro = 0
3183  END IF
3184 
3185  CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
3186  unit_nr_voro, qs_env, rho_elec_rspace)
3187 
3188  CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
3189 
3190  IF (unit_nr_voro > 0) THEN
3191  CALL cp_print_key_finished_output(unit_nr_voro, logger, ec_section, "PRINT%VORONOI")
3192  END IF
3193 
3194  END IF
3195 
3196  CALL timestop(handle)
3197 
3198  END SUBROUTINE ec_properties
3199 
3200 ! **************************************************************************************************
3201 !> \brief Solve the Harris functional by linear scaling density purification scheme,
3202 !> instead of the diagonalization performed in ec_diag_solver
3203 !>
3204 !> \param qs_env ...
3205 !> \param matrix_ks Harris Kohn-Sham matrix
3206 !> \param matrix_s Overlap matrix in Harris functional basis
3207 !> \par History
3208 !> 09.2020 created
3209 !> \author F.Belleflamme
3210 ! **************************************************************************************************
3211  SUBROUTINE ec_ls_init(qs_env, matrix_ks, matrix_s)
3212  TYPE(qs_environment_type), POINTER :: qs_env
3213  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
3214 
3215  CHARACTER(len=*), PARAMETER :: routinen = 'ec_ls_init'
3216 
3217  INTEGER :: handle, ispin, nspins
3218  TYPE(dft_control_type), POINTER :: dft_control
3219  TYPE(energy_correction_type), POINTER :: ec_env
3220  TYPE(ls_scf_env_type), POINTER :: ls_env
3221 
3222  CALL timeset(routinen, handle)
3223 
3224  CALL get_qs_env(qs_env=qs_env, &
3225  dft_control=dft_control, &
3226  ec_env=ec_env)
3227  nspins = dft_control%nspins
3228  ls_env => ec_env%ls_env
3229 
3230  ! create the matrix template for use in the ls procedures
3231  CALL matrix_ls_create(matrix_ls=ls_env%matrix_s, matrix_qs=matrix_s(1, 1)%matrix, &
3232  ls_mstruct=ls_env%ls_mstruct)
3233 
3234  IF (ALLOCATED(ls_env%matrix_p)) THEN
3235  DO ispin = 1, SIZE(ls_env%matrix_p)
3236  CALL dbcsr_release(ls_env%matrix_p(ispin))
3237  END DO
3238  ELSE
3239  ALLOCATE (ls_env%matrix_p(nspins))
3240  END IF
3241 
3242  DO ispin = 1, nspins
3243  CALL dbcsr_create(ls_env%matrix_p(ispin), template=ls_env%matrix_s, &
3244  matrix_type=dbcsr_type_no_symmetry)
3245  END DO
3246 
3247  ALLOCATE (ls_env%matrix_ks(nspins))
3248  DO ispin = 1, nspins
3249  CALL dbcsr_create(ls_env%matrix_ks(ispin), template=ls_env%matrix_s, &
3250  matrix_type=dbcsr_type_no_symmetry)
3251  END DO
3252 
3253  ! Set up S matrix and needed functions of S
3254  CALL ls_scf_init_matrix_s(matrix_s(1, 1)%matrix, ls_env)
3255 
3256  ! Bring KS matrix from QS to LS form
3257  ! EC KS-matrix already calculated
3258  DO ispin = 1, nspins
3259  CALL matrix_qs_to_ls(matrix_ls=ls_env%matrix_ks(ispin), &
3260  matrix_qs=matrix_ks(ispin, 1)%matrix, &
3261  ls_mstruct=ls_env%ls_mstruct, &
3262  covariant=.true.)
3263  END DO
3264 
3265  CALL timestop(handle)
3266 
3267  END SUBROUTINE ec_ls_init
3268 
3269 ! **************************************************************************************************
3270 !> \brief Solve the Harris functional by linear scaling density purification scheme,
3271 !> instead of the diagonalization performed in ec_diag_solver
3272 !>
3273 !> \param qs_env ...
3274 !> \param matrix_p Harris dentiy matrix, calculated here
3275 !> \param matrix_w Harris energy weighted density matrix, calculated here
3276 !> \param ec_ls_method which purification scheme should be used
3277 !> \par History
3278 !> 12.2019 created [JGH]
3279 !> 08.2020 refactoring [fbelle]
3280 !> \author Fabian Belleflamme
3281 ! **************************************************************************************************
3282 
3283  SUBROUTINE ec_ls_solver(qs_env, matrix_p, matrix_w, ec_ls_method)
3284  TYPE(qs_environment_type), POINTER :: qs_env
3285  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
3286  INTEGER, INTENT(IN) :: ec_ls_method
3287 
3288  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_ls_solver'
3289 
3290  INTEGER :: handle, ispin, nelectron_spin_real, &
3291  nspins
3292  INTEGER, DIMENSION(2) :: nmo
3293  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: wmat
3294  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_ks_deviation
3295  TYPE(dft_control_type), POINTER :: dft_control
3296  TYPE(energy_correction_type), POINTER :: ec_env
3297  TYPE(ls_scf_env_type), POINTER :: ls_env
3298  TYPE(mp_para_env_type), POINTER :: para_env
3299 
3300  CALL timeset(routinen, handle)
3301 
3302  NULLIFY (para_env)
3303  CALL get_qs_env(qs_env, &
3304  dft_control=dft_control, &
3305  para_env=para_env)
3306  nspins = dft_control%nspins
3307  ec_env => qs_env%ec_env
3308  ls_env => ec_env%ls_env
3309 
3310  nmo = 0
3311  CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
3312  IF (nspins == 1) nmo(1) = nmo(1)/2
3313  ls_env%homo_spin(:) = 0.0_dp
3314  ls_env%lumo_spin(:) = 0.0_dp
3315 
3316  ALLOCATE (matrix_ks_deviation(nspins))
3317  DO ispin = 1, nspins
3318  CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_env%matrix_ks(ispin))
3319  CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
3320  END DO
3321 
3322  ! F = S^-1/2 * F * S^-1/2
3323  IF (ls_env%has_s_preconditioner) THEN
3324  DO ispin = 1, nspins
3325  CALL apply_matrix_preconditioner(ls_env%matrix_ks(ispin), "forward", &
3326  ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3327 
3328  CALL dbcsr_filter(ls_env%matrix_ks(ispin), ls_env%eps_filter)
3329  END DO
3330  END IF
3331 
3332  DO ispin = 1, nspins
3333  nelectron_spin_real = ls_env%nelectron_spin(ispin)
3334  IF (ls_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
3335 
3336  SELECT CASE (ec_ls_method)
3337  CASE (ec_matrix_sign)
3338  CALL density_matrix_sign(ls_env%matrix_p(ispin), &
3339  ls_env%mu_spin(ispin), &
3340  ls_env%fixed_mu, &
3341  ls_env%sign_method, &
3342  ls_env%sign_order, &
3343  ls_env%matrix_ks(ispin), &
3344  ls_env%matrix_s, &
3345  ls_env%matrix_s_inv, &
3346  nelectron_spin_real, &
3347  ec_env%eps_default)
3348 
3349  CASE (ec_matrix_trs4)
3350  CALL density_matrix_trs4( &
3351  ls_env%matrix_p(ispin), &
3352  ls_env%matrix_ks(ispin), &
3353  ls_env%matrix_s_sqrt_inv, &
3354  nelectron_spin_real, &
3355  ec_env%eps_default, &
3356  ls_env%homo_spin(ispin), &
3357  ls_env%lumo_spin(ispin), &
3358  ls_env%mu_spin(ispin), &
3359  matrix_ks_deviation=matrix_ks_deviation(ispin), &
3360  dynamic_threshold=ls_env%dynamic_threshold, &
3361  eps_lanczos=ls_env%eps_lanczos, &
3362  max_iter_lanczos=ls_env%max_iter_lanczos)
3363 
3364  CASE (ec_matrix_tc2)
3365  CALL density_matrix_tc2( &
3366  ls_env%matrix_p(ispin), &
3367  ls_env%matrix_ks(ispin), &
3368  ls_env%matrix_s_sqrt_inv, &
3369  nelectron_spin_real, &
3370  ec_env%eps_default, &
3371  ls_env%homo_spin(ispin), &
3372  ls_env%lumo_spin(ispin), &
3373  non_monotonic=ls_env%non_monotonic, &
3374  eps_lanczos=ls_env%eps_lanczos, &
3375  max_iter_lanczos=ls_env%max_iter_lanczos)
3376 
3377  END SELECT
3378 
3379  END DO
3380 
3381  ! de-orthonormalize
3382  IF (ls_env%has_s_preconditioner) THEN
3383  DO ispin = 1, nspins
3384  ! P = S^-1/2 * P_tilde * S^-1/2 (forward)
3385  CALL apply_matrix_preconditioner(ls_env%matrix_p(ispin), "forward", &
3386  ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3387 
3388  CALL dbcsr_filter(ls_env%matrix_p(ispin), ls_env%eps_filter)
3389  END DO
3390  END IF
3391 
3392  ! Closed-shell
3393  IF (nspins == 1) CALL dbcsr_scale(ls_env%matrix_p(1), 2.0_dp)
3394 
3395  IF (ls_env%report_all_sparsities) CALL post_scf_sparsities(ls_env)
3396 
3397  ! ls_scf_dm_to_ks
3398  ! Density matrix from LS to EC
3399  DO ispin = 1, nspins
3400  CALL matrix_ls_to_qs(matrix_qs=matrix_p(ispin, 1)%matrix, &
3401  matrix_ls=ls_env%matrix_p(ispin), &
3402  ls_mstruct=ls_env%ls_mstruct, &
3403  covariant=.false.)
3404  END DO
3405 
3406  wmat => matrix_w(:, 1)
3407  CALL calculate_w_matrix_ls(wmat, ec_env%ls_env)
3408 
3409  ! clean up
3410  CALL dbcsr_release(ls_env%matrix_s)
3411  IF (ls_env%has_s_preconditioner) THEN
3412  CALL dbcsr_release(ls_env%matrix_bs_sqrt)
3413  CALL dbcsr_release(ls_env%matrix_bs_sqrt_inv)
3414  END IF
3415  IF (ls_env%needs_s_inv) THEN
3416  CALL dbcsr_release(ls_env%matrix_s_inv)
3417  END IF
3418  IF (ls_env%use_s_sqrt) THEN
3419  CALL dbcsr_release(ls_env%matrix_s_sqrt)
3420  CALL dbcsr_release(ls_env%matrix_s_sqrt_inv)
3421  END IF
3422 
3423  DO ispin = 1, SIZE(ls_env%matrix_ks)
3424  CALL dbcsr_release(ls_env%matrix_ks(ispin))
3425  END DO
3426  DEALLOCATE (ls_env%matrix_ks)
3427 
3428  DO ispin = 1, nspins
3429  CALL dbcsr_release(matrix_ks_deviation(ispin))
3430  END DO
3431  DEALLOCATE (matrix_ks_deviation)
3432 
3433  CALL timestop(handle)
3434 
3435  END SUBROUTINE ec_ls_solver
3436 
3437 ! **************************************************************************************************
3438 !> \brief Use OT-diagonalziation to obtain density matrix from Harris Kohn-Sham matrix
3439 !> Initial guess of density matrix is either the atomic block initial guess from SCF
3440 !> or the ground-state density matrix. The latter only works if the same basis is used
3441 !>
3442 !> \param qs_env ...
3443 !> \param ec_env ...
3444 !> \param matrix_ks Harris Kohn-Sham matrix
3445 !> \param matrix_s Overlap matrix in Harris functional basis
3446 !> \param matrix_p Harris dentiy matrix, calculated here
3447 !> \param matrix_w Harris energy weighted density matrix, calculated here
3448 !>
3449 !> \par History
3450 !> 09.2020 created
3451 !> \author F.Belleflamme
3452 ! **************************************************************************************************
3453  SUBROUTINE ec_ot_diag_solver(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
3454  TYPE(qs_environment_type), POINTER :: qs_env
3455  TYPE(energy_correction_type), POINTER :: ec_env
3456  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
3457  POINTER :: matrix_ks, matrix_s
3458  TYPE(dbcsr_p_type), DIMENSION(:, :), &
3459  INTENT(INOUT), POINTER :: matrix_p, matrix_w
3460 
3461  CHARACTER(len=*), PARAMETER :: routinen = 'ec_ot_diag_solver'
3462 
3463  INTEGER :: handle, homo, ikind, iounit, ispin, &
3464  max_iter, nao, nkind, nmo, nspins
3465  INTEGER, DIMENSION(2) :: nelectron_spin
3466  REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
3467  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3468  TYPE(cp_blacs_env_type), POINTER :: blacs_env
3469  TYPE(cp_fm_type) :: sv
3470  TYPE(cp_fm_type), POINTER :: mo_coeff
3471  TYPE(cp_logger_type), POINTER :: logger
3472  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
3473  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
3474  TYPE(dft_control_type), POINTER :: dft_control
3475  TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis
3476  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3477  TYPE(mp_para_env_type), POINTER :: para_env
3478  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3479  TYPE(preconditioner_type), POINTER :: local_preconditioner
3480  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3481  TYPE(qs_kind_type), POINTER :: qs_kind
3482  TYPE(qs_rho_type), POINTER :: rho
3483 
3484  CALL timeset(routinen, handle)
3485 
3486  logger => cp_get_default_logger()
3487  iounit = cp_logger_get_default_unit_nr(logger)
3488 
3489  cpassert(ASSOCIATED(qs_env))
3490  cpassert(ASSOCIATED(ec_env))
3491  cpassert(ASSOCIATED(matrix_ks))
3492  cpassert(ASSOCIATED(matrix_s))
3493  cpassert(ASSOCIATED(matrix_p))
3494  cpassert(ASSOCIATED(matrix_w))
3495 
3496  CALL get_qs_env(qs_env=qs_env, &
3497  atomic_kind_set=atomic_kind_set, &
3498  blacs_env=blacs_env, &
3499  dft_control=dft_control, &
3500  nelectron_spin=nelectron_spin, &
3501  para_env=para_env, &
3502  particle_set=particle_set, &
3503  qs_kind_set=qs_kind_set)
3504  nspins = dft_control%nspins
3505 
3506  ! Maximum number of OT iterations for diagonalization
3507  max_iter = 200
3508 
3509  ! If linear scaling, need to allocate and init MO set
3510  ! set it to qs_env%mos
3511  IF (dft_control%qs_control%do_ls_scf) THEN
3512  CALL ec_mos_init(qs_env, matrix_s(1, 1)%matrix)
3513  END IF
3514 
3515  CALL get_qs_env(qs_env, mos=mos)
3516 
3517  ! Inital guess to use
3518  NULLIFY (p_rmpv)
3519 
3520  ! Using ether ground-state DM or ATOMIC GUESS requires
3521  ! Harris functional to use the same basis set
3522  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
3523  CALL uppercase(ec_env%basis)
3524  ! Harris basis only differs from ground-state basis if explicitly added
3525  ! thus only two cases that need to be tested
3526  ! 1) explicit Harris basis present?
3527  IF (ec_env%basis == "HARRIS") THEN
3528  DO ikind = 1, nkind
3529  qs_kind => qs_kind_set(ikind)
3530  ! Basis sets of ground-state
3531  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
3532  ! Basis sets of energy correction
3533  CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
3534 
3535  IF (basis_set%name .NE. harris_basis%name) THEN
3536  cpabort("OT-Diag initial guess: Harris and ground state need to use the same basis")
3537  END IF
3538  END DO
3539  END IF
3540  ! 2) Harris uses MAOs
3541  IF (ec_env%mao) THEN
3542  cpabort("OT-Diag initial guess: not implemented for MAOs")
3543  END IF
3544 
3545  ! Initital guess obtained for OT Diagonalization
3546  SELECT CASE (ec_env%ec_initial_guess)
3547  CASE (ec_ot_atomic)
3548 
3549  p_rmpv => matrix_p(:, 1)
3550 
3551  CALL calculate_atomic_block_dm(p_rmpv, matrix_s(1, 1)%matrix, atomic_kind_set, qs_kind_set, &
3552  nspins, nelectron_spin, iounit, para_env)
3553 
3554  CASE (ec_ot_gs)
3555 
3556  CALL get_qs_env(qs_env, rho=rho)
3557  CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3558  p_rmpv => rho_ao(:, 1)
3559 
3560  CASE DEFAULT
3561  cpabort("Unknown inital guess for OT-Diagonalization (Harris functional)")
3562  END SELECT
3563 
3564  DO ispin = 1, nspins
3565  CALL get_mo_set(mo_set=mos(ispin), &
3566  mo_coeff=mo_coeff, &
3567  nmo=nmo, &
3568  nao=nao, &
3569  homo=homo)
3570 
3571  ! Calculate first MOs
3572  CALL cp_fm_set_all(mo_coeff, 0.0_dp)
3573  CALL cp_fm_init_random(mo_coeff, nmo)
3574 
3575  CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
3576  ! multiply times PS
3577  ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
3578  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff, sv, nmo)
3579  CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
3580  CALL cp_fm_release(sv)
3581  ! and ortho the result
3582  ! If DFBT or SE, then needs has_unit_metrix option
3583  CALL make_basis_sm(mo_coeff, nmo, matrix_s(1, 1)%matrix)
3584  END DO
3585 
3586  ! Preconditioner
3587  NULLIFY (local_preconditioner)
3588  ALLOCATE (local_preconditioner)
3589  CALL init_preconditioner(local_preconditioner, para_env=para_env, &
3590  blacs_env=blacs_env)
3591  DO ispin = 1, nspins
3592  CALL make_preconditioner(local_preconditioner, &
3593  precon_type=ot_precond_full_single_inverse, &
3594  solver_type=ot_precond_solver_default, &
3595  matrix_h=matrix_ks(ispin, 1)%matrix, &
3596  matrix_s=matrix_s(ispin, 1)%matrix, &
3597  convert_precond_to_dbcsr=.true., &
3598  mo_set=mos(ispin), energy_gap=0.2_dp)
3599 
3600  CALL get_mo_set(mos(ispin), &
3601  mo_coeff=mo_coeff, &
3602  eigenvalues=eigenvalues, &
3603  nmo=nmo, &
3604  homo=homo)
3605  CALL ot_eigensolver(matrix_h=matrix_ks(ispin, 1)%matrix, &
3606  matrix_s=matrix_s(1, 1)%matrix, &
3607  matrix_c_fm=mo_coeff, &
3608  preconditioner=local_preconditioner, &
3609  eps_gradient=ec_env%eps_default, &
3610  iter_max=max_iter, &
3611  silent=.false.)
3612  CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin, 1)%matrix, &
3613  evals_arg=eigenvalues, do_rotation=.true.)
3614 
3615  ! Deallocate preconditioner
3616  CALL destroy_preconditioner(local_preconditioner)
3617  DEALLOCATE (local_preconditioner)
3618 
3619  !fm->dbcsr
3620  CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
3621  mos(ispin)%mo_coeff_b)
3622  END DO
3623 
3624  ! Calculate density matrix from MOs
3625  DO ispin = 1, nspins
3626  CALL calculate_density_matrix(mos(ispin), matrix_p(ispin, 1)%matrix)
3627 
3628  CALL calculate_w_matrix(mos(ispin), matrix_w(ispin, 1)%matrix)
3629  END DO
3630 
3631  ! Get rid of MO environment again
3632  IF (dft_control%qs_control%do_ls_scf) THEN
3633  DO ispin = 1, nspins
3634  CALL deallocate_mo_set(mos(ispin))
3635  END DO
3636  IF (ASSOCIATED(qs_env%mos)) THEN
3637  DO ispin = 1, SIZE(qs_env%mos)
3638  CALL deallocate_mo_set(qs_env%mos(ispin))
3639  END DO
3640  DEALLOCATE (qs_env%mos)
3641  END IF
3642  END IF
3643 
3644  CALL timestop(handle)
3645 
3646  END SUBROUTINE ec_ot_diag_solver
3647 
3648 END MODULE energy_corrections
3649 
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
subroutine uppercase(string)
...
Definition: dumpdcd.F:1376
Contains ADMM methods which only require the density matrix.
subroutine, public admm_dm_calc_rho_aux(qs_env)
Entry methods: Calculates auxiliary density matrix from primary one.
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
subroutine, public admm_mo_calc_rho_aux(qs_env)
...
Definition: admm_methods.F:149
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public belleflamme2023
Definition: bibliography.F:43
Handles all functions related to the CELL.
Definition: cell_types.F:15
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition: core_ae.F:14
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Definition: core_ae.F:84
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition: core_ppl.F:18
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Definition: core_ppl.F:95
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: core_ppnl.F:15
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
Definition: core_ppnl.F:89
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 copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
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_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
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,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
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_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_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_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
Definition: cp_fm_types.F:452
subroutine, public cp_fm_to_fm_triangular(msource, mtarget, uplo)
copy just a triangular matrix
Definition: cp_fm_types.F:1428
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 ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
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,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
set of type/routines to handle the storage of results in force_envs
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
lower level routines for linear scaling SCF
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
subroutine, public apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv)
apply a preconditioner either forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs)) backward (rest...
subroutine, public density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, matrix_s_sqrt_inv)
compute the density matrix with a trace that is close to nelectron. take a mu as input,...
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged)
compute the density matrix using a trace-resetting algorithm
subroutine, public density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, non_monotonic, eps_lanczos, max_iter_lanczos)
compute the density matrix using a non monotonic trace conserving algorithm based on SIAM DOI....
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
Definition: dm_ls_scf_qs.F:15
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
Definition: dm_ls_scf_qs.F:307
subroutine, public matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
create a matrix for use (and as a template) in ls based on a qs template
Definition: dm_ls_scf_qs.F:97
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
Definition: dm_ls_scf_qs.F:213
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Routines for a linear scaling quickstep SCF run based on the density matrix.
Definition: dm_ls_scf.F:15
subroutine, public post_scf_sparsities(ls_scf_env)
Report on the sparsities of various interesting matrices.
Definition: dm_ls_scf.F:1027
subroutine, public calculate_w_matrix_ls(matrix_w, ls_scf_env)
Compute matrix_w as needed for the forces.
Definition: dm_ls_scf.F:1184
Calculates the energy contribution and the mo_derivative of a static electric field (nonperiodic)
subroutine, public ec_efield_local_operator(qs_env, ec_env, calculate_forces)
...
subroutine, public ec_efield_integrals(qs_env, ec_env, rpoint)
...
Types needed for a for a Energy Correction.
Definition: ec_env_types.F:14
Routines used for Harris functional Kohn-Sham calculation.
Definition: ec_methods.F:15
subroutine, public ec_mos_init(qs_env, matrix_s)
Allocate and initiate molecular orbitals environment.
Definition: ec_methods.F:162
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
Definition: ec_methods.F:88
Routines for an energy correction on top of a Kohn-Sham calculation.
subroutine, public energy_correction(qs_env, ec_init, calculate_forces)
Energy Correction to a Kohn-Sham simulation Available energy corrections: (1) Harris energy functiona...
Definition of the atomic potential types.
Routines to calculate EXX in RPA and energy correction methods.
Definition: hfx_exx.F:16
subroutine, public calculate_exx(qs_env, unit_nr, hfx_sections, x_data, do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, E_ex_from_GW, E_admm_from_GW, t3)
...
Definition: hfx_exx.F:106
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
Definition: hfx_exx.F:325
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ec_functional_harris
integer, parameter, public ec_functional_dc
integer, parameter, public vdw_pairpot_dftd3
integer, parameter, public ec_ot_atomic
integer, parameter, public ec_ot_diag
integer, parameter, public ec_ot_gs
integer, parameter, public ot_precond_solver_default
integer, parameter, public ot_precond_full_single_inverse
integer, parameter, public ec_matrix_tc2
integer, parameter, public ec_diagonalization
integer, parameter, public ec_matrix_trs4
integer, parameter, public ec_matrix_sign
integer, parameter, public xc_vdw_fun_pairpot
integer, parameter, public vdw_pairpot_dftd3bj
objects that represent the structure of input sections and the data contained in an input section
integer function, public section_get_ival(section_vals, keyword_name)
...
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
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
logical function, public section_get_lval(section_vals, keyword_name)
...
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
integer, parameter, public default_path_length
Definition: kinds.F:58
Calculate MAO's and analyze wavefunctions.
Definition: mao_basis.F:15
subroutine, public mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, molecular, max_iter, eps_grad, nmao_external, eps1_mao, iolevel, unit_nr)
...
Definition: mao_basis.F:75
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
Interface to the message passing library MPI.
Define the data structure for the molecule information.
Calculates the moment integrals <a|r^m|b>
Definition: moments_utils.F:15
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition: moments_utils.F:61
Define the data structure for the particle information.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public pascal
Definition: physcon.F:174
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
real(kind=dp), parameter, public debye
Definition: physcon.F:201
types of preconditioners
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
subroutine, public destroy_preconditioner(preconditioner_env)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
...
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
subroutine, public calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, ounit, para_env)
returns a block diagonal density matrix. Blocks correspond to the atomic densities.
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
Calculation of the energies concerning the core charge distribution.
subroutine, public calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, E_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
collects routines that calculate density matrices
Calculation of dispersion using pair potentials.
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition of disperson types for DFT calculations.
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.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Calculation of kinetic energy matrix and forces.
Definition: qs_kinetic.F:15
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition: qs_kinetic.F:101
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition: qs_mo_methods.F:87
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
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
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:558
Define the neighbor list data types and the corresponding functionality.
Generate the atomic neighbor lists.
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:120
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
Definition: qs_vxc.F:16
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition: qs_vxc.F:98
Calculate the CPKS equation and the resulting forces.
subroutine, public response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
...
subroutine, public response_calculation(qs_env, ec_env)
Initializes solver of linear response equation for energy correction.
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's)
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
generate the tasks lists used by collocate and integrate routines
subroutine, public generate_qs_task_list(ks_env, task_list, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
types for task lists
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
subroutine, public write_stress_tensor_components(virial, iw, cell)
...
subroutine, public write_stress_tensor(pv_virial, iw, cell, numerical)
Print stress tensor to output file.
pure real(kind=dp) function, public one_third_sum_diag(a)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Definition: virial_types.F:67
Interface for Voronoi Integration and output of BQB files.
subroutine, public entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
Does a Voronoi integration of density or stores the density to compressed BQB format.