1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines to somehow generate an initial guess 8!> \par History 9!> 2006.03 Moved here from qs_scf.F [Joost VandeVondele] 10! ************************************************************************************************** 11MODULE qs_initial_guess 12 USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind,& 15 get_atomic_kind_set 16 USE basis_set_types, ONLY: get_gto_basis_set,& 17 gto_basis_set_type 18 USE cp_control_types, ONLY: dft_control_type 19 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 20 copy_fm_to_dbcsr,& 21 cp_dbcsr_sm_fm_multiply,& 22 cp_fm_to_dbcsr_row_template 23 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 24 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 25 cp_fm_struct_get,& 26 cp_fm_struct_release,& 27 cp_fm_struct_type 28 USE cp_fm_types, ONLY: & 29 cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, cp_fm_p_type, & 30 cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type 31 USE cp_log_handling, ONLY: cp_get_default_logger,& 32 cp_logger_get_default_io_unit,& 33 cp_logger_type,& 34 cp_to_string 35 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 36 cp_print_key_unit_nr 37 USE cp_para_types, ONLY: cp_para_env_type 38 USE dbcsr_api, ONLY: & 39 dbcsr_add_on_diag, dbcsr_checksum, dbcsr_copy, dbcsr_dot, dbcsr_filter, dbcsr_get_diag, & 40 dbcsr_get_info, dbcsr_get_num_blocks, dbcsr_get_occupation, dbcsr_iterator_blocks_left, & 41 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 42 dbcsr_multiply, dbcsr_nfullrows_total, dbcsr_p_type, dbcsr_release, dbcsr_scale, & 43 dbcsr_set, dbcsr_set_diag, dbcsr_type, dbcsr_verify_matrix 44 USE external_potential_types, ONLY: all_potential_type,& 45 gth_potential_type,& 46 sgp_potential_type 47 USE input_constants, ONLY: atomic_guess,& 48 core_guess,& 49 history_guess,& 50 mopac_guess,& 51 no_guess,& 52 random_guess,& 53 restart_guess,& 54 sparse_guess 55 USE input_section_types, ONLY: section_vals_get_subs_vals,& 56 section_vals_type 57 USE kinds, ONLY: default_path_length,& 58 dp 59 USE kpoint_io, ONLY: read_kpoints_restart 60 USE kpoint_types, ONLY: kpoint_type 61 USE message_passing, ONLY: mp_bcast,& 62 mp_sum 63 USE particle_methods, ONLY: get_particle_set 64 USE particle_types, ONLY: particle_type 65 USE qs_dftb_utils, ONLY: get_dftb_atom_param 66 USE qs_environment_types, ONLY: get_qs_env,& 67 qs_environment_type 68 USE qs_kind_types, ONLY: get_qs_kind,& 69 get_qs_kind_set,& 70 qs_kind_type 71 USE qs_mo_io, ONLY: read_mo_set,& 72 wfn_restart_file_name 73 USE qs_mo_methods, ONLY: calculate_density_matrix,& 74 make_basis_lowdin,& 75 make_basis_simple,& 76 make_basis_sm 77 USE qs_mo_occupation, ONLY: set_mo_occupation 78 USE qs_mo_types, ONLY: get_mo_set,& 79 mo_set_p_type,& 80 mo_set_restrict 81 USE qs_mom_methods, ONLY: do_mom_guess 82 USE qs_rho_methods, ONLY: qs_rho_update_rho 83 USE qs_rho_types, ONLY: qs_rho_get,& 84 qs_rho_type 85 USE qs_scf_methods, ONLY: eigensolver,& 86 eigensolver_simple 87 USE qs_scf_types, ONLY: block_davidson_diag_method_nr,& 88 block_krylov_diag_method_nr,& 89 ot_diag_method_nr,& 90 qs_scf_env_type 91 USE qs_wf_history_methods, ONLY: wfi_update 92 USE scf_control_types, ONLY: scf_control_type 93 USE util, ONLY: sort 94 USE xtb_types, ONLY: get_xtb_atom_param,& 95 xtb_atom_type 96#include "./base/base_uses.f90" 97 98 IMPLICIT NONE 99 100 PRIVATE 101 102 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_initial_guess' 103 104 PUBLIC :: calculate_first_density_matrix, calculate_atomic_block_dm, calculate_mopac_dm 105 PUBLIC :: calculate_atomic_fock_matrix 106 107 TYPE atom_matrix_type 108 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat 109 END TYPE atom_matrix_type 110 111CONTAINS 112 113! ************************************************************************************************** 114!> \brief can use a variety of methods to come up with an initial 115!> density matrix and optionally an initial wavefunction 116!> \param scf_env SCF environment information 117!> \param qs_env QS environment 118!> \par History 119!> 03.2006 moved here from qs_scf [Joost VandeVondele] 120!> 06.2007 allow to skip the initial guess [jgh] 121!> 08.2014 kpoints [JGH] 122!> \note 123!> badly needs to be split in subroutines each doing one of the possible 124!> schemes 125! ************************************************************************************************** 126 SUBROUTINE calculate_first_density_matrix(scf_env, qs_env) 127 128 TYPE(qs_scf_env_type), POINTER :: scf_env 129 TYPE(qs_environment_type), POINTER :: qs_env 130 131 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_first_density_matrix', & 132 routineP = moduleN//':'//routineN 133 134 CHARACTER(LEN=default_path_length) :: file_name, filename 135 INTEGER :: atom_a, blk, density_guess, group, handle, homo, i, iatom, ic, icol, id_nr, & 136 ikind, irow, iseed(4), ispin, istart_col, istart_row, j, last_read, n, n_cols, n_rows, & 137 nao, natom, natoms, natoms_tmp, nblocks, nelectron, nmo, nmo_tmp, not_read, nsgf, nspin, & 138 nvec, ounit, qs_env_id, safe_density_guess, size_atomic_kind_set, z 139 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf 140 INTEGER, DIMENSION(2) :: nelectron_spin 141 INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nelec_kind, & 142 sort_kind 143 LOGICAL :: did_guess, do_kpoints, do_std_diag, & 144 exist, has_unit_metric, & 145 natom_mismatch, need_wm, ofgpw 146 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: buff, buff2 147 REAL(dp), DIMENSION(:, :), POINTER :: pdata 148 REAL(KIND=dp) :: checksum, eps, length, maxocc, occ, & 149 rscale, trps1, zeff 150 REAL(KIND=dp), DIMENSION(0:3) :: edftb 151 TYPE(atom_matrix_type), DIMENSION(:), POINTER :: pmat 152 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 153 TYPE(atomic_kind_type), POINTER :: atomic_kind 154 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: work1 155 TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, ao_mo_struct 156 TYPE(cp_fm_type), POINTER :: mo_coeff, moa, mob, ortho, sv, work2 157 TYPE(cp_logger_type), POINTER :: logger 158 TYPE(cp_para_env_type), POINTER :: para_env 159 TYPE(dbcsr_iterator_type) :: iter 160 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h_core_sparse, matrix_ks, p_rmpv, & 161 s_sparse 162 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h_kp, matrix_ks_kp, matrix_s_kp, & 163 rho_ao_kp 164 TYPE(dbcsr_type) :: mo_dbcsr, mo_tmp_dbcsr 165 TYPE(dft_control_type), POINTER :: dft_control 166 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 167 TYPE(kpoint_type), POINTER :: kpoints 168 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 169 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 170 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 171 TYPE(qs_kind_type), POINTER :: qs_kind 172 TYPE(qs_rho_type), POINTER :: rho 173 TYPE(scf_control_type), POINTER :: scf_control 174 TYPE(section_vals_type), POINTER :: dft_section, input, subsys_section 175 176 logger => cp_get_default_logger() 177 NULLIFY (atomic_kind, qs_kind, mo_coeff, sv, orb_basis_set, atomic_kind_set, & 178 qs_kind_set, particle_set, ortho, work2, work1, mo_array, s_sparse, & 179 scf_control, dft_control, p_rmpv, ortho, work2, work1, para_env, & 180 s_sparse, scf_control, dft_control, h_core_sparse, matrix_ks, rho) 181 NULLIFY (dft_section, input, subsys_section) 182 NULLIFY (matrix_s_kp, matrix_h_kp, matrix_ks_kp, rho_ao_kp) 183 NULLIFY (moa, mob) 184 NULLIFY (atom_list, elec_conf, kpoints) 185 edftb = 0.0_dp 186 187 CALL timeset(routineN, handle) 188 189 CALL get_qs_env(qs_env, & 190 atomic_kind_set=atomic_kind_set, & 191 qs_kind_set=qs_kind_set, & 192 particle_set=particle_set, & 193 mos=mo_array, & 194 matrix_s_kp=matrix_s_kp, & 195 matrix_h_kp=matrix_h_kp, & 196 matrix_ks_kp=matrix_ks_kp, & 197 input=input, & 198 scf_control=scf_control, & 199 id_nr=qs_env_id, & 200 dft_control=dft_control, & 201 has_unit_metric=has_unit_metric, & 202 do_kpoints=do_kpoints, & 203 kpoints=kpoints, & 204 rho=rho, & 205 nelectron_spin=nelectron_spin, & 206 para_env=para_env) 207 208 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 209 210 ! just initialize the first image, the other density are set to zero 211 DO ispin = 1, dft_control%nspins 212 DO ic = 1, SIZE(rho_ao_kp, 2) 213 CALL dbcsr_set(rho_ao_kp(ispin, ic)%matrix, 0.0_dp) 214 END DO 215 END DO 216 s_sparse => matrix_s_kp(:, 1) 217 h_core_sparse => matrix_h_kp(:, 1) 218 matrix_ks => matrix_ks_kp(:, 1) 219 p_rmpv => rho_ao_kp(:, 1) 220 221 work1 => scf_env%scf_work1 222 work2 => scf_env%scf_work2 223 ortho => scf_env%ortho 224 225 dft_section => section_vals_get_subs_vals(input, "DFT") 226 227 nspin = dft_control%nspins 228 ofgpw = dft_control%qs_control%ofgpw 229 density_guess = scf_control%density_guess 230 do_std_diag = .FALSE. 231 232 safe_density_guess = atomic_guess 233 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN 234 IF (density_guess == atomic_guess) density_guess = mopac_guess 235 safe_density_guess = mopac_guess 236 END IF 237 IF (dft_control%qs_control%xtb) THEN 238 IF (do_kpoints) THEN 239 IF (density_guess == atomic_guess) density_guess = mopac_guess 240 safe_density_guess = mopac_guess 241 ELSE 242 IF (density_guess == atomic_guess) density_guess = core_guess 243 safe_density_guess = core_guess 244 END IF 245 END IF 246 247 IF (scf_control%use_ot .AND. & 248 (.NOT. ((density_guess == random_guess) .OR. & 249 (density_guess == atomic_guess) .OR. & 250 (density_guess == core_guess) .OR. & 251 (density_guess == mopac_guess) .OR. & 252 (density_guess == sparse_guess) .OR. & 253 (((density_guess == restart_guess) .OR. & 254 (density_guess == history_guess)) .AND. & 255 (scf_control%level_shift == 0.0_dp))))) THEN 256 CALL cp_abort(__LOCATION__, & 257 "OT needs GUESS ATOMIC / CORE / RANDOM / SPARSE / RESTART / HISTORY RESTART: other options NYI") 258 END IF 259 260 ! if a restart was requested, check that the file exists, 261 ! if not we fall back to an atomic guess. No kidding, the file name should remain 262 ! in sync with read_mo_set 263 id_nr = 0 264 IF (density_guess == restart_guess) THEN 265 ! only check existence on I/O node, otherwise if file exists there but 266 ! not on compute nodes, everything goes crazy even though only I/O 267 ! node actually reads the file 268 IF (do_kpoints) THEN 269 IF (para_env%ionode) THEN 270 CALL wfn_restart_file_name(file_name, exist, dft_section, logger, kp=.TRUE.) 271 END IF 272 ELSE 273 IF (para_env%ionode) THEN 274 CALL wfn_restart_file_name(file_name, exist, dft_section, logger) 275 END IF 276 ENDIF 277 CALL mp_bcast(exist, para_env%source, para_env%group) 278 CALL mp_bcast(file_name, para_env%source, para_env%group) 279 IF (.NOT. exist) THEN 280 CALL cp_warn(__LOCATION__, & 281 "User requested to restart the wavefunction from the file named: "// & 282 TRIM(file_name)//". This file does not exist. Please check the existence of"// & 283 " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME."// & 284 " Calculation continues using ATOMIC GUESS. ") 285 density_guess = safe_density_guess 286 END IF 287 ELSE IF (density_guess == history_guess) THEN 288 IF (do_kpoints) THEN 289 CPABORT("calculate_first_density_matrix: history_guess not implemented for k-points") 290 ENDIF 291 IF (para_env%ionode) & 292 CALL wfn_restart_file_name(file_name, exist, dft_section, logger) 293 CALL mp_bcast(exist, para_env%source, para_env%group) 294 CALL mp_bcast(file_name, para_env%source, para_env%group) 295 nvec = qs_env%wf_history%memory_depth 296 not_read = nvec + 1 297 ! At this level we read the saved backup RESTART files.. 298 DO i = 1, nvec 299 j = i - 1 300 filename = TRIM(file_name) 301 IF (j /= 0) filename = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(j)) 302 IF (para_env%ionode) & 303 INQUIRE (FILE=filename, exist=exist) 304 CALL mp_bcast(exist, para_env%source, para_env%group) 305 IF ((.NOT. exist) .AND. (i < not_read)) THEN 306 not_read = i 307 END IF 308 END DO 309 IF (not_read == 1) THEN 310 density_guess = restart_guess 311 filename = TRIM(file_name) 312 IF (para_env%ionode) INQUIRE (FILE=filename, exist=exist) 313 CALL mp_bcast(exist, para_env%source, para_env%group) 314 IF (.NOT. exist) THEN 315 CALL cp_warn(__LOCATION__, & 316 "User requested to restart the wavefunction from a series of restart files named: "// & 317 TRIM(file_name)//" with extensions (.bak-n). These files do not exist."// & 318 " Even trying to switch to a plain restart wave-function failes because the"// & 319 " file named: "//TRIM(file_name)//" does not exist. Please check the existence of"// & 320 " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME. "// & 321 " Calculation continues using ATOMIC GUESS. ") 322 density_guess = safe_density_guess 323 END IF 324 END IF 325 last_read = not_read - 1 326 END IF 327 328 did_guess = .FALSE. 329 330 IF (density_guess == restart_guess) THEN 331 332 IF (do_kpoints) THEN 333 natoms = SIZE(particle_set) 334 CALL read_kpoints_restart(rho_ao_kp, kpoints, work1, & 335 natoms, para_env, id_nr, dft_section, natom_mismatch) 336 IF (natom_mismatch) density_guess = safe_density_guess 337 ELSE 338 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set) 339 CALL read_mo_set(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, & 340 id_nr=id_nr, multiplicity=dft_control%multiplicity, dft_section=dft_section, & 341 natom_mismatch=natom_mismatch) 342 343 IF (natom_mismatch) THEN 344 density_guess = safe_density_guess 345 ELSE 346 DO ispin = 1, nspin 347 IF (scf_control%level_shift /= 0.0_dp) THEN 348 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff) 349 CALL cp_fm_to_fm(mo_coeff, ortho) 350 END IF 351 352 ! make all nmo vectors present orthonormal 353 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 354 mo_coeff=mo_coeff, nmo=nmo, homo=homo) 355 356 IF (has_unit_metric) THEN 357 CALL make_basis_simple(mo_coeff, nmo) 358 ELSEIF (dft_control%smear) THEN 359 CALL make_basis_lowdin(vmatrix=mo_coeff, ncol=nmo, & 360 matrix_s=s_sparse(1)%matrix) 361 ELSE 362 ! ortho so that one can restart for different positions (basis sets?) 363 CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix) 364 ENDIF 365 ! only alpha spin is kept for restricted 366 IF (dft_control%restricted) EXIT 367 ENDDO 368 IF (dft_control%restricted) CALL mo_set_restrict(mo_array) 369 370 IF (.NOT. scf_control%diagonalization%mom) & 371 CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear) 372 373 DO ispin = 1, nspin 374 375 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr 376 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 377 mo_array(ispin)%mo_set%mo_coeff_b) !fm->dbcsr 378 ENDIF !fm->dbcsr 379 380 CALL calculate_density_matrix(mo_array(ispin)%mo_set, & 381 p_rmpv(ispin)%matrix) 382 ENDDO 383 ENDIF ! natom_mismatch 384 385 END IF 386 387 ! Maximum Overlap Method 388 IF (scf_control%diagonalization%mom) THEN 389 CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv) 390 ENDIF 391 392 did_guess = .TRUE. 393 END IF 394 395 IF (density_guess == history_guess) THEN 396 IF (not_read > 1) THEN 397 DO i = 1, last_read 398 j = last_read - i 399 CALL read_mo_set(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, & 400 id_nr=j, multiplicity=dft_control%multiplicity, & 401 dft_section=dft_section) 402 403 DO ispin = 1, nspin 404 IF (scf_control%level_shift /= 0.0_dp) THEN 405 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff) 406 CALL cp_fm_to_fm(mo_coeff, ortho) 407 END IF 408 409 ! make all nmo vectors present orthonormal 410 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo, homo=homo) 411 412 IF (has_unit_metric) THEN 413 CALL make_basis_simple(mo_coeff, nmo) 414 ELSE 415 ! ortho so that one can restart for different positions (basis sets?) 416 CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix) 417 ENDIF 418 ! only alpha spin is kept for restricted 419 IF (dft_control%restricted) EXIT 420 END DO 421 IF (dft_control%restricted) CALL mo_set_restrict(mo_array) 422 423 DO ispin = 1, nspin 424 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, & 425 smear=qs_env%scf_control%smear) 426 ENDDO 427 428 DO ispin = 1, nspin 429 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr 430 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 431 mo_array(ispin)%mo_set%mo_coeff_b) !fm->dbcsr 432 END IF !fm->dbcsr 433 CALL calculate_density_matrix(mo_array(ispin)%mo_set, p_rmpv(ispin)%matrix) 434 ENDDO 435 436 ! Write to extrapolation pipeline 437 CALL wfi_update(wf_history=qs_env%wf_history, qs_env=qs_env, dt=1.0_dp) 438 END DO 439 END IF 440 441 did_guess = .TRUE. 442 END IF 443 444 IF (density_guess == random_guess) THEN 445 446 DO ispin = 1, nspin 447 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 448 mo_coeff=mo_coeff, nmo=nmo) 449 CALL cp_fm_init_random(mo_coeff, nmo) 450 IF (has_unit_metric) THEN 451 CALL make_basis_simple(mo_coeff, nmo) 452 ELSE 453 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix) 454 ENDIF 455 ! only alpha spin is kept for restricted 456 IF (dft_control%restricted) EXIT 457 ENDDO 458 IF (dft_control%restricted) CALL mo_set_restrict(mo_array) 459 460 DO ispin = 1, nspin 461 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, & 462 smear=qs_env%scf_control%smear) 463 ENDDO 464 465 DO ispin = 1, nspin 466 467 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr 468 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 469 mo_array(ispin)%mo_set%mo_coeff_b) !fm->dbcsr 470 ENDIF !fm->dbcsr 471 472 CALL calculate_density_matrix(mo_array(ispin)%mo_set, p_rmpv(ispin)%matrix) 473 ENDDO 474 475 did_guess = .TRUE. 476 END IF 477 478 IF (density_guess == core_guess) THEN 479 480 IF (do_kpoints) THEN 481 CPABORT("calculate_first_density_matrix: core_guess not implemented for k-points") 482 ENDIF 483 484 IF (.NOT. ASSOCIATED(work1)) THEN 485 need_wm = .TRUE. 486 CPASSERT(.NOT. ASSOCIATED(work2)) 487 CPASSERT(.NOT. ASSOCIATED(ortho)) 488 ELSE 489 need_wm = .FALSE. 490 CPASSERT(ASSOCIATED(work2)) 491 END IF 492 493 IF (need_wm) THEN 494 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=moa) 495 CALL cp_fm_get_info(moa, matrix_struct=ao_mo_struct) 496 CALL cp_fm_struct_get(ao_mo_struct, nrow_global=nao, nrow_block=nblocks) 497 CALL cp_fm_struct_create(fmstruct=ao_ao_struct, & 498 nrow_block=nblocks, & 499 ncol_block=nblocks, & 500 nrow_global=nao, & 501 ncol_global=nao, & 502 template_fmstruct=ao_mo_struct) 503 ALLOCATE (work1(1)) 504 NULLIFY (work1(1)%matrix) 505 CALL cp_fm_create(work1(1)%matrix, ao_ao_struct) 506 CALL cp_fm_create(work2, ao_ao_struct) 507 CALL cp_fm_create(ortho, ao_ao_struct) 508 CALL copy_dbcsr_to_fm(matrix_s_kp(1, 1)%matrix, ortho) 509 CALL cp_fm_cholesky_decompose(ortho) 510 CALL cp_fm_struct_release(ao_ao_struct) 511 END IF 512 513 ispin = 1 514 ! Load core Hamiltonian into work matrix 515 CALL copy_dbcsr_to_fm(h_core_sparse(1)%matrix, work1(ispin)%matrix) 516 517 ! Diagonalize the core Hamiltonian matrix and retrieve a first set of 518 ! molecular orbitals (MOs) 519 IF (has_unit_metric) THEN 520 CALL eigensolver_simple(matrix_ks=work1(ispin)%matrix, & 521 mo_set=mo_array(ispin)%mo_set, & 522 work=work2, & 523 do_level_shift=.FALSE., & 524 level_shift=0.0_dp, & 525 use_jacobi=.FALSE., jacobi_threshold=0._dp) 526 ELSE 527 CALL eigensolver(matrix_ks_fm=work1(ispin)%matrix, & 528 mo_set=mo_array(ispin)%mo_set, & 529 ortho=ortho, & 530 work=work2, & 531 cholesky_method=scf_env%cholesky_method, & 532 do_level_shift=.FALSE., & 533 level_shift=0.0_dp, & 534 matrix_u_fm=null(), & 535 use_jacobi=.FALSE.) 536 END IF 537 538 ! Open shell case: copy alpha MOs to beta MOs 539 IF (nspin == 2) THEN 540 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=moa) 541 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff=mob, nmo=nmo) 542 CALL cp_fm_to_fm(moa, mob, nmo) 543 END IF 544 545 ! Build an initial density matrix (for each spin in the case of 546 ! an open shell calculation) from the first MOs set 547 DO ispin = 1, nspin 548 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, smear=scf_control%smear) 549 CALL calculate_density_matrix(mo_array(ispin)%mo_set, p_rmpv(ispin)%matrix) 550 END DO 551 552 ! release intermediate matrices 553 IF (need_wm) THEN 554 CALL cp_fm_release(ortho) 555 CALL cp_fm_release(work2) 556 CALL cp_fm_release(work1(1)%matrix) 557 DEALLOCATE (work1) 558 NULLIFY (work1, work2, ortho) 559 END IF 560 561 did_guess = .TRUE. 562 END IF 563 564 IF (density_guess == atomic_guess) THEN 565 566 subsys_section => section_vals_get_subs_vals(input, "SUBSYS") 567 ounit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log") 568 IF (ounit > 0) THEN 569 WRITE (UNIT=ounit, FMT="(/,(T2,A))") & 570 "Atomic guess: The first density matrix is obtained in terms of atomic orbitals", & 571 " and electronic configurations assigned to each atomic kind" 572 END IF 573 574 CALL calculate_atomic_block_dm(p_rmpv, s_sparse(1)%matrix, & 575 particle_set, atomic_kind_set, qs_kind_set, & 576 nspin, nelectron_spin, ounit, para_env) 577 578 DO ispin = 1, nspin 579 580 ! The orbital transformation method (OT) requires not only an 581 ! initial density matrix, but also an initial wavefunction (MO set) 582 IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN 583 ! get orbitals later 584 ELSE 585 IF (ASSOCIATED(scf_env%krylov_space)) do_std_diag = (scf_env%krylov_space%eps_std_diag > 0.0_dp) 586 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr .OR. & 587 (scf_env%method == block_krylov_diag_method_nr .AND. .NOT. do_std_diag) & 588 .OR. dft_control%do_admm .OR. scf_env%method == block_davidson_diag_method_nr) THEN 589 IF (dft_control%restricted .AND. (ispin == 2)) THEN 590 CALL mo_set_restrict(mo_array) 591 ELSE 592 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 593 mo_coeff=mo_coeff, & 594 nmo=nmo, nao=nao, homo=homo) 595 596 CALL cp_fm_set_all(mo_coeff, 0.0_dp) 597 CALL cp_fm_init_random(mo_coeff, nmo) 598 599 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV") 600 ! multiply times PS 601 IF (has_unit_metric) THEN 602 CALL cp_fm_to_fm(mo_coeff, sv) 603 ELSE 604 ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc)) 605 CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo) 606 END IF 607 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo) 608 609 CALL cp_fm_release(sv) 610 ! and ortho the result 611 IF (has_unit_metric) THEN 612 CALL make_basis_simple(mo_coeff, nmo) 613 ELSE 614 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix) 615 END IF 616 END IF 617 618 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, & 619 smear=qs_env%scf_control%smear) 620 621 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 622 mo_array(ispin)%mo_set%mo_coeff_b) !fm->dbcsr 623 624 CALL calculate_density_matrix(mo_array(ispin)%mo_set, & 625 p_rmpv(ispin)%matrix) 626 END IF 627 END IF 628 629 END DO 630 631 IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN 632 ! We fit a function to the square root of the density 633 CALL qs_rho_update_rho(rho, qs_env) 634 CPASSERT(1 == 0) 635! CALL cp_fm_create(sv,mo_coeff%matrix_struct,"SV") 636! DO ispin=1,nspin 637! CALL integrate_ppl_rspace(qs%rho%rho_r(ispin),qs_env) 638! CALL cp_cfm_solve(overlap,mos) 639! CALL get_mo_set(mo_set=mo_array(ispin)%mo_set,& 640! mo_coeff=mo_coeff, nmo=nmo, nao=nao) 641! CALL cp_fm_init_random(mo_coeff,nmo) 642! END DO 643! CALL cp_fm_release(sv) 644 END IF 645 646 IF (scf_control%diagonalization%mom) THEN 647 CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv) 648 ENDIF 649 650 CALL cp_print_key_finished_output(ounit, logger, subsys_section, & 651 "PRINT%KINDS") 652 653 did_guess = .TRUE. 654 END IF 655 656 IF (density_guess == sparse_guess) THEN 657 658 IF (ofgpw) CPABORT("SPARSE_GUESS not implemented for OFGPW") 659 IF (.NOT. scf_control%use_ot) CPABORT("OT needed!") 660 IF (do_kpoints) THEN 661 CPABORT("calculate_first_density_matrix: sparse_guess not implemented for k-points") 662 ENDIF 663 664 eps = 1.0E-5_dp 665 666 ounit = cp_logger_get_default_io_unit(logger) 667 group = para_env%group 668 natoms = SIZE(particle_set) 669 ALLOCATE (kind_of(natoms)) 670 ALLOCATE (first_sgf(natoms), last_sgf(natoms)) 671 672 checksum = dbcsr_checksum(s_sparse(1)%matrix) 673 i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL mp_sum(i, group) 674 IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum 675 CALL dbcsr_filter(s_sparse(1)%matrix, eps) 676 checksum = dbcsr_checksum(s_sparse(1)%matrix) 677 i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL mp_sum(i, group) 678 IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum 679 680 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, & 681 last_sgf=last_sgf) 682 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 683 684 ALLOCATE (pmat(SIZE(atomic_kind_set))) 685 686 rscale = 1._dp 687 IF (nspin == 2) rscale = 0.5_dp 688 DO ikind = 1, SIZE(atomic_kind_set) 689 atomic_kind => atomic_kind_set(ikind) 690 qs_kind => qs_kind_set(ikind) 691 NULLIFY (pmat(ikind)%mat) 692 CALL calculate_atomic_orbitals(atomic_kind, qs_kind, pmat=pmat(ikind)%mat) 693 NULLIFY (atomic_kind) 694 END DO 695 696 DO ispin = 1, nspin 697 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 698 maxocc=maxocc, & 699 nelectron=nelectron) 700 ! 701 CALL dbcsr_iterator_start(iter, p_rmpv(ispin)%matrix) 702 DO WHILE (dbcsr_iterator_blocks_left(iter)) 703 CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk) 704 ikind = kind_of(irow) 705 IF (icol .EQ. irow) THEN 706 IF (ispin == 1) THEN 707 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + & 708 pmat(ikind)%mat(:, :, 2)*rscale 709 ELSE 710 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - & 711 pmat(ikind)%mat(:, :, 2)*rscale 712 END IF 713 END IF 714 ENDDO 715 CALL dbcsr_iterator_stop(iter) 716 717 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix) 718 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix) 719 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix) 720 IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum 721 ! so far p needs to have the same sparsity as S 722 !CALL dbcsr_filter(p_rmpv(ispin)%matrix, eps) 723 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix) 724 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix) 725 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix) 726 IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum 727 728 CALL dbcsr_dot(p_rmpv(ispin)%matrix, s_sparse(1)%matrix, trps1) 729 rscale = REAL(nelectron, dp)/trps1 730 CALL dbcsr_scale(p_rmpv(ispin)%matrix, rscale) 731 732 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix) 733 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix) 734 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix) 735 IF (ounit > 0) WRITE (ounit, *) 'P occ', occ, ' checksum', checksum 736 ! 737 ! The orbital transformation method (OT) requires not only an 738 ! initial density matrix, but also an initial wavefunction (MO set) 739 IF (dft_control%restricted .AND. (ispin == 2)) THEN 740 CALL mo_set_restrict(mo_array) 741 ELSE 742 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 743 mo_coeff=mo_coeff, & 744 nmo=nmo, nao=nao, homo=homo) 745 CALL cp_fm_set_all(mo_coeff, 0.0_dp) 746 747 n = MAXVAL(last_sgf - first_sgf) + 1 748 size_atomic_kind_set = SIZE(atomic_kind_set) 749 750 ALLOCATE (buff(n, n), sort_kind(size_atomic_kind_set), & 751 nelec_kind(size_atomic_kind_set)) 752 ! 753 ! sort kind vs nbr electron 754 DO ikind = 1, size_atomic_kind_set 755 atomic_kind => atomic_kind_set(ikind) 756 qs_kind => qs_kind_set(ikind) 757 CALL get_atomic_kind(atomic_kind=atomic_kind, & 758 natom=natom, & 759 atom_list=atom_list, & 760 z=z) 761 CALL get_qs_kind(qs_kind, nsgf=nsgf, elec_conf=elec_conf, & 762 basis_set=orb_basis_set, zeff=zeff) 763 nelec_kind(ikind) = SUM(elec_conf) 764 ENDDO 765 CALL sort(nelec_kind, size_atomic_kind_set, sort_kind) 766 ! 767 ! a -very- naive sparse guess 768 nmo_tmp = nmo 769 natoms_tmp = natoms 770 istart_col = 1 771 iseed(1) = 4; iseed(2) = 3; iseed(3) = 2; iseed(4) = 1 ! set the seed for dlarnv 772 DO i = 1, size_atomic_kind_set 773 ikind = sort_kind(i) 774 atomic_kind => atomic_kind_set(ikind) 775 CALL get_atomic_kind(atomic_kind=atomic_kind, & 776 natom=natom, atom_list=atom_list) 777 DO iatom = 1, natom 778 ! 779 atom_a = atom_list(iatom) 780 istart_row = first_sgf(atom_a) 781 n_rows = last_sgf(atom_a) - first_sgf(atom_a) + 1 782 ! 783 ! compute the "potential" nbr of states for this atom 784 n_cols = MAX(INT(REAL(nmo_tmp, dp)/REAL(natoms_tmp, dp)), 1) 785 IF (n_cols .GT. n_rows) n_cols = n_rows 786 ! 787 nmo_tmp = nmo_tmp - n_cols 788 natoms_tmp = natoms_tmp - 1 789 IF (nmo_tmp .LT. 0 .OR. natoms_tmp .LT. 0) THEN 790 CPABORT("Wrong1!") 791 END IF 792 DO j = 1, n_cols 793 CALL dlarnv(1, iseed, n_rows, buff(1, j)) 794 ENDDO 795 CALL cp_fm_set_submatrix(mo_coeff, buff, istart_row, istart_col, & 796 n_rows, n_cols) 797 istart_col = istart_col + n_cols 798 ENDDO 799 ENDDO 800 801 IF (istart_col .LE. nmo) THEN 802 CPABORT("Wrong2!") 803 END IF 804 805 DEALLOCATE (buff, nelec_kind, sort_kind) 806 807 IF (.FALSE.) THEN 808 ALLOCATE (buff(nao, 1), buff2(nao, 1)) 809 DO i = 1, nmo 810 CALL cp_fm_get_submatrix(mo_coeff, buff, 1, i, nao, 1) 811 IF (SUM(buff**2) .LT. 1E-10_dp) THEN 812 WRITE (*, *) 'wrong', i, SUM(buff**2) 813 ENDIF 814 length = SQRT(DOT_PRODUCT(buff(:, 1), buff(:, 1))) 815 buff(:, :) = buff(:, :)/length 816 DO j = i + 1, nmo 817 CALL cp_fm_get_submatrix(mo_coeff, buff2, 1, j, nao, 1) 818 length = SQRT(DOT_PRODUCT(buff2(:, 1), buff2(:, 1))) 819 buff2(:, :) = buff2(:, :)/length 820 IF (ABS(DOT_PRODUCT(buff(:, 1), buff2(:, 1)) - 1.0_dp) .LT. 1E-10_dp) THEN 821 WRITE (*, *) 'wrong2', i, j, DOT_PRODUCT(buff(:, 1), buff2(:, 1)) 822 DO ikind = 1, nao 823 IF (ABS(mo_coeff%local_data(ikind, i)) .GT. 1e-10_dp) THEN 824 WRITE (*, *) 'c1', ikind, mo_coeff%local_data(ikind, i) 825 ENDIF 826 IF (ABS(mo_coeff%local_data(ikind, j)) .GT. 1e-10_dp) THEN 827 WRITE (*, *) 'c2', ikind, mo_coeff%local_data(ikind, j) 828 ENDIF 829 ENDDO 830 CPABORT("") 831 ENDIF 832 ENDDO 833 ENDDO 834 DEALLOCATE (buff, buff2) 835 836 ENDIF 837 ! 838 CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, mo_coeff, s_sparse(1)%matrix) 839 !CALL dbcsr_verify_matrix(mo_dbcsr) 840 checksum = dbcsr_checksum(mo_dbcsr) 841 842 occ = dbcsr_get_occupation(mo_dbcsr) 843 IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum 844 CALL dbcsr_filter(mo_dbcsr, eps) 845 !CALL dbcsr_verify_matrix(mo_dbcsr) 846 occ = dbcsr_get_occupation(mo_dbcsr) 847 checksum = dbcsr_checksum(mo_dbcsr) 848 IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum 849 ! 850 ! multiply times PS 851 IF (has_unit_metric) THEN 852 CPABORT("has_unit_metric will be removed soon") 853 END IF 854 ! 855 ! S*C 856 CALL dbcsr_copy(mo_tmp_dbcsr, mo_dbcsr, name="mo_tmp") 857 CALL dbcsr_multiply("N", "N", 1.0_dp, s_sparse(1)%matrix, mo_dbcsr, & 858 0.0_dp, mo_tmp_dbcsr, & 859 retain_sparsity=.TRUE.) 860 !CALL dbcsr_verify_matrix(mo_tmp_dbcsr) 861 checksum = dbcsr_checksum(mo_tmp_dbcsr) 862 occ = dbcsr_get_occupation(mo_tmp_dbcsr) 863 IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum 864 CALL dbcsr_filter(mo_tmp_dbcsr, eps) 865 !CALL dbcsr_verify_matrix(mo_tmp_dbcsr) 866 checksum = dbcsr_checksum(mo_tmp_dbcsr) 867 occ = dbcsr_get_occupation(mo_tmp_dbcsr) 868 IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum 869 ! 870 ! P*SC 871 ! the destroy is needed for the moment to avoid memory leaks ! 872 ! This one is not needed because _destroy takes care of zeroing. 873 CALL dbcsr_multiply("N", "N", 1.0_dp, p_rmpv(ispin)%matrix, & 874 mo_tmp_dbcsr, 0.0_dp, mo_dbcsr) 875 IF (.FALSE.) CALL dbcsr_verify_matrix(mo_dbcsr) 876 checksum = dbcsr_checksum(mo_dbcsr) 877 occ = dbcsr_get_occupation(mo_dbcsr) 878 IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum 879 CALL dbcsr_filter(mo_dbcsr, eps) 880 !CALL dbcsr_verify_matrix(mo_dbcsr) 881 checksum = dbcsr_checksum(mo_dbcsr) 882 occ = dbcsr_get_occupation(mo_dbcsr) 883 IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum 884 ! 885 CALL copy_dbcsr_to_fm(mo_dbcsr, mo_coeff) 886 887 CALL dbcsr_release(mo_dbcsr) 888 CALL dbcsr_release(mo_tmp_dbcsr) 889 890 ! and ortho the result 891 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix) 892 END IF 893 894 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, & 895 smear=qs_env%scf_control%smear) 896 897 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 898 mo_array(ispin)%mo_set%mo_coeff_b) !fm->dbcsr 899 900 CALL calculate_density_matrix(mo_array(ispin)%mo_set, & 901 p_rmpv(ispin)%matrix) 902 DO ikind = 1, SIZE(atomic_kind_set) 903 IF (ASSOCIATED(pmat(ikind)%mat)) THEN 904 DEALLOCATE (pmat(ikind)%mat) 905 END IF 906 END DO 907 END DO 908 909 DEALLOCATE (pmat) 910 911 DEALLOCATE (kind_of) 912 913 DEALLOCATE (first_sgf, last_sgf) 914 915 did_guess = .TRUE. 916 END IF 917 IF (density_guess == mopac_guess) THEN 918 919 CALL calculate_mopac_dm(p_rmpv, s_sparse(1)%matrix, has_unit_metric, dft_control, & 920 particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env) 921 922 DO ispin = 1, nspin 923 ! The orbital transformation method (OT) requires not only an 924 ! initial density matrix, but also an initial wavefunction (MO set) 925 IF (ASSOCIATED(scf_env%krylov_space)) do_std_diag = (scf_env%krylov_space%eps_std_diag > 0.0_dp) 926 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr .OR. & 927 (scf_env%method == block_krylov_diag_method_nr .AND. .NOT. do_std_diag)) THEN 928 IF (dft_control%restricted .AND. (ispin == 2)) THEN 929 CALL mo_set_restrict(mo_array) 930 ELSE 931 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 932 mo_coeff=mo_coeff, & 933 nmo=nmo, homo=homo) 934 CALL cp_fm_init_random(mo_coeff, nmo) 935 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV") 936 ! multiply times PS 937 IF (has_unit_metric) THEN 938 CALL cp_fm_to_fm(mo_coeff, sv) 939 ELSE 940 CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo) 941 END IF 942 ! here we could easily multiply with the diag that we actually have replicated already 943 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo) 944 CALL cp_fm_release(sv) 945 ! and ortho the result 946 IF (has_unit_metric) THEN 947 CALL make_basis_simple(mo_coeff, nmo) 948 ELSE 949 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix) 950 END IF 951 END IF 952 953 CALL set_mo_occupation(mo_set=mo_array(ispin)%mo_set, & 954 smear=qs_env%scf_control%smear) 955 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_set%mo_coeff, & 956 mo_array(ispin)%mo_set%mo_coeff_b) 957 958 CALL calculate_density_matrix(mo_array(ispin)%mo_set, & 959 p_rmpv(ispin)%matrix) 960 END IF 961 END DO 962 963 did_guess = .TRUE. 964 END IF 965 966 IF (density_guess == no_guess) THEN 967 did_guess = .TRUE. 968 END IF 969 970 IF (.NOT. did_guess) THEN 971 CPABORT("An invalid keyword for the initial density guess was specified") 972 END IF 973 974 CALL timestop(handle) 975 976 END SUBROUTINE calculate_first_density_matrix 977 978! ************************************************************************************************** 979!> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities. 980!> \param pmatrix ... 981!> \param matrix_s ... 982!> \param particle_set ... 983!> \param atomic_kind_set ... 984!> \param qs_kind_set ... 985!> \param nspin ... 986!> \param nelectron_spin ... 987!> \param ounit ... 988!> \param para_env ... 989! ************************************************************************************************** 990 SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, particle_set, atomic_kind_set, & 991 qs_kind_set, nspin, nelectron_spin, ounit, para_env) 992 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmatrix 993 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s 994 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 995 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 996 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 997 INTEGER, INTENT(IN) :: nspin 998 INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin 999 INTEGER, INTENT(IN) :: ounit 1000 TYPE(cp_para_env_type) :: para_env 1001 1002 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm', & 1003 routineP = moduleN//':'//routineN 1004 1005 INTEGER :: blk, group, handle, icol, ikind, irow, & 1006 ispin, natom, nc, nkind, nocc(2) 1007 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of 1008 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nok 1009 REAL(dp), DIMENSION(:, :), POINTER :: pdata 1010 REAL(KIND=dp) :: rds, rscale, trps1 1011 TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: pmat 1012 TYPE(atomic_kind_type), POINTER :: atomic_kind 1013 TYPE(dbcsr_iterator_type) :: iter 1014 TYPE(dbcsr_type), POINTER :: matrix_p 1015 TYPE(qs_kind_type), POINTER :: qs_kind 1016 1017 CALL timeset(routineN, handle) 1018 1019 natom = SIZE(particle_set) 1020 nkind = SIZE(atomic_kind_set) 1021 ALLOCATE (kind_of(natom)) 1022 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 1023 ALLOCATE (pmat(nkind)) 1024 ALLOCATE (nok(2, nkind)) 1025 1026 ! precompute the atomic blocks corresponding to spherical atoms 1027 DO ikind = 1, nkind 1028 atomic_kind => atomic_kind_set(ikind) 1029 qs_kind => qs_kind_set(ikind) 1030 NULLIFY (pmat(ikind)%mat) 1031 IF (ounit > 0) THEN 1032 WRITE (UNIT=ounit, FMT="(/,T2,A)") & 1033 "Guess for atomic kind: "//TRIM(atomic_kind%name) 1034 END IF 1035 CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, & 1036 pmat=pmat(ikind)%mat, nocc=nocc) 1037 nok(1:2, ikind) = nocc(1:2) 1038 END DO 1039 1040 rscale = 1.0_dp 1041 IF (nspin == 2) rscale = 0.5_dp 1042 1043 DO ispin = 1, nspin 1044 IF ((ounit > 0) .AND. (nspin > 1)) THEN 1045 WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin 1046 END IF 1047 1048 matrix_p => pmatrix(ispin)%matrix 1049 CALL dbcsr_set(matrix_p, 0.0_dp) 1050 1051 nocc(ispin) = 0 1052 CALL dbcsr_iterator_start(iter, matrix_p) 1053 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1054 CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk) 1055 ikind = kind_of(irow) 1056 IF (icol .EQ. irow) THEN 1057 IF (ispin == 1) THEN 1058 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + & 1059 pmat(ikind)%mat(:, :, 2)*rscale 1060 ELSE 1061 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - & 1062 pmat(ikind)%mat(:, :, 2)*rscale 1063 END IF 1064 nocc(ispin) = nocc(ispin) + nok(ispin, ikind) 1065 END IF 1066 ENDDO 1067 CALL dbcsr_iterator_stop(iter) 1068 1069 CALL dbcsr_dot(matrix_p, matrix_s, trps1) 1070 rds = 0.0_dp 1071 ! could be a ghost-atoms-only simulation 1072 IF (nelectron_spin(ispin) > 0) THEN 1073 rds = REAL(nelectron_spin(ispin), dp)/trps1 1074 END IF 1075 CALL dbcsr_scale(matrix_p, rds) 1076 1077 IF (ounit > 0) THEN 1078 IF (nspin > 1) THEN 1079 WRITE (UNIT=ounit, FMT="(T2,A,I1)") & 1080 "Re-scaling the density matrix to get the right number of electrons for spin ", ispin 1081 ELSE 1082 WRITE (UNIT=ounit, FMT="(T2,A)") & 1083 "Re-scaling the density matrix to get the right number of electrons" 1084 END IF 1085 WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor" 1086 WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds 1087 END IF 1088 1089 IF (nspin > 1) THEN 1090 group = para_env%group 1091 CALL mp_sum(nocc, group) 1092 IF (nelectron_spin(ispin) > nocc(ispin)) THEN 1093 rds = 0.99_dp 1094 CALL dbcsr_scale(matrix_p, rds) 1095 rds = (1.0_dp - rds)*nelectron_spin(ispin) 1096 CALL dbcsr_get_info(matrix_p, nfullcols_total=nc) 1097 rds = rds/REAL(nc, KIND=dp) 1098 CALL dbcsr_add_on_diag(matrix_p, rds) 1099 IF (ounit > 0) THEN 1100 WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") & 1101 "More MOs than initial guess orbitals detected", & 1102 "Add constant to diagonal elements ", rds 1103 END IF 1104 END IF 1105 END IF 1106 1107 END DO 1108 1109 DO ikind = 1, nkind 1110 IF (ASSOCIATED(pmat(ikind)%mat)) THEN 1111 DEALLOCATE (pmat(ikind)%mat) 1112 END IF 1113 END DO 1114 DEALLOCATE (pmat) 1115 1116 DEALLOCATE (kind_of, nok) 1117 1118 CALL timestop(handle) 1119 1120 END SUBROUTINE calculate_atomic_block_dm 1121 1122! ************************************************************************************************** 1123!> \brief returns a block diagonal fock matrix. 1124!> \param matrix_f ... 1125!> \param particle_set ... 1126!> \param atomic_kind_set ... 1127!> \param qs_kind_set ... 1128!> \param ounit ... 1129! ************************************************************************************************** 1130 SUBROUTINE calculate_atomic_fock_matrix(matrix_f, particle_set, atomic_kind_set, & 1131 qs_kind_set, ounit) 1132 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_f 1133 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1134 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1135 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1136 INTEGER, INTENT(IN) :: ounit 1137 1138 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_fock_matrix', & 1139 routineP = moduleN//':'//routineN 1140 1141 INTEGER :: handle, icol, ikind, irow 1142 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of 1143 REAL(dp), DIMENSION(:, :), POINTER :: block 1144 TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: fmat 1145 TYPE(atomic_kind_type), POINTER :: atomic_kind 1146 TYPE(dbcsr_iterator_type) :: iter 1147 TYPE(qs_kind_type), POINTER :: qs_kind 1148 1149 CALL timeset(routineN, handle) 1150 1151 ALLOCATE (kind_of(SIZE(particle_set))) 1152 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) 1153 ALLOCATE (fmat(SIZE(atomic_kind_set))) 1154 1155 ! precompute the atomic blocks for each atomic-kind 1156 DO ikind = 1, SIZE(atomic_kind_set) 1157 atomic_kind => atomic_kind_set(ikind) 1158 qs_kind => qs_kind_set(ikind) 1159 NULLIFY (fmat(ikind)%mat) 1160 IF (ounit > 0) WRITE (UNIT=ounit, FMT="(/,T2,A)") & 1161 "Calculating atomic Fock matrix for atomic kind: "//TRIM(atomic_kind%name) 1162 1163 !Currently only ispin=1 is supported 1164 CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, & 1165 fmat=fmat(ikind)%mat) 1166 END DO 1167 1168 ! zero result matrix 1169 CALL dbcsr_set(matrix_f, 0.0_dp) 1170 1171 ! copy precomputed blocks onto diagonal of result matrix 1172 CALL dbcsr_iterator_start(iter, matrix_f) 1173 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1174 CALL dbcsr_iterator_next_block(iter, irow, icol, block) 1175 ikind = kind_of(irow) 1176 IF (icol .EQ. irow) block(:, :) = fmat(ikind)%mat(:, :, 1) 1177 ENDDO 1178 CALL dbcsr_iterator_stop(iter) 1179 1180 ! cleanup 1181 DO ikind = 1, SIZE(atomic_kind_set) 1182 DEALLOCATE (fmat(ikind)%mat) 1183 END DO 1184 DEALLOCATE (fmat, kind_of) 1185 1186 CALL timestop(handle) 1187 1188 END SUBROUTINE calculate_atomic_fock_matrix 1189 1190! ************************************************************************************************** 1191!> \brief returns a block diagonal density matrix. Blocks correspond to the mopac initial guess. 1192!> \param pmat ... 1193!> \param matrix_s ... 1194!> \param has_unit_metric ... 1195!> \param dft_control ... 1196!> \param particle_set ... 1197!> \param atomic_kind_set ... 1198!> \param qs_kind_set ... 1199!> \param nspin ... 1200!> \param nelectron_spin ... 1201!> \param para_env ... 1202! ************************************************************************************************** 1203 SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, & 1204 dft_control, particle_set, atomic_kind_set, qs_kind_set, & 1205 nspin, nelectron_spin, para_env) 1206 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmat 1207 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s 1208 LOGICAL :: has_unit_metric 1209 TYPE(dft_control_type), POINTER :: dft_control 1210 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1211 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1212 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1213 INTEGER, INTENT(IN) :: nspin 1214 INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin 1215 TYPE(cp_para_env_type) :: para_env 1216 1217 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_mopac_dm', & 1218 routineP = moduleN//':'//routineN 1219 1220 INTEGER :: atom_a, group, handle, iatom, ikind, & 1221 iset, isgf, isgfa, ishell, ispin, la, & 1222 maxl, maxll, na, nao, natom, ncount, & 1223 nset, nsgf, z 1224 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf 1225 INTEGER, DIMENSION(25) :: laox, naox 1226 INTEGER, DIMENSION(5) :: occupation 1227 INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nshell 1228 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, l, last_sgfa 1229 LOGICAL :: has_pot 1230 REAL(KIND=dp) :: maxocc, my_sum, nelec, occ, paa, rscale, & 1231 trps1, trps2, yy, zeff 1232 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: econf, pdiag, sdiag 1233 REAL(KIND=dp), DIMENSION(0:3) :: edftb 1234 TYPE(all_potential_type), POINTER :: all_potential 1235 TYPE(dbcsr_type), POINTER :: matrix_p 1236 TYPE(gth_potential_type), POINTER :: gth_potential 1237 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 1238 TYPE(sgp_potential_type), POINTER :: sgp_potential 1239 TYPE(xtb_atom_type), POINTER :: xtb_kind 1240 1241 CALL timeset(routineN, handle) 1242 1243 DO ispin = 1, nspin 1244 matrix_p => pmat(ispin)%matrix 1245 CALL dbcsr_set(matrix_p, 0.0_dp) 1246 END DO 1247 1248 group = para_env%group 1249 natom = SIZE(particle_set) 1250 nao = dbcsr_nfullrows_total(pmat(1)%matrix) 1251 IF (nspin == 1) THEN 1252 maxocc = 2.0_dp 1253 ELSE 1254 maxocc = 1.0_dp 1255 ENDIF 1256 1257 ALLOCATE (first_sgf(natom)) 1258 1259 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf) 1260 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxl) 1261 1262 ALLOCATE (econf(0:maxl)) 1263 1264 ALLOCATE (pdiag(nao)) 1265 pdiag(:) = 0.0_dp 1266 1267 ALLOCATE (sdiag(nao)) 1268 1269 sdiag(:) = 0.0_dp 1270 IF (has_unit_metric) THEN 1271 sdiag(:) = 1.0_dp 1272 ELSE 1273 CALL dbcsr_get_diag(matrix_s, sdiag) 1274 CALL mp_sum(sdiag, group) 1275 END IF 1276 1277 ncount = 0 1278 trps1 = 0.0_dp 1279 trps2 = 0.0_dp 1280 pdiag(:) = 0.0_dp 1281 1282 IF (SUM(nelectron_spin) /= 0) THEN 1283 DO ikind = 1, SIZE(atomic_kind_set) 1284 1285 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list) 1286 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, & 1287 all_potential=all_potential, & 1288 gth_potential=gth_potential, & 1289 sgp_potential=sgp_potential) 1290 has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential) 1291 1292 IF (dft_control%qs_control%dftb) THEN 1293 CALL get_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter, & 1294 lmax=maxll, occupation=edftb) 1295 maxll = MIN(maxll, maxl) 1296 econf(0:maxl) = edftb(0:maxl) 1297 ELSEIF (dft_control%qs_control%xtb) THEN 1298 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) 1299 CALL get_xtb_atom_param(xtb_kind, natorb=nsgf, nao=naox, lao=laox, occupation=occupation) 1300 ELSEIF (has_pot) THEN 1301 CALL get_atomic_kind(atomic_kind_set(ikind), z=z) 1302 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, elec_conf=elec_conf, zeff=zeff) 1303 maxll = MIN(SIZE(elec_conf) - 1, maxl) 1304 econf(:) = 0.0_dp 1305 econf(0:maxll) = 0.5_dp*maxocc*REAL(elec_conf(0:maxll), dp) 1306 ELSE 1307 CYCLE 1308 END IF 1309 1310 ! MOPAC TYPE GUESS 1311 IF (dft_control%qs_control%dftb) THEN 1312 DO iatom = 1, natom 1313 atom_a = atom_list(iatom) 1314 isgfa = first_sgf(atom_a) 1315 DO la = 0, maxll 1316 SELECT CASE (la) 1317 CASE (0) 1318 pdiag(isgfa) = econf(0) 1319 CASE (1) 1320 pdiag(isgfa + 1) = econf(1)/3._dp 1321 pdiag(isgfa + 2) = econf(1)/3._dp 1322 pdiag(isgfa + 3) = econf(1)/3._dp 1323 CASE (2) 1324 pdiag(isgfa + 4) = econf(2)/5._dp 1325 pdiag(isgfa + 5) = econf(2)/5._dp 1326 pdiag(isgfa + 6) = econf(2)/5._dp 1327 pdiag(isgfa + 7) = econf(2)/5._dp 1328 pdiag(isgfa + 8) = econf(2)/5._dp 1329 CASE (3) 1330 pdiag(isgfa + 9) = econf(3)/7._dp 1331 pdiag(isgfa + 10) = econf(3)/7._dp 1332 pdiag(isgfa + 11) = econf(3)/7._dp 1333 pdiag(isgfa + 12) = econf(3)/7._dp 1334 pdiag(isgfa + 13) = econf(3)/7._dp 1335 pdiag(isgfa + 14) = econf(3)/7._dp 1336 pdiag(isgfa + 15) = econf(3)/7._dp 1337 CASE DEFAULT 1338 CPABORT("") 1339 END SELECT 1340 END DO 1341 END DO 1342 ELSEIF (dft_control%qs_control%xtb) THEN 1343 DO iatom = 1, natom 1344 atom_a = atom_list(iatom) 1345 isgfa = first_sgf(atom_a) 1346 DO isgf = 1, nsgf 1347 na = naox(isgf) 1348 la = laox(isgf) 1349 occ = REAL(occupation(na), dp)/REAL(2*la + 1, dp) 1350 pdiag(isgfa + isgf - 1) = occ 1351 END DO 1352 END DO 1353 ELSEIF (dft_control%qs_control%semi_empirical) THEN 1354 yy = REAL(dft_control%charge, KIND=dp)/REAL(nao, KIND=dp) 1355 DO iatom = 1, natom 1356 atom_a = atom_list(iatom) 1357 isgfa = first_sgf(atom_a) 1358 SELECT CASE (nsgf) 1359 CASE (1) ! s-basis 1360 pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc 1361 CASE (4) ! sp-basis 1362 IF (z == 1) THEN 1363 ! special case: hydrogen with sp basis 1364 pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc 1365 pdiag(isgfa + 1) = 0._dp 1366 pdiag(isgfa + 2) = 0._dp 1367 pdiag(isgfa + 3) = 0._dp 1368 ELSE 1369 pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1370 pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1371 pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1372 pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1373 END IF 1374 CASE (9) ! spd-basis 1375 IF (z < 21 .OR. z > 30 .AND. z < 39 .OR. z > 48 .AND. z < 57) THEN 1376 ! Main Group Element: The "d" shell is formally empty. 1377 pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1378 pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1379 pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1380 pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc 1381 pdiag(isgfa + 4) = (-yy)*0.5_dp*maxocc 1382 pdiag(isgfa + 5) = (-yy)*0.5_dp*maxocc 1383 pdiag(isgfa + 6) = (-yy)*0.5_dp*maxocc 1384 pdiag(isgfa + 7) = (-yy)*0.5_dp*maxocc 1385 pdiag(isgfa + 8) = (-yy)*0.5_dp*maxocc 1386 ELSE IF (z < 99) THEN 1387 my_sum = zeff - 9.0_dp*yy 1388 ! First, put 2 electrons in the 's' shell 1389 pdiag(isgfa) = (MAX(0.0_dp, MIN(my_sum, 2.0_dp)))*0.5_dp*maxocc 1390 my_sum = my_sum - 2.0_dp 1391 IF (my_sum > 0.0_dp) THEN 1392 ! Now put as many electrons as possible into the 'd' shell 1393 pdiag(isgfa + 4) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc 1394 pdiag(isgfa + 5) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc 1395 pdiag(isgfa + 6) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc 1396 pdiag(isgfa + 7) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc 1397 pdiag(isgfa + 8) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc 1398 my_sum = MAX(0.0_dp, my_sum - 10.0_dp) 1399 ! Put the remaining electrons in the 'p' shell 1400 pdiag(isgfa + 1) = (my_sum/3.0_dp)*0.5_dp*maxocc 1401 pdiag(isgfa + 2) = (my_sum/3.0_dp)*0.5_dp*maxocc 1402 pdiag(isgfa + 3) = (my_sum/3.0_dp)*0.5_dp*maxocc 1403 END IF 1404 END IF 1405 CASE DEFAULT 1406 CPABORT("") 1407 END SELECT 1408 END DO 1409 ELSE 1410 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 1411 nset=nset, & 1412 nshell=nshell, & 1413 l=l, & 1414 first_sgf=first_sgfa, & 1415 last_sgf=last_sgfa) 1416 1417 DO iset = 1, nset 1418 DO ishell = 1, nshell(iset) 1419 la = l(ishell, iset) 1420 nelec = maxocc*REAL(2*la + 1, dp) 1421 IF (econf(la) > 0.0_dp) THEN 1422 IF (econf(la) >= nelec) THEN 1423 paa = maxocc 1424 econf(la) = econf(la) - nelec 1425 ELSE 1426 paa = maxocc*econf(la)/nelec 1427 econf(la) = 0.0_dp 1428 ncount = ncount + NINT(nelec/maxocc) 1429 END IF 1430 DO isgfa = first_sgfa(ishell, iset), last_sgfa(ishell, iset) 1431 DO iatom = 1, natom 1432 atom_a = atom_list(iatom) 1433 isgf = first_sgf(atom_a) + isgfa - 1 1434 pdiag(isgf) = paa 1435 IF (paa == maxocc) THEN 1436 trps1 = trps1 + paa*sdiag(isgf) 1437 ELSE 1438 trps2 = trps2 + paa*sdiag(isgf) 1439 END IF 1440 END DO 1441 END DO 1442 END IF 1443 END DO ! ishell 1444 END DO ! iset 1445 END IF 1446 END DO ! ikind 1447 1448 IF (trps2 == 0.0_dp) THEN 1449 DO isgf = 1, nao 1450 IF (sdiag(isgf) > 0.0_dp) pdiag(isgf) = pdiag(isgf)/sdiag(isgf) 1451 END DO 1452 DO ispin = 1, nspin 1453 IF (nelectron_spin(ispin) /= 0) THEN 1454 matrix_p => pmat(ispin)%matrix 1455 CALL dbcsr_set_diag(matrix_p, pdiag) 1456 END IF 1457 END DO 1458 ELSE 1459 DO ispin = 1, nspin 1460 IF (nelectron_spin(ispin) /= 0) THEN 1461 rscale = (REAL(nelectron_spin(ispin), dp) - trps1)/trps2 1462 DO isgf = 1, nao 1463 IF (pdiag(isgf) < maxocc) pdiag(isgf) = rscale*pdiag(isgf) 1464 END DO 1465 matrix_p => pmat(ispin)%matrix 1466 CALL dbcsr_set_diag(matrix_p, pdiag) 1467 DO isgf = 1, nao 1468 IF (pdiag(isgf) < maxocc) pdiag(isgf) = pdiag(isgf)/rscale 1469 END DO 1470 END IF 1471 END DO 1472 END IF 1473 END IF 1474 1475 DEALLOCATE (econf) 1476 1477 DEALLOCATE (first_sgf) 1478 1479 DEALLOCATE (pdiag) 1480 1481 DEALLOCATE (sdiag) 1482 1483 CALL timestop(handle) 1484 1485 END SUBROUTINE calculate_mopac_dm 1486 1487END MODULE qs_initial_guess 1488