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 collects routines that perform operations directly related to MOs
10 : !> \note
11 : !> first version : most routines imported
12 : !> \author Joost VandeVondele (2003-08)
13 : ! **************************************************************************************************
14 : MODULE qs_mo_methods
15 : USE admm_types, ONLY: admm_type
16 : USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
17 : admm_uncorrect_for_eigenvalues
18 : USE cp_blacs_env, ONLY: cp_blacs_env_type
19 : USE cp_control_types, ONLY: hairy_probes_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
21 : dbcsr_get_info,&
22 : dbcsr_init_p,&
23 : dbcsr_multiply,&
24 : dbcsr_p_type,&
25 : dbcsr_release_p,&
26 : dbcsr_type,&
27 : dbcsr_type_no_symmetry
28 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd,&
29 : cp_dbcsr_syevx
30 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
31 : copy_fm_to_dbcsr,&
32 : cp_dbcsr_m_by_n_from_template,&
33 : cp_dbcsr_sm_fm_multiply
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_syrk,&
35 : cp_fm_triangular_multiply
36 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
37 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
38 : cp_fm_power,&
39 : cp_fm_syevx
40 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
41 : cp_fm_struct_release,&
42 : cp_fm_struct_type
43 : USE cp_fm_types, ONLY: cp_fm_create,&
44 : cp_fm_get_info,&
45 : cp_fm_get_submatrix,&
46 : cp_fm_release,&
47 : cp_fm_set_all,&
48 : cp_fm_to_fm,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
51 : USE kinds, ONLY: dp
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE parallel_gemm_api, ONLY: parallel_gemm
54 : USE physcon, ONLY: evolt
55 : USE qs_mo_occupation, ONLY: set_mo_occupation
56 : USE qs_mo_types, ONLY: get_mo_set,&
57 : mo_set_type
58 : USE scf_control_types, ONLY: scf_control_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 : PRIVATE
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_methods'
64 :
65 : PUBLIC :: make_basis_simple, make_basis_cholesky, make_basis_sv, make_basis_sm, &
66 : make_basis_lowdin, calculate_subspace_eigenvalues, &
67 : calculate_orthonormality, calculate_magnitude, make_mo_eig, &
68 : calculate_orbital_fock_elements, make_mo_set_canonical
69 :
70 : INTERFACE calculate_subspace_eigenvalues
71 : MODULE PROCEDURE subspace_eigenvalues_ks_fm
72 : MODULE PROCEDURE subspace_eigenvalues_ks_dbcsr
73 : END INTERFACE
74 :
75 : INTERFACE make_basis_sv
76 : MODULE PROCEDURE make_basis_sv_fm
77 : MODULE PROCEDURE make_basis_sv_dbcsr
78 : END INTERFACE
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief returns an S-orthonormal basis v (v^T S v ==1)
84 : !> \param vmatrix ...
85 : !> \param ncol ...
86 : !> \param matrix_s ...
87 : !> \par History
88 : !> 03.2006 created [Joost VandeVondele]
89 : ! **************************************************************************************************
90 31456 : SUBROUTINE make_basis_sm(vmatrix, ncol, matrix_s)
91 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
92 : INTEGER, INTENT(IN) :: ncol
93 : TYPE(dbcsr_type), POINTER :: matrix_s
94 :
95 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sm'
96 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
97 :
98 : INTEGER :: handle, n, ncol_global
99 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
100 : TYPE(cp_fm_type) :: overlap_vv, svmatrix
101 :
102 304 : IF (ncol == 0) RETURN
103 :
104 7788 : CALL timeset(routineN, handle)
105 :
106 7788 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
107 7788 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
108 :
109 7788 : CALL cp_fm_create(svmatrix, vmatrix%matrix_struct, "SV", set_zero=.TRUE.)
110 7788 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, svmatrix, ncol)
111 :
112 7788 : NULLIFY (fm_struct_tmp)
113 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
114 : para_env=vmatrix%matrix_struct%para_env, &
115 7788 : context=vmatrix%matrix_struct%context)
116 7788 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
117 7788 : CALL cp_fm_struct_release(fm_struct_tmp)
118 :
119 7788 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
120 7788 : CALL cp_fm_cholesky_decompose(overlap_vv)
121 7788 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
122 :
123 7788 : CALL cp_fm_release(overlap_vv)
124 7788 : CALL cp_fm_release(svmatrix)
125 :
126 7788 : CALL timestop(handle)
127 :
128 8092 : END SUBROUTINE make_basis_sm
129 :
130 : ! **************************************************************************************************
131 : !> \brief returns an S-orthonormal basis v and the corresponding matrix S*v as well
132 : !> \param vmatrix ...
133 : !> \param ncol ...
134 : !> \param svmatrix ...
135 : !> \par History
136 : !> 03.2006 created [Joost VandeVondele]
137 : ! **************************************************************************************************
138 0 : SUBROUTINE make_basis_sv_fm(vmatrix, ncol, svmatrix)
139 :
140 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
141 : INTEGER, INTENT(IN) :: ncol
142 : TYPE(cp_fm_type), INTENT(IN) :: svmatrix
143 :
144 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sv_fm'
145 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
146 :
147 : INTEGER :: handle, n, ncol_global
148 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
149 : TYPE(cp_fm_type) :: overlap_vv
150 :
151 0 : IF (ncol == 0) RETURN
152 :
153 0 : CALL timeset(routineN, handle)
154 0 : NULLIFY (fm_struct_tmp)
155 :
156 0 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
157 0 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
158 :
159 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
160 : para_env=vmatrix%matrix_struct%para_env, &
161 0 : context=vmatrix%matrix_struct%context)
162 0 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
163 0 : CALL cp_fm_struct_release(fm_struct_tmp)
164 :
165 0 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
166 0 : CALL cp_fm_cholesky_decompose(overlap_vv)
167 0 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
168 0 : CALL cp_fm_triangular_multiply(overlap_vv, svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
169 :
170 0 : CALL cp_fm_release(overlap_vv)
171 :
172 0 : CALL timestop(handle)
173 :
174 0 : END SUBROUTINE make_basis_sv_fm
175 :
176 : ! **************************************************************************************************
177 : !> \brief ...
178 : !> \param vmatrix ...
179 : !> \param ncol ...
180 : !> \param svmatrix ...
181 : !> \param para_env ...
182 : !> \param blacs_env ...
183 : ! **************************************************************************************************
184 3686 : SUBROUTINE make_basis_sv_dbcsr(vmatrix, ncol, svmatrix, para_env, blacs_env)
185 :
186 : TYPE(dbcsr_type) :: vmatrix
187 : INTEGER, INTENT(IN) :: ncol
188 : TYPE(dbcsr_type) :: svmatrix
189 : TYPE(mp_para_env_type), POINTER :: para_env
190 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
191 :
192 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sv_dbcsr'
193 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
194 :
195 : INTEGER :: handle, n, ncol_global
196 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
197 : TYPE(cp_fm_type) :: fm_svmatrix, fm_vmatrix, overlap_vv
198 :
199 122 : IF (ncol == 0) RETURN
200 :
201 594 : CALL timeset(routineN, handle)
202 :
203 : !CALL cp_fm_get_info(matrix=vmatrix,nrow_global=n,ncol_global=ncol_global)
204 594 : CALL dbcsr_get_info(vmatrix, nfullrows_total=n, nfullcols_total=ncol_global)
205 594 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
206 :
207 : CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=ncol, &
208 594 : ncol_global=ncol, para_env=para_env)
209 594 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, name="fm_overlap_vv")
210 594 : CALL cp_fm_struct_release(fm_struct_tmp)
211 :
212 : CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=n, &
213 594 : ncol_global=ncol_global, para_env=para_env)
214 594 : CALL cp_fm_create(fm_vmatrix, fm_struct_tmp, name="fm_vmatrix")
215 594 : CALL cp_fm_create(fm_svmatrix, fm_struct_tmp, name="fm_svmatrix")
216 594 : CALL cp_fm_struct_release(fm_struct_tmp)
217 :
218 594 : CALL copy_dbcsr_to_fm(vmatrix, fm_vmatrix)
219 594 : CALL copy_dbcsr_to_fm(svmatrix, fm_svmatrix)
220 :
221 594 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, fm_vmatrix, fm_svmatrix, rzero, overlap_vv)
222 594 : CALL cp_fm_cholesky_decompose(overlap_vv)
223 594 : CALL cp_fm_triangular_multiply(overlap_vv, fm_vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
224 594 : CALL cp_fm_triangular_multiply(overlap_vv, fm_svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
225 :
226 594 : CALL copy_fm_to_dbcsr(fm_vmatrix, vmatrix)
227 594 : CALL copy_fm_to_dbcsr(fm_svmatrix, svmatrix)
228 :
229 594 : CALL cp_fm_release(overlap_vv)
230 594 : CALL cp_fm_release(fm_vmatrix)
231 594 : CALL cp_fm_release(fm_svmatrix)
232 :
233 594 : CALL timestop(handle)
234 :
235 716 : END SUBROUTINE make_basis_sv_dbcsr
236 :
237 : ! **************************************************************************************************
238 : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
239 : !> the cholesky decomposed form of S is passed as an argument
240 : !> \param vmatrix ...
241 : !> \param ncol ...
242 : !> \param ortho cholesky decomposed S matrix
243 : !> \par History
244 : !> 03.2006 created [Joost VandeVondele]
245 : !> \note
246 : !> if the cholesky decomposed S matrix is not available
247 : !> use make_basis_sm since this is much faster than computing the
248 : !> cholesky decomposition of S
249 : ! **************************************************************************************************
250 33012 : SUBROUTINE make_basis_cholesky(vmatrix, ncol, ortho)
251 :
252 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
253 : INTEGER, INTENT(IN) :: ncol
254 : TYPE(cp_fm_type), INTENT(IN) :: ortho
255 :
256 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_cholesky'
257 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
258 :
259 : INTEGER :: handle, n, ncol_global
260 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
261 : TYPE(cp_fm_type) :: overlap_vv
262 :
263 84 : IF (ncol == 0) RETURN
264 :
265 8232 : CALL timeset(routineN, handle)
266 8232 : NULLIFY (fm_struct_tmp)
267 :
268 8232 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
269 8232 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
270 :
271 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
272 : para_env=vmatrix%matrix_struct%para_env, &
273 8232 : context=vmatrix%matrix_struct%context)
274 8232 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
275 8232 : CALL cp_fm_struct_release(fm_struct_tmp)
276 :
277 8232 : CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol)
278 8232 : CALL cp_fm_syrk('U', 'T', n, rone, vmatrix, 1, 1, rzero, overlap_vv)
279 8232 : CALL cp_fm_cholesky_decompose(overlap_vv)
280 8232 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
281 8232 : CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol, invert_tr=.TRUE.)
282 :
283 8232 : CALL cp_fm_release(overlap_vv)
284 :
285 8232 : CALL timestop(handle)
286 :
287 8316 : END SUBROUTINE make_basis_cholesky
288 :
289 : ! **************************************************************************************************
290 : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
291 : !> a Loedwin transformation is applied to keep the rotated vectors as close
292 : !> as possible to the original ones
293 : !> \param vmatrix ...
294 : !> \param ncol ...
295 : !> \param matrix_s ...
296 : !> \param
297 : !> \par History
298 : !> 05.2009 created [MI]
299 : !> \note
300 : ! **************************************************************************************************
301 12540 : SUBROUTINE make_basis_lowdin(vmatrix, ncol, matrix_s)
302 :
303 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
304 : INTEGER, INTENT(IN) :: ncol
305 : TYPE(dbcsr_type) :: matrix_s
306 :
307 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_lowdin'
308 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
309 :
310 : INTEGER :: handle, n, ncol_global, ndep
311 : REAL(dp) :: threshold
312 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
313 : TYPE(cp_fm_type) :: csc, sc, work
314 :
315 0 : IF (ncol == 0) RETURN
316 :
317 2090 : CALL timeset(routineN, handle)
318 2090 : NULLIFY (fm_struct_tmp)
319 2090 : threshold = 1.0E-7_dp
320 2090 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
321 2090 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
322 :
323 2090 : CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
324 2090 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
325 :
326 2090 : NULLIFY (fm_struct_tmp)
327 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
328 : para_env=vmatrix%matrix_struct%para_env, &
329 2090 : context=vmatrix%matrix_struct%context)
330 2090 : CALL cp_fm_create(csc, fm_struct_tmp, "csc")
331 2090 : CALL cp_fm_set_all(csc, 0.0_dp)
332 2090 : CALL cp_fm_create(work, fm_struct_tmp, "work")
333 2090 : CALL cp_fm_set_all(work, 0.0_dp)
334 2090 : CALL cp_fm_struct_release(fm_struct_tmp)
335 :
336 2090 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
337 2090 : CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
338 2090 : CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
339 2090 : CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
340 :
341 2090 : CALL cp_fm_release(csc)
342 2090 : CALL cp_fm_release(sc)
343 2090 : CALL cp_fm_release(work)
344 :
345 2090 : CALL timestop(handle)
346 :
347 2090 : END SUBROUTINE make_basis_lowdin
348 :
349 : ! **************************************************************************************************
350 : !> \brief given a set of vectors, return an orthogonal (C^T C == 1) set
351 : !> spanning the same space (notice, only for cases where S==1)
352 : !> \param vmatrix ...
353 : !> \param ncol ...
354 : !> \par History
355 : !> 03.2006 created [Joost VandeVondele]
356 : ! **************************************************************************************************
357 14582 : SUBROUTINE make_basis_simple(vmatrix, ncol)
358 :
359 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
360 : INTEGER, INTENT(IN) :: ncol
361 :
362 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_simple'
363 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
364 :
365 : INTEGER :: handle, n, ncol_global
366 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
367 : TYPE(cp_fm_type) :: overlap_vv
368 :
369 70 : IF (ncol == 0) RETURN
370 :
371 3628 : CALL timeset(routineN, handle)
372 :
373 3628 : NULLIFY (fm_struct_tmp)
374 :
375 3628 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
376 3628 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
377 :
378 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
379 : para_env=vmatrix%matrix_struct%para_env, &
380 3628 : context=vmatrix%matrix_struct%context)
381 3628 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
382 3628 : CALL cp_fm_struct_release(fm_struct_tmp)
383 :
384 3628 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, vmatrix, rzero, overlap_vv)
385 3628 : CALL cp_fm_cholesky_decompose(overlap_vv)
386 3628 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
387 :
388 3628 : CALL cp_fm_release(overlap_vv)
389 :
390 3628 : CALL timestop(handle)
391 :
392 3698 : END SUBROUTINE make_basis_simple
393 :
394 : ! **************************************************************************************************
395 : !> \brief computes ritz values of a set of orbitals given a ks_matrix
396 : !> rotates the orbitals into eigenstates depending on do_rotation
397 : !> writes the evals to the screen depending on ionode/scr
398 : !> \param orbitals S-orthonormal orbitals
399 : !> \param ks_matrix Kohn-Sham matrix
400 : !> \param evals_arg optional, filled with the evals
401 : !> \param ionode , scr if present write to unit scr where ionode
402 : !> \param scr ...
403 : !> \param do_rotation optional rotate orbitals (default=.TRUE.)
404 : !> note that rotating the orbitals is slower
405 : !> \param co_rotate an optional set of orbitals rotated by the same rotation matrix
406 : !> \param co_rotate_dbcsr ...
407 : !> \par History
408 : !> 08.2004 documented and added do_rotation [Joost VandeVondele]
409 : !> 09.2008 only compute eigenvalues if rotation is not needed
410 : ! **************************************************************************************************
411 8511 : SUBROUTINE subspace_eigenvalues_ks_fm(orbitals, ks_matrix, evals_arg, ionode, scr, &
412 : do_rotation, co_rotate, co_rotate_dbcsr)
413 :
414 : TYPE(cp_fm_type), INTENT(IN) :: orbitals
415 : TYPE(dbcsr_type), POINTER :: ks_matrix
416 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: evals_arg
417 : LOGICAL, INTENT(IN), OPTIONAL :: ionode
418 : INTEGER, INTENT(IN), OPTIONAL :: scr
419 : LOGICAL, INTENT(IN), OPTIONAL :: do_rotation
420 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: co_rotate
421 : TYPE(dbcsr_type), OPTIONAL, POINTER :: co_rotate_dbcsr
422 :
423 : CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_fm'
424 :
425 : INTEGER :: handle, i, j, n, ncol_global, nrow_global
426 : LOGICAL :: compute_evecs, do_rotation_local
427 8511 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
428 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
429 : TYPE(cp_fm_type) :: e_vectors, h_block, weighted_vectors, &
430 : weighted_vectors2
431 :
432 8511 : CALL timeset(routineN, handle)
433 :
434 8511 : do_rotation_local = .TRUE.
435 8511 : IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
436 :
437 8511 : NULLIFY (fm_struct_tmp)
438 : CALL cp_fm_get_info(matrix=orbitals, &
439 : ncol_global=ncol_global, &
440 8511 : nrow_global=nrow_global)
441 :
442 : IF (do_rotation_local) THEN
443 : compute_evecs = .TRUE.
444 : ELSE
445 : ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
446 : compute_evecs = .FALSE.
447 : ! this is not the case, so lets compute evecs always
448 : compute_evecs = .TRUE.
449 : END IF
450 :
451 8511 : IF (ncol_global > 0) THEN
452 :
453 25167 : ALLOCATE (evals(ncol_global))
454 :
455 8389 : CALL cp_fm_create(weighted_vectors, orbitals%matrix_struct, "weighted_vectors")
456 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
457 : para_env=orbitals%matrix_struct%para_env, &
458 8389 : context=orbitals%matrix_struct%context)
459 8389 : CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
460 8389 : IF (compute_evecs) THEN
461 8389 : CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e vectors")
462 : END IF
463 8389 : CALL cp_fm_struct_release(fm_struct_tmp)
464 :
465 : ! h subblock and diag
466 8389 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, orbitals, weighted_vectors, ncol_global)
467 :
468 : CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 1.0_dp, &
469 8389 : orbitals, weighted_vectors, 0.0_dp, h_block)
470 :
471 : ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
472 : IF (compute_evecs) THEN
473 8389 : CALL choose_eigv_solver(h_block, e_vectors, evals)
474 : ELSE
475 : CALL cp_fm_syevx(h_block, eigenvalues=evals)
476 : END IF
477 :
478 : ! rotate the orbitals
479 8389 : IF (do_rotation_local) THEN
480 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
481 8077 : orbitals, e_vectors, 0.0_dp, weighted_vectors)
482 8077 : CALL cp_fm_to_fm(weighted_vectors, orbitals)
483 8077 : IF (PRESENT(co_rotate)) THEN
484 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
485 10 : co_rotate, e_vectors, 0.0_dp, weighted_vectors)
486 10 : CALL cp_fm_to_fm(weighted_vectors, co_rotate)
487 : END IF
488 8077 : IF (PRESENT(co_rotate_dbcsr)) THEN
489 362 : IF (ASSOCIATED(co_rotate_dbcsr)) THEN
490 358 : CALL cp_fm_create(weighted_vectors2, orbitals%matrix_struct, "weighted_vectors")
491 358 : CALL copy_dbcsr_to_fm(co_rotate_dbcsr, weighted_vectors2)
492 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
493 358 : weighted_vectors2, e_vectors, 0.0_dp, weighted_vectors)
494 358 : CALL copy_fm_to_dbcsr(weighted_vectors, co_rotate_dbcsr)
495 358 : CALL cp_fm_release(weighted_vectors2)
496 : END IF
497 : END IF
498 : END IF
499 :
500 : ! give output
501 8389 : IF (PRESENT(evals_arg)) THEN
502 8357 : n = MIN(SIZE(evals_arg), SIZE(evals))
503 80199 : evals_arg(1:n) = evals(1:n)
504 : END IF
505 :
506 8389 : IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
507 182 : IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
508 182 : IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
509 182 : IF (ionode) THEN
510 228 : DO i = 1, ncol_global, 4
511 137 : j = MIN(3, ncol_global - i)
512 91 : SELECT CASE (j)
513 : CASE (3)
514 105 : WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
515 : CASE (2)
516 2 : WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
517 : CASE (1)
518 7 : WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
519 : CASE (0)
520 137 : WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
521 : END SELECT
522 : END DO
523 : END IF
524 : END IF
525 :
526 8389 : CALL cp_fm_release(weighted_vectors)
527 8389 : CALL cp_fm_release(h_block)
528 : IF (compute_evecs) THEN
529 8389 : CALL cp_fm_release(e_vectors)
530 : END IF
531 :
532 25167 : DEALLOCATE (evals)
533 :
534 : END IF
535 :
536 8511 : CALL timestop(handle)
537 :
538 17022 : END SUBROUTINE subspace_eigenvalues_ks_fm
539 :
540 : ! **************************************************************************************************
541 : !> \brief ...
542 : !> \param orbitals ...
543 : !> \param ks_matrix ...
544 : !> \param evals_arg ...
545 : !> \param ionode ...
546 : !> \param scr ...
547 : !> \param do_rotation ...
548 : !> \param co_rotate ...
549 : !> \param para_env ...
550 : !> \param blacs_env ...
551 : ! **************************************************************************************************
552 6802 : SUBROUTINE subspace_eigenvalues_ks_dbcsr(orbitals, ks_matrix, evals_arg, ionode, scr, &
553 : do_rotation, co_rotate, para_env, blacs_env)
554 :
555 : TYPE(dbcsr_type), POINTER :: orbitals, ks_matrix
556 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: evals_arg
557 : LOGICAL, INTENT(IN), OPTIONAL :: ionode
558 : INTEGER, INTENT(IN), OPTIONAL :: scr
559 : LOGICAL, INTENT(IN), OPTIONAL :: do_rotation
560 : TYPE(dbcsr_type), OPTIONAL, POINTER :: co_rotate
561 : TYPE(mp_para_env_type), POINTER :: para_env
562 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
563 :
564 : CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_dbcsr'
565 :
566 : INTEGER :: handle, i, j, ncol_global, nrow_global
567 : LOGICAL :: compute_evecs, do_rotation_local
568 6802 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
569 : TYPE(dbcsr_type), POINTER :: e_vectors, h_block, weighted_vectors
570 :
571 6802 : CALL timeset(routineN, handle)
572 :
573 6802 : do_rotation_local = .TRUE.
574 6802 : IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
575 :
576 6802 : NULLIFY (e_vectors, h_block, weighted_vectors)
577 :
578 : CALL dbcsr_get_info(matrix=orbitals, &
579 : nfullcols_total=ncol_global, &
580 6802 : nfullrows_total=nrow_global)
581 :
582 : IF (do_rotation_local) THEN
583 : compute_evecs = .TRUE.
584 : ELSE
585 : ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
586 : compute_evecs = .FALSE.
587 : ! this is not the case, so lets compute evecs always
588 : compute_evecs = .TRUE.
589 : END IF
590 :
591 6802 : IF (ncol_global > 0) THEN
592 :
593 19746 : ALLOCATE (evals(ncol_global))
594 :
595 6582 : CALL dbcsr_init_p(weighted_vectors)
596 6582 : CALL dbcsr_copy(weighted_vectors, orbitals, name="weighted_vectors")
597 :
598 6582 : CALL dbcsr_init_p(h_block)
599 : CALL cp_dbcsr_m_by_n_from_template(h_block, template=orbitals, m=ncol_global, n=ncol_global, &
600 6582 : sym=dbcsr_type_no_symmetry)
601 :
602 : !!!!!!!!!!!!!!XXXXXXXXXXXXXXXXXXX!!!!!!!!!!!!!
603 : IF (compute_evecs) THEN
604 6582 : CALL dbcsr_init_p(e_vectors)
605 : CALL cp_dbcsr_m_by_n_from_template(e_vectors, template=orbitals, m=ncol_global, n=ncol_global, &
606 6582 : sym=dbcsr_type_no_symmetry)
607 : END IF
608 :
609 : ! h subblock and diag
610 : CALL dbcsr_multiply('N', 'N', 1.0_dp, ks_matrix, orbitals, &
611 6582 : 0.0_dp, weighted_vectors)
612 : !CALL cp_dbcsr_sm_fm_multiply(ks_matrix,orbitals,weighted_vectors, ncol_global)
613 :
614 6582 : CALL dbcsr_multiply('T', 'N', 1.0_dp, orbitals, weighted_vectors, 0.0_dp, h_block)
615 : !CALL parallel_gemm('T','N',ncol_global,ncol_global,nrow_global,1.0_dp, &
616 : ! orbitals,weighted_vectors,0.0_dp,h_block)
617 :
618 : ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
619 : IF (compute_evecs) THEN
620 6582 : CALL cp_dbcsr_syevd(h_block, e_vectors, evals, para_env=para_env, blacs_env=blacs_env)
621 : ELSE
622 : CALL cp_dbcsr_syevx(h_block, eigenvalues=evals, para_env=para_env, blacs_env=blacs_env)
623 : END IF
624 :
625 : ! rotate the orbitals
626 6582 : IF (do_rotation_local) THEN
627 3222 : CALL dbcsr_multiply('N', 'N', 1.0_dp, orbitals, e_vectors, 0.0_dp, weighted_vectors)
628 : !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
629 : ! orbitals,e_vectors,0.0_dp,weighted_vectors)
630 3222 : CALL dbcsr_copy(orbitals, weighted_vectors)
631 : !CALL cp_fm_to_fm(weighted_vectors,orbitals)
632 3222 : IF (PRESENT(co_rotate)) THEN
633 3202 : IF (ASSOCIATED(co_rotate)) THEN
634 3202 : CALL dbcsr_multiply('N', 'N', 1.0_dp, co_rotate, e_vectors, 0.0_dp, weighted_vectors)
635 : !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
636 : ! co_rotate,e_vectors,0.0_dp,weighted_vectors)
637 3202 : CALL dbcsr_copy(co_rotate, weighted_vectors)
638 : !CALL cp_fm_to_fm(weighted_vectors,co_rotate)
639 : END IF
640 : END IF
641 : END IF
642 :
643 : ! give output
644 6582 : IF (PRESENT(evals_arg)) THEN
645 19944 : evals_arg(:) = evals(:)
646 : END IF
647 :
648 6582 : IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
649 0 : IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
650 0 : IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
651 0 : IF (ionode) THEN
652 0 : DO i = 1, ncol_global, 4
653 0 : j = MIN(3, ncol_global - i)
654 0 : SELECT CASE (j)
655 : CASE (3)
656 0 : WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
657 : CASE (2)
658 0 : WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
659 : CASE (1)
660 0 : WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
661 : CASE (0)
662 0 : WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
663 : END SELECT
664 : END DO
665 : END IF
666 : END IF
667 :
668 6582 : CALL dbcsr_release_p(weighted_vectors)
669 6582 : CALL dbcsr_release_p(h_block)
670 : IF (compute_evecs) THEN
671 6582 : CALL dbcsr_release_p(e_vectors)
672 : END IF
673 :
674 6582 : DEALLOCATE (evals)
675 :
676 : END IF
677 :
678 6802 : CALL timestop(handle)
679 :
680 6802 : END SUBROUTINE subspace_eigenvalues_ks_dbcsr
681 :
682 : ! computes the effective orthonormality of a set of mos given an s-matrix
683 : ! orthonormality is the max deviation from unity of the C^T S C
684 : ! **************************************************************************************************
685 : !> \brief ...
686 : !> \param orthonormality ...
687 : !> \param mo_array ...
688 : !> \param matrix_s ...
689 : ! **************************************************************************************************
690 1196 : SUBROUTINE calculate_orthonormality(orthonormality, mo_array, matrix_s)
691 : REAL(KIND=dp) :: orthonormality
692 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
693 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
694 :
695 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_orthonormality'
696 :
697 : INTEGER :: handle, i, ispin, j, k, n, ncol_local, &
698 : nrow_local, nspin
699 1196 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
700 : REAL(KIND=dp) :: alpha, max_alpha
701 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
702 : TYPE(cp_fm_type) :: overlap, svec
703 :
704 1196 : NULLIFY (tmp_fm_struct)
705 :
706 1196 : CALL timeset(routineN, handle)
707 :
708 1196 : nspin = SIZE(mo_array)
709 1196 : max_alpha = 0.0_dp
710 :
711 2404 : DO ispin = 1, nspin
712 1208 : IF (PRESENT(matrix_s)) THEN
713 : ! get S*C
714 1150 : CALL cp_fm_create(svec, mo_array(ispin)%mo_coeff%matrix_struct)
715 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
716 1150 : nrow_global=n, ncol_global=k)
717 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_array(ispin)%mo_coeff, &
718 1150 : svec, k)
719 : ! get C^T (S*C)
720 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
721 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
722 1150 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
723 1150 : CALL cp_fm_create(overlap, tmp_fm_struct)
724 1150 : CALL cp_fm_struct_release(tmp_fm_struct)
725 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
726 1150 : svec, 0.0_dp, overlap)
727 1150 : CALL cp_fm_release(svec)
728 : ELSE
729 : ! orthogonal basis C^T C
730 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
731 58 : nrow_global=n, ncol_global=k)
732 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
733 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
734 58 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
735 58 : CALL cp_fm_create(overlap, tmp_fm_struct)
736 58 : CALL cp_fm_struct_release(tmp_fm_struct)
737 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
738 58 : mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
739 : END IF
740 : CALL cp_fm_get_info(overlap, nrow_local=nrow_local, ncol_local=ncol_local, &
741 1208 : row_indices=row_indices, col_indices=col_indices)
742 8305 : DO i = 1, nrow_local
743 425678 : DO j = 1, ncol_local
744 417373 : alpha = overlap%local_data(i, j)
745 417373 : IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
746 424470 : max_alpha = MAX(max_alpha, ABS(alpha))
747 : END DO
748 : END DO
749 3612 : CALL cp_fm_release(overlap)
750 : END DO
751 1196 : CALL mo_array(1)%mo_coeff%matrix_struct%para_env%max(max_alpha)
752 1196 : orthonormality = max_alpha
753 : ! write(*,*) "max deviation from orthonormalization ",orthonormality
754 :
755 1196 : CALL timestop(handle)
756 :
757 1196 : END SUBROUTINE calculate_orthonormality
758 :
759 : ! computes the minimum/maximum magnitudes of C^T C. This could be useful
760 : ! to detect problems in the case of nearly singular overlap matrices.
761 : ! in this case, we expect the ratio of min/max to be large
762 : ! this routine is only similar to mo_orthonormality if S==1
763 : ! **************************************************************************************************
764 : !> \brief ...
765 : !> \param mo_array ...
766 : !> \param mo_mag_min ...
767 : !> \param mo_mag_max ...
768 : ! **************************************************************************************************
769 1176 : SUBROUTINE calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
770 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
771 : REAL(KIND=dp) :: mo_mag_min, mo_mag_max
772 :
773 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_magnitude'
774 :
775 : INTEGER :: handle, ispin, k, n, nspin
776 1176 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
777 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
778 : TYPE(cp_fm_type) :: evecs, overlap
779 :
780 1176 : NULLIFY (tmp_fm_struct)
781 :
782 1176 : CALL timeset(routineN, handle)
783 :
784 1176 : nspin = SIZE(mo_array)
785 1176 : mo_mag_min = HUGE(0.0_dp)
786 1176 : mo_mag_max = -HUGE(0.0_dp)
787 2352 : DO ispin = 1, nspin
788 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
789 1176 : nrow_global=n, ncol_global=k)
790 3528 : ALLOCATE (evals(k))
791 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
792 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
793 1176 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
794 1176 : CALL cp_fm_create(overlap, tmp_fm_struct)
795 1176 : CALL cp_fm_create(evecs, tmp_fm_struct)
796 1176 : CALL cp_fm_struct_release(tmp_fm_struct)
797 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
798 1176 : mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
799 1176 : CALL choose_eigv_solver(overlap, evecs, evals)
800 15238 : mo_mag_min = MIN(MINVAL(evals), mo_mag_min)
801 15238 : mo_mag_max = MAX(MAXVAL(evals), mo_mag_max)
802 1176 : CALL cp_fm_release(overlap)
803 1176 : CALL cp_fm_release(evecs)
804 4704 : DEALLOCATE (evals)
805 : END DO
806 1176 : CALL timestop(handle)
807 :
808 2352 : END SUBROUTINE calculate_magnitude
809 :
810 : ! **************************************************************************************************
811 : !> \brief Calculate KS eigenvalues starting from OF MOS
812 : !> \param mos ...
813 : !> \param nspins ...
814 : !> \param ks_rmpv ...
815 : !> \param scf_control ...
816 : !> \param mo_derivs ...
817 : !> \param admm_env ...
818 : !> \param hairy_probes ...
819 : !> \param probe ...
820 : !> \par History
821 : !> 02.2013 moved from qs_scf_post_gpw
822 : !>
823 : ! **************************************************************************************************
824 204 : SUBROUTINE make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env, &
825 : hairy_probes, probe)
826 :
827 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
828 : INTEGER, INTENT(IN) :: nspins
829 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv
830 : TYPE(scf_control_type), POINTER :: scf_control
831 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
832 : TYPE(admm_type), OPTIONAL, POINTER :: admm_env
833 : LOGICAL, OPTIONAL :: hairy_probes
834 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
835 : POINTER :: probe
836 :
837 : CHARACTER(len=*), PARAMETER :: routineN = 'make_mo_eig'
838 :
839 : INTEGER :: handle, homo, ispin, nmo, output_unit
840 204 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
841 : TYPE(cp_fm_type), POINTER :: mo_coeff
842 : TYPE(dbcsr_type), POINTER :: mo_coeff_deriv
843 :
844 204 : CALL timeset(routineN, handle)
845 :
846 204 : NULLIFY (mo_coeff_deriv, mo_coeff, mo_eigenvalues)
847 204 : output_unit = cp_logger_get_default_io_unit()
848 :
849 448 : DO ispin = 1, nspins
850 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
851 244 : eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
852 244 : IF (output_unit > 0) WRITE (output_unit, *) " "
853 244 : IF (output_unit > 0) WRITE (output_unit, *) " Eigenvalues of the occupied subspace spin ", ispin
854 : ! IF (homo /= nmo) THEN
855 : ! IF (output_unit>0) WRITE(output_unit,*)" and ",nmo-homo," added MO eigenvalues"
856 : ! END IF
857 244 : IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
858 :
859 244 : IF (scf_control%use_ot) THEN
860 : ! Also rotate the OT derivs, since they are needed for force calculations
861 128 : IF (ASSOCIATED(mo_derivs)) THEN
862 128 : mo_coeff_deriv => mo_derivs(ispin)%matrix
863 : ELSE
864 0 : mo_coeff_deriv => NULL()
865 : END IF
866 :
867 : ! If we do ADMM, we add have to modify the Kohn-Sham matrix
868 128 : IF (PRESENT(admm_env)) THEN
869 0 : CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
870 : END IF
871 :
872 : CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, &
873 : scr=output_unit, ionode=output_unit > 0, do_rotation=.TRUE., &
874 128 : co_rotate_dbcsr=mo_coeff_deriv)
875 :
876 : ! If we do ADMM, we restore the original Kohn-Sham matrix
877 128 : IF (PRESENT(admm_env)) THEN
878 0 : CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
879 : END IF
880 : ELSE
881 567 : IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo)
882 : END IF
883 244 : IF (.NOT. scf_control%diagonalization%mom) THEN
884 244 : IF (PRESENT(hairy_probes) .EQV. .TRUE.) THEN
885 0 : scf_control%smear%do_smear = .FALSE.
886 : CALL set_mo_occupation(mo_set=mos(ispin), &
887 : smear=scf_control%smear, &
888 0 : probe=probe)
889 : ELSE
890 244 : CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear)
891 : END IF
892 : END IF
893 244 : IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
894 570 : "Fermi Energy [eV] :", mos(ispin)%mu*evolt
895 : END DO
896 :
897 204 : CALL timestop(handle)
898 :
899 204 : END SUBROUTINE make_mo_eig
900 :
901 : ! **************************************************************************************************
902 : !> \brief Orbital energies <phi_i|F|phi_i> of a set of orbitals and the largest off-diagonal Fock
903 : !> matrix element max_{i/=j} |<phi_i|F|phi_j>| within the set. Unlike
904 : !> calculate_subspace_eigenvalues, the orbitals are not rotated and the energies belong to
905 : !> the given orbitals. The orbitals are canonical if the off-diagonal elements vanish.
906 : !> \param orbitals S-orthonormal orbitals (columns)
907 : !> \param ks_matrix Kohn-Sham or Fock matrix
908 : !> \param orbital_energies diagonal elements; filled up to MIN(SIZE, number of orbitals)
909 : !> \param max_offdiag largest off-diagonal element in absolute value
910 : ! **************************************************************************************************
911 252 : SUBROUTINE calculate_orbital_fock_elements(orbitals, ks_matrix, orbital_energies, max_offdiag)
912 : TYPE(cp_fm_type), INTENT(IN) :: orbitals
913 : TYPE(dbcsr_type), POINTER :: ks_matrix
914 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: orbital_energies
915 : REAL(KIND=dp), INTENT(OUT) :: max_offdiag
916 :
917 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_orbital_fock_elements'
918 :
919 : INTEGER :: handle, i, ncol_global, nrow_global
920 252 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fock
921 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
922 : TYPE(cp_fm_type) :: h_block, weighted_vectors
923 :
924 252 : CALL timeset(routineN, handle)
925 :
926 1496 : orbital_energies = 0.0_dp
927 252 : max_offdiag = 0.0_dp
928 252 : CALL cp_fm_get_info(matrix=orbitals, ncol_global=ncol_global, nrow_global=nrow_global)
929 252 : IF (ncol_global > 0) THEN
930 252 : NULLIFY (fm_struct_tmp)
931 252 : CALL cp_fm_create(weighted_vectors, orbitals%matrix_struct, "weighted_vectors")
932 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
933 : para_env=orbitals%matrix_struct%para_env, &
934 252 : context=orbitals%matrix_struct%context)
935 252 : CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
936 252 : CALL cp_fm_struct_release(fm_struct_tmp)
937 :
938 252 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, orbitals, weighted_vectors, ncol_global)
939 : CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 1.0_dp, &
940 252 : orbitals, weighted_vectors, 0.0_dp, h_block)
941 :
942 1008 : ALLOCATE (fock(ncol_global, ncol_global))
943 252 : CALL cp_fm_get_submatrix(h_block, fock)
944 1480 : DO i = 1, MIN(ncol_global, SIZE(orbital_energies))
945 1480 : orbital_energies(i) = fock(i, i)
946 : END DO
947 1824 : DO i = 1, ncol_global
948 1824 : fock(i, i) = 0.0_dp
949 : END DO
950 16892 : max_offdiag = MAXVAL(ABS(fock))
951 252 : DEALLOCATE (fock)
952 :
953 252 : CALL cp_fm_release(weighted_vectors)
954 756 : CALL cp_fm_release(h_block)
955 : END IF
956 :
957 252 : CALL timestop(handle)
958 :
959 504 : END SUBROUTINE calculate_orbital_fock_elements
960 :
961 : ! **************************************************************************************************
962 : !> \brief Rotates the orbitals of an MO set into canonical orbitals (eigenvectors of the Fock
963 : !> matrix in the space of the orbitals), separately within the occupied and within the
964 : !> unoccupied orbitals, so that the occupied space (the determinant) is unchanged.
965 : !> The orbitals are rotated in place: the eigenvalues, the DBCSR copy of the coefficients
966 : !> (if used) and the orbital derivatives (if given) are rotated consistently, C' = C U and
967 : !> D' = D U. Nothing is done if the orbitals are already canonical, i.e. if no
968 : !> off-diagonal Fock matrix element exceeds eps_canonical; orbitals of a converged
969 : !> diagonalization SCF are therefore never changed.
970 : !> Requires integer occupations (0 or the maximum occupation); fractionally occupied
971 : !> orbitals and ROKS orbitals must not be mixed this way.
972 : !> \param mo_set the MO set (e.g. the SCF orbitals after an OT optimization)
973 : !> \param ks_matrix Kohn-Sham or Fock matrix
974 : !> \param eps_canonical largest off-diagonal Fock matrix element accepted as canonical
975 : !> \param rotated whether the orbitals were rotated
976 : !> \param max_offdiag largest off-diagonal Fock matrix element before the rotation
977 : !> \param mo_derivs derivatives of the energy with respect to the orbitals (OT), rotated as well
978 : ! **************************************************************************************************
979 10 : SUBROUTINE make_mo_set_canonical(mo_set, ks_matrix, eps_canonical, rotated, max_offdiag, mo_derivs)
980 : TYPE(mo_set_type), INTENT(INOUT) :: mo_set
981 : TYPE(dbcsr_type), POINTER :: ks_matrix
982 : REAL(KIND=dp), INTENT(IN) :: eps_canonical
983 : LOGICAL, INTENT(OUT) :: rotated
984 : REAL(KIND=dp), INTENT(OUT) :: max_offdiag
985 : TYPE(dbcsr_type), OPTIONAL, POINTER :: mo_derivs
986 :
987 : CHARACTER(len=*), PARAMETER :: routineN = 'make_mo_set_canonical'
988 : REAL(KIND=dp), PARAMETER :: eps_occ = 1.0E-10_dp
989 :
990 : INTEGER :: first, handle, i, iblock, n, nmo, nocc, &
991 : nrow_global
992 : LOGICAL :: do_derivs
993 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: energies
994 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
995 : TYPE(cp_fm_type) :: block_coeff, block_derivs, derivs
996 :
997 10 : CALL timeset(routineN, handle)
998 :
999 10 : rotated = .FALSE.
1000 10 : nmo = mo_set%nmo
1001 10 : IF (mo_set%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mo_set%mo_coeff_b, mo_set%mo_coeff)
1002 30 : ALLOCATE (energies(nmo))
1003 10 : CALL calculate_orbital_fock_elements(mo_set%mo_coeff, ks_matrix, energies, max_offdiag)
1004 10 : DEALLOCATE (energies)
1005 :
1006 10 : IF (max_offdiag > eps_canonical) THEN
1007 50 : DO i = 1, nmo
1008 40 : IF (ABS(mo_set%occupation_numbers(i)) > eps_occ .AND. &
1009 10 : ABS(mo_set%occupation_numbers(i) - mo_set%maxocc) > eps_occ) THEN
1010 0 : CPABORT("Canonicalization of orbitals with fractional occupations is not possible")
1011 : END IF
1012 : END DO
1013 10 : do_derivs = .FALSE.
1014 10 : IF (PRESENT(mo_derivs)) do_derivs = ASSOCIATED(mo_derivs)
1015 : ! occupied orbitals: those with nonzero occupation, they come first
1016 50 : nocc = COUNT(mo_set%occupation_numbers(1:nmo) > eps_occ)
1017 10 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, nrow_global=nrow_global)
1018 10 : IF (do_derivs) THEN
1019 10 : CALL cp_fm_create(derivs, mo_set%mo_coeff%matrix_struct, name="orbital derivatives")
1020 10 : CALL copy_dbcsr_to_fm(mo_derivs, derivs)
1021 : END IF
1022 30 : DO iblock = 1, 2
1023 20 : IF (iblock == 1) THEN
1024 10 : first = 1
1025 10 : n = nocc
1026 : ELSE
1027 10 : first = nocc + 1
1028 10 : n = nmo - nocc
1029 : END IF
1030 20 : IF (n <= 0) CYCLE
1031 10 : NULLIFY (fm_struct_tmp)
1032 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, ncol_global=n, &
1033 10 : template_fmstruct=mo_set%mo_coeff%matrix_struct)
1034 10 : CALL cp_fm_create(block_coeff, fm_struct_tmp, name="canonicalized orbitals")
1035 10 : CALL cp_fm_to_fm(mo_set%mo_coeff, block_coeff, n, first, 1)
1036 10 : IF (do_derivs) THEN
1037 10 : CALL cp_fm_create(block_derivs, fm_struct_tmp, name="rotated orbital derivatives")
1038 10 : CALL cp_fm_to_fm(derivs, block_derivs, n, first, 1)
1039 : CALL calculate_subspace_eigenvalues(block_coeff, ks_matrix, &
1040 : evals_arg=mo_set%eigenvalues(first:first + n - 1), &
1041 10 : do_rotation=.TRUE., co_rotate=block_derivs)
1042 10 : CALL cp_fm_to_fm(block_derivs, derivs, n, 1, first)
1043 10 : CALL cp_fm_release(block_derivs)
1044 : ELSE
1045 : CALL calculate_subspace_eigenvalues(block_coeff, ks_matrix, &
1046 : evals_arg=mo_set%eigenvalues(first:first + n - 1), &
1047 0 : do_rotation=.TRUE.)
1048 : END IF
1049 10 : CALL cp_fm_struct_release(fm_struct_tmp)
1050 10 : CALL cp_fm_to_fm(block_coeff, mo_set%mo_coeff, n, 1, first)
1051 40 : CALL cp_fm_release(block_coeff)
1052 : END DO
1053 10 : IF (do_derivs) THEN
1054 10 : CALL copy_fm_to_dbcsr(derivs, mo_derivs)
1055 10 : CALL cp_fm_release(derivs)
1056 : END IF
1057 10 : IF (mo_set%use_mo_coeff_b) CALL copy_fm_to_dbcsr(mo_set%mo_coeff, mo_set%mo_coeff_b)
1058 10 : rotated = .TRUE.
1059 : END IF
1060 :
1061 10 : CALL timestop(handle)
1062 :
1063 20 : END SUBROUTINE make_mo_set_canonical
1064 :
1065 : END MODULE qs_mo_methods
|
