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