(git:ed6f26b)
Loading...
Searching...
No Matches
qs_initial_guess.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to somehow generate an initial guess
10!> \par History
11!> 2006.03 Moved here from qs_scf.F [Joost VandeVondele]
12! **************************************************************************************************
21 USE cp_dbcsr_api, ONLY: &
27 dbcsr_dot,&
39 USE cp_fm_types, ONLY: &
51 USE hfx_types, ONLY: hfx_type
52 USE input_constants, ONLY: &
55 USE input_cp2k_hfx, ONLY: ri_mo
59 USE kinds, ONLY: default_path_length,&
60 dp
62 USE kpoint_types, ONLY: kpoint_type
72 USE qs_kind_types, ONLY: get_qs_kind,&
81 USE qs_mo_types, ONLY: get_mo_set,&
87 USE qs_rho_types, ONLY: qs_rho_get,&
89 USE qs_scf_methods, ONLY: eigensolver,&
98 USE util, ONLY: sort
99 USE xtb_types, ONLY: get_xtb_atom_param,&
101#include "./base/base_uses.f90"
102
103 IMPLICIT NONE
104
105 PRIVATE
106
107 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_initial_guess'
108
111
112 TYPE atom_matrix_type
113 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat => null()
114 END TYPE atom_matrix_type
115
116CONTAINS
117
118! **************************************************************************************************
119!> \brief can use a variety of methods to come up with an initial
120!> density matrix and optionally an initial wavefunction
121!> \param scf_env SCF environment information
122!> \param qs_env QS environment
123!> \par History
124!> 03.2006 moved here from qs_scf [Joost VandeVondele]
125!> 06.2007 allow to skip the initial guess [jgh]
126!> 08.2014 kpoints [JGH]
127!> 10.2019 tot_corr_zeff, switch_surf_dip [SGh]
128!> \note
129!> badly needs to be split in subroutines each doing one of the possible
130!> schemes
131! **************************************************************************************************
132 SUBROUTINE calculate_first_density_matrix(scf_env, qs_env)
133
134 TYPE(qs_scf_env_type), POINTER :: scf_env
135 TYPE(qs_environment_type), POINTER :: qs_env
136
137 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_first_density_matrix'
138
139 CHARACTER(LEN=default_path_length) :: file_name, filename
140 INTEGER :: atom_a, density_guess, handle, homo, i, iatom, ic, icol, id_nr, ikind, irow, &
141 iseed(4), ispin, istart_col, istart_row, j, last_read, n, n_cols, n_rows, nao, natom, &
142 natoms, natoms_tmp, nblocks, nelectron, nmo, nmo_tmp, not_read, nsgf, nspin, nvec, ounit, &
143 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, do_std_diag, exist, has_unit_metric, &
149 natom_mismatch, need_mos, need_wm, ofgpw, owns_ortho, print_history_log, print_log
150 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: buff, buff2
151 REAL(dp), DIMENSION(:, :), POINTER :: pdata
152 REAL(kind=dp) :: checksum, eps, length, maxocc, occ, &
153 rscale, tot_corr_zeff, trps1, zeff
154 REAL(kind=dp), DIMENSION(0:3) :: edftb
155 TYPE(atom_matrix_type), DIMENSION(:), POINTER :: pmat
156 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
157 TYPE(atomic_kind_type), POINTER :: atomic_kind
158 TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, ao_mo_struct
159 TYPE(cp_fm_type) :: sv
160 TYPE(cp_fm_type), DIMENSION(:), POINTER :: work1
161 TYPE(cp_fm_type), POINTER :: mo_coeff, moa, mob, ortho, work2
162 TYPE(cp_logger_type), POINTER :: logger
163 TYPE(dbcsr_iterator_type) :: iter
164 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h_core_sparse, matrix_ks, p_rmpv, &
165 s_sparse
166 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h_kp, matrix_ks_kp, matrix_s_kp, &
167 rho_ao_kp
168 TYPE(dbcsr_type) :: mo_dbcsr, mo_tmp_dbcsr
169 TYPE(dft_control_type), POINTER :: dft_control
170 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
171 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
172 TYPE(kpoint_type), POINTER :: kpoints
173 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array, mos_last_converged
174 TYPE(mp_para_env_type), POINTER :: para_env
175 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
177 TYPE(qs_kind_type), POINTER :: qs_kind
178 TYPE(qs_rho_type), POINTER :: rho
179 TYPE(scf_control_type), POINTER :: scf_control
180 TYPE(section_vals_type), POINTER :: dft_section, input, subsys_section
181
182 logger => cp_get_default_logger()
183 NULLIFY (atomic_kind, qs_kind, mo_coeff, orb_basis_set, atomic_kind_set, &
184 qs_kind_set, particle_set, ortho, work2, work1, mo_array, s_sparse, &
185 scf_control, dft_control, p_rmpv, para_env, h_core_sparse, matrix_ks, rho, &
186 mos_last_converged)
187 NULLIFY (dft_section, input, subsys_section)
188 NULLIFY (matrix_s_kp, matrix_h_kp, matrix_ks_kp, rho_ao_kp)
189 NULLIFY (moa, mob)
190 NULLIFY (atom_list, elec_conf, kpoints)
191 edftb = 0.0_dp
192 tot_corr_zeff = 0.0_dp
193
194 CALL timeset(routinen, handle)
195
196 CALL get_qs_env(qs_env, &
197 atomic_kind_set=atomic_kind_set, &
198 qs_kind_set=qs_kind_set, &
199 particle_set=particle_set, &
200 mos=mo_array, &
201 matrix_s_kp=matrix_s_kp, &
202 matrix_h_kp=matrix_h_kp, &
203 matrix_ks_kp=matrix_ks_kp, &
204 input=input, &
205 scf_control=scf_control, &
206 dft_control=dft_control, &
207 has_unit_metric=has_unit_metric, &
208 do_kpoints=do_kpoints, &
209 kpoints=kpoints, &
210 rho=rho, &
211 nelectron_spin=nelectron_spin, &
212 para_env=para_env, &
213 x_data=x_data)
214
215 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
216
217 IF (dft_control%switch_surf_dip) THEN
218 CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
219 END IF
220
221 ! just initialize the first image, the other density are set to zero
222 DO ispin = 1, dft_control%nspins
223 DO ic = 1, SIZE(rho_ao_kp, 2)
224 CALL dbcsr_set(rho_ao_kp(ispin, ic)%matrix, 0.0_dp)
225 END DO
226 END DO
227 s_sparse => matrix_s_kp(:, 1)
228 h_core_sparse => matrix_h_kp(:, 1)
229 matrix_ks => matrix_ks_kp(:, 1)
230 p_rmpv => rho_ao_kp(:, 1)
231
232 work1 => scf_env%scf_work1
233 work2 => scf_env%scf_work2
234 ortho => scf_env%ortho
235
236 dft_section => section_vals_get_subs_vals(input, "DFT")
237
238 nspin = dft_control%nspins
239 ofgpw = dft_control%qs_control%ofgpw
240 density_guess = scf_control%density_guess
241 do_std_diag = .false.
242
243 do_hfx_ri_mo = .false.
244 IF (ASSOCIATED(x_data)) THEN
245 IF (x_data(1, 1)%do_hfx_ri) THEN
246 IF (x_data(1, 1)%ri_data%flavor == ri_mo) do_hfx_ri_mo = .true.
247 END IF
248 END IF
249
250 IF (ASSOCIATED(scf_env%krylov_space)) do_std_diag = (scf_env%krylov_space%eps_std_diag > 0.0_dp)
251
252 need_mos = scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr .OR. &
253 (scf_env%method == block_krylov_diag_method_nr .AND. .NOT. do_std_diag) &
254 .OR. dft_control%do_admm .OR. scf_env%method == block_davidson_diag_method_nr &
255 .OR. do_hfx_ri_mo
256
257 safe_density_guess = atomic_guess
258 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN
259 IF (density_guess == atomic_guess) density_guess = mopac_guess
260 safe_density_guess = mopac_guess
261 END IF
262 IF (dft_control%qs_control%xtb) THEN
263 IF (do_kpoints) THEN
264 IF (density_guess == atomic_guess) density_guess = mopac_guess
265 safe_density_guess = mopac_guess
266 ELSE
267 IF (density_guess == atomic_guess) density_guess = core_guess
268 safe_density_guess = core_guess
269 END IF
270 END IF
271
272 IF (scf_control%use_ot .AND. &
273 (.NOT. ((density_guess == random_guess) .OR. &
274 (density_guess == atomic_guess) .OR. &
275 (density_guess == core_guess) .OR. &
276 (density_guess == mopac_guess) .OR. &
277 (density_guess == eht_guess) .OR. &
278 (density_guess == sparse_guess) .OR. &
279 (((density_guess == restart_guess) .OR. &
280 (density_guess == history_guess)) .AND. &
281 (scf_control%level_shift == 0.0_dp))))) THEN
282 CALL cp_abort(__location__, &
283 "OT needs GUESS ATOMIC / CORE / RANDOM / SPARSE / RESTART / HISTORY RESTART: other options NYI")
284 END IF
285
286 ! if a restart was requested, check that the file exists,
287 ! if not we fall back to an atomic guess. No kidding, the file name should remain
288 ! in sync with read_mo_set_from_restart
289 id_nr = 0
290 IF (density_guess == restart_guess) THEN
291 ! only check existence on I/O node, otherwise if file exists there but
292 ! not on compute nodes, everything goes crazy even though only I/O
293 ! node actually reads the file
294 IF (do_kpoints) THEN
295 IF (para_env%is_source()) THEN
296 CALL wfn_restart_file_name(file_name, exist, dft_section, logger, kp=.true.)
297 END IF
298 ELSE
299 IF (para_env%is_source()) THEN
300 CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
301 END IF
302 END IF
303 CALL para_env%bcast(exist)
304 CALL para_env%bcast(file_name)
305 IF (.NOT. exist) THEN
306 CALL cp_warn(__location__, &
307 "User requested to restart the wavefunction from the file named: "// &
308 trim(file_name)//". This file does not exist. Please check the existence of"// &
309 " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME."// &
310 " Calculation continues using ATOMIC GUESS. ")
311 density_guess = safe_density_guess
312 END IF
313 ELSE IF (density_guess == history_guess) THEN
314 IF (do_kpoints) THEN
315 cpabort("calculate_first_density_matrix: history_guess not implemented for k-points")
316 END IF
317 IF (para_env%is_source()) THEN
318 CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
319 END IF
320 CALL para_env%bcast(exist)
321 CALL para_env%bcast(file_name)
322 nvec = qs_env%wf_history%memory_depth
323 not_read = nvec + 1
324 ! At this level we read the saved backup RESTART files..
325 DO i = 1, nvec
326 j = i - 1
327 filename = trim(file_name)
328 IF (j /= 0) THEN
329 filename = trim(file_name)//".bak-"//adjustl(cp_to_string(j))
330 END IF
331 IF (para_env%is_source()) &
332 INQUIRE (file=filename, exist=exist)
333 CALL para_env%bcast(exist)
334 IF ((.NOT. exist) .AND. (i < not_read)) THEN
335 not_read = i
336 END IF
337 END DO
338 IF (not_read == 1) THEN
339 density_guess = restart_guess
340 filename = trim(file_name)
341 IF (para_env%is_source()) INQUIRE (file=filename, exist=exist)
342 CALL para_env%bcast(exist)
343 IF (.NOT. exist) THEN
344 CALL cp_warn(__location__, &
345 "User requested to restart the wavefunction from a series of restart files named: "// &
346 trim(file_name)//" with extensions (.bak-n). These files do not exist."// &
347 " Even trying to switch to a plain restart wave-function failes because the"// &
348 " file named: "//trim(file_name)//" does not exist. Please check the existence of"// &
349 " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME."// &
350 " Calculation continues using ATOMIC GUESS. ")
351 density_guess = safe_density_guess
352 END IF
353 END IF
354 last_read = not_read - 1
355 END IF
356
357 did_guess = .false.
358
359 IF (dft_control%correct_el_density_dip) THEN
360 tot_corr_zeff = qs_env%total_zeff_corr
361 !WRITE(*,*) "tot_corr_zeff = ", tot_corr_zeff
362 IF ((abs(tot_corr_zeff) > 0.0_dp) .AND. (density_guess /= restart_guess)) THEN
363 CALL cp_warn(__location__, &
364 "Use SCF_GUESS RESTART in conjunction with "// &
365 "CORE_CORRECTION /= 0.0 and SURFACE_DIPOLE_CORRECTION TRUE. "// &
366 "It is always advisable to perform SURFACE_DIPOLE_CORRECTION "// &
367 "after a simulation without the surface dipole correction "// &
368 "and using the ensuing wavefunction restart file. ")
369 END IF
370 END IF
371
372 ounit = -1
373 print_log = .false.
374 print_history_log = .false.
375 IF (para_env%is_source()) THEN
376 CALL section_vals_val_get(dft_section, &
377 "SCF%PRINT%RESTART%LOG_PRINT_KEY", &
378 l_val=print_log)
379 CALL section_vals_val_get(dft_section, &
380 "SCF%PRINT%RESTART_HISTORY%LOG_PRINT_KEY", &
381 l_val=print_history_log)
382 IF (print_log .OR. print_history_log) THEN
383 ounit = cp_logger_get_default_io_unit(logger)
384 END IF
385 END IF
386
387 IF (density_guess == restart_guess) THEN
388 IF (ounit > 0) THEN
389 WRITE (unit=ounit, fmt="(/,T2,A)") &
390 "WFN_RESTART| Reading restart file"
391 END IF
392 IF (do_kpoints) THEN
393 natoms = SIZE(particle_set)
394 CALL read_kpoints_restart(rho_ao_kp, kpoints, work1, &
395 natoms, para_env, id_nr, dft_section, natom_mismatch)
396 IF (natom_mismatch) density_guess = safe_density_guess
397 ELSE
398 CALL read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, &
399 id_nr=id_nr, multiplicity=dft_control%multiplicity, dft_section=dft_section, &
400 natom_mismatch=natom_mismatch, out_unit=ounit)
401
402 IF (natom_mismatch) THEN
403 density_guess = safe_density_guess
404 ELSE
405 DO ispin = 1, nspin
406 IF (scf_control%level_shift /= 0.0_dp) THEN
407 CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
408 CALL cp_fm_to_fm(mo_coeff, ortho)
409 END IF
410
411 ! make all nmo vectors present orthonormal
412 CALL get_mo_set(mo_set=mo_array(ispin), &
413 mo_coeff=mo_coeff, nmo=nmo, homo=homo)
414
415 IF (has_unit_metric) THEN
416 CALL make_basis_simple(mo_coeff, nmo)
417 ELSEIF (dft_control%smear) THEN
418 CALL make_basis_lowdin(vmatrix=mo_coeff, ncol=nmo, &
419 matrix_s=s_sparse(1)%matrix)
420 ELSE
421 ! ortho so that one can restart for different positions (basis sets?)
422 CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
423 END IF
424 ! only alpha spin is kept for restricted
425 IF (dft_control%restricted) EXIT
426 END DO
427 IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
428
429 IF (.NOT. scf_control%diagonalization%mom) THEN
430 IF (dft_control%correct_surf_dip) THEN
431 IF (abs(tot_corr_zeff) > 0.0_dp) THEN
432 CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear, &
433 tot_zeff_corr=tot_corr_zeff)
434 ELSE
435 CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
436 END IF
437 ELSE
438 CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
439 END IF
440 END IF
441
442 DO ispin = 1, nspin
443
444 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
445 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
446 mo_array(ispin)%mo_coeff_b) !fm->dbcsr
447 END IF !fm->dbcsr
448
449 CALL calculate_density_matrix(mo_array(ispin), &
450 p_rmpv(ispin)%matrix)
451 END DO
452 END IF ! natom_mismatch
453
454 END IF
455
456 ! Maximum Overlap Method
457 IF (scf_control%diagonalization%mom) THEN
458 CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv)
459 END IF
460
461 did_guess = .true.
462 END IF
463
464 IF (density_guess == history_guess) THEN
465 IF (not_read > 1) THEN
466 IF (ounit > 0) THEN
467 WRITE (unit=ounit, fmt="(/,T2,A)") &
468 "WFN_RESTART| Reading restart file history"
469 END IF
470 DO i = 1, last_read
471 j = last_read - i
472 CALL read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, &
473 id_nr=j, multiplicity=dft_control%multiplicity, &
474 dft_section=dft_section, out_unit=ounit)
475
476 DO ispin = 1, nspin
477 IF (scf_control%level_shift /= 0.0_dp) THEN
478 CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
479 CALL cp_fm_to_fm(mo_coeff, ortho)
480 END IF
481
482 ! make all nmo vectors present orthonormal
483 CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, nmo=nmo, homo=homo)
484
485 IF (has_unit_metric) THEN
486 CALL make_basis_simple(mo_coeff, nmo)
487 ELSE
488 ! ortho so that one can restart for different positions (basis sets?)
489 CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
490 END IF
491 ! only alpha spin is kept for restricted
492 IF (dft_control%restricted) EXIT
493 END DO
494 IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
495
496 DO ispin = 1, nspin
497 CALL set_mo_occupation(mo_set=mo_array(ispin), &
498 smear=qs_env%scf_control%smear)
499 END DO
500
501 DO ispin = 1, nspin
502 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
503 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
504 mo_array(ispin)%mo_coeff_b) !fm->dbcsr
505 END IF !fm->dbcsr
506 CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
507 END DO
508
509 ! Write to extrapolation pipeline
510 CALL wfi_update(wf_history=qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
511 END DO
512 END IF
513
514 did_guess = .true.
515 END IF
516
517 IF (density_guess == random_guess) THEN
518
519 DO ispin = 1, nspin
520 CALL get_mo_set(mo_set=mo_array(ispin), &
521 mo_coeff=mo_coeff, nmo=nmo)
522 CALL cp_fm_init_random(mo_coeff, nmo)
523 IF (has_unit_metric) THEN
524 CALL make_basis_simple(mo_coeff, nmo)
525 ELSE
526 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
527 END IF
528 ! only alpha spin is kept for restricted
529 IF (dft_control%restricted) EXIT
530 END DO
531 IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
532
533 DO ispin = 1, nspin
534 CALL set_mo_occupation(mo_set=mo_array(ispin), &
535 smear=qs_env%scf_control%smear)
536 END DO
537
538 DO ispin = 1, nspin
539
540 IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
541 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
542 mo_array(ispin)%mo_coeff_b) !fm->dbcsr
543 END IF !fm->dbcsr
544
545 CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
546 END DO
547
548 did_guess = .true.
549 END IF
550
551 IF (density_guess == core_guess) THEN
552
553 IF (do_kpoints) THEN
554 cpabort("calculate_first_density_matrix: core_guess not implemented for k-points")
555 END IF
556
557 owns_ortho = .false.
558 IF (.NOT. ASSOCIATED(work1)) THEN
559 need_wm = .true.
560 cpassert(.NOT. ASSOCIATED(work2))
561 cpassert(.NOT. ASSOCIATED(ortho))
562 ELSE
563 need_wm = .false.
564 cpassert(ASSOCIATED(work2))
565 IF (.NOT. ASSOCIATED(ortho)) THEN
566 ALLOCATE (ortho)
567 owns_ortho = .true.
568 END IF
569 END IF
570
571 IF (need_wm) THEN
572 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
573 CALL cp_fm_get_info(moa, matrix_struct=ao_mo_struct)
574 CALL cp_fm_struct_get(ao_mo_struct, nrow_global=nao, nrow_block=nblocks)
575 CALL cp_fm_struct_create(fmstruct=ao_ao_struct, &
576 nrow_block=nblocks, &
577 ncol_block=nblocks, &
578 nrow_global=nao, &
579 ncol_global=nao, &
580 template_fmstruct=ao_mo_struct)
581 ALLOCATE (work1(1))
582 ALLOCATE (work2, ortho)
583 CALL cp_fm_create(work1(1), ao_ao_struct)
584 CALL cp_fm_create(work2, ao_ao_struct)
585 CALL cp_fm_create(ortho, ao_ao_struct)
586 CALL copy_dbcsr_to_fm(matrix_s_kp(1, 1)%matrix, ortho)
587 CALL cp_fm_cholesky_decompose(ortho)
588 CALL cp_fm_struct_release(ao_ao_struct)
589 END IF
590
591 ispin = 1
592 ! Load core Hamiltonian into work matrix
593 CALL copy_dbcsr_to_fm(h_core_sparse(1)%matrix, work1(ispin))
594
595 ! Diagonalize the core Hamiltonian matrix and retrieve a first set of
596 ! molecular orbitals (MOs)
597 IF (has_unit_metric) THEN
598 CALL eigensolver_simple(matrix_ks=work1(ispin), &
599 mo_set=mo_array(ispin), &
600 work=work2, &
601 do_level_shift=.false., &
602 level_shift=0.0_dp, &
603 use_jacobi=.false., jacobi_threshold=0._dp)
604 ELSE
605 CALL eigensolver(matrix_ks_fm=work1(ispin), &
606 mo_set=mo_array(ispin), &
607 ortho=ortho, &
608 work=work2, &
609 cholesky_method=scf_env%cholesky_method, &
610 do_level_shift=.false., &
611 level_shift=0.0_dp, &
612 use_jacobi=.false.)
613 END IF
614
615 ! Open shell case: copy alpha MOs to beta MOs
616 IF (nspin == 2) THEN
617 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
618 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mob, nmo=nmo)
619 CALL cp_fm_to_fm(moa, mob, nmo)
620 END IF
621
622 ! Build an initial density matrix (for each spin in the case of
623 ! an open shell calculation) from the first MOs set
624 DO ispin = 1, nspin
625 CALL set_mo_occupation(mo_set=mo_array(ispin), smear=scf_control%smear)
626 CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
627 END DO
628
629 ! release intermediate matrices
630 IF (need_wm) THEN
631 CALL cp_fm_release(ortho)
632 CALL cp_fm_release(work2)
633 CALL cp_fm_release(work1(1))
634 DEALLOCATE (ortho, work2)
635 DEALLOCATE (work1)
636 NULLIFY (work1, work2, ortho)
637 ELSE IF (owns_ortho) THEN
638 DEALLOCATE (ortho)
639 END IF
640
641 did_guess = .true.
642 END IF
643
644 IF (density_guess == atomic_guess) THEN
645
646 subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
647 ounit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
648 IF (ounit > 0) THEN
649 WRITE (unit=ounit, fmt="(/,(T2,A))") &
650 "Atomic guess: The first density matrix is obtained in terms of atomic orbitals", &
651 " and electronic configurations assigned to each atomic kind"
652 END IF
653
654 CALL calculate_atomic_block_dm(p_rmpv, s_sparse(1)%matrix, atomic_kind_set, qs_kind_set, &
655 nspin, nelectron_spin, ounit, para_env)
656
657 DO ispin = 1, nspin
658
659 ! The orbital transformation method (OT) requires not only an
660 ! initial density matrix, but also an initial wavefunction (MO set)
661 IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN
662 ! get orbitals later
663 ELSE
664 IF (need_mos) THEN
665
666 IF (dft_control%restricted .AND. (ispin == 2)) THEN
667 CALL mo_set_restrict(mo_array)
668 ELSE
669 CALL get_mo_set(mo_set=mo_array(ispin), &
670 mo_coeff=mo_coeff, &
671 nmo=nmo, nao=nao, homo=homo)
672
673 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
674 CALL cp_fm_init_random(mo_coeff, nmo)
675
676 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
677 ! multiply times PS
678 IF (has_unit_metric) THEN
679 CALL cp_fm_to_fm(mo_coeff, sv)
680 ELSE
681 ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
682 CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
683 END IF
684 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
685
686 CALL cp_fm_release(sv)
687 ! and ortho the result
688 IF (has_unit_metric) THEN
689 CALL make_basis_simple(mo_coeff, nmo)
690 ELSE
691 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
692 END IF
693 END IF
694
695 CALL set_mo_occupation(mo_set=mo_array(ispin), &
696 smear=qs_env%scf_control%smear)
697
698 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
699 mo_array(ispin)%mo_coeff_b) !fm->dbcsr
700
701 CALL calculate_density_matrix(mo_array(ispin), &
702 p_rmpv(ispin)%matrix)
703 END IF
704 ! adjust el_density in case surface_dipole_correction is switched
705 ! on and CORE_CORRECTION is non-zero
706 IF (scf_env%method == general_diag_method_nr) THEN
707 IF (dft_control%correct_surf_dip) THEN
708 IF (abs(tot_corr_zeff) > 0.0_dp) THEN
709 CALL get_mo_set(mo_set=mo_array(ispin), &
710 mo_coeff=mo_coeff, &
711 nmo=nmo, nao=nao, homo=homo)
712
713 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
714 CALL cp_fm_init_random(mo_coeff, nmo)
715
716 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
717 ! multiply times PS
718 IF (has_unit_metric) THEN
719 CALL cp_fm_to_fm(mo_coeff, sv)
720 ELSE
721 ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
722 CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
723 END IF
724 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
725
726 CALL cp_fm_release(sv)
727 ! and ortho the result
728 IF (has_unit_metric) THEN
729 CALL make_basis_simple(mo_coeff, nmo)
730 ELSE
731 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
732 END IF
733
734 CALL set_mo_occupation(mo_set=mo_array(ispin), smear=qs_env%scf_control%smear, &
735 tot_zeff_corr=tot_corr_zeff)
736
737 CALL calculate_density_matrix(mo_array(ispin), &
738 p_rmpv(ispin)%matrix)
739 END IF
740 END IF
741 END IF
742
743 END IF
744
745 END DO
746
747 IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN
748 ! We fit a function to the square root of the density
749 CALL qs_rho_update_rho(rho, qs_env)
750 cpassert(1 == 0)
751! CALL cp_fm_create(sv,mo_coeff%matrix_struct,"SV")
752! DO ispin=1,nspin
753! CALL integrate_ppl_rspace(qs%rho%rho_r(ispin),qs_env)
754! CALL cp_cfm_solve(overlap,mos)
755! CALL get_mo_set(mo_set=mo_array(ispin),&
756! mo_coeff=mo_coeff, nmo=nmo, nao=nao)
757! CALL cp_fm_init_random(mo_coeff,nmo)
758! END DO
759! CALL cp_fm_release(sv)
760 END IF
761
762 IF (scf_control%diagonalization%mom) THEN
763 CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv)
764 END IF
765
766 CALL cp_print_key_finished_output(ounit, logger, subsys_section, &
767 "PRINT%KINDS")
768
769 did_guess = .true.
770 END IF
771
772 IF (density_guess == sparse_guess) THEN
773
774 IF (ofgpw) cpabort("SPARSE_GUESS not implemented for OFGPW")
775 IF (.NOT. scf_control%use_ot) cpabort("OT needed!")
776 IF (do_kpoints) THEN
777 cpabort("calculate_first_density_matrix: sparse_guess not implemented for k-points")
778 END IF
779
780 eps = 1.0e-5_dp
781
782 ounit = cp_logger_get_default_io_unit(logger)
783 natoms = SIZE(particle_set)
784 ALLOCATE (kind_of(natoms))
785 ALLOCATE (first_sgf(natoms), last_sgf(natoms))
786
787 checksum = dbcsr_checksum(s_sparse(1)%matrix)
788 i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
789 IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum
790 CALL dbcsr_filter(s_sparse(1)%matrix, eps)
791 checksum = dbcsr_checksum(s_sparse(1)%matrix)
792 i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
793 IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum
794
795 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, &
796 last_sgf=last_sgf)
797 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
798
799 ALLOCATE (pmat(SIZE(atomic_kind_set)))
800
801 rscale = 1._dp
802 IF (nspin == 2) rscale = 0.5_dp
803 DO ikind = 1, SIZE(atomic_kind_set)
804 atomic_kind => atomic_kind_set(ikind)
805 qs_kind => qs_kind_set(ikind)
806 NULLIFY (pmat(ikind)%mat)
807 CALL calculate_atomic_orbitals(atomic_kind, qs_kind, pmat=pmat(ikind)%mat)
808 NULLIFY (atomic_kind)
809 END DO
810
811 DO ispin = 1, nspin
812 CALL get_mo_set(mo_set=mo_array(ispin), &
813 maxocc=maxocc, &
814 nelectron=nelectron)
815 !
816 CALL dbcsr_iterator_start(iter, p_rmpv(ispin)%matrix)
817 DO WHILE (dbcsr_iterator_blocks_left(iter))
818 CALL dbcsr_iterator_next_block(iter, irow, icol, pdata)
819 ikind = kind_of(irow)
820 IF (icol .EQ. irow) THEN
821 IF (ispin == 1) THEN
822 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
823 pmat(ikind)%mat(:, :, 2)*rscale
824 ELSE
825 pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
826 pmat(ikind)%mat(:, :, 2)*rscale
827 END IF
828 END IF
829 END DO
830 CALL dbcsr_iterator_stop(iter)
831
832 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
833 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
834 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
835 IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum
836 ! so far p needs to have the same sparsity as S
837 !CALL dbcsr_filter(p_rmpv(ispin)%matrix, eps)
838 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
839 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
840 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
841 IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum
842
843 CALL dbcsr_dot(p_rmpv(ispin)%matrix, s_sparse(1)%matrix, trps1)
844 rscale = real(nelectron, dp)/trps1
845 CALL dbcsr_scale(p_rmpv(ispin)%matrix, rscale)
846
847 !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
848 checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
849 occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
850 IF (ounit > 0) WRITE (ounit, *) 'P occ', occ, ' checksum', checksum
851 !
852 ! The orbital transformation method (OT) requires not only an
853 ! initial density matrix, but also an initial wavefunction (MO set)
854 IF (dft_control%restricted .AND. (ispin == 2)) THEN
855 CALL mo_set_restrict(mo_array)
856 ELSE
857 CALL get_mo_set(mo_set=mo_array(ispin), &
858 mo_coeff=mo_coeff, &
859 nmo=nmo, nao=nao, homo=homo)
860 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
861
862 n = maxval(last_sgf - first_sgf) + 1
863 size_atomic_kind_set = SIZE(atomic_kind_set)
864
865 ALLOCATE (buff(n, n), sort_kind(size_atomic_kind_set), &
866 nelec_kind(size_atomic_kind_set))
867 !
868 ! sort kind vs nbr electron
869 DO ikind = 1, size_atomic_kind_set
870 atomic_kind => atomic_kind_set(ikind)
871 qs_kind => qs_kind_set(ikind)
872 CALL get_atomic_kind(atomic_kind=atomic_kind, &
873 natom=natom, &
874 atom_list=atom_list, &
875 z=z)
876 CALL get_qs_kind(qs_kind, nsgf=nsgf, elec_conf=elec_conf, &
877 basis_set=orb_basis_set, zeff=zeff)
878 nelec_kind(ikind) = sum(elec_conf)
879 END DO
880 CALL sort(nelec_kind, size_atomic_kind_set, sort_kind)
881 !
882 ! a -very- naive sparse guess
883 nmo_tmp = nmo
884 natoms_tmp = natoms
885 istart_col = 1
886 iseed(1) = 4; iseed(2) = 3; iseed(3) = 2; iseed(4) = 1 ! set the seed for dlarnv
887 DO i = 1, size_atomic_kind_set
888 ikind = sort_kind(i)
889 atomic_kind => atomic_kind_set(ikind)
890 CALL get_atomic_kind(atomic_kind=atomic_kind, &
891 natom=natom, atom_list=atom_list)
892 DO iatom = 1, natom
893 !
894 atom_a = atom_list(iatom)
895 istart_row = first_sgf(atom_a)
896 n_rows = last_sgf(atom_a) - first_sgf(atom_a) + 1
897 !
898 ! compute the "potential" nbr of states for this atom
899 n_cols = max(int(real(nmo_tmp, dp)/real(natoms_tmp, dp)), 1)
900 IF (n_cols .GT. n_rows) n_cols = n_rows
901 !
902 nmo_tmp = nmo_tmp - n_cols
903 natoms_tmp = natoms_tmp - 1
904 IF (nmo_tmp .LT. 0 .OR. natoms_tmp .LT. 0) THEN
905 cpabort("Wrong1!")
906 END IF
907 DO j = 1, n_cols
908 CALL dlarnv(1, iseed, n_rows, buff(1, j))
909 END DO
910 CALL cp_fm_set_submatrix(mo_coeff, buff, istart_row, istart_col, &
911 n_rows, n_cols)
912 istart_col = istart_col + n_cols
913 END DO
914 END DO
915
916 IF (istart_col .LE. nmo) THEN
917 cpabort("Wrong2!")
918 END IF
919
920 DEALLOCATE (buff, nelec_kind, sort_kind)
921
922 IF (.false.) THEN
923 ALLOCATE (buff(nao, 1), buff2(nao, 1))
924 DO i = 1, nmo
925 CALL cp_fm_get_submatrix(mo_coeff, buff, 1, i, nao, 1)
926 IF (sum(buff**2) .LT. 1e-10_dp) THEN
927 WRITE (*, *) 'wrong', i, sum(buff**2)
928 END IF
929 length = sqrt(dot_product(buff(:, 1), buff(:, 1)))
930 buff(:, :) = buff(:, :)/length
931 DO j = i + 1, nmo
932 CALL cp_fm_get_submatrix(mo_coeff, buff2, 1, j, nao, 1)
933 length = sqrt(dot_product(buff2(:, 1), buff2(:, 1)))
934 buff2(:, :) = buff2(:, :)/length
935 IF (abs(dot_product(buff(:, 1), buff2(:, 1)) - 1.0_dp) .LT. 1e-10_dp) THEN
936 WRITE (*, *) 'wrong2', i, j, dot_product(buff(:, 1), buff2(:, 1))
937 DO ikind = 1, nao
938 IF (abs(mo_coeff%local_data(ikind, i)) .GT. 1e-10_dp) THEN
939 WRITE (*, *) 'c1', ikind, mo_coeff%local_data(ikind, i)
940 END IF
941 IF (abs(mo_coeff%local_data(ikind, j)) .GT. 1e-10_dp) THEN
942 WRITE (*, *) 'c2', ikind, mo_coeff%local_data(ikind, j)
943 END IF
944 END DO
945 cpabort("")
946 END IF
947 END DO
948 END DO
949 DEALLOCATE (buff, buff2)
950
951 END IF
952 !
953 CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, mo_coeff, s_sparse(1)%matrix)
954 !CALL dbcsr_verify_matrix(mo_dbcsr)
955 checksum = dbcsr_checksum(mo_dbcsr)
956
957 occ = dbcsr_get_occupation(mo_dbcsr)
958 IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
959 CALL dbcsr_filter(mo_dbcsr, eps)
960 !CALL dbcsr_verify_matrix(mo_dbcsr)
961 occ = dbcsr_get_occupation(mo_dbcsr)
962 checksum = dbcsr_checksum(mo_dbcsr)
963 IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
964 !
965 ! multiply times PS
966 IF (has_unit_metric) THEN
967 cpabort("has_unit_metric will be removed soon")
968 END IF
969 !
970 ! S*C
971 CALL dbcsr_copy(mo_tmp_dbcsr, mo_dbcsr, name="mo_tmp")
972 CALL dbcsr_multiply("N", "N", 1.0_dp, s_sparse(1)%matrix, mo_dbcsr, &
973 0.0_dp, mo_tmp_dbcsr, &
974 retain_sparsity=.true.)
975 !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
976 checksum = dbcsr_checksum(mo_tmp_dbcsr)
977 occ = dbcsr_get_occupation(mo_tmp_dbcsr)
978 IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum
979 CALL dbcsr_filter(mo_tmp_dbcsr, eps)
980 !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
981 checksum = dbcsr_checksum(mo_tmp_dbcsr)
982 occ = dbcsr_get_occupation(mo_tmp_dbcsr)
983 IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum
984 !
985 ! P*SC
986 ! the destroy is needed for the moment to avoid memory leaks !
987 ! This one is not needed because _destroy takes care of zeroing.
988 CALL dbcsr_multiply("N", "N", 1.0_dp, p_rmpv(ispin)%matrix, &
989 mo_tmp_dbcsr, 0.0_dp, mo_dbcsr)
990 IF (.false.) CALL dbcsr_verify_matrix(mo_dbcsr)
991 checksum = dbcsr_checksum(mo_dbcsr)
992 occ = dbcsr_get_occupation(mo_dbcsr)
993 IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
994 CALL dbcsr_filter(mo_dbcsr, eps)
995 !CALL dbcsr_verify_matrix(mo_dbcsr)
996 checksum = dbcsr_checksum(mo_dbcsr)
997 occ = dbcsr_get_occupation(mo_dbcsr)
998 IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
999 !
1000 CALL copy_dbcsr_to_fm(mo_dbcsr, mo_coeff)
1001
1002 CALL dbcsr_release(mo_dbcsr)
1003 CALL dbcsr_release(mo_tmp_dbcsr)
1004
1005 ! and ortho the result
1006 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
1007 END IF
1008
1009 CALL set_mo_occupation(mo_set=mo_array(ispin), &
1010 smear=qs_env%scf_control%smear)
1011
1012 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1013 mo_array(ispin)%mo_coeff_b) !fm->dbcsr
1014
1015 CALL calculate_density_matrix(mo_array(ispin), &
1016 p_rmpv(ispin)%matrix)
1017 DO ikind = 1, SIZE(atomic_kind_set)
1018 IF (ASSOCIATED(pmat(ikind)%mat)) THEN
1019 DEALLOCATE (pmat(ikind)%mat)
1020 END IF
1021 END DO
1022 END DO
1023
1024 DEALLOCATE (pmat)
1025
1026 DEALLOCATE (kind_of)
1027
1028 DEALLOCATE (first_sgf, last_sgf)
1029
1030 did_guess = .true.
1031 END IF
1032 IF (density_guess == mopac_guess) THEN
1033
1034 CALL calculate_mopac_dm(p_rmpv, s_sparse(1)%matrix, has_unit_metric, dft_control, &
1035 particle_set, atomic_kind_set, qs_kind_set, &
1036 nspin, nelectron_spin, para_env)
1037
1038 DO ispin = 1, nspin
1039 ! The orbital transformation method (OT) requires not only an
1040 ! initial density matrix, but also an initial wavefunction (MO set)
1041 IF (need_mos) THEN
1042 IF (dft_control%restricted .AND. (ispin == 2)) THEN
1043 CALL mo_set_restrict(mo_array)
1044 ELSE
1045 CALL get_mo_set(mo_set=mo_array(ispin), &
1046 mo_coeff=mo_coeff, &
1047 nmo=nmo, homo=homo)
1048 CALL cp_fm_init_random(mo_coeff, nmo)
1049 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
1050 ! multiply times PS
1051 IF (has_unit_metric) THEN
1052 CALL cp_fm_to_fm(mo_coeff, sv)
1053 ELSE
1054 CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
1055 END IF
1056 ! here we could easily multiply with the diag that we actually have replicated already
1057 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
1058 CALL cp_fm_release(sv)
1059 ! and ortho the result
1060 IF (has_unit_metric) THEN
1061 CALL make_basis_simple(mo_coeff, nmo)
1062 ELSE
1063 CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
1064 END IF
1065 END IF
1066
1067 CALL set_mo_occupation(mo_set=mo_array(ispin), &
1068 smear=qs_env%scf_control%smear)
1069 CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1070 mo_array(ispin)%mo_coeff_b)
1071
1072 CALL calculate_density_matrix(mo_array(ispin), &
1073 p_rmpv(ispin)%matrix)
1074 END IF
1075 END DO
1076
1077 did_guess = .true.
1078 END IF
1079 !
1080 ! EHT guess (gfn0-xTB)
1081 IF (density_guess == eht_guess) THEN
1082 CALL calculate_eht_guess(qs_env, mo_array)
1083 DO ispin = 1, nspin
1084 CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
1085 END DO
1086 did_guess = .true.
1087 END IF
1088 ! switch_surf_dip [SGh]
1089 IF (dft_control%switch_surf_dip) THEN
1090 DO ispin = 1, nspin
1091 CALL reassign_allocated_mos(mos_last_converged(ispin), &
1092 mo_array(ispin))
1093 END DO
1094 END IF
1095
1096 IF (density_guess == no_guess) THEN
1097 did_guess = .true.
1098 END IF
1099
1100 IF (.NOT. did_guess) THEN
1101 cpabort("An invalid keyword for the initial density guess was specified")
1102 END IF
1103
1104 CALL timestop(handle)
1105
1106 END SUBROUTINE calculate_first_density_matrix
1107
1108! **************************************************************************************************
1109!> \brief returns a block diagonal fock matrix.
1110!> \param matrix_f ...
1111!> \param atomic_kind_set ...
1112!> \param qs_kind_set ...
1113!> \param ounit ...
1114! **************************************************************************************************
1115 SUBROUTINE calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
1116 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_f
1117 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1118 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1119 INTEGER, INTENT(IN) :: ounit
1120
1121 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atomic_fock_matrix'
1122
1123 INTEGER :: handle, icol, ikind, irow
1124 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1125 REAL(dp), DIMENSION(:, :), POINTER :: block
1126 TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: fmat
1127 TYPE(atomic_kind_type), POINTER :: atomic_kind
1128 TYPE(dbcsr_iterator_type) :: iter
1129 TYPE(qs_kind_type), POINTER :: qs_kind
1130
1131 CALL timeset(routinen, handle)
1132
1133 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1134 ALLOCATE (fmat(SIZE(atomic_kind_set)))
1135
1136 ! precompute the atomic blocks for each atomic-kind
1137 DO ikind = 1, SIZE(atomic_kind_set)
1138 atomic_kind => atomic_kind_set(ikind)
1139 qs_kind => qs_kind_set(ikind)
1140 NULLIFY (fmat(ikind)%mat)
1141 IF (ounit > 0) WRITE (unit=ounit, fmt="(/,T2,A)") &
1142 "Calculating atomic Fock matrix for atomic kind: "//trim(atomic_kind%name)
1143
1144 !Currently only ispin=1 is supported
1145 CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
1146 fmat=fmat(ikind)%mat)
1147 END DO
1148
1149 ! zero result matrix
1150 CALL dbcsr_set(matrix_f, 0.0_dp)
1151
1152 ! copy precomputed blocks onto diagonal of result matrix
1153 CALL dbcsr_iterator_start(iter, matrix_f)
1154 DO WHILE (dbcsr_iterator_blocks_left(iter))
1155 CALL dbcsr_iterator_next_block(iter, irow, icol, block)
1156 ikind = kind_of(irow)
1157 IF (icol .EQ. irow) block(:, :) = fmat(ikind)%mat(:, :, 1)
1158 END DO
1159 CALL dbcsr_iterator_stop(iter)
1160
1161 ! cleanup
1162 DO ikind = 1, SIZE(atomic_kind_set)
1163 DEALLOCATE (fmat(ikind)%mat)
1164 END DO
1165 DEALLOCATE (fmat)
1166
1167 CALL timestop(handle)
1168
1169 END SUBROUTINE calculate_atomic_fock_matrix
1170
1171! **************************************************************************************************
1172!> \brief returns a block diagonal density matrix. Blocks correspond to the mopac initial guess.
1173!> \param pmat ...
1174!> \param matrix_s ...
1175!> \param has_unit_metric ...
1176!> \param dft_control ...
1177!> \param particle_set ...
1178!> \param atomic_kind_set ...
1179!> \param qs_kind_set ...
1180!> \param nspin ...
1181!> \param nelectron_spin ...
1182!> \param para_env ...
1183! **************************************************************************************************
1184 SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
1185 dft_control, particle_set, atomic_kind_set, qs_kind_set, &
1186 nspin, nelectron_spin, para_env)
1187 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmat
1188 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
1189 LOGICAL :: has_unit_metric
1190 TYPE(dft_control_type), POINTER :: dft_control
1191 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1192 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1193 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1194 INTEGER, INTENT(IN) :: nspin
1195 INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
1196 TYPE(mp_para_env_type) :: para_env
1197
1198 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_mopac_dm'
1199
1200 INTEGER :: atom_a, handle, iatom, ikind, iset, &
1201 isgf, isgfa, ishell, ispin, la, maxl, &
1202 maxll, na, nao, natom, ncount, nset, &
1203 nsgf, z
1204 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
1205 INTEGER, DIMENSION(25) :: laox, naox
1206 INTEGER, DIMENSION(5) :: occupation
1207 INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nshell
1208 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, l, last_sgfa
1209 LOGICAL :: has_pot
1210 REAL(kind=dp) :: maxocc, my_sum, nelec, occ, paa, rscale, &
1211 trps1, trps2, yy, zeff
1212 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: econf, pdiag, sdiag
1213 REAL(kind=dp), DIMENSION(0:3) :: edftb
1214 TYPE(all_potential_type), POINTER :: all_potential
1215 TYPE(dbcsr_type), POINTER :: matrix_p
1216 TYPE(gth_potential_type), POINTER :: gth_potential
1217 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1218 TYPE(sgp_potential_type), POINTER :: sgp_potential
1219 TYPE(xtb_atom_type), POINTER :: xtb_kind
1220
1221 CALL timeset(routinen, handle)
1222
1223 DO ispin = 1, nspin
1224 matrix_p => pmat(ispin)%matrix
1225 CALL dbcsr_set(matrix_p, 0.0_dp)
1226 END DO
1227
1228 natom = SIZE(particle_set)
1229 CALL dbcsr_get_info(pmat(1)%matrix, nfullrows_total=nao)
1230 IF (nspin == 1) THEN
1231 maxocc = 2.0_dp
1232 ELSE
1233 maxocc = 1.0_dp
1234 END IF
1235
1236 ALLOCATE (first_sgf(natom))
1237
1238 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1239 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxl)
1240
1241 ALLOCATE (econf(0:maxl))
1242
1243 ALLOCATE (pdiag(nao))
1244 pdiag(:) = 0.0_dp
1245
1246 ALLOCATE (sdiag(nao))
1247
1248 sdiag(:) = 0.0_dp
1249 IF (has_unit_metric) THEN
1250 sdiag(:) = 1.0_dp
1251 ELSE
1252 CALL dbcsr_get_diag(matrix_s, sdiag)
1253 CALL para_env%sum(sdiag)
1254 END IF
1255
1256 ncount = 0
1257 trps1 = 0.0_dp
1258 trps2 = 0.0_dp
1259 pdiag(:) = 0.0_dp
1260
1261 IF (sum(nelectron_spin) /= 0) THEN
1262 DO ikind = 1, SIZE(atomic_kind_set)
1263
1264 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
1265 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
1266 all_potential=all_potential, &
1267 gth_potential=gth_potential, &
1268 sgp_potential=sgp_potential)
1269 has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
1270
1271 IF (dft_control%qs_control%dftb) THEN
1272 CALL get_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter, &
1273 lmax=maxll, occupation=edftb)
1274 maxll = min(maxll, maxl)
1275 econf(0:maxl) = edftb(0:maxl)
1276 ELSEIF (dft_control%qs_control%xtb) THEN
1277 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1278 CALL get_xtb_atom_param(xtb_kind, z=z, natorb=nsgf, nao=naox, lao=laox, occupation=occupation)
1279 ELSEIF (has_pot) THEN
1280 CALL get_atomic_kind(atomic_kind_set(ikind), z=z)
1281 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, elec_conf=elec_conf, zeff=zeff)
1282 maxll = min(SIZE(elec_conf) - 1, maxl)
1283 econf(:) = 0.0_dp
1284 econf(0:maxll) = 0.5_dp*maxocc*real(elec_conf(0:maxll), dp)
1285 ELSE
1286 cycle
1287 END IF
1288
1289 ! MOPAC TYPE GUESS
1290 IF (dft_control%qs_control%dftb) THEN
1291 DO iatom = 1, natom
1292 atom_a = atom_list(iatom)
1293 isgfa = first_sgf(atom_a)
1294 DO la = 0, maxll
1295 SELECT CASE (la)
1296 CASE (0)
1297 pdiag(isgfa) = econf(0)
1298 CASE (1)
1299 pdiag(isgfa + 1) = econf(1)/3._dp
1300 pdiag(isgfa + 2) = econf(1)/3._dp
1301 pdiag(isgfa + 3) = econf(1)/3._dp
1302 CASE (2)
1303 pdiag(isgfa + 4) = econf(2)/5._dp
1304 pdiag(isgfa + 5) = econf(2)/5._dp
1305 pdiag(isgfa + 6) = econf(2)/5._dp
1306 pdiag(isgfa + 7) = econf(2)/5._dp
1307 pdiag(isgfa + 8) = econf(2)/5._dp
1308 CASE (3)
1309 pdiag(isgfa + 9) = econf(3)/7._dp
1310 pdiag(isgfa + 10) = econf(3)/7._dp
1311 pdiag(isgfa + 11) = econf(3)/7._dp
1312 pdiag(isgfa + 12) = econf(3)/7._dp
1313 pdiag(isgfa + 13) = econf(3)/7._dp
1314 pdiag(isgfa + 14) = econf(3)/7._dp
1315 pdiag(isgfa + 15) = econf(3)/7._dp
1316 CASE DEFAULT
1317 cpabort("")
1318 END SELECT
1319 END DO
1320 END DO
1321 ELSEIF (dft_control%qs_control%xtb) THEN
1322 DO iatom = 1, natom
1323 atom_a = atom_list(iatom)
1324 isgfa = first_sgf(atom_a)
1325 IF (z == 1 .AND. nsgf == 2) THEN
1326 ! Hydrogen 2s basis
1327 pdiag(isgfa) = 1.0_dp
1328 pdiag(isgfa + 1) = 0.0_dp
1329 ELSE
1330 DO isgf = 1, nsgf
1331 na = naox(isgf)
1332 la = laox(isgf)
1333 occ = real(occupation(la + 1), dp)/real(2*la + 1, dp)
1334 pdiag(isgfa + isgf - 1) = occ
1335 END DO
1336 END IF
1337 END DO
1338 ELSEIF (dft_control%qs_control%semi_empirical) THEN
1339 yy = real(dft_control%charge, kind=dp)/real(nao, kind=dp)
1340 DO iatom = 1, natom
1341 atom_a = atom_list(iatom)
1342 isgfa = first_sgf(atom_a)
1343 SELECT CASE (nsgf)
1344 CASE (1) ! s-basis
1345 pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1346 CASE (4) ! sp-basis
1347 IF (z == 1) THEN
1348 ! special case: hydrogen with sp basis
1349 pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1350 pdiag(isgfa + 1) = 0._dp
1351 pdiag(isgfa + 2) = 0._dp
1352 pdiag(isgfa + 3) = 0._dp
1353 ELSE
1354 pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1355 pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1356 pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1357 pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1358 END IF
1359 CASE (9) ! spd-basis
1360 IF (z < 21 .OR. z > 30 .AND. z < 39 .OR. z > 48 .AND. z < 57) THEN
1361 ! Main Group Element: The "d" shell is formally empty.
1362 pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1363 pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1364 pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1365 pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1366 pdiag(isgfa + 4) = (-yy)*0.5_dp*maxocc
1367 pdiag(isgfa + 5) = (-yy)*0.5_dp*maxocc
1368 pdiag(isgfa + 6) = (-yy)*0.5_dp*maxocc
1369 pdiag(isgfa + 7) = (-yy)*0.5_dp*maxocc
1370 pdiag(isgfa + 8) = (-yy)*0.5_dp*maxocc
1371 ELSE IF (z < 99) THEN
1372 my_sum = zeff - 9.0_dp*yy
1373 ! First, put 2 electrons in the 's' shell
1374 pdiag(isgfa) = (max(0.0_dp, min(my_sum, 2.0_dp)))*0.5_dp*maxocc
1375 my_sum = my_sum - 2.0_dp
1376 IF (my_sum > 0.0_dp) THEN
1377 ! Now put as many electrons as possible into the 'd' shell
1378 pdiag(isgfa + 4) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1379 pdiag(isgfa + 5) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1380 pdiag(isgfa + 6) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1381 pdiag(isgfa + 7) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1382 pdiag(isgfa + 8) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1383 my_sum = max(0.0_dp, my_sum - 10.0_dp)
1384 ! Put the remaining electrons in the 'p' shell
1385 pdiag(isgfa + 1) = (my_sum/3.0_dp)*0.5_dp*maxocc
1386 pdiag(isgfa + 2) = (my_sum/3.0_dp)*0.5_dp*maxocc
1387 pdiag(isgfa + 3) = (my_sum/3.0_dp)*0.5_dp*maxocc
1388 END IF
1389 END IF
1390 CASE DEFAULT
1391 cpabort("")
1392 END SELECT
1393 END DO
1394 ELSE
1395 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1396 nset=nset, &
1397 nshell=nshell, &
1398 l=l, &
1399 first_sgf=first_sgfa, &
1400 last_sgf=last_sgfa)
1401
1402 DO iset = 1, nset
1403 DO ishell = 1, nshell(iset)
1404 la = l(ishell, iset)
1405 nelec = maxocc*real(2*la + 1, dp)
1406 IF (econf(la) > 0.0_dp) THEN
1407 IF (econf(la) >= nelec) THEN
1408 paa = maxocc
1409 econf(la) = econf(la) - nelec
1410 ELSE
1411 paa = maxocc*econf(la)/nelec
1412 econf(la) = 0.0_dp
1413 ncount = ncount + nint(nelec/maxocc)
1414 END IF
1415 DO isgfa = first_sgfa(ishell, iset), last_sgfa(ishell, iset)
1416 DO iatom = 1, natom
1417 atom_a = atom_list(iatom)
1418 isgf = first_sgf(atom_a) + isgfa - 1
1419 pdiag(isgf) = paa
1420 IF (paa == maxocc) THEN
1421 trps1 = trps1 + paa*sdiag(isgf)
1422 ELSE
1423 trps2 = trps2 + paa*sdiag(isgf)
1424 END IF
1425 END DO
1426 END DO
1427 END IF
1428 END DO ! ishell
1429 END DO ! iset
1430 END IF
1431 END DO ! ikind
1432
1433 IF (trps2 == 0.0_dp) THEN
1434 DO isgf = 1, nao
1435 IF (sdiag(isgf) > 0.0_dp) pdiag(isgf) = pdiag(isgf)/sdiag(isgf)
1436 END DO
1437 DO ispin = 1, nspin
1438 IF (nelectron_spin(ispin) /= 0) THEN
1439 matrix_p => pmat(ispin)%matrix
1440 CALL dbcsr_set_diag(matrix_p, pdiag)
1441 END IF
1442 END DO
1443 ELSE
1444 DO ispin = 1, nspin
1445 IF (nelectron_spin(ispin) /= 0) THEN
1446 rscale = (real(nelectron_spin(ispin), dp) - trps1)/trps2
1447 DO isgf = 1, nao
1448 IF (pdiag(isgf) < maxocc) pdiag(isgf) = rscale*pdiag(isgf)
1449 END DO
1450 matrix_p => pmat(ispin)%matrix
1451 CALL dbcsr_set_diag(matrix_p, pdiag)
1452 DO isgf = 1, nao
1453 IF (pdiag(isgf) < maxocc) pdiag(isgf) = pdiag(isgf)/rscale
1454 END DO
1455 END IF
1456 END DO
1457 END IF
1458 END IF
1459
1460 DEALLOCATE (econf)
1461
1462 DEALLOCATE (first_sgf)
1463
1464 DEALLOCATE (pdiag)
1465
1466 DEALLOCATE (sdiag)
1467
1468 CALL timestop(handle)
1469
1470 END SUBROUTINE calculate_mopac_dm
1471
1472END MODULE qs_initial_guess
calculate the orbitals for a given atomic kind type
subroutine, public calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, density, wavefunction, wfninfo, confine, xc_section, nocc)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_verify_matrix(matrix, verbosity, local)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_filter(matrix, eps)
...
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
integer function, public dbcsr_get_num_blocks(matrix)
...
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
real(kind=dp) function, public dbcsr_checksum(matrix, pos)
Calculates the checksum of a DBCSR matrix.
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition of the atomic potential types.
Types and set/get functions for HFX.
Definition hfx_types.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public core_guess
integer, parameter, public mopac_guess
integer, parameter, public no_guess
integer, parameter, public atomic_guess
integer, parameter, public history_guess
integer, parameter, public random_guess
integer, parameter, public sparse_guess
integer, parameter, public restart_guess
integer, parameter, public eht_guess
function that builds the hartree fock exchange section of the input
integer, parameter, public ri_mo
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_path_length
Definition kinds.F:58
Restart file for k point calculations.
Definition kpoint_io.F:13
subroutine, public read_kpoints_restart(denmat, kpoints, fmwork, natom, para_env, id_nr, dft_section, natom_mismatch)
...
Definition kpoint_io.F:297
Types and basic routines needed for a kpoint calculation.
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
subroutine, public calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, ounit, para_env)
returns a block diagonal density matrix. Blocks correspond to the atomic densities.
collects routines that calculate density matrices
Working with the DFTB parameter types.
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Generate an initial guess (dm and orb) from EHT calculation.
subroutine, public calculate_eht_guess(qs_env, mo_array)
EHT MO guess calclulation.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Routines to somehow generate an initial guess.
subroutine, public calculate_first_density_matrix(scf_env, qs_env)
can use a variety of methods to come up with an initial density matrix and optionally an initial wave...
subroutine, public calculate_mopac_dm(pmat, matrix_s, has_unit_metric, dft_control, particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env)
returns a block diagonal density matrix. Blocks correspond to the mopac initial guess.
subroutine, public calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
returns a block diagonal fock matrix.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
subroutine, public wfn_restart_file_name(filename, exist, section, logger, kp, xas, rtp)
...
Definition qs_mo_io.F:428
subroutine, public read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, id_nr, multiplicity, dft_section, natom_mismatch, cdft, out_unit)
...
Definition qs_mo_io.F:495
collects routines that perform operations directly related to MOs
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public mo_set_restrict(mo_array, convert_dbcsr)
make the beta orbitals explicitly equal to the alpha orbitals effectively copying the orbital data
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public reassign_allocated_mos(mo_set_new, mo_set_old)
reassign an already allocated mo_set
methods for deltaSCF calculations
subroutine, public do_mom_guess(nspins, mos, scf_control, p_rmpv)
initial guess for the maximum overlap method
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
module that contains the definitions of the scf types
integer, parameter, public ot_diag_method_nr
integer, parameter, public block_davidson_diag_method_nr
integer, parameter, public block_krylov_diag_method_nr
integer, parameter, public general_diag_method_nr
Storage of past states of the qs_env. Methods to interpolate (or actually normally extrapolate) the n...
subroutine, public wfi_update(wf_history, qs_env, dt)
updates the snapshot buffer, taking a new snapshot
parameters that control an scf iteration
All kind of helpful little routines.
Definition util.F:14
Definition of the xTB parameter types.
Definition xtb_types.F:20
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
Provides all information about an atomic kind.
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:509
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.