(git:ed6f26b)
Loading...
Searching...
No Matches
minbas_wfn_analysis.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate localized minimal basis and analyze wavefunctions
10!> \par History
11!> 12.2016 created [JGH]
12!> \author JGH
13! **************************************************************************************************
18 USE bibliography, ONLY: lu2004,&
19 cite_reference
20 USE cell_types, ONLY: cell_type
23 USE cp_dbcsr_api, ONLY: &
27 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
28 dbcsr_type_symmetric
39 USE cp_fm_types, ONLY: cp_fm_create,&
55 USE kinds, ONLY: default_path_length,&
56 dp
60 USE mulliken, ONLY: compute_bond_order,&
66 USE pw_env_types, ONLY: pw_env_get,&
69 USE pw_types, ONLY: pw_c1d_gs_type,&
75 USE qs_ks_types, ONLY: get_ks_env,&
78 USE qs_mo_types, ONLY: allocate_mo_set,&
85#include "./base/base_uses.f90"
86
87 IMPLICIT NONE
88 PRIVATE
89
90 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_wfn_analysis'
91
92 PUBLIC :: minbas_analysis
93
94! **************************************************************************************************
95
96CONTAINS
97
98! **************************************************************************************************
99!> \brief ...
100!> \param qs_env ...
101!> \param input_section ...
102!> \param unit_nr ...
103! **************************************************************************************************
104 SUBROUTINE minbas_analysis(qs_env, input_section, unit_nr)
105 TYPE(qs_environment_type), POINTER :: qs_env
106 TYPE(section_vals_type), POINTER :: input_section
107 INTEGER, INTENT(IN) :: unit_nr
108
109 CHARACTER(len=*), PARAMETER :: routinen = 'minbas_analysis'
110
111 INTEGER :: handle, homo, i, ispin, nao, natom, &
112 nimages, nmao, nmo, nspin
113 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: ecount
114 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
115 LOGICAL :: do_bondorder, explicit, full_ortho, occeq
116 REAL(kind=dp) :: alpha, amax, eps_filter, filter_eps, &
117 trace
118 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: border, fnorm, mcharge, prmao
119 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers
120 TYPE(cp_blacs_env_type), POINTER :: blacs_env
121 TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c
122 TYPE(cp_fm_type) :: fm1, fm2, fm3, fm4
123 TYPE(cp_fm_type), POINTER :: fm_mos
124 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
125 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, pqmat, quambo, sqmat
126 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
127 TYPE(dbcsr_type) :: psmat, sinv, smao, smaox, spmat
128 TYPE(dbcsr_type), POINTER :: smat
129 TYPE(dft_control_type), POINTER :: dft_control
130 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mbas
131 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
132 TYPE(mp_para_env_type), POINTER :: para_env
133 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
134 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
135 TYPE(qs_ks_env_type), POINTER :: ks_env
136 TYPE(section_vals_type), POINTER :: molden_section
137
138 ! only do MINBAS analysis if explicitly requested
139 CALL section_vals_get(input_section, explicit=explicit)
140 IF (.NOT. explicit) RETURN
141
142 ! k-points?
143 CALL get_qs_env(qs_env, dft_control=dft_control)
144 nspin = dft_control%nspins
145 nimages = dft_control%nimages
146 IF (nimages > 1) THEN
147 IF (unit_nr > 0) THEN
148 WRITE (unit=unit_nr, fmt="(T2,A)") &
149 "K-Points: Localized Minimal Basis Analysis not available."
150 END IF
151 END IF
152 IF (nimages > 1) RETURN
153
154 CALL timeset(routinen, handle)
155
156 IF (unit_nr > 0) THEN
157 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
158 WRITE (unit=unit_nr, fmt="(T26,A)") "LOCALIZED MINIMAL BASIS ANALYSIS"
159 WRITE (unit=unit_nr, fmt="(T18,A)") "W.C. Lu et al, J. Chem. Phys. 120, 2629 (2004)"
160 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
161 END IF
162 CALL cite_reference(lu2004)
163
164 ! input options
165 CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter)
166 CALL section_vals_val_get(input_section, "FULL_ORTHOGONALIZATION", l_val=full_ortho)
167 CALL section_vals_val_get(input_section, "BOND_ORDER", l_val=do_bondorder)
168
169 ! generate MAOs and QUAMBOs
170 CALL get_qs_env(qs_env, mos=mos)
171 NULLIFY (quambo, mao_coef)
172 CALL minbas_calculation(qs_env, mos, quambo, mao=mao_coef, iounit=unit_nr, &
173 full_ortho=full_ortho, eps_filter=eps_filter)
174 IF (ASSOCIATED(quambo)) THEN
175 CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=nmo)
176 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
177 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
178 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
179 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
180 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
181 CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
182 nmao = sum(col_blk_sizes)
183
184 NULLIFY (pqmat, sqmat)
185 CALL dbcsr_allocate_matrix_set(sqmat, nspin)
186 CALL dbcsr_allocate_matrix_set(pqmat, nspin)
187 DO ispin = 1, nspin
188 ALLOCATE (sqmat(ispin)%matrix)
189 CALL dbcsr_create(matrix=sqmat(ispin)%matrix, &
190 name="SQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
191 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
192 ALLOCATE (pqmat(ispin)%matrix)
193 CALL dbcsr_create(matrix=pqmat(ispin)%matrix, &
194 name="PQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
195 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
196 END DO
197 DEALLOCATE (row_blk_sizes, col_blk_sizes)
198
199 ! Start wfn analysis
200 IF (unit_nr > 0) THEN
201 WRITE (unit_nr, '(/,T2,A)') 'Localized Minimal Basis Wavefunction Analysis'
202 END IF
203
204 ! localization of basis
205 DO ispin = 1, nspin
206 amax = dbcsr_get_occupation(quambo(ispin)%matrix)
207 IF (unit_nr > 0) THEN
208 WRITE (unit_nr, '(/,T2,A,I2,T69,F10.3,A2,/)') &
209 'Occupation of Basis Function Representation (Spin) ', ispin, amax*100._dp, ' %'
210 END IF
211 END DO
212
213 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
214 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
215 CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, &
216 para_env=para_env, context=blacs_env)
217 CALL cp_fm_create(fm1, fm_struct_a)
218 CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmao, ncol_global=nmo, &
219 para_env=para_env, context=blacs_env)
220 CALL cp_fm_create(fm2, fm_struct_b)
221 CALL cp_fm_create(fm3, fm_struct_b)
222 CALL cp_fm_struct_create(fm_struct_c, nrow_global=nmo, ncol_global=nmo, &
223 para_env=para_env, context=blacs_env)
224 CALL cp_fm_create(fm4, fm_struct_c)
225 ALLOCATE (fnorm(nmo, nspin), ecount(natom, 3, nspin), prmao(natom, nspin))
226 ecount = 0
227 prmao = 0.0_dp
228 DO ispin = 1, nspin
229 CALL dbcsr_create(smao, name="S*QM", template=mao_coef(1)%matrix)
230 smat => matrix_s(1, 1)%matrix
231 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, quambo(ispin)%matrix, 0.0_dp, smao)
232 ! calculate atomic extend of basis
233 CALL pm_extend(quambo(ispin)%matrix, smao, ecount(:, :, ispin))
234 CALL dbcsr_create(sinv, name="QM*S*QM", template=sqmat(ispin)%matrix)
235 CALL dbcsr_multiply("T", "N", 1.0_dp, quambo(ispin)%matrix, smao, 0.0_dp, sqmat(ispin)%matrix)
236 ! atomic MAO projection
237 CALL project_mao(mao_coef(ispin)%matrix, smao, sqmat(ispin)%matrix, prmao(:, ispin))
238 ! invert overlap
239 CALL invert_hotelling(sinv, sqmat(ispin)%matrix, 1.e-6_dp, silent=.true.)
240 CALL dbcsr_create(smaox, name="S*QM*SINV", template=smao)
241 CALL dbcsr_multiply("N", "N", 1.0_dp, smao, sinv, 0.0_dp, smaox)
242 CALL copy_dbcsr_to_fm(smaox, fm1)
243 CALL get_mo_set(mos(ispin), mo_coeff=fm_mos, homo=homo)
244 CALL parallel_gemm("T", "N", nmao, nmo, nao, 1.0_dp, fm1, fm_mos, 0.0_dp, fm2)
245 CALL cp_dbcsr_sm_fm_multiply(sqmat(ispin)%matrix, fm2, fm3, nmo)
246 CALL parallel_gemm("T", "N", nmo, nmo, nmao, 1.0_dp, fm2, fm3, 0.0_dp, fm4)
247 CALL cp_fm_get_diag(fm4, fnorm(1:nmo, ispin))
248 ! fm2 are the projected MOs (in MAO basis); orthogonalize the occupied subspace
249 CALL make_basis_lowdin(vmatrix=fm2, ncol=homo, matrix_s=sqmat(ispin)%matrix)
250 ! pmat
251 CALL get_mo_set(mos(ispin), occupation_numbers=occupation_numbers, maxocc=alpha)
252 occeq = all(occupation_numbers(1:homo) == alpha)
253 CALL dbcsr_copy(pqmat(ispin)%matrix, sqmat(ispin)%matrix)
254 CALL dbcsr_set(pqmat(ispin)%matrix, 0.0_dp)
255 IF (occeq) THEN
256 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, &
257 ncol=homo, alpha=alpha, keep_sparsity=.false.)
258 ELSE
259 CALL cp_fm_to_fm(fm2, fm3)
260 CALL cp_fm_column_scale(fm3, occupation_numbers(1:homo))
261 alpha = 1.0_dp
262 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, &
263 matrix_g=fm3, ncol=homo, alpha=alpha, keep_sparsity=.true.)
264 END IF
265
266 CALL dbcsr_release(smao)
267 CALL dbcsr_release(smaox)
268 CALL dbcsr_release(sinv)
269 END DO
270 ! Basis extension
271 CALL para_env%sum(ecount)
272 IF (unit_nr > 0) THEN
273 IF (nspin == 1) THEN
274 WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 '
275 DO i = 1, natom
276 WRITE (unit_nr, '(T2,I8,T20,I10,T40,I10,T60,I10)') i, ecount(i, 1:3, 1)
277 END DO
278 ELSE
279 WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 '
280 DO i = 1, natom
281 WRITE (unit_nr, '(T2,I8,T20,2I6,T40,2I6,T60,2I6)') &
282 i, ecount(i, 1, 1:2), ecount(i, 2, 1:2), ecount(i, 3, 1:2)
283 END DO
284 END IF
285 END IF
286 ! MAO projection
287 CALL para_env%sum(prmao)
288 IF (unit_nr > 0) THEN
289 DO ispin = 1, nspin
290 WRITE (unit_nr, '(/,T2,A,I2)') 'Projection on same atom MAO orbitals: Spin ', ispin
291 DO i = 1, natom, 2
292 IF (i < natom) THEN
293 WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6,T42,A,I8,T60,A,F10.6)') &
294 " Atom:", i, "Projection:", prmao(i, ispin), " Atom:", i + 1, "Projection:", prmao(i + 1, ispin)
295 ELSE
296 WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6)') " Atom:", i, "Projection:", prmao(i, ispin)
297 END IF
298 END DO
299 END DO
300 END IF
301 ! MO expansion completness
302 DO ispin = 1, nspin
303 CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo)
304 IF (unit_nr > 0) THEN
305 WRITE (unit_nr, '(/,T2,A,I2)') 'MO expansion in Localized Minimal Basis: Spin ', ispin
306 WRITE (unit_nr, '(T64,A)') 'Occupied Orbitals'
307 WRITE (unit_nr, '(8F10.6)') fnorm(1:homo, ispin)
308 WRITE (unit_nr, '(T65,A)') 'Virtual Orbitals'
309 WRITE (unit_nr, '(8F10.6)') fnorm(homo + 1:nmo, ispin)
310 END IF
311 END DO
312 ! Mulliken population
313 IF (unit_nr > 0) THEN
314 WRITE (unit_nr, '(/,T2,A)') 'Mulliken Population Analysis '
315 END IF
316 ALLOCATE (mcharge(natom, nspin))
317 DO ispin = 1, nspin
318 CALL dbcsr_dot(pqmat(ispin)%matrix, sqmat(ispin)%matrix, trace)
319 IF (unit_nr > 0) THEN
320 WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(PS) Spin ', ispin, trace
321 END IF
322 CALL mulliken_charges(pqmat(ispin)%matrix, sqmat(ispin)%matrix, para_env, mcharge(:, ispin))
323 END DO
324 CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Minimal Basis Mulliken Charges", &
325 electronic_charges=mcharge)
326 ! Mayer bond orders
327 IF (do_bondorder) THEN
328 ALLOCATE (border(natom, natom))
329 border = 0.0_dp
330 CALL dbcsr_create(psmat, name="PS", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
331 CALL dbcsr_create(spmat, name="SP", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
332 filter_eps = 1.e-6_dp
333 DO ispin = 1, nspin
334 CALL dbcsr_multiply("N", "N", 1.0_dp, pqmat(ispin)%matrix, sqmat(ispin)%matrix, 0.0_dp, psmat, &
335 filter_eps=filter_eps)
336 CALL dbcsr_multiply("N", "N", 1.0_dp, sqmat(ispin)%matrix, pqmat(ispin)%matrix, 0.0_dp, spmat, &
337 filter_eps=filter_eps)
338 CALL compute_bond_order(psmat, spmat, border)
339 END DO
340 CALL para_env%sum(border)
341 border = border*real(nspin, kind=dp)
342 CALL dbcsr_release(psmat)
343 CALL dbcsr_release(spmat)
344 CALL print_bond_orders(particle_set, unit_nr, border)
345 DEALLOCATE (border)
346 END IF
347
348 ! for printing purposes we now copy the QUAMBOs into MO format
349 ALLOCATE (mbas(nspin))
350 DO ispin = 1, nspin
351 CALL allocate_mo_set(mbas(ispin), nao, nmao, nmao, 0.0_dp, 1.0_dp, 0.0_dp)
352 CALL set_mo_set(mbas(ispin), homo=nmao)
353 ALLOCATE (mbas(ispin)%eigenvalues(nmao))
354 mbas(ispin)%eigenvalues = 0.0_dp
355 ALLOCATE (mbas(ispin)%occupation_numbers(nmao))
356 mbas(ispin)%occupation_numbers = 1.0_dp
357 CALL cp_fm_create(mbas(ispin)%mo_coeff, fm_struct_a)
358 CALL copy_dbcsr_to_fm(quambo(ispin)%matrix, mbas(ispin)%mo_coeff)
359 END DO
360
361 ! Print basis functions: cube files
362 DO ispin = 1, nspin
363 CALL get_mo_set(mbas(ispin), mo_coeff=fm_mos)
364 CALL post_minbas_cubes(qs_env, input_section, fm_mos, ispin)
365 END DO
366 ! Print basis functions: molden format
367 molden_section => section_vals_get_subs_vals(input_section, "MINBAS_MOLDEN")
368 CALL write_mos_molden(mbas, qs_kind_set, particle_set, molden_section)
369 DO ispin = 1, nspin
370 CALL deallocate_mo_set(mbas(ispin))
371 END DO
372 DEALLOCATE (mbas)
373
374 DEALLOCATE (fnorm, ecount, prmao, mcharge)
375 CALL cp_fm_release(fm1)
376 CALL cp_fm_release(fm2)
377 CALL cp_fm_release(fm3)
378 CALL cp_fm_release(fm4)
379 CALL cp_fm_struct_release(fm_struct_a)
380 CALL cp_fm_struct_release(fm_struct_b)
381 CALL cp_fm_struct_release(fm_struct_c)
382
383 ! clean up
386 CALL dbcsr_deallocate_matrix_set(mao_coef)
387 CALL dbcsr_deallocate_matrix_set(quambo)
388
389 END IF
390
391 IF (unit_nr > 0) THEN
392 WRITE (unit_nr, '(/,T2,A)') &
393 '!--------------------------END OF MINBAS ANALYSIS-----------------------------!'
394 END IF
395
396 CALL timestop(handle)
397
398 END SUBROUTINE minbas_analysis
399
400! **************************************************************************************************
401!> \brief ...
402!> \param quambo ...
403!> \param smao ...
404!> \param ecount ...
405! **************************************************************************************************
406 SUBROUTINE pm_extend(quambo, smao, ecount)
407 TYPE(dbcsr_type) :: quambo, smao
408 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: ecount
409
410 INTEGER :: iatom, jatom, n
411 LOGICAL :: found
412 REAL(kind=dp) :: wij
413 REAL(kind=dp), DIMENSION(:, :), POINTER :: qblock, sblock
414 TYPE(dbcsr_iterator_type) :: dbcsr_iter
415
416 CALL dbcsr_iterator_start(dbcsr_iter, quambo)
417 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
418 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock)
419 CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, block=sblock, found=found)
420 IF (found) THEN
421 n = SIZE(qblock, 2)
422 wij = abs(sum(qblock*sblock))/real(n, kind=dp)
423 IF (wij > 0.1_dp) THEN
424 ecount(jatom, 1) = ecount(jatom, 1) + 1
425 ELSEIF (wij > 0.01_dp) THEN
426 ecount(jatom, 2) = ecount(jatom, 2) + 1
427 ELSEIF (wij > 0.001_dp) THEN
428 ecount(jatom, 3) = ecount(jatom, 3) + 1
429 END IF
430 END IF
431 END DO
432 CALL dbcsr_iterator_stop(dbcsr_iter)
433
434 END SUBROUTINE pm_extend
435
436! **************************************************************************************************
437!> \brief ...
438!> \param mao ...
439!> \param smao ...
440!> \param sovl ...
441!> \param prmao ...
442! **************************************************************************************************
443 SUBROUTINE project_mao(mao, smao, sovl, prmao)
444 TYPE(dbcsr_type) :: mao, smao, sovl
445 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: prmao
446
447 INTEGER :: i, iatom, jatom, n
448 LOGICAL :: found
449 REAL(kind=dp) :: wi
450 REAL(kind=dp), DIMENSION(:, :), POINTER :: qblock, sblock, so
451 TYPE(dbcsr_iterator_type) :: dbcsr_iter
452
453 CALL dbcsr_iterator_start(dbcsr_iter, mao)
454 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
455 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock)
456 cpassert(iatom == jatom)
457 CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, block=sblock, found=found)
458 IF (found) THEN
459 CALL dbcsr_get_block_p(matrix=sovl, row=iatom, col=jatom, block=so, found=found)
460 n = SIZE(qblock, 2)
461 DO i = 1, n
462 wi = sum(qblock(:, i)*sblock(:, i))
463 prmao(iatom) = prmao(iatom) + wi/so(i, i)
464 END DO
465 END IF
466 END DO
467 CALL dbcsr_iterator_stop(dbcsr_iter)
468
469 END SUBROUTINE project_mao
470
471! **************************************************************************************************
472!> \brief Computes and prints the Cube Files for the minimal basis set
473!> \param qs_env ...
474!> \param print_section ...
475!> \param minbas_coeff ...
476!> \param ispin ...
477! **************************************************************************************************
478 SUBROUTINE post_minbas_cubes(qs_env, print_section, minbas_coeff, ispin)
479 TYPE(qs_environment_type), POINTER :: qs_env
480 TYPE(section_vals_type), POINTER :: print_section
481 TYPE(cp_fm_type), INTENT(IN) :: minbas_coeff
482 INTEGER, INTENT(IN) :: ispin
483
484 CHARACTER(LEN=default_path_length) :: filename, title
485 INTEGER :: i, i_rep, ivec, iw, j, n_rep, natom
486 INTEGER, DIMENSION(:), POINTER :: blk_sizes, first_bas, ilist, stride
487 LOGICAL :: explicit, mpi_io
488 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
489 TYPE(cell_type), POINTER :: cell
490 TYPE(cp_logger_type), POINTER :: logger
491 TYPE(dft_control_type), POINTER :: dft_control
492 TYPE(particle_list_type), POINTER :: particles
493 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
494 TYPE(pw_c1d_gs_type) :: wf_g
495 TYPE(pw_env_type), POINTER :: pw_env
496 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
497 TYPE(pw_r3d_rs_type) :: wf_r
498 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
499 TYPE(qs_subsys_type), POINTER :: subsys
500 TYPE(section_vals_type), POINTER :: minbas_section
501
502 minbas_section => section_vals_get_subs_vals(print_section, "MINBAS_CUBE")
503 CALL section_vals_get(minbas_section, explicit=explicit)
504 IF (.NOT. explicit) RETURN
505
506 logger => cp_get_default_logger()
507 stride => section_get_ivals(print_section, "MINBAS_CUBE%STRIDE")
508
509 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
510 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
511 CALL qs_subsys_get(subsys, particles=particles)
512
513 CALL get_qs_env(qs_env=qs_env, natom=natom)
514 ALLOCATE (blk_sizes(natom), first_bas(0:natom))
515 CALL get_particle_set(particle_set, qs_kind_set, nmao=blk_sizes)
516 first_bas(0) = 0
517 DO i = 1, natom
518 first_bas(i) = first_bas(i - 1) + blk_sizes(i)
519 END DO
520
521 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
522 CALL auxbas_pw_pool%create_pw(wf_r)
523 CALL auxbas_pw_pool%create_pw(wf_g)
524
525 ! loop over list of atoms
526 CALL section_vals_val_get(minbas_section, "ATOM_LIST", n_rep_val=n_rep)
527 IF (n_rep == 0) THEN
528 DO i = 1, natom
529 DO ivec = first_bas(i - 1) + 1, first_bas(i)
530 WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin
531 WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", i, " spin ", ispin
532 mpi_io = .true.
533 iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", &
534 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
535 mpi_io=mpi_io)
536 CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
537 cell, dft_control, particle_set, pw_env)
538 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
539 CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io)
540 END DO
541 END DO
542 ELSE
543 DO i_rep = 1, n_rep
544 CALL section_vals_val_get(minbas_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=ilist)
545 DO i = 1, SIZE(ilist, 1)
546 j = ilist(i)
547 DO ivec = first_bas(j - 1) + 1, first_bas(j)
548 WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin
549 WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", j, " spin ", ispin
550 mpi_io = .true.
551 iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", &
552 middle_name=trim(filename), file_position="REWIND", log_filename=.false., &
553 mpi_io=mpi_io)
554 CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
555 cell, dft_control, particle_set, pw_env)
556 CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
557 CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io)
558 END DO
559 END DO
560 END DO
561 END IF
562 DEALLOCATE (blk_sizes, first_bas)
563 CALL auxbas_pw_pool%give_back_pw(wf_r)
564 CALL auxbas_pw_pool%give_back_pw(wf_g)
565
566 END SUBROUTINE post_minbas_cubes
567
568END MODULE minbas_wfn_analysis
simple routine to print charges for all atomic charge methods (currently mulliken,...
subroutine, public print_bond_orders(particle_set, scr, bond_orders)
Print Mayer bond orders.
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public lu2004
Handles all functions related to the CELL.
Definition cell_types.F:15
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
basic linear algebra operations for full matrices
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
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
integer, parameter, public default_path_length
Definition kinds.F:58
Interface to the message passing library MPI.
Calculate localized minimal basis.
subroutine, public minbas_calculation(qs_env, mos, quambo, mao, iounit, full_ortho, eps_filter)
...
Calculate localized minimal basis and analyze wavefunctions.
subroutine, public minbas_analysis(qs_env, input_section, unit_nr)
...
Functions handling the MOLDEN format. Split from mode_selective.
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section)
Write out the MOs in molden format for visualisation.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
subroutine, public compute_bond_order(psmat, spmat, bond_order)
compute Mayer bond orders for a single spin channel for complete result sum up over all spins and mul...
Definition mulliken.F:657
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
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.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
collects routines that perform operations directly related to MOs
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...