Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Provide various population analyses and print the requested output
9 : !> information
10 : !>
11 : !> \author Matthias Krack (MK)
12 : !> \date 09.07.2010
13 : !> \version 1.0
14 : ! **************************************************************************************************
15 :
16 : MODULE population_analyses
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind,&
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: get_gto_basis_set,&
21 : gto_basis_set_type
22 : USE cp_blacs_env, ONLY: cp_blacs_env_type
23 : USE cp_dbcsr_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_setname, 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 4568 : SUBROUTINE mulliken_population_analysis(qs_env, output_unit, print_level)
252 :
253 : TYPE(qs_environment_type), POINTER :: qs_env
254 : INTEGER, INTENT(IN) :: output_unit, print_level
255 :
256 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_population_analysis'
257 :
258 : CHARACTER(LEN=default_string_length) :: headline
259 : INTEGER :: blk, handle, iatom, ic, isgf, ispin, &
260 : jatom, jsgf, natom, nsgf, nspin, sgfa, &
261 : sgfb
262 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom
263 : LOGICAL :: found, print_gop
264 : REAL(KIND=dp) :: ps
265 4568 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop, p_block, ps_block, s_block
266 4568 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
267 : TYPE(dbcsr_iterator_type) :: iter
268 4568 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
269 : TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s
270 : TYPE(mp_para_env_type), POINTER :: para_env
271 4568 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
272 4568 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
273 : TYPE(qs_rho_type), POINTER :: rho
274 :
275 4568 : CALL timeset(routineN, handle)
276 :
277 4568 : NULLIFY (atomic_kind_set)
278 4568 : NULLIFY (qs_kind_set)
279 4568 : NULLIFY (matrix_p)
280 4568 : NULLIFY (matrix_s)
281 : NULLIFY (orbpop)
282 4568 : NULLIFY (particle_set)
283 4568 : NULLIFY (ps_block)
284 4568 : NULLIFY (p_block)
285 4568 : NULLIFY (rho)
286 4568 : NULLIFY (sm_p)
287 4568 : NULLIFY (sm_ps)
288 4568 : NULLIFY (sm_s)
289 4568 : NULLIFY (s_block)
290 4568 : NULLIFY (para_env)
291 :
292 : CALL get_qs_env(qs_env=qs_env, &
293 : atomic_kind_set=atomic_kind_set, &
294 : qs_kind_set=qs_kind_set, &
295 : matrix_s_kp=matrix_s, &
296 : particle_set=particle_set, &
297 : rho=rho, &
298 4568 : para_env=para_env)
299 :
300 4568 : CPASSERT(ASSOCIATED(atomic_kind_set))
301 4568 : CPASSERT(ASSOCIATED(qs_kind_set))
302 4568 : CPASSERT(ASSOCIATED(particle_set))
303 4568 : CPASSERT(ASSOCIATED(rho))
304 4568 : CPASSERT(ASSOCIATED(matrix_s))
305 :
306 4568 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
307 4568 : nspin = SIZE(matrix_p, 1)
308 :
309 : ! Get the total number of contracted spherical Gaussian basis functions
310 4568 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
311 4568 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
312 13704 : ALLOCATE (first_sgf_atom(natom))
313 24164 : first_sgf_atom(:) = 0
314 :
315 4568 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
316 :
317 : ! Provide an array to store the orbital populations for each spin
318 18272 : ALLOCATE (orbpop(nsgf, nspin))
319 163128 : orbpop(:, :) = 0.0_dp
320 :
321 : ! Write headline
322 4568 : IF (output_unit > 0) THEN
323 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
324 2298 : '!-----------------------------------------------------------------------------!'
325 2298 : WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
326 : END IF
327 :
328 : ! Create a DBCSR work matrix, if needed
329 4568 : IF (print_level > 2) THEN
330 2 : sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
331 2 : ALLOCATE (sm_ps)
332 2 : headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
333 : CALL dbcsr_copy(matrix_b=sm_ps, &
334 : matrix_a=sm_s, &
335 2 : name=TRIM(headline))
336 : END IF
337 :
338 : ! Build Mulliken population matrix for each spin
339 9886 : DO ispin = 1, nspin
340 14690 : DO ic = 1, SIZE(matrix_s, 2)
341 9372 : IF (print_level > 2) THEN
342 4 : CALL dbcsr_set(sm_ps, 0.0_dp)
343 : END IF
344 9372 : sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
345 9372 : sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
346 : ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
347 : ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
348 9372 : CALL dbcsr_iterator_start(iter, sm_s)
349 94964 : DO WHILE (dbcsr_iterator_blocks_left(iter))
350 85592 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block, blk)
351 85592 : IF (.NOT. (ASSOCIATED(s_block))) CYCLE
352 : CALL dbcsr_get_block_p(matrix=sm_p, &
353 : row=iatom, &
354 : col=jatom, &
355 : block=p_block, &
356 85592 : found=found)
357 85592 : IF (print_level > 2) THEN
358 : CALL dbcsr_get_block_p(matrix=sm_ps, &
359 : row=iatom, &
360 : col=jatom, &
361 : block=ps_block, &
362 12 : found=found)
363 12 : CPASSERT(ASSOCIATED(ps_block))
364 : END IF
365 :
366 85592 : sgfb = first_sgf_atom(jatom)
367 463100 : DO jsgf = 1, SIZE(s_block, 2)
368 4848141 : DO isgf = 1, SIZE(s_block, 1)
369 4470633 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
370 4470633 : IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
371 4848141 : orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
372 : END DO
373 463100 : sgfb = sgfb + 1
374 : END DO
375 94964 : IF (iatom /= jatom) THEN
376 69666 : sgfa = first_sgf_atom(iatom)
377 365022 : DO isgf = 1, SIZE(s_block, 1)
378 3224583 : DO jsgf = 1, SIZE(s_block, 2)
379 2929227 : ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
380 3224583 : orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
381 : END DO
382 365022 : sgfa = sgfa + 1
383 : END DO
384 : END IF
385 : END DO
386 24062 : CALL dbcsr_iterator_stop(iter)
387 : END DO
388 :
389 9886 : IF (print_level > 2) THEN
390 : ! Write the full Mulliken net AO and overlap population matrix
391 4 : IF (nspin > 1) THEN
392 4 : IF (ispin == 1) THEN
393 2 : CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Alpha Spin")
394 : ELSE
395 2 : CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Beta Spin")
396 : END IF
397 : END IF
398 : CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, &
399 4 : output_unit=output_unit)
400 : END IF
401 : END DO
402 :
403 321688 : CALL para_env%sum(orbpop)
404 :
405 : ! Write atomic populations and charges
406 4568 : IF (output_unit > 0) THEN
407 2298 : print_gop = (print_level > 1) ! Print also orbital populations
408 2298 : CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
409 : END IF
410 :
411 : ! Save the Mulliken charges to results
412 4568 : CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
413 :
414 : ! Release local working storage
415 4568 : IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
416 4568 : IF (ASSOCIATED(orbpop)) THEN
417 4568 : DEALLOCATE (orbpop)
418 : END IF
419 4568 : IF (ALLOCATED(first_sgf_atom)) THEN
420 4568 : DEALLOCATE (first_sgf_atom)
421 : END IF
422 :
423 4568 : IF (output_unit > 0) THEN
424 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
425 2298 : '!-----------------------------------------------------------------------------!'
426 : END IF
427 :
428 4568 : CALL timestop(handle)
429 :
430 13704 : END SUBROUTINE mulliken_population_analysis
431 :
432 : ! **************************************************************************************************
433 : !> \brief Save the Mulliken charges in results
434 : !>
435 : !> \param orbpop ...
436 : !> \param atomic_kind_set ...
437 : !> \param qs_kind_set ...
438 : !> \param particle_set ...
439 : !> \param qs_env ...
440 : !> \date 27.05.2022
441 : !> \author Bo Thomsen (BT)
442 : !> \version 1.0
443 : ! **************************************************************************************************
444 4568 : SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
445 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
446 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
447 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
448 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
449 : TYPE(qs_environment_type), POINTER :: qs_env
450 :
451 : CHARACTER(LEN=default_string_length) :: description
452 : INTEGER :: iao, iatom, ikind, iset, isgf, ishell, &
453 : iso, l, natom, nset, nsgf, nspin
454 4568 : INTEGER, DIMENSION(:), POINTER :: nshell
455 4568 : INTEGER, DIMENSION(:, :), POINTER :: lshell
456 : REAL(KIND=dp) :: zeff
457 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop
458 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_save
459 : TYPE(cp_result_type), POINTER :: results
460 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
461 :
462 4568 : NULLIFY (lshell)
463 4568 : NULLIFY (nshell)
464 4568 : NULLIFY (orb_basis_set)
465 :
466 0 : CPASSERT(ASSOCIATED(orbpop))
467 4568 : CPASSERT(ASSOCIATED(atomic_kind_set))
468 4568 : CPASSERT(ASSOCIATED(particle_set))
469 :
470 4568 : nspin = SIZE(orbpop, 2)
471 :
472 4568 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
473 4568 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
474 4568 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
475 4568 : NULLIFY (results)
476 4568 : CALL get_qs_env(qs_env, results=results)
477 13704 : ALLOCATE (charges_save(natom))
478 :
479 4568 : iao = 1
480 24164 : DO iatom = 1, natom
481 19596 : sumorbpop(:) = 0.0_dp
482 19596 : NULLIFY (orb_basis_set)
483 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
484 19596 : kind_number=ikind)
485 19596 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
486 24164 : 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 19596 : l=lshell)
491 19596 : isgf = 1
492 53182 : DO iset = 1, nset
493 113916 : DO ishell = 1, nshell(iset)
494 60734 : l = lshell(ishell, iset)
495 223870 : DO iso = 1, nso(l)
496 129550 : IF (nspin == 1) THEN
497 105858 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
498 : ELSE
499 71076 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
500 23692 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
501 : END IF
502 129550 : isgf = isgf + 1
503 190284 : iao = iao + 1
504 : END DO
505 : END DO
506 : END DO
507 19596 : IF (nspin == 1) THEN
508 16736 : charges_save(iatom) = zeff - sumorbpop(1)
509 : ELSE
510 2860 : charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
511 : END IF
512 : END IF ! atom has an orbital basis
513 : END DO ! next atom iatom
514 :
515 : ! Store charges in results
516 4568 : description = "[MULLIKEN-CHARGES]"
517 4568 : CALL cp_results_erase(results=results, description=description)
518 : CALL put_results(results=results, description=description, &
519 4568 : values=charges_save)
520 :
521 4568 : DEALLOCATE (charges_save)
522 :
523 9136 : END SUBROUTINE save_mulliken_charges
524 :
525 : ! **************************************************************************************************
526 : !> \brief Write atomic orbital populations and net atomic charges
527 : !>
528 : !> \param orbpop ...
529 : !> \param atomic_kind_set ...
530 : !> \param qs_kind_set ...
531 : !> \param particle_set ...
532 : !> \param output_unit ...
533 : !> \param print_orbital_contributions ...
534 : !> \date 07.07.2010
535 : !> \author Matthias Krack (MK)
536 : !> \version 1.0
537 : ! **************************************************************************************************
538 7017 : SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
539 : print_orbital_contributions)
540 :
541 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop
542 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
543 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
544 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
545 : INTEGER, INTENT(IN) :: output_unit
546 : LOGICAL, INTENT(IN) :: print_orbital_contributions
547 :
548 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_orbpop'
549 :
550 : CHARACTER(LEN=2) :: element_symbol
551 2339 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol
552 : INTEGER :: handle, iao, iatom, ikind, iset, isgf, &
553 : ishell, iso, l, natom, nset, nsgf, &
554 : nspin
555 2339 : INTEGER, DIMENSION(:), POINTER :: nshell
556 2339 : INTEGER, DIMENSION(:, :), POINTER :: lshell
557 : REAL(KIND=dp) :: zeff
558 : REAL(KIND=dp), DIMENSION(3) :: sumorbpop, totsumorbpop
559 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
560 :
561 2339 : CALL timeset(routineN, handle)
562 :
563 2339 : NULLIFY (lshell)
564 2339 : NULLIFY (nshell)
565 2339 : NULLIFY (orb_basis_set)
566 2339 : NULLIFY (sgf_symbol)
567 :
568 2339 : CPASSERT(ASSOCIATED(orbpop))
569 2339 : CPASSERT(ASSOCIATED(atomic_kind_set))
570 2339 : CPASSERT(ASSOCIATED(particle_set))
571 :
572 2339 : nspin = SIZE(orbpop, 2)
573 :
574 2339 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
575 2339 : CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
576 :
577 : ! Select and write headline
578 2339 : IF (nspin == 1) THEN
579 1962 : IF (print_orbital_contributions) THEN
580 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
581 0 : "# Orbital AO symbol Orbital population Net charge"
582 : ELSE
583 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
584 1962 : "# Atom Element Kind Atomic population Net charge"
585 : END IF
586 : ELSE
587 377 : IF (print_orbital_contributions) THEN
588 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
589 2 : "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment"
590 : ELSE
591 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
592 375 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
593 : END IF
594 : END IF
595 :
596 2339 : totsumorbpop(:) = 0.0_dp
597 :
598 2339 : iao = 1
599 12276 : DO iatom = 1, natom
600 9937 : sumorbpop(:) = 0.0_dp
601 9937 : NULLIFY (orb_basis_set)
602 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
603 : element_symbol=element_symbol, &
604 9937 : kind_number=ikind)
605 9937 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
606 12276 : IF (ASSOCIATED(orb_basis_set)) THEN
607 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
608 : nset=nset, &
609 : nshell=nshell, &
610 : l=lshell, &
611 9937 : sgf_symbol=sgf_symbol)
612 9937 : isgf = 1
613 27181 : DO iset = 1, nset
614 58129 : DO ishell = 1, nshell(iset)
615 30948 : l = lshell(ishell, iset)
616 114264 : DO iso = 1, nso(l)
617 66072 : IF (nspin == 1) THEN
618 54161 : sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
619 54161 : IF (print_orbital_contributions) THEN
620 0 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
621 : WRITE (UNIT=output_unit, &
622 : FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
623 0 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
624 : END IF
625 : ELSE
626 35733 : sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
627 11911 : sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
628 11911 : IF (print_orbital_contributions) THEN
629 78 : IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
630 : WRITE (UNIT=output_unit, &
631 : FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
632 78 : iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
633 156 : orbpop(iao, 1) - orbpop(iao, 2)
634 : END IF
635 : END IF
636 66072 : isgf = isgf + 1
637 97020 : iao = iao + 1
638 : END DO
639 : END DO
640 : END DO
641 9937 : IF (nspin == 1) THEN
642 8502 : totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
643 8502 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
644 : WRITE (UNIT=output_unit, &
645 : FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
646 8502 : iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
647 : ELSE
648 4305 : totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
649 1435 : totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
650 : WRITE (UNIT=output_unit, &
651 : FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
652 1435 : iatom, element_symbol, ikind, sumorbpop(1:2), &
653 2870 : zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
654 : END IF
655 : END IF ! atom has an orbital basis
656 : END DO ! next atom iatom
657 :
658 : ! Write total sums
659 2339 : IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
660 2339 : IF (nspin == 1) THEN
661 : WRITE (UNIT=output_unit, &
662 : FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
663 1962 : "# Total charge", totsumorbpop(1), totsumorbpop(3)
664 : ELSE
665 : WRITE (UNIT=output_unit, &
666 : FMT="(T2,A,T28,4(1X,F12.6),/)") &
667 377 : "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
668 754 : totsumorbpop(1) - totsumorbpop(2)
669 : END IF
670 :
671 2339 : IF (output_unit > 0) CALL m_flush(output_unit)
672 :
673 2339 : CALL timestop(handle)
674 :
675 2339 : END SUBROUTINE write_orbpop
676 :
677 : END MODULE population_analyses
|
