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