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