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