1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Does all kind of post scf calculations for semi-empirical 8!> \par History 9!> Started printing preliminary stuff for MO_CUBES and MO requires some 10!> more work to complete all other functionalities 11!> \author Teodoro Laino (07.2008) 12! ************************************************************************************************** 13MODULE qs_scf_post_se 14! 15 USE ai_moments, ONLY: moment 16 USE atomic_kind_types, ONLY: atomic_kind_type,& 17 get_atomic_kind 18 USE basis_set_types, ONLY: gto_basis_set_type 19 USE cell_types, ONLY: cell_type,& 20 pbc 21 USE cp_control_types, ONLY: dft_control_type 22 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 23 USE cp_log_handling, ONLY: cp_get_default_logger,& 24 cp_logger_get_default_io_unit,& 25 cp_logger_type 26 USE cp_output_handling, ONLY: cp_p_file,& 27 cp_print_key_finished_output,& 28 cp_print_key_should_output,& 29 cp_print_key_unit_nr 30 USE cp_para_types, ONLY: cp_para_env_type 31 USE cp_result_methods, ONLY: cp_results_erase,& 32 put_results 33 USE cp_result_types, ONLY: cp_result_type 34 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 35 dbcsr_p_type 36 USE input_section_types, ONLY: section_get_ival,& 37 section_vals_get,& 38 section_vals_get_subs_vals,& 39 section_vals_type,& 40 section_vals_val_get 41 USE kinds, ONLY: default_string_length,& 42 dp 43 USE mathconstants, ONLY: twopi 44 USE message_passing, ONLY: mp_sum 45 USE moments_utils, ONLY: get_reference_point 46 USE orbital_pointers, ONLY: coset,& 47 ncoset 48 USE particle_list_types, ONLY: particle_list_type 49 USE particle_types, ONLY: particle_type 50 USE physcon, ONLY: debye 51 USE qs_environment_types, ONLY: get_qs_env,& 52 qs_environment_type 53 USE qs_kind_types, ONLY: get_qs_kind,& 54 qs_kind_type 55 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 56 USE qs_ks_types, ONLY: qs_ks_did_change 57 USE qs_mo_io, ONLY: write_mo_set 58 USE qs_mo_types, ONLY: mo_set_p_type 59 USE qs_rho_types, ONLY: qs_rho_get,& 60 qs_rho_type 61 USE qs_subsys_types, ONLY: qs_subsys_get,& 62 qs_subsys_type 63 USE semi_empirical_types, ONLY: get_se_param,& 64 semi_empirical_type 65#include "./base/base_uses.f90" 66 67 IMPLICIT NONE 68 PRIVATE 69 70 ! Global parameters 71 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_se' 72 PUBLIC :: scf_post_calculation_se 73 74CONTAINS 75 76! ************************************************************************************************** 77!> \brief collects possible post - scf calculations and prints info / computes properties. 78!> specific for Semi-empirical calculations 79!> \param qs_env the qs_env in which the qs_env lives 80!> \par History 81!> 07.2008 created [tlaino] - Split from qs_scf_post (general) 82!> \author tlaino 83!> \note 84!> this function changes mo_eigenvectors and mo_eigenvalues, depending on the print keys. 85!> In particular, MO_CUBES causes the MOs to be rotated to make them eigenstates of the KS 86!> matrix, and mo_eigenvalues is updated accordingly. This can, for unconverged wavefunctions, 87!> change afterwards slightly the forces (hence small numerical differences between MD 88!> with and without the debug print level). Ideally this should not happen... 89! ************************************************************************************************** 90 SUBROUTINE scf_post_calculation_se(qs_env) 91 92 TYPE(qs_environment_type), POINTER :: qs_env 93 94 CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_se', & 95 routineP = moduleN//':'//routineN 96 97 INTEGER :: handle, output_unit 98 LOGICAL :: explicit, my_localized_wfn 99 TYPE(cp_logger_type), POINTER :: logger 100 TYPE(cp_para_env_type), POINTER :: para_env 101 TYPE(dft_control_type), POINTER :: dft_control 102 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 103 TYPE(particle_list_type), POINTER :: particles 104 TYPE(qs_rho_type), POINTER :: rho 105 TYPE(qs_subsys_type), POINTER :: subsys 106 TYPE(section_vals_type), POINTER :: input, print_key, wfn_mix_section 107 108 CALL timeset(routineN, handle) 109 110 ! Writes the data that is already available in qs_env 111 CALL write_available_results(qs_env) 112 113 my_localized_wfn = .FALSE. 114 NULLIFY (dft_control, mos, rho, & 115 subsys, particles, input, print_key, para_env) 116 117 logger => cp_get_default_logger() 118 output_unit = cp_logger_get_default_io_unit(logger) 119 120 CPASSERT(ASSOCIATED(qs_env)) 121 ! Here we start with data that needs a postprocessing... 122 CALL get_qs_env(qs_env, & 123 dft_control=dft_control, & 124 mos=mos, & 125 rho=rho, & 126 input=input, & 127 subsys=subsys, & 128 para_env=para_env) 129 CALL qs_subsys_get(subsys, particles=particles) 130 131 ! Compute Atomic Charges 132 CALL qs_scf_post_charges(input, logger, qs_env, rho, para_env) 133 134 ! Moments of charge distribution 135 CALL qs_scf_post_moments(input, logger, qs_env) 136 137 ! MO_CUBES 138 print_key => section_vals_get_subs_vals(section_vals=input, & 139 subsection_name="DFT%PRINT%MO_CUBES") 140 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 141 CPWARN("Printing of MO cube files not implemented for Semi-Empirical method.") 142 END IF 143 144 ! STM 145 print_key => section_vals_get_subs_vals(section_vals=input, & 146 subsection_name="DFT%PRINT%STM") 147 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 148 CPWARN("STM not implemented for Semi-Empirical method.") 149 END IF 150 151 ! DFT+U 152 print_key => section_vals_get_subs_vals(section_vals=input, & 153 subsection_name="DFT%PRINT%PLUS_U") 154 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 155 CPWARN("DFT+U not available for Semi-Empirical method.") 156 END IF 157 158 ! Kinetic Energy 159 print_key => section_vals_get_subs_vals(section_vals=input, & 160 subsection_name="DFT%PRINT%KINETIC_ENERGY") 161 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 162 CPWARN("Kinetic energy not available for Semi-Empirical method.") 163 END IF 164 165 ! Wavefunction mixing 166 wfn_mix_section => section_vals_get_subs_vals(input, "DFT%PRINT%WFN_MIX") 167 CALL section_vals_get(wfn_mix_section, explicit=explicit) 168 IF (explicit .AND. .NOT. qs_env%run_rtp) THEN 169 CPWARN("Wavefunction mixing not implemented for Semi-Empirical method.") 170 END IF 171 172 ! Print coherent X-ray diffraction spectrum 173 print_key => section_vals_get_subs_vals(section_vals=input, & 174 subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM") 175 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 176 CPWARN("XRAY_DIFFRACTION_SPECTRUM not implemented for Semi-Empirical calculations!!") 177 END IF 178 179 ! Calculation of Electric Field Gradients 180 print_key => section_vals_get_subs_vals(section_vals=input, & 181 subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT") 182 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 183 CPWARN("ELECTRIC_FIELD_GRADIENT not implemented for Semi-Empirical calculations!!") 184 END IF 185 186 ! Calculation of EPR Hyperfine Coupling Tensors 187 print_key => section_vals_get_subs_vals(section_vals=input, & 188 subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR") 189 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), & 190 cp_p_file)) THEN 191 CPWARN("HYPERFINE_COUPLING_TENSOR not implemented for Semi-Empirical calculations!!") 192 END IF 193 194 CALL timestop(handle) 195 196 END SUBROUTINE scf_post_calculation_se 197 198! ************************************************************************************************** 199!> \brief Computes and prints electric dipole moments 200!> We use the approximation for NDDO from 201!> Pople and Beveridge, Approximate Molecular Orbital Theory, 202!> Mc Graw Hill 1970 203!> mu = \sum_A [ Q_A * R_a + Tr(P_A*D_A) ] 204!> \param input ... 205!> \param logger ... 206!> \param qs_env the qs_env in which the qs_env lives 207! ************************************************************************************************** 208 SUBROUTINE qs_scf_post_moments(input, logger, qs_env) 209 TYPE(section_vals_type), POINTER :: input 210 TYPE(cp_logger_type), POINTER :: logger 211 TYPE(qs_environment_type), POINTER :: qs_env 212 213 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_moments', & 214 routineP = moduleN//':'//routineN 215 216 CHARACTER(LEN=default_string_length) :: description, dipole_type 217 COMPLEX(KIND=dp) :: dzeta, zeta 218 COMPLEX(KIND=dp), DIMENSION(3) :: dggamma, dzphase, ggamma, zphase 219 INTEGER :: i, iat, iatom, ikind, ix, j, nat, natom, & 220 natorb, nkind, nspin, reference, & 221 unit_nr 222 LOGICAL :: do_berry, found 223 REAL(KIND=dp) :: charge_tot, ci(3), dci(3), dipole(3), dipole_deriv(3), drcc(3), dria(3), & 224 dtheta, gvec(3), q, rcc(3), ria(3), tcharge(2), theta, tmp(3), via(3), zeff 225 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ncharge 226 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mom 227 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 228 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock 229 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 230 TYPE(cell_type), POINTER :: cell 231 TYPE(cp_result_type), POINTER :: results 232 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p 233 TYPE(dft_control_type), POINTER :: dft_control 234 TYPE(gto_basis_set_type), POINTER :: basis_set 235 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 236 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 237 TYPE(qs_rho_type), POINTER :: rho 238 TYPE(section_vals_type), POINTER :: print_key 239 TYPE(semi_empirical_type), POINTER :: se_kind 240 241 NULLIFY (results) 242 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MOMENTS") 243 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 244 ! Dipole Moments 245 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MOMENTS", & 246 extension=".data", middle_name="se_dipole", log_filename=.FALSE.) 247 248 ! Reference point 249 reference = section_get_ival(print_key, keyword_name="REFERENCE") 250 NULLIFY (ref_point) 251 description = '[DIPOLE]' 252 CALL section_vals_val_get(print_key, "REF_POINT", r_vals=ref_point) 253 CALL section_vals_val_get(print_key, "PERIODIC", l_val=do_berry) 254 CALL get_reference_point(rcc, drcc, qs_env=qs_env, reference=reference, & 255 ref_point=ref_point) 256 ! 257 NULLIFY (particle_set) 258 CALL get_qs_env(qs_env=qs_env, & 259 rho=rho, & 260 cell=cell, & 261 atomic_kind_set=atomic_kind_set, & 262 natom=natom, & 263 qs_kind_set=qs_kind_set, & 264 particle_set=particle_set, & 265 results=results, & 266 dft_control=dft_control) 267 268 CALL qs_rho_get(rho, rho_ao=matrix_p) 269 nspin = SIZE(matrix_p) 270 nkind = SIZE(atomic_kind_set) 271 ! net charges 272 ALLOCATE (ncharge(natom)) 273 ncharge = 0.0_dp 274 DO ikind = 1, nkind 275 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 276 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind) 277 CALL get_se_param(se_kind, zeff=zeff, natorb=natorb) 278 DO iatom = 1, nat 279 iat = atomic_kind_set(ikind)%atom_list(iatom) 280 tcharge = 0.0_dp 281 DO i = 1, nspin 282 CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, & 283 block=pblock, found=found) 284 IF (found) THEN 285 DO j = 1, natorb 286 tcharge(i) = tcharge(i) + pblock(j, j) 287 END DO 288 END IF 289 END DO 290 ncharge(iat) = zeff - SUM(tcharge) 291 END DO 292 END DO 293 ! Contributions from net atomic charges 294 ! Dipole deriv will be the derivative of the Dipole(dM/dt=\sum e_j v_j) 295 dipole_deriv = 0.0_dp 296 dipole = 0.0_dp 297 IF (do_berry) THEN 298 dipole_type = "[BERRY PHASE]" 299 rcc = pbc(rcc, cell) 300 charge_tot = 0._dp 301 charge_tot = SUM(ncharge) 302 ria = twopi*MATMUL(cell%h_inv, rcc) 303 zphase = CMPLX(COS(ria), SIN(ria), dp)**charge_tot 304 305 dria = twopi*MATMUL(cell%h_inv, drcc) 306 dzphase = charge_tot*CMPLX(-SIN(ria), COS(ria), dp)**(charge_tot - 1.0_dp)*dria 307 308 ggamma = CMPLX(1.0_dp, 0.0_dp, KIND=dp) 309 dggamma = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 310 DO ikind = 1, SIZE(atomic_kind_set) 311 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 312 DO i = 1, nat 313 iat = atomic_kind_set(ikind)%atom_list(i) 314 ria = particle_set(iat)%r(:) 315 ria = pbc(ria, cell) 316 via = particle_set(iat)%v(:) 317 q = ncharge(iat) 318 DO j = 1, 3 319 gvec = twopi*cell%h_inv(j, :) 320 theta = SUM(ria(:)*gvec(:)) 321 dtheta = SUM(via(:)*gvec(:)) 322 zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-q) 323 dzeta = -q*CMPLX(-SIN(theta), COS(theta), KIND=dp)**(-q - 1.0_dp)*dtheta 324 dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta 325 ggamma(j) = ggamma(j)*zeta 326 END DO 327 ENDDO 328 END DO 329 dggamma = dggamma*zphase + ggamma*dzphase 330 ggamma = ggamma*zphase 331 IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN 332 tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp) 333 ci = ATAN(tmp) 334 dci = (1.0_dp/(1.0_dp + tmp**2))* & 335 (AIMAG(dggamma)*REAL(ggamma, KIND=dp) - AIMAG(ggamma)* & 336 REAL(dggamma, KIND=dp))/(REAL(ggamma, KIND=dp))**2 337 dipole = MATMUL(cell%hmat, ci)/twopi 338 dipole_deriv = MATMUL(cell%hmat, dci)/twopi 339 END IF 340 ELSE 341 dipole_type = "[Non Periodic]" 342 DO i = 1, natom 343 ! no pbc(particle_set(i)%r(:),cell) so that the total dipole is the sum of the molecular dipoles 344 ria = particle_set(i)%r(:) 345 q = ncharge(i) 346 dipole = dipole - q*(ria - rcc) 347 dipole_deriv(:) = dipole_deriv(:) - q*(particle_set(i)%v(:) - drcc) 348 END DO 349 END IF 350 ! Contributions from atomic polarization 351 ! No contribution to dipole derivatives 352 DO ikind = 1, nkind 353 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 354 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set) 355 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind) 356 CALL get_se_param(se_kind, natorb=natorb) 357 ALLOCATE (mom(natorb, natorb, 3)) 358 mom = 0.0_dp 359 CALL atomic_moments(mom, basis_set) 360 DO iatom = 1, nat 361 iat = atomic_kind_set(ikind)%atom_list(iatom) 362 DO i = 1, nspin 363 CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, & 364 block=pblock, found=found) 365 IF (found) THEN 366 CPASSERT(natorb == SIZE(pblock, 1)) 367 ix = coset(1, 0, 0) - 1 368 dipole(1) = dipole(1) + SUM(pblock*mom(:, :, ix)) 369 ix = coset(0, 1, 0) - 1 370 dipole(2) = dipole(2) + SUM(pblock*mom(:, :, ix)) 371 ix = coset(0, 0, 1) - 1 372 dipole(3) = dipole(3) + SUM(pblock*mom(:, :, ix)) 373 END IF 374 END DO 375 END DO 376 DEALLOCATE (mom) 377 END DO 378 CALL cp_results_erase(results=results, description=description) 379 CALL put_results(results=results, description=description, & 380 values=dipole(1:3)) 381 IF (unit_nr > 0) THEN 382 WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//"(A.U.)|", dipole 383 WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//"(Debye)|", dipole*debye 384 WRITE (unit_nr, '(1X,A,T48,3F11.6)') "DIPOLE "//TRIM(dipole_type)//" DERIVATIVE(A.U.)|", dipole_deriv 385 END IF 386 CALL cp_print_key_finished_output(unit_nr, logger, print_key) 387 END IF 388 389 END SUBROUTINE qs_scf_post_moments 390 391! ************************************************************************************************** 392!> \brief Computes the dipole integrals for an atom (a|x|b), a,b on atom A 393!> \param mom ... 394!> \param basis_set ... 395! ************************************************************************************************** 396 SUBROUTINE atomic_moments(mom, basis_set) 397 REAL(KIND=dp), DIMENSION(:, :, :) :: mom 398 TYPE(gto_basis_set_type), POINTER :: basis_set 399 400 CHARACTER(len=*), PARAMETER :: routineN = 'atomic_moments', routineP = moduleN//':'//routineN 401 402 INTEGER :: i, iset, jset, ncoa, ncob, nm, nset, & 403 sgfa, sgfb 404 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgf, nsgf 405 INTEGER, DIMENSION(:, :), POINTER :: first_sgf 406 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work 407 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab 408 REAL(KIND=dp), DIMENSION(3) :: rac, rbc 409 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf, sphi, zet 410 411 rac = 0.0_dp 412 rbc = 0.0_dp 413 414 first_sgf => basis_set%first_sgf 415 la_max => basis_set%lmax 416 la_min => basis_set%lmin 417 npgf => basis_set%npgf 418 nset = basis_set%nset 419 nsgf => basis_set%nsgf_set 420 rpgf => basis_set%pgf_radius 421 sphi => basis_set%sphi 422 zet => basis_set%zet 423 424 nm = 0 425 DO iset = 1, nset 426 ncoa = npgf(iset)*ncoset(la_max(iset)) 427 nm = MAX(nm, ncoa) 428 END DO 429 ALLOCATE (mab(nm, nm, 4), work(nm, nm)) 430 431 DO iset = 1, nset 432 ncoa = npgf(iset)*ncoset(la_max(iset)) 433 sgfa = first_sgf(1, iset) 434 DO jset = 1, nset 435 ncob = npgf(jset)*ncoset(la_max(jset)) 436 sgfb = first_sgf(1, jset) 437 !*** Calculate the primitive integrals *** 438 CALL moment(la_max(iset), npgf(iset), zet(:, iset), rpgf(:, iset), la_min(iset), & 439 la_max(jset), npgf(jset), zet(:, jset), rpgf(:, jset), 1, rac, rbc, mab) 440 !*** Contraction step *** 441 DO i = 1, 3 442 CALL dgemm("N", "N", ncoa, nsgf(jset), ncob, 1.0_dp, mab(1, 1, i), SIZE(mab, 1), & 443 sphi(1, sgfb), SIZE(sphi, 1), 0.0_dp, work(1, 1), SIZE(work, 1)) 444 CALL dgemm("T", "N", nsgf(iset), nsgf(jset), ncoa, 1.0_dp, sphi(1, sgfa), SIZE(sphi, 1), & 445 work(1, 1), SIZE(work, 1), 1.0_dp, mom(sgfa, sgfb, i), SIZE(mom, 1)) 446 END DO 447 END DO 448 END DO 449 DEALLOCATE (mab, work) 450 451 END SUBROUTINE atomic_moments 452! ************************************************************************************************** 453!> \brief Computes and Prints Atomic Charges with several methods 454!> \param input ... 455!> \param logger ... 456!> \param qs_env the qs_env in which the qs_env lives 457!> \param rho ... 458!> \param para_env ... 459! ************************************************************************************************** 460 SUBROUTINE qs_scf_post_charges(input, logger, qs_env, rho, & 461 para_env) 462 TYPE(section_vals_type), POINTER :: input 463 TYPE(cp_logger_type), POINTER :: logger 464 TYPE(qs_environment_type), POINTER :: qs_env 465 TYPE(qs_rho_type), POINTER :: rho 466 TYPE(cp_para_env_type), POINTER :: para_env 467 468 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_charges', & 469 routineP = moduleN//':'//routineN 470 471 CHARACTER(LEN=2) :: ana 472 CHARACTER(LEN=default_string_length) :: aname 473 INTEGER :: i, iat, iatom, ikind, j, nat, natom, & 474 natorb, nkind, nspin, unit_nr 475 LOGICAL :: found 476 REAL(KIND=dp) :: zeff 477 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge 478 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: charges 479 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock 480 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 481 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p 482 TYPE(dft_control_type), POINTER :: dft_control 483 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 484 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 485 TYPE(section_vals_type), POINTER :: print_key 486 TYPE(semi_empirical_type), POINTER :: se_kind 487 488 NULLIFY (particle_set) 489 CALL get_qs_env(qs_env=qs_env, & 490 atomic_kind_set=atomic_kind_set, & 491 natom=natom, & 492 qs_kind_set=qs_kind_set, & 493 particle_set=particle_set, & 494 dft_control=dft_control) 495 496 ! Compute the mulliken charges 497 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MULLIKEN") 498 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 499 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MULLIKEN", extension=".mulliken", log_filename=.FALSE.) 500 CALL qs_rho_get(rho, rho_ao=matrix_p) 501 nspin = SIZE(matrix_p) 502 ALLOCATE (charges(natom, nspin), mcharge(natom)) 503 charges = 0.0_dp 504 mcharge = 0.0_dp 505 ! calculate atomic charges 506 nkind = SIZE(atomic_kind_set) 507 DO ikind = 1, nkind 508 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 509 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind) 510 CALL get_se_param(se_kind, zeff=zeff, natorb=natorb) 511 DO iatom = 1, nat 512 iat = atomic_kind_set(ikind)%atom_list(iatom) 513 DO i = 1, nspin 514 CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, & 515 block=pblock, found=found) 516 IF (found) THEN 517 DO j = 1, natorb 518 charges(iat, i) = charges(iat, i) + pblock(j, j) 519 END DO 520 END IF 521 END DO 522 mcharge(iat) = zeff - SUM(charges(iat, 1:nspin)) 523 END DO 524 END DO 525 ! 526 CALL mp_sum(charges, para_env%group) 527 CALL mp_sum(mcharge, para_env%group) 528 ! 529 IF (unit_nr > 0) THEN 530 WRITE (UNIT=unit_nr, FMT="(/,/,T2,A)") "POPULATION ANALYSIS" 531 IF (nspin == 1) THEN 532 WRITE (UNIT=unit_nr, FMT="(/,T2,A,T70,A)") & 533 " # Atom Element Kind Atomic population", " Net charge" 534 DO ikind = 1, nkind 535 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 536 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind) 537 CALL get_se_param(se_kind, name=aname) 538 ana = ADJUSTR(TRIM(ADJUSTL(aname))) 539 DO iatom = 1, nat 540 iat = atomic_kind_set(ikind)%atom_list(iatom) 541 WRITE (UNIT=unit_nr, & 542 FMT="(T2,I7,6X,A2,3X,I6,T39,F12.6,T69,F12.6)") & 543 iat, ana, ikind, charges(iat, 1), mcharge(iat) 544 END DO 545 END DO 546 WRITE (UNIT=unit_nr, & 547 FMT="(T2,A,T39,F12.6,T69,F12.6,/)") & 548 "# Total charge", SUM(charges(:, 1)), SUM(mcharge(:)) 549 ELSE 550 WRITE (UNIT=unit_nr, FMT="(/,T2,A)") & 551 "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment" 552 DO ikind = 1, nkind 553 CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat) 554 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind) 555 CALL get_se_param(se_kind, name=aname) 556 ana = ADJUSTR(TRIM(ADJUSTL(aname))) 557 DO iatom = 1, nat 558 iat = atomic_kind_set(ikind)%atom_list(iatom) 559 WRITE (UNIT=unit_nr, & 560 FMT="(T2,I6,5X,A2,2X,I6,T29,4(1X,F12.6))") & 561 iat, ana, ikind, charges(iat, 1:2), mcharge(iat), charges(iat, 1) - charges(iat, 2) 562 END DO 563 END DO 564 WRITE (UNIT=unit_nr, & 565 FMT="(T2,A,T29,4(1X,F12.6),/)") & 566 "# Total charge and spin", SUM(charges(:, 1)), SUM(charges(:, 2)), SUM(mcharge(:)) 567 END IF 568 END IF 569 570 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN") 571 572 DEALLOCATE (charges, mcharge) 573 END IF 574 575 ! Compute the Lowdin charges 576 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%LOWDIN") 577 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 578 CPWARN("Lowdin charges not available for semi-empirical calculations!") 579 END IF 580 581 ! Hirshfeld charges 582 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%HIRSHFELD") 583 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 584 CPWARN("Hirshfeld charges not available for semi-empirical calculations!") 585 END IF 586 587 ! MAO 588 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MAO_ANALYSIS") 589 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 590 CPWARN("MAO analysis not available for semi-empirical calculations!") 591 END IF 592 593 END SUBROUTINE qs_scf_post_charges 594 595! ************************************************************************************************** 596!> \brief Write QS results always available (if switched on through the print_keys) 597!> \param qs_env the qs_env in which the qs_env lives 598! ************************************************************************************************** 599 SUBROUTINE write_available_results(qs_env) 600 TYPE(qs_environment_type), POINTER :: qs_env 601 602 CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results', & 603 routineP = moduleN//':'//routineN 604 605 INTEGER :: after, handle, ispin, iw, output_unit 606 LOGICAL :: omit_headers 607 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 608 TYPE(cell_type), POINTER :: cell 609 TYPE(cp_logger_type), POINTER :: logger 610 TYPE(cp_para_env_type), POINTER :: para_env 611 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, rho_ao 612 TYPE(dft_control_type), POINTER :: dft_control 613 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 614 TYPE(particle_list_type), POINTER :: particles 615 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 616 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 617 TYPE(qs_rho_type), POINTER :: rho 618 TYPE(qs_subsys_type), POINTER :: subsys 619 TYPE(section_vals_type), POINTER :: dft_section, input 620 621 CALL timeset(routineN, handle) 622 NULLIFY (cell, dft_control, mos, atomic_kind_set, particle_set, rho, & 623 ks_rmpv, dft_section, input, & 624 particles, subsys, para_env, rho_ao) 625 logger => cp_get_default_logger() 626 output_unit = cp_logger_get_default_io_unit(logger) 627 628 CPASSERT(ASSOCIATED(qs_env)) 629 CALL get_qs_env(qs_env, & 630 dft_control=dft_control, & 631 mos=mos, & 632 atomic_kind_set=atomic_kind_set, & 633 qs_kind_set=qs_kind_set, & 634 particle_set=particle_set, & 635 rho=rho, & 636 matrix_ks=ks_rmpv, & 637 input=input, & 638 cell=cell, & 639 subsys=subsys, & 640 para_env=para_env) 641 CALL qs_subsys_get(subsys, particles=particles) 642 643 CALL qs_rho_get(rho, rho_ao=rho_ao) 644 ! *** if the dft_section tells you to do so, write last wavefunction to screen 645 dft_section => section_vals_get_subs_vals(input, "DFT") 646 IF (dft_control%nspins == 2) THEN 647 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 648 dft_section, spin="ALPHA", last=.TRUE.) 649 CALL write_mo_set(mos(2)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 650 dft_section, spin="BETA", last=.TRUE.) 651 ELSE 652 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 653 dft_section, last=.TRUE.) 654 END IF 655 656 ! *** at the end of scf print out the projected dos per kind 657 IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%PDOS") & 658 , cp_p_file)) THEN 659 CPWARN("PDOS not implemented for Semi-Empirical calculations!!") 660 ENDIF 661 662 ! Print the total density (electronic + core charge) 663 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 664 "DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN 665 CPWARN("TOT_DENSITY_CUBE not implemented for Semi-Empirical calculations!!") 666 END IF 667 668 ! Write cube file with electron density 669 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 670 "DFT%PRINT%E_DENSITY_CUBE"), cp_p_file)) THEN 671 CPWARN("E_DENSITY_CUBE not implemented for Semi-Empirical calculations!!") 672 END IF ! print key 673 674 ! Write cube file with EFIELD 675 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 676 "DFT%PRINT%EFIELD_CUBE"), cp_p_file)) THEN 677 CPWARN("EFIELD_CUBE not implemented for Semi-Empirical calculations!!") 678 END IF ! print key 679 680 ! Write cube file with ELF 681 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 682 "DFT%PRINT%ELF_CUBE"), cp_p_file)) THEN 683 CPWARN("ELF function not implemented for Semi-Empirical calculations!!") 684 END IF ! print key 685 686 ! Print the hartree potential 687 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 688 "DFT%PRINT%V_HARTREE_CUBE"), cp_p_file)) THEN 689 CPWARN("V_HARTREE_CUBE not implemented for Semi-Empirical calculations!!") 690 ENDIF 691 692 ! Print the XC potential 693 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 694 "DFT%PRINT%V_XC_CUBE"), cp_p_file)) THEN 695 CPWARN("V_XC_CUBE not available for Semi-Empirical calculations!!") 696 ENDIF 697 698 ! Write the density matrix 699 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 700 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 701 "DFT%PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN 702 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/DENSITY", & 703 extension=".Log") 704 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 705 after = MIN(MAX(after, 1), 16) 706 DO ispin = 1, dft_control%nspins 707 CALL cp_dbcsr_write_sparse_matrix(rho_ao(ispin)%matrix, 4, after, qs_env, & 708 para_env, output_unit=iw, omit_headers=omit_headers) 709 END DO 710 CALL cp_print_key_finished_output(iw, logger, input, & 711 "DFT%PRINT%AO_MATRICES/DENSITY") 712 END IF 713 714 ! The Kohn-Sham matrix itself 715 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 716 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN 717 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.) 718 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 719 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", & 720 extension=".Log") 721 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 722 after = MIN(MAX(after, 1), 16) 723 CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(1)%matrix, 4, after, qs_env, & 724 para_env, output_unit=iw, omit_headers=omit_headers) 725 CALL cp_print_key_finished_output(iw, logger, input, & 726 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX") 727 END IF 728 729 CALL timestop(handle) 730 731 END SUBROUTINE write_available_results 732 733END MODULE qs_scf_post_se 734