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