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