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