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