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