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