(git:0de0cc2)
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 ! **************************************************************************************************
17  USE basis_set_types, ONLY: gto_basis_set_p_type
18  USE bibliography, ONLY: ehrhardt1985,&
20  cite_reference
21  USE cp_blacs_env, ONLY: cp_blacs_env_type
22  USE cp_control_types, ONLY: dft_control_type
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
35  section_vals_type,&
38  USE kinds, ONLY: dp
39  USE kpoint_types, ONLY: kpoint_type
40  USE mao_methods, ONLY: mao_basis_analysis,&
41  mao_build_q,&
43  USE mao_optimizer, ONLY: mao_optimize
44  USE mathlib, ONLY: invmat_symm
45  USE message_passing, ONLY: mp_para_env_type
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
57  neighbor_list_iterator_p_type,&
59  neighbor_list_set_p_type,&
63  USE qs_rho_types, ONLY: qs_rho_get,&
64  qs_rho_type
65 #include "./base/base_uses.f90"
66 
67  IMPLICIT NONE
68  PRIVATE
69 
70  TYPE block_type
71  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: mat
72  END TYPE block_type
73 
74  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_wfn_analysis'
75 
76  PUBLIC :: mao_analysis
77 
78 ! **************************************************************************************************
79 
80 CONTAINS
81 
82 ! **************************************************************************************************
83 !> \brief ...
84 !> \param qs_env ...
85 !> \param input_section ...
86 !> \param unit_nr ...
87 ! **************************************************************************************************
88  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 
872 END MODULE mao_wfn_analysis
for(int lxp=0;lxp<=lp;lxp++)
Define the atomic kind types and their sub types.
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.
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public ehrhardt1985
Definition: bibliography.F:43
integer, save, public heinzmann1976
Definition: bibliography.F:43
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to (sca)lapack for the Cholesky based procedures.
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,...
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
DBCSR operations in CP2K.
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Routines useful for iterative matrix calculations.
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 ...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
Calculate MAO's and analyze wavefunctions.
Definition: mao_basis.F:15
Calculate MAO's and analyze wavefunctions.
Definition: mao_methods.F:15
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
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
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
Calculate MAO's and analyze wavefunctions.
Definition: mao_optimizer.F:15
subroutine, public mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iolevel, iw)
...
Definition: mao_optimizer.F:55
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_analysis(qs_env, input_section, unit_nr)
...
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:579
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
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.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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.
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
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Generate the atomic neighbor lists.
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
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
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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