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