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 110 : 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 110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
97 110 : 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 110 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
102 110 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
103 110 : 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 110 : POINTER :: sab_nl
109 110 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 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 110 : CALL timeset(routineN, handle)
115 :
116 110 : 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 110 : blacs_env=blacs_env)
128 :
129 110 : CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
130 110 : nspin = SIZE(matrixkp_p, 1)
131 110 : nimg = SIZE(matrixkp_p, 2)
132 :
133 : ! Get the total number of contracted spherical Gaussian basis functions
134 110 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
135 : ! Provide an array to store the orbital populations for each spin
136 440 : ALLOCATE (orbpop(nsgf, nspin))
137 110 : orbpop(:, :) = 0.0_dp
138 :
139 : ! Write headline
140 110 : IF (output_unit > 0) THEN
141 55 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "LOWDIN POPULATION ANALYSIS"
142 : END IF
143 :
144 110 : IF (do_kpoints) THEN
145 :
146 6 : 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 6 : nrow_global=nsgf, ncol_global=nsgf)
151 30 : ALLOCATE (fmwork(4))
152 30 : DO i = 1, 4
153 30 : CALL cp_fm_create(fmwork(i), matrix_struct=fmstruct, name="work")
154 : END DO
155 6 : CALL cp_fm_struct_release(fmstruct)
156 :
157 : ! S^1/2
158 6 : CALL diag_kp_smat(matrixkp_s, kpoints, fmwork)
159 : ! Lowdin P Matrix Transform
160 6 : CALL lowdin_kp_trans(kpoints, orbpop)
161 :
162 6 : 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) THEN
191 : CALL cp_warn(__LOCATION__, &
192 : "Overlap matrix exhibits linear dependencies. At least some "// &
193 0 : "eigenvalues have been quenched.")
194 : END IF
195 :
196 : ! Build Lowdin population matrix for each spin
197 212 : DO ispin = 1, nspin
198 108 : 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 108 : 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 108 : matrix_c=fm_work2)
211 108 : IF (print_level > 2) THEN
212 : ! Write the full Lowdin population matrix
213 4 : IF (nspin > 1) THEN
214 4 : IF (ispin == 1) THEN
215 2 : fm_work2%name = TRIM(headline)//" FOR ALPHA SPIN"
216 : ELSE
217 2 : 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 4 : output_unit=output_unit)
222 : END IF
223 212 : CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
224 : END DO ! next spin ispin
225 :
226 : ! Release local working storage
227 104 : CALL cp_fm_release(matrix=fm_s_half)
228 104 : CALL cp_fm_release(matrix=fm_work1)
229 104 : CALL cp_fm_release(matrix=fm_work2)
230 :
231 : END IF
232 :
233 : ! Write atomic populations and charges
234 110 : IF (output_unit > 0) THEN
235 55 : print_gop = (print_level > 1) ! Print also orbital populations
236 55 : CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
237 : END IF
238 :
239 110 : DEALLOCATE (orbpop)
240 :
241 110 : CALL timestop(handle)
242 :
243 330 : 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 5302 : 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 5302 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
266 : LOGICAL :: found, print_gop
267 : REAL(KIND=dp) :: ps
268 5302 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbpop
269 5302 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, ps_block, s_block
270 5302 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
271 : TYPE(dbcsr_iterator_type) :: iter
272 5302 : 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 5302 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
276 5302 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
277 : TYPE(qs_rho_type), POINTER :: rho
278 :
279 5302 : CALL timeset(routineN, handle)
280 :
281 5302 : NULLIFY (sm_p, sm_ps, sm_s)
282 5302 : 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 5302 : para_env=para_env)
291 :
292 5302 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
293 5302 : nspin = SIZE(matrix_p, 1)
294 :
295 : ! Get the total number of contracted spherical Gaussian basis functions
296 5302 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
297 5302 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
298 15906 : ALLOCATE (first_sgf_atom(natom))
299 5302 : first_sgf_atom(:) = 0
300 :
301 5302 : 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 21208 : ALLOCATE (orbpop(nsgf, nspin))
305 5302 : orbpop(:, :) = 0.0_dp
306 :
307 : ! Write headline
308 5302 : IF (output_unit > 0) THEN
309 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
310 2666 : '!-----------------------------------------------------------------------------!'
311 2666 : WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
312 : END IF
313 :
314 : ! Create a DBCSR work matrix, if needed
315 5302 : IF (print_level > 2) THEN
316 2 : sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
317 2 : ALLOCATE (sm_ps)
318 2 : headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
319 2 : IF (nspin > 1) THEN
320 2 : IF (ispin == 1) THEN
321 0 : headline = TRIM(headline)//" For Alpha Spin"
322 : ELSE
323 2 : headline = TRIM(headline)//" For Beta Spin"
324 : END IF
325 : END IF
326 2 : 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 11376 : DO ispin = 1, nspin
331 24168 : DO ic = 1, SIZE(matrix_s, 2)
332 18094 : IF (print_level > 2) THEN
333 4 : CALL dbcsr_set(sm_ps, 0.0_dp)
334 : END IF
335 18094 : sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
336 18094 : 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 18094 : CALL dbcsr_iterator_start(iter, sm_s)
340 176164 : DO WHILE (dbcsr_iterator_blocks_left(iter))
341 158070 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block)
342 158070 : 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 158070 : found=found)
348 158070 : 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 12 : found=found)
354 12 : CPASSERT(ASSOCIATED(ps_block))
355 : END IF
356 :
357 158070 : sgfb = first_sgf_atom(jatom)
358 862126 : DO jsgf = 1, SIZE(s_block, 2)
359 6636451 : DO isgf = 1, SIZE(s_block, 1)
360 5932395 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
361 5932395 : IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
362 6636451 : orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
363 : END DO
364 862126 : sgfb = sgfb + 1
365 : END DO
366 176164 : IF (iatom /= jatom) THEN
367 121901 : sgfa = first_sgf_atom(iatom)
368 650578 : DO isgf = 1, SIZE(s_block, 1)
369 4127591 : DO jsgf = 1, SIZE(s_block, 2)
370 3598914 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
371 4127591 : orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
372 : END DO
373 650578 : sgfa = sgfa + 1
374 : END DO
375 : END IF
376 : END DO
377 42262 : CALL dbcsr_iterator_stop(iter)
378 : END DO
379 :
380 11376 : IF (print_level > 2) THEN
381 : ! Write the full Mulliken net AO and overlap population matrix
382 4 : 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 5302 : CALL para_env%sum(orbpop)
387 :
388 : ! Write atomic populations and charges
389 5302 : IF (output_unit > 0) THEN
390 2666 : print_gop = (print_level > 1) ! Print also orbital populations
391 2666 : 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 5302 : CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
396 :
397 : ! Release local working storage
398 5302 : IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
399 5302 : DEALLOCATE (orbpop)
400 5302 : DEALLOCATE (first_sgf_atom)
401 :
402 5302 : IF (output_unit > 0) THEN
403 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
404 2666 : '!-----------------------------------------------------------------------------!'
405 : END IF
406 :
407 5302 : CALL timestop(handle)
408 :
409 15906 : 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 5302 : 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 5302 : INTEGER, DIMENSION(:), POINTER :: nshell
437 5302 : 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 5302 : nspin = SIZE(orbpop, 2)
445 :
446 5302 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
447 5302 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
448 5302 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
449 5302 : NULLIFY (results)
450 5302 : CALL get_qs_env(qs_env, results=results)
451 15906 : ALLOCATE (all_sumorbpop(natom))
452 10604 : ALLOCATE (charges_save(natom))
453 :
454 5302 : iao = 1
455 27956 : DO iatom = 1, natom
456 22654 : sumorbpop(:) = 0.0_dp
457 22654 : NULLIFY (orb_basis_set)
458 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
459 22654 : kind_number=ikind)
460 22654 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
461 27956 : 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 22654 : l=lshell)
466 22654 : isgf = 1
467 61022 : DO iset = 1, nset
468 128806 : DO ishell = 1, nshell(iset)
469 67784 : l = lshell(ishell, iset)
470 246256 : DO iso = 1, nso(l)
471 140104 : IF (nspin == 1) THEN
472 117684 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
473 : ELSE
474 67260 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
475 22420 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
476 : END IF
477 140104 : isgf = isgf + 1
478 207888 : iao = iao + 1
479 : END DO
480 : END DO
481 : END DO
482 22654 : IF (nspin == 1) THEN
483 19742 : charges_save(iatom) = zeff - sumorbpop(1)
484 19742 : all_sumorbpop(iatom) = sumorbpop(1)
485 : ELSE
486 2912 : charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
487 2912 : 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 5302 : description = "[MULLIKEN-ORBPOP]"
494 5302 : CALL cp_results_erase(results=results, description=description)
495 : CALL put_results(results=results, description=description, &
496 5302 : values=orbpop)
497 :
498 : ! Store sum orbital population in results
499 5302 : description = "[MULLIKEN-SUMORBPOP]"
500 5302 : CALL cp_results_erase(results=results, description=description)
501 : CALL put_results(results=results, description=description, &
502 5302 : values=all_sumorbpop)
503 :
504 : ! Store charges in results
505 5302 : description = "[MULLIKEN-CHARGES]"
506 5302 : CALL cp_results_erase(results=results, description=description)
507 : CALL put_results(results=results, description=description, &
508 5302 : values=charges_save)
509 :
510 5302 : DEALLOCATE (all_sumorbpop)
511 5302 : DEALLOCATE (charges_save)
512 :
513 5302 : 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 2721 : 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 2721 : 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 2721 : INTEGER, DIMENSION(:), POINTER :: nshell
546 2721 : 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 2721 : CALL timeset(routineN, handle)
552 :
553 2721 : nspin = SIZE(orbpop, 2)
554 :
555 2721 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
556 2721 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
557 :
558 : ! Select and write headline
559 2721 : IF (nspin == 1) THEN
560 2331 : IF (print_orbital_contributions) THEN
561 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
562 0 : "# Orbital AO symbol Orbital population Net charge"
563 : ELSE
564 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
565 2331 : "# Atom Element Kind Atomic population Net charge"
566 : END IF
567 : ELSE
568 390 : IF (print_orbital_contributions) THEN
569 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
570 4 : "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
571 : ELSE
572 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
573 386 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
574 : END IF
575 : END IF
576 :
577 2721 : totsumorbpop(:) = 0.0_dp
578 :
579 2721 : iao = 1
580 14246 : DO iatom = 1, natom
581 11525 : sumorbpop(:) = 0.0_dp
582 11525 : NULLIFY (orb_basis_set)
583 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
584 : element_symbol=element_symbol, &
585 11525 : kind_number=ikind)
586 11525 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
587 14246 : 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 11525 : sgf_symbol=sgf_symbol)
593 11525 : isgf = 1
594 31280 : DO iset = 1, nset
595 65893 : DO ishell = 1, nshell(iset)
596 34613 : l = lshell(ishell, iset)
597 126013 : DO iso = 1, nso(l)
598 71645 : IF (nspin == 1) THEN
599 60292 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
600 60292 : IF (print_orbital_contributions) THEN
601 0 : 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 0 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
605 : END IF
606 : ELSE
607 34059 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
608 11353 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
609 11353 : IF (print_orbital_contributions) THEN
610 156 : 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 156 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
614 312 : orbpop(iao, 1) - orbpop(iao, 2)
615 : END IF
616 : END IF
617 71645 : isgf = isgf + 1
618 106258 : iao = iao + 1
619 : END DO
620 : END DO
621 : END DO
622 11525 : IF (nspin == 1) THEN
623 10058 : totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
624 10058 : 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 10058 : iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
628 : ELSE
629 4401 : totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
630 1467 : 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 1467 : iatom, element_symbol, ikind, sumorbpop(1:2), &
634 2934 : 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 2721 : IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
641 2721 : IF (nspin == 1) THEN
642 : WRITE (UNIT=output_unit, &
643 : FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
644 2331 : "# Total charge", totsumorbpop(1), totsumorbpop(3)
645 : ELSE
646 : WRITE (UNIT=output_unit, &
647 : FMT="(T2,A,T28,4(1X,F12.6),/)") &
648 390 : "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
649 780 : totsumorbpop(1) - totsumorbpop(2)
650 : END IF
651 :
652 2721 : IF (output_unit > 0) CALL m_flush(output_unit)
653 :
654 2721 : CALL timestop(handle)
655 :
656 2721 : END SUBROUTINE write_orbpop
657 :
658 : END MODULE population_analyses
|
