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