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