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