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