Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate Energy Decomposition analysis
10 : !> \par History
11 : !> 07.2023 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ed_analysis
15 : USE admm_types, ONLY: admm_type
16 : USE atomic_kind_types, ONLY: atomic_kind_type
17 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
18 : gto_basis_set_type
19 : USE bibliography, ONLY: Eriksen2020,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type
22 : USE core_ae, ONLY: build_erfc
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, &
26 : dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
27 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
28 : dbcsr_reserve_diag_blocks
29 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
30 : cp_dbcsr_sm_fm_multiply
31 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE cp_fm_types, ONLY: &
35 : cp_fm_create, cp_fm_get_diag, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, &
36 : cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_set_submatrix, cp_fm_to_fm, &
37 : cp_fm_type
38 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
39 : USE hartree_local_types, ONLY: ecoul_1center_type
40 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
41 : USE iao_analysis, ONLY: iao_calculate_dmat,&
42 : iao_charges,&
43 : iao_wfn_analysis,&
44 : print_center_spread
45 : USE iao_types, ONLY: iao_env_type,&
46 : iao_set_default
47 : USE input_constants, ONLY: do_admm_aux_exch_func_none
48 : USE input_section_types, ONLY: section_vals_get,&
49 : section_vals_get_subs_vals,&
50 : section_vals_type,&
51 : section_vals_val_get
52 : USE kinds, ONLY: dp
53 : USE mathconstants, ONLY: pi
54 : USE message_passing, ONLY: mp_comm_type,&
55 : mp_para_env_type
56 : USE mulliken, ONLY: mulliken_charges
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE particle_methods, ONLY: get_particle_set
59 : USE particle_types, ONLY: particle_type
60 : USE physcon, ONLY: angstrom
61 : USE pw_env_types, ONLY: pw_env_get,&
62 : pw_env_type
63 : USE pw_methods, ONLY: pw_axpy,&
64 : pw_scale,&
65 : pw_transfer
66 : USE pw_poisson_methods, ONLY: pw_poisson_solve
67 : USE pw_poisson_types, ONLY: pw_poisson_type
68 : USE pw_pool_types, ONLY: pw_pool_type
69 : USE pw_types, ONLY: pw_c1d_gs_type,&
70 : pw_r3d_rs_type
71 : USE qs_collocate_density, ONLY: calculate_rho_core
72 : USE qs_core_energies, ONLY: calculate_ecore_alpha,&
73 : calculate_ecore_overlap,&
74 : calculate_ecore_self
75 : USE qs_core_hamiltonian, ONLY: core_matrices
76 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
77 : USE qs_dispersion_types, ONLY: qs_dispersion_type
78 : USE qs_energy_types, ONLY: qs_energy_type
79 : USE qs_environment_types, ONLY: get_qs_env,&
80 : qs_environment_type
81 : USE qs_gcp_method, ONLY: calculate_gcp_pairpot
82 : USE qs_gcp_types, ONLY: qs_gcp_type
83 : USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
84 : integrate_v_gaussian_rspace,&
85 : integrate_v_rspace
86 : USE qs_kind_types, ONLY: get_qs_kind,&
87 : qs_kind_type
88 : USE qs_ks_atom, ONLY: update_ks_atom
89 : USE qs_ks_types, ONLY: qs_ks_env_type
90 : USE qs_local_rho_types, ONLY: local_rho_type
91 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
92 : USE qs_mo_types, ONLY: deallocate_mo_set,&
93 : duplicate_mo_set,&
94 : get_mo_set,&
95 : mo_set_type
96 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
97 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
98 : USE qs_rho_atom_types, ONLY: rho_atom_type,&
99 : zero_rho_atom_integrals
100 : USE qs_rho_types, ONLY: qs_rho_get,&
101 : qs_rho_type
102 : USE qs_vxc, ONLY: qs_xc_density
103 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
104 : USE xc_derivatives, ONLY: xc_functionals_get_needs
105 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
106 : #include "./base/base_uses.f90"
107 :
108 : IMPLICIT NONE
109 : PRIVATE
110 :
111 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ed_analysis'
112 :
113 : PUBLIC :: edmf_analysis
114 :
115 : ! **************************************************************************************************
116 :
117 : CONTAINS
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param qs_env ...
122 : !> \param input_section ...
123 : !> \param unit_nr ...
124 : ! **************************************************************************************************
125 58 : SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr)
126 : TYPE(qs_environment_type), POINTER :: qs_env
127 : TYPE(section_vals_type), POINTER :: input_section
128 : INTEGER, INTENT(IN) :: unit_nr
129 :
130 : CHARACTER(len=*), PARAMETER :: routineN = 'edmf_analysis'
131 :
132 : INTEGER :: handle, iatom, ikind, iorb, ispin, jorb, &
133 : nao, natom, nimages, nkind, no, norb, &
134 : nref, nspin
135 : INTEGER, DIMENSION(2) :: nocc
136 58 : INTEGER, DIMENSION(:), POINTER :: refbas_blk_sizes
137 : LOGICAL :: detailed_ener, do_hfx, ewald_correction, explicit, gapw, gapw_xc, &
138 : ref_orb_canonical, skip_localize, uniform_occupation, uocc
139 : REAL(KIND=dp) :: ateps, checksum, e1, e2, e_pot, ealpha, &
140 : ecc, egcp, ehfx, ekts, evdw, focc, &
141 : sum_energy
142 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: amval, atcore, ate1h, ate1xc, atecc, &
143 58 : ateks, atener, atewald, odiag
144 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atdet, atmul, mcharge, mweight
145 58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: bcenter
146 58 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
147 : TYPE(cell_type), POINTER :: cell
148 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
149 : TYPE(cp_fm_type) :: cvec, cvec2, smo
150 58 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: c_iao_coef, ciao, fij_mat, orb_weight, &
151 58 : rotmat
152 : TYPE(cp_fm_type), POINTER :: cloc, moref
153 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
154 : TYPE(dbcsr_p_type) :: dve_mat
155 58 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: exc_mat, ks_mat, vhxc_mat
156 58 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: core_mat, matrix_h, matrix_hfx, &
157 58 : matrix_ks, matrix_p, matrix_s, matrix_t
158 58 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: math, matp
159 : TYPE(dbcsr_type) :: dkmat, dmat
160 : TYPE(dbcsr_type), POINTER :: smat
161 : TYPE(dft_control_type), POINTER :: dft_control
162 58 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: ref_basis_set_list
163 : TYPE(gto_basis_set_type), POINTER :: refbasis
164 : TYPE(iao_env_type) :: iao_env
165 58 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_loc
166 : TYPE(mp_comm_type) :: group
167 : TYPE(mp_para_env_type), POINTER :: para_env
168 58 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
169 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
170 : TYPE(qs_energy_type), POINTER :: energy
171 : TYPE(qs_gcp_type), POINTER :: gcp_env
172 58 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
173 : TYPE(qs_kind_type), POINTER :: qs_kind
174 : TYPE(qs_rho_type), POINTER :: rho
175 : TYPE(section_vals_type), POINTER :: hfx_sections, input
176 :
177 58 : CALL section_vals_get(input_section, explicit=explicit)
178 58 : IF (.NOT. explicit) RETURN
179 :
180 30 : CALL timeset(routineN, handle)
181 :
182 30 : IF (unit_nr > 0) THEN
183 : WRITE (UNIT=unit_nr, FMT="(/,4(T2,A))") &
184 15 : "!-----------------------------------------------------------------------------!", &
185 15 : "! ENERGY DECOMPOSITION ANALYSIS !", &
186 15 : "! Janus J Eriksen, JCP 153 214109 (2020) !", &
187 30 : "!-----------------------------------------------------------------------------!"
188 : END IF
189 30 : CALL cite_reference(Eriksen2020)
190 :
191 : ! input options
192 30 : CALL section_vals_val_get(input_section, "REFERENCE_ORB_CANONICAL", l_val=ref_orb_canonical)
193 30 : CALL section_vals_val_get(input_section, "DETAILED_ENERGY", l_val=detailed_ener)
194 30 : CALL section_vals_val_get(input_section, "SKIP_LOCALIZATION", l_val=skip_localize)
195 30 : CALL section_vals_val_get(input_section, "EWALD_ALPHA_PARAMETER", r_val=ealpha)
196 30 : ewald_correction = (ealpha /= 0.0_dp)
197 : ! k-points?
198 30 : CALL get_qs_env(qs_env, dft_control=dft_control)
199 30 : nimages = dft_control%nimages
200 30 : IF (nimages > 1) THEN
201 0 : IF (unit_nr > 0) THEN
202 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
203 0 : "K-Points: MAO's determined and analyzed using Gamma-Point only."
204 : END IF
205 : END IF
206 30 : gapw = dft_control%qs_control%gapw
207 30 : gapw_xc = dft_control%qs_control%gapw_xc
208 30 : IF (ewald_correction) THEN
209 2 : IF (gapw .OR. gapw_xc) THEN
210 0 : CALL cp_warn(__LOCATION__, "Ewald Correction for GAPW and GAPW_XC not available.")
211 0 : CPABORT("Ewald Correction")
212 : END IF
213 : END IF
214 30 : CALL get_qs_env(qs_env, input=input)
215 30 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
216 30 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
217 :
218 30 : CALL get_qs_env(qs_env, mos=mos)
219 30 : nspin = dft_control%nspins
220 122 : ALLOCATE (mos_loc(nspin))
221 :
222 : ! do we have a uniform occupation
223 30 : uocc = .TRUE.
224 62 : DO ispin = 1, nspin
225 32 : CALL get_mo_set(mos(ispin), uniform_occupation=uniform_occupation)
226 62 : IF (.NOT. uniform_occupation) uocc = .FALSE.
227 : END DO
228 30 : IF (unit_nr > 0) THEN
229 15 : IF (uocc) THEN
230 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
231 14 : "MO's have a uniform occupation pattern"
232 : ELSE
233 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
234 1 : "MO's have varying occupations"
235 : END IF
236 : END IF
237 :
238 : ! perform IAO analysis
239 30 : CALL iao_set_default(iao_env)
240 30 : iao_env%do_iao = .TRUE.
241 30 : iao_env%do_charges = .TRUE.
242 30 : IF (skip_localize) THEN
243 2 : iao_env%do_bondorbitals = .FALSE.
244 2 : iao_env%do_center = .FALSE.
245 : ELSE
246 28 : iao_env%do_bondorbitals = .TRUE.
247 28 : iao_env%do_center = .TRUE.
248 : END IF
249 30 : iao_env%eps_occ = 1.0E-4_dp
250 30 : CALL get_qs_env(qs_env, cell=cell)
251 36 : iao_env%pos_periodic = .NOT. ALL(cell%perd == 0)
252 30 : no = 0
253 30 : nocc = 0
254 62 : DO ispin = 1, nspin
255 32 : CALL duplicate_mo_set(mos_loc(ispin), mos(ispin))
256 32 : CALL get_mo_set(mos_loc(ispin), nmo=norb)
257 32 : no = MAX(no, norb)
258 32 : nocc(ispin) = norb
259 62 : IF (ref_orb_canonical) THEN
260 32 : CALL get_mo_set(mos_loc(ispin), mo_coeff=moref)
261 32 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
262 : CALL calculate_subspace_eigenvalues(moref, matrix_ks(ispin)%matrix, &
263 32 : do_rotation=.TRUE.)
264 : END IF
265 : END DO
266 120 : ALLOCATE (bcenter(5, no, nspin))
267 1154 : bcenter = 0.0_dp
268 : CALL iao_wfn_analysis(qs_env, iao_env, unit_nr, &
269 30 : c_iao_coef=c_iao_coef, mos=mos_loc, bond_centers=bcenter)
270 :
271 : ! output bond centers
272 30 : IF (iao_env%do_bondorbitals .AND. iao_env%do_center) THEN
273 28 : CALL print_center_spread(bcenter, nocc, iounit=unit_nr)
274 : END IF
275 :
276 : ! Calculate orbital rotation matrix
277 30 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
278 30 : smat => matrix_s(1)%matrix
279 122 : ALLOCATE (rotmat(nspin))
280 62 : DO ispin = 1, nspin
281 32 : CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
282 32 : CALL get_mo_set(mos(ispin), mo_coeff=moref)
283 32 : CALL cp_fm_get_info(cloc, nrow_global=nao, ncol_global=norb)
284 32 : CALL cp_fm_create(smo, cloc%matrix_struct)
285 32 : NULLIFY (fm_struct)
286 : CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
287 32 : template_fmstruct=cloc%matrix_struct)
288 32 : CALL cp_fm_create(rotmat(ispin), fm_struct)
289 32 : CALL cp_fm_struct_release(fm_struct)
290 : ! ROTMAT = Cref(T)*S*Cloc
291 32 : CALL cp_dbcsr_sm_fm_multiply(smat, cloc, smo, ncol=norb)
292 : CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, moref, &
293 32 : smo, 0.0_dp, rotmat(ispin))
294 94 : CALL cp_fm_release(smo)
295 : END DO
296 :
297 : ! calculate occupation matrix
298 30 : IF (.NOT. uocc) THEN
299 6 : ALLOCATE (fij_mat(nspin))
300 4 : DO ispin = 1, nspin
301 2 : CALL cp_fm_create(fij_mat(ispin), rotmat(ispin)%matrix_struct)
302 2 : CALL cp_fm_set_all(fij_mat(ispin), 0.0_dp)
303 2 : CALL cp_fm_create(smo, rotmat(ispin)%matrix_struct)
304 : ! fii = f
305 2 : CALL get_mo_set(mos(ispin), nmo=norb, occupation_numbers=occupation_numbers)
306 54 : DO iorb = 1, norb
307 54 : CALL cp_fm_set_element(fij_mat(ispin), iorb, iorb, occupation_numbers(iorb))
308 : END DO
309 : ! fij = U(T)*f*U
310 : CALL parallel_gemm('N', 'N', norb, norb, norb, 1.0_dp, fij_mat(ispin), &
311 2 : rotmat(ispin), 0.0_dp, smo)
312 : CALL parallel_gemm('T', 'N', norb, norb, norb, 1.0_dp, rotmat(ispin), &
313 2 : smo, 0.0_dp, fij_mat(ispin))
314 6 : CALL cp_fm_release(smo)
315 : END DO
316 : END IF
317 :
318 : ! localized orbitals in IAO basis => CIAO
319 92 : ALLOCATE (ciao(nspin))
320 62 : DO ispin = 1, nspin
321 32 : CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
322 32 : CALL cp_fm_get_info(cloc, nrow_global=nao, ncol_global=norb)
323 32 : CALL cp_fm_get_info(c_iao_coef(ispin), ncol_global=nref)
324 32 : CALL cp_fm_create(smo, cloc%matrix_struct)
325 32 : NULLIFY (fm_struct)
326 : CALL cp_fm_struct_create(fm_struct, nrow_global=nref, &
327 32 : template_fmstruct=cloc%matrix_struct)
328 32 : CALL cp_fm_create(ciao(ispin), fm_struct)
329 32 : CALL cp_fm_struct_release(fm_struct)
330 : ! CIAO = A(T)*S*C
331 32 : CALL cp_dbcsr_sm_fm_multiply(smat, cloc, smo, ncol=norb)
332 : CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, c_iao_coef(ispin), &
333 32 : smo, 0.0_dp, ciao(ispin))
334 94 : CALL cp_fm_release(smo)
335 : END DO
336 :
337 : ! Reference basis set
338 30 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
339 30 : nkind = SIZE(qs_kind_set)
340 148 : ALLOCATE (ref_basis_set_list(nkind))
341 88 : DO ikind = 1, nkind
342 58 : qs_kind => qs_kind_set(ikind)
343 58 : NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
344 58 : NULLIFY (refbasis)
345 58 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type="MIN")
346 88 : IF (ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
347 : END DO
348 30 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, natom=natom)
349 90 : ALLOCATE (refbas_blk_sizes(natom))
350 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=refbas_blk_sizes, &
351 30 : basis=ref_basis_set_list)
352 :
353 : ! Atomic orbital weights
354 92 : ALLOCATE (orb_weight(nspin))
355 90 : ALLOCATE (mcharge(natom, 1))
356 62 : DO ispin = 1, nspin
357 32 : CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
358 32 : NULLIFY (fm_struct)
359 : CALL cp_fm_struct_create(fm_struct, nrow_global=natom, &
360 32 : template_fmstruct=cloc%matrix_struct)
361 32 : CALL cp_fm_create(orb_weight(ispin), fm_struct)
362 32 : CALL cp_fm_struct_release(fm_struct)
363 62 : CALL cp_fm_set_all(orb_weight(ispin), 0.0_dp)
364 : END DO
365 : !
366 30 : CALL dbcsr_get_info(smat, distribution=dbcsr_dist)
367 : CALL dbcsr_create(matrix=dmat, name="DMAT", &
368 : dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
369 30 : row_blk_size=refbas_blk_sizes, col_blk_size=refbas_blk_sizes)
370 30 : CALL dbcsr_reserve_diag_blocks(dmat)
371 : !
372 30 : NULLIFY (fm_struct)
373 : CALL cp_fm_struct_create(fm_struct, ncol_global=1, &
374 30 : template_fmstruct=ciao(1)%matrix_struct)
375 30 : CALL cp_fm_create(cvec, fm_struct)
376 30 : CALL cp_fm_create(cvec2, fm_struct)
377 30 : CALL cp_fm_struct_release(fm_struct)
378 : !
379 62 : DO ispin = 1, nspin
380 : CALL get_mo_set(mos_loc(ispin), &
381 : mo_coeff=cloc, nmo=norb, &
382 32 : occupation_numbers=occupation_numbers)
383 62 : IF (uocc) THEN
384 : ! uniform occupation
385 158 : DO iorb = 1, norb
386 128 : CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
387 128 : focc = occupation_numbers(iorb)
388 128 : CALL dbcsr_set(dmat, 0.0_dp)
389 128 : CALL cp_dbcsr_plus_fm_fm_t(dmat, cvec, ncol=1, alpha=focc)
390 128 : CALL iao_charges(dmat, mcharge(:, 1))
391 : CALL cp_fm_set_submatrix(orb_weight(ispin), mcharge, start_row=1, &
392 128 : start_col=iorb, n_rows=natom, n_cols=1)
393 560 : checksum = SUM(mcharge(:, 1))
394 158 : IF (ABS(focc - checksum) > 1.E-6_dp) THEN
395 0 : CALL cp_warn(__LOCATION__, "Sum of atomic orbital weights is incorrect")
396 0 : IF (unit_nr > 0) THEN
397 : WRITE (UNIT=unit_nr, FMT="(T2,A,F10.6,T40,A,F10.6)") &
398 0 : "Orbital occupation:", focc, &
399 0 : "Sum of atomic orbital weights:", checksum
400 : END IF
401 : END IF
402 : END DO
403 : ELSE
404 : ! non-diagonal occupation matrix
405 6 : ALLOCATE (odiag(norb))
406 2 : CALL cp_fm_get_diag(fij_mat(ispin), odiag)
407 54 : DO iorb = 1, norb
408 52 : IF (odiag(iorb) < 1.E-8_dp) CYCLE
409 44 : CALL dbcsr_set(dmat, 0.0_dp)
410 44 : CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
411 1188 : DO jorb = 1, norb
412 1144 : CALL cp_fm_get_element(fij_mat(ispin), iorb, jorb, focc)
413 1144 : IF (focc < 1.E-12_dp) CYCLE
414 534 : CALL cp_fm_to_fm(ciao(ispin), cvec2, ncol=1, source_start=jorb, target_start=1)
415 1722 : CALL cp_dbcsr_plus_fm_fm_t(dmat, cvec, cvec2, 1, alpha=focc, symmetry_mode=1)
416 : END DO
417 44 : CALL iao_charges(dmat, mcharge(:, 1))
418 396 : checksum = SUM(mcharge(:, 1))
419 44 : focc = odiag(iorb)
420 44 : IF (ABS(focc - checksum) > 1.E-6_dp) THEN
421 0 : CALL cp_warn(__LOCATION__, "Sum of atomic orbital weights is incorrect")
422 0 : IF (unit_nr > 0) THEN
423 : WRITE (UNIT=unit_nr, FMT="(T2,A,F10.6,T40,A,F10.6)") &
424 0 : "Orbital occupation:", focc, &
425 0 : "Sum of atomic orbital weights:", checksum
426 : END IF
427 : END IF
428 396 : mcharge(:, 1) = mcharge(:, 1)/focc
429 : CALL cp_fm_set_submatrix(orb_weight(ispin), mcharge, start_row=1, &
430 54 : start_col=iorb, n_rows=natom, n_cols=1)
431 : END DO
432 2 : DEALLOCATE (odiag)
433 : END IF
434 : END DO
435 30 : DEALLOCATE (mcharge)
436 30 : CALL dbcsr_release(dmat)
437 30 : CALL cp_fm_release(cvec)
438 30 : CALL cp_fm_release(cvec2)
439 :
440 30 : CALL get_qs_env(qs_env, para_env=para_env)
441 30 : group = para_env
442 : ! energy arrays
443 180 : ALLOCATE (atener(natom), ateks(natom), atecc(natom), atcore(natom))
444 120 : ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
445 136 : atener = 0.0_dp
446 136 : ateks = 0.0_dp
447 136 : atecc = 0.0_dp
448 136 : atcore = 0.0_dp
449 136 : ate1xc = 0.0_dp
450 136 : ate1h = 0.0_dp
451 136 : atewald = 0.0_dp
452 30 : IF (detailed_ener) THEN
453 30 : ALLOCATE (atdet(natom, 10), atmul(natom, 10), amval(natom))
454 246 : atdet = 0.0_dp
455 246 : atmul = 0.0_dp
456 : END IF
457 : ! atom dependent density matrix
458 30 : CALL dbcsr_create(dkmat, template=smat)
459 30 : CALL dbcsr_copy(dkmat, smat)
460 30 : CALL dbcsr_set(dkmat, 0.0_dp)
461 : ! KS matrix + correction
462 30 : CALL get_qs_env(qs_env, matrix_h=matrix_h, matrix_ks=matrix_ks, kinetic=matrix_t)
463 244 : ALLOCATE (ks_mat(nspin), core_mat(1), vhxc_mat(nspin), exc_mat(1))
464 62 : DO ispin = 1, nspin
465 32 : ALLOCATE (ks_mat(ispin)%matrix)
466 32 : CALL dbcsr_create(ks_mat(ispin)%matrix, template=matrix_h(1)%matrix)
467 32 : CALL dbcsr_copy(ks_mat(ispin)%matrix, matrix_h(1)%matrix)
468 32 : CALL dbcsr_add(ks_mat(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
469 32 : CALL dbcsr_scale(ks_mat(ispin)%matrix, 0.5_dp)
470 : !
471 32 : ALLOCATE (vhxc_mat(ispin)%matrix)
472 32 : CALL dbcsr_create(vhxc_mat(ispin)%matrix, template=smat)
473 32 : CALL dbcsr_copy(vhxc_mat(ispin)%matrix, smat)
474 62 : CALL dbcsr_set(vhxc_mat(ispin)%matrix, 0.0_dp)
475 : END DO
476 : !
477 30 : ALLOCATE (core_mat(1)%matrix)
478 30 : CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
479 30 : CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
480 30 : CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
481 : !
482 30 : ALLOCATE (exc_mat(1)%matrix)
483 30 : CALL dbcsr_create(exc_mat(1)%matrix, template=smat)
484 30 : CALL dbcsr_copy(exc_mat(1)%matrix, smat)
485 30 : CALL dbcsr_set(exc_mat(1)%matrix, 0.0_dp)
486 : !
487 30 : CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
488 : !
489 30 : CALL get_qs_env(qs_env, rho=rho)
490 30 : CALL qs_rho_get(rho, rho_ao=matrix_p)
491 30 : math(1:1, 1:1) => core_mat(1:1)
492 30 : matp(1:nspin, 1:1) => matrix_p(1:nspin)
493 30 : CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
494 30 : CALL group%sum(atcore)
495 : !
496 30 : IF (ewald_correction) THEN
497 2 : ALLOCATE (dve_mat%matrix)
498 2 : CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
499 2 : CALL dbcsr_copy(dve_mat%matrix, matrix_h(1)%matrix)
500 2 : CALL dbcsr_set(dve_mat%matrix, 0.0_dp)
501 2 : CALL vh_ewald_correction(qs_env, ealpha, dve_mat, atewald)
502 2 : CALL group%sum(atewald)
503 : END IF
504 : !
505 62 : DO ispin = 1, nspin
506 32 : CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
507 62 : CALL dbcsr_add(ks_mat(ispin)%matrix, core_mat(1)%matrix, 1.0_dp, -0.5_dp)
508 : END DO
509 : !
510 30 : IF (detailed_ener .AND. do_hfx) THEN
511 6 : ALLOCATE (matrix_hfx(nspin))
512 4 : DO ispin = 1, nspin
513 2 : ALLOCATE (matrix_hfx(nspin)%matrix)
514 2 : CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=smat)
515 2 : CALL dbcsr_copy(matrix_hfx(ispin)%matrix, smat)
516 4 : CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
517 : END DO
518 2 : CALL get_hfx_ks_matrix(qs_env, matrix_hfx)
519 : END IF
520 : ! Loop over spins and atoms
521 62 : DO ispin = 1, nspin
522 32 : CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc, nmo=norb)
523 96 : ALLOCATE (mweight(1, norb))
524 144 : DO iatom = 1, natom
525 : CALL cp_fm_get_submatrix(orb_weight(ispin), mweight, start_row=iatom, &
526 112 : start_col=1, n_rows=1, n_cols=norb)
527 112 : IF (uocc) THEN
528 96 : CALL iao_calculate_dmat(cloc, dkmat, mweight(1, :), .FALSE.)
529 : ELSE
530 16 : CALL iao_calculate_dmat(cloc, dkmat, mweight(1, :), fij_mat(ispin))
531 : END IF
532 112 : CALL dbcsr_dot(dkmat, ks_mat(ispin)%matrix, ecc)
533 112 : ateks(iatom) = ateks(iatom) + ecc
534 : !
535 112 : IF (detailed_ener) THEN
536 18 : CALL dbcsr_dot(dkmat, matrix_t(1)%matrix, e1)
537 18 : atdet(iatom, 1) = atdet(iatom, 1) + e1
538 18 : CALL dbcsr_dot(dkmat, matrix_h(1)%matrix, e2)
539 18 : atdet(iatom, 2) = atdet(iatom, 2) + (e2 - e1)
540 18 : IF (do_hfx) THEN
541 6 : CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 0.5_dp)
542 6 : CALL dbcsr_dot(dkmat, matrix_hfx(ispin)%matrix, ehfx)
543 6 : atdet(iatom, 5) = atdet(iatom, 5) + ehfx
544 6 : CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 2.0_dp)
545 : ELSE
546 12 : ehfx = 0.0_dp
547 : END IF
548 18 : CALL dbcsr_dot(dkmat, exc_mat(1)%matrix, e1)
549 18 : atdet(iatom, 3) = atdet(iatom, 3) + ecc - e2 - e1 - ehfx
550 18 : atdet(iatom, 4) = atdet(iatom, 4) + e1
551 : END IF
552 256 : IF (ewald_correction) THEN
553 6 : CALL dbcsr_dot(dkmat, dve_mat%matrix, ecc)
554 6 : atewald(iatom) = atewald(iatom) + ecc
555 : END IF
556 : !
557 : END DO
558 94 : DEALLOCATE (mweight)
559 : END DO
560 : !
561 30 : IF (detailed_ener) THEN
562 : ! Mulliken
563 6 : CALL get_qs_env(qs_env, rho=rho, para_env=para_env)
564 6 : CALL qs_rho_get(rho, rho_ao=matrix_p)
565 : ! kinetic energy
566 12 : DO ispin = 1, nspin
567 : CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_t(1)%matrix, &
568 6 : para_env, amval)
569 24 : atmul(1:natom, 1) = atmul(1:natom, 1) + amval(1:natom)
570 : CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_h(1)%matrix, &
571 6 : para_env, amval)
572 24 : atmul(1:natom, 2) = atmul(1:natom, 2) + amval(1:natom)
573 24 : atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
574 : CALL mulliken_charges(matrix_p(ispin)%matrix, ks_mat(ispin)%matrix, &
575 6 : para_env, amval)
576 24 : atmul(1:natom, 3) = atmul(1:natom, 3) + amval(1:natom)
577 6 : IF (do_hfx) THEN
578 : CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_hfx(ispin)%matrix, &
579 2 : para_env, amval)
580 8 : atmul(1:natom, 5) = atmul(1:natom, 5) + 0.5_dp*amval(1:natom)
581 8 : atmul(1:natom, 3) = atmul(1:natom, 3) - 0.5_dp*amval(1:natom)
582 : END IF
583 : CALL mulliken_charges(matrix_p(ispin)%matrix, exc_mat(1)%matrix, &
584 6 : para_env, amval)
585 24 : atmul(1:natom, 4) = atmul(1:natom, 4) + amval(1:natom)
586 30 : atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
587 : END DO
588 24 : atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
589 24 : atmul(1:natom, 3) = atmul(1:natom, 3) + 0.5_dp*atmul(1:natom, 2)
590 24 : atmul(1:natom, 2) = 0.5_dp*atmul(1:natom, 2) + 0.5_dp*atcore(1:natom)
591 : END IF
592 : !
593 30 : CALL dbcsr_release(dkmat)
594 62 : DO ispin = 1, nspin
595 32 : CALL dbcsr_release(ks_mat(ispin)%matrix)
596 32 : CALL dbcsr_release(vhxc_mat(ispin)%matrix)
597 32 : DEALLOCATE (ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix)
598 62 : CALL deallocate_mo_set(mos_loc(ispin))
599 : END DO
600 30 : CALL dbcsr_release(core_mat(1)%matrix)
601 30 : DEALLOCATE (core_mat(1)%matrix)
602 30 : CALL dbcsr_release(exc_mat(1)%matrix)
603 30 : DEALLOCATE (exc_mat(1)%matrix)
604 30 : DEALLOCATE (ks_mat, core_mat, vhxc_mat, exc_mat)
605 30 : DEALLOCATE (mos_loc)
606 30 : DEALLOCATE (refbas_blk_sizes)
607 30 : DEALLOCATE (ref_basis_set_list)
608 30 : IF (detailed_ener .AND. do_hfx) THEN
609 4 : DO ispin = 1, nspin
610 2 : CALL dbcsr_release(matrix_hfx(ispin)%matrix)
611 4 : DEALLOCATE (matrix_hfx(ispin)%matrix)
612 : END DO
613 2 : DEALLOCATE (matrix_hfx)
614 : END IF
615 30 : IF (ewald_correction) THEN
616 2 : CALL dbcsr_release(dve_mat%matrix)
617 2 : DEALLOCATE (dve_mat%matrix)
618 : END IF
619 30 : CALL cp_fm_release(orb_weight)
620 30 : CALL cp_fm_release(ciao)
621 30 : CALL cp_fm_release(rotmat)
622 30 : CALL cp_fm_release(c_iao_coef)
623 30 : IF (.NOT. uocc) THEN
624 2 : CALL cp_fm_release(fij_mat)
625 : END IF
626 :
627 : ! KS energy
628 136 : atener(1:natom) = ateks(1:natom)
629 : ! 1/2 of VPP contribution Tr[VPP(K)*P]
630 136 : atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
631 : ! core energy corrections
632 30 : CALL group%sum(atecc)
633 136 : atener(1:natom) = atener(1:natom) + atecc(1:natom)
634 48 : IF (detailed_ener) atdet(1:natom, 6) = atdet(1:natom, 6) + atecc(1:natom)
635 : ! one center terms (GAPW)
636 30 : CALL group%sum(ate1xc)
637 30 : CALL group%sum(ate1h)
638 136 : atener(1:natom) = atener(1:natom) + ate1xc(1:natom) + ate1h(1:natom)
639 30 : IF (detailed_ener) THEN
640 6 : IF (gapw .OR. gapw_xc) atdet(1:natom, 8) = atdet(1:natom, 8) + ate1xc(1:natom)
641 0 : IF (gapw) atdet(1:natom, 9) = atdet(1:natom, 9) + ate1h(1:natom)
642 : END IF
643 : ! core correction
644 136 : atecc(1:natom) = 0.0_dp
645 30 : CALL calculate_ecore_overlap(qs_env, para_env, .FALSE., atecc=atecc)
646 30 : CALL group%sum(atecc)
647 136 : atener(1:natom) = atener(1:natom) + atecc(1:natom)
648 48 : IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
649 30 : IF (ewald_correction) THEN
650 8 : atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
651 : END IF
652 : ! e self
653 136 : atecc(1:natom) = 0.0_dp
654 30 : CALL calculate_ecore_self(qs_env, atecc=atecc)
655 30 : CALL group%sum(atecc)
656 136 : atener(1:natom) = atener(1:natom) + atecc(1:natom)
657 48 : IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
658 30 : IF (ewald_correction) THEN
659 8 : atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
660 2 : CALL calculate_ecore_alpha(qs_env, ealpha, atewald)
661 : END IF
662 : ! vdW pair-potentials
663 136 : atecc(1:natom) = 0.0_dp
664 30 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
665 30 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, evdw, .FALSE., atevdw=atecc)
666 : ! Pair potential gCP energy
667 30 : CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
668 30 : IF (ASSOCIATED(gcp_env)) THEN
669 30 : CALL calculate_gcp_pairpot(qs_env, gcp_env, egcp, .FALSE., ategcp=atecc)
670 : END IF
671 30 : CALL group%sum(atecc)
672 136 : atener(1:natom) = atener(1:natom) + atecc(1:natom)
673 48 : IF (detailed_ener) atdet(1:natom, 10) = atdet(1:natom, 10) + atecc(1:natom)
674 : ! distribute the entropic energy
675 30 : CALL get_qs_env(qs_env, energy=energy)
676 30 : ekts = energy%kts/REAL(natom, KIND=dp)
677 136 : atener(1:natom) = atener(1:natom) + ekts
678 : ! 0.5 Vpp(at)*D + 0.5 * Vpp*D(at)
679 30 : IF (detailed_ener) THEN
680 24 : atdet(1:natom, 3) = atdet(1:natom, 3) + 0.5_dp*atdet(1:natom, 2)
681 24 : atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
682 : END IF
683 : !
684 30 : IF (detailed_ener) THEN
685 6 : IF (unit_nr > 0) THEN
686 3 : WRITE (unit_nr, FMT="(/,T2,A)") "Detailed IAO Atomic Energy Components "
687 3 : CALL write_atdet(unit_nr, atdet(:, 1), "Kinetic Energy")
688 3 : CALL write_atdet(unit_nr, atdet(:, 2), "Short-Range Core and/or PP Energy")
689 3 : CALL write_atdet(unit_nr, atdet(:, 3), "Hartree Energy (long-ranged part)")
690 3 : CALL write_atdet(unit_nr, atdet(:, 5), "Exact Exchange Energy")
691 3 : CALL write_atdet(unit_nr, atdet(:, 4), "Exchange-Correlation Energy")
692 3 : CALL write_atdet(unit_nr, atdet(:, 6), "Atomic Core Hartree: Int(nc V(n+nc) dx")
693 3 : CALL write_atdet(unit_nr, atdet(:, 7), "Core Self Energy Corrections")
694 3 : IF (gapw) THEN
695 0 : CALL write_atdet(unit_nr, atdet(:, 9), "GAPW: 1-center Hartree Energy")
696 : END IF
697 3 : IF (gapw .OR. gapw_xc) THEN
698 0 : CALL write_atdet(unit_nr, atdet(:, 8), "GAPW: 1-center XC Energy")
699 : END IF
700 3 : IF (ABS(evdw) > 1.E-14_dp) THEN
701 0 : CALL write_atdet(unit_nr, atdet(:, 10), "vdW Pairpotential Energy")
702 : END IF
703 12 : DO iatom = 1, natom
704 102 : atecc(iatom) = SUM(atdet(iatom, 1:10)) - atener(iatom)
705 : END DO
706 3 : CALL write_atdet(unit_nr, atecc(:), "Missing Energy")
707 : !
708 3 : WRITE (unit_nr, FMT="(/,T2,A)") "Detailed Mulliken Atomic Energy Components "
709 3 : CALL write_atdet(unit_nr, atmul(:, 1), "Kinetic Energy")
710 3 : CALL write_atdet(unit_nr, atmul(:, 2), "Short-Range Core and/or PP Energy")
711 3 : CALL write_atdet(unit_nr, atmul(:, 3), "Hartree Energy (long-ranged part)")
712 3 : CALL write_atdet(unit_nr, atmul(:, 5), "Exact Exchange Energy")
713 3 : CALL write_atdet(unit_nr, atmul(:, 4), "Exchange-Correlation Energy")
714 3 : CALL write_atdet(unit_nr, atdet(:, 6), "Atomic Core Hartree: Int(nc V(n+nc) dx")
715 3 : CALL write_atdet(unit_nr, atdet(:, 7), "Core Self Energy Corrections")
716 3 : IF (gapw) THEN
717 0 : CALL write_atdet(unit_nr, atdet(:, 9), "GAPW: 1-center Hartree Energy")
718 : END IF
719 3 : IF (gapw .OR. gapw_xc) THEN
720 0 : CALL write_atdet(unit_nr, atdet(:, 8), "GAPW: 1-center XC Energy")
721 : END IF
722 3 : IF (ABS(evdw) > 1.E-14_dp) THEN
723 0 : CALL write_atdet(unit_nr, atdet(:, 10), "vdW Pairpotential Energy")
724 : END IF
725 : END IF
726 : END IF
727 :
728 30 : IF (unit_nr > 0) THEN
729 15 : e_pot = energy%total
730 15 : ateps = 1.E-6_dp
731 15 : CALL write_atener(unit_nr, particle_set, atener, "Atomic Energy Decomposition")
732 68 : sum_energy = SUM(atener(:))
733 15 : checksum = ABS(e_pot - sum_energy)
734 : WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
735 15 : "Potential energy (Atomic):", sum_energy, &
736 15 : "Potential energy (Total) :", e_pot, &
737 30 : "Difference :", checksum
738 15 : CPASSERT((checksum < ateps*ABS(e_pot)))
739 : !
740 15 : IF (ewald_correction) THEN
741 1 : WRITE (UNIT=unit_nr, FMT="(/,(T2,A,F10.3,A))") "Universal Ewald Parameter: ", ealpha, " [Angstrom]"
742 1 : CALL write_atener(unit_nr, particle_set, atewald, "Change in Atomic Energy Decomposition")
743 4 : sum_energy = SUM(atewald(:))
744 : WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
745 1 : "Total Change in Potential energy:", sum_energy
746 : END IF
747 : END IF
748 :
749 30 : IF (detailed_ener) THEN
750 6 : DEALLOCATE (atdet, atmul, amval)
751 : END IF
752 :
753 30 : IF (unit_nr > 0) THEN
754 : WRITE (UNIT=unit_nr, FMT="(/,T2,A)") &
755 15 : "!--------------------------- END OF ED ANALYSIS ------------------------------!"
756 : END IF
757 30 : DEALLOCATE (bcenter)
758 30 : DEALLOCATE (atener, ateks, atecc, atcore, ate1xc, ate1h, atewald)
759 :
760 30 : CALL timestop(handle)
761 :
762 266 : END SUBROUTINE edmf_analysis
763 :
764 : ! **************************************************************************************************
765 : !> \brief ...
766 : !> \param qs_env ...
767 : !> \param vhxc_mat ...
768 : !> \param exc_mat ...
769 : !> \param atecc ...
770 : !> \param ate1xc ...
771 : !> \param ate1h ...
772 : ! **************************************************************************************************
773 30 : SUBROUTINE vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
774 : TYPE(qs_environment_type), POINTER :: qs_env
775 : TYPE(dbcsr_p_type), DIMENSION(:) :: vhxc_mat, exc_mat
776 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: atecc, ate1xc, ate1h
777 :
778 : CHARACTER(len=*), PARAMETER :: routineN = 'vhxc_correction'
779 :
780 : INTEGER :: handle, iatom, ispin, natom, nspins
781 : LOGICAL :: do_admm_corr, gapw, gapw_xc
782 : REAL(KIND=dp) :: eh1, exc1
783 : TYPE(admm_type), POINTER :: admm_env
784 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
785 30 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
786 : TYPE(dft_control_type), POINTER :: dft_control
787 30 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
788 : TYPE(local_rho_type), POINTER :: local_rho_set
789 : TYPE(mp_para_env_type), POINTER :: para_env
790 : TYPE(pw_env_type), POINTER :: pw_env
791 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
792 : TYPE(pw_r3d_rs_type) :: xc_den
793 30 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vtau, vxc
794 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
795 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
796 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
797 : TYPE(qs_ks_env_type), POINTER :: ks_env
798 : TYPE(qs_rho_type), POINTER :: rho_struct
799 30 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
800 : TYPE(section_vals_type), POINTER :: xc_fun_section, xc_section
801 : TYPE(xc_rho_cflags_type) :: needs
802 :
803 30 : CALL timeset(routineN, handle)
804 :
805 30 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
806 :
807 30 : nspins = dft_control%nspins
808 30 : gapw = dft_control%qs_control%gapw
809 30 : gapw_xc = dft_control%qs_control%gapw_xc
810 30 : do_admm_corr = .FALSE.
811 30 : IF (dft_control%do_admm) THEN
812 0 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) do_admm_corr = .TRUE.
813 : END IF
814 : IF (do_admm_corr) THEN
815 0 : CALL get_qs_env(qs_env, admm_env=admm_env)
816 0 : xc_section => admm_env%xc_section_aux
817 : ELSE
818 30 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
819 : END IF
820 30 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
821 30 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
822 :
823 30 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
824 30 : CALL auxbas_pw_pool%create_pw(xc_den)
825 122 : ALLOCATE (vxc(nspins))
826 62 : DO ispin = 1, nspins
827 62 : CALL auxbas_pw_pool%create_pw(vxc(ispin))
828 : END DO
829 30 : IF (needs%tau .OR. needs%tau_spin) THEN
830 24 : ALLOCATE (vtau(nspins))
831 16 : DO ispin = 1, nspins
832 38 : CALL auxbas_pw_pool%create_pw(vtau(ispin))
833 : END DO
834 : END IF
835 :
836 : ! Nuclear charge correction
837 30 : CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
838 30 : CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
839 30 : IF (gapw .OR. gapw_xc) THEN
840 : CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, &
841 : rho_atom_set=rho_atom_set, ecoul_1c=ecoul_1c, &
842 8 : natom=natom, para_env=para_env)
843 8 : CALL zero_rho_atom_integrals(rho_atom_set)
844 8 : CALL calculate_vxc_atom(qs_env, .FALSE., exc1)
845 8 : IF (gapw) THEN
846 6 : CALL Vh_1c_gg_integrals(qs_env, eh1, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
847 6 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
848 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
849 6 : local_rho_set=local_rho_set, atener=atecc)
850 : END IF
851 : END IF
852 :
853 8 : IF (gapw_xc) THEN
854 2 : CALL get_qs_env(qs_env, rho_xc=rho_struct, dispersion_env=dispersion_env)
855 : ELSE
856 28 : CALL get_qs_env(qs_env, rho=rho_struct, dispersion_env=dispersion_env)
857 : END IF
858 : !
859 30 : CPASSERT(.NOT. do_admm_corr)
860 : !
861 30 : IF (needs%tau .OR. needs%tau_spin) THEN
862 : CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
863 8 : xc_den=xc_den, vxc=vxc, vtau=vtau)
864 : ELSE
865 : CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
866 22 : xc_den=xc_den, vxc=vxc)
867 : END IF
868 62 : DO ispin = 1, nspins
869 32 : CALL pw_scale(vxc(ispin), -0.5_dp)
870 32 : CALL pw_axpy(xc_den, vxc(ispin))
871 32 : CALL pw_scale(vxc(ispin), vxc(ispin)%pw_grid%dvol)
872 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vxc(ispin), &
873 : hmat=vhxc_mat(ispin), calculate_forces=.FALSE., &
874 32 : gapw=(gapw .OR. gapw_xc))
875 62 : IF (needs%tau .OR. needs%tau_spin) THEN
876 8 : CALL pw_scale(vtau(ispin), -0.5_dp*vtau(ispin)%pw_grid%dvol)
877 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vtau(ispin), &
878 : hmat=vhxc_mat(ispin), calculate_forces=.FALSE., &
879 8 : compute_tau=.TRUE., gapw=(gapw .OR. gapw_xc))
880 : END IF
881 : END DO
882 30 : CALL pw_scale(xc_den, xc_den%pw_grid%dvol)
883 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=xc_den, &
884 : hmat=exc_mat(1), calculate_forces=.FALSE., &
885 30 : gapw=(gapw .OR. gapw_xc))
886 :
887 30 : IF (gapw .OR. gapw_xc) THEN
888 : ! remove one-center potential matrix part
889 8 : CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
890 8 : CALL update_ks_atom(qs_env, vhxc_mat, matrix_p, forces=.FALSE., kscale=-0.5_dp)
891 : !
892 32 : DO iatom = 1, natom
893 : ate1xc(iatom) = ate1xc(iatom) + &
894 32 : rho_atom_set(iatom)%exc_h - rho_atom_set(iatom)%exc_s
895 : END DO
896 8 : IF (gapw) THEN
897 24 : DO iatom = 1, natom
898 : ate1h(iatom) = ate1h(iatom) + &
899 : ecoul_1c(iatom)%ecoul_1_h - ecoul_1c(iatom)%ecoul_1_s + &
900 24 : ecoul_1c(iatom)%ecoul_1_z - ecoul_1c(iatom)%ecoul_1_0
901 : END DO
902 : END IF
903 : END IF
904 :
905 30 : CALL auxbas_pw_pool%give_back_pw(xc_den)
906 62 : DO ispin = 1, nspins
907 62 : CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
908 : END DO
909 30 : IF (needs%tau .OR. needs%tau_spin) THEN
910 16 : DO ispin = 1, nspins
911 38 : CALL auxbas_pw_pool%give_back_pw(vtau(ispin))
912 : END DO
913 : END IF
914 :
915 30 : CALL timestop(handle)
916 :
917 60 : END SUBROUTINE vhxc_correction
918 :
919 : ! **************************************************************************************************
920 : !> \brief ...
921 : !> \param qs_env ...
922 : !> \param ealpha ...
923 : !> \param vh_mat ...
924 : !> \param atewd ...
925 : ! **************************************************************************************************
926 2 : SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd)
927 : TYPE(qs_environment_type), POINTER :: qs_env
928 : REAL(KIND=dp), INTENT(IN) :: ealpha
929 : TYPE(dbcsr_p_type) :: vh_mat
930 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: atewd
931 :
932 : CHARACTER(len=*), PARAMETER :: routineN = 'vh_ewald_correction'
933 :
934 : INTEGER :: handle, ikind, ispin, natom, nkind
935 : REAL(KIND=dp) :: eps_core_charge, rhotot, zeff
936 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, atcore, atecc, ccore
937 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
938 : TYPE(dbcsr_p_type) :: e_mat
939 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
940 : TYPE(dft_control_type), POINTER :: dft_control
941 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
942 2 : POINTER :: sab_orb, sac_ae
943 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
944 : TYPE(pw_c1d_gs_type) :: rho_tot_ewd_g, v_hewd_gspace
945 : TYPE(pw_env_type), POINTER :: pw_env
946 : TYPE(pw_poisson_type), POINTER :: poisson_env
947 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
948 : TYPE(pw_r3d_rs_type) :: rho_tot_ewd_r, v_hewd_rspace
949 2 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
950 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
951 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
952 : TYPE(qs_rho_type), POINTER :: rho
953 :
954 2 : CALL timeset(routineN, handle)
955 :
956 2 : natom = SIZE(atewd)
957 8 : ALLOCATE (atecc(natom), atcore(natom))
958 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, qs_kind_set=qs_kind_set, &
959 2 : dft_control=dft_control)
960 8 : ALLOCATE (alpha(nkind), ccore(nkind))
961 6 : DO ikind = 1, nkind
962 4 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
963 4 : alpha(ikind) = ealpha
964 6 : ccore(ikind) = zeff*(ealpha/pi)**1.5_dp
965 : END DO
966 : !
967 2 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
968 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
969 2 : poisson_env=poisson_env)
970 : !
971 2 : CALL auxbas_pw_pool%create_pw(v_hewd_gspace)
972 2 : CALL auxbas_pw_pool%create_pw(v_hewd_rspace)
973 2 : CALL auxbas_pw_pool%create_pw(rho_tot_ewd_g)
974 2 : CALL auxbas_pw_pool%create_pw(rho_tot_ewd_r)
975 : rhotot = 0.0_dp
976 2 : CALL calculate_rho_core(rho_tot_ewd_r, rhotot, qs_env, calpha=alpha, ccore=ccore)
977 : ! Get the total density in g-space [ions + electrons]
978 2 : CALL get_qs_env(qs_env=qs_env, rho=rho)
979 2 : CALL qs_rho_get(rho, rho_r=rho_r)
980 4 : DO ispin = 1, dft_control%nspins
981 4 : CALL pw_axpy(rho_r(ispin), rho_tot_ewd_r)
982 : END DO
983 2 : CALL pw_transfer(rho_tot_ewd_r, rho_tot_ewd_g)
984 : !
985 2 : CALL pw_poisson_solve(poisson_env, density=rho_tot_ewd_g, vhartree=v_hewd_gspace)
986 2 : CALL pw_transfer(v_hewd_gspace, v_hewd_rspace)
987 2 : CALL pw_scale(v_hewd_rspace, v_hewd_rspace%pw_grid%dvol)
988 8 : atecc = 0.0_dp
989 2 : CALL integrate_v_gaussian_rspace(v_hewd_rspace, qs_env, alpha, ccore, atecc)
990 8 : atewd(1:natom) = atewd(1:natom) + atecc(1:natom)
991 : !
992 8 : atecc = 0.0_dp
993 2 : CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
994 2 : CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
995 8 : atewd(1:natom) = atewd(1:natom) - atecc(1:natom)
996 : !
997 2 : CALL pw_axpy(v_hartree_rspace, v_hewd_rspace, -1.0_dp)
998 2 : CALL dbcsr_set(vh_mat%matrix, 0.0_dp)
999 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hewd_rspace, hmat=vh_mat, &
1000 2 : calculate_forces=.FALSE.)
1001 : !
1002 2 : CALL dbcsr_scale(vh_mat%matrix, 0.5_dp)
1003 : !
1004 : ! delta erfc
1005 2 : CALL qs_rho_get(rho, rho_ao=matrix_p)
1006 2 : eps_core_charge = dft_control%qs_control%eps_core_charge
1007 2 : ALLOCATE (e_mat%matrix)
1008 2 : CALL dbcsr_create(e_mat%matrix, template=vh_mat%matrix)
1009 2 : CALL dbcsr_copy(e_mat%matrix, vh_mat%matrix)
1010 : !
1011 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
1012 2 : particle_set=particle_set, sab_orb=sab_orb, sac_ppl=sac_ae)
1013 2 : CALL dbcsr_set(e_mat%matrix, 0.0_dp)
1014 8 : atcore = 0.0_dp
1015 : CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
1016 2 : alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
1017 8 : atewd(1:natom) = atewd(1:natom) + 0.5_dp*atcore(1:natom)
1018 2 : CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, 0.5_dp)
1019 6 : DO ikind = 1, nkind
1020 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha(ikind), &
1021 6 : ccore_charge=ccore(ikind))
1022 : END DO
1023 2 : CALL dbcsr_set(e_mat%matrix, 0.0_dp)
1024 8 : atcore = 0.0_dp
1025 : CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
1026 2 : alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
1027 8 : atewd(1:natom) = atewd(1:natom) - 0.5_dp*atcore(1:natom)
1028 2 : CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, -0.5_dp)
1029 : !
1030 2 : CALL dbcsr_release(e_mat%matrix)
1031 2 : DEALLOCATE (e_mat%matrix)
1032 2 : CALL auxbas_pw_pool%give_back_pw(v_hewd_gspace)
1033 2 : CALL auxbas_pw_pool%give_back_pw(v_hewd_rspace)
1034 2 : CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_g)
1035 2 : CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_r)
1036 2 : DEALLOCATE (atecc, atcore)
1037 2 : DEALLOCATE (alpha, ccore)
1038 :
1039 2 : CALL timestop(handle)
1040 :
1041 4 : END SUBROUTINE vh_ewald_correction
1042 :
1043 : ! **************************************************************************************************
1044 : !> \brief ...
1045 : !> \param qs_env ...
1046 : !> \param matrix_hfx ...
1047 : ! **************************************************************************************************
1048 2 : SUBROUTINE get_hfx_ks_matrix(qs_env, matrix_hfx)
1049 : TYPE(qs_environment_type), POINTER :: qs_env
1050 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hfx
1051 :
1052 : CHARACTER(len=*), PARAMETER :: routineN = 'get_hfx_ks_matrix'
1053 :
1054 : INTEGER :: handle
1055 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
1056 : TYPE(qs_rho_type), POINTER :: rho
1057 :
1058 2 : CALL timeset(routineN, handle)
1059 :
1060 2 : CALL get_qs_env(qs_env, rho=rho)
1061 2 : CALL qs_rho_get(rho, rho_ao=matrix_p)
1062 : CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, update_energy=.FALSE., &
1063 2 : recalc_integrals=.FALSE.)
1064 :
1065 2 : CALL timestop(handle)
1066 :
1067 2 : END SUBROUTINE get_hfx_ks_matrix
1068 : ! **************************************************************************************************
1069 : !> \brief Write the atomic coordinates, types, and energies
1070 : !> \param iounit ...
1071 : !> \param particle_set ...
1072 : !> \param atener ...
1073 : !> \param label ...
1074 : !> \date 05.06.2023
1075 : !> \author JGH
1076 : !> \version 1.0
1077 : ! **************************************************************************************************
1078 16 : SUBROUTINE write_atener(iounit, particle_set, atener, label)
1079 :
1080 : INTEGER, INTENT(IN) :: iounit
1081 : TYPE(particle_type), DIMENSION(:) :: particle_set
1082 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: atener
1083 : CHARACTER(LEN=*), INTENT(IN) :: label
1084 :
1085 : INTEGER :: i, natom
1086 :
1087 16 : IF (iounit > 0) THEN
1088 16 : WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
1089 : WRITE (UNIT=iounit, FMT="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
1090 16 : "Atom Element", "X", "Y", "Z", "Energy[a.u.]"
1091 16 : natom = SIZE(atener)
1092 72 : DO i = 1, natom
1093 56 : WRITE (UNIT=iounit, FMT="(I6,T12,A2,T24,3F12.6,F21.10)") i, &
1094 56 : TRIM(ADJUSTL(particle_set(i)%atomic_kind%element_symbol)), &
1095 296 : particle_set(i)%r(1:3)*angstrom, atener(i)
1096 : END DO
1097 16 : WRITE (UNIT=iounit, FMT="(A)") ""
1098 : END IF
1099 :
1100 16 : END SUBROUTINE write_atener
1101 :
1102 : ! **************************************************************************************************
1103 : !> \brief Write the one component of atomic energies
1104 : !> \param iounit ...
1105 : !> \param atener ...
1106 : !> \param label ...
1107 : !> \date 18.08.2023
1108 : !> \author JGH
1109 : !> \version 1.0
1110 : ! **************************************************************************************************
1111 45 : SUBROUTINE write_atdet(iounit, atener, label)
1112 : INTEGER, INTENT(IN) :: iounit
1113 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: atener
1114 : CHARACTER(LEN=*), INTENT(IN) :: label
1115 :
1116 : INTEGER :: natom
1117 :
1118 45 : IF (iounit > 0) THEN
1119 45 : natom = SIZE(atener)
1120 45 : WRITE (UNIT=iounit, FMT="(T2,A)") TRIM(label)
1121 45 : WRITE (UNIT=iounit, FMT="(5G16.8)") atener(1:natom)
1122 : END IF
1123 :
1124 45 : END SUBROUTINE write_atdet
1125 :
1126 : END MODULE ed_analysis
|
