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