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