Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \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 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 102 : 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 102 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
93 102 : 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 102 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
98 102 : 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 102 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
102 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 102 : CALL timeset(routineN, handle)
107 :
108 102 : NULLIFY (atomic_kind_set)
109 102 : NULLIFY (qs_kind_set)
110 102 : NULLIFY (fmstruct)
111 102 : NULLIFY (matrix_p)
112 102 : NULLIFY (matrixkp_p)
113 102 : NULLIFY (matrixkp_s)
114 102 : NULLIFY (orbpop)
115 102 : NULLIFY (particle_set)
116 102 : NULLIFY (rho)
117 102 : NULLIFY (scf_control)
118 102 : NULLIFY (sm_p)
119 102 : NULLIFY (sm_s)
120 : NULLIFY (orbpop)
121 102 : NULLIFY (para_env)
122 102 : 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 102 : blacs_env=blacs_env)
133 :
134 102 : CPASSERT(ASSOCIATED(atomic_kind_set))
135 102 : CPASSERT(ASSOCIATED(qs_kind_set))
136 102 : CPASSERT(ASSOCIATED(matrixkp_s))
137 102 : CPASSERT(ASSOCIATED(particle_set))
138 102 : CPASSERT(ASSOCIATED(rho))
139 102 : CPASSERT(ASSOCIATED(scf_control))
140 :
141 102 : IF (SIZE(matrixkp_s, 2) > 1) THEN
142 :
143 0 : CPWARN("Lowdin population analysis not implemented for k-points.")
144 :
145 : ELSE
146 :
147 102 : sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
148 102 : CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
149 :
150 102 : matrix_p => matrixkp_p(:, 1)
151 102 : nspin = SIZE(matrix_p, 1)
152 :
153 : ! Get the total number of contracted spherical Gaussian basis functions
154 102 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
155 :
156 : ! Provide an array to store the orbital populations for each spin
157 408 : ALLOCATE (orbpop(nsgf, nspin))
158 3046 : orbpop(:, :) = 0.0_dp
159 :
160 : ! Write headline
161 102 : IF (output_unit > 0) THEN
162 51 : 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 102 : ncol_global=nsgf)
171 : CALL cp_fm_create(matrix=fm_s_half, &
172 : matrix_struct=fmstruct, &
173 102 : name="S^(1/2) MATRIX")
174 : CALL cp_fm_create(matrix=fm_work1, &
175 : matrix_struct=fmstruct, &
176 102 : name="FULL WORK MATRIX 1")
177 102 : headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
178 : CALL cp_fm_create(matrix=fm_work2, &
179 : matrix_struct=fmstruct, &
180 102 : name=TRIM(headline))
181 102 : CALL cp_fm_struct_release(fmstruct=fmstruct)
182 :
183 : ! Build full S^(1/2) matrix (computationally expensive)
184 102 : CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
185 102 : CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
186 102 : IF (ndep /= 0) &
187 : CALL cp_warn(__LOCATION__, &
188 : "Overlap matrix exhibits linear dependencies. At least some "// &
189 0 : "eigenvalues have been quenched.")
190 :
191 : ! Build Lowdin population matrix for each spin
192 208 : DO ispin = 1, nspin
193 106 : 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 106 : 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 106 : matrix_c=fm_work2)
206 106 : IF (print_level > 2) THEN
207 : ! Write the full Lowdin population matrix
208 4 : IF (nspin > 1) THEN
209 4 : IF (ispin == 1) THEN
210 2 : fm_work2%name = TRIM(headline)//" FOR ALPHA SPIN"
211 : ELSE
212 2 : 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 4 : output_unit=output_unit)
217 : END IF
218 208 : CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
219 : END DO ! next spin ispin
220 :
221 : ! Write atomic populations and charges
222 102 : IF (output_unit > 0) THEN
223 51 : print_gop = (print_level > 1) ! Print also orbital populations
224 51 : 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 102 : CALL cp_fm_release(matrix=fm_s_half)
229 102 : CALL cp_fm_release(matrix=fm_work1)
230 102 : CALL cp_fm_release(matrix=fm_work2)
231 306 : IF (ASSOCIATED(orbpop)) THEN
232 102 : DEALLOCATE (orbpop)
233 : END IF
234 :
235 : END IF
236 :
237 102 : CALL timestop(handle)
238 :
239 102 : 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 4780 : 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 :: handle, iatom, ic, isgf, ispin, jatom, &
260 : jsgf, natom, nsgf, nspin, sgfa, sgfb
261 4780 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
262 : LOGICAL :: found, print_gop
263 : REAL(KIND=dp) :: ps
264 4780 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop, p_block, ps_block, s_block
265 4780 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
266 : TYPE(dbcsr_iterator_type) :: iter
267 4780 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
268 : TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s
269 : TYPE(mp_para_env_type), POINTER :: para_env
270 4780 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
271 4780 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
272 : TYPE(qs_rho_type), POINTER :: rho
273 :
274 4780 : CALL timeset(routineN, handle)
275 :
276 4780 : NULLIFY (atomic_kind_set)
277 4780 : NULLIFY (qs_kind_set)
278 4780 : NULLIFY (matrix_p)
279 4780 : NULLIFY (matrix_s)
280 4780 : NULLIFY (orbpop)
281 4780 : NULLIFY (particle_set)
282 4780 : NULLIFY (ps_block)
283 4780 : NULLIFY (p_block)
284 4780 : NULLIFY (rho)
285 4780 : NULLIFY (sm_p)
286 4780 : NULLIFY (sm_ps)
287 4780 : NULLIFY (sm_s)
288 4780 : NULLIFY (s_block)
289 4780 : NULLIFY (para_env)
290 :
291 : CALL get_qs_env(qs_env=qs_env, &
292 : atomic_kind_set=atomic_kind_set, &
293 : qs_kind_set=qs_kind_set, &
294 : matrix_s_kp=matrix_s, &
295 : particle_set=particle_set, &
296 : rho=rho, &
297 4780 : para_env=para_env)
298 :
299 4780 : CPASSERT(ASSOCIATED(atomic_kind_set))
300 4780 : CPASSERT(ASSOCIATED(qs_kind_set))
301 4780 : CPASSERT(ASSOCIATED(particle_set))
302 4780 : CPASSERT(ASSOCIATED(rho))
303 4780 : CPASSERT(ASSOCIATED(matrix_s))
304 :
305 4780 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
306 4780 : nspin = SIZE(matrix_p, 1)
307 :
308 : ! Get the total number of contracted spherical Gaussian basis functions
309 4780 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
310 4780 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
311 14340 : ALLOCATE (first_sgf_atom(natom))
312 25088 : first_sgf_atom(:) = 0
313 :
314 4780 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
315 :
316 : ! Provide an array to store the orbital populations for each spin
317 19120 : ALLOCATE (orbpop(nsgf, nspin))
318 168140 : orbpop(:, :) = 0.0_dp
319 :
320 : ! Write headline
321 4780 : IF (output_unit > 0) THEN
322 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
323 2404 : '!-----------------------------------------------------------------------------!'
324 2404 : WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
325 : END IF
326 :
327 : ! Create a DBCSR work matrix, if needed
328 4780 : IF (print_level > 2) THEN
329 2 : sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
330 2 : ALLOCATE (sm_ps)
331 2 : headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
332 2 : IF (nspin > 1) THEN
333 2 : IF (ispin == 1) THEN
334 0 : headline = TRIM(headline)//" For Alpha Spin"
335 : ELSE
336 2 : headline = TRIM(headline)//" For Beta Spin"
337 : END IF
338 : END IF
339 2 : CALL dbcsr_copy(matrix_b=sm_ps, matrix_a=sm_s, name=TRIM(headline))
340 : END IF
341 :
342 : ! Build Mulliken population matrix for each spin
343 10316 : DO ispin = 1, nspin
344 15750 : DO ic = 1, SIZE(matrix_s, 2)
345 10214 : IF (print_level > 2) THEN
346 4 : CALL dbcsr_set(sm_ps, 0.0_dp)
347 : END IF
348 10214 : sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
349 10214 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
350 : ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
351 : ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
352 10214 : CALL dbcsr_iterator_start(iter, sm_s)
353 97929 : DO WHILE (dbcsr_iterator_blocks_left(iter))
354 87715 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block)
355 87715 : IF (.NOT. (ASSOCIATED(s_block))) CYCLE
356 : CALL dbcsr_get_block_p(matrix=sm_p, &
357 : row=iatom, &
358 : col=jatom, &
359 : block=p_block, &
360 87715 : found=found)
361 87715 : IF (print_level > 2) THEN
362 : CALL dbcsr_get_block_p(matrix=sm_ps, &
363 : row=iatom, &
364 : col=jatom, &
365 : block=ps_block, &
366 12 : found=found)
367 12 : CPASSERT(ASSOCIATED(ps_block))
368 : END IF
369 :
370 87715 : sgfb = first_sgf_atom(jatom)
371 478968 : DO jsgf = 1, SIZE(s_block, 2)
372 5012028 : DO isgf = 1, SIZE(s_block, 1)
373 4620775 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
374 4620775 : IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
375 5012028 : orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
376 : END DO
377 478968 : sgfb = sgfb + 1
378 : END DO
379 97929 : IF (iatom /= jatom) THEN
380 70663 : sgfa = first_sgf_atom(iatom)
381 375310 : DO isgf = 1, SIZE(s_block, 1)
382 3292599 : DO jsgf = 1, SIZE(s_block, 2)
383 2987952 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
384 3292599 : orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
385 : END DO
386 375310 : sgfa = sgfa + 1
387 : END DO
388 : END IF
389 : END DO
390 25964 : CALL dbcsr_iterator_stop(iter)
391 : END DO
392 :
393 10316 : IF (print_level > 2) THEN
394 : ! Write the full Mulliken net AO and overlap population matrix
395 4 : CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, output_unit=output_unit)
396 : END IF
397 : END DO
398 :
399 331500 : CALL para_env%sum(orbpop)
400 :
401 : ! Write atomic populations and charges
402 4780 : IF (output_unit > 0) THEN
403 2404 : print_gop = (print_level > 1) ! Print also orbital populations
404 2404 : CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
405 : END IF
406 :
407 : ! Save the Mulliken charges to results
408 4780 : CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
409 :
410 : ! Release local working storage
411 4780 : IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
412 4780 : IF (ASSOCIATED(orbpop)) THEN
413 4780 : DEALLOCATE (orbpop)
414 : END IF
415 4780 : IF (ALLOCATED(first_sgf_atom)) THEN
416 4780 : DEALLOCATE (first_sgf_atom)
417 : END IF
418 :
419 4780 : IF (output_unit > 0) THEN
420 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
421 2404 : '!-----------------------------------------------------------------------------!'
422 : END IF
423 :
424 4780 : CALL timestop(handle)
425 :
426 14340 : END SUBROUTINE mulliken_population_analysis
427 :
428 : ! **************************************************************************************************
429 : !> \brief Save the Mulliken atomic orbital populations and charges in results
430 : !>
431 : !> \param orbpop ...
432 : !> \param atomic_kind_set ...
433 : !> \param qs_kind_set ...
434 : !> \param particle_set ...
435 : !> \param qs_env ...
436 : !> \par History
437 : !> 27.05.2022 BT
438 : !> 16.07.2025 RK
439 : !> \author Bo Thomsen (BT)
440 : !> Rangsiman Ketkaew (RK)
441 : !> \version 1.0
442 : ! **************************************************************************************************
443 4780 : SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
444 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
445 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
446 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
447 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
448 : TYPE(qs_environment_type), POINTER :: qs_env
449 :
450 : CHARACTER(LEN=default_string_length) :: description
451 : INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
452 : iso, l, natom, nset, nsgf, nspin
453 4780 : INTEGER, DIMENSION(:), POINTER :: nshell
454 4780 : INTEGER, DIMENSION(:, :), POINTER :: lshell
455 : REAL(KIND=dp) :: zeff
456 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop
457 : REAL(KIND=dp), DIMENSION(:), POINTER :: all_sumorbpop, charges_save
458 : TYPE(cp_result_type), POINTER :: results
459 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
460 :
461 4780 : NULLIFY (lshell)
462 4780 : NULLIFY (nshell)
463 4780 : NULLIFY (orb_basis_set)
464 :
465 0 : CPASSERT(ASSOCIATED(orbpop))
466 4780 : CPASSERT(ASSOCIATED(atomic_kind_set))
467 4780 : CPASSERT(ASSOCIATED(particle_set))
468 :
469 4780 : nspin = SIZE(orbpop, 2)
470 :
471 4780 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
472 4780 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
473 4780 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
474 4780 : NULLIFY (results)
475 4780 : CALL get_qs_env(qs_env, results=results)
476 14340 : ALLOCATE (all_sumorbpop(natom))
477 9560 : ALLOCATE (charges_save(natom))
478 :
479 4780 : iao = 1
480 25088 : DO iatom = 1, natom
481 20308 : sumorbpop(:) = 0.0_dp
482 20308 : NULLIFY (orb_basis_set)
483 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
484 20308 : kind_number=ikind)
485 20308 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
486 25088 : 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 20308 : l=lshell)
491 20308 : isgf = 1
492 55142 : DO iset = 1, nset
493 118122 : DO ishell = 1, nshell(iset)
494 62980 : l = lshell(ishell, iset)
495 231854 : DO iso = 1, nso(l)
496 134040 : IF (nspin == 1) THEN
497 110256 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
498 : ELSE
499 71352 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
500 23784 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
501 : END IF
502 134040 : isgf = isgf + 1
503 197020 : iao = iao + 1
504 : END DO
505 : END DO
506 : END DO
507 20308 : IF (nspin == 1) THEN
508 17436 : charges_save(iatom) = zeff - sumorbpop(1)
509 17436 : all_sumorbpop(iatom) = sumorbpop(1)
510 : ELSE
511 2872 : charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
512 2872 : all_sumorbpop(iatom) = sumorbpop(1) + sumorbpop(2)
513 : END IF
514 : END IF ! atom has an orbital basis
515 : END DO ! next atom iatom
516 :
517 : ! Store atomic orbital populations in results
518 4780 : description = "[MULLIKEN-ORBPOP]"
519 4780 : CALL cp_results_erase(results=results, description=description)
520 : CALL put_results(results=results, description=description, &
521 4780 : values=orbpop)
522 :
523 : ! Store sum orbital population in results
524 4780 : description = "[MULLIKEN-SUMORBPOP]"
525 4780 : CALL cp_results_erase(results=results, description=description)
526 : CALL put_results(results=results, description=description, &
527 4780 : values=all_sumorbpop)
528 :
529 : ! Store charges in results
530 4780 : description = "[MULLIKEN-CHARGES]"
531 4780 : CALL cp_results_erase(results=results, description=description)
532 : CALL put_results(results=results, description=description, &
533 4780 : values=charges_save)
534 :
535 4780 : DEALLOCATE (all_sumorbpop)
536 4780 : DEALLOCATE (charges_save)
537 :
538 9560 : END SUBROUTINE save_mulliken_charges
539 :
540 : ! **************************************************************************************************
541 : !> \brief Write atomic orbital populations and net atomic charges
542 : !>
543 : !> \param orbpop ...
544 : !> \param atomic_kind_set ...
545 : !> \param qs_kind_set ...
546 : !> \param particle_set ...
547 : !> \param output_unit ...
548 : !> \param print_orbital_contributions ...
549 : !> \date 07.07.2010
550 : !> \author Matthias Krack (MK)
551 : !> \version 1.0
552 : ! **************************************************************************************************
553 7365 : SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
554 : print_orbital_contributions)
555 :
556 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
557 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
558 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
559 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
560 : INTEGER, INTENT(IN) :: output_unit
561 : LOGICAL, INTENT(IN) :: print_orbital_contributions
562 :
563 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_orbpop'
564 :
565 : CHARACTER(LEN=2) :: element_symbol
566 2455 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
567 : INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
568 : ishell, iso, l, natom, nset, nsgf, &
569 : nspin
570 2455 : INTEGER, DIMENSION(:), POINTER :: nshell
571 2455 : INTEGER, DIMENSION(:, :), POINTER :: lshell
572 : REAL(KIND=dp) :: zeff
573 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
574 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
575 :
576 2455 : CALL timeset(routineN, handle)
577 :
578 2455 : NULLIFY (lshell)
579 2455 : NULLIFY (nshell)
580 2455 : NULLIFY (orb_basis_set)
581 2455 : NULLIFY (sgf_symbol)
582 :
583 2455 : CPASSERT(ASSOCIATED(orbpop))
584 2455 : CPASSERT(ASSOCIATED(atomic_kind_set))
585 2455 : CPASSERT(ASSOCIATED(particle_set))
586 :
587 2455 : nspin = SIZE(orbpop, 2)
588 :
589 2455 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
590 2455 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
591 :
592 : ! Select and write headline
593 2455 : IF (nspin == 1) THEN
594 2075 : IF (print_orbital_contributions) THEN
595 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
596 0 : "# Orbital AO symbol Orbital population Net charge"
597 : ELSE
598 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
599 2075 : "# Atom Element Kind Atomic population Net charge"
600 : END IF
601 : ELSE
602 380 : IF (print_orbital_contributions) THEN
603 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
604 2 : "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
605 : ELSE
606 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
607 378 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
608 : END IF
609 : END IF
610 :
611 2455 : totsumorbpop(:) = 0.0_dp
612 :
613 2455 : iao = 1
614 12791 : DO iatom = 1, natom
615 10336 : sumorbpop(:) = 0.0_dp
616 10336 : NULLIFY (orb_basis_set)
617 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
618 : element_symbol=element_symbol, &
619 10336 : kind_number=ikind)
620 10336 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
621 12791 : IF (ASSOCIATED(orb_basis_set)) THEN
622 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
623 : nset=nset, &
624 : nshell=nshell, &
625 : l=lshell, &
626 10336 : sgf_symbol=sgf_symbol)
627 10336 : isgf = 1
628 28299 : DO iset = 1, nset
629 60465 : DO ishell = 1, nshell(iset)
630 32166 : l = lshell(ishell, iset)
631 118643 : DO iso = 1, nso(l)
632 68514 : IF (nspin == 1) THEN
633 56557 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
634 56557 : IF (print_orbital_contributions) THEN
635 0 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
636 : WRITE (UNIT=output_unit, &
637 : FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
638 0 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
639 : END IF
640 : ELSE
641 35871 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
642 11957 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
643 11957 : IF (print_orbital_contributions) THEN
644 78 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
645 : WRITE (UNIT=output_unit, &
646 : FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
647 78 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
648 156 : orbpop(iao, 1) - orbpop(iao, 2)
649 : END IF
650 : END IF
651 68514 : isgf = isgf + 1
652 100680 : iao = iao + 1
653 : END DO
654 : END DO
655 : END DO
656 10336 : IF (nspin == 1) THEN
657 8895 : totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
658 8895 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
659 : WRITE (UNIT=output_unit, &
660 : FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
661 8895 : iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
662 : ELSE
663 4323 : totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
664 1441 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
665 : WRITE (UNIT=output_unit, &
666 : FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
667 1441 : iatom, element_symbol, ikind, sumorbpop(1:2), &
668 2882 : zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
669 : END IF
670 : END IF ! atom has an orbital basis
671 : END DO ! next atom iatom
672 :
673 : ! Write total sums
674 2455 : IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
675 2455 : IF (nspin == 1) THEN
676 : WRITE (UNIT=output_unit, &
677 : FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
678 2075 : "# Total charge", totsumorbpop(1), totsumorbpop(3)
679 : ELSE
680 : WRITE (UNIT=output_unit, &
681 : FMT="(T2,A,T28,4(1X,F12.6),/)") &
682 380 : "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
683 760 : totsumorbpop(1) - totsumorbpop(2)
684 : END IF
685 :
686 2455 : IF (output_unit > 0) CALL m_flush(output_unit)
687 :
688 2455 : CALL timestop(handle)
689 :
690 2455 : END SUBROUTINE write_orbpop
691 :
692 : END MODULE population_analyses
|
