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