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