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