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