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