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