Line data Source code
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 :
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_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 :
73 : CONTAINS
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 108 : 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 108 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
97 108 : 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 108 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
102 108 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
103 108 : 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 108 : POINTER :: sab_nl
109 108 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 108 : 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 108 : CALL timeset(routineN, handle)
115 :
116 108 : 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 108 : blacs_env=blacs_env)
128 :
129 108 : CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
130 108 : nspin = SIZE(matrixkp_p, 1)
131 108 : nimg = SIZE(matrixkp_p, 2)
132 :
133 : ! Get the total number of contracted spherical Gaussian basis functions
134 108 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
135 : ! Provide an array to store the orbital populations for each spin
136 432 : ALLOCATE (orbpop(nsgf, nspin))
137 3238 : orbpop(:, :) = 0.0_dp
138 :
139 : ! Write headline
140 108 : IF (output_unit > 0) THEN
141 54 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "LOWDIN POPULATION ANALYSIS"
142 : END IF
143 :
144 108 : IF (do_kpoints) THEN
145 :
146 4 : 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 4 : nrow_global=nsgf, ncol_global=nsgf)
151 20 : ALLOCATE (fmwork(4))
152 20 : DO i = 1, 4
153 20 : CALL cp_fm_create(fmwork(i), matrix_struct=fmstruct, name="work")
154 : END DO
155 4 : CALL cp_fm_struct_release(fmstruct)
156 :
157 : ! S^1/2
158 4 : CALL diag_kp_smat(matrixkp_s, kpoints, fmwork)
159 : ! Lowdin P Matrix Transform
160 4 : CALL lowdin_kp_trans(kpoints, orbpop)
161 :
162 4 : CALL cp_fm_release(fmwork)
163 :
164 : ELSE
165 :
166 104 : sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
167 104 : 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 104 : ncol_global=nsgf)
175 : CALL cp_fm_create(matrix=fm_s_half, &
176 : matrix_struct=fmstruct, &
177 104 : name="S^(1/2) MATRIX")
178 : CALL cp_fm_create(matrix=fm_work1, &
179 : matrix_struct=fmstruct, &
180 104 : name="FULL WORK MATRIX 1")
181 104 : headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
182 : CALL cp_fm_create(matrix=fm_work2, &
183 : matrix_struct=fmstruct, &
184 104 : name=TRIM(headline))
185 104 : CALL cp_fm_struct_release(fmstruct=fmstruct)
186 :
187 : ! Build full S^(1/2) matrix (computationally expensive)
188 104 : CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
189 104 : CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
190 104 : IF (ndep /= 0) &
191 : CALL cp_warn(__LOCATION__, &
192 : "Overlap matrix exhibits linear dependencies. At least some "// &
193 0 : "eigenvalues have been quenched.")
194 :
195 : ! Build Lowdin population matrix for each spin
196 212 : DO ispin = 1, nspin
197 108 : sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format
198 : ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
199 108 : CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
200 : CALL parallel_gemm(transa="N", &
201 : transb="N", &
202 : m=nsgf, &
203 : n=nsgf, &
204 : k=nsgf, &
205 : alpha=1.0_dp, &
206 : matrix_a=fm_s_half, &
207 : matrix_b=fm_work1, &
208 : beta=0.0_dp, &
209 108 : matrix_c=fm_work2)
210 108 : IF (print_level > 2) THEN
211 : ! Write the full Lowdin population matrix
212 4 : IF (nspin > 1) THEN
213 4 : IF (ispin == 1) THEN
214 2 : fm_work2%name = TRIM(headline)//" FOR ALPHA SPIN"
215 : ELSE
216 2 : fm_work2%name = TRIM(headline)//" FOR BETA SPIN"
217 : END IF
218 : END IF
219 : CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
220 4 : output_unit=output_unit)
221 : END IF
222 212 : CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
223 : END DO ! next spin ispin
224 :
225 : ! Release local working storage
226 104 : CALL cp_fm_release(matrix=fm_s_half)
227 104 : CALL cp_fm_release(matrix=fm_work1)
228 104 : CALL cp_fm_release(matrix=fm_work2)
229 :
230 : END IF
231 :
232 : ! Write atomic populations and charges
233 108 : IF (output_unit > 0) THEN
234 54 : print_gop = (print_level > 1) ! Print also orbital populations
235 54 : CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
236 : END IF
237 :
238 108 : DEALLOCATE (orbpop)
239 :
240 108 : CALL timestop(handle)
241 :
242 324 : END SUBROUTINE lowdin_population_analysis
243 :
244 : ! **************************************************************************************************
245 : !> \brief Perform a Mulliken population analysis
246 : !>
247 : !> \param qs_env ...
248 : !> \param output_unit ...
249 : !> \param print_level ...
250 : !> \date 10.07.2010
251 : !> \author Matthias Krack (MK)
252 : !> \version 1.0
253 : ! **************************************************************************************************
254 5030 : SUBROUTINE mulliken_population_analysis(qs_env, output_unit, print_level)
255 :
256 : TYPE(qs_environment_type), POINTER :: qs_env
257 : INTEGER, INTENT(IN) :: output_unit, print_level
258 :
259 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_population_analysis'
260 :
261 : CHARACTER(LEN=default_string_length) :: headline
262 : INTEGER :: handle, iatom, ic, isgf, ispin, jatom, &
263 : jsgf, natom, nsgf, nspin, sgfa, sgfb
264 5030 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
265 : LOGICAL :: found, print_gop
266 : REAL(KIND=dp) :: ps
267 5030 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
268 5030 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, ps_block, s_block
269 5030 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
270 : TYPE(dbcsr_iterator_type) :: iter
271 5030 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
272 : TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s
273 : TYPE(mp_para_env_type), POINTER :: para_env
274 5030 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
275 5030 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
276 : TYPE(qs_rho_type), POINTER :: rho
277 :
278 5030 : CALL timeset(routineN, handle)
279 :
280 5030 : NULLIFY (sm_p, sm_ps, sm_s)
281 5030 : NULLIFY (p_block, s_block, ps_block)
282 :
283 : CALL get_qs_env(qs_env=qs_env, &
284 : atomic_kind_set=atomic_kind_set, &
285 : qs_kind_set=qs_kind_set, &
286 : matrix_s_kp=matrix_s, &
287 : particle_set=particle_set, &
288 : rho=rho, &
289 5030 : para_env=para_env)
290 :
291 5030 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
292 5030 : nspin = SIZE(matrix_p, 1)
293 :
294 : ! Get the total number of contracted spherical Gaussian basis functions
295 5030 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
296 5030 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
297 15090 : ALLOCATE (first_sgf_atom(natom))
298 26182 : first_sgf_atom(:) = 0
299 :
300 5030 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
301 :
302 : ! Provide an array to store the orbital populations for each spin
303 20120 : ALLOCATE (orbpop(nsgf, nspin))
304 175836 : orbpop(:, :) = 0.0_dp
305 :
306 : ! Write headline
307 5030 : IF (output_unit > 0) THEN
308 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
309 2529 : '!-----------------------------------------------------------------------------!'
310 2529 : WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
311 : END IF
312 :
313 : ! Create a DBCSR work matrix, if needed
314 5030 : IF (print_level > 2) THEN
315 2 : sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
316 2 : ALLOCATE (sm_ps)
317 2 : headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
318 2 : IF (nspin > 1) THEN
319 2 : IF (ispin == 1) THEN
320 0 : headline = TRIM(headline)//" For Alpha Spin"
321 : ELSE
322 2 : headline = TRIM(headline)//" For Beta Spin"
323 : END IF
324 : END IF
325 2 : CALL dbcsr_copy(matrix_b=sm_ps, matrix_a=sm_s, name=TRIM(headline))
326 : END IF
327 :
328 : ! Build Mulliken population matrix for each spin
329 10824 : DO ispin = 1, nspin
330 17320 : DO ic = 1, SIZE(matrix_s, 2)
331 11526 : IF (print_level > 2) THEN
332 4 : CALL dbcsr_set(sm_ps, 0.0_dp)
333 : END IF
334 11526 : sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
335 11526 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
336 : ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
337 : ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
338 11526 : CALL dbcsr_iterator_start(iter, sm_s)
339 105530 : DO WHILE (dbcsr_iterator_blocks_left(iter))
340 94004 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block)
341 94004 : IF (.NOT. (ASSOCIATED(s_block))) CYCLE
342 : CALL dbcsr_get_block_p(matrix=sm_p, &
343 : row=iatom, &
344 : col=jatom, &
345 : block=p_block, &
346 94004 : found=found)
347 94004 : IF (print_level > 2) THEN
348 : CALL dbcsr_get_block_p(matrix=sm_ps, &
349 : row=iatom, &
350 : col=jatom, &
351 : block=ps_block, &
352 12 : found=found)
353 12 : CPASSERT(ASSOCIATED(ps_block))
354 : END IF
355 :
356 94004 : sgfb = first_sgf_atom(jatom)
357 555906 : DO jsgf = 1, SIZE(s_block, 2)
358 5972682 : DO isgf = 1, SIZE(s_block, 1)
359 5510780 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
360 5510780 : IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
361 5972682 : orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
362 : END DO
363 555906 : sgfb = sgfb + 1
364 : END DO
365 105530 : IF (iatom /= jatom) THEN
366 74412 : sgfa = first_sgf_atom(iatom)
367 422459 : DO isgf = 1, SIZE(s_block, 1)
368 3861376 : DO jsgf = 1, SIZE(s_block, 2)
369 3513329 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
370 3861376 : orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
371 : END DO
372 422459 : sgfa = sgfa + 1
373 : END DO
374 : END IF
375 : END DO
376 28846 : CALL dbcsr_iterator_stop(iter)
377 : END DO
378 :
379 10824 : IF (print_level > 2) THEN
380 : ! Write the full Mulliken net AO and overlap population matrix
381 4 : CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, output_unit=output_unit)
382 : END IF
383 : END DO
384 :
385 5030 : CALL para_env%sum(orbpop)
386 :
387 : ! Write atomic populations and charges
388 5030 : IF (output_unit > 0) THEN
389 2529 : print_gop = (print_level > 1) ! Print also orbital populations
390 2529 : CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
391 : END IF
392 :
393 : ! Save the Mulliken charges to results
394 5030 : CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
395 :
396 : ! Release local working storage
397 5030 : IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
398 5030 : DEALLOCATE (orbpop)
399 5030 : DEALLOCATE (first_sgf_atom)
400 :
401 5030 : IF (output_unit > 0) THEN
402 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
403 2529 : '!-----------------------------------------------------------------------------!'
404 : END IF
405 :
406 5030 : CALL timestop(handle)
407 :
408 15090 : END SUBROUTINE mulliken_population_analysis
409 :
410 : ! **************************************************************************************************
411 : !> \brief Save the Mulliken atomic orbital populations and charges in results
412 : !>
413 : !> \param orbpop ...
414 : !> \param atomic_kind_set ...
415 : !> \param qs_kind_set ...
416 : !> \param particle_set ...
417 : !> \param qs_env ...
418 : !> \par History
419 : !> 27.05.2022 BT
420 : !> 16.07.2025 RK
421 : !> \author Bo Thomsen (BT)
422 : !> Rangsiman Ketkaew (RK)
423 : !> \version 1.0
424 : ! **************************************************************************************************
425 5030 : SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
426 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: orbpop
427 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
428 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
429 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
430 : TYPE(qs_environment_type), POINTER :: qs_env
431 :
432 : CHARACTER(LEN=default_string_length) :: description
433 : INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
434 : iso, l, natom, nset, nsgf, nspin
435 5030 : INTEGER, DIMENSION(:), POINTER :: nshell
436 5030 : INTEGER, DIMENSION(:, :), POINTER :: lshell
437 : REAL(KIND=dp) :: zeff
438 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: all_sumorbpop, charges_save
439 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop
440 : TYPE(cp_result_type), POINTER :: results
441 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
442 :
443 5030 : nspin = SIZE(orbpop, 2)
444 :
445 5030 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
446 5030 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
447 5030 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
448 5030 : NULLIFY (results)
449 5030 : CALL get_qs_env(qs_env, results=results)
450 15090 : ALLOCATE (all_sumorbpop(natom))
451 10060 : ALLOCATE (charges_save(natom))
452 :
453 5030 : iao = 1
454 26182 : DO iatom = 1, natom
455 21152 : sumorbpop(:) = 0.0_dp
456 21152 : NULLIFY (orb_basis_set)
457 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
458 21152 : kind_number=ikind)
459 21152 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
460 26182 : IF (ASSOCIATED(orb_basis_set)) THEN
461 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
462 : nset=nset, &
463 : nshell=nshell, &
464 21152 : l=lshell)
465 21152 : isgf = 1
466 57538 : DO iset = 1, nset
467 123790 : DO ishell = 1, nshell(iset)
468 66252 : l = lshell(ishell, iset)
469 243718 : DO iso = 1, nso(l)
470 141080 : IF (nspin == 1) THEN
471 117148 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
472 : ELSE
473 71796 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
474 23932 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
475 : END IF
476 141080 : isgf = isgf + 1
477 207332 : iao = iao + 1
478 : END DO
479 : END DO
480 : END DO
481 21152 : IF (nspin == 1) THEN
482 18260 : charges_save(iatom) = zeff - sumorbpop(1)
483 18260 : all_sumorbpop(iatom) = sumorbpop(1)
484 : ELSE
485 2892 : charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
486 2892 : all_sumorbpop(iatom) = sumorbpop(1) + sumorbpop(2)
487 : END IF
488 : END IF ! atom has an orbital basis
489 : END DO ! next atom iatom
490 :
491 : ! Store atomic orbital populations in results
492 5030 : description = "[MULLIKEN-ORBPOP]"
493 5030 : CALL cp_results_erase(results=results, description=description)
494 : CALL put_results(results=results, description=description, &
495 5030 : values=orbpop)
496 :
497 : ! Store sum orbital population in results
498 5030 : description = "[MULLIKEN-SUMORBPOP]"
499 5030 : CALL cp_results_erase(results=results, description=description)
500 : CALL put_results(results=results, description=description, &
501 5030 : values=all_sumorbpop)
502 :
503 : ! Store charges in results
504 5030 : description = "[MULLIKEN-CHARGES]"
505 5030 : CALL cp_results_erase(results=results, description=description)
506 : CALL put_results(results=results, description=description, &
507 5030 : values=charges_save)
508 :
509 5030 : DEALLOCATE (all_sumorbpop)
510 5030 : DEALLOCATE (charges_save)
511 :
512 5030 : END SUBROUTINE save_mulliken_charges
513 :
514 : ! **************************************************************************************************
515 : !> \brief Write atomic orbital populations and net atomic charges
516 : !>
517 : !> \param orbpop ...
518 : !> \param atomic_kind_set ...
519 : !> \param qs_kind_set ...
520 : !> \param particle_set ...
521 : !> \param output_unit ...
522 : !> \param print_orbital_contributions ...
523 : !> \date 07.07.2010
524 : !> \author Matthias Krack (MK)
525 : !> \version 1.0
526 : ! **************************************************************************************************
527 2583 : SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
528 : print_orbital_contributions)
529 :
530 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: orbpop
531 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
532 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
533 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
534 : INTEGER, INTENT(IN) :: output_unit
535 : LOGICAL, INTENT(IN) :: print_orbital_contributions
536 :
537 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_orbpop'
538 :
539 : CHARACTER(LEN=2) :: element_symbol
540 2583 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
541 : INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
542 : ishell, iso, l, natom, nset, nsgf, &
543 : nspin
544 2583 : INTEGER, DIMENSION(:), POINTER :: nshell
545 2583 : INTEGER, DIMENSION(:, :), POINTER :: lshell
546 : REAL(KIND=dp) :: zeff
547 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
548 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
549 :
550 2583 : CALL timeset(routineN, handle)
551 :
552 2583 : nspin = SIZE(orbpop, 2)
553 :
554 2583 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
555 2583 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
556 :
557 : ! Select and write headline
558 2583 : IF (nspin == 1) THEN
559 2198 : IF (print_orbital_contributions) THEN
560 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
561 0 : "# Orbital AO symbol Orbital population Net charge"
562 : ELSE
563 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
564 2198 : "# Atom Element Kind Atomic population Net charge"
565 : END IF
566 : ELSE
567 385 : IF (print_orbital_contributions) THEN
568 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
569 3 : "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
570 : ELSE
571 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
572 382 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
573 : END IF
574 : END IF
575 :
576 2583 : totsumorbpop(:) = 0.0_dp
577 :
578 2583 : iao = 1
579 13352 : DO iatom = 1, natom
580 10769 : sumorbpop(:) = 0.0_dp
581 10769 : NULLIFY (orb_basis_set)
582 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
583 : element_symbol=element_symbol, &
584 10769 : kind_number=ikind)
585 10769 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
586 13352 : IF (ASSOCIATED(orb_basis_set)) THEN
587 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
588 : nset=nset, &
589 : nshell=nshell, &
590 : l=lshell, &
591 10769 : sgf_symbol=sgf_symbol)
592 10769 : isgf = 1
593 29523 : DO iset = 1, nset
594 63349 : DO ishell = 1, nshell(iset)
595 33826 : l = lshell(ishell, iset)
596 124664 : DO iso = 1, nso(l)
597 72084 : IF (nspin == 1) THEN
598 60014 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
599 60014 : IF (print_orbital_contributions) THEN
600 0 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
601 : WRITE (UNIT=output_unit, &
602 : FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
603 0 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
604 : END IF
605 : ELSE
606 36210 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
607 12070 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
608 12070 : IF (print_orbital_contributions) THEN
609 117 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
610 : WRITE (UNIT=output_unit, &
611 : FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
612 117 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
613 234 : orbpop(iao, 1) - orbpop(iao, 2)
614 : END IF
615 : END IF
616 72084 : isgf = isgf + 1
617 105910 : iao = iao + 1
618 : END DO
619 : END DO
620 : END DO
621 10769 : IF (nspin == 1) THEN
622 9315 : totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
623 9315 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
624 : WRITE (UNIT=output_unit, &
625 : FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
626 9315 : iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
627 : ELSE
628 4362 : totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
629 1454 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
630 : WRITE (UNIT=output_unit, &
631 : FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
632 1454 : iatom, element_symbol, ikind, sumorbpop(1:2), &
633 2908 : zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
634 : END IF
635 : END IF ! atom has an orbital basis
636 : END DO ! next atom iatom
637 :
638 : ! Write total sums
639 2583 : IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
640 2583 : IF (nspin == 1) THEN
641 : WRITE (UNIT=output_unit, &
642 : FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
643 2198 : "# Total charge", totsumorbpop(1), totsumorbpop(3)
644 : ELSE
645 : WRITE (UNIT=output_unit, &
646 : FMT="(T2,A,T28,4(1X,F12.6),/)") &
647 385 : "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
648 770 : totsumorbpop(1) - totsumorbpop(2)
649 : END IF
650 :
651 2583 : IF (output_unit > 0) CALL m_flush(output_unit)
652 :
653 2583 : CALL timestop(handle)
654 :
655 2583 : END SUBROUTINE write_orbpop
656 :
657 : END MODULE population_analyses
|
