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 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 82 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
93 82 : 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 82 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
98 82 : 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 82 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
102 82 : 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 82 : CALL timeset(routineN, handle)
107 :
108 82 : NULLIFY (atomic_kind_set)
109 82 : NULLIFY (qs_kind_set)
110 82 : NULLIFY (fmstruct)
111 82 : NULLIFY (matrix_p)
112 82 : NULLIFY (matrixkp_p)
113 82 : NULLIFY (matrixkp_s)
114 82 : NULLIFY (orbpop)
115 82 : NULLIFY (particle_set)
116 82 : NULLIFY (rho)
117 82 : NULLIFY (scf_control)
118 82 : NULLIFY (sm_p)
119 82 : NULLIFY (sm_s)
120 : NULLIFY (orbpop)
121 82 : NULLIFY (para_env)
122 82 : 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 82 : blacs_env=blacs_env)
133 :
134 82 : CPASSERT(ASSOCIATED(atomic_kind_set))
135 82 : CPASSERT(ASSOCIATED(qs_kind_set))
136 82 : CPASSERT(ASSOCIATED(matrixkp_s))
137 82 : CPASSERT(ASSOCIATED(particle_set))
138 82 : CPASSERT(ASSOCIATED(rho))
139 82 : CPASSERT(ASSOCIATED(scf_control))
140 :
141 82 : IF (SIZE(matrixkp_s, 2) > 1) THEN
142 :
143 0 : CPWARN("Lowdin population analysis not implemented for k-points.")
144 :
145 : ELSE
146 :
147 82 : sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
148 82 : CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
149 :
150 82 : matrix_p => matrixkp_p(:, 1)
151 82 : nspin = SIZE(matrix_p, 1)
152 :
153 : ! Get the total number of contracted spherical Gaussian basis functions
154 82 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
155 :
156 : ! Provide an array to store the orbital populations for each spin
157 328 : ALLOCATE (orbpop(nsgf, nspin))
158 2612 : orbpop(:, :) = 0.0_dp
159 :
160 : ! Write headline
161 82 : IF (output_unit > 0) THEN
162 41 : 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 82 : ncol_global=nsgf)
171 : CALL cp_fm_create(matrix=fm_s_half, &
172 : matrix_struct=fmstruct, &
173 82 : name="S^(1/2) MATRIX")
174 : CALL cp_fm_create(matrix=fm_work1, &
175 : matrix_struct=fmstruct, &
176 82 : name="FULL WORK MATRIX 1")
177 82 : headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
178 : CALL cp_fm_create(matrix=fm_work2, &
179 : matrix_struct=fmstruct, &
180 82 : name=TRIM(headline))
181 82 : CALL cp_fm_struct_release(fmstruct=fmstruct)
182 :
183 : ! Build full S^(1/2) matrix (computationally expensive)
184 82 : CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
185 82 : CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
186 82 : 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 168 : DO ispin = 1, nspin
193 86 : 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 86 : 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 86 : matrix_c=fm_work2)
206 86 : 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 168 : CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
219 : END DO ! next spin ispin
220 :
221 : ! Write atomic populations and charges
222 82 : IF (output_unit > 0) THEN
223 41 : print_gop = (print_level > 1) ! Print also orbital populations
224 41 : 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 82 : CALL cp_fm_release(matrix=fm_s_half)
229 82 : CALL cp_fm_release(matrix=fm_work1)
230 82 : CALL cp_fm_release(matrix=fm_work2)
231 246 : IF (ASSOCIATED(orbpop)) THEN
232 82 : DEALLOCATE (orbpop)
233 : END IF
234 :
235 : END IF
236 :
237 82 : CALL timestop(handle)
238 :
239 82 : 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 4582 : 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 4582 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
262 : LOGICAL :: found, print_gop
263 : REAL(KIND=dp) :: ps
264 4582 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop, p_block, ps_block, s_block
265 4582 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
266 : TYPE(dbcsr_iterator_type) :: iter
267 4582 : 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 4582 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
271 4582 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
272 : TYPE(qs_rho_type), POINTER :: rho
273 :
274 4582 : CALL timeset(routineN, handle)
275 :
276 4582 : NULLIFY (atomic_kind_set)
277 4582 : NULLIFY (qs_kind_set)
278 4582 : NULLIFY (matrix_p)
279 4582 : NULLIFY (matrix_s)
280 4582 : NULLIFY (orbpop)
281 4582 : NULLIFY (particle_set)
282 4582 : NULLIFY (ps_block)
283 4582 : NULLIFY (p_block)
284 4582 : NULLIFY (rho)
285 4582 : NULLIFY (sm_p)
286 4582 : NULLIFY (sm_ps)
287 4582 : NULLIFY (sm_s)
288 4582 : NULLIFY (s_block)
289 4582 : 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 4582 : para_env=para_env)
298 :
299 4582 : CPASSERT(ASSOCIATED(atomic_kind_set))
300 4582 : CPASSERT(ASSOCIATED(qs_kind_set))
301 4582 : CPASSERT(ASSOCIATED(particle_set))
302 4582 : CPASSERT(ASSOCIATED(rho))
303 4582 : CPASSERT(ASSOCIATED(matrix_s))
304 :
305 4582 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
306 4582 : nspin = SIZE(matrix_p, 1)
307 :
308 : ! Get the total number of contracted spherical Gaussian basis functions
309 4582 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
310 4582 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
311 13746 : ALLOCATE (first_sgf_atom(natom))
312 24216 : first_sgf_atom(:) = 0
313 :
314 4582 : 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 18328 : ALLOCATE (orbpop(nsgf, nspin))
318 163266 : orbpop(:, :) = 0.0_dp
319 :
320 : ! Write headline
321 4582 : IF (output_unit > 0) THEN
322 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
323 2305 : '!-----------------------------------------------------------------------------!'
324 2305 : WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
325 : END IF
326 :
327 : ! Create a DBCSR work matrix, if needed
328 4582 : 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 9914 : DO ispin = 1, nspin
344 14718 : DO ic = 1, SIZE(matrix_s, 2)
345 9386 : IF (print_level > 2) THEN
346 4 : CALL dbcsr_set(sm_ps, 0.0_dp)
347 : END IF
348 9386 : sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
349 9386 : 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 9386 : CALL dbcsr_iterator_start(iter, sm_s)
353 94909 : DO WHILE (dbcsr_iterator_blocks_left(iter))
354 85523 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block)
355 85523 : 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 85523 : found=found)
361 85523 : 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 85523 : sgfb = first_sgf_atom(jatom)
371 462006 : DO jsgf = 1, SIZE(s_block, 2)
372 4835135 : DO isgf = 1, SIZE(s_block, 1)
373 4458652 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
374 4458652 : IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
375 4835135 : orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
376 : END DO
377 462006 : sgfb = sgfb + 1
378 : END DO
379 94909 : IF (iatom /= jatom) THEN
380 69578 : sgfa = first_sgf_atom(iatom)
381 364046 : DO isgf = 1, SIZE(s_block, 1)
382 3211959 : DO jsgf = 1, SIZE(s_block, 2)
383 2917491 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
384 3211959 : orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
385 : END DO
386 364046 : sgfa = sgfa + 1
387 : END DO
388 : END IF
389 : END DO
390 24104 : CALL dbcsr_iterator_stop(iter)
391 : END DO
392 :
393 9914 : 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 321950 : CALL para_env%sum(orbpop)
400 :
401 : ! Write atomic populations and charges
402 4582 : IF (output_unit > 0) THEN
403 2305 : print_gop = (print_level > 1) ! Print also orbital populations
404 2305 : 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 4582 : CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
409 :
410 : ! Release local working storage
411 4582 : IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
412 4582 : IF (ASSOCIATED(orbpop)) THEN
413 4582 : DEALLOCATE (orbpop)
414 : END IF
415 4582 : IF (ALLOCATED(first_sgf_atom)) THEN
416 4582 : DEALLOCATE (first_sgf_atom)
417 : END IF
418 :
419 4582 : IF (output_unit > 0) THEN
420 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
421 2305 : '!-----------------------------------------------------------------------------!'
422 : END IF
423 :
424 4582 : CALL timestop(handle)
425 :
426 13746 : END SUBROUTINE mulliken_population_analysis
427 :
428 : ! **************************************************************************************************
429 : !> \brief Save the Mulliken 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 : !> \date 27.05.2022
437 : !> \author Bo Thomsen (BT)
438 : !> \version 1.0
439 : ! **************************************************************************************************
440 4582 : SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
441 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
442 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
443 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
444 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
445 : TYPE(qs_environment_type), POINTER :: qs_env
446 :
447 : CHARACTER(LEN=default_string_length) :: description
448 : INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
449 : iso, l, natom, nset, nsgf, nspin
450 4582 : INTEGER, DIMENSION(:), POINTER :: nshell
451 4582 : INTEGER, DIMENSION(:, :), POINTER :: lshell
452 : REAL(KIND=dp) :: zeff
453 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop
454 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_save
455 : TYPE(cp_result_type), POINTER :: results
456 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
457 :
458 4582 : NULLIFY (lshell)
459 4582 : NULLIFY (nshell)
460 4582 : NULLIFY (orb_basis_set)
461 :
462 0 : CPASSERT(ASSOCIATED(orbpop))
463 4582 : CPASSERT(ASSOCIATED(atomic_kind_set))
464 4582 : CPASSERT(ASSOCIATED(particle_set))
465 :
466 4582 : nspin = SIZE(orbpop, 2)
467 :
468 4582 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
469 4582 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
470 4582 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
471 4582 : NULLIFY (results)
472 4582 : CALL get_qs_env(qs_env, results=results)
473 13746 : ALLOCATE (charges_save(natom))
474 :
475 4582 : iao = 1
476 24216 : DO iatom = 1, natom
477 19634 : sumorbpop(:) = 0.0_dp
478 19634 : NULLIFY (orb_basis_set)
479 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
480 19634 : kind_number=ikind)
481 19634 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
482 24216 : IF (ASSOCIATED(orb_basis_set)) THEN
483 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
484 : nset=nset, &
485 : nshell=nshell, &
486 19634 : l=lshell)
487 19634 : isgf = 1
488 53238 : DO iset = 1, nset
489 114066 : DO ishell = 1, nshell(iset)
490 60828 : l = lshell(ishell, iset)
491 224092 : DO iso = 1, nso(l)
492 129660 : IF (nspin == 1) THEN
493 105968 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
494 : ELSE
495 71076 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
496 23692 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
497 : END IF
498 129660 : isgf = isgf + 1
499 190488 : iao = iao + 1
500 : END DO
501 : END DO
502 : END DO
503 19634 : IF (nspin == 1) THEN
504 16774 : charges_save(iatom) = zeff - sumorbpop(1)
505 : ELSE
506 2860 : charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
507 : END IF
508 : END IF ! atom has an orbital basis
509 : END DO ! next atom iatom
510 :
511 : ! Store charges in results
512 4582 : description = "[MULLIKEN-CHARGES]"
513 4582 : CALL cp_results_erase(results=results, description=description)
514 : CALL put_results(results=results, description=description, &
515 4582 : values=charges_save)
516 :
517 4582 : DEALLOCATE (charges_save)
518 :
519 9164 : END SUBROUTINE save_mulliken_charges
520 :
521 : ! **************************************************************************************************
522 : !> \brief Write atomic orbital populations and net atomic charges
523 : !>
524 : !> \param orbpop ...
525 : !> \param atomic_kind_set ...
526 : !> \param qs_kind_set ...
527 : !> \param particle_set ...
528 : !> \param output_unit ...
529 : !> \param print_orbital_contributions ...
530 : !> \date 07.07.2010
531 : !> \author Matthias Krack (MK)
532 : !> \version 1.0
533 : ! **************************************************************************************************
534 7038 : SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
535 : print_orbital_contributions)
536 :
537 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
538 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
539 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
540 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
541 : INTEGER, INTENT(IN) :: output_unit
542 : LOGICAL, INTENT(IN) :: print_orbital_contributions
543 :
544 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_orbpop'
545 :
546 : CHARACTER(LEN=2) :: element_symbol
547 2346 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
548 : INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
549 : ishell, iso, l, natom, nset, nsgf, &
550 : nspin
551 2346 : INTEGER, DIMENSION(:), POINTER :: nshell
552 2346 : INTEGER, DIMENSION(:, :), POINTER :: lshell
553 : REAL(KIND=dp) :: zeff
554 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
555 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
556 :
557 2346 : CALL timeset(routineN, handle)
558 :
559 2346 : NULLIFY (lshell)
560 2346 : NULLIFY (nshell)
561 2346 : NULLIFY (orb_basis_set)
562 2346 : NULLIFY (sgf_symbol)
563 :
564 2346 : CPASSERT(ASSOCIATED(orbpop))
565 2346 : CPASSERT(ASSOCIATED(atomic_kind_set))
566 2346 : CPASSERT(ASSOCIATED(particle_set))
567 :
568 2346 : nspin = SIZE(orbpop, 2)
569 :
570 2346 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
571 2346 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
572 :
573 : ! Select and write headline
574 2346 : IF (nspin == 1) THEN
575 1969 : IF (print_orbital_contributions) THEN
576 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
577 0 : "# Orbital AO symbol Orbital population Net charge"
578 : ELSE
579 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
580 1969 : "# Atom Element Kind Atomic population Net charge"
581 : END IF
582 : ELSE
583 377 : IF (print_orbital_contributions) THEN
584 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
585 2 : "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
586 : ELSE
587 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
588 375 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
589 : END IF
590 : END IF
591 :
592 2346 : totsumorbpop(:) = 0.0_dp
593 :
594 2346 : iao = 1
595 12302 : DO iatom = 1, natom
596 9956 : sumorbpop(:) = 0.0_dp
597 9956 : NULLIFY (orb_basis_set)
598 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
599 : element_symbol=element_symbol, &
600 9956 : kind_number=ikind)
601 9956 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
602 12302 : IF (ASSOCIATED(orb_basis_set)) THEN
603 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
604 : nset=nset, &
605 : nshell=nshell, &
606 : l=lshell, &
607 9956 : sgf_symbol=sgf_symbol)
608 9956 : isgf = 1
609 27209 : DO iset = 1, nset
610 58204 : DO ishell = 1, nshell(iset)
611 30995 : l = lshell(ishell, iset)
612 114375 : DO iso = 1, nso(l)
613 66127 : IF (nspin == 1) THEN
614 54216 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
615 54216 : IF (print_orbital_contributions) THEN
616 0 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
617 : WRITE (UNIT=output_unit, &
618 : FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
619 0 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
620 : END IF
621 : ELSE
622 35733 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
623 11911 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
624 11911 : IF (print_orbital_contributions) THEN
625 78 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
626 : WRITE (UNIT=output_unit, &
627 : FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
628 78 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
629 156 : orbpop(iao, 1) - orbpop(iao, 2)
630 : END IF
631 : END IF
632 66127 : isgf = isgf + 1
633 97122 : iao = iao + 1
634 : END DO
635 : END DO
636 : END DO
637 9956 : IF (nspin == 1) THEN
638 8521 : totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
639 8521 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
640 : WRITE (UNIT=output_unit, &
641 : FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
642 8521 : iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
643 : ELSE
644 4305 : totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
645 1435 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
646 : WRITE (UNIT=output_unit, &
647 : FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
648 1435 : iatom, element_symbol, ikind, sumorbpop(1:2), &
649 2870 : zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
650 : END IF
651 : END IF ! atom has an orbital basis
652 : END DO ! next atom iatom
653 :
654 : ! Write total sums
655 2346 : IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
656 2346 : IF (nspin == 1) THEN
657 : WRITE (UNIT=output_unit, &
658 : FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
659 1969 : "# Total charge", totsumorbpop(1), totsumorbpop(3)
660 : ELSE
661 : WRITE (UNIT=output_unit, &
662 : FMT="(T2,A,T28,4(1X,F12.6),/)") &
663 377 : "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
664 754 : totsumorbpop(1) - totsumorbpop(2)
665 : END IF
666 :
667 2346 : IF (output_unit > 0) CALL m_flush(output_unit)
668 :
669 2346 : CALL timestop(handle)
670 :
671 2346 : END SUBROUTINE write_orbpop
672 :
673 : END MODULE population_analyses
|
