Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate MAO's and analyze wavefunctions
10 : !> \par History
11 : !> 03.2016 created [JGH]
12 : !> 12.2016 split into four modules [JGH]
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE mao_wfn_analysis
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE basis_set_types, ONLY: gto_basis_set_p_type
18 : USE bibliography, ONLY: Ehrhardt1985,&
19 : Heinzmann1976,&
20 : cite_reference
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_get_block_p, &
25 : dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
26 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
27 : dbcsr_p_type, dbcsr_release, dbcsr_replicate_all, dbcsr_type, dbcsr_type_no_symmetry, &
28 : dbcsr_type_symmetric
29 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
30 : cp_dbcsr_cholesky_restore
31 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
32 : dbcsr_get_block_diag,&
33 : dbcsr_reserve_diag_blocks
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
36 : dbcsr_deallocate_matrix_set
37 : USE input_section_types, ONLY: section_vals_get,&
38 : section_vals_type,&
39 : section_vals_val_get
40 : USE iterate_matrix, ONLY: invert_Hotelling
41 : USE kinds, ONLY: dp
42 : USE kpoint_types, ONLY: kpoint_type
43 : USE mao_io, ONLY: mao_write_pao_restart
44 : USE mao_methods, ONLY: mao_basis_analysis,&
45 : mao_build_q,&
46 : mao_reference_basis
47 : USE mao_optimizer, ONLY: mao_optimize
48 : USE mathlib, ONLY: invmat_symm
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE particle_methods, ONLY: get_particle_set
51 : USE particle_types, ONLY: particle_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_kind_types, ONLY: get_qs_kind,&
55 : qs_kind_type
56 : USE qs_ks_types, ONLY: get_ks_env,&
57 : qs_ks_env_type
58 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
59 : neighbor_list_iterate,&
60 : neighbor_list_iterator_create,&
61 : neighbor_list_iterator_p_type,&
62 : neighbor_list_iterator_release,&
63 : neighbor_list_set_p_type,&
64 : release_neighbor_list_sets
65 : USE qs_neighbor_lists, ONLY: setup_neighbor_list
66 : USE qs_overlap, ONLY: build_overlap_matrix_simple
67 : USE qs_rho_types, ONLY: qs_rho_get,&
68 : qs_rho_type
69 : #include "./base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 : PRIVATE
73 :
74 : TYPE block_type
75 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: mat
76 : END TYPE block_type
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_wfn_analysis'
79 :
80 : PUBLIC :: mao_analysis
81 :
82 : ! **************************************************************************************************
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief ...
88 : !> \param qs_env ...
89 : !> \param input_section ...
90 : !> \param unit_nr ...
91 : ! **************************************************************************************************
92 38 : SUBROUTINE mao_analysis(qs_env, input_section, unit_nr)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : TYPE(section_vals_type), POINTER :: input_section
95 : INTEGER, INTENT(IN) :: unit_nr
96 :
97 : CHARACTER(len=*), PARAMETER :: routineN = 'mao_analysis'
98 :
99 : CHARACTER(len=2) :: element_symbol, esa, esb, esc
100 : INTEGER :: fall, handle, ia, iab, iabc, iatom, ib, ic, icol, ikind, irow, ispin, jatom, &
101 : mao_basis, max_iter, me, na, nab, nabc, natom, nb, nc, nimages, nspin, ssize
102 38 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, mao_blk_sizes, &
103 38 : orb_blk, row_blk_sizes
104 : LOGICAL :: analyze_ua, explicit, fo, for, fos, &
105 : found, neglect_abc, print_basis, &
106 : print_pao
107 : REAL(KIND=dp) :: deltaq, electra(2), eps_ab, eps_abc, eps_filter, eps_fun, eps_grad, epsx, &
108 : senabc, senmax, threshold, total_charge, total_spin, ua_charge(2), zeff
109 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: occnumA, occnumABC, qab, qmatab, qmatac, &
110 38 : qmatbc, raq, sab, selnABC, sinv, &
111 38 : smatab, smatac, smatbc, uaq
112 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: occnumAB, selnAB
113 38 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, cmao, diag, qblka, qblkb, qblkc, &
114 38 : rblkl, rblku, sblk, sblka, sblkb, sblkc
115 38 : TYPE(block_type), ALLOCATABLE, DIMENSION(:) :: rowblock
116 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
117 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
118 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
119 38 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, mao_dmat, mao_qmat, mao_smat, &
120 38 : matrix_q, matrix_smm, matrix_smo
121 38 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
122 : TYPE(dbcsr_type) :: amat, axmat, cgmat, cholmat, crumat, &
123 : qmat, qmat_diag, rumat, smat_diag, &
124 : sumat, tmat
125 : TYPE(dft_control_type), POINTER :: dft_control
126 38 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
127 : TYPE(kpoint_type), POINTER :: kpoints
128 : TYPE(mp_para_env_type), POINTER :: para_env
129 : TYPE(neighbor_list_iterator_p_type), &
130 38 : DIMENSION(:), POINTER :: nl_iterator
131 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
132 38 : POINTER :: sab_all, sab_orb, smm_list, smo_list
133 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
134 38 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135 : TYPE(qs_ks_env_type), POINTER :: ks_env
136 : TYPE(qs_rho_type), POINTER :: rho
137 :
138 : ! only do MAO analysis if explicitely requested
139 :
140 38 : CALL section_vals_get(input_section, explicit=explicit)
141 38 : IF (.NOT. explicit) RETURN
142 :
143 10 : CALL timeset(routineN, handle)
144 :
145 10 : IF (unit_nr > 0) THEN
146 5 : WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
147 5 : WRITE (UNIT=unit_nr, FMT="(T36,A)") "MAO ANALYSIS"
148 5 : WRITE (UNIT=unit_nr, FMT="(T12,A)") "Claus Ehrhardt and Reinhart Ahlrichs, TCA 68:231-245 (1985)"
149 5 : WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
150 : END IF
151 10 : CALL cite_reference(Heinzmann1976)
152 10 : CALL cite_reference(Ehrhardt1985)
153 :
154 : ! input options
155 10 : CALL section_vals_val_get(input_section, "REFERENCE_BASIS", i_val=mao_basis)
156 10 : CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter)
157 10 : CALL section_vals_val_get(input_section, "EPS_FUNCTION", r_val=eps_fun)
158 10 : CALL section_vals_val_get(input_section, "EPS_GRAD", r_val=eps_grad)
159 10 : CALL section_vals_val_get(input_section, "MAX_ITER", i_val=max_iter)
160 10 : CALL section_vals_val_get(input_section, "PRINT_BASIS", l_val=print_basis)
161 10 : CALL section_vals_val_get(input_section, "PRINT_PAO", l_val=print_pao)
162 10 : CALL section_vals_val_get(input_section, "NEGLECT_ABC", l_val=neglect_abc)
163 10 : CALL section_vals_val_get(input_section, "AB_THRESHOLD", r_val=eps_ab)
164 10 : CALL section_vals_val_get(input_section, "ABC_THRESHOLD", r_val=eps_abc)
165 10 : CALL section_vals_val_get(input_section, "ANALYZE_UNASSIGNED_CHARGE", l_val=analyze_ua)
166 :
167 : ! k-points?
168 10 : CALL get_qs_env(qs_env, dft_control=dft_control)
169 10 : nimages = dft_control%nimages
170 10 : IF (nimages > 1) THEN
171 0 : IF (unit_nr > 0) THEN
172 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
173 0 : "K-Points: MAO's determined and analyzed using Gamma-Point only."
174 : END IF
175 : END IF
176 :
177 : ! Reference basis set
178 10 : NULLIFY (mao_basis_set_list, orb_basis_set_list)
179 : CALL mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, &
180 10 : unit_nr, print_basis)
181 :
182 : ! neighbor lists
183 10 : NULLIFY (smm_list, smo_list)
184 10 : CALL setup_neighbor_list(smm_list, mao_basis_set_list, qs_env=qs_env)
185 10 : CALL setup_neighbor_list(smo_list, mao_basis_set_list, orb_basis_set_list, qs_env=qs_env)
186 :
187 : ! overlap matrices
188 10 : NULLIFY (matrix_smm, matrix_smo)
189 10 : CALL get_qs_env(qs_env, ks_env=ks_env)
190 : CALL build_overlap_matrix_simple(ks_env, matrix_smm, &
191 10 : mao_basis_set_list, mao_basis_set_list, smm_list)
192 : CALL build_overlap_matrix_simple(ks_env, matrix_smo, &
193 10 : mao_basis_set_list, orb_basis_set_list, smo_list)
194 :
195 : ! get reference density matrix and overlap matrix
196 10 : CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s)
197 10 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
198 10 : nspin = SIZE(matrix_p, 1)
199 : !
200 : ! Q matrix
201 10 : IF (nimages == 1) THEN
202 10 : CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter)
203 : ELSE
204 0 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints)
205 : CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, &
206 0 : nimages=nimages, kpoints=kpoints, matrix_ks=matrix_ks, sab_orb=sab_orb)
207 : END IF
208 :
209 : ! check for extended basis sets
210 10 : fall = 0
211 10 : CALL neighbor_list_iterator_create(nl_iterator, smm_list)
212 97 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
213 87 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
214 87 : IF (iatom <= jatom) THEN
215 53 : irow = iatom
216 53 : icol = jatom
217 : ELSE
218 34 : irow = jatom
219 34 : icol = iatom
220 : END IF
221 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
222 87 : row=irow, col=icol, block=block, found=found)
223 97 : IF (.NOT. found) fall = fall + 1
224 : END DO
225 10 : CALL neighbor_list_iterator_release(nl_iterator)
226 :
227 10 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
228 10 : CALL para_env%sum(fall)
229 10 : IF (unit_nr > 0 .AND. fall > 0) THEN
230 : WRITE (UNIT=unit_nr, FMT="(/,T2,A,/,T2,A,/)") &
231 0 : "Warning: Extended MAO basis used with original basis filtered density matrix", &
232 0 : "Warning: Possible errors can be controlled with EPS_PGF_ORB"
233 : END IF
234 :
235 : ! MAO matrices
236 10 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
237 10 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
238 10 : NULLIFY (mao_coef)
239 10 : CALL dbcsr_allocate_matrix_set(mao_coef, nspin)
240 40 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
241 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
242 10 : basis=mao_basis_set_list)
243 10 : CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
244 : ! check if MAOs have been specified
245 58 : DO iab = 1, natom
246 48 : IF (col_blk_sizes(iab) < 0) &
247 10 : CPABORT("Number of MAOs has to be specified in KIND section for all elements")
248 : END DO
249 22 : DO ispin = 1, nspin
250 : ! coeficients
251 12 : ALLOCATE (mao_coef(ispin)%matrix)
252 : CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, &
253 : name="MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
254 12 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
255 22 : CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
256 : END DO
257 10 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
258 :
259 : ! optimize MAOs
260 10 : epsx = 1000.0_dp
261 : CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, epsx, &
262 10 : 3, unit_nr)
263 :
264 : ! Analyze the MAO basis
265 : CALL mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
266 10 : qs_kind_set, unit_nr, para_env)
267 :
268 : ! Calculate the overlap and density matrix in the new MAO basis
269 10 : NULLIFY (mao_dmat, mao_smat, mao_qmat)
270 10 : CALL dbcsr_allocate_matrix_set(mao_qmat, nspin)
271 10 : CALL dbcsr_allocate_matrix_set(mao_dmat, nspin)
272 10 : CALL dbcsr_allocate_matrix_set(mao_smat, nspin)
273 10 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
274 22 : DO ispin = 1, nspin
275 12 : ALLOCATE (mao_dmat(ispin)%matrix)
276 : CALL dbcsr_create(mao_dmat(ispin)%matrix, name="MAO density", dist=dbcsr_dist, &
277 : matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
278 12 : col_blk_size=col_blk_sizes)
279 12 : ALLOCATE (mao_smat(ispin)%matrix)
280 : CALL dbcsr_create(mao_smat(ispin)%matrix, name="MAO overlap", dist=dbcsr_dist, &
281 : matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
282 12 : col_blk_size=col_blk_sizes)
283 12 : ALLOCATE (mao_qmat(ispin)%matrix)
284 : CALL dbcsr_create(mao_qmat(ispin)%matrix, name="MAO covar density", dist=dbcsr_dist, &
285 : matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
286 22 : col_blk_size=col_blk_sizes)
287 : END DO
288 10 : CALL dbcsr_create(amat, name="MAO overlap", template=mao_dmat(1)%matrix)
289 10 : CALL dbcsr_create(tmat, name="MAO Overlap Inverse", template=amat)
290 10 : CALL dbcsr_create(qmat, name="MAO covar density", template=amat)
291 10 : CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
292 10 : CALL dbcsr_create(axmat, name="TEMP", template=amat, matrix_type=dbcsr_type_no_symmetry)
293 22 : DO ispin = 1, nspin
294 : ! calculate MAO overlap matrix
295 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smm(1)%matrix, mao_coef(ispin)%matrix, &
296 12 : 0.0_dp, cgmat)
297 12 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, amat)
298 : ! calculate inverse of MAO overlap
299 12 : threshold = 1.e-8_dp
300 12 : CALL invert_Hotelling(tmat, amat, threshold, norm_convergence=1.e-4_dp, silent=.TRUE.)
301 12 : CALL dbcsr_copy(mao_smat(ispin)%matrix, amat)
302 : ! calculate q-matrix q = C*Q*C
303 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_q(ispin)%matrix, mao_coef(ispin)%matrix, &
304 12 : 0.0_dp, cgmat, filter_eps=eps_filter)
305 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, &
306 12 : 0.0_dp, qmat, filter_eps=eps_filter)
307 12 : CALL dbcsr_copy(mao_qmat(ispin)%matrix, qmat)
308 : ! calculate density matrix
309 12 : CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, tmat, 0.0_dp, axmat, filter_eps=eps_filter)
310 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, axmat, 0.0_dp, mao_dmat(ispin)%matrix, &
311 22 : filter_eps=eps_filter)
312 : END DO
313 10 : CALL dbcsr_release(amat)
314 10 : CALL dbcsr_release(tmat)
315 10 : CALL dbcsr_release(qmat)
316 10 : CALL dbcsr_release(cgmat)
317 10 : CALL dbcsr_release(axmat)
318 :
319 : ! calculate unassigned charge : n - Tr PS
320 22 : DO ispin = 1, nspin
321 12 : CALL dbcsr_dot(mao_dmat(ispin)%matrix, mao_smat(ispin)%matrix, ua_charge(ispin))
322 22 : ua_charge(ispin) = electra(ispin) - ua_charge(ispin)
323 : END DO
324 10 : IF (unit_nr > 0) THEN
325 5 : WRITE (unit_nr, *)
326 11 : DO ispin = 1, nspin
327 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
328 11 : "Unassigned charge", "Spin ", ispin, "delta charge =", ua_charge(ispin)
329 : END DO
330 : END IF
331 :
332 : ! occupation numbers: single atoms
333 : ! We use S_A = 1
334 : ! At the gamma point we use an effective MIC
335 10 : CALL get_qs_env(qs_env, natom=natom)
336 40 : ALLOCATE (occnumA(natom, nspin))
337 10 : occnumA = 0.0_dp
338 22 : DO ispin = 1, nspin
339 76 : DO iatom = 1, natom
340 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
341 54 : row=iatom, col=iatom, block=block, found=found)
342 66 : IF (found) THEN
343 81 : DO iab = 1, SIZE(block, 1)
344 81 : occnumA(iatom, ispin) = occnumA(iatom, ispin) + block(iab, iab)
345 : END DO
346 : END IF
347 : END DO
348 : END DO
349 10 : CALL para_env%sum(occnumA)
350 :
351 : ! occupation numbers: atom pairs
352 50 : ALLOCATE (occnumAB(natom, natom, nspin))
353 10 : occnumAB = 0.0_dp
354 22 : DO ispin = 1, nspin
355 12 : CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
356 12 : CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix)
357 : ! replicate the diagonal blocks of the density and overlap matrices
358 12 : CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
359 12 : CALL dbcsr_replicate_all(qmat_diag)
360 12 : CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
361 12 : CALL dbcsr_replicate_all(smat_diag)
362 66 : DO ia = 1, natom
363 174 : DO ib = ia + 1, natom
364 108 : iab = 0
365 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
366 108 : row=ia, col=ib, block=block, found=found)
367 108 : IF (found) iab = 1
368 108 : CALL para_env%sum(iab)
369 108 : CPASSERT(iab <= 1)
370 270 : IF (iab == 0 .AND. para_env%is_source()) THEN
371 : ! AB block is not available N_AB = N_A + N_B
372 : ! Do this only on the "source" processor
373 0 : occnumAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin)
374 0 : occnumAB(ib, ia, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin)
375 108 : ELSE IF (found) THEN
376 : ! owner of AB block performs calculation
377 54 : na = SIZE(block, 1)
378 54 : nb = SIZE(block, 2)
379 54 : nab = na + nb
380 432 : ALLOCATE (sab(nab, nab), qab(nab, nab), sinv(nab, nab))
381 : ! qmat
382 327 : qab(1:na, na + 1:nab) = block(1:na, 1:nb)
383 375 : qab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb))
384 54 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=diag, found=fo)
385 54 : CPASSERT(fo)
386 630 : qab(1:na, 1:na) = diag(1:na, 1:na)
387 54 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=diag, found=fo)
388 54 : CPASSERT(fo)
389 342 : qab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
390 : ! smat
391 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, &
392 54 : row=ia, col=ib, block=block, found=fo)
393 54 : CPASSERT(fo)
394 327 : sab(1:na, na + 1:nab) = block(1:na, 1:nb)
395 375 : sab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb))
396 54 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=diag, found=fo)
397 54 : CPASSERT(fo)
398 630 : sab(1:na, 1:na) = diag(1:na, 1:na)
399 54 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=diag, found=fo)
400 54 : CPASSERT(fo)
401 342 : sab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
402 : ! inv smat
403 1296 : sinv(1:nab, 1:nab) = sab(1:nab, 1:nab)
404 54 : CALL invmat_symm(sinv)
405 : ! Tr(Q*Sinv)
406 1296 : occnumAB(ia, ib, ispin) = SUM(qab*sinv)
407 54 : occnumAB(ib, ia, ispin) = occnumAB(ia, ib, ispin)
408 : !
409 324 : DEALLOCATE (sab, qab, sinv)
410 : END IF
411 : END DO
412 : END DO
413 12 : CALL dbcsr_release(qmat_diag)
414 22 : CALL dbcsr_release(smat_diag)
415 : END DO
416 10 : CALL para_env%sum(occnumAB)
417 :
418 : ! calculate shared electron numbers (AB)
419 50 : ALLOCATE (selnAB(natom, natom, nspin))
420 10 : selnAB = 0.0_dp
421 22 : DO ispin = 1, nspin
422 76 : DO ia = 1, natom
423 174 : DO ib = ia + 1, natom
424 108 : selnAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) - occnumAB(ia, ib, ispin)
425 162 : selnAB(ib, ia, ispin) = selnAB(ia, ib, ispin)
426 : END DO
427 : END DO
428 : END DO
429 :
430 10 : IF (.NOT. neglect_abc) THEN
431 : ! calculate N_ABC
432 8 : nabc = (natom*(natom - 1)*(natom - 2))/6
433 32 : ALLOCATE (occnumABC(nabc, nspin))
434 142 : occnumABC = -1.0_dp
435 18 : DO ispin = 1, nspin
436 10 : CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
437 10 : CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix)
438 : ! replicate the diagonal blocks of the density and overlap matrices
439 10 : CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
440 10 : CALL dbcsr_replicate_all(qmat_diag)
441 10 : CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
442 10 : CALL dbcsr_replicate_all(smat_diag)
443 10 : iabc = 0
444 58 : DO ia = 1, natom
445 48 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=qblka, found=fo)
446 48 : CPASSERT(fo)
447 48 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=sblka, found=fo)
448 48 : CPASSERT(fo)
449 48 : na = SIZE(qblka, 1)
450 256 : DO ib = ia + 1, natom
451 : ! screen with SEN(AB)
452 102 : IF (selnAB(ia, ib, ispin) < eps_abc) THEN
453 14 : iabc = iabc + (natom - ib)
454 14 : CYCLE
455 : END IF
456 88 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=qblkb, found=fo)
457 88 : CPASSERT(fo)
458 88 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=sblkb, found=fo)
459 88 : CPASSERT(fo)
460 88 : nb = SIZE(qblkb, 1)
461 88 : nab = na + nb
462 528 : ALLOCATE (qmatab(na, nb), smatab(na, nb))
463 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ib, &
464 88 : block=block, found=found)
465 88 : qmatab = 0.0_dp
466 320 : IF (found) qmatab(1:na, 1:nb) = block(1:na, 1:nb)
467 88 : CALL para_env%sum(qmatab)
468 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ib, &
469 88 : block=block, found=found)
470 88 : smatab = 0.0_dp
471 320 : IF (found) smatab(1:na, 1:nb) = block(1:na, 1:nb)
472 88 : CALL para_env%sum(smatab)
473 202 : DO ic = ib + 1, natom
474 : ! screen with SEN(AB)
475 114 : IF ((selnAB(ia, ic, ispin) < eps_abc) .OR. (selnAB(ib, ic, ispin) < eps_abc)) THEN
476 24 : iabc = iabc + 1
477 24 : CYCLE
478 : END IF
479 90 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ic, col=ic, block=qblkc, found=fo)
480 90 : CPASSERT(fo)
481 90 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ic, col=ic, block=sblkc, found=fo)
482 90 : CPASSERT(fo)
483 90 : nc = SIZE(qblkc, 1)
484 540 : ALLOCATE (qmatac(na, nc), smatac(na, nc))
485 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ic, &
486 90 : block=block, found=found)
487 90 : qmatac = 0.0_dp
488 348 : IF (found) qmatac(1:na, 1:nc) = block(1:na, 1:nc)
489 90 : CALL para_env%sum(qmatac)
490 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ic, &
491 90 : block=block, found=found)
492 90 : smatac = 0.0_dp
493 348 : IF (found) smatac(1:na, 1:nc) = block(1:na, 1:nc)
494 90 : CALL para_env%sum(smatac)
495 540 : ALLOCATE (qmatbc(nb, nc), smatbc(nb, nc))
496 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ib, col=ic, &
497 90 : block=block, found=found)
498 90 : qmatbc = 0.0_dp
499 258 : IF (found) qmatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
500 90 : CALL para_env%sum(qmatbc)
501 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ib, col=ic, &
502 90 : block=block, found=found)
503 90 : smatbc = 0.0_dp
504 258 : IF (found) smatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
505 90 : CALL para_env%sum(smatbc)
506 : !
507 90 : nabc = na + nb + nc
508 720 : ALLOCATE (sab(nabc, nabc), sinv(nabc, nabc), qab(nabc, nabc))
509 : !
510 1242 : qab(1:na, 1:na) = qblka(1:na, 1:na)
511 702 : qab(na + 1:nab, na + 1:nab) = qblkb(1:nb, 1:nb)
512 522 : qab(nab + 1:nabc, nab + 1:nabc) = qblkc(1:nc, 1:nc)
513 648 : qab(1:na, na + 1:nab) = qmatab(1:na, 1:nb)
514 738 : qab(na + 1:nab, 1:na) = TRANSPOSE(qmatab(1:na, 1:nb))
515 606 : qab(1:na, nab + 1:nabc) = qmatac(1:na, 1:nc)
516 726 : qab(nab + 1:nabc, 1:na) = TRANSPOSE(qmatac(1:na, 1:nc))
517 426 : qab(na + 1:nab, nab + 1:nabc) = qmatbc(1:nb, 1:nc)
518 456 : qab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(qmatbc(1:nb, 1:nc))
519 : !
520 1242 : sab(1:na, 1:na) = sblka(1:na, 1:na)
521 702 : sab(na + 1:nab, na + 1:nab) = sblkb(1:nb, 1:nb)
522 522 : sab(nab + 1:nabc, nab + 1:nabc) = sblkc(1:nc, 1:nc)
523 648 : sab(1:na, na + 1:nab) = smatab(1:na, 1:nb)
524 738 : sab(na + 1:nab, 1:na) = TRANSPOSE(smatab(1:na, 1:nb))
525 606 : sab(1:na, nab + 1:nabc) = smatac(1:na, 1:nc)
526 726 : sab(nab + 1:nabc, 1:na) = TRANSPOSE(smatac(1:na, 1:nc))
527 426 : sab(na + 1:nab, nab + 1:nabc) = smatbc(1:nb, 1:nc)
528 456 : sab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(smatbc(1:nb, 1:nc))
529 : ! inv smat
530 4254 : sinv(1:nabc, 1:nabc) = sab(1:nabc, 1:nabc)
531 90 : CALL invmat_symm(sinv)
532 : ! Tr(Q*Sinv)
533 90 : iabc = iabc + 1
534 90 : me = MOD(iabc, para_env%num_pe)
535 90 : IF (me == para_env%mepos) THEN
536 2127 : occnumABC(iabc, ispin) = SUM(qab*sinv)
537 : ELSE
538 45 : occnumABC(iabc, ispin) = 0.0_dp
539 : END IF
540 : !
541 90 : DEALLOCATE (sab, sinv, qab)
542 90 : DEALLOCATE (qmatac, smatac)
543 718 : DEALLOCATE (qmatbc, smatbc)
544 : END DO
545 488 : DEALLOCATE (qmatab, smatab)
546 : END DO
547 : END DO
548 10 : CALL dbcsr_release(qmat_diag)
549 18 : CALL dbcsr_release(smat_diag)
550 : END DO
551 8 : CALL para_env%sum(occnumABC)
552 : END IF
553 :
554 10 : IF (.NOT. neglect_abc) THEN
555 : ! calculate shared electron numbers (ABC)
556 8 : nabc = (natom*(natom - 1)*(natom - 2))/6
557 32 : ALLOCATE (selnABC(nabc, nspin))
558 8 : selnABC = 0.0_dp
559 18 : DO ispin = 1, nspin
560 10 : iabc = 0
561 66 : DO ia = 1, natom
562 160 : DO ib = ia + 1, natom
563 274 : DO ic = ib + 1, natom
564 124 : iabc = iabc + 1
565 226 : IF (occnumABC(iabc, ispin) >= 0.0_dp) THEN
566 : selnABC(iabc, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) + occnumA(ic, ispin) - &
567 : occnumAB(ia, ib, ispin) - occnumAB(ia, ic, ispin) - occnumAB(ib, ic, ispin) + &
568 90 : occnumABC(iabc, ispin)
569 : END IF
570 : END DO
571 : END DO
572 : END DO
573 : END DO
574 : END IF
575 :
576 : ! calculate atomic charge
577 40 : ALLOCATE (raq(natom, nspin))
578 10 : raq = 0.0_dp
579 22 : DO ispin = 1, nspin
580 66 : DO ia = 1, natom
581 54 : raq(ia, ispin) = occnumA(ia, ispin)
582 336 : DO ib = 1, natom
583 324 : raq(ia, ispin) = raq(ia, ispin) - 0.5_dp*selnAB(ia, ib, ispin)
584 : END DO
585 : END DO
586 22 : IF (.NOT. neglect_abc) THEN
587 10 : iabc = 0
588 58 : DO ia = 1, natom
589 160 : DO ib = ia + 1, natom
590 274 : DO ic = ib + 1, natom
591 124 : iabc = iabc + 1
592 124 : raq(ia, ispin) = raq(ia, ispin) + selnABC(iabc, ispin)/3._dp
593 124 : raq(ib, ispin) = raq(ib, ispin) + selnABC(iabc, ispin)/3._dp
594 226 : raq(ic, ispin) = raq(ic, ispin) + selnABC(iabc, ispin)/3._dp
595 : END DO
596 : END DO
597 : END DO
598 : END IF
599 : END DO
600 :
601 : ! calculate unassigned charge (from sum over atomic charges)
602 22 : DO ispin = 1, nspin
603 66 : deltaq = (electra(ispin) - SUM(raq(1:natom, ispin))) - ua_charge(ispin)
604 22 : IF (unit_nr > 0) THEN
605 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
606 6 : "Cutoff error on charge", "Spin ", ispin, "error charge =", deltaq
607 : END IF
608 : END DO
609 :
610 : ! analyze unassigned charge
611 40 : ALLOCATE (uaq(natom, nspin))
612 10 : uaq = 0.0_dp
613 10 : IF (analyze_ua) THEN
614 8 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
615 8 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, sab_all=sab_all)
616 : CALL dbcsr_get_info(mao_coef(1)%matrix, row_blk_size=mao_blk_sizes, &
617 8 : col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
618 8 : CALL dbcsr_get_info(matrix_s(1, 1)%matrix, row_blk_size=row_blk_sizes)
619 8 : CALL dbcsr_create(amat, name="temp", template=matrix_s(1, 1)%matrix)
620 8 : CALL dbcsr_create(tmat, name="temp", template=mao_coef(1)%matrix)
621 : ! replicate diagonal of smm matrix
622 8 : CALL dbcsr_get_block_diag(matrix_smm(1)%matrix, smat_diag)
623 8 : CALL dbcsr_replicate_all(smat_diag)
624 :
625 32 : ALLOCATE (orb_blk(natom), mao_blk(natom))
626 50 : DO ia = 1, natom
627 510 : orb_blk = row_blk_sizes
628 510 : mao_blk = row_blk_sizes
629 42 : mao_blk(ia) = col_blk_sizes(ia)
630 : CALL dbcsr_create(sumat, name="Smat", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
631 42 : row_blk_size=mao_blk, col_blk_size=mao_blk)
632 42 : CALL cp_dbcsr_alloc_block_from_nbl(sumat, sab_orb)
633 : CALL dbcsr_create(cholmat, name="Cholesky matrix", dist=dbcsr_dist, &
634 42 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=mao_blk, col_blk_size=mao_blk)
635 : CALL dbcsr_create(rumat, name="Rmat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
636 42 : row_blk_size=orb_blk, col_blk_size=mao_blk)
637 42 : CALL cp_dbcsr_alloc_block_from_nbl(rumat, sab_orb, .TRUE.)
638 : CALL dbcsr_create(crumat, name="Rmat*Umat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
639 42 : row_blk_size=orb_blk, col_blk_size=mao_blk)
640 : ! replicate row and col of smo matrix
641 360 : ALLOCATE (rowblock(natom))
642 276 : DO ib = 1, natom
643 234 : na = mao_blk_sizes(ia)
644 234 : nb = row_blk_sizes(ib)
645 936 : ALLOCATE (rowblock(ib)%mat(na, nb))
646 20396 : rowblock(ib)%mat = 0.0_dp
647 : CALL dbcsr_get_block_p(matrix=matrix_smo(1)%matrix, row=ia, col=ib, &
648 234 : block=block, found=found)
649 10315 : IF (found) rowblock(ib)%mat(1:na, 1:nb) = block(1:na, 1:nb)
650 510 : CALL para_env%sum(rowblock(ib)%mat)
651 : END DO
652 : !
653 90 : DO ispin = 1, nspin
654 48 : CALL dbcsr_copy(tmat, mao_coef(ispin)%matrix)
655 48 : CALL dbcsr_replicate_all(tmat)
656 48 : CALL dbcsr_iterator_start(dbcsr_iter, matrix_s(1, 1)%matrix)
657 462 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
658 414 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, block)
659 414 : CALL dbcsr_get_block_p(matrix=sumat, row=iatom, col=jatom, block=sblk, found=fos)
660 414 : CPASSERT(fos)
661 414 : CALL dbcsr_get_block_p(matrix=rumat, row=iatom, col=jatom, block=rblku, found=for)
662 414 : CPASSERT(for)
663 414 : CALL dbcsr_get_block_p(matrix=rumat, row=jatom, col=iatom, block=rblkl, found=for)
664 414 : CPASSERT(for)
665 414 : CALL dbcsr_get_block_p(matrix=tmat, row=ia, col=ia, block=cmao, found=found)
666 414 : CPASSERT(found)
667 462 : IF (iatom /= ia .AND. jatom /= ia) THEN
668 : ! copy original overlap matrix
669 24864 : sblk = block
670 24864 : rblku = block
671 26008 : rblkl = TRANSPOSE(block)
672 126 : ELSE IF (iatom /= ia) THEN
673 3435 : rblkl = TRANSPOSE(block)
674 51390 : sblk = MATMUL(TRANSPOSE(rowblock(iatom)%mat), cmao)
675 1267 : rblku = sblk
676 75 : ELSE IF (jatom /= ia) THEN
677 3083 : rblku = block
678 45327 : sblk = MATMUL(TRANSPOSE(cmao), rowblock(jatom)%mat)
679 1203 : rblkl = TRANSPOSE(sblk)
680 : ELSE
681 24 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=block, found=found)
682 24 : CPASSERT(found)
683 202604 : sblk = MATMUL(TRANSPOSE(cmao), MATMUL(block, cmao))
684 72928 : rblku = MATMUL(TRANSPOSE(rowblock(ia)%mat), cmao)
685 : END IF
686 : END DO
687 48 : CALL dbcsr_iterator_stop(dbcsr_iter)
688 : ! Cholesky decomposition of SUMAT = U'U
689 48 : CALL dbcsr_desymmetrize(sumat, cholmat)
690 48 : CALL cp_dbcsr_cholesky_decompose(cholmat, para_env=para_env, blacs_env=blacs_env)
691 : ! T = R*inv(U)
692 300 : ssize = SUM(mao_blk)
693 : CALL cp_dbcsr_cholesky_restore(rumat, ssize, cholmat, crumat, op="SOLVE", pos="RIGHT", &
694 48 : transa="N", para_env=para_env, blacs_env=blacs_env)
695 : ! A = T*transpose(T)
696 : CALL dbcsr_multiply("N", "T", 1.0_dp, crumat, crumat, 0.0_dp, amat, &
697 48 : filter_eps=eps_filter)
698 : ! Tr(P*A)
699 48 : CALL dbcsr_dot(matrix_p(ispin, 1)%matrix, amat, uaq(ia, ispin))
700 138 : uaq(ia, ispin) = uaq(ia, ispin) - electra(ispin)
701 : END DO
702 : !
703 42 : CALL dbcsr_release(sumat)
704 42 : CALL dbcsr_release(cholmat)
705 42 : CALL dbcsr_release(rumat)
706 42 : CALL dbcsr_release(crumat)
707 : !
708 276 : DO ib = 1, natom
709 276 : DEALLOCATE (rowblock(ib)%mat)
710 : END DO
711 284 : DEALLOCATE (rowblock)
712 : END DO
713 8 : CALL dbcsr_release(smat_diag)
714 8 : CALL dbcsr_release(amat)
715 8 : CALL dbcsr_release(tmat)
716 16 : DEALLOCATE (orb_blk, mao_blk)
717 : END IF
718 : !
719 76 : raq(1:natom, 1:nspin) = raq(1:natom, 1:nspin) - uaq(1:natom, 1:nspin)
720 22 : DO ispin = 1, nspin
721 66 : deltaq = electra(ispin) - SUM(raq(1:natom, ispin))
722 22 : IF (unit_nr > 0) THEN
723 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
724 6 : "Charge/Atom redistributed", "Spin ", ispin, "delta charge =", &
725 12 : (deltaq + ua_charge(ispin))/REAL(natom, KIND=dp)
726 : END IF
727 : END DO
728 :
729 : ! output charges
730 10 : IF (unit_nr > 0) THEN
731 5 : IF (nspin == 1) THEN
732 4 : WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO atomic charges ", "Atom", "Charge"
733 : ELSE
734 1 : WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO atomic charges ", "Atom", "Charge", "Spin Charge"
735 : END IF
736 11 : DO ispin = 1, nspin
737 33 : deltaq = electra(ispin) - SUM(raq(1:natom, ispin))
738 38 : raq(:, ispin) = raq(:, ispin) + deltaq/REAL(natom, KIND=dp)
739 : END DO
740 5 : total_charge = 0.0_dp
741 5 : total_spin = 0.0_dp
742 29 : DO iatom = 1, natom
743 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
744 24 : element_symbol=element_symbol, kind_number=ikind)
745 24 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
746 29 : IF (nspin == 1) THEN
747 21 : WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, zeff - raq(iatom, 1)
748 21 : total_charge = total_charge + (zeff - raq(iatom, 1))
749 : ELSE
750 3 : WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
751 6 : zeff - raq(iatom, 1) - raq(iatom, 2), raq(iatom, 1) - raq(iatom, 2)
752 3 : total_charge = total_charge + (zeff - raq(iatom, 1) - raq(iatom, 2))
753 3 : total_spin = total_spin + (raq(iatom, 1) - raq(iatom, 2))
754 : END IF
755 : END DO
756 5 : IF (nspin == 1) THEN
757 4 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
758 : ELSE
759 1 : WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin
760 : END IF
761 : END IF
762 :
763 10 : IF (analyze_ua) THEN
764 : ! output unassigned charges
765 8 : IF (unit_nr > 0) THEN
766 4 : IF (nspin == 1) THEN
767 3 : WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO hypervalent charges ", "Atom", "Charge"
768 : ELSE
769 1 : WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO hypervalent charges ", "Atom", &
770 2 : "Charge", "Spin Charge"
771 : END IF
772 4 : total_charge = 0.0_dp
773 4 : total_spin = 0.0_dp
774 25 : DO iatom = 1, natom
775 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
776 21 : element_symbol=element_symbol)
777 25 : IF (nspin == 1) THEN
778 18 : WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, uaq(iatom, 1)
779 18 : total_charge = total_charge + uaq(iatom, 1)
780 : ELSE
781 3 : WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
782 6 : uaq(iatom, 1) + uaq(iatom, 2), uaq(iatom, 1) - uaq(iatom, 2)
783 3 : total_charge = total_charge + uaq(iatom, 1) + uaq(iatom, 2)
784 3 : total_spin = total_spin + uaq(iatom, 1) - uaq(iatom, 2)
785 : END IF
786 : END DO
787 4 : IF (nspin == 1) THEN
788 3 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
789 : ELSE
790 1 : WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin
791 : END IF
792 : END IF
793 : END IF
794 :
795 : ! output shared electron numbers AB
796 10 : IF (unit_nr > 0) THEN
797 5 : IF (nspin == 1) THEN
798 4 : WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T78,A)") "Shared electron numbers ", "Atom", "Atom", "SEN"
799 : ELSE
800 1 : WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T51,A,T63,A,T71,A)") "Shared electron numbers ", "Atom", "Atom", &
801 2 : "SEN(1)", "SEN(2)", "SEN(total)"
802 : END IF
803 29 : DO ia = 1, natom
804 80 : DO ib = ia + 1, natom
805 51 : CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
806 51 : CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
807 75 : IF (nspin == 1) THEN
808 48 : IF (selnAB(ia, ib, 1) > eps_ab) THEN
809 31 : WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T69,F12.6)") ia, esa, ib, esb, selnAB(ia, ib, 1)
810 : END IF
811 : ELSE
812 3 : IF ((selnAB(ia, ib, 1) + selnAB(ia, ib, 2)) > eps_ab) THEN
813 3 : WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T45,3F12.6)") ia, esa, ib, esb, &
814 6 : selnAB(ia, ib, 1), selnAB(ia, ib, 2), (selnAB(ia, ib, 1) + selnAB(ia, ib, 2))
815 : END IF
816 : END IF
817 : END DO
818 : END DO
819 : END IF
820 :
821 10 : IF (.NOT. neglect_abc) THEN
822 : ! output shared electron numbers ABC
823 8 : IF (unit_nr > 0) THEN
824 4 : WRITE (unit_nr, "(/,T2,A,T40,A,T49,A,T58,A,T78,A)") "Shared electron numbers ABC", &
825 8 : "Atom", "Atom", "Atom", "SEN"
826 4 : senmax = 0.0_dp
827 4 : iabc = 0
828 25 : DO ia = 1, natom
829 73 : DO ib = ia + 1, natom
830 130 : DO ic = ib + 1, natom
831 61 : iabc = iabc + 1
832 123 : senabc = SUM(selnABC(iabc, :))
833 61 : senmax = MAX(senmax, senabc)
834 109 : IF (senabc > eps_abc) THEN
835 43 : CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
836 43 : CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
837 43 : CALL get_atomic_kind(atomic_kind=particle_set(ic)%atomic_kind, element_symbol=esc)
838 : WRITE (unit_nr, "(T35,I6,' ',A2,T44,I6,' ',A2,T53,I6,' ',A2,T69,F12.6)") &
839 43 : ia, esa, ib, esb, ic, esc, senabc
840 : END IF
841 : END DO
842 : END DO
843 : END DO
844 4 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Maximum SEN value calculated", senmax
845 : END IF
846 : END IF
847 :
848 10 : IF (print_pao) THEN
849 4 : CALL mao_write_pao_restart(mao_coef, qs_env)
850 : END IF
851 :
852 10 : IF (unit_nr > 0) THEN
853 : WRITE (unit_nr, '(/,T2,A)') &
854 5 : '!---------------------------END OF MAO ANALYSIS-------------------------------!'
855 : END IF
856 :
857 : ! Deallocate temporary arrays
858 10 : DEALLOCATE (occnumA, occnumAB, selnAB, raq, uaq)
859 10 : IF (.NOT. neglect_abc) THEN
860 8 : DEALLOCATE (occnumABC, selnABC)
861 : END IF
862 :
863 : ! Deallocate the neighbor list structure
864 10 : CALL release_neighbor_list_sets(smm_list)
865 10 : CALL release_neighbor_list_sets(smo_list)
866 :
867 10 : DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
868 :
869 10 : IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm)
870 10 : IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo)
871 10 : IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q)
872 :
873 10 : IF (ASSOCIATED(mao_coef)) CALL dbcsr_deallocate_matrix_set(mao_coef)
874 10 : IF (ASSOCIATED(mao_dmat)) CALL dbcsr_deallocate_matrix_set(mao_dmat)
875 10 : IF (ASSOCIATED(mao_smat)) CALL dbcsr_deallocate_matrix_set(mao_smat)
876 10 : IF (ASSOCIATED(mao_qmat)) CALL dbcsr_deallocate_matrix_set(mao_qmat)
877 :
878 10 : CALL timestop(handle)
879 :
880 106 : END SUBROUTINE mao_analysis
881 :
882 24 : END MODULE mao_wfn_analysis
|
