(git:34ef472)
population_analyses.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 !> \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 
16 MODULE 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_operations, ONLY: copy_dbcsr_to_fm,&
24  cp_dbcsr_sm_fm_multiply
25  USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,&
26  write_fm_with_basis_info
27  USE cp_fm_diag, ONLY: cp_fm_power
28  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
29  cp_fm_struct_release,&
30  cp_fm_struct_type
31  USE cp_fm_types, ONLY: cp_fm_create,&
32  cp_fm_get_diag,&
33  cp_fm_release,&
34  cp_fm_type
35  USE cp_result_methods, ONLY: cp_results_erase,&
36  put_results
37  USE cp_result_types, ONLY: cp_result_type
38  USE dbcsr_api, ONLY: &
39  dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
40  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
41  dbcsr_p_type, dbcsr_set, dbcsr_setname, dbcsr_type
42  USE kinds, ONLY: default_string_length,&
43  dp
44  USE machine, ONLY: m_flush
45  USE message_passing, ONLY: mp_para_env_type
46  USE orbital_pointers, ONLY: nso
47  USE parallel_gemm_api, ONLY: parallel_gemm
48  USE particle_methods, ONLY: get_particle_set
49  USE particle_types, ONLY: particle_type
50  USE qs_environment_types, ONLY: get_qs_env,&
51  qs_environment_type
52  USE qs_kind_types, ONLY: get_qs_kind,&
53  get_qs_kind_set,&
54  qs_kind_type
55  USE qs_rho_types, ONLY: qs_rho_get,&
56  qs_rho_type
57  USE scf_control_types, ONLY: scf_control_type
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61 
62  PRIVATE
63 
64  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'population_analyses'
65 
66  PUBLIC :: lowdin_population_analysis, &
67  mulliken_population_analysis
68 
69 CONTAINS
70 
71 ! **************************************************************************************************
72 !> \brief Perform a Lowdin population analysis based on a symmetric
73 !> orthogonalisation of the density matrix using S^(1/2)
74 !>
75 !> \param qs_env ...
76 !> \param output_unit ...
77 !> \param print_level ...
78 !> \date 06.07.2010
79 !> \author Matthias Krack (MK)
80 !> \version 1.0
81 ! **************************************************************************************************
82  SUBROUTINE lowdin_population_analysis(qs_env, output_unit, print_level)
83 
84  TYPE(qs_environment_type), POINTER :: qs_env
85  INTEGER, INTENT(IN) :: output_unit, print_level
86 
87  CHARACTER(LEN=*), PARAMETER :: routinen = 'lowdin_population_analysis'
88 
89  CHARACTER(LEN=default_string_length) :: headline
90  INTEGER :: handle, ispin, ndep, nsgf, nspin
91  LOGICAL :: print_gop
92  REAL(kind=dp), DIMENSION(:, :), POINTER :: orbpop
93  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
94  TYPE(cp_blacs_env_type), POINTER :: blacs_env
95  TYPE(cp_fm_struct_type), POINTER :: fmstruct
96  TYPE(cp_fm_type) :: fm_s_half, fm_work1, fm_work2
97  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
98  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_p, matrixkp_s
99  TYPE(dbcsr_type), POINTER :: sm_p, sm_s
100  TYPE(mp_para_env_type), POINTER :: para_env
101  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
102  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
103  TYPE(qs_rho_type), POINTER :: rho
104  TYPE(scf_control_type), POINTER :: scf_control
105 
106  CALL timeset(routinen, handle)
107 
108  NULLIFY (atomic_kind_set)
109  NULLIFY (qs_kind_set)
110  NULLIFY (fmstruct)
111  NULLIFY (matrix_p)
112  NULLIFY (matrixkp_p)
113  NULLIFY (matrixkp_s)
114  NULLIFY (orbpop)
115  NULLIFY (particle_set)
116  NULLIFY (rho)
117  NULLIFY (scf_control)
118  NULLIFY (sm_p)
119  NULLIFY (sm_s)
120  NULLIFY (orbpop)
121  NULLIFY (para_env)
122  NULLIFY (blacs_env)
123 
124  CALL get_qs_env(qs_env=qs_env, &
125  atomic_kind_set=atomic_kind_set, &
126  qs_kind_set=qs_kind_set, &
127  matrix_s_kp=matrixkp_s, &
128  particle_set=particle_set, &
129  rho=rho, &
130  scf_control=scf_control, &
131  para_env=para_env, &
132  blacs_env=blacs_env)
133 
134  cpassert(ASSOCIATED(atomic_kind_set))
135  cpassert(ASSOCIATED(qs_kind_set))
136  cpassert(ASSOCIATED(matrixkp_s))
137  cpassert(ASSOCIATED(particle_set))
138  cpassert(ASSOCIATED(rho))
139  cpassert(ASSOCIATED(scf_control))
140 
141  IF (SIZE(matrixkp_s, 2) > 1) THEN
142 
143  cpwarn("Lowdin population analysis not implemented for k-points.")
144 
145  ELSE
146 
147  sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
148  CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
149 
150  matrix_p => matrixkp_p(:, 1)
151  nspin = SIZE(matrix_p, 1)
152 
153  ! Get the total number of contracted spherical Gaussian basis functions
154  CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
155 
156  ! Provide an array to store the orbital populations for each spin
157  ALLOCATE (orbpop(nsgf, nspin))
158  orbpop(:, :) = 0.0_dp
159 
160  ! Write headline
161  IF (output_unit > 0) THEN
162  WRITE (unit=output_unit, fmt="(/,/,T2,A)") "LOWDIN POPULATION ANALYSIS"
163  END IF
164 
165  ! Provide full size work matrices
166  CALL cp_fm_struct_create(fmstruct=fmstruct, &
167  para_env=para_env, &
168  context=blacs_env, &
169  nrow_global=nsgf, &
170  ncol_global=nsgf)
171  CALL cp_fm_create(matrix=fm_s_half, &
172  matrix_struct=fmstruct, &
173  name="S^(1/2) MATRIX")
174  CALL cp_fm_create(matrix=fm_work1, &
175  matrix_struct=fmstruct, &
176  name="FULL WORK MATRIX 1")
177  headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
178  CALL cp_fm_create(matrix=fm_work2, &
179  matrix_struct=fmstruct, &
180  name=trim(headline))
181  CALL cp_fm_struct_release(fmstruct=fmstruct)
182 
183  ! Build full S^(1/2) matrix (computationally expensive)
184  CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
185  CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
186  IF (ndep /= 0) &
187  CALL cp_warn(__location__, &
188  "Overlap matrix exhibits linear dependencies. At least some "// &
189  "eigenvalues have been quenched.")
190 
191  ! Build Lowdin population matrix for each spin
192  DO ispin = 1, nspin
193  sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format
194  ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
195  CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
196  CALL parallel_gemm(transa="N", &
197  transb="N", &
198  m=nsgf, &
199  n=nsgf, &
200  k=nsgf, &
201  alpha=1.0_dp, &
202  matrix_a=fm_s_half, &
203  matrix_b=fm_work1, &
204  beta=0.0_dp, &
205  matrix_c=fm_work2)
206  IF (print_level > 2) THEN
207  ! Write the full Lowdin population matrix
208  IF (nspin > 1) THEN
209  IF (ispin == 1) THEN
210  fm_work2%name = trim(headline)//" FOR ALPHA SPIN"
211  ELSE
212  fm_work2%name = trim(headline)//" FOR BETA SPIN"
213  END IF
214  END IF
215  CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
216  output_unit=output_unit)
217  END IF
218  CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
219  END DO ! next spin ispin
220 
221  ! Write atomic populations and charges
222  IF (output_unit > 0) THEN
223  print_gop = (print_level > 1) ! Print also orbital populations
224  CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
225  END IF
226 
227  ! Release local working storage
228  CALL cp_fm_release(matrix=fm_s_half)
229  CALL cp_fm_release(matrix=fm_work1)
230  CALL cp_fm_release(matrix=fm_work2)
231  IF (ASSOCIATED(orbpop)) THEN
232  DEALLOCATE (orbpop)
233  END IF
234 
235  END IF
236 
237  CALL timestop(handle)
238 
239  END SUBROUTINE lowdin_population_analysis
240 
241 ! **************************************************************************************************
242 !> \brief Perform a Mulliken population analysis
243 !>
244 !> \param qs_env ...
245 !> \param output_unit ...
246 !> \param print_level ...
247 !> \date 10.07.2010
248 !> \author Matthias Krack (MK)
249 !> \version 1.0
250 ! **************************************************************************************************
251  SUBROUTINE mulliken_population_analysis(qs_env, output_unit, print_level)
252 
253  TYPE(qs_environment_type), POINTER :: qs_env
254  INTEGER, INTENT(IN) :: output_unit, print_level
255 
256  CHARACTER(LEN=*), PARAMETER :: routinen = 'mulliken_population_analysis'
257 
258  CHARACTER(LEN=default_string_length) :: headline
259  INTEGER :: blk, handle, iatom, ic, isgf, ispin, &
260  jatom, jsgf, natom, nsgf, nspin, sgfa, &
261  sgfb
262  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
263  LOGICAL :: found, print_gop
264  REAL(kind=dp) :: ps
265  REAL(kind=dp), DIMENSION(:, :), POINTER :: orbpop, p_block, ps_block, s_block
266  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
267  TYPE(dbcsr_iterator_type) :: iter
268  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
269  TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s
270  TYPE(mp_para_env_type), POINTER :: para_env
271  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
272  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
273  TYPE(qs_rho_type), POINTER :: rho
274 
275  CALL timeset(routinen, handle)
276 
277  NULLIFY (atomic_kind_set)
278  NULLIFY (qs_kind_set)
279  NULLIFY (matrix_p)
280  NULLIFY (matrix_s)
281  NULLIFY (orbpop)
282  NULLIFY (particle_set)
283  NULLIFY (ps_block)
284  NULLIFY (p_block)
285  NULLIFY (rho)
286  NULLIFY (sm_p)
287  NULLIFY (sm_ps)
288  NULLIFY (sm_s)
289  NULLIFY (s_block)
290  NULLIFY (para_env)
291 
292  CALL get_qs_env(qs_env=qs_env, &
293  atomic_kind_set=atomic_kind_set, &
294  qs_kind_set=qs_kind_set, &
295  matrix_s_kp=matrix_s, &
296  particle_set=particle_set, &
297  rho=rho, &
298  para_env=para_env)
299 
300  cpassert(ASSOCIATED(atomic_kind_set))
301  cpassert(ASSOCIATED(qs_kind_set))
302  cpassert(ASSOCIATED(particle_set))
303  cpassert(ASSOCIATED(rho))
304  cpassert(ASSOCIATED(matrix_s))
305 
306  CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
307  nspin = SIZE(matrix_p, 1)
308 
309  ! Get the total number of contracted spherical Gaussian basis functions
310  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
311  CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
312  ALLOCATE (first_sgf_atom(natom))
313  first_sgf_atom(:) = 0
314 
315  CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
316 
317  ! Provide an array to store the orbital populations for each spin
318  ALLOCATE (orbpop(nsgf, nspin))
319  orbpop(:, :) = 0.0_dp
320 
321  ! Write headline
322  IF (output_unit > 0) THEN
323  WRITE (unit=output_unit, fmt="(/,T2,A)") &
324  '!-----------------------------------------------------------------------------!'
325  WRITE (unit=output_unit, fmt="(T22,A)") "Mulliken Population Analysis"
326  END IF
327 
328  ! Create a DBCSR work matrix, if needed
329  IF (print_level > 2) THEN
330  sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
331  ALLOCATE (sm_ps)
332  headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
333  CALL dbcsr_copy(matrix_b=sm_ps, &
334  matrix_a=sm_s, &
335  name=trim(headline))
336  END IF
337 
338  ! Build Mulliken population matrix for each spin
339  DO ispin = 1, nspin
340  DO ic = 1, SIZE(matrix_s, 2)
341  IF (print_level > 2) THEN
342  CALL dbcsr_set(sm_ps, 0.0_dp)
343  END IF
344  sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
345  sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
346  ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
347  ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
348  CALL dbcsr_iterator_start(iter, sm_s)
349  DO WHILE (dbcsr_iterator_blocks_left(iter))
350  CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block, blk)
351  IF (.NOT. (ASSOCIATED(s_block))) cycle
352  CALL dbcsr_get_block_p(matrix=sm_p, &
353  row=iatom, &
354  col=jatom, &
355  block=p_block, &
356  found=found)
357  IF (print_level > 2) THEN
358  CALL dbcsr_get_block_p(matrix=sm_ps, &
359  row=iatom, &
360  col=jatom, &
361  block=ps_block, &
362  found=found)
363  cpassert(ASSOCIATED(ps_block))
364  END IF
365 
366  sgfb = first_sgf_atom(jatom)
367  DO jsgf = 1, SIZE(s_block, 2)
368  DO isgf = 1, SIZE(s_block, 1)
369  ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
370  IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
371  orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
372  END DO
373  sgfb = sgfb + 1
374  END DO
375  IF (iatom /= jatom) THEN
376  sgfa = first_sgf_atom(iatom)
377  DO isgf = 1, SIZE(s_block, 1)
378  DO jsgf = 1, SIZE(s_block, 2)
379  ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
380  orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
381  END DO
382  sgfa = sgfa + 1
383  END DO
384  END IF
385  END DO
386  CALL dbcsr_iterator_stop(iter)
387  END DO
388 
389  IF (print_level > 2) THEN
390  ! Write the full Mulliken net AO and overlap population matrix
391  IF (nspin > 1) THEN
392  IF (ispin == 1) THEN
393  CALL dbcsr_setname(sm_ps, trim(headline)//" For Alpha Spin")
394  ELSE
395  CALL dbcsr_setname(sm_ps, trim(headline)//" For Beta Spin")
396  END IF
397  END IF
398  CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, &
399  output_unit=output_unit)
400  END IF
401  END DO
402 
403  CALL para_env%sum(orbpop)
404 
405  ! Write atomic populations and charges
406  IF (output_unit > 0) THEN
407  print_gop = (print_level > 1) ! Print also orbital populations
408  CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
409  END IF
410 
411  ! Save the Mulliken charges to results
412  CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
413 
414  ! Release local working storage
415  IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
416  IF (ASSOCIATED(orbpop)) THEN
417  DEALLOCATE (orbpop)
418  END IF
419  IF (ALLOCATED(first_sgf_atom)) THEN
420  DEALLOCATE (first_sgf_atom)
421  END IF
422 
423  IF (output_unit > 0) THEN
424  WRITE (unit=output_unit, fmt="(T2,A)") &
425  '!-----------------------------------------------------------------------------!'
426  END IF
427 
428  CALL timestop(handle)
429 
430  END SUBROUTINE mulliken_population_analysis
431 
432 ! **************************************************************************************************
433 !> \brief Save the Mulliken charges in results
434 !>
435 !> \param orbpop ...
436 !> \param atomic_kind_set ...
437 !> \param qs_kind_set ...
438 !> \param particle_set ...
439 !> \param qs_env ...
440 !> \date 27.05.2022
441 !> \author Bo Thomsen (BT)
442 !> \version 1.0
443 ! **************************************************************************************************
444  SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
445  REAL(kind=dp), DIMENSION(:, :), POINTER :: orbpop
446  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
447  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
448  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
449  TYPE(qs_environment_type), POINTER :: qs_env
450 
451  CHARACTER(LEN=default_string_length) :: description
452  INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
453  iso, l, natom, nset, nsgf, nspin
454  INTEGER, DIMENSION(:), POINTER :: nshell
455  INTEGER, DIMENSION(:, :), POINTER :: lshell
456  REAL(kind=dp) :: zeff
457  REAL(kind=dp), DIMENSION(3) :: sumorbpop
458  REAL(kind=dp), DIMENSION(:), POINTER :: charges_save
459  TYPE(cp_result_type), POINTER :: results
460  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
461 
462  NULLIFY (lshell)
463  NULLIFY (nshell)
464  NULLIFY (orb_basis_set)
465 
466  cpassert(ASSOCIATED(orbpop))
467  cpassert(ASSOCIATED(atomic_kind_set))
468  cpassert(ASSOCIATED(particle_set))
469 
470  nspin = SIZE(orbpop, 2)
471 
472  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
473  CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
474  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
475  NULLIFY (results)
476  CALL get_qs_env(qs_env, results=results)
477  ALLOCATE (charges_save(natom))
478 
479  iao = 1
480  DO iatom = 1, natom
481  sumorbpop(:) = 0.0_dp
482  NULLIFY (orb_basis_set)
483  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
484  kind_number=ikind)
485  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
486  IF (ASSOCIATED(orb_basis_set)) THEN
487  CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
488  nset=nset, &
489  nshell=nshell, &
490  l=lshell)
491  isgf = 1
492  DO iset = 1, nset
493  DO ishell = 1, nshell(iset)
494  l = lshell(ishell, iset)
495  DO iso = 1, nso(l)
496  IF (nspin == 1) THEN
497  sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
498  ELSE
499  sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
500  sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
501  END IF
502  isgf = isgf + 1
503  iao = iao + 1
504  END DO
505  END DO
506  END DO
507  IF (nspin == 1) THEN
508  charges_save(iatom) = zeff - sumorbpop(1)
509  ELSE
510  charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
511  END IF
512  END IF ! atom has an orbital basis
513  END DO ! next atom iatom
514 
515  ! Store charges in results
516  description = "[MULLIKEN-CHARGES]"
517  CALL cp_results_erase(results=results, description=description)
518  CALL put_results(results=results, description=description, &
519  values=charges_save)
520 
521  DEALLOCATE (charges_save)
522 
523  END SUBROUTINE save_mulliken_charges
524 
525 ! **************************************************************************************************
526 !> \brief Write atomic orbital populations and net atomic charges
527 !>
528 !> \param orbpop ...
529 !> \param atomic_kind_set ...
530 !> \param qs_kind_set ...
531 !> \param particle_set ...
532 !> \param output_unit ...
533 !> \param print_orbital_contributions ...
534 !> \date 07.07.2010
535 !> \author Matthias Krack (MK)
536 !> \version 1.0
537 ! **************************************************************************************************
538  SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
539  print_orbital_contributions)
540 
541  REAL(kind=dp), DIMENSION(:, :), POINTER :: orbpop
542  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
543  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
544  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
545  INTEGER, INTENT(IN) :: output_unit
546  LOGICAL, INTENT(IN) :: print_orbital_contributions
547 
548  CHARACTER(LEN=*), PARAMETER :: routinen = 'write_orbpop'
549 
550  CHARACTER(LEN=2) :: element_symbol
551  CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
552  INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
553  ishell, iso, l, natom, nset, nsgf, &
554  nspin
555  INTEGER, DIMENSION(:), POINTER :: nshell
556  INTEGER, DIMENSION(:, :), POINTER :: lshell
557  REAL(kind=dp) :: zeff
558  REAL(kind=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
559  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
560 
561  CALL timeset(routinen, handle)
562 
563  NULLIFY (lshell)
564  NULLIFY (nshell)
565  NULLIFY (orb_basis_set)
566  NULLIFY (sgf_symbol)
567 
568  cpassert(ASSOCIATED(orbpop))
569  cpassert(ASSOCIATED(atomic_kind_set))
570  cpassert(ASSOCIATED(particle_set))
571 
572  nspin = SIZE(orbpop, 2)
573 
574  CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
575  CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
576 
577  ! Select and write headline
578  IF (nspin == 1) THEN
579  IF (print_orbital_contributions) THEN
580  WRITE (unit=output_unit, fmt="(/,T2,A)") &
581  "# Orbital AO symbol Orbital population Net charge"
582  ELSE
583  WRITE (unit=output_unit, fmt="(/,T2,A)") &
584  "# Atom Element Kind Atomic population Net charge"
585  END IF
586  ELSE
587  IF (print_orbital_contributions) THEN
588  WRITE (unit=output_unit, fmt="(/,T2,A)") &
589  "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
590  ELSE
591  WRITE (unit=output_unit, fmt="(/,T2,A)") &
592  "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
593  END IF
594  END IF
595 
596  totsumorbpop(:) = 0.0_dp
597 
598  iao = 1
599  DO iatom = 1, natom
600  sumorbpop(:) = 0.0_dp
601  NULLIFY (orb_basis_set)
602  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
603  element_symbol=element_symbol, &
604  kind_number=ikind)
605  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
606  IF (ASSOCIATED(orb_basis_set)) THEN
607  CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
608  nset=nset, &
609  nshell=nshell, &
610  l=lshell, &
611  sgf_symbol=sgf_symbol)
612  isgf = 1
613  DO iset = 1, nset
614  DO ishell = 1, nshell(iset)
615  l = lshell(ishell, iset)
616  DO iso = 1, nso(l)
617  IF (nspin == 1) THEN
618  sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
619  IF (print_orbital_contributions) THEN
620  IF (isgf == 1) WRITE (unit=output_unit, fmt="(A)") ""
621  WRITE (unit=output_unit, &
622  fmt="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
623  iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
624  END IF
625  ELSE
626  sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
627  sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
628  IF (print_orbital_contributions) THEN
629  IF (isgf == 1) WRITE (unit=output_unit, fmt="(A)") ""
630  WRITE (unit=output_unit, &
631  fmt="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
632  iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
633  orbpop(iao, 1) - orbpop(iao, 2)
634  END IF
635  END IF
636  isgf = isgf + 1
637  iao = iao + 1
638  END DO
639  END DO
640  END DO
641  IF (nspin == 1) THEN
642  totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
643  totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
644  WRITE (unit=output_unit, &
645  fmt="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
646  iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
647  ELSE
648  totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
649  totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
650  WRITE (unit=output_unit, &
651  fmt="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
652  iatom, element_symbol, ikind, sumorbpop(1:2), &
653  zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
654  END IF
655  END IF ! atom has an orbital basis
656  END DO ! next atom iatom
657 
658  ! Write total sums
659  IF (print_orbital_contributions) WRITE (unit=output_unit, fmt="(A)") ""
660  IF (nspin == 1) THEN
661  WRITE (unit=output_unit, &
662  fmt="(T2,A,T42,F12.6,T68,F12.6,/)") &
663  "# Total charge", totsumorbpop(1), totsumorbpop(3)
664  ELSE
665  WRITE (unit=output_unit, &
666  fmt="(T2,A,T28,4(1X,F12.6),/)") &
667  "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
668  totsumorbpop(1) - totsumorbpop(2)
669  END IF
670 
671  IF (output_unit > 0) CALL m_flush(output_unit)
672 
673  CALL timestop(handle)
674 
675  END SUBROUTINE write_orbpop
676 
677 END MODULE population_analyses
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)
...
Definition: basis_set_types.F:615
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
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:744
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:208
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:78
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)
...
Definition: cp_dbcsr_output.F:163
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:896
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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:570
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
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:362
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
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:106
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)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
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:83
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:252
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_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, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
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
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17