Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_release,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
49 : USE kinds, ONLY: dp
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE parallel_gemm_api, ONLY: parallel_gemm
52 : USE physcon, ONLY: evolt
53 : USE qs_mo_occupation, ONLY: set_mo_occupation
54 : USE qs_mo_types, ONLY: get_mo_set,&
55 : mo_set_type
56 : USE scf_control_types, ONLY: scf_control_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 : PRIVATE
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_methods'
62 :
63 : PUBLIC :: make_basis_simple, make_basis_cholesky, make_basis_sv, make_basis_sm, &
64 : make_basis_lowdin, calculate_subspace_eigenvalues, &
65 : calculate_orthonormality, calculate_magnitude, make_mo_eig
66 :
67 : INTERFACE calculate_subspace_eigenvalues
68 : MODULE PROCEDURE subspace_eigenvalues_ks_fm
69 : MODULE PROCEDURE subspace_eigenvalues_ks_dbcsr
70 : END INTERFACE
71 :
72 : INTERFACE make_basis_sv
73 : MODULE PROCEDURE make_basis_sv_fm
74 : MODULE PROCEDURE make_basis_sv_dbcsr
75 : END INTERFACE
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief returns an S-orthonormal basis v (v^T S v ==1)
81 : !> \param vmatrix ...
82 : !> \param ncol ...
83 : !> \param matrix_s ...
84 : !> \par History
85 : !> 03.2006 created [Joost VandeVondele]
86 : ! **************************************************************************************************
87 27376 : SUBROUTINE make_basis_sm(vmatrix, ncol, matrix_s)
88 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
89 : INTEGER, INTENT(IN) :: ncol
90 : TYPE(dbcsr_type), POINTER :: matrix_s
91 :
92 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sm'
93 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
94 :
95 : INTEGER :: handle, n, ncol_global
96 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
97 : TYPE(cp_fm_type) :: overlap_vv, svmatrix
98 :
99 128 : IF (ncol == 0) RETURN
100 :
101 6812 : CALL timeset(routineN, handle)
102 :
103 6812 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
104 6812 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
105 :
106 6812 : CALL cp_fm_create(svmatrix, vmatrix%matrix_struct, "SV", set_zero=.TRUE.)
107 6812 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, svmatrix, ncol)
108 :
109 6812 : NULLIFY (fm_struct_tmp)
110 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
111 : para_env=vmatrix%matrix_struct%para_env, &
112 6812 : context=vmatrix%matrix_struct%context)
113 6812 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
114 6812 : CALL cp_fm_struct_release(fm_struct_tmp)
115 :
116 6812 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
117 6812 : CALL cp_fm_cholesky_decompose(overlap_vv)
118 6812 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
119 :
120 6812 : CALL cp_fm_release(overlap_vv)
121 6812 : CALL cp_fm_release(svmatrix)
122 :
123 6812 : CALL timestop(handle)
124 :
125 6940 : END SUBROUTINE make_basis_sm
126 :
127 : ! **************************************************************************************************
128 : !> \brief returns an S-orthonormal basis v and the corresponding matrix S*v as well
129 : !> \param vmatrix ...
130 : !> \param ncol ...
131 : !> \param svmatrix ...
132 : !> \par History
133 : !> 03.2006 created [Joost VandeVondele]
134 : ! **************************************************************************************************
135 0 : SUBROUTINE make_basis_sv_fm(vmatrix, ncol, svmatrix)
136 :
137 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
138 : INTEGER, INTENT(IN) :: ncol
139 : TYPE(cp_fm_type), INTENT(IN) :: svmatrix
140 :
141 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sv_fm'
142 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
143 :
144 : INTEGER :: handle, n, ncol_global
145 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
146 : TYPE(cp_fm_type) :: overlap_vv
147 :
148 0 : IF (ncol == 0) RETURN
149 :
150 0 : CALL timeset(routineN, handle)
151 0 : NULLIFY (fm_struct_tmp)
152 :
153 0 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
154 0 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
155 :
156 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
157 : para_env=vmatrix%matrix_struct%para_env, &
158 0 : context=vmatrix%matrix_struct%context)
159 0 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
160 0 : CALL cp_fm_struct_release(fm_struct_tmp)
161 :
162 0 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
163 0 : CALL cp_fm_cholesky_decompose(overlap_vv)
164 0 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
165 0 : CALL cp_fm_triangular_multiply(overlap_vv, svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
166 :
167 0 : CALL cp_fm_release(overlap_vv)
168 :
169 0 : CALL timestop(handle)
170 :
171 0 : END SUBROUTINE make_basis_sv_fm
172 :
173 : ! **************************************************************************************************
174 : !> \brief ...
175 : !> \param vmatrix ...
176 : !> \param ncol ...
177 : !> \param svmatrix ...
178 : !> \param para_env ...
179 : !> \param blacs_env ...
180 : ! **************************************************************************************************
181 2838 : SUBROUTINE make_basis_sv_dbcsr(vmatrix, ncol, svmatrix, para_env, blacs_env)
182 :
183 : TYPE(dbcsr_type) :: vmatrix
184 : INTEGER, INTENT(IN) :: ncol
185 : TYPE(dbcsr_type) :: svmatrix
186 : TYPE(mp_para_env_type), POINTER :: para_env
187 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
188 :
189 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_sv_dbcsr'
190 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
191 :
192 : INTEGER :: handle, n, ncol_global
193 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
194 : TYPE(cp_fm_type) :: fm_svmatrix, fm_vmatrix, overlap_vv
195 :
196 78 : IF (ncol == 0) RETURN
197 :
198 460 : CALL timeset(routineN, handle)
199 :
200 : !CALL cp_fm_get_info(matrix=vmatrix,nrow_global=n,ncol_global=ncol_global)
201 460 : CALL dbcsr_get_info(vmatrix, nfullrows_total=n, nfullcols_total=ncol_global)
202 460 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
203 :
204 : CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=ncol, &
205 460 : ncol_global=ncol, para_env=para_env)
206 460 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, name="fm_overlap_vv")
207 460 : CALL cp_fm_struct_release(fm_struct_tmp)
208 :
209 : CALL cp_fm_struct_create(fm_struct_tmp, context=blacs_env, nrow_global=n, &
210 460 : ncol_global=ncol_global, para_env=para_env)
211 460 : CALL cp_fm_create(fm_vmatrix, fm_struct_tmp, name="fm_vmatrix")
212 460 : CALL cp_fm_create(fm_svmatrix, fm_struct_tmp, name="fm_svmatrix")
213 460 : CALL cp_fm_struct_release(fm_struct_tmp)
214 :
215 460 : CALL copy_dbcsr_to_fm(vmatrix, fm_vmatrix)
216 460 : CALL copy_dbcsr_to_fm(svmatrix, fm_svmatrix)
217 :
218 460 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, fm_vmatrix, fm_svmatrix, rzero, overlap_vv)
219 460 : CALL cp_fm_cholesky_decompose(overlap_vv)
220 460 : CALL cp_fm_triangular_multiply(overlap_vv, fm_vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
221 460 : CALL cp_fm_triangular_multiply(overlap_vv, fm_svmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
222 :
223 460 : CALL copy_fm_to_dbcsr(fm_vmatrix, vmatrix)
224 460 : CALL copy_fm_to_dbcsr(fm_svmatrix, svmatrix)
225 :
226 460 : CALL cp_fm_release(overlap_vv)
227 460 : CALL cp_fm_release(fm_vmatrix)
228 460 : CALL cp_fm_release(fm_svmatrix)
229 :
230 460 : CALL timestop(handle)
231 :
232 538 : END SUBROUTINE make_basis_sv_dbcsr
233 :
234 : ! **************************************************************************************************
235 : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
236 : !> the cholesky decomposed form of S is passed as an argument
237 : !> \param vmatrix ...
238 : !> \param ncol ...
239 : !> \param ortho cholesky decomposed S matrix
240 : !> \par History
241 : !> 03.2006 created [Joost VandeVondele]
242 : !> \note
243 : !> if the cholesky decomposed S matrix is not available
244 : !> use make_basis_sm since this is much faster than computing the
245 : !> cholesky decomposition of S
246 : ! **************************************************************************************************
247 24632 : SUBROUTINE make_basis_cholesky(vmatrix, ncol, ortho)
248 :
249 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
250 : INTEGER, INTENT(IN) :: ncol
251 : TYPE(cp_fm_type), INTENT(IN) :: ortho
252 :
253 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_cholesky'
254 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
255 :
256 : INTEGER :: handle, n, ncol_global
257 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
258 : TYPE(cp_fm_type) :: overlap_vv
259 :
260 80 : IF (ncol == 0) RETURN
261 :
262 6138 : CALL timeset(routineN, handle)
263 6138 : NULLIFY (fm_struct_tmp)
264 :
265 6138 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
266 6138 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
267 :
268 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
269 : para_env=vmatrix%matrix_struct%para_env, &
270 6138 : context=vmatrix%matrix_struct%context)
271 6138 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
272 6138 : CALL cp_fm_struct_release(fm_struct_tmp)
273 :
274 6138 : CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol)
275 6138 : CALL cp_fm_syrk('U', 'T', n, rone, vmatrix, 1, 1, rzero, overlap_vv)
276 6138 : CALL cp_fm_cholesky_decompose(overlap_vv)
277 6138 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
278 6138 : CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol, invert_tr=.TRUE.)
279 :
280 6138 : CALL cp_fm_release(overlap_vv)
281 :
282 6138 : CALL timestop(handle)
283 :
284 6218 : END SUBROUTINE make_basis_cholesky
285 :
286 : ! **************************************************************************************************
287 : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
288 : !> a Loedwin transformation is applied to keep the rotated vectors as close
289 : !> as possible to the original ones
290 : !> \param vmatrix ...
291 : !> \param ncol ...
292 : !> \param matrix_s ...
293 : !> \param
294 : !> \par History
295 : !> 05.2009 created [MI]
296 : !> \note
297 : ! **************************************************************************************************
298 2868 : SUBROUTINE make_basis_lowdin(vmatrix, ncol, matrix_s)
299 :
300 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
301 : INTEGER, INTENT(IN) :: ncol
302 : TYPE(dbcsr_type) :: matrix_s
303 :
304 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_lowdin'
305 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
306 :
307 : INTEGER :: handle, n, ncol_global, ndep
308 : REAL(dp) :: threshold
309 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
310 : TYPE(cp_fm_type) :: csc, sc, work
311 :
312 0 : IF (ncol == 0) RETURN
313 :
314 478 : CALL timeset(routineN, handle)
315 478 : NULLIFY (fm_struct_tmp)
316 478 : threshold = 1.0E-7_dp
317 478 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
318 478 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
319 :
320 478 : CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
321 478 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
322 :
323 478 : NULLIFY (fm_struct_tmp)
324 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
325 : para_env=vmatrix%matrix_struct%para_env, &
326 478 : context=vmatrix%matrix_struct%context)
327 478 : CALL cp_fm_create(csc, fm_struct_tmp, "csc")
328 478 : CALL cp_fm_create(work, fm_struct_tmp, "work")
329 478 : CALL cp_fm_struct_release(fm_struct_tmp)
330 :
331 478 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
332 478 : CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
333 478 : CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
334 478 : CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
335 :
336 478 : CALL cp_fm_release(csc)
337 478 : CALL cp_fm_release(sc)
338 478 : CALL cp_fm_release(work)
339 :
340 478 : CALL timestop(handle)
341 :
342 478 : END SUBROUTINE make_basis_lowdin
343 :
344 : ! **************************************************************************************************
345 : !> \brief given a set of vectors, return an orthogonal (C^T C == 1) set
346 : !> spanning the same space (notice, only for cases where S==1)
347 : !> \param vmatrix ...
348 : !> \param ncol ...
349 : !> \par History
350 : !> 03.2006 created [Joost VandeVondele]
351 : ! **************************************************************************************************
352 14526 : SUBROUTINE make_basis_simple(vmatrix, ncol)
353 :
354 : TYPE(cp_fm_type), INTENT(IN) :: vmatrix
355 : INTEGER, INTENT(IN) :: ncol
356 :
357 : CHARACTER(LEN=*), PARAMETER :: routineN = 'make_basis_simple'
358 : REAL(KIND=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
359 :
360 : INTEGER :: handle, n, ncol_global
361 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
362 : TYPE(cp_fm_type) :: overlap_vv
363 :
364 70 : IF (ncol == 0) RETURN
365 :
366 3614 : CALL timeset(routineN, handle)
367 :
368 3614 : NULLIFY (fm_struct_tmp)
369 :
370 3614 : CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
371 3614 : IF (ncol > ncol_global) CPABORT("Wrong ncol value")
372 :
373 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
374 : para_env=vmatrix%matrix_struct%para_env, &
375 3614 : context=vmatrix%matrix_struct%context)
376 3614 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
377 3614 : CALL cp_fm_struct_release(fm_struct_tmp)
378 :
379 3614 : CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, vmatrix, rzero, overlap_vv)
380 3614 : CALL cp_fm_cholesky_decompose(overlap_vv)
381 3614 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
382 :
383 3614 : CALL cp_fm_release(overlap_vv)
384 :
385 3614 : CALL timestop(handle)
386 :
387 3684 : END SUBROUTINE make_basis_simple
388 :
389 : ! **************************************************************************************************
390 : !> \brief computes ritz values of a set of orbitals given a ks_matrix
391 : !> rotates the orbitals into eigenstates depending on do_rotation
392 : !> writes the evals to the screen depending on ionode/scr
393 : !> \param orbitals S-orthonormal orbitals
394 : !> \param ks_matrix Kohn-Sham matrix
395 : !> \param evals_arg optional, filled with the evals
396 : !> \param ionode , scr if present write to unit scr where ionode
397 : !> \param scr ...
398 : !> \param do_rotation optional rotate orbitals (default=.TRUE.)
399 : !> note that rotating the orbitals is slower
400 : !> \param co_rotate an optional set of orbitals rotated by the same rotation matrix
401 : !> \param co_rotate_dbcsr ...
402 : !> \par History
403 : !> 08.2004 documented and added do_rotation [Joost VandeVondele]
404 : !> 09.2008 only compute eigenvalues if rotation is not needed
405 : ! **************************************************************************************************
406 7515 : SUBROUTINE subspace_eigenvalues_ks_fm(orbitals, ks_matrix, evals_arg, ionode, scr, &
407 : do_rotation, co_rotate, co_rotate_dbcsr)
408 :
409 : TYPE(cp_fm_type), INTENT(IN) :: orbitals
410 : TYPE(dbcsr_type), POINTER :: ks_matrix
411 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: evals_arg
412 : LOGICAL, INTENT(IN), OPTIONAL :: ionode
413 : INTEGER, INTENT(IN), OPTIONAL :: scr
414 : LOGICAL, INTENT(IN), OPTIONAL :: do_rotation
415 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: co_rotate
416 : TYPE(dbcsr_type), OPTIONAL, POINTER :: co_rotate_dbcsr
417 :
418 : CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_fm'
419 :
420 : INTEGER :: handle, i, j, n, ncol_global, nrow_global
421 : LOGICAL :: compute_evecs, do_rotation_local
422 7515 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
423 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
424 : TYPE(cp_fm_type) :: e_vectors, h_block, weighted_vectors, &
425 : weighted_vectors2
426 :
427 7515 : CALL timeset(routineN, handle)
428 :
429 7515 : do_rotation_local = .TRUE.
430 7515 : IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
431 :
432 7515 : NULLIFY (fm_struct_tmp)
433 : CALL cp_fm_get_info(matrix=orbitals, &
434 : ncol_global=ncol_global, &
435 7515 : nrow_global=nrow_global)
436 :
437 : IF (do_rotation_local) THEN
438 : compute_evecs = .TRUE.
439 : ELSE
440 : ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
441 : compute_evecs = .FALSE.
442 : ! this is not the case, so lets compute evecs always
443 : compute_evecs = .TRUE.
444 : END IF
445 :
446 7515 : IF (ncol_global > 0) THEN
447 :
448 22311 : ALLOCATE (evals(ncol_global))
449 :
450 7437 : CALL cp_fm_create(weighted_vectors, orbitals%matrix_struct, "weighted_vectors")
451 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
452 : para_env=orbitals%matrix_struct%para_env, &
453 7437 : context=orbitals%matrix_struct%context)
454 7437 : CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
455 7437 : IF (compute_evecs) THEN
456 7437 : CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e vectors")
457 : END IF
458 7437 : CALL cp_fm_struct_release(fm_struct_tmp)
459 :
460 : ! h subblock and diag
461 7437 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, orbitals, weighted_vectors, ncol_global)
462 :
463 : CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 1.0_dp, &
464 7437 : orbitals, weighted_vectors, 0.0_dp, h_block)
465 :
466 : ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
467 : IF (compute_evecs) THEN
468 7437 : CALL choose_eigv_solver(h_block, e_vectors, evals)
469 : ELSE
470 : CALL cp_fm_syevx(h_block, eigenvalues=evals)
471 : END IF
472 :
473 : ! rotate the orbitals
474 7437 : IF (do_rotation_local) THEN
475 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
476 7153 : orbitals, e_vectors, 0.0_dp, weighted_vectors)
477 7153 : CALL cp_fm_to_fm(weighted_vectors, orbitals)
478 7153 : IF (PRESENT(co_rotate)) THEN
479 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
480 0 : co_rotate, e_vectors, 0.0_dp, weighted_vectors)
481 0 : CALL cp_fm_to_fm(weighted_vectors, co_rotate)
482 : END IF
483 7153 : IF (PRESENT(co_rotate_dbcsr)) THEN
484 136 : IF (ASSOCIATED(co_rotate_dbcsr)) THEN
485 136 : CALL cp_fm_create(weighted_vectors2, orbitals%matrix_struct, "weighted_vectors")
486 136 : CALL copy_dbcsr_to_fm(co_rotate_dbcsr, weighted_vectors2)
487 : CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
488 136 : weighted_vectors2, e_vectors, 0.0_dp, weighted_vectors)
489 136 : CALL copy_fm_to_dbcsr(weighted_vectors, co_rotate_dbcsr)
490 136 : CALL cp_fm_release(weighted_vectors2)
491 : END IF
492 : END IF
493 : END IF
494 :
495 : ! give output
496 7437 : IF (PRESENT(evals_arg)) THEN
497 7405 : n = MIN(SIZE(evals_arg), SIZE(evals))
498 77057 : evals_arg(1:n) = evals(1:n)
499 : END IF
500 :
501 7437 : IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
502 180 : IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
503 180 : IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
504 180 : IF (ionode) THEN
505 225 : DO i = 1, ncol_global, 4
506 135 : j = MIN(3, ncol_global - i)
507 90 : SELECT CASE (j)
508 : CASE (3)
509 104 : WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
510 : CASE (2)
511 2 : WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
512 : CASE (1)
513 7 : WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
514 : CASE (0)
515 135 : WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
516 : END SELECT
517 : END DO
518 : END IF
519 : END IF
520 :
521 7437 : CALL cp_fm_release(weighted_vectors)
522 7437 : CALL cp_fm_release(h_block)
523 : IF (compute_evecs) THEN
524 7437 : CALL cp_fm_release(e_vectors)
525 : END IF
526 :
527 22311 : DEALLOCATE (evals)
528 :
529 : END IF
530 :
531 7515 : CALL timestop(handle)
532 :
533 15030 : END SUBROUTINE subspace_eigenvalues_ks_fm
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param orbitals ...
538 : !> \param ks_matrix ...
539 : !> \param evals_arg ...
540 : !> \param ionode ...
541 : !> \param scr ...
542 : !> \param do_rotation ...
543 : !> \param co_rotate ...
544 : !> \param para_env ...
545 : !> \param blacs_env ...
546 : ! **************************************************************************************************
547 5628 : SUBROUTINE subspace_eigenvalues_ks_dbcsr(orbitals, ks_matrix, evals_arg, ionode, scr, &
548 : do_rotation, co_rotate, para_env, blacs_env)
549 :
550 : TYPE(dbcsr_type), POINTER :: orbitals, ks_matrix
551 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: evals_arg
552 : LOGICAL, INTENT(IN), OPTIONAL :: ionode
553 : INTEGER, INTENT(IN), OPTIONAL :: scr
554 : LOGICAL, INTENT(IN), OPTIONAL :: do_rotation
555 : TYPE(dbcsr_type), OPTIONAL, POINTER :: co_rotate
556 : TYPE(mp_para_env_type), POINTER :: para_env
557 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
558 :
559 : CHARACTER(len=*), PARAMETER :: routineN = 'subspace_eigenvalues_ks_dbcsr'
560 :
561 : INTEGER :: handle, i, j, ncol_global, nrow_global
562 : LOGICAL :: compute_evecs, do_rotation_local
563 5628 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
564 : TYPE(dbcsr_type), POINTER :: e_vectors, h_block, weighted_vectors
565 :
566 5628 : CALL timeset(routineN, handle)
567 :
568 5628 : do_rotation_local = .TRUE.
569 5628 : IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
570 :
571 5628 : NULLIFY (e_vectors, h_block, weighted_vectors)
572 :
573 : CALL dbcsr_get_info(matrix=orbitals, &
574 : nfullcols_total=ncol_global, &
575 5628 : nfullrows_total=nrow_global)
576 :
577 : IF (do_rotation_local) THEN
578 : compute_evecs = .TRUE.
579 : ELSE
580 : ! this would be the logical choice if syevx computing only evals were faster than syevd computing evecs and evals.
581 : compute_evecs = .FALSE.
582 : ! this is not the case, so lets compute evecs always
583 : compute_evecs = .TRUE.
584 : END IF
585 :
586 5628 : IF (ncol_global > 0) THEN
587 :
588 16248 : ALLOCATE (evals(ncol_global))
589 :
590 5416 : CALL dbcsr_init_p(weighted_vectors)
591 5416 : CALL dbcsr_copy(weighted_vectors, orbitals, name="weighted_vectors")
592 :
593 5416 : CALL dbcsr_init_p(h_block)
594 : CALL cp_dbcsr_m_by_n_from_template(h_block, template=orbitals, m=ncol_global, n=ncol_global, &
595 5416 : sym=dbcsr_type_no_symmetry)
596 :
597 : !!!!!!!!!!!!!!XXXXXXXXXXXXXXXXXXX!!!!!!!!!!!!!
598 : IF (compute_evecs) THEN
599 5416 : CALL dbcsr_init_p(e_vectors)
600 : CALL cp_dbcsr_m_by_n_from_template(e_vectors, template=orbitals, m=ncol_global, n=ncol_global, &
601 5416 : sym=dbcsr_type_no_symmetry)
602 : END IF
603 :
604 : ! h subblock and diag
605 : CALL dbcsr_multiply('N', 'N', 1.0_dp, ks_matrix, orbitals, &
606 5416 : 0.0_dp, weighted_vectors)
607 : !CALL cp_dbcsr_sm_fm_multiply(ks_matrix,orbitals,weighted_vectors, ncol_global)
608 :
609 5416 : CALL dbcsr_multiply('T', 'N', 1.0_dp, orbitals, weighted_vectors, 0.0_dp, h_block)
610 : !CALL parallel_gemm('T','N',ncol_global,ncol_global,nrow_global,1.0_dp, &
611 : ! orbitals,weighted_vectors,0.0_dp,h_block)
612 :
613 : ! if eigenvectors are required, go for syevd, otherwise only compute eigenvalues
614 : IF (compute_evecs) THEN
615 5416 : CALL cp_dbcsr_syevd(h_block, e_vectors, evals, para_env=para_env, blacs_env=blacs_env)
616 : ELSE
617 : CALL cp_dbcsr_syevx(h_block, eigenvalues=evals, para_env=para_env, blacs_env=blacs_env)
618 : END IF
619 :
620 : ! rotate the orbitals
621 5416 : IF (do_rotation_local) THEN
622 2654 : CALL dbcsr_multiply('N', 'N', 1.0_dp, orbitals, e_vectors, 0.0_dp, weighted_vectors)
623 : !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
624 : ! orbitals,e_vectors,0.0_dp,weighted_vectors)
625 2654 : CALL dbcsr_copy(orbitals, weighted_vectors)
626 : !CALL cp_fm_to_fm(weighted_vectors,orbitals)
627 2654 : IF (PRESENT(co_rotate)) THEN
628 2636 : IF (ASSOCIATED(co_rotate)) THEN
629 2636 : CALL dbcsr_multiply('N', 'N', 1.0_dp, co_rotate, e_vectors, 0.0_dp, weighted_vectors)
630 : !CALL parallel_gemm('N','N',nrow_global,ncol_global,ncol_global,1.0_dp, &
631 : ! co_rotate,e_vectors,0.0_dp,weighted_vectors)
632 2636 : CALL dbcsr_copy(co_rotate, weighted_vectors)
633 : !CALL cp_fm_to_fm(weighted_vectors,co_rotate)
634 : END IF
635 : END IF
636 : END IF
637 :
638 : ! give output
639 5416 : IF (PRESENT(evals_arg)) THEN
640 16440 : evals_arg(:) = evals(:)
641 : END IF
642 :
643 5416 : IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
644 0 : IF (.NOT. PRESENT(ionode)) CPABORT("IONODE?")
645 0 : IF (.NOT. PRESENT(scr)) CPABORT("SCR?")
646 0 : IF (ionode) THEN
647 0 : DO i = 1, ncol_global, 4
648 0 : j = MIN(3, ncol_global - i)
649 0 : SELECT CASE (j)
650 : CASE (3)
651 0 : WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
652 : CASE (2)
653 0 : WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
654 : CASE (1)
655 0 : WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
656 : CASE (0)
657 0 : WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
658 : END SELECT
659 : END DO
660 : END IF
661 : END IF
662 :
663 5416 : CALL dbcsr_release_p(weighted_vectors)
664 5416 : CALL dbcsr_release_p(h_block)
665 : IF (compute_evecs) THEN
666 5416 : CALL dbcsr_release_p(e_vectors)
667 : END IF
668 :
669 5416 : DEALLOCATE (evals)
670 :
671 : END IF
672 :
673 5628 : CALL timestop(handle)
674 :
675 5628 : END SUBROUTINE subspace_eigenvalues_ks_dbcsr
676 :
677 : ! computes the effective orthonormality of a set of mos given an s-matrix
678 : ! orthonormality is the max deviation from unity of the C^T S C
679 : ! **************************************************************************************************
680 : !> \brief ...
681 : !> \param orthonormality ...
682 : !> \param mo_array ...
683 : !> \param matrix_s ...
684 : ! **************************************************************************************************
685 868 : SUBROUTINE calculate_orthonormality(orthonormality, mo_array, matrix_s)
686 : REAL(KIND=dp) :: orthonormality
687 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
688 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_orthonormality'
691 :
692 : INTEGER :: handle, i, ispin, j, k, n, ncol_local, &
693 : nrow_local, nspin
694 868 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
695 : REAL(KIND=dp) :: alpha, max_alpha
696 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
697 : TYPE(cp_fm_type) :: overlap, svec
698 :
699 868 : NULLIFY (tmp_fm_struct)
700 :
701 868 : CALL timeset(routineN, handle)
702 :
703 868 : nspin = SIZE(mo_array)
704 868 : max_alpha = 0.0_dp
705 :
706 1748 : DO ispin = 1, nspin
707 880 : IF (PRESENT(matrix_s)) THEN
708 : ! get S*C
709 822 : CALL cp_fm_create(svec, mo_array(ispin)%mo_coeff%matrix_struct)
710 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
711 822 : nrow_global=n, ncol_global=k)
712 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_array(ispin)%mo_coeff, &
713 822 : svec, k)
714 : ! get C^T (S*C)
715 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
716 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
717 822 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
718 822 : CALL cp_fm_create(overlap, tmp_fm_struct)
719 822 : CALL cp_fm_struct_release(tmp_fm_struct)
720 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
721 822 : svec, 0.0_dp, overlap)
722 822 : CALL cp_fm_release(svec)
723 : ELSE
724 : ! orthogonal basis C^T C
725 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
726 58 : nrow_global=n, ncol_global=k)
727 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
728 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
729 58 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
730 58 : CALL cp_fm_create(overlap, tmp_fm_struct)
731 58 : CALL cp_fm_struct_release(tmp_fm_struct)
732 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
733 58 : mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
734 : END IF
735 : CALL cp_fm_get_info(overlap, nrow_local=nrow_local, ncol_local=ncol_local, &
736 880 : row_indices=row_indices, col_indices=col_indices)
737 5031 : DO i = 1, nrow_local
738 324298 : DO j = 1, ncol_local
739 319267 : alpha = overlap%local_data(i, j)
740 319267 : IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
741 323418 : max_alpha = MAX(max_alpha, ABS(alpha))
742 : END DO
743 : END DO
744 2628 : CALL cp_fm_release(overlap)
745 : END DO
746 868 : CALL mo_array(1)%mo_coeff%matrix_struct%para_env%max(max_alpha)
747 868 : orthonormality = max_alpha
748 : ! write(*,*) "max deviation from orthonormalization ",orthonormality
749 :
750 868 : CALL timestop(handle)
751 :
752 868 : END SUBROUTINE calculate_orthonormality
753 :
754 : ! computes the minimum/maximum magnitudes of C^T C. This could be useful
755 : ! to detect problems in the case of nearly singular overlap matrices.
756 : ! in this case, we expect the ratio of min/max to be large
757 : ! this routine is only similar to mo_orthonormality if S==1
758 : ! **************************************************************************************************
759 : !> \brief ...
760 : !> \param mo_array ...
761 : !> \param mo_mag_min ...
762 : !> \param mo_mag_max ...
763 : ! **************************************************************************************************
764 848 : SUBROUTINE calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
765 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
766 : REAL(KIND=dp) :: mo_mag_min, mo_mag_max
767 :
768 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_magnitude'
769 :
770 : INTEGER :: handle, ispin, k, n, nspin
771 848 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
772 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
773 : TYPE(cp_fm_type) :: evecs, overlap
774 :
775 848 : NULLIFY (tmp_fm_struct)
776 :
777 848 : CALL timeset(routineN, handle)
778 :
779 848 : nspin = SIZE(mo_array)
780 848 : mo_mag_min = HUGE(0.0_dp)
781 848 : mo_mag_max = -HUGE(0.0_dp)
782 1696 : DO ispin = 1, nspin
783 : CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
784 848 : nrow_global=n, ncol_global=k)
785 2544 : ALLOCATE (evals(k))
786 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=k, ncol_global=k, &
787 : para_env=mo_array(ispin)%mo_coeff%matrix_struct%para_env, &
788 848 : context=mo_array(ispin)%mo_coeff%matrix_struct%context)
789 848 : CALL cp_fm_create(overlap, tmp_fm_struct)
790 848 : CALL cp_fm_create(evecs, tmp_fm_struct)
791 848 : CALL cp_fm_struct_release(tmp_fm_struct)
792 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, mo_array(ispin)%mo_coeff, &
793 848 : mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
794 848 : CALL choose_eigv_solver(overlap, evecs, evals)
795 9866 : mo_mag_min = MIN(MINVAL(evals), mo_mag_min)
796 9866 : mo_mag_max = MAX(MAXVAL(evals), mo_mag_max)
797 848 : CALL cp_fm_release(overlap)
798 848 : CALL cp_fm_release(evecs)
799 3392 : DEALLOCATE (evals)
800 : END DO
801 848 : CALL timestop(handle)
802 :
803 1696 : END SUBROUTINE calculate_magnitude
804 :
805 : ! **************************************************************************************************
806 : !> \brief Calculate KS eigenvalues starting from OF MOS
807 : !> \param mos ...
808 : !> \param nspins ...
809 : !> \param ks_rmpv ...
810 : !> \param scf_control ...
811 : !> \param mo_derivs ...
812 : !> \param admm_env ...
813 : !> \param hairy_probes ...
814 : !> \param probe ...
815 : !> \par History
816 : !> 02.2013 moved from qs_scf_post_gpw
817 : !>
818 : ! **************************************************************************************************
819 176 : SUBROUTINE make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env, &
820 : hairy_probes, probe)
821 :
822 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
823 : INTEGER, INTENT(IN) :: nspins
824 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv
825 : TYPE(scf_control_type), POINTER :: scf_control
826 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
827 : TYPE(admm_type), OPTIONAL, POINTER :: admm_env
828 : LOGICAL, OPTIONAL :: hairy_probes
829 : TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
830 : POINTER :: probe
831 :
832 : CHARACTER(len=*), PARAMETER :: routineN = 'make_mo_eig'
833 :
834 : INTEGER :: handle, homo, ispin, nmo, output_unit
835 176 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
836 : TYPE(cp_fm_type), POINTER :: mo_coeff
837 : TYPE(dbcsr_type), POINTER :: mo_coeff_deriv
838 :
839 176 : CALL timeset(routineN, handle)
840 :
841 176 : NULLIFY (mo_coeff_deriv, mo_coeff, mo_eigenvalues)
842 176 : output_unit = cp_logger_get_default_io_unit()
843 :
844 390 : DO ispin = 1, nspins
845 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
846 214 : eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
847 214 : IF (output_unit > 0) WRITE (output_unit, *) " "
848 214 : IF (output_unit > 0) WRITE (output_unit, *) " Eigenvalues of the occupied subspace spin ", ispin
849 : ! IF (homo /= nmo) THEN
850 : ! IF (output_unit>0) WRITE(output_unit,*)" and ",nmo-homo," added MO eigenvalues"
851 : ! END IF
852 214 : IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
853 :
854 214 : IF (scf_control%use_ot) THEN
855 : ! Also rotate the OT derivs, since they are needed for force calculations
856 128 : IF (ASSOCIATED(mo_derivs)) THEN
857 128 : mo_coeff_deriv => mo_derivs(ispin)%matrix
858 : ELSE
859 0 : mo_coeff_deriv => NULL()
860 : END IF
861 :
862 : ! If we do ADMM, we add have to modify the Kohn-Sham matrix
863 128 : IF (PRESENT(admm_env)) THEN
864 0 : CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
865 : END IF
866 :
867 : CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, &
868 : scr=output_unit, ionode=output_unit > 0, do_rotation=.TRUE., &
869 128 : co_rotate_dbcsr=mo_coeff_deriv)
870 :
871 : ! If we do ADMM, we restore the original Kohn-Sham matrix
872 128 : IF (PRESENT(admm_env)) THEN
873 0 : CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
874 : END IF
875 : ELSE
876 86 : IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo)
877 : END IF
878 214 : IF (.NOT. scf_control%diagonalization%mom) THEN
879 214 : IF (PRESENT(hairy_probes) .EQV. .TRUE.) THEN
880 0 : scf_control%smear%do_smear = .FALSE.
881 : CALL set_mo_occupation(mo_set=mos(ispin), &
882 : smear=scf_control%smear, &
883 0 : probe=probe)
884 : ELSE
885 214 : CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear)
886 : END IF
887 : END IF
888 214 : IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
889 497 : "Fermi Energy [eV] :", mos(ispin)%mu*evolt
890 : END DO
891 :
892 176 : CALL timestop(handle)
893 :
894 176 : END SUBROUTINE make_mo_eig
895 :
896 : END MODULE qs_mo_methods
|
