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