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 58 : IF (col_blk_sizes(iab) < 0) THEN
247 0 : CPABORT("Number of MAOs has to be specified in KIND section for all elements")
248 : END IF
249 : END DO
250 22 : DO ispin = 1, nspin
251 : ! coeficients
252 12 : ALLOCATE (mao_coef(ispin)%matrix)
253 : CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, &
254 : name="MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
255 12 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
256 22 : CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
257 : END DO
258 10 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
259 :
260 : ! optimize MAOs
261 10 : epsx = 1000.0_dp
262 : CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, epsx, &
263 10 : 3, unit_nr)
264 :
265 : ! Analyze the MAO basis
266 : CALL mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
267 10 : qs_kind_set, unit_nr, para_env)
268 :
269 : ! Calculate the overlap and density matrix in the new MAO basis
270 10 : NULLIFY (mao_dmat, mao_smat, mao_qmat)
271 10 : CALL dbcsr_allocate_matrix_set(mao_qmat, nspin)
272 10 : CALL dbcsr_allocate_matrix_set(mao_dmat, nspin)
273 10 : CALL dbcsr_allocate_matrix_set(mao_smat, nspin)
274 10 : CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
275 22 : DO ispin = 1, nspin
276 12 : ALLOCATE (mao_dmat(ispin)%matrix)
277 : CALL dbcsr_create(mao_dmat(ispin)%matrix, name="MAO density", 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_smat(ispin)%matrix)
281 : CALL dbcsr_create(mao_smat(ispin)%matrix, name="MAO overlap", dist=dbcsr_dist, &
282 : matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
283 12 : col_blk_size=col_blk_sizes)
284 12 : ALLOCATE (mao_qmat(ispin)%matrix)
285 : CALL dbcsr_create(mao_qmat(ispin)%matrix, name="MAO covar density", dist=dbcsr_dist, &
286 : matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
287 22 : col_blk_size=col_blk_sizes)
288 : END DO
289 10 : CALL dbcsr_create(amat, name="MAO overlap", template=mao_dmat(1)%matrix)
290 10 : CALL dbcsr_create(tmat, name="MAO Overlap Inverse", template=amat)
291 10 : CALL dbcsr_create(qmat, name="MAO covar density", template=amat)
292 10 : CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
293 10 : CALL dbcsr_create(axmat, name="TEMP", template=amat, matrix_type=dbcsr_type_no_symmetry)
294 22 : DO ispin = 1, nspin
295 : ! calculate MAO overlap matrix
296 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smm(1)%matrix, mao_coef(ispin)%matrix, &
297 12 : 0.0_dp, cgmat)
298 12 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, amat)
299 : ! calculate inverse of MAO overlap
300 12 : threshold = 1.e-8_dp
301 12 : CALL invert_Hotelling(tmat, amat, threshold, norm_convergence=1.e-4_dp, silent=.TRUE.)
302 12 : CALL dbcsr_copy(mao_smat(ispin)%matrix, amat)
303 : ! calculate q-matrix q = C*Q*C
304 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_q(ispin)%matrix, mao_coef(ispin)%matrix, &
305 12 : 0.0_dp, cgmat, filter_eps=eps_filter)
306 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, &
307 12 : 0.0_dp, qmat, filter_eps=eps_filter)
308 12 : CALL dbcsr_copy(mao_qmat(ispin)%matrix, qmat)
309 : ! calculate density matrix
310 12 : CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, tmat, 0.0_dp, axmat, filter_eps=eps_filter)
311 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, axmat, 0.0_dp, mao_dmat(ispin)%matrix, &
312 22 : filter_eps=eps_filter)
313 : END DO
314 10 : CALL dbcsr_release(amat)
315 10 : CALL dbcsr_release(tmat)
316 10 : CALL dbcsr_release(qmat)
317 10 : CALL dbcsr_release(cgmat)
318 10 : CALL dbcsr_release(axmat)
319 :
320 : ! calculate unassigned charge : n - Tr PS
321 22 : DO ispin = 1, nspin
322 12 : CALL dbcsr_dot(mao_dmat(ispin)%matrix, mao_smat(ispin)%matrix, ua_charge(ispin))
323 22 : ua_charge(ispin) = electra(ispin) - ua_charge(ispin)
324 : END DO
325 10 : IF (unit_nr > 0) THEN
326 5 : WRITE (unit_nr, *)
327 11 : DO ispin = 1, nspin
328 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
329 11 : "Unassigned charge", "Spin ", ispin, "delta charge =", ua_charge(ispin)
330 : END DO
331 : END IF
332 :
333 : ! occupation numbers: single atoms
334 : ! We use S_A = 1
335 : ! At the gamma point we use an effective MIC
336 10 : CALL get_qs_env(qs_env, natom=natom)
337 40 : ALLOCATE (occnumA(natom, nspin))
338 10 : occnumA = 0.0_dp
339 22 : DO ispin = 1, nspin
340 76 : DO iatom = 1, natom
341 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
342 54 : row=iatom, col=iatom, block=block, found=found)
343 66 : IF (found) THEN
344 81 : DO iab = 1, SIZE(block, 1)
345 81 : occnumA(iatom, ispin) = occnumA(iatom, ispin) + block(iab, iab)
346 : END DO
347 : END IF
348 : END DO
349 : END DO
350 10 : CALL para_env%sum(occnumA)
351 :
352 : ! occupation numbers: atom pairs
353 50 : ALLOCATE (occnumAB(natom, natom, nspin))
354 10 : occnumAB = 0.0_dp
355 22 : DO ispin = 1, nspin
356 12 : CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
357 12 : CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix)
358 : ! replicate the diagonal blocks of the density and overlap matrices
359 12 : CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
360 12 : CALL dbcsr_replicate_all(qmat_diag)
361 12 : CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
362 12 : CALL dbcsr_replicate_all(smat_diag)
363 66 : DO ia = 1, natom
364 174 : DO ib = ia + 1, natom
365 108 : iab = 0
366 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
367 108 : row=ia, col=ib, block=block, found=found)
368 108 : IF (found) iab = 1
369 108 : CALL para_env%sum(iab)
370 108 : CPASSERT(iab <= 1)
371 270 : IF (iab == 0 .AND. para_env%is_source()) THEN
372 : ! AB block is not available N_AB = N_A + N_B
373 : ! Do this only on the "source" processor
374 0 : occnumAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin)
375 0 : occnumAB(ib, ia, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin)
376 108 : ELSE IF (found) THEN
377 : ! owner of AB block performs calculation
378 54 : na = SIZE(block, 1)
379 54 : nb = SIZE(block, 2)
380 54 : nab = na + nb
381 432 : ALLOCATE (sab(nab, nab), qab(nab, nab), sinv(nab, nab))
382 : ! qmat
383 327 : qab(1:na, na + 1:nab) = block(1:na, 1:nb)
384 375 : qab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb))
385 54 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=diag, found=fo)
386 54 : CPASSERT(fo)
387 630 : qab(1:na, 1:na) = diag(1:na, 1:na)
388 54 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=diag, found=fo)
389 54 : CPASSERT(fo)
390 342 : qab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
391 : ! smat
392 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, &
393 54 : row=ia, col=ib, block=block, found=fo)
394 54 : CPASSERT(fo)
395 327 : sab(1:na, na + 1:nab) = block(1:na, 1:nb)
396 375 : sab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb))
397 54 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=diag, found=fo)
398 54 : CPASSERT(fo)
399 630 : sab(1:na, 1:na) = diag(1:na, 1:na)
400 54 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=diag, found=fo)
401 54 : CPASSERT(fo)
402 342 : sab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
403 : ! inv smat
404 1296 : sinv(1:nab, 1:nab) = sab(1:nab, 1:nab)
405 54 : CALL invmat_symm(sinv)
406 : ! Tr(Q*Sinv)
407 1296 : occnumAB(ia, ib, ispin) = SUM(qab*sinv)
408 54 : occnumAB(ib, ia, ispin) = occnumAB(ia, ib, ispin)
409 : !
410 324 : DEALLOCATE (sab, qab, sinv)
411 : END IF
412 : END DO
413 : END DO
414 12 : CALL dbcsr_release(qmat_diag)
415 22 : CALL dbcsr_release(smat_diag)
416 : END DO
417 10 : CALL para_env%sum(occnumAB)
418 :
419 : ! calculate shared electron numbers (AB)
420 50 : ALLOCATE (selnAB(natom, natom, nspin))
421 10 : selnAB = 0.0_dp
422 22 : DO ispin = 1, nspin
423 76 : DO ia = 1, natom
424 174 : DO ib = ia + 1, natom
425 108 : selnAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) - occnumAB(ia, ib, ispin)
426 162 : selnAB(ib, ia, ispin) = selnAB(ia, ib, ispin)
427 : END DO
428 : END DO
429 : END DO
430 :
431 10 : IF (.NOT. neglect_abc) THEN
432 : ! calculate N_ABC
433 8 : nabc = (natom*(natom - 1)*(natom - 2))/6
434 32 : ALLOCATE (occnumABC(nabc, nspin))
435 142 : occnumABC = -1.0_dp
436 18 : DO ispin = 1, nspin
437 10 : CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
438 10 : CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix)
439 : ! replicate the diagonal blocks of the density and overlap matrices
440 10 : CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
441 10 : CALL dbcsr_replicate_all(qmat_diag)
442 10 : CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
443 10 : CALL dbcsr_replicate_all(smat_diag)
444 10 : iabc = 0
445 58 : DO ia = 1, natom
446 48 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=qblka, found=fo)
447 48 : CPASSERT(fo)
448 48 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=sblka, found=fo)
449 48 : CPASSERT(fo)
450 48 : na = SIZE(qblka, 1)
451 256 : DO ib = ia + 1, natom
452 : ! screen with SEN(AB)
453 102 : IF (selnAB(ia, ib, ispin) < eps_abc) THEN
454 14 : iabc = iabc + (natom - ib)
455 14 : CYCLE
456 : END IF
457 88 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=qblkb, found=fo)
458 88 : CPASSERT(fo)
459 88 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=sblkb, found=fo)
460 88 : CPASSERT(fo)
461 88 : nb = SIZE(qblkb, 1)
462 88 : nab = na + nb
463 528 : ALLOCATE (qmatab(na, nb), smatab(na, nb))
464 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ib, &
465 88 : block=block, found=found)
466 88 : qmatab = 0.0_dp
467 320 : IF (found) qmatab(1:na, 1:nb) = block(1:na, 1:nb)
468 88 : CALL para_env%sum(qmatab)
469 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ib, &
470 88 : block=block, found=found)
471 88 : smatab = 0.0_dp
472 320 : IF (found) smatab(1:na, 1:nb) = block(1:na, 1:nb)
473 88 : CALL para_env%sum(smatab)
474 202 : DO ic = ib + 1, natom
475 : ! screen with SEN(AB)
476 114 : IF ((selnAB(ia, ic, ispin) < eps_abc) .OR. (selnAB(ib, ic, ispin) < eps_abc)) THEN
477 24 : iabc = iabc + 1
478 24 : CYCLE
479 : END IF
480 90 : CALL dbcsr_get_block_p(matrix=qmat_diag, row=ic, col=ic, block=qblkc, found=fo)
481 90 : CPASSERT(fo)
482 90 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ic, col=ic, block=sblkc, found=fo)
483 90 : CPASSERT(fo)
484 90 : nc = SIZE(qblkc, 1)
485 540 : ALLOCATE (qmatac(na, nc), smatac(na, nc))
486 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ic, &
487 90 : block=block, found=found)
488 90 : qmatac = 0.0_dp
489 348 : IF (found) qmatac(1:na, 1:nc) = block(1:na, 1:nc)
490 90 : CALL para_env%sum(qmatac)
491 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ic, &
492 90 : block=block, found=found)
493 90 : smatac = 0.0_dp
494 348 : IF (found) smatac(1:na, 1:nc) = block(1:na, 1:nc)
495 90 : CALL para_env%sum(smatac)
496 540 : ALLOCATE (qmatbc(nb, nc), smatbc(nb, nc))
497 : CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ib, col=ic, &
498 90 : block=block, found=found)
499 90 : qmatbc = 0.0_dp
500 258 : IF (found) qmatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
501 90 : CALL para_env%sum(qmatbc)
502 : CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ib, col=ic, &
503 90 : block=block, found=found)
504 90 : smatbc = 0.0_dp
505 258 : IF (found) smatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
506 90 : CALL para_env%sum(smatbc)
507 : !
508 90 : nabc = na + nb + nc
509 720 : ALLOCATE (sab(nabc, nabc), sinv(nabc, nabc), qab(nabc, nabc))
510 : !
511 1242 : qab(1:na, 1:na) = qblka(1:na, 1:na)
512 702 : qab(na + 1:nab, na + 1:nab) = qblkb(1:nb, 1:nb)
513 522 : qab(nab + 1:nabc, nab + 1:nabc) = qblkc(1:nc, 1:nc)
514 648 : qab(1:na, na + 1:nab) = qmatab(1:na, 1:nb)
515 738 : qab(na + 1:nab, 1:na) = TRANSPOSE(qmatab(1:na, 1:nb))
516 606 : qab(1:na, nab + 1:nabc) = qmatac(1:na, 1:nc)
517 726 : qab(nab + 1:nabc, 1:na) = TRANSPOSE(qmatac(1:na, 1:nc))
518 426 : qab(na + 1:nab, nab + 1:nabc) = qmatbc(1:nb, 1:nc)
519 456 : qab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(qmatbc(1:nb, 1:nc))
520 : !
521 1242 : sab(1:na, 1:na) = sblka(1:na, 1:na)
522 702 : sab(na + 1:nab, na + 1:nab) = sblkb(1:nb, 1:nb)
523 522 : sab(nab + 1:nabc, nab + 1:nabc) = sblkc(1:nc, 1:nc)
524 648 : sab(1:na, na + 1:nab) = smatab(1:na, 1:nb)
525 738 : sab(na + 1:nab, 1:na) = TRANSPOSE(smatab(1:na, 1:nb))
526 606 : sab(1:na, nab + 1:nabc) = smatac(1:na, 1:nc)
527 726 : sab(nab + 1:nabc, 1:na) = TRANSPOSE(smatac(1:na, 1:nc))
528 426 : sab(na + 1:nab, nab + 1:nabc) = smatbc(1:nb, 1:nc)
529 456 : sab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(smatbc(1:nb, 1:nc))
530 : ! inv smat
531 4254 : sinv(1:nabc, 1:nabc) = sab(1:nabc, 1:nabc)
532 90 : CALL invmat_symm(sinv)
533 : ! Tr(Q*Sinv)
534 90 : iabc = iabc + 1
535 90 : me = MOD(iabc, para_env%num_pe)
536 90 : IF (me == para_env%mepos) THEN
537 2127 : occnumABC(iabc, ispin) = SUM(qab*sinv)
538 : ELSE
539 45 : occnumABC(iabc, ispin) = 0.0_dp
540 : END IF
541 : !
542 90 : DEALLOCATE (sab, sinv, qab)
543 90 : DEALLOCATE (qmatac, smatac)
544 718 : DEALLOCATE (qmatbc, smatbc)
545 : END DO
546 488 : DEALLOCATE (qmatab, smatab)
547 : END DO
548 : END DO
549 10 : CALL dbcsr_release(qmat_diag)
550 18 : CALL dbcsr_release(smat_diag)
551 : END DO
552 8 : CALL para_env%sum(occnumABC)
553 : END IF
554 :
555 10 : IF (.NOT. neglect_abc) THEN
556 : ! calculate shared electron numbers (ABC)
557 8 : nabc = (natom*(natom - 1)*(natom - 2))/6
558 32 : ALLOCATE (selnABC(nabc, nspin))
559 8 : selnABC = 0.0_dp
560 18 : DO ispin = 1, nspin
561 10 : iabc = 0
562 66 : DO ia = 1, natom
563 160 : DO ib = ia + 1, natom
564 274 : DO ic = ib + 1, natom
565 124 : iabc = iabc + 1
566 226 : IF (occnumABC(iabc, ispin) >= 0.0_dp) THEN
567 : selnABC(iabc, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) + occnumA(ic, ispin) - &
568 : occnumAB(ia, ib, ispin) - occnumAB(ia, ic, ispin) - occnumAB(ib, ic, ispin) + &
569 90 : occnumABC(iabc, ispin)
570 : END IF
571 : END DO
572 : END DO
573 : END DO
574 : END DO
575 : END IF
576 :
577 : ! calculate atomic charge
578 40 : ALLOCATE (raq(natom, nspin))
579 10 : raq = 0.0_dp
580 22 : DO ispin = 1, nspin
581 66 : DO ia = 1, natom
582 54 : raq(ia, ispin) = occnumA(ia, ispin)
583 336 : DO ib = 1, natom
584 324 : raq(ia, ispin) = raq(ia, ispin) - 0.5_dp*selnAB(ia, ib, ispin)
585 : END DO
586 : END DO
587 22 : IF (.NOT. neglect_abc) THEN
588 10 : iabc = 0
589 58 : DO ia = 1, natom
590 160 : DO ib = ia + 1, natom
591 274 : DO ic = ib + 1, natom
592 124 : iabc = iabc + 1
593 124 : raq(ia, ispin) = raq(ia, ispin) + selnABC(iabc, ispin)/3._dp
594 124 : raq(ib, ispin) = raq(ib, ispin) + selnABC(iabc, ispin)/3._dp
595 226 : raq(ic, ispin) = raq(ic, ispin) + selnABC(iabc, ispin)/3._dp
596 : END DO
597 : END DO
598 : END DO
599 : END IF
600 : END DO
601 :
602 : ! calculate unassigned charge (from sum over atomic charges)
603 22 : DO ispin = 1, nspin
604 66 : deltaq = (electra(ispin) - SUM(raq(1:natom, ispin))) - ua_charge(ispin)
605 22 : IF (unit_nr > 0) THEN
606 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
607 6 : "Cutoff error on charge", "Spin ", ispin, "error charge =", deltaq
608 : END IF
609 : END DO
610 :
611 : ! analyze unassigned charge
612 40 : ALLOCATE (uaq(natom, nspin))
613 10 : uaq = 0.0_dp
614 10 : IF (analyze_ua) THEN
615 8 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
616 8 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, sab_all=sab_all)
617 : CALL dbcsr_get_info(mao_coef(1)%matrix, row_blk_size=mao_blk_sizes, &
618 8 : col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
619 8 : CALL dbcsr_get_info(matrix_s(1, 1)%matrix, row_blk_size=row_blk_sizes)
620 8 : CALL dbcsr_create(amat, name="temp", template=matrix_s(1, 1)%matrix)
621 8 : CALL dbcsr_create(tmat, name="temp", template=mao_coef(1)%matrix)
622 : ! replicate diagonal of smm matrix
623 8 : CALL dbcsr_get_block_diag(matrix_smm(1)%matrix, smat_diag)
624 8 : CALL dbcsr_replicate_all(smat_diag)
625 :
626 32 : ALLOCATE (orb_blk(natom), mao_blk(natom))
627 50 : DO ia = 1, natom
628 510 : orb_blk = row_blk_sizes
629 510 : mao_blk = row_blk_sizes
630 42 : mao_blk(ia) = col_blk_sizes(ia)
631 : CALL dbcsr_create(sumat, name="Smat", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
632 42 : row_blk_size=mao_blk, col_blk_size=mao_blk)
633 42 : CALL cp_dbcsr_alloc_block_from_nbl(sumat, sab_orb)
634 : CALL dbcsr_create(cholmat, name="Cholesky matrix", dist=dbcsr_dist, &
635 42 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=mao_blk, col_blk_size=mao_blk)
636 : CALL dbcsr_create(rumat, name="Rmat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
637 42 : row_blk_size=orb_blk, col_blk_size=mao_blk)
638 42 : CALL cp_dbcsr_alloc_block_from_nbl(rumat, sab_orb, .TRUE.)
639 : CALL dbcsr_create(crumat, name="Rmat*Umat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
640 42 : row_blk_size=orb_blk, col_blk_size=mao_blk)
641 : ! replicate row and col of smo matrix
642 360 : ALLOCATE (rowblock(natom))
643 276 : DO ib = 1, natom
644 234 : na = mao_blk_sizes(ia)
645 234 : nb = row_blk_sizes(ib)
646 936 : ALLOCATE (rowblock(ib)%mat(na, nb))
647 20396 : rowblock(ib)%mat = 0.0_dp
648 : CALL dbcsr_get_block_p(matrix=matrix_smo(1)%matrix, row=ia, col=ib, &
649 234 : block=block, found=found)
650 10315 : IF (found) rowblock(ib)%mat(1:na, 1:nb) = block(1:na, 1:nb)
651 510 : CALL para_env%sum(rowblock(ib)%mat)
652 : END DO
653 : !
654 90 : DO ispin = 1, nspin
655 48 : CALL dbcsr_copy(tmat, mao_coef(ispin)%matrix)
656 48 : CALL dbcsr_replicate_all(tmat)
657 48 : CALL dbcsr_iterator_start(dbcsr_iter, matrix_s(1, 1)%matrix)
658 462 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
659 414 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, block)
660 414 : CALL dbcsr_get_block_p(matrix=sumat, row=iatom, col=jatom, block=sblk, found=fos)
661 414 : CPASSERT(fos)
662 414 : CALL dbcsr_get_block_p(matrix=rumat, row=iatom, col=jatom, block=rblku, found=for)
663 414 : CPASSERT(for)
664 414 : CALL dbcsr_get_block_p(matrix=rumat, row=jatom, col=iatom, block=rblkl, found=for)
665 414 : CPASSERT(for)
666 414 : CALL dbcsr_get_block_p(matrix=tmat, row=ia, col=ia, block=cmao, found=found)
667 414 : CPASSERT(found)
668 462 : IF (iatom /= ia .AND. jatom /= ia) THEN
669 : ! copy original overlap matrix
670 24864 : sblk = block
671 24864 : rblku = block
672 26008 : rblkl = TRANSPOSE(block)
673 126 : ELSE IF (iatom /= ia) THEN
674 3435 : rblkl = TRANSPOSE(block)
675 51390 : sblk = MATMUL(TRANSPOSE(rowblock(iatom)%mat), cmao)
676 1267 : rblku = sblk
677 75 : ELSE IF (jatom /= ia) THEN
678 3083 : rblku = block
679 45327 : sblk = MATMUL(TRANSPOSE(cmao), rowblock(jatom)%mat)
680 1203 : rblkl = TRANSPOSE(sblk)
681 : ELSE
682 24 : CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=block, found=found)
683 24 : CPASSERT(found)
684 202604 : sblk = MATMUL(TRANSPOSE(cmao), MATMUL(block, cmao))
685 72928 : rblku = MATMUL(TRANSPOSE(rowblock(ia)%mat), cmao)
686 : END IF
687 : END DO
688 48 : CALL dbcsr_iterator_stop(dbcsr_iter)
689 : ! Cholesky decomposition of SUMAT = U'U
690 48 : CALL dbcsr_desymmetrize(sumat, cholmat)
691 48 : CALL cp_dbcsr_cholesky_decompose(cholmat, para_env=para_env, blacs_env=blacs_env)
692 : ! T = R*inv(U)
693 300 : ssize = SUM(mao_blk)
694 : CALL cp_dbcsr_cholesky_restore(rumat, ssize, cholmat, crumat, op="SOLVE", pos="RIGHT", &
695 48 : transa="N", para_env=para_env, blacs_env=blacs_env)
696 : ! A = T*transpose(T)
697 : CALL dbcsr_multiply("N", "T", 1.0_dp, crumat, crumat, 0.0_dp, amat, &
698 48 : filter_eps=eps_filter)
699 : ! Tr(P*A)
700 48 : CALL dbcsr_dot(matrix_p(ispin, 1)%matrix, amat, uaq(ia, ispin))
701 138 : uaq(ia, ispin) = uaq(ia, ispin) - electra(ispin)
702 : END DO
703 : !
704 42 : CALL dbcsr_release(sumat)
705 42 : CALL dbcsr_release(cholmat)
706 42 : CALL dbcsr_release(rumat)
707 42 : CALL dbcsr_release(crumat)
708 : !
709 276 : DO ib = 1, natom
710 276 : DEALLOCATE (rowblock(ib)%mat)
711 : END DO
712 284 : DEALLOCATE (rowblock)
713 : END DO
714 8 : CALL dbcsr_release(smat_diag)
715 8 : CALL dbcsr_release(amat)
716 8 : CALL dbcsr_release(tmat)
717 16 : DEALLOCATE (orb_blk, mao_blk)
718 : END IF
719 : !
720 76 : raq(1:natom, 1:nspin) = raq(1:natom, 1:nspin) - uaq(1:natom, 1:nspin)
721 22 : DO ispin = 1, nspin
722 66 : deltaq = electra(ispin) - SUM(raq(1:natom, ispin))
723 22 : IF (unit_nr > 0) THEN
724 : WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") &
725 6 : "Charge/Atom redistributed", "Spin ", ispin, "delta charge =", &
726 12 : (deltaq + ua_charge(ispin))/REAL(natom, KIND=dp)
727 : END IF
728 : END DO
729 :
730 : ! output charges
731 10 : IF (unit_nr > 0) THEN
732 5 : IF (nspin == 1) THEN
733 4 : WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO atomic charges ", "Atom", "Charge"
734 : ELSE
735 1 : WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO atomic charges ", "Atom", "Charge", "Spin Charge"
736 : END IF
737 11 : DO ispin = 1, nspin
738 33 : deltaq = electra(ispin) - SUM(raq(1:natom, ispin))
739 38 : raq(:, ispin) = raq(:, ispin) + deltaq/REAL(natom, KIND=dp)
740 : END DO
741 5 : total_charge = 0.0_dp
742 5 : total_spin = 0.0_dp
743 29 : DO iatom = 1, natom
744 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
745 24 : element_symbol=element_symbol, kind_number=ikind)
746 24 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
747 29 : IF (nspin == 1) THEN
748 21 : WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, zeff - raq(iatom, 1)
749 21 : total_charge = total_charge + (zeff - raq(iatom, 1))
750 : ELSE
751 3 : WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
752 6 : zeff - raq(iatom, 1) - raq(iatom, 2), raq(iatom, 1) - raq(iatom, 2)
753 3 : total_charge = total_charge + (zeff - raq(iatom, 1) - raq(iatom, 2))
754 3 : total_spin = total_spin + (raq(iatom, 1) - raq(iatom, 2))
755 : END IF
756 : END DO
757 5 : IF (nspin == 1) THEN
758 4 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
759 : ELSE
760 1 : WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin
761 : END IF
762 : END IF
763 :
764 10 : IF (analyze_ua) THEN
765 : ! output unassigned charges
766 8 : IF (unit_nr > 0) THEN
767 4 : IF (nspin == 1) THEN
768 3 : WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO hypervalent charges ", "Atom", "Charge"
769 : ELSE
770 1 : WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO hypervalent charges ", "Atom", &
771 2 : "Charge", "Spin Charge"
772 : END IF
773 4 : total_charge = 0.0_dp
774 4 : total_spin = 0.0_dp
775 25 : DO iatom = 1, natom
776 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
777 21 : element_symbol=element_symbol)
778 25 : IF (nspin == 1) THEN
779 18 : WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, uaq(iatom, 1)
780 18 : total_charge = total_charge + uaq(iatom, 1)
781 : ELSE
782 3 : WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
783 6 : uaq(iatom, 1) + uaq(iatom, 2), uaq(iatom, 1) - uaq(iatom, 2)
784 3 : total_charge = total_charge + uaq(iatom, 1) + uaq(iatom, 2)
785 3 : total_spin = total_spin + uaq(iatom, 1) - uaq(iatom, 2)
786 : END IF
787 : END DO
788 4 : IF (nspin == 1) THEN
789 3 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
790 : ELSE
791 1 : WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin
792 : END IF
793 : END IF
794 : END IF
795 :
796 : ! output shared electron numbers AB
797 10 : IF (unit_nr > 0) THEN
798 5 : IF (nspin == 1) THEN
799 4 : WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T78,A)") "Shared electron numbers ", "Atom", "Atom", "SEN"
800 : ELSE
801 1 : WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T51,A,T63,A,T71,A)") "Shared electron numbers ", "Atom", "Atom", &
802 2 : "SEN(1)", "SEN(2)", "SEN(total)"
803 : END IF
804 29 : DO ia = 1, natom
805 80 : DO ib = ia + 1, natom
806 51 : CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
807 51 : CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
808 75 : IF (nspin == 1) THEN
809 48 : IF (selnAB(ia, ib, 1) > eps_ab) THEN
810 31 : WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T69,F12.6)") ia, esa, ib, esb, selnAB(ia, ib, 1)
811 : END IF
812 : ELSE
813 3 : IF ((selnAB(ia, ib, 1) + selnAB(ia, ib, 2)) > eps_ab) THEN
814 3 : WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T45,3F12.6)") ia, esa, ib, esb, &
815 6 : selnAB(ia, ib, 1), selnAB(ia, ib, 2), (selnAB(ia, ib, 1) + selnAB(ia, ib, 2))
816 : END IF
817 : END IF
818 : END DO
819 : END DO
820 : END IF
821 :
822 10 : IF (.NOT. neglect_abc) THEN
823 : ! output shared electron numbers ABC
824 8 : IF (unit_nr > 0) THEN
825 4 : WRITE (unit_nr, "(/,T2,A,T40,A,T49,A,T58,A,T78,A)") "Shared electron numbers ABC", &
826 8 : "Atom", "Atom", "Atom", "SEN"
827 4 : senmax = 0.0_dp
828 4 : iabc = 0
829 25 : DO ia = 1, natom
830 73 : DO ib = ia + 1, natom
831 130 : DO ic = ib + 1, natom
832 61 : iabc = iabc + 1
833 123 : senabc = SUM(selnABC(iabc, :))
834 61 : senmax = MAX(senmax, senabc)
835 109 : IF (senabc > eps_abc) THEN
836 43 : CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
837 43 : CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
838 43 : CALL get_atomic_kind(atomic_kind=particle_set(ic)%atomic_kind, element_symbol=esc)
839 : WRITE (unit_nr, "(T35,I6,' ',A2,T44,I6,' ',A2,T53,I6,' ',A2,T69,F12.6)") &
840 43 : ia, esa, ib, esb, ic, esc, senabc
841 : END IF
842 : END DO
843 : END DO
844 : END DO
845 4 : WRITE (unit_nr, "(T2,A,T69,F12.6)") "Maximum SEN value calculated", senmax
846 : END IF
847 : END IF
848 :
849 10 : IF (print_pao) THEN
850 4 : CALL mao_write_pao_restart(mao_coef, qs_env)
851 : END IF
852 :
853 10 : IF (unit_nr > 0) THEN
854 : WRITE (unit_nr, '(/,T2,A)') &
855 5 : '!---------------------------END OF MAO ANALYSIS-------------------------------!'
856 : END IF
857 :
858 : ! Deallocate temporary arrays
859 10 : DEALLOCATE (occnumA, occnumAB, selnAB, raq, uaq)
860 10 : IF (.NOT. neglect_abc) THEN
861 8 : DEALLOCATE (occnumABC, selnABC)
862 : END IF
863 :
864 : ! Deallocate the neighbor list structure
865 10 : CALL release_neighbor_list_sets(smm_list)
866 10 : CALL release_neighbor_list_sets(smo_list)
867 :
868 10 : DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
869 :
870 10 : IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm)
871 10 : IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo)
872 10 : IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q)
873 :
874 10 : IF (ASSOCIATED(mao_coef)) CALL dbcsr_deallocate_matrix_set(mao_coef)
875 10 : IF (ASSOCIATED(mao_dmat)) CALL dbcsr_deallocate_matrix_set(mao_dmat)
876 10 : IF (ASSOCIATED(mao_smat)) CALL dbcsr_deallocate_matrix_set(mao_smat)
877 10 : IF (ASSOCIATED(mao_qmat)) CALL dbcsr_deallocate_matrix_set(mao_qmat)
878 :
879 10 : CALL timestop(handle)
880 :
881 106 : END SUBROUTINE mao_analysis
882 :
883 24 : END MODULE mao_wfn_analysis
|
