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