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