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