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