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