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