1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6!> \brief Provide various population analyses and print the requested output 7!> information 8!> 9!> \author Matthias Krack (MK) 10!> \date 09.07.2010 11!> \version 1.0 12! ************************************************************************************************** 13 14MODULE population_analyses 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind,& 17 get_atomic_kind_set 18 USE basis_set_types, ONLY: get_gto_basis_set,& 19 gto_basis_set_type 20 USE cp_blacs_env, ONLY: cp_blacs_env_type 21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 22 cp_dbcsr_sm_fm_multiply 23 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,& 24 write_fm_with_basis_info 25 USE cp_fm_diag, ONLY: cp_fm_power 26 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 27 cp_fm_struct_release,& 28 cp_fm_struct_type 29 USE cp_fm_types, ONLY: cp_fm_create,& 30 cp_fm_get_diag,& 31 cp_fm_release,& 32 cp_fm_type 33 USE cp_gemm_interface, ONLY: cp_gemm 34 USE cp_para_types, ONLY: cp_para_env_type 35 USE dbcsr_api, ONLY: & 36 dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_iterator_blocks_left, & 37 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 38 dbcsr_p_type, dbcsr_set, dbcsr_setname, dbcsr_type 39 USE kinds, ONLY: default_string_length,& 40 dp 41 USE machine, ONLY: m_flush 42 USE message_passing, ONLY: mp_sum 43 USE orbital_pointers, ONLY: nso 44 USE particle_methods, ONLY: get_particle_set 45 USE particle_types, ONLY: particle_type 46 USE qs_environment_types, ONLY: get_qs_env,& 47 qs_environment_type 48 USE qs_kind_types, ONLY: get_qs_kind,& 49 get_qs_kind_set,& 50 qs_kind_type 51 USE qs_rho_types, ONLY: qs_rho_get,& 52 qs_rho_type 53 USE scf_control_types, ONLY: scf_control_type 54#include "./base/base_uses.f90" 55 56 IMPLICIT NONE 57 58 PRIVATE 59 60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'population_analyses' 61 62 PUBLIC :: lowdin_population_analysis, & 63 mulliken_population_analysis 64 65CONTAINS 66 67! ************************************************************************************************** 68!> \brief Perform a Lowdin population analysis based on a symmetric 69!> orthogonalisation of the density matrix using S^(1/2) 70!> 71!> \param qs_env ... 72!> \param output_unit ... 73!> \param print_level ... 74!> \date 06.07.2010 75!> \author Matthias Krack (MK) 76!> \version 1.0 77! ************************************************************************************************** 78 SUBROUTINE lowdin_population_analysis(qs_env, output_unit, print_level) 79 80 TYPE(qs_environment_type), POINTER :: qs_env 81 INTEGER, INTENT(IN) :: output_unit, print_level 82 83 CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin_population_analysis', & 84 routineP = moduleN//':'//routineN 85 86 CHARACTER(LEN=default_string_length) :: headline 87 INTEGER :: handle, ispin, ndep, nsgf, nspin 88 LOGICAL :: print_gop 89 REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop 90 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 91 TYPE(cp_blacs_env_type), POINTER :: blacs_env 92 TYPE(cp_fm_struct_type), POINTER :: fmstruct 93 TYPE(cp_fm_type), POINTER :: fm_s_half, fm_work1, fm_work2 94 TYPE(cp_para_env_type), POINTER :: para_env 95 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p 96 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_p, matrixkp_s 97 TYPE(dbcsr_type), POINTER :: sm_p, sm_s 98 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 99 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 100 TYPE(qs_rho_type), POINTER :: rho 101 TYPE(scf_control_type), POINTER :: scf_control 102 103 CALL timeset(routineN, handle) 104 105 NULLIFY (atomic_kind_set) 106 NULLIFY (qs_kind_set) 107 NULLIFY (fmstruct) 108 NULLIFY (fm_s_half) 109 NULLIFY (fm_work1) 110 NULLIFY (fm_work2) 111 NULLIFY (matrix_p) 112 NULLIFY (matrixkp_p) 113 NULLIFY (matrixkp_s) 114 NULLIFY (orbpop) 115 NULLIFY (particle_set) 116 NULLIFY (rho) 117 NULLIFY (scf_control) 118 NULLIFY (sm_p) 119 NULLIFY (sm_s) 120 NULLIFY (orbpop) 121 NULLIFY (para_env) 122 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 blacs_env=blacs_env) 133 134 CPASSERT(ASSOCIATED(atomic_kind_set)) 135 CPASSERT(ASSOCIATED(qs_kind_set)) 136 CPASSERT(ASSOCIATED(matrixkp_s)) 137 CPASSERT(ASSOCIATED(particle_set)) 138 CPASSERT(ASSOCIATED(rho)) 139 CPASSERT(ASSOCIATED(scf_control)) 140 141 IF (SIZE(matrixkp_s, 2) > 1) THEN 142 143 CPWARN("Lowdin population analysis not implemented for k-points.") 144 145 ELSE 146 147 sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format 148 CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format 149 150 matrix_p => matrixkp_p(:, 1) 151 nspin = SIZE(matrix_p, 1) 152 153 ! Get the total number of contracted spherical Gaussian basis functions 154 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf) 155 156 ! Provide an array to store the orbital populations for each spin 157 ALLOCATE (orbpop(nsgf, nspin)) 158 orbpop(:, :) = 0.0_dp 159 160 ! Write headline 161 IF (output_unit > 0) THEN 162 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 ncol_global=nsgf) 171 CALL cp_fm_create(matrix=fm_s_half, & 172 matrix_struct=fmstruct, & 173 name="S^(1/2) MATRIX") 174 CALL cp_fm_create(matrix=fm_work1, & 175 matrix_struct=fmstruct, & 176 name="FULL WORK MATRIX 1") 177 headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX" 178 CALL cp_fm_create(matrix=fm_work2, & 179 matrix_struct=fmstruct, & 180 name=TRIM(headline)) 181 CALL cp_fm_struct_release(fmstruct=fmstruct) 182 183 ! Build full S^(1/2) matrix (computationally expensive) 184 CALL copy_dbcsr_to_fm(sm_s, fm_s_half) 185 CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep) 186 IF (ndep /= 0) & 187 CALL cp_warn(__LOCATION__, & 188 "Overlap matrix exhibits linear dependencies. At least some "// & 189 "eigenvalues have been quenched.") 190 191 ! Build Lowdin population matrix for each spin 192 DO ispin = 1, nspin 193 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 CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf) 196 CALL cp_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 matrix_c=fm_work2) 206 IF (print_level > 2) THEN 207 ! Write the full Lowdin population matrix 208 IF (nspin > 1) THEN 209 IF (ispin == 1) THEN 210 fm_work2%name = TRIM(headline)//" FOR ALPHA SPIN" 211 ELSE 212 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 output_unit=output_unit) 217 END IF 218 CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin)) 219 END DO ! next spin ispin 220 221 ! Write atomic populations and charges 222 IF (output_unit > 0) THEN 223 print_gop = (print_level > 1) ! Print also orbital populations 224 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 CALL cp_fm_release(matrix=fm_s_half) 229 CALL cp_fm_release(matrix=fm_work1) 230 CALL cp_fm_release(matrix=fm_work2) 231 IF (ASSOCIATED(orbpop)) THEN 232 DEALLOCATE (orbpop) 233 END IF 234 235 END IF 236 237 CALL timestop(handle) 238 239 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 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 routineP = moduleN//':'//routineN 258 259 CHARACTER(LEN=default_string_length) :: headline 260 INTEGER :: blk, handle, iatom, ic, isgf, ispin, & 261 jatom, jsgf, natom, nsgf, nspin, sgfa, & 262 sgfb 263 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom 264 LOGICAL :: found, print_gop 265 REAL(KIND=dp) :: ps 266 REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop, p_block, ps_block, s_block 267 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 268 TYPE(cp_para_env_type), POINTER :: para_env 269 TYPE(dbcsr_iterator_type) :: iter 270 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 271 TYPE(dbcsr_type), POINTER :: sm_p, sm_ps, sm_s 272 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 273 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 274 TYPE(qs_rho_type), POINTER :: rho 275 276 CALL timeset(routineN, handle) 277 278 NULLIFY (atomic_kind_set) 279 NULLIFY (qs_kind_set) 280 NULLIFY (matrix_p) 281 NULLIFY (matrix_s) 282 NULLIFY (orbpop) 283 NULLIFY (particle_set) 284 NULLIFY (ps_block) 285 NULLIFY (p_block) 286 NULLIFY (rho) 287 NULLIFY (sm_p) 288 NULLIFY (sm_ps) 289 NULLIFY (sm_s) 290 NULLIFY (s_block) 291 NULLIFY (para_env) 292 293 CALL get_qs_env(qs_env=qs_env, & 294 atomic_kind_set=atomic_kind_set, & 295 qs_kind_set=qs_kind_set, & 296 matrix_s_kp=matrix_s, & 297 particle_set=particle_set, & 298 rho=rho, & 299 para_env=para_env) 300 301 CPASSERT(ASSOCIATED(atomic_kind_set)) 302 CPASSERT(ASSOCIATED(qs_kind_set)) 303 CPASSERT(ASSOCIATED(particle_set)) 304 CPASSERT(ASSOCIATED(rho)) 305 CPASSERT(ASSOCIATED(matrix_s)) 306 307 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format 308 nspin = SIZE(matrix_p, 1) 309 310 ! Get the total number of contracted spherical Gaussian basis functions 311 CALL get_atomic_kind_set(atomic_kind_set, natom=natom) 312 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf) 313 ALLOCATE (first_sgf_atom(natom)) 314 first_sgf_atom(:) = 0 315 316 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom) 317 318 ! Provide an array to store the orbital populations for each spin 319 ALLOCATE (orbpop(nsgf, nspin)) 320 orbpop(:, :) = 0.0_dp 321 322 ! Write headline 323 IF (output_unit > 0) THEN 324 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 325 '!-----------------------------------------------------------------------------!' 326 WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis" 327 END IF 328 329 ! Create a DBCSR work matrix, if needed 330 IF (print_level > 2) THEN 331 sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format 332 ALLOCATE (sm_ps) 333 headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX" 334 CALL dbcsr_copy(matrix_b=sm_ps, & 335 matrix_a=sm_s, & 336 name=TRIM(headline)) 337 END IF 338 339 ! Build Mulliken population matrix for each spin 340 DO ispin = 1, nspin 341 DO ic = 1, SIZE(matrix_s, 2) 342 IF (print_level > 2) THEN 343 CALL dbcsr_set(sm_ps, 0.0_dp) 344 END IF 345 sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format 346 sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format 347 ! Calculate Hadamard product of P and S as sparse matrix (Mulliken) 348 ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps) 349 CALL dbcsr_iterator_start(iter, sm_s) 350 DO WHILE (dbcsr_iterator_blocks_left(iter)) 351 CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block, blk) 352 IF (.NOT. (ASSOCIATED(s_block))) CYCLE 353 CALL dbcsr_get_block_p(matrix=sm_p, & 354 row=iatom, & 355 col=jatom, & 356 block=p_block, & 357 found=found) 358 IF (print_level > 2) THEN 359 CALL dbcsr_get_block_p(matrix=sm_ps, & 360 row=iatom, & 361 col=jatom, & 362 block=ps_block, & 363 found=found) 364 CPASSERT(ASSOCIATED(ps_block)) 365 END IF 366 367 sgfb = first_sgf_atom(jatom) 368 DO jsgf = 1, SIZE(s_block, 2) 369 DO isgf = 1, SIZE(s_block, 1) 370 ps = p_block(isgf, jsgf)*s_block(isgf, jsgf) 371 IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps 372 orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps 373 END DO 374 sgfb = sgfb + 1 375 END DO 376 IF (iatom /= jatom) THEN 377 sgfa = first_sgf_atom(iatom) 378 DO isgf = 1, SIZE(s_block, 1) 379 DO jsgf = 1, SIZE(s_block, 2) 380 ps = p_block(isgf, jsgf)*s_block(isgf, jsgf) 381 orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps 382 END DO 383 sgfa = sgfa + 1 384 END DO 385 END IF 386 END DO 387 CALL dbcsr_iterator_stop(iter) 388 END DO 389 390 IF (print_level > 2) THEN 391 ! Write the full Mulliken net AO and overlap population matrix 392 IF (nspin > 1) THEN 393 IF (ispin == 1) THEN 394 CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Alpha Spin") 395 ELSE 396 CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Beta Spin") 397 END IF 398 END IF 399 CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, & 400 output_unit=output_unit) 401 END IF 402 END DO 403 404 CALL mp_sum(orbpop, para_env%group) 405 406 ! Write atomic populations and charges 407 IF (output_unit > 0) THEN 408 print_gop = (print_level > 1) ! Print also orbital populations 409 CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop) 410 END IF 411 412 ! Release local working storage 413 IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps) 414 IF (ASSOCIATED(orbpop)) THEN 415 DEALLOCATE (orbpop) 416 END IF 417 IF (ALLOCATED(first_sgf_atom)) THEN 418 DEALLOCATE (first_sgf_atom) 419 END IF 420 421 IF (output_unit > 0) THEN 422 WRITE (UNIT=output_unit, FMT="(T2,A)") & 423 '!-----------------------------------------------------------------------------!' 424 END IF 425 426 CALL timestop(handle) 427 428 END SUBROUTINE mulliken_population_analysis 429 430! ************************************************************************************************** 431!> \brief Write atomic orbital populations and net atomic charges 432!> 433!> \param orbpop ... 434!> \param atomic_kind_set ... 435!> \param qs_kind_set ... 436!> \param particle_set ... 437!> \param output_unit ... 438!> \param print_orbital_contributions ... 439!> \date 07.07.2010 440!> \author Matthias Krack (MK) 441!> \version 1.0 442! ************************************************************************************************** 443 SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, & 444 print_orbital_contributions) 445 446 REAL(KIND=dp), DIMENSION(:, :), POINTER :: orbpop 447 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 448 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 449 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 450 INTEGER, INTENT(IN) :: output_unit 451 LOGICAL, INTENT(IN) :: print_orbital_contributions 452 453 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_orbpop', routineP = moduleN//':'//routineN 454 455 CHARACTER(LEN=2) :: element_symbol 456 CHARACTER(LEN=6), DIMENSION(:), POINTER :: sgf_symbol 457 INTEGER :: handle, iao, iatom, ikind, iset, isgf, & 458 ishell, iso, l, natom, nset, nsgf, & 459 nspin 460 INTEGER, DIMENSION(:), POINTER :: nshell 461 INTEGER, DIMENSION(:, :), POINTER :: lshell 462 REAL(KIND=dp) :: zeff 463 REAL(KIND=dp), DIMENSION(3) :: sumorbpop, totsumorbpop 464 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 465 466 CALL timeset(routineN, handle) 467 468 NULLIFY (lshell) 469 NULLIFY (nshell) 470 NULLIFY (orb_basis_set) 471 NULLIFY (sgf_symbol) 472 473 CPASSERT(ASSOCIATED(orbpop)) 474 CPASSERT(ASSOCIATED(atomic_kind_set)) 475 CPASSERT(ASSOCIATED(particle_set)) 476 477 nspin = SIZE(orbpop, 2) 478 479 CALL get_atomic_kind_set(atomic_kind_set, natom=natom) 480 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf) 481 482 ! Select and write headline 483 IF (nspin == 1) THEN 484 IF (print_orbital_contributions) THEN 485 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 486 "# Orbital AO symbol Orbital population Net charge" 487 ELSE 488 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 489 "# Atom Element Kind Atomic population Net charge" 490 END IF 491 ELSE 492 IF (print_orbital_contributions) THEN 493 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 494 "# Orbital AO symbol Orbital population (alpha,beta) Net charge Spin moment" 495 ELSE 496 WRITE (UNIT=output_unit, FMT="(/,T2,A)") & 497 "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment" 498 END IF 499 END IF 500 501 totsumorbpop(:) = 0.0_dp 502 503 iao = 1 504 DO iatom = 1, natom 505 sumorbpop(:) = 0.0_dp 506 NULLIFY (orb_basis_set) 507 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 508 element_symbol=element_symbol, & 509 kind_number=ikind) 510 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff) 511 IF (ASSOCIATED(orb_basis_set)) THEN 512 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 513 nset=nset, & 514 nshell=nshell, & 515 l=lshell, & 516 sgf_symbol=sgf_symbol) 517 isgf = 1 518 DO iset = 1, nset 519 DO ishell = 1, nshell(iset) 520 l = lshell(ishell, iset) 521 DO iso = 1, nso(l) 522 IF (nspin == 1) THEN 523 sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1) 524 IF (print_orbital_contributions) THEN 525 IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") "" 526 WRITE (UNIT=output_unit, & 527 FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") & 528 iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1) 529 END IF 530 ELSE 531 sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2) 532 sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2) 533 IF (print_orbital_contributions) THEN 534 IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") "" 535 WRITE (UNIT=output_unit, & 536 FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") & 537 iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), & 538 orbpop(iao, 1) - orbpop(iao, 2) 539 END IF 540 END IF 541 isgf = isgf + 1 542 iao = iao + 1 543 END DO 544 END DO 545 END DO 546 IF (nspin == 1) THEN 547 totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1) 548 totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) 549 WRITE (UNIT=output_unit, & 550 FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") & 551 iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1) 552 ELSE 553 totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2) 554 totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2) 555 WRITE (UNIT=output_unit, & 556 FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") & 557 iatom, element_symbol, ikind, sumorbpop(1:2), & 558 zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3) 559 END IF 560 END IF ! atom has an orbital basis 561 END DO ! next atom iatom 562 563 ! Write total sums 564 IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") "" 565 IF (nspin == 1) THEN 566 WRITE (UNIT=output_unit, & 567 FMT="(T2,A,T42,F12.6,T68,F12.6,/)") & 568 "# Total charge", totsumorbpop(1), totsumorbpop(3) 569 ELSE 570 WRITE (UNIT=output_unit, & 571 FMT="(T2,A,T28,4(1X,F12.6),/)") & 572 "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), & 573 totsumorbpop(1) - totsumorbpop(2) 574 END IF 575 576 IF (output_unit > 0) CALL m_flush(output_unit) 577 578 CALL timestop(handle) 579 580 END SUBROUTINE write_orbpop 581 582END MODULE population_analyses 583