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 for an energy correction on top of a Kohn-Sham calculation
10 : !> \par History
11 : !> 03.2026 created from energy_correction
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ec_diag_solver
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE basis_set_types, ONLY: gto_basis_set_type
17 : USE cp_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_scale, &
21 : dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
22 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 : copy_fm_to_dbcsr,&
24 : cp_dbcsr_sm_fm_multiply
25 : USE cp_fm_diag, ONLY: cp_fm_geeig
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27 : cp_fm_struct_release,&
28 : cp_fm_struct_type
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_init_random,&
31 : cp_fm_release,&
32 : cp_fm_set_all,&
33 : cp_fm_type
34 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr
35 : USE dm_ls_scf, ONLY: calculate_w_matrix_ls,&
36 : post_scf_sparsities
37 : USE dm_ls_scf_methods, ONLY: apply_matrix_preconditioner,&
38 : density_matrix_sign,&
39 : density_matrix_tc2,&
40 : density_matrix_trs4,&
41 : ls_scf_init_matrix_s
42 : USE dm_ls_scf_qs, ONLY: matrix_ls_create,&
43 : matrix_ls_to_qs,&
44 : matrix_qs_to_ls
45 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
46 : USE ec_env_types, ONLY: energy_correction_type
47 : USE ec_methods, ONLY: ec_mos_init
48 : USE input_constants, ONLY: ec_matrix_sign,&
49 : ec_matrix_tc2,&
50 : ec_matrix_trs4,&
51 : ec_ot_atomic,&
52 : ec_ot_gs,&
53 : ot_precond_full_single_inverse,&
54 : ot_precond_solver_default
55 : USE kinds, ONLY: dp
56 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
57 : kpoint_density_transform,&
58 : kpoint_set_mo_occupation
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE particle_types, ONLY: particle_type
61 : USE preconditioner, ONLY: make_preconditioner
62 : USE preconditioner_types, ONLY: destroy_preconditioner,&
63 : init_preconditioner,&
64 : preconditioner_type
65 : USE qs_atomic_block, ONLY: calculate_atomic_block_dm
66 : USE qs_density_matrices, ONLY: calculate_density_matrix,&
67 : calculate_w_matrix
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
70 : USE qs_kind_types, ONLY: get_qs_kind,&
71 : qs_kind_type
72 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,&
73 : make_basis_sm
74 : USE qs_mo_occupation, ONLY: set_mo_occupation
75 : USE qs_mo_types, ONLY: allocate_mo_set,&
76 : deallocate_mo_set,&
77 : get_mo_set,&
78 : init_mo_set,&
79 : mo_set_type
80 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
81 : USE qs_rho_types, ONLY: qs_rho_get,&
82 : qs_rho_type
83 : USE qs_scf_diagonalization, ONLY: diag_kp_basic
84 : USE string_utilities, ONLY: uppercase
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : ! Global parameters
92 :
93 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_diag_solver'
94 :
95 : PUBLIC :: ec_diag_solver_gamma, ec_diag_solver_kp, ec_ot_diag_solver, ec_ls_init, ec_ls_solver
96 :
97 : ! **************************************************************************************************
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief Solve KS equation using diagonalization
103 : !> \param qs_env ...
104 : !> \param ec_env ...
105 : !> \param matrix_ks ...
106 : !> \param matrix_s ...
107 : !> \param matrix_p ...
108 : !> \param matrix_w ...
109 : !> \par History
110 : !> 03.2014 created [JGH]
111 : !> \author JGH
112 : ! **************************************************************************************************
113 312 : SUBROUTINE ec_diag_solver_gamma(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
114 :
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 : TYPE(energy_correction_type) :: ec_env
117 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
118 :
119 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_diag_solver_gamma'
120 :
121 : INTEGER :: handle, ispin, nmo(2), nsize, nspins
122 : REAL(KIND=dp) :: eps_filter, flexible_electron_count, &
123 : focc(2), n_el_f
124 312 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigvals
125 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
126 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
127 : TYPE(cp_fm_type) :: fm_k, fm_s, fm_w
128 : TYPE(cp_fm_type), POINTER :: fm_mos
129 : TYPE(dbcsr_type), POINTER :: ref_matrix
130 : TYPE(dft_control_type), POINTER :: dft_control
131 312 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: moset
132 : TYPE(mp_para_env_type), POINTER :: para_env
133 :
134 312 : CALL timeset(routineN, handle)
135 :
136 312 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
137 312 : eps_filter = dft_control%qs_control%eps_filter_matrix
138 312 : nspins = dft_control%nspins
139 :
140 312 : nmo = 0
141 312 : CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
142 936 : focc = 1._dp
143 1248 : IF (nspins == 1) focc = 2._dp
144 :
145 312 : CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
146 :
147 312 : NULLIFY (blacs_env, para_env)
148 312 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
149 312 : NULLIFY (fm_struct)
150 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
151 312 : ncol_global=nsize, para_env=para_env)
152 312 : CALL cp_fm_create(fm_k, fm_struct)
153 312 : CALL cp_fm_create(fm_s, fm_struct)
154 312 : CALL cp_fm_create(fm_w, fm_struct)
155 312 : CALL cp_fm_struct_release(fm_struct)
156 :
157 : ! MO sets
158 312 : flexible_electron_count = 0.0_dp
159 312 : IF (ec_env%smear%do_smear) flexible_electron_count = 0.0001_dp
160 :
161 1248 : ALLOCATE (moset(nspins))
162 624 : DO ispin = 1, nspins
163 312 : n_el_f = nmo(ispin)
164 : CALL allocate_mo_set(moset(ispin), nsize, nsize, nmo(ispin), n_el_f, &
165 312 : focc(ispin), flexible_electron_count)
166 624 : CALL init_mo_set(moset(ispin), fm_ref=fm_w, name="MO")
167 : END DO
168 :
169 624 : DO ispin = 1, nspins
170 312 : ref_matrix => matrix_s(1, 1)%matrix
171 312 : CALL copy_dbcsr_to_fm(ref_matrix, fm_s)
172 312 : ref_matrix => matrix_ks(ispin, 1)%matrix
173 312 : CALL copy_dbcsr_to_fm(ref_matrix, fm_k)
174 312 : CALL get_mo_set(moset(ispin), eigenvalues=eigvals, mo_coeff=fm_mos)
175 624 : CALL cp_fm_geeig(fm_k, fm_s, fm_mos, eigvals, fm_w)
176 : END DO
177 : !
178 312 : CALL set_mo_occupation(moset, ec_env%smear)
179 : !
180 624 : DO ispin = 1, nspins
181 312 : ec_env%ekTS = ec_env%ekTS + moset(ispin)%kTS
182 312 : CALL calculate_density_matrix(moset(ispin), matrix_p(ispin, 1)%matrix)
183 624 : CALL calculate_w_matrix(moset(ispin), matrix_w(ispin, 1)%matrix)
184 : END DO
185 : !
186 312 : CALL cp_fm_release(fm_k)
187 312 : CALL cp_fm_release(fm_s)
188 312 : CALL cp_fm_release(fm_w)
189 624 : DO ispin = 1, nspins
190 624 : CALL deallocate_mo_set(moset(ispin))
191 : END DO
192 312 : DEALLOCATE (moset)
193 :
194 312 : CALL timestop(handle)
195 :
196 1248 : END SUBROUTINE ec_diag_solver_gamma
197 :
198 : ! **************************************************************************************************
199 : !> \brief Solve Kpoint-KS equation using diagonalization
200 : !> \param qs_env ...
201 : !> \param ec_env ...
202 : !> \param matrix_ks ...
203 : !> \param matrix_s ...
204 : !> \param matrix_p ...
205 : !> \param matrix_w ...
206 : !> \par History
207 : !> 03.2026 created [JGH]
208 : !> \author JGH
209 : ! **************************************************************************************************
210 10 : SUBROUTINE ec_diag_solver_kp(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
211 :
212 : TYPE(qs_environment_type), POINTER :: qs_env
213 : TYPE(energy_correction_type) :: ec_env
214 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
215 :
216 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_diag_solver_kp'
217 :
218 : INTEGER :: handle, i, ispin, nsize
219 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
220 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
221 10 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
222 : TYPE(dbcsr_type), POINTER :: tempmat
223 10 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
224 : TYPE(mp_para_env_type), POINTER :: para_env
225 :
226 10 : CALL timeset(routineN, handle)
227 :
228 10 : CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
229 :
230 10 : NULLIFY (blacs_env, para_env)
231 10 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
232 10 : NULLIFY (fm_struct)
233 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
234 10 : ncol_global=nsize, para_env=para_env)
235 50 : ALLOCATE (fmwork(4))
236 50 : DO i = 1, 4
237 50 : CALL cp_fm_create(fmwork(i), fm_struct)
238 : END DO
239 10 : CALL cp_fm_struct_release(fm_struct)
240 :
241 : ! Diagonalization
242 10 : CALL diag_kp_basic(matrix_ks, matrix_s, ec_env%kpoints, fmwork)
243 : ! MO occupations
244 10 : CALL kpoint_set_mo_occupation(ec_env%kpoints, ec_env%smear)
245 : ! density matrices
246 10 : CALL kpoint_density_matrices(ec_env%kpoints)
247 10 : CALL kpoint_density_matrices(ec_env%kpoints, .TRUE.)
248 : ! density matrices in real space
249 10 : tempmat => matrix_s(1, 1)%matrix
250 : CALL kpoint_density_transform(ec_env%kpoints, matrix_p, .FALSE., tempmat, &
251 10 : ec_env%sab_orb, fmwork)
252 : CALL kpoint_density_transform(ec_env%kpoints, matrix_w, .TRUE., tempmat, &
253 10 : ec_env%sab_orb, fmwork)
254 :
255 10 : ec_env%ekTS = 0.0_dp
256 10 : mos => ec_env%kpoints%kp_env(1)%kpoint_env%mos
257 20 : DO ispin = 1, SIZE(mos, 2)
258 20 : ec_env%ekTS = ec_env%ekTS + mos(1, ispin)%kTS
259 : END DO
260 :
261 10 : CALL cp_fm_release(fmwork)
262 :
263 10 : CALL timestop(handle)
264 :
265 10 : END SUBROUTINE ec_diag_solver_kp
266 :
267 : ! **************************************************************************************************
268 : !> \brief Use OT-diagonalziation to obtain density matrix from Harris Kohn-Sham matrix
269 : !> Initial guess of density matrix is either the atomic block initial guess from SCF
270 : !> or the ground-state density matrix. The latter only works if the same basis is used
271 : !>
272 : !> \param qs_env ...
273 : !> \param ec_env ...
274 : !> \param matrix_ks Harris Kohn-Sham matrix
275 : !> \param matrix_s Overlap matrix in Harris functional basis
276 : !> \param matrix_p Harris dentiy matrix, calculated here
277 : !> \param matrix_w Harris energy weighted density matrix, calculated here
278 : !>
279 : !> \par History
280 : !> 09.2020 created
281 : !> \author F.Belleflamme
282 : ! **************************************************************************************************
283 4 : SUBROUTINE ec_ot_diag_solver(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
284 : TYPE(qs_environment_type), POINTER :: qs_env
285 : TYPE(energy_correction_type), POINTER :: ec_env
286 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
287 : POINTER :: matrix_ks, matrix_s
288 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
289 : INTENT(INOUT), POINTER :: matrix_p, matrix_w
290 :
291 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_ot_diag_solver'
292 :
293 : INTEGER :: handle, homo, ikind, iounit, ispin, &
294 : max_iter, nao, nkind, nmo, nspins
295 : INTEGER, DIMENSION(2) :: nelectron_spin
296 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
297 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
298 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
299 : TYPE(cp_fm_type) :: sv
300 : TYPE(cp_fm_type), POINTER :: mo_coeff
301 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
302 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
303 : TYPE(dft_control_type), POINTER :: dft_control
304 : TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis
305 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
306 : TYPE(mp_para_env_type), POINTER :: para_env
307 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
308 : TYPE(preconditioner_type), POINTER :: local_preconditioner
309 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
310 : TYPE(qs_kind_type), POINTER :: qs_kind
311 : TYPE(qs_rho_type), POINTER :: rho
312 :
313 4 : CALL timeset(routineN, handle)
314 :
315 4 : iounit = cp_logger_get_default_unit_nr(local=.FALSE.)
316 :
317 4 : CPASSERT(ASSOCIATED(qs_env))
318 4 : CPASSERT(ASSOCIATED(ec_env))
319 4 : CPASSERT(ASSOCIATED(matrix_ks))
320 4 : CPASSERT(ASSOCIATED(matrix_s))
321 4 : CPASSERT(ASSOCIATED(matrix_p))
322 4 : CPASSERT(ASSOCIATED(matrix_w))
323 :
324 : CALL get_qs_env(qs_env=qs_env, &
325 : atomic_kind_set=atomic_kind_set, &
326 : blacs_env=blacs_env, &
327 : dft_control=dft_control, &
328 : nelectron_spin=nelectron_spin, &
329 : para_env=para_env, &
330 : particle_set=particle_set, &
331 4 : qs_kind_set=qs_kind_set)
332 4 : nspins = dft_control%nspins
333 :
334 : ! Maximum number of OT iterations for diagonalization
335 4 : max_iter = 200
336 :
337 : ! If linear scaling, need to allocate and init MO set
338 : ! set it to qs_env%mos
339 4 : IF (dft_control%qs_control%do_ls_scf) THEN
340 0 : CALL ec_mos_init(qs_env, matrix_s(1, 1)%matrix)
341 : END IF
342 :
343 4 : CALL get_qs_env(qs_env, mos=mos)
344 :
345 : ! Inital guess to use
346 4 : NULLIFY (p_rmpv)
347 :
348 : ! Using ether ground-state DM or ATOMIC GUESS requires
349 : ! Harris functional to use the same basis set
350 4 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
351 4 : CALL uppercase(ec_env%basis)
352 : ! Harris basis only differs from ground-state basis if explicitly added
353 : ! thus only two cases that need to be tested
354 : ! 1) explicit Harris basis present?
355 4 : IF (ec_env%basis == "HARRIS") THEN
356 12 : DO ikind = 1, nkind
357 8 : qs_kind => qs_kind_set(ikind)
358 : ! Basis sets of ground-state
359 8 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
360 : ! Basis sets of energy correction
361 8 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
362 :
363 12 : IF (basis_set%name /= harris_basis%name) THEN
364 0 : CPABORT("OT-Diag initial guess: Harris and ground state need to use the same basis")
365 : END IF
366 : END DO
367 : END IF
368 : ! 2) Harris uses MAOs
369 4 : IF (ec_env%mao) THEN
370 0 : CPABORT("OT-Diag initial guess: not implemented for MAOs")
371 : END IF
372 :
373 : ! Initital guess obtained for OT Diagonalization
374 6 : SELECT CASE (ec_env%ec_initial_guess)
375 : CASE (ec_ot_atomic)
376 :
377 2 : p_rmpv => matrix_p(:, 1)
378 :
379 : CALL calculate_atomic_block_dm(p_rmpv, matrix_s(1, 1)%matrix, atomic_kind_set, qs_kind_set, &
380 2 : nspins, nelectron_spin, iounit, para_env)
381 :
382 : CASE (ec_ot_gs)
383 :
384 2 : CALL get_qs_env(qs_env, rho=rho)
385 2 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
386 2 : p_rmpv => rho_ao(:, 1)
387 :
388 : CASE DEFAULT
389 4 : CPABORT("Unknown inital guess for OT-Diagonalization (Harris functional)")
390 : END SELECT
391 :
392 8 : DO ispin = 1, nspins
393 : CALL get_mo_set(mo_set=mos(ispin), &
394 : mo_coeff=mo_coeff, &
395 : nmo=nmo, &
396 : nao=nao, &
397 4 : homo=homo)
398 :
399 : ! Calculate first MOs
400 4 : CALL cp_fm_set_all(mo_coeff, 0.0_dp)
401 4 : CALL cp_fm_init_random(mo_coeff, nmo)
402 :
403 4 : CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
404 : ! multiply times PS
405 : ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
406 4 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff, sv, nmo)
407 4 : CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
408 4 : CALL cp_fm_release(sv)
409 : ! and ortho the result
410 : ! If DFBT or SE, then needs has_unit_metrix option
411 16 : CALL make_basis_sm(mo_coeff, nmo, matrix_s(1, 1)%matrix)
412 : END DO
413 :
414 : ! Preconditioner
415 : NULLIFY (local_preconditioner)
416 4 : ALLOCATE (local_preconditioner)
417 : CALL init_preconditioner(local_preconditioner, para_env=para_env, &
418 4 : blacs_env=blacs_env)
419 8 : DO ispin = 1, nspins
420 : CALL make_preconditioner(local_preconditioner, &
421 : precon_type=ot_precond_full_single_inverse, &
422 : solver_type=ot_precond_solver_default, &
423 : matrix_h=matrix_ks(ispin, 1)%matrix, &
424 : matrix_s=matrix_s(ispin, 1)%matrix, &
425 : convert_precond_to_dbcsr=.TRUE., &
426 4 : mo_set=mos(ispin), energy_gap=0.2_dp)
427 :
428 : CALL get_mo_set(mos(ispin), &
429 : mo_coeff=mo_coeff, &
430 : eigenvalues=eigenvalues, &
431 : nmo=nmo, &
432 4 : homo=homo)
433 : CALL ot_eigensolver(matrix_h=matrix_ks(ispin, 1)%matrix, &
434 : matrix_s=matrix_s(1, 1)%matrix, &
435 : matrix_c_fm=mo_coeff, &
436 : preconditioner=local_preconditioner, &
437 : eps_gradient=ec_env%eps_default, &
438 : iter_max=max_iter, &
439 4 : silent=.FALSE.)
440 : CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin, 1)%matrix, &
441 4 : evals_arg=eigenvalues, do_rotation=.TRUE.)
442 :
443 : ! Deallocate preconditioner
444 4 : CALL destroy_preconditioner(local_preconditioner)
445 4 : DEALLOCATE (local_preconditioner)
446 :
447 : !fm->dbcsr
448 : CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
449 12 : mos(ispin)%mo_coeff_b)
450 : END DO
451 :
452 : ! Calculate density matrix from MOs
453 8 : DO ispin = 1, nspins
454 4 : CALL calculate_density_matrix(mos(ispin), matrix_p(ispin, 1)%matrix)
455 :
456 8 : CALL calculate_w_matrix(mos(ispin), matrix_w(ispin, 1)%matrix)
457 : END DO
458 :
459 : ! Get rid of MO environment again
460 4 : IF (dft_control%qs_control%do_ls_scf) THEN
461 0 : DO ispin = 1, nspins
462 0 : CALL deallocate_mo_set(mos(ispin))
463 : END DO
464 0 : IF (ASSOCIATED(qs_env%mos)) THEN
465 0 : DO ispin = 1, SIZE(qs_env%mos)
466 0 : CALL deallocate_mo_set(qs_env%mos(ispin))
467 : END DO
468 0 : DEALLOCATE (qs_env%mos)
469 : END IF
470 : END IF
471 :
472 4 : CALL timestop(handle)
473 :
474 4 : END SUBROUTINE ec_ot_diag_solver
475 :
476 : ! **************************************************************************************************
477 : !> \brief Solve the Harris functional by linear scaling density purification scheme,
478 : !> instead of the diagonalization performed in ec_diag_solver
479 : !>
480 : !> \param qs_env ...
481 : !> \param matrix_ks Harris Kohn-Sham matrix
482 : !> \param matrix_s Overlap matrix in Harris functional basis
483 : !> \par History
484 : !> 09.2020 created
485 : !> \author F.Belleflamme
486 : ! **************************************************************************************************
487 30 : SUBROUTINE ec_ls_init(qs_env, matrix_ks, matrix_s)
488 : TYPE(qs_environment_type), POINTER :: qs_env
489 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
490 :
491 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_ls_init'
492 :
493 : INTEGER :: handle, ispin, nspins
494 : TYPE(dft_control_type), POINTER :: dft_control
495 : TYPE(energy_correction_type), POINTER :: ec_env
496 : TYPE(ls_scf_env_type), POINTER :: ls_env
497 :
498 30 : CALL timeset(routineN, handle)
499 :
500 : CALL get_qs_env(qs_env=qs_env, &
501 : dft_control=dft_control, &
502 30 : ec_env=ec_env)
503 30 : nspins = dft_control%nspins
504 30 : ls_env => ec_env%ls_env
505 :
506 : ! create the matrix template for use in the ls procedures
507 : CALL matrix_ls_create(matrix_ls=ls_env%matrix_s, matrix_qs=matrix_s(1, 1)%matrix, &
508 30 : ls_mstruct=ls_env%ls_mstruct)
509 :
510 30 : IF (ALLOCATED(ls_env%matrix_p)) THEN
511 16 : DO ispin = 1, SIZE(ls_env%matrix_p)
512 16 : CALL dbcsr_release(ls_env%matrix_p(ispin))
513 : END DO
514 : ELSE
515 88 : ALLOCATE (ls_env%matrix_p(nspins))
516 : END IF
517 :
518 60 : DO ispin = 1, nspins
519 : CALL dbcsr_create(ls_env%matrix_p(ispin), template=ls_env%matrix_s, &
520 60 : matrix_type=dbcsr_type_no_symmetry)
521 : END DO
522 :
523 120 : ALLOCATE (ls_env%matrix_ks(nspins))
524 60 : DO ispin = 1, nspins
525 : CALL dbcsr_create(ls_env%matrix_ks(ispin), template=ls_env%matrix_s, &
526 60 : matrix_type=dbcsr_type_no_symmetry)
527 : END DO
528 :
529 : ! Set up S matrix and needed functions of S
530 30 : CALL ls_scf_init_matrix_s(matrix_s(1, 1)%matrix, ls_env)
531 :
532 : ! Bring KS matrix from QS to LS form
533 : ! EC KS-matrix already calculated
534 60 : DO ispin = 1, nspins
535 : CALL matrix_qs_to_ls(matrix_ls=ls_env%matrix_ks(ispin), &
536 : matrix_qs=matrix_ks(ispin, 1)%matrix, &
537 : ls_mstruct=ls_env%ls_mstruct, &
538 60 : covariant=.TRUE.)
539 : END DO
540 :
541 30 : CALL timestop(handle)
542 :
543 30 : END SUBROUTINE ec_ls_init
544 :
545 : ! **************************************************************************************************
546 : !> \brief Solve the Harris functional by linear scaling density purification scheme,
547 : !> instead of the diagonalization performed in ec_diag_solver
548 : !>
549 : !> \param qs_env ...
550 : !> \param matrix_p Harris dentiy matrix, calculated here
551 : !> \param matrix_w Harris energy weighted density matrix, calculated here
552 : !> \param ec_ls_method which purification scheme should be used
553 : !> \par History
554 : !> 12.2019 created [JGH]
555 : !> 08.2020 refactoring [fbelle]
556 : !> \author Fabian Belleflamme
557 : ! **************************************************************************************************
558 :
559 30 : SUBROUTINE ec_ls_solver(qs_env, matrix_p, matrix_w, ec_ls_method)
560 : TYPE(qs_environment_type), POINTER :: qs_env
561 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
562 : INTEGER, INTENT(IN) :: ec_ls_method
563 :
564 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ls_solver'
565 :
566 : INTEGER :: handle, ispin, nelectron_spin_real, &
567 : nspins
568 : INTEGER, DIMENSION(2) :: nmo
569 30 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: wmat
570 30 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_ks_deviation
571 : TYPE(dft_control_type), POINTER :: dft_control
572 : TYPE(energy_correction_type), POINTER :: ec_env
573 : TYPE(ls_scf_env_type), POINTER :: ls_env
574 : TYPE(mp_para_env_type), POINTER :: para_env
575 :
576 30 : CALL timeset(routineN, handle)
577 :
578 30 : NULLIFY (para_env)
579 : CALL get_qs_env(qs_env, &
580 : dft_control=dft_control, &
581 30 : para_env=para_env)
582 30 : nspins = dft_control%nspins
583 30 : ec_env => qs_env%ec_env
584 30 : ls_env => ec_env%ls_env
585 :
586 30 : nmo = 0
587 30 : CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
588 30 : IF (nspins == 1) nmo(1) = nmo(1)/2
589 90 : ls_env%homo_spin(:) = 0.0_dp
590 90 : ls_env%lumo_spin(:) = 0.0_dp
591 :
592 120 : ALLOCATE (matrix_ks_deviation(nspins))
593 60 : DO ispin = 1, nspins
594 30 : CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_env%matrix_ks(ispin))
595 60 : CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
596 : END DO
597 :
598 : ! F = S^-1/2 * F * S^-1/2
599 30 : IF (ls_env%has_s_preconditioner) THEN
600 60 : DO ispin = 1, nspins
601 : CALL apply_matrix_preconditioner(ls_env%matrix_ks(ispin), "forward", &
602 30 : ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
603 :
604 60 : CALL dbcsr_filter(ls_env%matrix_ks(ispin), ls_env%eps_filter)
605 : END DO
606 : END IF
607 :
608 60 : DO ispin = 1, nspins
609 30 : nelectron_spin_real = ls_env%nelectron_spin(ispin)
610 30 : IF (ls_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
611 :
612 30 : SELECT CASE (ec_ls_method)
613 : CASE (ec_matrix_sign)
614 : CALL density_matrix_sign(ls_env%matrix_p(ispin), &
615 : ls_env%mu_spin(ispin), &
616 : ls_env%fixed_mu, &
617 : ls_env%sign_method, &
618 : ls_env%sign_order, &
619 : ls_env%matrix_ks(ispin), &
620 : ls_env%matrix_s, &
621 : ls_env%matrix_s_inv, &
622 : nelectron_spin_real, &
623 2 : ec_env%eps_default)
624 :
625 : CASE (ec_matrix_trs4)
626 : CALL density_matrix_trs4( &
627 : ls_env%matrix_p(ispin), &
628 : ls_env%matrix_ks(ispin), &
629 : ls_env%matrix_s_sqrt_inv, &
630 : nelectron_spin_real, &
631 : ec_env%eps_default, &
632 : ls_env%homo_spin(ispin), &
633 : ls_env%lumo_spin(ispin), &
634 : ls_env%mu_spin(ispin), &
635 : matrix_ks_deviation=matrix_ks_deviation(ispin), &
636 : dynamic_threshold=ls_env%dynamic_threshold, &
637 : eps_lanczos=ls_env%eps_lanczos, &
638 26 : max_iter_lanczos=ls_env%max_iter_lanczos)
639 :
640 : CASE (ec_matrix_tc2)
641 : CALL density_matrix_tc2( &
642 : ls_env%matrix_p(ispin), &
643 : ls_env%matrix_ks(ispin), &
644 : ls_env%matrix_s_sqrt_inv, &
645 : nelectron_spin_real, &
646 : ec_env%eps_default, &
647 : ls_env%homo_spin(ispin), &
648 : ls_env%lumo_spin(ispin), &
649 : non_monotonic=ls_env%non_monotonic, &
650 : eps_lanczos=ls_env%eps_lanczos, &
651 30 : max_iter_lanczos=ls_env%max_iter_lanczos)
652 :
653 : END SELECT
654 :
655 : END DO
656 :
657 : ! de-orthonormalize
658 30 : IF (ls_env%has_s_preconditioner) THEN
659 60 : DO ispin = 1, nspins
660 : ! P = S^-1/2 * P_tilde * S^-1/2 (forward)
661 : CALL apply_matrix_preconditioner(ls_env%matrix_p(ispin), "forward", &
662 30 : ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
663 :
664 60 : CALL dbcsr_filter(ls_env%matrix_p(ispin), ls_env%eps_filter)
665 : END DO
666 : END IF
667 :
668 : ! Closed-shell
669 30 : IF (nspins == 1) CALL dbcsr_scale(ls_env%matrix_p(1), 2.0_dp)
670 :
671 30 : IF (ls_env%report_all_sparsities) CALL post_scf_sparsities(ls_env)
672 :
673 : ! ls_scf_dm_to_ks
674 : ! Density matrix from LS to EC
675 60 : DO ispin = 1, nspins
676 : CALL matrix_ls_to_qs(matrix_qs=matrix_p(ispin, 1)%matrix, &
677 : matrix_ls=ls_env%matrix_p(ispin), &
678 : ls_mstruct=ls_env%ls_mstruct, &
679 60 : covariant=.FALSE.)
680 : END DO
681 :
682 30 : wmat => matrix_w(:, 1)
683 30 : CALL calculate_w_matrix_ls(wmat, ec_env%ls_env)
684 :
685 : ! clean up
686 30 : CALL dbcsr_release(ls_env%matrix_s)
687 30 : IF (ls_env%has_s_preconditioner) THEN
688 30 : CALL dbcsr_release(ls_env%matrix_bs_sqrt)
689 30 : CALL dbcsr_release(ls_env%matrix_bs_sqrt_inv)
690 : END IF
691 30 : IF (ls_env%needs_s_inv) THEN
692 2 : CALL dbcsr_release(ls_env%matrix_s_inv)
693 : END IF
694 30 : IF (ls_env%use_s_sqrt) THEN
695 30 : CALL dbcsr_release(ls_env%matrix_s_sqrt)
696 30 : CALL dbcsr_release(ls_env%matrix_s_sqrt_inv)
697 : END IF
698 :
699 60 : DO ispin = 1, SIZE(ls_env%matrix_ks)
700 60 : CALL dbcsr_release(ls_env%matrix_ks(ispin))
701 : END DO
702 30 : DEALLOCATE (ls_env%matrix_ks)
703 :
704 60 : DO ispin = 1, nspins
705 60 : CALL dbcsr_release(matrix_ks_deviation(ispin))
706 : END DO
707 30 : DEALLOCATE (matrix_ks_deviation)
708 :
709 30 : CALL timestop(handle)
710 :
711 30 : END SUBROUTINE ec_ls_solver
712 :
713 : END MODULE ec_diag_solver
714 :
|
