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