(git:c51e757)
Loading...
Searching...
No Matches
population_analyses.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Provide various population analyses and print the requested output
9!> information
10!>
11!> \author Matthias Krack (MK)
12!> \date 09.07.2010
13!> \version 1.0
14! **************************************************************************************************
15
16MODULE population_analyses
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind,&
19 get_atomic_kind_set
20 USE basis_set_types, ONLY: get_gto_basis_set,&
21 gto_basis_set_type
22 USE cp_blacs_env, ONLY: cp_blacs_env_type
23 USE cp_dbcsr_api, ONLY: &
24 dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
25 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
26 dbcsr_p_type, dbcsr_set, dbcsr_type
27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 cp_dbcsr_sm_fm_multiply
29 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,&
30 write_fm_with_basis_info
31 USE cp_fm_diag, ONLY: cp_fm_power
32 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 cp_fm_struct_release,&
34 cp_fm_struct_type
35 USE cp_fm_types, ONLY: cp_fm_create,&
36 cp_fm_get_diag,&
37 cp_fm_release,&
38 cp_fm_type
39 USE cp_result_methods, ONLY: cp_results_erase,&
40 put_results
41 USE cp_result_types, ONLY: cp_result_type
42 USE kinds, ONLY: default_string_length,&
43 dp
44 USE kpoint_methods, ONLY: lowdin_kp_trans
45 USE kpoint_types, ONLY: kpoint_type
46 USE machine, ONLY: m_flush
47 USE message_passing, ONLY: mp_para_env_type
48 USE orbital_pointers, ONLY: nso
49 USE parallel_gemm_api, ONLY: parallel_gemm
50 USE particle_methods, ONLY: get_particle_set
51 USE particle_types, ONLY: particle_type
52 USE qs_environment_types, ONLY: get_qs_env,&
53 qs_environment_type
54 USE qs_kind_types, ONLY: get_qs_kind,&
55 get_qs_kind_set,&
56 qs_kind_type
57 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
58 USE qs_rho_types, ONLY: qs_rho_get,&
59 qs_rho_type
60 USE qs_scf_diagonalization, ONLY: diag_kp_smat
61 USE scf_control_types, ONLY: scf_control_type
62#include "./base/base_uses.f90"
63
64 IMPLICIT NONE
65
66 PRIVATE
67
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'population_analyses'
69
70 PUBLIC :: lowdin_population_analysis, &
71 mulliken_population_analysis
72
73CONTAINS
74
75! **************************************************************************************************
76!> \brief Perform a Lowdin population analysis based on a symmetric
77!> orthogonalisation of the density matrix using S^(1/2)
78!>
79!> \param qs_env ...
80!> \param output_unit ...
81!> \param print_level ...
82!> \date 06.07.2010
83!> \author Matthias Krack (MK)
84!> \version 1.0
85! **************************************************************************************************
86 SUBROUTINE lowdin_population_analysis(qs_env, output_unit, print_level)
87
88 TYPE(qs_environment_type), POINTER :: qs_env
89 INTEGER, INTENT(IN) :: output_unit, print_level
90
91 CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin_population_analysis'
92
93 CHARACTER(LEN=default_string_length) :: headline
94 INTEGER :: handle, i, ispin, ndep, nimg, nsgf, nspin
95 LOGICAL :: do_kpoints, print_gop
96 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
97 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
98 TYPE(cp_blacs_env_type), POINTER :: blacs_env
99 TYPE(cp_fm_struct_type), POINTER :: fmstruct
100 TYPE(cp_fm_type) :: fm_s_half, fm_work1, fm_work2
101 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
102 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
103 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_p, matrixkp_s
104 TYPE(dbcsr_type), POINTER :: sm_p, sm_s
105 TYPE(kpoint_type), POINTER :: kpoints
106 TYPE(mp_para_env_type), POINTER :: para_env
107 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108 POINTER :: sab_nl
109 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
111 TYPE(qs_rho_type), POINTER :: rho
112 TYPE(scf_control_type), POINTER :: scf_control
113
114 CALL timeset(routinen, handle)
115
116 NULLIFY (sm_p, sm_s)
117
118 CALL get_qs_env(qs_env=qs_env, &
119 atomic_kind_set=atomic_kind_set, &
120 qs_kind_set=qs_kind_set, &
121 do_kpoints=do_kpoints, &
122 matrix_s_kp=matrixkp_s, &
123 particle_set=particle_set, &
124 rho=rho, &
125 scf_control=scf_control, &
126 para_env=para_env, &
127 blacs_env=blacs_env)
128
129 CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
130 nspin = SIZE(matrixkp_p, 1)
131 nimg = SIZE(matrixkp_p, 2)
132
133 ! Get the total number of contracted spherical Gaussian basis functions
134 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
135 ! Provide an array to store the orbital populations for each spin
136 ALLOCATE (orbpop(nsgf, nspin))
137 orbpop(:, :) = 0.0_dp
138
139 ! Write headline
140 IF (output_unit > 0) THEN
141 WRITE (unit=output_unit, fmt="(/,/,T2,A)") "LOWDIN POPULATION ANALYSIS"
142 END IF
143
144 IF (do_kpoints) THEN
145
146 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints, sab_orb=sab_nl)
147
148 ! work matrices
149 CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, &
150 nrow_global=nsgf, ncol_global=nsgf)
151 ALLOCATE (fmwork(4))
152 DO i = 1, 4
153 CALL cp_fm_create(fmwork(i), matrix_struct=fmstruct, name="work")
154 END DO
155 CALL cp_fm_struct_release(fmstruct)
156
157 ! S^1/2
158 CALL diag_kp_smat(matrixkp_s, kpoints, fmwork)
159 ! Lowdin P Matrix Transform
160 CALL lowdin_kp_trans(kpoints, orbpop)
161
162 CALL cp_fm_release(fmwork)
163
164 ELSE
165
166 sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
167 matrix_p => matrixkp_p(:, 1)
168
169 ! Provide full size work matrices
170 CALL cp_fm_struct_create(fmstruct=fmstruct, &
171 para_env=para_env, &
172 context=blacs_env, &
173 nrow_global=nsgf, &
174 ncol_global=nsgf)
175 CALL cp_fm_create(matrix=fm_s_half, &
176 matrix_struct=fmstruct, &
177 name="S^(1/2) MATRIX")
178 CALL cp_fm_create(matrix=fm_work1, &
179 matrix_struct=fmstruct, &
180 name="FULL WORK MATRIX 1")
181 headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
182 CALL cp_fm_create(matrix=fm_work2, &
183 matrix_struct=fmstruct, &
184 name=trim(headline))
185 CALL cp_fm_struct_release(fmstruct=fmstruct)
186
187 ! Build full S^(1/2) matrix (computationally expensive)
188 CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
189 CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
190 IF (ndep /= 0) THEN
191 CALL cp_warn(__location__, &
192 "Overlap matrix exhibits linear dependencies. At least some "// &
193 "eigenvalues have been quenched.")
194 END IF
195
196 ! Build Lowdin population matrix for each spin
197 DO ispin = 1, nspin
198 sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format
199 ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
200 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
201 CALL parallel_gemm(transa="N", &
202 transb="N", &
203 m=nsgf, &
204 n=nsgf, &
205 k=nsgf, &
206 alpha=1.0_dp, &
207 matrix_a=fm_s_half, &
208 matrix_b=fm_work1, &
209 beta=0.0_dp, &
210 matrix_c=fm_work2)
211 IF (print_level > 2) THEN
212 ! Write the full Lowdin population matrix
213 IF (nspin > 1) THEN
214 IF (ispin == 1) THEN
215 fm_work2%name = trim(headline)//" FOR ALPHA SPIN"
216 ELSE
217 fm_work2%name = trim(headline)//" FOR BETA SPIN"
218 END IF
219 END IF
220 CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
221 output_unit=output_unit)
222 END IF
223 CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
224 END DO ! next spin ispin
225
226 ! Release local working storage
227 CALL cp_fm_release(matrix=fm_s_half)
228 CALL cp_fm_release(matrix=fm_work1)
229 CALL cp_fm_release(matrix=fm_work2)
230
231 END IF
232
233 ! Write atomic populations and charges
234 IF (output_unit > 0) THEN
235 print_gop = (print_level > 1) ! Print also orbital populations
236 CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
237 END IF
238
239 DEALLOCATE (orbpop)
240
241 CALL timestop(handle)
242
243 END SUBROUTINE lowdin_population_analysis
244
245! **************************************************************************************************
246!> \brief Perform a Mulliken population analysis
247!>
248!> \param qs_env ...
249!> \param output_unit ...
250!> \param print_level ...
251!> \date 10.07.2010
252!> \author Matthias Krack (MK)
253!> \version 1.0
254! **************************************************************************************************
255 SUBROUTINE mulliken_population_analysis(qs_env, output_unit, print_level)
256
257 TYPE(qs_environment_type), POINTER :: qs_env
258 INTEGER, INTENT(IN) :: output_unit, print_level
259
260 CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken_population_analysis'
261
262 CHARACTER(LEN=default_string_length) :: headline
263 INTEGER :: handle, iatom, ic, isgf, ispin, jatom, &
264 jsgf, natom, nsgf, nspin, sgfa, sgfb
265 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
266 LOGICAL :: found, print_gop
267 REAL(kind=dp) :: ps
268 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
269 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, ps_block, s_block
270 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
271 TYPE(dbcsr_iterator_type) :: iter
272 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
273 TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s
274 TYPE(mp_para_env_type), POINTER :: para_env
275 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
276 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
277 TYPE(qs_rho_type), POINTER :: rho
278
279 CALL timeset(routinen, handle)
280
281 NULLIFY (sm_p, sm_ps, sm_s)
282 NULLIFY (p_block, s_block, ps_block)
283
284 CALL get_qs_env(qs_env=qs_env, &
285 atomic_kind_set=atomic_kind_set, &
286 qs_kind_set=qs_kind_set, &
287 matrix_s_kp=matrix_s, &
288 particle_set=particle_set, &
289 rho=rho, &
290 para_env=para_env)
291
292 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
293 nspin = SIZE(matrix_p, 1)
294
295 ! Get the total number of contracted spherical Gaussian basis functions
296 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
297 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
298 ALLOCATE (first_sgf_atom(natom))
299 first_sgf_atom(:) = 0
300
301 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
302
303 ! Provide an array to store the orbital populations for each spin
304 ALLOCATE (orbpop(nsgf, nspin))
305 orbpop(:, :) = 0.0_dp
306
307 ! Write headline
308 IF (output_unit > 0) THEN
309 WRITE (unit=output_unit, fmt="(/,T2,A)") &
310 '!-----------------------------------------------------------------------------!'
311 WRITE (unit=output_unit, fmt="(T22,A)") "Mulliken Population Analysis"
312 END IF
313
314 ! Create a DBCSR work matrix, if needed
315 IF (print_level > 2) THEN
316 sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
317 ALLOCATE (sm_ps)
318 headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
319 IF (nspin > 1) THEN
320 IF (ispin == 1) THEN
321 headline = trim(headline)//" For Alpha Spin"
322 ELSE
323 headline = trim(headline)//" For Beta Spin"
324 END IF
325 END IF
326 CALL dbcsr_copy(matrix_b=sm_ps, matrix_a=sm_s, name=trim(headline))
327 END IF
328
329 ! Build Mulliken population matrix for each spin
330 DO ispin = 1, nspin
331 DO ic = 1, SIZE(matrix_s, 2)
332 IF (print_level > 2) THEN
333 CALL dbcsr_set(sm_ps, 0.0_dp)
334 END IF
335 sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
336 sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
337 ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
338 ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
339 CALL dbcsr_iterator_start(iter, sm_s)
340 DO WHILE (dbcsr_iterator_blocks_left(iter))
341 CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block)
342 IF (.NOT. (ASSOCIATED(s_block))) cycle
343 CALL dbcsr_get_block_p(matrix=sm_p, &
344 row=iatom, &
345 col=jatom, &
346 block=p_block, &
347 found=found)
348 IF (print_level > 2) THEN
349 CALL dbcsr_get_block_p(matrix=sm_ps, &
350 row=iatom, &
351 col=jatom, &
352 block=ps_block, &
353 found=found)
354 cpassert(ASSOCIATED(ps_block))
355 END IF
356
357 sgfb = first_sgf_atom(jatom)
358 DO jsgf = 1, SIZE(s_block, 2)
359 DO isgf = 1, SIZE(s_block, 1)
360 ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
361 IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
362 orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
363 END DO
364 sgfb = sgfb + 1
365 END DO
366 IF (iatom /= jatom) THEN
367 sgfa = first_sgf_atom(iatom)
368 DO isgf = 1, SIZE(s_block, 1)
369 DO jsgf = 1, SIZE(s_block, 2)
370 ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
371 orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
372 END DO
373 sgfa = sgfa + 1
374 END DO
375 END IF
376 END DO
377 CALL dbcsr_iterator_stop(iter)
378 END DO
379
380 IF (print_level > 2) THEN
381 ! Write the full Mulliken net AO and overlap population matrix
382 CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, output_unit=output_unit)
383 END IF
384 END DO
385
386 CALL para_env%sum(orbpop)
387
388 ! Write atomic populations and charges
389 IF (output_unit > 0) THEN
390 print_gop = (print_level > 1) ! Print also orbital populations
391 CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
392 END IF
393
394 ! Save the Mulliken charges to results
395 CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
396
397 ! Release local working storage
398 IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
399 DEALLOCATE (orbpop)
400 DEALLOCATE (first_sgf_atom)
401
402 IF (output_unit > 0) THEN
403 WRITE (unit=output_unit, fmt="(T2,A)") &
404 '!-----------------------------------------------------------------------------!'
405 END IF
406
407 CALL timestop(handle)
408
409 END SUBROUTINE mulliken_population_analysis
410
411! **************************************************************************************************
412!> \brief Save the Mulliken atomic orbital populations and charges in results
413!>
414!> \param orbpop ...
415!> \param atomic_kind_set ...
416!> \param qs_kind_set ...
417!> \param particle_set ...
418!> \param qs_env ...
419!> \par History
420!> 27.05.2022 BT
421!> 16.07.2025 RK
422!> \author Bo Thomsen (BT)
423!> Rangsiman Ketkaew (RK)
424!> \version 1.0
425! **************************************************************************************************
426 SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
427 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: orbpop
428 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
429 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
430 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
431 TYPE(qs_environment_type), POINTER :: qs_env
432
433 CHARACTER(LEN=default_string_length) :: description
434 INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
435 iso, l, natom, nset, nsgf, nspin
436 INTEGER, DIMENSION(:), POINTER :: nshell
437 INTEGER, DIMENSION(:, :), POINTER :: lshell
438 REAL(kind=dp) :: zeff
439 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: all_sumorbpop, charges_save
440 REAL(kind=dp), DIMENSION(3) :: sumorbpop
441 TYPE(cp_result_type), POINTER :: results
442 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
443
444 nspin = SIZE(orbpop, 2)
445
446 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
447 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
448 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
449 NULLIFY (results)
450 CALL get_qs_env(qs_env, results=results)
451 ALLOCATE (all_sumorbpop(natom))
452 ALLOCATE (charges_save(natom))
453
454 iao = 1
455 DO iatom = 1, natom
456 sumorbpop(:) = 0.0_dp
457 NULLIFY (orb_basis_set)
458 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
459 kind_number=ikind)
460 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
461 IF (ASSOCIATED(orb_basis_set)) THEN
462 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
463 nset=nset, &
464 nshell=nshell, &
465 l=lshell)
466 isgf = 1
467 DO iset = 1, nset
468 DO ishell = 1, nshell(iset)
469 l = lshell(ishell, iset)
470 DO iso = 1, nso(l)
471 IF (nspin == 1) THEN
472 sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
473 ELSE
474 sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
475 sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
476 END IF
477 isgf = isgf + 1
478 iao = iao + 1
479 END DO
480 END DO
481 END DO
482 IF (nspin == 1) THEN
483 charges_save(iatom) = zeff - sumorbpop(1)
484 all_sumorbpop(iatom) = sumorbpop(1)
485 ELSE
486 charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
487 all_sumorbpop(iatom) = sumorbpop(1) + sumorbpop(2)
488 END IF
489 END IF ! atom has an orbital basis
490 END DO ! next atom iatom
491
492 ! Store atomic orbital populations in results
493 description = "[MULLIKEN-ORBPOP]"
494 CALL cp_results_erase(results=results, description=description)
495 CALL put_results(results=results, description=description, &
496 values=orbpop)
497
498 ! Store sum orbital population in results
499 description = "[MULLIKEN-SUMORBPOP]"
500 CALL cp_results_erase(results=results, description=description)
501 CALL put_results(results=results, description=description, &
502 values=all_sumorbpop)
503
504 ! Store charges in results
505 description = "[MULLIKEN-CHARGES]"
506 CALL cp_results_erase(results=results, description=description)
507 CALL put_results(results=results, description=description, &
508 values=charges_save)
509
510 DEALLOCATE (all_sumorbpop)
511 DEALLOCATE (charges_save)
512
513 END SUBROUTINE save_mulliken_charges
514
515! **************************************************************************************************
516!> \brief Write atomic orbital populations and net atomic charges
517!>
518!> \param orbpop ...
519!> \param atomic_kind_set ...
520!> \param qs_kind_set ...
521!> \param particle_set ...
522!> \param output_unit ...
523!> \param print_orbital_contributions ...
524!> \date 07.07.2010
525!> \author Matthias Krack (MK)
526!> \version 1.0
527! **************************************************************************************************
528 SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
529 print_orbital_contributions)
530
531 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: orbpop
532 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
533 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
534 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
535 INTEGER, INTENT(IN) :: output_unit
536 LOGICAL, INTENT(IN) :: print_orbital_contributions
537
538 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_orbpop'
539
540 CHARACTER(LEN=2) :: element_symbol
541 CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
542 INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
543 ishell, iso, l, natom, nset, nsgf, &
544 nspin
545 INTEGER, DIMENSION(:), POINTER :: nshell
546 INTEGER, DIMENSION(:, :), POINTER :: lshell
547 REAL(kind=dp) :: zeff
548 REAL(kind=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
549 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
550
551 CALL timeset(routinen, handle)
552
553 nspin = SIZE(orbpop, 2)
554
555 CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
556 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
557
558 ! Select and write headline
559 IF (nspin == 1) THEN
560 IF (print_orbital_contributions) THEN
561 WRITE (unit=output_unit, fmt="(/,T2,A)") &
562 "# Orbital AO symbol Orbital population Net charge"
563 ELSE
564 WRITE (unit=output_unit, fmt="(/,T2,A)") &
565 "# Atom Element Kind Atomic population Net charge"
566 END IF
567 ELSE
568 IF (print_orbital_contributions) THEN
569 WRITE (unit=output_unit, fmt="(/,T2,A)") &
570 "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
571 ELSE
572 WRITE (unit=output_unit, fmt="(/,T2,A)") &
573 "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
574 END IF
575 END IF
576
577 totsumorbpop(:) = 0.0_dp
578
579 iao = 1
580 DO iatom = 1, natom
581 sumorbpop(:) = 0.0_dp
582 NULLIFY (orb_basis_set)
583 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
584 element_symbol=element_symbol, &
585 kind_number=ikind)
586 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
587 IF (ASSOCIATED(orb_basis_set)) THEN
588 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
589 nset=nset, &
590 nshell=nshell, &
591 l=lshell, &
592 sgf_symbol=sgf_symbol)
593 isgf = 1
594 DO iset = 1, nset
595 DO ishell = 1, nshell(iset)
596 l = lshell(ishell, iset)
597 DO iso = 1, nso(l)
598 IF (nspin == 1) THEN
599 sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
600 IF (print_orbital_contributions) THEN
601 IF (isgf == 1) WRITE (unit=output_unit, fmt="(A)") ""
602 WRITE (unit=output_unit, &
603 fmt="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
604 iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
605 END IF
606 ELSE
607 sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
608 sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
609 IF (print_orbital_contributions) THEN
610 IF (isgf == 1) WRITE (unit=output_unit, fmt="(A)") ""
611 WRITE (unit=output_unit, &
612 fmt="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
613 iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
614 orbpop(iao, 1) - orbpop(iao, 2)
615 END IF
616 END IF
617 isgf = isgf + 1
618 iao = iao + 1
619 END DO
620 END DO
621 END DO
622 IF (nspin == 1) THEN
623 totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
624 totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
625 WRITE (unit=output_unit, &
626 fmt="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
627 iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
628 ELSE
629 totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
630 totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
631 WRITE (unit=output_unit, &
632 fmt="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
633 iatom, element_symbol, ikind, sumorbpop(1:2), &
634 zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
635 END IF
636 END IF ! atom has an orbital basis
637 END DO ! next atom iatom
638
639 ! Write total sums
640 IF (print_orbital_contributions) WRITE (unit=output_unit, fmt="(A)") ""
641 IF (nspin == 1) THEN
642 WRITE (unit=output_unit, &
643 fmt="(T2,A,T42,F12.6,T68,F12.6,/)") &
644 "# Total charge", totsumorbpop(1), totsumorbpop(3)
645 ELSE
646 WRITE (unit=output_unit, &
647 fmt="(T2,A,T28,4(1X,F12.6),/)") &
648 "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
649 totsumorbpop(1) - totsumorbpop(2)
650 END IF
651
652 IF (output_unit > 0) CALL m_flush(output_unit)
653
654 CALL timestop(handle)
655
656 END SUBROUTINE write_orbpop
657
658END MODULE population_analyses
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_result_methods::put_results
Definition cp_result_methods.F:39
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
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
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::write_fm_with_basis_info
subroutine, public write_fm_with_basis_info(blacs_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, output_unit, omit_headers)
Print a spherical matrix of blacs type.
Definition cp_dbcsr_output.F:79
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1127
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
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
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition cp_result_methods.F:365
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::lowdin_kp_trans
subroutine, public lowdin_kp_trans(kpoint, pmat_diag)
Calculate Lowdin transformation of density matrix S^1/2 P S^1/2 Integrate diagonal elements over k-po...
Definition kpoint_methods.F:1610
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
population_analyses
Provide various population analyses and print the requested output information.
Definition population_analyses.F:16
population_analyses::lowdin_population_analysis
subroutine, public lowdin_population_analysis(qs_env, output_unit, print_level)
Perform a Lowdin population analysis based on a symmetric orthogonalisation of the density matrix usi...
Definition population_analyses.F:87
population_analyses::mulliken_population_analysis
subroutine, public mulliken_population_analysis(qs_env, output_unit, print_level)
Perform a Mulliken population analysis.
Definition population_analyses.F:256
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:486
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:979
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
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
qs_scf_diagonalization
Different diagonalization schemes that can be used for the iterative solution of the eigenvalue probl...
Definition qs_scf_diagonalization.F:15
qs_scf_diagonalization::diag_kp_smat
subroutine, public diag_kp_smat(matrix_s, kpoints, fmwork)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:2135
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
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_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:185
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
scf_control_types::scf_control_type
Definition scf_control_types.F:126