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