(git:07a6c39)
Loading...
Searching...
No Matches
mao_wfn_analysis.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
84CONTAINS
85
86! **************************************************************************************************
87!> \brief ...
88!> \param qs_env ...
89!> \param input_section ...
90!> \param unit_nr ...
91! **************************************************************************************************
92 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 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, mao_blk_sizes, &
103 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: occnuma, occnumabc, qab, qmatab, qmatac, &
110 qmatbc, raq, sab, selnabc, sinv, &
111 smatab, smatac, smatbc, uaq
112 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: occnumab, selnab
113 REAL(kind=dp), DIMENSION(:, :), POINTER :: block, cmao, diag, qblka, qblkb, qblkc, &
114 rblkl, rblku, sblk, sblka, sblkb, sblkc
115 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, mao_dmat, mao_qmat, mao_smat, &
120 matrix_q, matrix_smm, matrix_smo
121 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 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 DIMENSION(:), POINTER :: nl_iterator
131 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
132 POINTER :: sab_all, sab_orb, smm_list, smo_list
133 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
134 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 CALL section_vals_get(input_section, explicit=explicit)
141 IF (.NOT. explicit) RETURN
142
143 CALL timeset(routinen, handle)
144
145 IF (unit_nr > 0) THEN
146 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
147 WRITE (unit=unit_nr, fmt="(T36,A)") "MAO ANALYSIS"
148 WRITE (unit=unit_nr, fmt="(T12,A)") "Claus Ehrhardt and Reinhart Ahlrichs, TCA 68:231-245 (1985)"
149 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
150 END IF
151 CALL cite_reference(heinzmann1976)
152 CALL cite_reference(ehrhardt1985)
153
154 ! input options
155 CALL section_vals_val_get(input_section, "REFERENCE_BASIS", i_val=mao_basis)
156 CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter)
157 CALL section_vals_val_get(input_section, "EPS_FUNCTION", r_val=eps_fun)
158 CALL section_vals_val_get(input_section, "EPS_GRAD", r_val=eps_grad)
159 CALL section_vals_val_get(input_section, "MAX_ITER", i_val=max_iter)
160 CALL section_vals_val_get(input_section, "PRINT_BASIS", l_val=print_basis)
161 CALL section_vals_val_get(input_section, "PRINT_PAO", l_val=print_pao)
162 CALL section_vals_val_get(input_section, "NEGLECT_ABC", l_val=neglect_abc)
163 CALL section_vals_val_get(input_section, "AB_THRESHOLD", r_val=eps_ab)
164 CALL section_vals_val_get(input_section, "ABC_THRESHOLD", r_val=eps_abc)
165 CALL section_vals_val_get(input_section, "ANALYZE_UNASSIGNED_CHARGE", l_val=analyze_ua)
166
167 ! k-points?
168 CALL get_qs_env(qs_env, dft_control=dft_control)
169 nimages = dft_control%nimages
170 IF (nimages > 1) THEN
171 IF (unit_nr > 0) THEN
172 WRITE (unit=unit_nr, fmt="(T2,A)") &
173 "K-Points: MAO's determined and analyzed using Gamma-Point only."
174 END IF
175 END IF
176
177 ! Reference basis set
178 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 unit_nr, print_basis)
181
182 ! neighbor lists
183 NULLIFY (smm_list, smo_list)
184 CALL setup_neighbor_list(smm_list, mao_basis_set_list, qs_env=qs_env)
185 CALL setup_neighbor_list(smo_list, mao_basis_set_list, orb_basis_set_list, qs_env=qs_env)
186
187 ! overlap matrices
188 NULLIFY (matrix_smm, matrix_smo)
189 CALL get_qs_env(qs_env, ks_env=ks_env)
190 CALL build_overlap_matrix_simple(ks_env, matrix_smm, &
191 mao_basis_set_list, mao_basis_set_list, smm_list)
192 CALL build_overlap_matrix_simple(ks_env, matrix_smo, &
193 mao_basis_set_list, orb_basis_set_list, smo_list)
194
195 ! get reference density matrix and overlap matrix
196 CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s)
197 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
198 nspin = SIZE(matrix_p, 1)
199 !
200 ! Q matrix
201 IF (nimages == 1) THEN
202 CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter)
203 ELSE
204 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 nimages=nimages, kpoints=kpoints, matrix_ks=matrix_ks, sab_orb=sab_orb)
207 END IF
208
209 ! check for extended basis sets
210 fall = 0
211 CALL neighbor_list_iterator_create(nl_iterator, smm_list)
212 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
213 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
214 IF (iatom <= jatom) THEN
215 irow = iatom
216 icol = jatom
217 ELSE
218 irow = jatom
219 icol = iatom
220 END IF
221 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
222 row=irow, col=icol, block=block, found=found)
223 IF (.NOT. found) fall = fall + 1
224 END DO
225 CALL neighbor_list_iterator_release(nl_iterator)
226
227 CALL get_qs_env(qs_env=qs_env, para_env=para_env)
228 CALL para_env%sum(fall)
229 IF (unit_nr > 0 .AND. fall > 0) THEN
230 WRITE (unit=unit_nr, fmt="(/,T2,A,/,T2,A,/)") &
231 "Warning: Extended MAO basis used with original basis filtered density matrix", &
232 "Warning: Possible errors can be controlled with EPS_PGF_ORB"
233 END IF
234
235 ! MAO matrices
236 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
237 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
238 NULLIFY (mao_coef)
239 CALL dbcsr_allocate_matrix_set(mao_coef, nspin)
240 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
241 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
242 basis=mao_basis_set_list)
243 CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
244 ! check if MAOs have been specified
245 DO iab = 1, natom
246 IF (col_blk_sizes(iab) < 0) THEN
247 cpabort("Number of MAOs has to be specified in KIND section for all elements")
248 END IF
249 END DO
250 DO ispin = 1, nspin
251 ! coeficients
252 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 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
256 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
257 END DO
258 DEALLOCATE (row_blk_sizes, col_blk_sizes)
259
260 ! optimize MAOs
261 epsx = 1000.0_dp
262 CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, epsx, &
263 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 qs_kind_set, unit_nr, para_env)
268
269 ! Calculate the overlap and density matrix in the new MAO basis
270 NULLIFY (mao_dmat, mao_smat, mao_qmat)
271 CALL dbcsr_allocate_matrix_set(mao_qmat, nspin)
272 CALL dbcsr_allocate_matrix_set(mao_dmat, nspin)
273 CALL dbcsr_allocate_matrix_set(mao_smat, nspin)
274 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
275 DO ispin = 1, nspin
276 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 col_blk_size=col_blk_sizes)
280 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 col_blk_size=col_blk_sizes)
284 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 col_blk_size=col_blk_sizes)
288 END DO
289 CALL dbcsr_create(amat, name="MAO overlap", template=mao_dmat(1)%matrix)
290 CALL dbcsr_create(tmat, name="MAO Overlap Inverse", template=amat)
291 CALL dbcsr_create(qmat, name="MAO covar density", template=amat)
292 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
293 CALL dbcsr_create(axmat, name="TEMP", template=amat, matrix_type=dbcsr_type_no_symmetry)
294 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 0.0_dp, cgmat)
298 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, amat)
299 ! calculate inverse of MAO overlap
300 threshold = 1.e-8_dp
301 CALL invert_hotelling(tmat, amat, threshold, norm_convergence=1.e-4_dp, silent=.true.)
302 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 0.0_dp, cgmat, filter_eps=eps_filter)
306 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, &
307 0.0_dp, qmat, filter_eps=eps_filter)
308 CALL dbcsr_copy(mao_qmat(ispin)%matrix, qmat)
309 ! calculate density matrix
310 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 filter_eps=eps_filter)
313 END DO
314 CALL dbcsr_release(amat)
315 CALL dbcsr_release(tmat)
316 CALL dbcsr_release(qmat)
317 CALL dbcsr_release(cgmat)
318 CALL dbcsr_release(axmat)
319
320 ! calculate unassigned charge : n - Tr PS
321 DO ispin = 1, nspin
322 CALL dbcsr_dot(mao_dmat(ispin)%matrix, mao_smat(ispin)%matrix, ua_charge(ispin))
323 ua_charge(ispin) = electra(ispin) - ua_charge(ispin)
324 END DO
325 IF (unit_nr > 0) THEN
326 WRITE (unit_nr, *)
327 DO ispin = 1, nspin
328 WRITE (unit=unit_nr, fmt="(T2,A,T32,A,i2,T55,A,F12.8)") &
329 "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 CALL get_qs_env(qs_env, natom=natom)
337 ALLOCATE (occnuma(natom, nspin))
338 occnuma = 0.0_dp
339 DO ispin = 1, nspin
340 DO iatom = 1, natom
341 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
342 row=iatom, col=iatom, block=block, found=found)
343 IF (found) THEN
344 DO iab = 1, SIZE(block, 1)
345 occnuma(iatom, ispin) = occnuma(iatom, ispin) + block(iab, iab)
346 END DO
347 END IF
348 END DO
349 END DO
350 CALL para_env%sum(occnuma)
351
352 ! occupation numbers: atom pairs
353 ALLOCATE (occnumab(natom, natom, nspin))
354 occnumab = 0.0_dp
355 DO ispin = 1, nspin
356 CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
357 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 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
360 CALL dbcsr_replicate_all(qmat_diag)
361 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
362 CALL dbcsr_replicate_all(smat_diag)
363 DO ia = 1, natom
364 DO ib = ia + 1, natom
365 iab = 0
366 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
367 row=ia, col=ib, block=block, found=found)
368 IF (found) iab = 1
369 CALL para_env%sum(iab)
370 cpassert(iab <= 1)
371 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 occnumab(ia, ib, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin)
375 occnumab(ib, ia, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin)
376 ELSE IF (found) THEN
377 ! owner of AB block performs calculation
378 na = SIZE(block, 1)
379 nb = SIZE(block, 2)
380 nab = na + nb
381 ALLOCATE (sab(nab, nab), qab(nab, nab), sinv(nab, nab))
382 ! qmat
383 qab(1:na, na + 1:nab) = block(1:na, 1:nb)
384 qab(na + 1:nab, 1:na) = transpose(block(1:na, 1:nb))
385 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=diag, found=fo)
386 cpassert(fo)
387 qab(1:na, 1:na) = diag(1:na, 1:na)
388 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=diag, found=fo)
389 cpassert(fo)
390 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 row=ia, col=ib, block=block, found=fo)
394 cpassert(fo)
395 sab(1:na, na + 1:nab) = block(1:na, 1:nb)
396 sab(na + 1:nab, 1:na) = transpose(block(1:na, 1:nb))
397 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=diag, found=fo)
398 cpassert(fo)
399 sab(1:na, 1:na) = diag(1:na, 1:na)
400 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=diag, found=fo)
401 cpassert(fo)
402 sab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
403 ! inv smat
404 sinv(1:nab, 1:nab) = sab(1:nab, 1:nab)
405 CALL invmat_symm(sinv)
406 ! Tr(Q*Sinv)
407 occnumab(ia, ib, ispin) = sum(qab*sinv)
408 occnumab(ib, ia, ispin) = occnumab(ia, ib, ispin)
409 !
410 DEALLOCATE (sab, qab, sinv)
411 END IF
412 END DO
413 END DO
414 CALL dbcsr_release(qmat_diag)
415 CALL dbcsr_release(smat_diag)
416 END DO
417 CALL para_env%sum(occnumab)
418
419 ! calculate shared electron numbers (AB)
420 ALLOCATE (selnab(natom, natom, nspin))
421 selnab = 0.0_dp
422 DO ispin = 1, nspin
423 DO ia = 1, natom
424 DO ib = ia + 1, natom
425 selnab(ia, ib, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin) - occnumab(ia, ib, ispin)
426 selnab(ib, ia, ispin) = selnab(ia, ib, ispin)
427 END DO
428 END DO
429 END DO
430
431 IF (.NOT. neglect_abc) THEN
432 ! calculate N_ABC
433 nabc = (natom*(natom - 1)*(natom - 2))/6
434 ALLOCATE (occnumabc(nabc, nspin))
435 occnumabc = -1.0_dp
436 DO ispin = 1, nspin
437 CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix)
438 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 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
441 CALL dbcsr_replicate_all(qmat_diag)
442 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
443 CALL dbcsr_replicate_all(smat_diag)
444 iabc = 0
445 DO ia = 1, natom
446 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=qblka, found=fo)
447 cpassert(fo)
448 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=sblka, found=fo)
449 cpassert(fo)
450 na = SIZE(qblka, 1)
451 DO ib = ia + 1, natom
452 ! screen with SEN(AB)
453 IF (selnab(ia, ib, ispin) < eps_abc) THEN
454 iabc = iabc + (natom - ib)
455 cycle
456 END IF
457 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=qblkb, found=fo)
458 cpassert(fo)
459 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=sblkb, found=fo)
460 cpassert(fo)
461 nb = SIZE(qblkb, 1)
462 nab = na + nb
463 ALLOCATE (qmatab(na, nb), smatab(na, nb))
464 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ib, &
465 block=block, found=found)
466 qmatab = 0.0_dp
467 IF (found) qmatab(1:na, 1:nb) = block(1:na, 1:nb)
468 CALL para_env%sum(qmatab)
469 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ib, &
470 block=block, found=found)
471 smatab = 0.0_dp
472 IF (found) smatab(1:na, 1:nb) = block(1:na, 1:nb)
473 CALL para_env%sum(smatab)
474 DO ic = ib + 1, natom
475 ! screen with SEN(AB)
476 IF ((selnab(ia, ic, ispin) < eps_abc) .OR. (selnab(ib, ic, ispin) < eps_abc)) THEN
477 iabc = iabc + 1
478 cycle
479 END IF
480 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ic, col=ic, block=qblkc, found=fo)
481 cpassert(fo)
482 CALL dbcsr_get_block_p(matrix=smat_diag, row=ic, col=ic, block=sblkc, found=fo)
483 cpassert(fo)
484 nc = SIZE(qblkc, 1)
485 ALLOCATE (qmatac(na, nc), smatac(na, nc))
486 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ic, &
487 block=block, found=found)
488 qmatac = 0.0_dp
489 IF (found) qmatac(1:na, 1:nc) = block(1:na, 1:nc)
490 CALL para_env%sum(qmatac)
491 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ic, &
492 block=block, found=found)
493 smatac = 0.0_dp
494 IF (found) smatac(1:na, 1:nc) = block(1:na, 1:nc)
495 CALL para_env%sum(smatac)
496 ALLOCATE (qmatbc(nb, nc), smatbc(nb, nc))
497 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ib, col=ic, &
498 block=block, found=found)
499 qmatbc = 0.0_dp
500 IF (found) qmatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
501 CALL para_env%sum(qmatbc)
502 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ib, col=ic, &
503 block=block, found=found)
504 smatbc = 0.0_dp
505 IF (found) smatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
506 CALL para_env%sum(smatbc)
507 !
508 nabc = na + nb + nc
509 ALLOCATE (sab(nabc, nabc), sinv(nabc, nabc), qab(nabc, nabc))
510 !
511 qab(1:na, 1:na) = qblka(1:na, 1:na)
512 qab(na + 1:nab, na + 1:nab) = qblkb(1:nb, 1:nb)
513 qab(nab + 1:nabc, nab + 1:nabc) = qblkc(1:nc, 1:nc)
514 qab(1:na, na + 1:nab) = qmatab(1:na, 1:nb)
515 qab(na + 1:nab, 1:na) = transpose(qmatab(1:na, 1:nb))
516 qab(1:na, nab + 1:nabc) = qmatac(1:na, 1:nc)
517 qab(nab + 1:nabc, 1:na) = transpose(qmatac(1:na, 1:nc))
518 qab(na + 1:nab, nab + 1:nabc) = qmatbc(1:nb, 1:nc)
519 qab(nab + 1:nabc, na + 1:nab) = transpose(qmatbc(1:nb, 1:nc))
520 !
521 sab(1:na, 1:na) = sblka(1:na, 1:na)
522 sab(na + 1:nab, na + 1:nab) = sblkb(1:nb, 1:nb)
523 sab(nab + 1:nabc, nab + 1:nabc) = sblkc(1:nc, 1:nc)
524 sab(1:na, na + 1:nab) = smatab(1:na, 1:nb)
525 sab(na + 1:nab, 1:na) = transpose(smatab(1:na, 1:nb))
526 sab(1:na, nab + 1:nabc) = smatac(1:na, 1:nc)
527 sab(nab + 1:nabc, 1:na) = transpose(smatac(1:na, 1:nc))
528 sab(na + 1:nab, nab + 1:nabc) = smatbc(1:nb, 1:nc)
529 sab(nab + 1:nabc, na + 1:nab) = transpose(smatbc(1:nb, 1:nc))
530 ! inv smat
531 sinv(1:nabc, 1:nabc) = sab(1:nabc, 1:nabc)
532 CALL invmat_symm(sinv)
533 ! Tr(Q*Sinv)
534 iabc = iabc + 1
535 me = mod(iabc, para_env%num_pe)
536 IF (me == para_env%mepos) THEN
537 occnumabc(iabc, ispin) = sum(qab*sinv)
538 ELSE
539 occnumabc(iabc, ispin) = 0.0_dp
540 END IF
541 !
542 DEALLOCATE (sab, sinv, qab)
543 DEALLOCATE (qmatac, smatac)
544 DEALLOCATE (qmatbc, smatbc)
545 END DO
546 DEALLOCATE (qmatab, smatab)
547 END DO
548 END DO
549 CALL dbcsr_release(qmat_diag)
550 CALL dbcsr_release(smat_diag)
551 END DO
552 CALL para_env%sum(occnumabc)
553 END IF
554
555 IF (.NOT. neglect_abc) THEN
556 ! calculate shared electron numbers (ABC)
557 nabc = (natom*(natom - 1)*(natom - 2))/6
558 ALLOCATE (selnabc(nabc, nspin))
559 selnabc = 0.0_dp
560 DO ispin = 1, nspin
561 iabc = 0
562 DO ia = 1, natom
563 DO ib = ia + 1, natom
564 DO ic = ib + 1, natom
565 iabc = iabc + 1
566 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 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 ALLOCATE (raq(natom, nspin))
579 raq = 0.0_dp
580 DO ispin = 1, nspin
581 DO ia = 1, natom
582 raq(ia, ispin) = occnuma(ia, ispin)
583 DO ib = 1, natom
584 raq(ia, ispin) = raq(ia, ispin) - 0.5_dp*selnab(ia, ib, ispin)
585 END DO
586 END DO
587 IF (.NOT. neglect_abc) THEN
588 iabc = 0
589 DO ia = 1, natom
590 DO ib = ia + 1, natom
591 DO ic = ib + 1, natom
592 iabc = iabc + 1
593 raq(ia, ispin) = raq(ia, ispin) + selnabc(iabc, ispin)/3._dp
594 raq(ib, ispin) = raq(ib, ispin) + selnabc(iabc, ispin)/3._dp
595 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 DO ispin = 1, nspin
604 deltaq = (electra(ispin) - sum(raq(1:natom, ispin))) - ua_charge(ispin)
605 IF (unit_nr > 0) THEN
606 WRITE (unit=unit_nr, fmt="(T2,A,T32,A,i2,T55,A,F12.8)") &
607 "Cutoff error on charge", "Spin ", ispin, "error charge =", deltaq
608 END IF
609 END DO
610
611 ! analyze unassigned charge
612 ALLOCATE (uaq(natom, nspin))
613 uaq = 0.0_dp
614 IF (analyze_ua) THEN
615 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
616 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 col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
619 CALL dbcsr_get_info(matrix_s(1, 1)%matrix, row_blk_size=row_blk_sizes)
620 CALL dbcsr_create(amat, name="temp", template=matrix_s(1, 1)%matrix)
621 CALL dbcsr_create(tmat, name="temp", template=mao_coef(1)%matrix)
622 ! replicate diagonal of smm matrix
623 CALL dbcsr_get_block_diag(matrix_smm(1)%matrix, smat_diag)
624 CALL dbcsr_replicate_all(smat_diag)
625
626 ALLOCATE (orb_blk(natom), mao_blk(natom))
627 DO ia = 1, natom
628 orb_blk = row_blk_sizes
629 mao_blk = row_blk_sizes
630 mao_blk(ia) = col_blk_sizes(ia)
631 CALL dbcsr_create(sumat, name="Smat", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
632 row_blk_size=mao_blk, col_blk_size=mao_blk)
633 CALL cp_dbcsr_alloc_block_from_nbl(sumat, sab_orb)
634 CALL dbcsr_create(cholmat, name="Cholesky matrix", dist=dbcsr_dist, &
635 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 row_blk_size=orb_blk, col_blk_size=mao_blk)
638 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 row_blk_size=orb_blk, col_blk_size=mao_blk)
641 ! replicate row and col of smo matrix
642 ALLOCATE (rowblock(natom))
643 DO ib = 1, natom
644 na = mao_blk_sizes(ia)
645 nb = row_blk_sizes(ib)
646 ALLOCATE (rowblock(ib)%mat(na, nb))
647 rowblock(ib)%mat = 0.0_dp
648 CALL dbcsr_get_block_p(matrix=matrix_smo(1)%matrix, row=ia, col=ib, &
649 block=block, found=found)
650 IF (found) rowblock(ib)%mat(1:na, 1:nb) = block(1:na, 1:nb)
651 CALL para_env%sum(rowblock(ib)%mat)
652 END DO
653 !
654 DO ispin = 1, nspin
655 CALL dbcsr_copy(tmat, mao_coef(ispin)%matrix)
656 CALL dbcsr_replicate_all(tmat)
657 CALL dbcsr_iterator_start(dbcsr_iter, matrix_s(1, 1)%matrix)
658 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
659 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, block)
660 CALL dbcsr_get_block_p(matrix=sumat, row=iatom, col=jatom, block=sblk, found=fos)
661 cpassert(fos)
662 CALL dbcsr_get_block_p(matrix=rumat, row=iatom, col=jatom, block=rblku, found=for)
663 cpassert(for)
664 CALL dbcsr_get_block_p(matrix=rumat, row=jatom, col=iatom, block=rblkl, found=for)
665 cpassert(for)
666 CALL dbcsr_get_block_p(matrix=tmat, row=ia, col=ia, block=cmao, found=found)
667 cpassert(found)
668 IF (iatom /= ia .AND. jatom /= ia) THEN
669 ! copy original overlap matrix
670 sblk = block
671 rblku = block
672 rblkl = transpose(block)
673 ELSE IF (iatom /= ia) THEN
674 rblkl = transpose(block)
675 sblk = matmul(transpose(rowblock(iatom)%mat), cmao)
676 rblku = sblk
677 ELSE IF (jatom /= ia) THEN
678 rblku = block
679 sblk = matmul(transpose(cmao), rowblock(jatom)%mat)
680 rblkl = transpose(sblk)
681 ELSE
682 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=block, found=found)
683 cpassert(found)
684 sblk = matmul(transpose(cmao), matmul(block, cmao))
685 rblku = matmul(transpose(rowblock(ia)%mat), cmao)
686 END IF
687 END DO
688 CALL dbcsr_iterator_stop(dbcsr_iter)
689 ! Cholesky decomposition of SUMAT = U'U
690 CALL dbcsr_desymmetrize(sumat, cholmat)
691 CALL cp_dbcsr_cholesky_decompose(cholmat, para_env=para_env, blacs_env=blacs_env)
692 ! T = R*inv(U)
693 ssize = sum(mao_blk)
694 CALL cp_dbcsr_cholesky_restore(rumat, ssize, cholmat, crumat, op="SOLVE", pos="RIGHT", &
695 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 filter_eps=eps_filter)
699 ! Tr(P*A)
700 CALL dbcsr_dot(matrix_p(ispin, 1)%matrix, amat, uaq(ia, ispin))
701 uaq(ia, ispin) = uaq(ia, ispin) - electra(ispin)
702 END DO
703 !
704 CALL dbcsr_release(sumat)
705 CALL dbcsr_release(cholmat)
706 CALL dbcsr_release(rumat)
707 CALL dbcsr_release(crumat)
708 !
709 DO ib = 1, natom
710 DEALLOCATE (rowblock(ib)%mat)
711 END DO
712 DEALLOCATE (rowblock)
713 END DO
714 CALL dbcsr_release(smat_diag)
715 CALL dbcsr_release(amat)
716 CALL dbcsr_release(tmat)
717 DEALLOCATE (orb_blk, mao_blk)
718 END IF
719 !
720 raq(1:natom, 1:nspin) = raq(1:natom, 1:nspin) - uaq(1:natom, 1:nspin)
721 DO ispin = 1, nspin
722 deltaq = electra(ispin) - sum(raq(1:natom, ispin))
723 IF (unit_nr > 0) THEN
724 WRITE (unit=unit_nr, fmt="(T2,A,T32,A,i2,T55,A,F12.8)") &
725 "Charge/Atom redistributed", "Spin ", ispin, "delta charge =", &
726 (deltaq + ua_charge(ispin))/real(natom, kind=dp)
727 END IF
728 END DO
729
730 ! output charges
731 IF (unit_nr > 0) THEN
732 IF (nspin == 1) THEN
733 WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO atomic charges ", "Atom", "Charge"
734 ELSE
735 WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO atomic charges ", "Atom", "Charge", "Spin Charge"
736 END IF
737 DO ispin = 1, nspin
738 deltaq = electra(ispin) - sum(raq(1:natom, ispin))
739 raq(:, ispin) = raq(:, ispin) + deltaq/real(natom, kind=dp)
740 END DO
741 total_charge = 0.0_dp
742 total_spin = 0.0_dp
743 DO iatom = 1, natom
744 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
745 element_symbol=element_symbol, kind_number=ikind)
746 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
747 IF (nspin == 1) THEN
748 WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, zeff - raq(iatom, 1)
749 total_charge = total_charge + (zeff - raq(iatom, 1))
750 ELSE
751 WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
752 zeff - raq(iatom, 1) - raq(iatom, 2), raq(iatom, 1) - raq(iatom, 2)
753 total_charge = total_charge + (zeff - raq(iatom, 1) - raq(iatom, 2))
754 total_spin = total_spin + (raq(iatom, 1) - raq(iatom, 2))
755 END IF
756 END DO
757 IF (nspin == 1) THEN
758 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
759 ELSE
760 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 IF (analyze_ua) THEN
765 ! output unassigned charges
766 IF (unit_nr > 0) THEN
767 IF (nspin == 1) THEN
768 WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO hypervalent charges ", "Atom", "Charge"
769 ELSE
770 WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO hypervalent charges ", "Atom", &
771 "Charge", "Spin Charge"
772 END IF
773 total_charge = 0.0_dp
774 total_spin = 0.0_dp
775 DO iatom = 1, natom
776 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
777 element_symbol=element_symbol)
778 IF (nspin == 1) THEN
779 WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, uaq(iatom, 1)
780 total_charge = total_charge + uaq(iatom, 1)
781 ELSE
782 WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
783 uaq(iatom, 1) + uaq(iatom, 2), uaq(iatom, 1) - uaq(iatom, 2)
784 total_charge = total_charge + uaq(iatom, 1) + uaq(iatom, 2)
785 total_spin = total_spin + uaq(iatom, 1) - uaq(iatom, 2)
786 END IF
787 END DO
788 IF (nspin == 1) THEN
789 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge
790 ELSE
791 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 IF (unit_nr > 0) THEN
798 IF (nspin == 1) THEN
799 WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T78,A)") "Shared electron numbers ", "Atom", "Atom", "SEN"
800 ELSE
801 WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T51,A,T63,A,T71,A)") "Shared electron numbers ", "Atom", "Atom", &
802 "SEN(1)", "SEN(2)", "SEN(total)"
803 END IF
804 DO ia = 1, natom
805 DO ib = ia + 1, natom
806 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
807 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
808 IF (nspin == 1) THEN
809 IF (selnab(ia, ib, 1) > eps_ab) THEN
810 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 IF ((selnab(ia, ib, 1) + selnab(ia, ib, 2)) > eps_ab) THEN
814 WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T45,3F12.6)") ia, esa, ib, esb, &
815 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 IF (.NOT. neglect_abc) THEN
823 ! output shared electron numbers ABC
824 IF (unit_nr > 0) THEN
825 WRITE (unit_nr, "(/,T2,A,T40,A,T49,A,T58,A,T78,A)") "Shared electron numbers ABC", &
826 "Atom", "Atom", "Atom", "SEN"
827 senmax = 0.0_dp
828 iabc = 0
829 DO ia = 1, natom
830 DO ib = ia + 1, natom
831 DO ic = ib + 1, natom
832 iabc = iabc + 1
833 senabc = sum(selnabc(iabc, :))
834 senmax = max(senmax, senabc)
835 IF (senabc > eps_abc) THEN
836 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
837 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
838 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 ia, esa, ib, esb, ic, esc, senabc
841 END IF
842 END DO
843 END DO
844 END DO
845 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Maximum SEN value calculated", senmax
846 END IF
847 END IF
848
849 IF (print_pao) THEN
850 CALL mao_write_pao_restart(mao_coef, qs_env)
851 END IF
852
853 IF (unit_nr > 0) THEN
854 WRITE (unit_nr, '(/,T2,A)') &
855 '!---------------------------END OF MAO ANALYSIS-------------------------------!'
856 END IF
857
858 ! Deallocate temporary arrays
859 DEALLOCATE (occnuma, occnumab, selnab, raq, uaq)
860 IF (.NOT. neglect_abc) THEN
861 DEALLOCATE (occnumabc, selnabc)
862 END IF
863
864 ! Deallocate the neighbor list structure
865 CALL release_neighbor_list_sets(smm_list)
866 CALL release_neighbor_list_sets(smo_list)
867
868 DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
869
870 IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm)
871 IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo)
872 IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q)
873
874 IF (ASSOCIATED(mao_coef)) CALL dbcsr_deallocate_matrix_set(mao_coef)
875 IF (ASSOCIATED(mao_dmat)) CALL dbcsr_deallocate_matrix_set(mao_dmat)
876 IF (ASSOCIATED(mao_smat)) CALL dbcsr_deallocate_matrix_set(mao_smat)
877 IF (ASSOCIATED(mao_qmat)) CALL dbcsr_deallocate_matrix_set(mao_qmat)
878
879 CALL timestop(handle)
880
881 END SUBROUTINE mao_analysis
882
883END MODULE mao_wfn_analysis
for
for(int lxp=0;lxp<=lp;lxp++)
Definition grid_cpu_collint.h:66
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::ehrhardt1985
integer, save, public ehrhardt1985
Definition bibliography.F:36
bibliography::heinzmann1976
integer, save, public heinzmann1976
Definition bibliography.F:36
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_replicate_all
subroutine, public dbcsr_replicate_all(matrix)
...
Definition cp_dbcsr_api.F:1156
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_restore
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
Definition cp_dbcsr_cholesky.F:176
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_reserve_diag_blocks
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
Definition cp_dbcsr_contrib.F:255
cp_dbcsr_contrib::dbcsr_get_block_diag
subroutine, public dbcsr_get_block_diag(matrix, diag)
Copies the diagonal blocks of matrix into diag.
Definition cp_dbcsr_contrib.F:450
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition iterate_matrix.F:473
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
mao_basis
Calculate MAO's and analyze wavefunctions.
Definition mao_basis.F:15
mao_io
Routines for writing PAO restart files from MAO.
Definition mao_io.F:11
mao_io::mao_write_pao_restart
subroutine, public mao_write_pao_restart(mao_coef, qs_env)
Writes restart file.
Definition mao_io.F:54
mao_methods
Calculate MAO's and analyze wavefunctions.
Definition mao_methods.F:15
mao_methods::mao_reference_basis
subroutine, public mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, iunit, print_basis)
Define the MAO reference basis set.
Definition mao_methods.F:684
mao_methods::mao_basis_analysis
subroutine, public mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, qs_kind_set, unit_nr, para_env)
Analyze the MAO basis, projection on angular functions.
Definition mao_methods.F:795
mao_methods::mao_build_q
subroutine, public mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap.
Definition mao_methods.F:890
mao_optimizer
Calculate MAO's and analyze wavefunctions.
Definition mao_optimizer.F:15
mao_optimizer::mao_optimize
subroutine, public mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iolevel, iw)
...
Definition mao_optimizer.F:56
mao_wfn_analysis
Calculate MAO's and analyze wavefunctions.
Definition mao_wfn_analysis.F:15
mao_wfn_analysis::mao_analysis
subroutine, public mao_analysis(qs_env, input_section, unit_nr)
...
Definition mao_wfn_analysis.F:93
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::invmat_symm
subroutine, public invmat_symm(a, potrf, uplo)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:588
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1633
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1096
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition qs_neighbor_lists.F:19
qs_neighbor_lists::setup_neighbor_list
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Definition qs_neighbor_lists.F:1529
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix_simple
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl, lcart)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:575
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62