Line data Source code
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 : ! **************************************************************************************************
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_copy, dbcsr_filter, dbcsr_get_info, dbcsr_get_num_blocks, dbcsr_get_occupation, &
23 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
24 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
25 : dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_verify_matrix
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum,&
27 : dbcsr_dot,&
28 : dbcsr_get_diag,&
29 : dbcsr_set_diag
30 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
31 : copy_fm_to_dbcsr,&
32 : cp_dbcsr_sm_fm_multiply,&
33 : cp_fm_to_dbcsr_row_template
34 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_get,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: &
40 : cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, cp_fm_release, &
41 : cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
42 : USE cp_log_handling, ONLY: cp_get_default_logger,&
43 : cp_logger_get_default_io_unit,&
44 : cp_logger_type,&
45 : cp_to_string
46 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
47 : cp_print_key_unit_nr
48 : USE external_potential_types, ONLY: all_potential_type,&
49 : gth_potential_type,&
50 : sgp_potential_type
51 : USE hfx_types, ONLY: hfx_type
52 : USE input_constants, ONLY: &
53 : atomic_guess, core_guess, eht_guess, history_guess, mopac_guess, no_guess, random_guess, &
54 : restart_guess, sparse_guess
55 : USE input_cp2k_hfx, ONLY: ri_mo
56 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
57 : section_vals_type,&
58 : section_vals_val_get
59 : USE kinds, ONLY: default_path_length,&
60 : dp
61 : USE kpoint_io, ONLY: read_kpoints_restart
62 : USE kpoint_types, ONLY: kpoint_type
63 : USE message_passing, ONLY: mp_para_env_type
64 : USE particle_methods, ONLY: get_particle_set
65 : USE particle_types, ONLY: particle_type
66 : USE qs_atomic_block, ONLY: calculate_atomic_block_dm
67 : USE qs_density_matrices, ONLY: calculate_density_matrix
68 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
69 : USE qs_eht_guess, ONLY: calculate_eht_guess
70 : USE qs_environment_types, ONLY: get_qs_env,&
71 : qs_environment_type
72 : USE qs_kind_types, ONLY: get_qs_kind,&
73 : get_qs_kind_set,&
74 : qs_kind_type
75 : USE qs_mo_io, ONLY: read_mo_set_from_restart,&
76 : wfn_restart_file_name
77 : USE qs_mo_methods, ONLY: make_basis_lowdin,&
78 : make_basis_simple,&
79 : make_basis_sm
80 : USE qs_mo_occupation, ONLY: set_mo_occupation
81 : USE qs_mo_types, ONLY: get_mo_set,&
82 : mo_set_restrict,&
83 : mo_set_type,&
84 : reassign_allocated_mos
85 : USE qs_mom_methods, ONLY: do_mom_guess
86 : USE qs_rho_methods, ONLY: qs_rho_update_rho
87 : USE qs_rho_types, ONLY: qs_rho_get,&
88 : qs_rho_type
89 : USE qs_scf_methods, ONLY: eigensolver,&
90 : eigensolver_simple
91 : USE qs_scf_types, ONLY: block_davidson_diag_method_nr,&
92 : block_krylov_diag_method_nr,&
93 : general_diag_method_nr,&
94 : ot_diag_method_nr,&
95 : qs_scf_env_type
96 : USE qs_wf_history_methods, ONLY: wfi_update
97 : USE scf_control_types, ONLY: scf_control_type
98 : USE util, ONLY: sort
99 : USE xtb_types, ONLY: get_xtb_atom_param,&
100 : xtb_atom_type
101 : #include "./base/base_uses.f90"
102 :
103 : IMPLICIT NONE
104 :
105 : PRIVATE
106 :
107 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_initial_guess'
108 :
109 : PUBLIC :: calculate_first_density_matrix, calculate_mopac_dm
110 : PUBLIC :: calculate_atomic_fock_matrix
111 :
112 : TYPE atom_matrix_type
113 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat => NULL()
114 : END TYPE atom_matrix_type
115 :
116 : CONTAINS
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 6751 : 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 6751 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
145 : INTEGER, DIMENSION(2) :: nelectron_spin
146 6751 : INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nelec_kind, &
147 6751 : 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 6751 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: buff, buff2
151 6751 : 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 6751 : TYPE(atom_matrix_type), DIMENSION(:), POINTER :: pmat
156 6751 : 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 6751 : 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 6751 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h_core_sparse, matrix_ks, p_rmpv, &
165 6751 : s_sparse
166 6751 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h_kp, matrix_ks_kp, matrix_s_kp, &
167 6751 : 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 6751 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
172 : TYPE(kpoint_type), POINTER :: kpoints
173 6751 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array, mos_last_converged
174 : TYPE(mp_para_env_type), POINTER :: para_env
175 6751 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 6751 : 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 13502 : logger => cp_get_default_logger()
183 6751 : NULLIFY (atomic_kind, qs_kind, mo_coeff, orb_basis_set, atomic_kind_set, &
184 6751 : qs_kind_set, particle_set, ortho, work2, work1, mo_array, s_sparse, &
185 6751 : scf_control, dft_control, p_rmpv, para_env, h_core_sparse, matrix_ks, rho, &
186 6751 : mos_last_converged)
187 6751 : NULLIFY (dft_section, input, subsys_section)
188 6751 : NULLIFY (matrix_s_kp, matrix_h_kp, matrix_ks_kp, rho_ao_kp)
189 6751 : NULLIFY (moa, mob)
190 6751 : NULLIFY (atom_list, elec_conf, kpoints)
191 : edftb = 0.0_dp
192 6751 : tot_corr_zeff = 0.0_dp
193 :
194 6751 : 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 6751 : x_data=x_data)
214 :
215 6751 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
216 :
217 6751 : IF (dft_control%switch_surf_dip) THEN
218 2 : 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 15035 : DO ispin = 1, dft_control%nspins
223 90527 : DO ic = 1, SIZE(rho_ao_kp, 2)
224 83776 : CALL dbcsr_set(rho_ao_kp(ispin, ic)%matrix, 0.0_dp)
225 : END DO
226 : END DO
227 6751 : s_sparse => matrix_s_kp(:, 1)
228 6751 : h_core_sparse => matrix_h_kp(:, 1)
229 6751 : matrix_ks => matrix_ks_kp(:, 1)
230 6751 : p_rmpv => rho_ao_kp(:, 1)
231 :
232 6751 : work1 => scf_env%scf_work1
233 6751 : work2 => scf_env%scf_work2
234 6751 : ortho => scf_env%ortho
235 :
236 6751 : dft_section => section_vals_get_subs_vals(input, "DFT")
237 :
238 6751 : nspin = dft_control%nspins
239 6751 : ofgpw = dft_control%qs_control%ofgpw
240 6751 : density_guess = scf_control%density_guess
241 6751 : do_std_diag = .FALSE.
242 :
243 6751 : do_hfx_ri_mo = .FALSE.
244 6751 : IF (ASSOCIATED(x_data)) THEN
245 1160 : IF (x_data(1, 1)%do_hfx_ri) THEN
246 114 : IF (x_data(1, 1)%ri_data%flavor == ri_mo) do_hfx_ri_mo = .TRUE.
247 : END IF
248 : END IF
249 :
250 6751 : 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 6751 : .OR. do_hfx_ri_mo
256 :
257 6751 : safe_density_guess = atomic_guess
258 6751 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN
259 802 : IF (density_guess == atomic_guess) density_guess = mopac_guess
260 : safe_density_guess = mopac_guess
261 : END IF
262 6751 : IF (dft_control%qs_control%xtb) THEN
263 1092 : IF (do_kpoints) THEN
264 220 : IF (density_guess == atomic_guess) density_guess = mopac_guess
265 : safe_density_guess = mopac_guess
266 : ELSE
267 872 : IF (density_guess == atomic_guess) density_guess = core_guess
268 : safe_density_guess = core_guess
269 : END IF
270 : END IF
271 :
272 6751 : 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 0 : "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 6751 : id_nr = 0
290 6751 : 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 568 : IF (do_kpoints) THEN
295 8 : IF (para_env%is_source()) THEN
296 4 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger, kp=.TRUE.)
297 : END IF
298 : ELSE
299 560 : IF (para_env%is_source()) THEN
300 294 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
301 : END IF
302 : END IF
303 568 : CALL para_env%bcast(exist)
304 568 : CALL para_env%bcast(file_name)
305 568 : 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 90 : " Calculation continues using ATOMIC GUESS. ")
311 90 : density_guess = safe_density_guess
312 : END IF
313 6183 : ELSE IF (density_guess == history_guess) THEN
314 2 : IF (do_kpoints) THEN
315 0 : CPABORT("calculate_first_density_matrix: history_guess not implemented for k-points")
316 : END IF
317 2 : IF (para_env%is_source()) THEN
318 1 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
319 : END IF
320 2 : CALL para_env%bcast(exist)
321 2 : CALL para_env%bcast(file_name)
322 2 : nvec = qs_env%wf_history%memory_depth
323 2 : not_read = nvec + 1
324 : ! At this level we read the saved backup RESTART files..
325 6 : DO i = 1, nvec
326 4 : j = i - 1
327 4 : filename = TRIM(file_name)
328 4 : IF (j /= 0) THEN
329 2 : filename = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(j))
330 : END IF
331 4 : IF (para_env%is_source()) &
332 2 : INQUIRE (FILE=filename, exist=exist)
333 4 : CALL para_env%bcast(exist)
334 6 : IF ((.NOT. exist) .AND. (i < not_read)) THEN
335 : not_read = i
336 : END IF
337 : END DO
338 2 : IF (not_read == 1) THEN
339 0 : density_guess = restart_guess
340 0 : filename = TRIM(file_name)
341 0 : IF (para_env%is_source()) INQUIRE (FILE=filename, exist=exist)
342 0 : CALL para_env%bcast(exist)
343 0 : 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 0 : " Calculation continues using ATOMIC GUESS. ")
351 0 : density_guess = safe_density_guess
352 : END IF
353 : END IF
354 2 : last_read = not_read - 1
355 : END IF
356 :
357 6751 : did_guess = .FALSE.
358 :
359 6751 : IF (dft_control%correct_el_density_dip) THEN
360 4 : tot_corr_zeff = qs_env%total_zeff_corr
361 : !WRITE(*,*) "tot_corr_zeff = ", tot_corr_zeff
362 4 : 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 4 : "and using the ensuing wavefunction restart file. ")
369 : END IF
370 : END IF
371 :
372 6751 : ounit = -1
373 6751 : print_log = .FALSE.
374 6751 : print_history_log = .FALSE.
375 6751 : IF (para_env%is_source()) THEN
376 : CALL section_vals_val_get(dft_section, &
377 : "SCF%PRINT%RESTART%LOG_PRINT_KEY", &
378 3413 : l_val=print_log)
379 : CALL section_vals_val_get(dft_section, &
380 : "SCF%PRINT%RESTART_HISTORY%LOG_PRINT_KEY", &
381 3413 : l_val=print_history_log)
382 3413 : IF (print_log .OR. print_history_log) THEN
383 13 : ounit = cp_logger_get_default_io_unit(logger)
384 : END IF
385 : END IF
386 :
387 6751 : IF (density_guess == restart_guess) THEN
388 478 : IF (ounit > 0) THEN
389 : WRITE (UNIT=ounit, FMT="(/,T2,A)") &
390 4 : "WFN_RESTART| Reading restart file"
391 : END IF
392 478 : IF (do_kpoints) THEN
393 6 : natoms = SIZE(particle_set)
394 : CALL read_kpoints_restart(rho_ao_kp, kpoints, work1, &
395 6 : natoms, para_env, id_nr, dft_section, natom_mismatch)
396 6 : 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 472 : natom_mismatch=natom_mismatch, out_unit=ounit)
401 :
402 472 : IF (natom_mismatch) THEN
403 : density_guess = safe_density_guess
404 : ELSE
405 1256 : DO ispin = 1, nspin
406 804 : IF (scf_control%level_shift /= 0.0_dp) THEN
407 0 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
408 0 : 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 804 : mo_coeff=mo_coeff, nmo=nmo, homo=homo)
414 :
415 804 : IF (has_unit_metric) THEN
416 4 : CALL make_basis_simple(mo_coeff, nmo)
417 800 : ELSEIF (dft_control%smear) THEN
418 : CALL make_basis_lowdin(vmatrix=mo_coeff, ncol=nmo, &
419 104 : matrix_s=s_sparse(1)%matrix)
420 : ELSE
421 : ! ortho so that one can restart for different positions (basis sets?)
422 696 : CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
423 : END IF
424 : ! only alpha spin is kept for restricted
425 2060 : IF (dft_control%restricted) EXIT
426 : END DO
427 472 : IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
428 :
429 472 : IF (.NOT. scf_control%diagonalization%mom) THEN
430 456 : IF (dft_control%correct_surf_dip) THEN
431 0 : IF (ABS(tot_corr_zeff) > 0.0_dp) THEN
432 : CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear, &
433 0 : tot_zeff_corr=tot_corr_zeff)
434 : ELSE
435 0 : CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
436 : END IF
437 : ELSE
438 456 : CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
439 : END IF
440 : END IF
441 :
442 1296 : DO ispin = 1, nspin
443 :
444 824 : 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 564 : mo_array(ispin)%mo_coeff_b) !fm->dbcsr
447 : END IF !fm->dbcsr
448 :
449 : CALL calculate_density_matrix(mo_array(ispin), &
450 1296 : p_rmpv(ispin)%matrix)
451 : END DO
452 : END IF ! natom_mismatch
453 :
454 : END IF
455 :
456 : ! Maximum Overlap Method
457 478 : IF (scf_control%diagonalization%mom) THEN
458 16 : CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv)
459 : END IF
460 :
461 : did_guess = .TRUE.
462 : END IF
463 :
464 6751 : IF (density_guess == history_guess) THEN
465 2 : IF (not_read > 1) THEN
466 2 : IF (ounit > 0) THEN
467 : WRITE (UNIT=ounit, FMT="(/,T2,A)") &
468 1 : "WFN_RESTART| Reading restart file history"
469 : END IF
470 6 : DO i = 1, last_read
471 4 : 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 4 : dft_section=dft_section, out_unit=ounit)
475 :
476 8 : DO ispin = 1, nspin
477 4 : IF (scf_control%level_shift /= 0.0_dp) THEN
478 0 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
479 0 : CALL cp_fm_to_fm(mo_coeff, ortho)
480 : END IF
481 :
482 : ! make all nmo vectors present orthonormal
483 4 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, nmo=nmo, homo=homo)
484 :
485 4 : IF (has_unit_metric) THEN
486 0 : CALL make_basis_simple(mo_coeff, nmo)
487 : ELSE
488 : ! ortho so that one can restart for different positions (basis sets?)
489 4 : CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
490 : END IF
491 : ! only alpha spin is kept for restricted
492 12 : IF (dft_control%restricted) EXIT
493 : END DO
494 4 : IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
495 :
496 8 : DO ispin = 1, nspin
497 : CALL set_mo_occupation(mo_set=mo_array(ispin), &
498 8 : smear=qs_env%scf_control%smear)
499 : END DO
500 :
501 8 : DO ispin = 1, nspin
502 4 : 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 4 : mo_array(ispin)%mo_coeff_b) !fm->dbcsr
505 : END IF !fm->dbcsr
506 8 : CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
507 : END DO
508 :
509 : ! Write to extrapolation pipeline
510 6 : 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 6751 : IF (density_guess == random_guess) THEN
518 :
519 52 : DO ispin = 1, nspin
520 : CALL get_mo_set(mo_set=mo_array(ispin), &
521 30 : mo_coeff=mo_coeff, nmo=nmo)
522 30 : CALL cp_fm_init_random(mo_coeff, nmo)
523 30 : IF (has_unit_metric) THEN
524 2 : CALL make_basis_simple(mo_coeff, nmo)
525 : ELSE
526 28 : CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
527 : END IF
528 : ! only alpha spin is kept for restricted
529 82 : IF (dft_control%restricted) EXIT
530 : END DO
531 22 : IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
532 :
533 52 : DO ispin = 1, nspin
534 : CALL set_mo_occupation(mo_set=mo_array(ispin), &
535 52 : smear=qs_env%scf_control%smear)
536 : END DO
537 :
538 52 : DO ispin = 1, nspin
539 :
540 30 : 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 22 : mo_array(ispin)%mo_coeff_b) !fm->dbcsr
543 : END IF !fm->dbcsr
544 :
545 52 : CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
546 : END DO
547 :
548 : did_guess = .TRUE.
549 : END IF
550 :
551 6751 : IF (density_guess == core_guess) THEN
552 :
553 158 : IF (do_kpoints) THEN
554 0 : CPABORT("calculate_first_density_matrix: core_guess not implemented for k-points")
555 : END IF
556 :
557 158 : owns_ortho = .FALSE.
558 158 : IF (.NOT. ASSOCIATED(work1)) THEN
559 42 : need_wm = .TRUE.
560 42 : CPASSERT(.NOT. ASSOCIATED(work2))
561 42 : CPASSERT(.NOT. ASSOCIATED(ortho))
562 : ELSE
563 116 : need_wm = .FALSE.
564 116 : CPASSERT(ASSOCIATED(work2))
565 116 : IF (.NOT. ASSOCIATED(ortho)) THEN
566 6 : ALLOCATE (ortho)
567 6 : owns_ortho = .TRUE.
568 : END IF
569 : END IF
570 :
571 158 : IF (need_wm) THEN
572 42 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
573 42 : CALL cp_fm_get_info(moa, matrix_struct=ao_mo_struct)
574 42 : 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 42 : template_fmstruct=ao_mo_struct)
581 84 : ALLOCATE (work1(1))
582 42 : ALLOCATE (work2, ortho)
583 42 : CALL cp_fm_create(work1(1), ao_ao_struct)
584 42 : CALL cp_fm_create(work2, ao_ao_struct)
585 42 : CALL cp_fm_create(ortho, ao_ao_struct)
586 42 : CALL copy_dbcsr_to_fm(matrix_s_kp(1, 1)%matrix, ortho)
587 42 : CALL cp_fm_cholesky_decompose(ortho)
588 126 : CALL cp_fm_struct_release(ao_ao_struct)
589 : END IF
590 :
591 158 : ispin = 1
592 : ! Load core Hamiltonian into work matrix
593 158 : 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 158 : 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 6 : 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 152 : use_jacobi=.FALSE.)
613 : END IF
614 :
615 : ! Open shell case: copy alpha MOs to beta MOs
616 158 : IF (nspin == 2) THEN
617 22 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
618 22 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mob, nmo=nmo)
619 22 : 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 338 : DO ispin = 1, nspin
625 180 : CALL set_mo_occupation(mo_set=mo_array(ispin), smear=scf_control%smear)
626 338 : CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
627 : END DO
628 :
629 : ! release intermediate matrices
630 158 : IF (need_wm) THEN
631 42 : CALL cp_fm_release(ortho)
632 42 : CALL cp_fm_release(work2)
633 42 : CALL cp_fm_release(work1(1))
634 42 : DEALLOCATE (ortho, work2)
635 42 : DEALLOCATE (work1)
636 42 : NULLIFY (work1, work2, ortho)
637 116 : ELSE IF (owns_ortho) THEN
638 6 : DEALLOCATE (ortho)
639 : END IF
640 :
641 : did_guess = .TRUE.
642 : END IF
643 :
644 6751 : IF (density_guess == atomic_guess) THEN
645 :
646 4309 : subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
647 4309 : ounit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
648 4309 : IF (ounit > 0) THEN
649 : WRITE (UNIT=ounit, FMT="(/,(T2,A))") &
650 985 : "Atomic guess: The first density matrix is obtained in terms of atomic orbitals", &
651 1970 : " 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 4309 : nspin, nelectron_spin, ounit, para_env)
656 :
657 9563 : 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 9563 : IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN
662 : ! get orbitals later
663 : ELSE
664 5254 : IF (need_mos) THEN
665 :
666 2074 : IF (dft_control%restricted .AND. (ispin == 2)) THEN
667 22 : CALL mo_set_restrict(mo_array)
668 : ELSE
669 : CALL get_mo_set(mo_set=mo_array(ispin), &
670 : mo_coeff=mo_coeff, &
671 2052 : nmo=nmo, nao=nao, homo=homo)
672 :
673 2052 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
674 2052 : CALL cp_fm_init_random(mo_coeff, nmo)
675 :
676 2052 : CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
677 : ! multiply times PS
678 2052 : IF (has_unit_metric) THEN
679 0 : CALL cp_fm_to_fm(mo_coeff, sv)
680 : ELSE
681 : ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
682 2052 : CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
683 : END IF
684 2052 : CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
685 :
686 2052 : CALL cp_fm_release(sv)
687 : ! and ortho the result
688 2052 : IF (has_unit_metric) THEN
689 0 : CALL make_basis_simple(mo_coeff, nmo)
690 : ELSE
691 2052 : 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 2074 : smear=qs_env%scf_control%smear)
697 :
698 : CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
699 2074 : mo_array(ispin)%mo_coeff_b) !fm->dbcsr
700 :
701 : CALL calculate_density_matrix(mo_array(ispin), &
702 2074 : 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 5254 : IF (scf_env%method == general_diag_method_nr) THEN
707 3424 : IF (dft_control%correct_surf_dip) THEN
708 8 : 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 6 : nmo=nmo, nao=nao, homo=homo)
712 :
713 6 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
714 6 : CALL cp_fm_init_random(mo_coeff, nmo)
715 :
716 6 : CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
717 : ! multiply times PS
718 6 : IF (has_unit_metric) THEN
719 0 : CALL cp_fm_to_fm(mo_coeff, sv)
720 : ELSE
721 : ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
722 6 : CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
723 : END IF
724 6 : CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
725 :
726 6 : CALL cp_fm_release(sv)
727 : ! and ortho the result
728 6 : IF (has_unit_metric) THEN
729 0 : CALL make_basis_simple(mo_coeff, nmo)
730 : ELSE
731 6 : 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 6 : tot_zeff_corr=tot_corr_zeff)
736 :
737 : CALL calculate_density_matrix(mo_array(ispin), &
738 6 : p_rmpv(ispin)%matrix)
739 : END IF
740 : END IF
741 : END IF
742 :
743 : END IF
744 :
745 : END DO
746 :
747 4309 : 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 0 : CALL qs_rho_update_rho(rho, qs_env)
750 4309 : 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 4309 : IF (scf_control%diagonalization%mom) THEN
763 4 : 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 4309 : "PRINT%KINDS")
768 :
769 4309 : did_guess = .TRUE.
770 : END IF
771 :
772 6751 : IF (density_guess == sparse_guess) THEN
773 :
774 0 : IF (ofgpw) CPABORT("SPARSE_GUESS not implemented for OFGPW")
775 0 : IF (.NOT. scf_control%use_ot) CPABORT("OT needed!")
776 0 : IF (do_kpoints) THEN
777 0 : CPABORT("calculate_first_density_matrix: sparse_guess not implemented for k-points")
778 : END IF
779 :
780 0 : eps = 1.0E-5_dp
781 :
782 0 : ounit = cp_logger_get_default_io_unit(logger)
783 0 : natoms = SIZE(particle_set)
784 0 : ALLOCATE (kind_of(natoms))
785 0 : ALLOCATE (first_sgf(natoms), last_sgf(natoms))
786 :
787 0 : checksum = dbcsr_checksum(s_sparse(1)%matrix)
788 0 : i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
789 0 : IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum
790 0 : CALL dbcsr_filter(s_sparse(1)%matrix, eps)
791 0 : checksum = dbcsr_checksum(s_sparse(1)%matrix)
792 0 : i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
793 0 : 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 0 : last_sgf=last_sgf)
797 0 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
798 :
799 0 : ALLOCATE (pmat(SIZE(atomic_kind_set)))
800 :
801 0 : rscale = 1._dp
802 0 : IF (nspin == 2) rscale = 0.5_dp
803 0 : DO ikind = 1, SIZE(atomic_kind_set)
804 0 : atomic_kind => atomic_kind_set(ikind)
805 0 : qs_kind => qs_kind_set(ikind)
806 0 : NULLIFY (pmat(ikind)%mat)
807 0 : CALL calculate_atomic_orbitals(atomic_kind, qs_kind, pmat=pmat(ikind)%mat)
808 0 : NULLIFY (atomic_kind)
809 : END DO
810 :
811 0 : DO ispin = 1, nspin
812 : CALL get_mo_set(mo_set=mo_array(ispin), &
813 : maxocc=maxocc, &
814 0 : nelectron=nelectron)
815 : !
816 0 : CALL dbcsr_iterator_start(iter, p_rmpv(ispin)%matrix)
817 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
818 0 : CALL dbcsr_iterator_next_block(iter, irow, icol, pdata)
819 0 : ikind = kind_of(irow)
820 0 : IF (icol .EQ. irow) THEN
821 0 : IF (ispin == 1) THEN
822 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
823 0 : pmat(ikind)%mat(:, :, 2)*rscale
824 : ELSE
825 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
826 0 : pmat(ikind)%mat(:, :, 2)*rscale
827 : END IF
828 : END IF
829 : END DO
830 0 : CALL dbcsr_iterator_stop(iter)
831 :
832 : !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
833 0 : checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
834 0 : occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
835 0 : 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 0 : checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
840 0 : occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
841 0 : IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum
842 :
843 0 : CALL dbcsr_dot(p_rmpv(ispin)%matrix, s_sparse(1)%matrix, trps1)
844 0 : rscale = REAL(nelectron, dp)/trps1
845 0 : CALL dbcsr_scale(p_rmpv(ispin)%matrix, rscale)
846 :
847 : !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
848 0 : checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
849 0 : occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
850 0 : 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 0 : IF (dft_control%restricted .AND. (ispin == 2)) THEN
855 0 : CALL mo_set_restrict(mo_array)
856 : ELSE
857 : CALL get_mo_set(mo_set=mo_array(ispin), &
858 : mo_coeff=mo_coeff, &
859 0 : nmo=nmo, nao=nao, homo=homo)
860 0 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
861 :
862 0 : n = MAXVAL(last_sgf - first_sgf) + 1
863 0 : size_atomic_kind_set = SIZE(atomic_kind_set)
864 :
865 0 : ALLOCATE (buff(n, n), sort_kind(size_atomic_kind_set), &
866 0 : nelec_kind(size_atomic_kind_set))
867 : !
868 : ! sort kind vs nbr electron
869 0 : DO ikind = 1, size_atomic_kind_set
870 0 : atomic_kind => atomic_kind_set(ikind)
871 0 : qs_kind => qs_kind_set(ikind)
872 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
873 : natom=natom, &
874 : atom_list=atom_list, &
875 0 : z=z)
876 : CALL get_qs_kind(qs_kind, nsgf=nsgf, elec_conf=elec_conf, &
877 0 : basis_set=orb_basis_set, zeff=zeff)
878 0 : nelec_kind(ikind) = SUM(elec_conf)
879 : END DO
880 0 : CALL sort(nelec_kind, size_atomic_kind_set, sort_kind)
881 : !
882 : ! a -very- naive sparse guess
883 0 : nmo_tmp = nmo
884 0 : natoms_tmp = natoms
885 0 : istart_col = 1
886 0 : iseed(1) = 4; iseed(2) = 3; iseed(3) = 2; iseed(4) = 1 ! set the seed for dlarnv
887 0 : DO i = 1, size_atomic_kind_set
888 0 : ikind = sort_kind(i)
889 0 : atomic_kind => atomic_kind_set(ikind)
890 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
891 0 : natom=natom, atom_list=atom_list)
892 0 : DO iatom = 1, natom
893 : !
894 0 : atom_a = atom_list(iatom)
895 0 : istart_row = first_sgf(atom_a)
896 0 : n_rows = last_sgf(atom_a) - first_sgf(atom_a) + 1
897 : !
898 : ! compute the "potential" nbr of states for this atom
899 0 : n_cols = MAX(INT(REAL(nmo_tmp, dp)/REAL(natoms_tmp, dp)), 1)
900 0 : IF (n_cols .GT. n_rows) n_cols = n_rows
901 : !
902 0 : nmo_tmp = nmo_tmp - n_cols
903 0 : natoms_tmp = natoms_tmp - 1
904 0 : IF (nmo_tmp .LT. 0 .OR. natoms_tmp .LT. 0) THEN
905 0 : CPABORT("Wrong1!")
906 : END IF
907 0 : DO j = 1, n_cols
908 0 : 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 0 : n_rows, n_cols)
912 0 : istart_col = istart_col + n_cols
913 : END DO
914 : END DO
915 :
916 0 : IF (istart_col .LE. nmo) THEN
917 0 : CPABORT("Wrong2!")
918 : END IF
919 :
920 0 : 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 0 : CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, mo_coeff, s_sparse(1)%matrix)
954 : !CALL dbcsr_verify_matrix(mo_dbcsr)
955 0 : checksum = dbcsr_checksum(mo_dbcsr)
956 :
957 0 : occ = dbcsr_get_occupation(mo_dbcsr)
958 0 : IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
959 0 : CALL dbcsr_filter(mo_dbcsr, eps)
960 : !CALL dbcsr_verify_matrix(mo_dbcsr)
961 0 : occ = dbcsr_get_occupation(mo_dbcsr)
962 0 : checksum = dbcsr_checksum(mo_dbcsr)
963 0 : IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
964 : !
965 : ! multiply times PS
966 0 : IF (has_unit_metric) THEN
967 0 : CPABORT("has_unit_metric will be removed soon")
968 : END IF
969 : !
970 : ! S*C
971 0 : 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 0 : retain_sparsity=.TRUE.)
975 : !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
976 0 : checksum = dbcsr_checksum(mo_tmp_dbcsr)
977 0 : occ = dbcsr_get_occupation(mo_tmp_dbcsr)
978 0 : IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum
979 0 : CALL dbcsr_filter(mo_tmp_dbcsr, eps)
980 : !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
981 0 : checksum = dbcsr_checksum(mo_tmp_dbcsr)
982 0 : occ = dbcsr_get_occupation(mo_tmp_dbcsr)
983 0 : 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 0 : mo_tmp_dbcsr, 0.0_dp, mo_dbcsr)
990 : IF (.FALSE.) CALL dbcsr_verify_matrix(mo_dbcsr)
991 0 : checksum = dbcsr_checksum(mo_dbcsr)
992 0 : occ = dbcsr_get_occupation(mo_dbcsr)
993 0 : IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
994 0 : CALL dbcsr_filter(mo_dbcsr, eps)
995 : !CALL dbcsr_verify_matrix(mo_dbcsr)
996 0 : checksum = dbcsr_checksum(mo_dbcsr)
997 0 : occ = dbcsr_get_occupation(mo_dbcsr)
998 0 : IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
999 : !
1000 0 : CALL copy_dbcsr_to_fm(mo_dbcsr, mo_coeff)
1001 :
1002 0 : CALL dbcsr_release(mo_dbcsr)
1003 0 : CALL dbcsr_release(mo_tmp_dbcsr)
1004 :
1005 : ! and ortho the result
1006 0 : 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 0 : smear=qs_env%scf_control%smear)
1011 :
1012 : CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1013 0 : mo_array(ispin)%mo_coeff_b) !fm->dbcsr
1014 :
1015 : CALL calculate_density_matrix(mo_array(ispin), &
1016 0 : p_rmpv(ispin)%matrix)
1017 0 : DO ikind = 1, SIZE(atomic_kind_set)
1018 0 : IF (ASSOCIATED(pmat(ikind)%mat)) THEN
1019 0 : DEALLOCATE (pmat(ikind)%mat)
1020 : END IF
1021 : END DO
1022 : END DO
1023 :
1024 0 : DEALLOCATE (pmat)
1025 :
1026 0 : DEALLOCATE (kind_of)
1027 :
1028 0 : DEALLOCATE (first_sgf, last_sgf)
1029 :
1030 0 : did_guess = .TRUE.
1031 : END IF
1032 6751 : 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 842 : nspin, nelectron_spin, para_env)
1037 :
1038 1754 : 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 1754 : IF (need_mos) THEN
1042 224 : IF (dft_control%restricted .AND. (ispin == 2)) THEN
1043 2 : CALL mo_set_restrict(mo_array)
1044 : ELSE
1045 : CALL get_mo_set(mo_set=mo_array(ispin), &
1046 : mo_coeff=mo_coeff, &
1047 222 : nmo=nmo, homo=homo)
1048 222 : CALL cp_fm_init_random(mo_coeff, nmo)
1049 222 : CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
1050 : ! multiply times PS
1051 222 : IF (has_unit_metric) THEN
1052 178 : CALL cp_fm_to_fm(mo_coeff, sv)
1053 : ELSE
1054 44 : 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 222 : CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
1058 222 : CALL cp_fm_release(sv)
1059 : ! and ortho the result
1060 222 : IF (has_unit_metric) THEN
1061 178 : CALL make_basis_simple(mo_coeff, nmo)
1062 : ELSE
1063 44 : 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 224 : smear=qs_env%scf_control%smear)
1069 : CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1070 224 : mo_array(ispin)%mo_coeff_b)
1071 :
1072 : CALL calculate_density_matrix(mo_array(ispin), &
1073 224 : 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 6751 : IF (density_guess == eht_guess) THEN
1082 4 : CALL calculate_eht_guess(qs_env, mo_array)
1083 8 : DO ispin = 1, nspin
1084 8 : 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 6751 : IF (dft_control%switch_surf_dip) THEN
1090 4 : DO ispin = 1, nspin
1091 : CALL reassign_allocated_mos(mos_last_converged(ispin), &
1092 4 : mo_array(ispin))
1093 : END DO
1094 : END IF
1095 :
1096 6751 : IF (density_guess == no_guess) THEN
1097 : did_guess = .TRUE.
1098 : END IF
1099 :
1100 5815 : IF (.NOT. did_guess) THEN
1101 0 : CPABORT("An invalid keyword for the initial density guess was specified")
1102 : END IF
1103 :
1104 6751 : CALL timestop(handle)
1105 :
1106 13502 : 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 98 : 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 98 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1125 98 : REAL(dp), DIMENSION(:, :), POINTER :: block
1126 98 : 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 98 : CALL timeset(routineN, handle)
1132 :
1133 98 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1134 432 : ALLOCATE (fmat(SIZE(atomic_kind_set)))
1135 :
1136 : ! precompute the atomic blocks for each atomic-kind
1137 236 : DO ikind = 1, SIZE(atomic_kind_set)
1138 138 : atomic_kind => atomic_kind_set(ikind)
1139 138 : qs_kind => qs_kind_set(ikind)
1140 138 : NULLIFY (fmat(ikind)%mat)
1141 138 : IF (ounit > 0) WRITE (UNIT=ounit, FMT="(/,T2,A)") &
1142 69 : "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 236 : fmat=fmat(ikind)%mat)
1147 : END DO
1148 :
1149 : ! zero result matrix
1150 98 : CALL dbcsr_set(matrix_f, 0.0_dp)
1151 :
1152 : ! copy precomputed blocks onto diagonal of result matrix
1153 98 : CALL dbcsr_iterator_start(iter, matrix_f)
1154 217 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1155 119 : CALL dbcsr_iterator_next_block(iter, irow, icol, block)
1156 119 : ikind = kind_of(irow)
1157 6937 : IF (icol .EQ. irow) block(:, :) = fmat(ikind)%mat(:, :, 1)
1158 : END DO
1159 98 : CALL dbcsr_iterator_stop(iter)
1160 :
1161 : ! cleanup
1162 236 : DO ikind = 1, SIZE(atomic_kind_set)
1163 236 : DEALLOCATE (fmat(ikind)%mat)
1164 : END DO
1165 98 : DEALLOCATE (fmat)
1166 :
1167 98 : CALL timestop(handle)
1168 :
1169 294 : 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 918 : SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
1185 : dft_control, particle_set, atomic_kind_set, qs_kind_set, &
1186 918 : 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 918 : INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nshell
1208 918 : 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 918 : 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 918 : CALL timeset(routineN, handle)
1222 :
1223 1906 : DO ispin = 1, nspin
1224 988 : matrix_p => pmat(ispin)%matrix
1225 1906 : CALL dbcsr_set(matrix_p, 0.0_dp)
1226 : END DO
1227 :
1228 918 : natom = SIZE(particle_set)
1229 918 : CALL dbcsr_get_info(pmat(1)%matrix, nfullrows_total=nao)
1230 918 : IF (nspin == 1) THEN
1231 : maxocc = 2.0_dp
1232 : ELSE
1233 70 : maxocc = 1.0_dp
1234 : END IF
1235 :
1236 2754 : ALLOCATE (first_sgf(natom))
1237 :
1238 918 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1239 918 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxl)
1240 :
1241 2754 : ALLOCATE (econf(0:maxl))
1242 :
1243 2754 : ALLOCATE (pdiag(nao))
1244 33922 : pdiag(:) = 0.0_dp
1245 :
1246 1836 : ALLOCATE (sdiag(nao))
1247 :
1248 33922 : sdiag(:) = 0.0_dp
1249 918 : IF (has_unit_metric) THEN
1250 12560 : sdiag(:) = 1.0_dp
1251 : ELSE
1252 556 : CALL dbcsr_get_diag(matrix_s, sdiag)
1253 556 : CALL para_env%sum(sdiag)
1254 : END IF
1255 :
1256 : ncount = 0
1257 33922 : trps1 = 0.0_dp
1258 33922 : trps2 = 0.0_dp
1259 33922 : pdiag(:) = 0.0_dp
1260 :
1261 2754 : IF (SUM(nelectron_spin) /= 0) THEN
1262 3038 : DO ikind = 1, SIZE(atomic_kind_set)
1263 :
1264 2134 : 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 2134 : sgp_potential=sgp_potential)
1269 2134 : has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
1270 :
1271 2134 : IF (dft_control%qs_control%dftb) THEN
1272 : CALL get_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter, &
1273 588 : lmax=maxll, occupation=edftb)
1274 588 : maxll = MIN(maxll, maxl)
1275 1764 : econf(0:maxl) = edftb(0:maxl)
1276 1546 : ELSEIF (dft_control%qs_control%xtb) THEN
1277 510 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1278 510 : CALL get_xtb_atom_param(xtb_kind, z=z, natorb=nsgf, nao=naox, lao=laox, occupation=occupation)
1279 1036 : ELSEIF (has_pot) THEN
1280 1036 : CALL get_atomic_kind(atomic_kind_set(ikind), z=z)
1281 1036 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, elec_conf=elec_conf, zeff=zeff)
1282 1036 : maxll = MIN(SIZE(elec_conf) - 1, maxl)
1283 3324 : econf(:) = 0.0_dp
1284 3262 : 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 5172 : IF (dft_control%qs_control%dftb) THEN
1291 4038 : DO iatom = 1, natom
1292 3450 : atom_a = atom_list(iatom)
1293 3450 : isgfa = first_sgf(atom_a)
1294 8738 : DO la = 0, maxll
1295 3450 : SELECT CASE (la)
1296 : CASE (0)
1297 3450 : pdiag(isgfa) = econf(0)
1298 : CASE (1)
1299 1250 : pdiag(isgfa + 1) = econf(1)/3._dp
1300 1250 : pdiag(isgfa + 2) = econf(1)/3._dp
1301 1250 : pdiag(isgfa + 3) = econf(1)/3._dp
1302 : CASE (2)
1303 0 : pdiag(isgfa + 4) = econf(2)/5._dp
1304 0 : pdiag(isgfa + 5) = econf(2)/5._dp
1305 0 : pdiag(isgfa + 6) = econf(2)/5._dp
1306 0 : pdiag(isgfa + 7) = econf(2)/5._dp
1307 0 : pdiag(isgfa + 8) = econf(2)/5._dp
1308 : CASE (3)
1309 0 : pdiag(isgfa + 9) = econf(3)/7._dp
1310 0 : pdiag(isgfa + 10) = econf(3)/7._dp
1311 0 : pdiag(isgfa + 11) = econf(3)/7._dp
1312 0 : pdiag(isgfa + 12) = econf(3)/7._dp
1313 0 : pdiag(isgfa + 13) = econf(3)/7._dp
1314 0 : pdiag(isgfa + 14) = econf(3)/7._dp
1315 0 : pdiag(isgfa + 15) = econf(3)/7._dp
1316 : CASE DEFAULT
1317 4700 : CPABORT("")
1318 : END SELECT
1319 : END DO
1320 : END DO
1321 1546 : ELSEIF (dft_control%qs_control%xtb) THEN
1322 3042 : DO iatom = 1, natom
1323 2532 : atom_a = atom_list(iatom)
1324 2532 : isgfa = first_sgf(atom_a)
1325 3042 : IF (z == 1 .AND. nsgf == 2) THEN
1326 : ! Hydrogen 2s basis
1327 1482 : pdiag(isgfa) = 1.0_dp
1328 1482 : pdiag(isgfa + 1) = 0.0_dp
1329 : ELSE
1330 6312 : DO isgf = 1, nsgf
1331 5262 : na = naox(isgf)
1332 5262 : la = laox(isgf)
1333 5262 : occ = REAL(occupation(la + 1), dp)/REAL(2*la + 1, dp)
1334 6312 : pdiag(isgfa + isgf - 1) = occ
1335 : END DO
1336 : END IF
1337 : END DO
1338 1036 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
1339 962 : yy = REAL(dft_control%charge, KIND=dp)/REAL(nao, KIND=dp)
1340 5482 : DO iatom = 1, natom
1341 4520 : atom_a = atom_list(iatom)
1342 4520 : isgfa = first_sgf(atom_a)
1343 962 : SELECT CASE (nsgf)
1344 : CASE (1) ! s-basis
1345 2188 : pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1346 : CASE (4) ! sp-basis
1347 2206 : IF (z == 1) THEN
1348 : ! special case: hydrogen with sp basis
1349 136 : pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1350 136 : pdiag(isgfa + 1) = 0._dp
1351 136 : pdiag(isgfa + 2) = 0._dp
1352 136 : pdiag(isgfa + 3) = 0._dp
1353 : ELSE
1354 2070 : pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1355 2070 : pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1356 2070 : pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1357 2070 : pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1358 : END IF
1359 : CASE (9) ! spd-basis
1360 126 : 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 92 : pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1363 92 : pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1364 92 : pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1365 92 : pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1366 92 : pdiag(isgfa + 4) = (-yy)*0.5_dp*maxocc
1367 92 : pdiag(isgfa + 5) = (-yy)*0.5_dp*maxocc
1368 92 : pdiag(isgfa + 6) = (-yy)*0.5_dp*maxocc
1369 92 : pdiag(isgfa + 7) = (-yy)*0.5_dp*maxocc
1370 92 : pdiag(isgfa + 8) = (-yy)*0.5_dp*maxocc
1371 34 : ELSE IF (z < 99) THEN
1372 34 : my_sum = zeff - 9.0_dp*yy
1373 : ! First, put 2 electrons in the 's' shell
1374 34 : pdiag(isgfa) = (MAX(0.0_dp, MIN(my_sum, 2.0_dp)))*0.5_dp*maxocc
1375 34 : my_sum = my_sum - 2.0_dp
1376 34 : IF (my_sum > 0.0_dp) THEN
1377 : ! Now put as many electrons as possible into the 'd' shell
1378 30 : pdiag(isgfa + 4) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1379 30 : pdiag(isgfa + 5) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1380 30 : pdiag(isgfa + 6) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1381 30 : pdiag(isgfa + 7) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1382 30 : pdiag(isgfa + 8) = (MAX(0.0_dp, MIN(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1383 30 : my_sum = MAX(0.0_dp, my_sum - 10.0_dp)
1384 : ! Put the remaining electrons in the 'p' shell
1385 30 : pdiag(isgfa + 1) = (my_sum/3.0_dp)*0.5_dp*maxocc
1386 30 : pdiag(isgfa + 2) = (my_sum/3.0_dp)*0.5_dp*maxocc
1387 30 : pdiag(isgfa + 3) = (my_sum/3.0_dp)*0.5_dp*maxocc
1388 : END IF
1389 : END IF
1390 : CASE DEFAULT
1391 4520 : 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 74 : last_sgf=last_sgfa)
1401 :
1402 202 : DO iset = 1, nset
1403 504 : DO ishell = 1, nshell(iset)
1404 302 : la = l(ishell, iset)
1405 302 : nelec = maxocc*REAL(2*la + 1, dp)
1406 430 : IF (econf(la) > 0.0_dp) THEN
1407 140 : IF (econf(la) >= nelec) THEN
1408 66 : paa = maxocc
1409 66 : econf(la) = econf(la) - nelec
1410 : ELSE
1411 74 : paa = maxocc*econf(la)/nelec
1412 74 : econf(la) = 0.0_dp
1413 : ncount = ncount + NINT(nelec/maxocc)
1414 : END IF
1415 412 : DO isgfa = first_sgfa(ishell, iset), last_sgfa(ishell, iset)
1416 2624 : DO iatom = 1, natom
1417 2212 : atom_a = atom_list(iatom)
1418 2212 : isgf = first_sgf(atom_a) + isgfa - 1
1419 2212 : pdiag(isgf) = paa
1420 2484 : IF (paa == maxocc) THEN
1421 550 : trps1 = trps1 + paa*sdiag(isgf)
1422 : ELSE
1423 1662 : 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 904 : IF (trps2 == 0.0_dp) THEN
1434 28420 : DO isgf = 1, nao
1435 28420 : IF (sdiag(isgf) > 0.0_dp) pdiag(isgf) = pdiag(isgf)/sdiag(isgf)
1436 : END DO
1437 1760 : DO ispin = 1, nspin
1438 1760 : IF (nelectron_spin(ispin) /= 0) THEN
1439 908 : matrix_p => pmat(ispin)%matrix
1440 908 : CALL dbcsr_set_diag(matrix_p, pdiag)
1441 : END IF
1442 : END DO
1443 : ELSE
1444 116 : DO ispin = 1, nspin
1445 116 : IF (nelectron_spin(ispin) /= 0) THEN
1446 60 : rscale = (REAL(nelectron_spin(ispin), dp) - trps1)/trps2
1447 5856 : DO isgf = 1, nao
1448 5856 : IF (pdiag(isgf) < maxocc) pdiag(isgf) = rscale*pdiag(isgf)
1449 : END DO
1450 60 : matrix_p => pmat(ispin)%matrix
1451 60 : CALL dbcsr_set_diag(matrix_p, pdiag)
1452 5856 : DO isgf = 1, nao
1453 5856 : 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 918 : DEALLOCATE (econf)
1461 :
1462 918 : DEALLOCATE (first_sgf)
1463 :
1464 918 : DEALLOCATE (pdiag)
1465 :
1466 918 : DEALLOCATE (sdiag)
1467 :
1468 918 : CALL timestop(handle)
1469 :
1470 2754 : END SUBROUTINE calculate_mopac_dm
1471 :
1472 0 : END MODULE qs_initial_guess
|
