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