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