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