(git:08825b9)
Loading...
Searching...
No Matches
qs_mo_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
80CONTAINS
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 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 IF (ncol == 0) RETURN
103
104 CALL timeset(routinen, handle)
105
106 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
107 IF (ncol > ncol_global) cpabort("Wrong ncol value")
108
109 CALL cp_fm_create(svmatrix, vmatrix%matrix_struct, "SV", set_zero=.true.)
110 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, svmatrix, ncol)
111
112 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 context=vmatrix%matrix_struct%context)
116 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
117 CALL cp_fm_struct_release(fm_struct_tmp)
118
119 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
120 CALL cp_fm_cholesky_decompose(overlap_vv)
121 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
122
123 CALL cp_fm_release(overlap_vv)
124 CALL cp_fm_release(svmatrix)
125
126 CALL timestop(handle)
127
128 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 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 IF (ncol == 0) RETURN
152
153 CALL timeset(routinen, handle)
154 NULLIFY (fm_struct_tmp)
155
156 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
157 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 context=vmatrix%matrix_struct%context)
162 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
163 CALL cp_fm_struct_release(fm_struct_tmp)
164
165 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
166 CALL cp_fm_cholesky_decompose(overlap_vv)
167 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
168 CALL cp_fm_triangular_multiply(overlap_vv, svmatrix, n_cols=ncol, side='R', invert_tr=.true.)
169
170 CALL cp_fm_release(overlap_vv)
171
172 CALL timestop(handle)
173
174 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 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 IF (ncol == 0) RETURN
200
201 CALL timeset(routinen, handle)
202
203 !CALL cp_fm_get_info(matrix=vmatrix,nrow_global=n,ncol_global=ncol_global)
204 CALL dbcsr_get_info(vmatrix, nfullrows_total=n, nfullcols_total=ncol_global)
205 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 ncol_global=ncol, para_env=para_env)
209 CALL cp_fm_create(overlap_vv, fm_struct_tmp, name="fm_overlap_vv")
210 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 ncol_global=ncol_global, para_env=para_env)
214 CALL cp_fm_create(fm_vmatrix, fm_struct_tmp, name="fm_vmatrix")
215 CALL cp_fm_create(fm_svmatrix, fm_struct_tmp, name="fm_svmatrix")
216 CALL cp_fm_struct_release(fm_struct_tmp)
217
218 CALL copy_dbcsr_to_fm(vmatrix, fm_vmatrix)
219 CALL copy_dbcsr_to_fm(svmatrix, fm_svmatrix)
220
221 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, fm_vmatrix, fm_svmatrix, rzero, overlap_vv)
222 CALL cp_fm_cholesky_decompose(overlap_vv)
223 CALL cp_fm_triangular_multiply(overlap_vv, fm_vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
224 CALL cp_fm_triangular_multiply(overlap_vv, fm_svmatrix, n_cols=ncol, side='R', invert_tr=.true.)
225
226 CALL copy_fm_to_dbcsr(fm_vmatrix, vmatrix)
227 CALL copy_fm_to_dbcsr(fm_svmatrix, svmatrix)
228
229 CALL cp_fm_release(overlap_vv)
230 CALL cp_fm_release(fm_vmatrix)
231 CALL cp_fm_release(fm_svmatrix)
232
233 CALL timestop(handle)
234
235 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 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 IF (ncol == 0) RETURN
264
265 CALL timeset(routinen, handle)
266 NULLIFY (fm_struct_tmp)
267
268 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
269 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 context=vmatrix%matrix_struct%context)
274 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
275 CALL cp_fm_struct_release(fm_struct_tmp)
276
277 CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol)
278 CALL cp_fm_syrk('U', 'T', n, rone, vmatrix, 1, 1, rzero, overlap_vv)
279 CALL cp_fm_cholesky_decompose(overlap_vv)
280 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
281 CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol, invert_tr=.true.)
282
283 CALL cp_fm_release(overlap_vv)
284
285 CALL timestop(handle)
286
287 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 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 IF (ncol == 0) RETURN
316
317 CALL timeset(routinen, handle)
318 NULLIFY (fm_struct_tmp)
319 threshold = 1.0e-7_dp
320 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
321 IF (ncol > ncol_global) cpabort("Wrong ncol value")
322
323 CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
324 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
325
326 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 context=vmatrix%matrix_struct%context)
330 CALL cp_fm_create(csc, fm_struct_tmp, "csc")
331 CALL cp_fm_set_all(csc, 0.0_dp)
332 CALL cp_fm_create(work, fm_struct_tmp, "work")
333 CALL cp_fm_set_all(work, 0.0_dp)
334 CALL cp_fm_struct_release(fm_struct_tmp)
335
336 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
337 CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
338 CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
339 CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
340
341 CALL cp_fm_release(csc)
342 CALL cp_fm_release(sc)
343 CALL cp_fm_release(work)
344
345 CALL timestop(handle)
346
347 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 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 IF (ncol == 0) RETURN
370
371 CALL timeset(routinen, handle)
372
373 NULLIFY (fm_struct_tmp)
374
375 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
376 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 context=vmatrix%matrix_struct%context)
381 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
382 CALL cp_fm_struct_release(fm_struct_tmp)
383
384 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, vmatrix, rzero, overlap_vv)
385 CALL cp_fm_cholesky_decompose(overlap_vv)
386 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
387
388 CALL cp_fm_release(overlap_vv)
389
390 CALL timestop(handle)
391
392 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 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 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 CALL timeset(routinen, handle)
433
434 do_rotation_local = .true.
435 IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
436
437 NULLIFY (fm_struct_tmp)
438 CALL cp_fm_get_info(matrix=orbitals, &
439 ncol_global=ncol_global, &
440 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 IF (ncol_global > 0) THEN
452
453 ALLOCATE (evals(ncol_global))
454
455 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 context=orbitals%matrix_struct%context)
459 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
460 IF (compute_evecs) THEN
461 CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e vectors")
462 END IF
463 CALL cp_fm_struct_release(fm_struct_tmp)
464
465 ! h subblock and diag
466 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 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 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 IF (do_rotation_local) THEN
480 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
481 orbitals, e_vectors, 0.0_dp, weighted_vectors)
482 CALL cp_fm_to_fm(weighted_vectors, orbitals)
483 IF (PRESENT(co_rotate)) THEN
484 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
485 co_rotate, e_vectors, 0.0_dp, weighted_vectors)
486 CALL cp_fm_to_fm(weighted_vectors, co_rotate)
487 END IF
488 IF (PRESENT(co_rotate_dbcsr)) THEN
489 IF (ASSOCIATED(co_rotate_dbcsr)) THEN
490 CALL cp_fm_create(weighted_vectors2, orbitals%matrix_struct, "weighted_vectors")
491 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 weighted_vectors2, e_vectors, 0.0_dp, weighted_vectors)
494 CALL copy_fm_to_dbcsr(weighted_vectors, co_rotate_dbcsr)
495 CALL cp_fm_release(weighted_vectors2)
496 END IF
497 END IF
498 END IF
499
500 ! give output
501 IF (PRESENT(evals_arg)) THEN
502 n = min(SIZE(evals_arg), SIZE(evals))
503 evals_arg(1:n) = evals(1:n)
504 END IF
505
506 IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
507 IF (.NOT. PRESENT(ionode)) cpabort("IONODE?")
508 IF (.NOT. PRESENT(scr)) cpabort("SCR?")
509 IF (ionode) THEN
510 DO i = 1, ncol_global, 4
511 j = min(3, ncol_global - i)
512 SELECT CASE (j)
513 CASE (3)
514 WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
515 CASE (2)
516 WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
517 CASE (1)
518 WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
519 CASE (0)
520 WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
521 END SELECT
522 END DO
523 END IF
524 END IF
525
526 CALL cp_fm_release(weighted_vectors)
527 CALL cp_fm_release(h_block)
528 IF (compute_evecs) THEN
529 CALL cp_fm_release(e_vectors)
530 END IF
531
532 DEALLOCATE (evals)
533
534 END IF
535
536 CALL timestop(handle)
537
538 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 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 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
569 TYPE(dbcsr_type), POINTER :: e_vectors, h_block, weighted_vectors
570
571 CALL timeset(routinen, handle)
572
573 do_rotation_local = .true.
574 IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
575
576 NULLIFY (e_vectors, h_block, weighted_vectors)
577
578 CALL dbcsr_get_info(matrix=orbitals, &
579 nfullcols_total=ncol_global, &
580 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 IF (ncol_global > 0) THEN
592
593 ALLOCATE (evals(ncol_global))
594
595 CALL dbcsr_init_p(weighted_vectors)
596 CALL dbcsr_copy(weighted_vectors, orbitals, name="weighted_vectors")
597
598 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 sym=dbcsr_type_no_symmetry)
601
602!!!!!!!!!!!!!!XXXXXXXXXXXXXXXXXXX!!!!!!!!!!!!!
603 IF (compute_evecs) THEN
604 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 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 0.0_dp, weighted_vectors)
612 !CALL cp_dbcsr_sm_fm_multiply(ks_matrix,orbitals,weighted_vectors, ncol_global)
613
614 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 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 IF (do_rotation_local) THEN
627 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 CALL dbcsr_copy(orbitals, weighted_vectors)
631 !CALL cp_fm_to_fm(weighted_vectors,orbitals)
632 IF (PRESENT(co_rotate)) THEN
633 IF (ASSOCIATED(co_rotate)) THEN
634 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 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 IF (PRESENT(evals_arg)) THEN
645 evals_arg(:) = evals(:)
646 END IF
647
648 IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
649 IF (.NOT. PRESENT(ionode)) cpabort("IONODE?")
650 IF (.NOT. PRESENT(scr)) cpabort("SCR?")
651 IF (ionode) THEN
652 DO i = 1, ncol_global, 4
653 j = min(3, ncol_global - i)
654 SELECT CASE (j)
655 CASE (3)
656 WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
657 CASE (2)
658 WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
659 CASE (1)
660 WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
661 CASE (0)
662 WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
663 END SELECT
664 END DO
665 END IF
666 END IF
667
668 CALL dbcsr_release_p(weighted_vectors)
669 CALL dbcsr_release_p(h_block)
670 IF (compute_evecs) THEN
671 CALL dbcsr_release_p(e_vectors)
672 END IF
673
674 DEALLOCATE (evals)
675
676 END IF
677
678 CALL timestop(handle)
679
680 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 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 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 NULLIFY (tmp_fm_struct)
705
706 CALL timeset(routinen, handle)
707
708 nspin = SIZE(mo_array)
709 max_alpha = 0.0_dp
710
711 DO ispin = 1, nspin
712 IF (PRESENT(matrix_s)) THEN
713 ! get S*C
714 CALL cp_fm_create(svec, mo_array(ispin)%mo_coeff%matrix_struct)
715 CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
716 nrow_global=n, ncol_global=k)
717 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_array(ispin)%mo_coeff, &
718 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
723 CALL cp_fm_create(overlap, tmp_fm_struct)
724 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 svec, 0.0_dp, overlap)
727 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 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
735 CALL cp_fm_create(overlap, tmp_fm_struct)
736 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 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 row_indices=row_indices, col_indices=col_indices)
742 DO i = 1, nrow_local
743 DO j = 1, ncol_local
744 alpha = overlap%local_data(i, j)
745 IF (row_indices(i) == col_indices(j)) alpha = alpha - 1.0_dp
746 max_alpha = max(max_alpha, abs(alpha))
747 END DO
748 END DO
749 CALL cp_fm_release(overlap)
750 END DO
751 CALL mo_array(1)%mo_coeff%matrix_struct%para_env%max(max_alpha)
752 orthonormality = max_alpha
753 ! write(*,*) "max deviation from orthonormalization ",orthonormality
754
755 CALL timestop(handle)
756
757 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 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 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 NULLIFY (tmp_fm_struct)
781
782 CALL timeset(routinen, handle)
783
784 nspin = SIZE(mo_array)
785 mo_mag_min = huge(0.0_dp)
786 mo_mag_max = -huge(0.0_dp)
787 DO ispin = 1, nspin
788 CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
789 nrow_global=n, ncol_global=k)
790 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
794 CALL cp_fm_create(overlap, tmp_fm_struct)
795 CALL cp_fm_create(evecs, tmp_fm_struct)
796 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 mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
799 CALL choose_eigv_solver(overlap, evecs, evals)
800 mo_mag_min = min(minval(evals), mo_mag_min)
801 mo_mag_max = max(maxval(evals), mo_mag_max)
802 CALL cp_fm_release(overlap)
803 CALL cp_fm_release(evecs)
804 DEALLOCATE (evals)
805 END DO
806 CALL timestop(handle)
807
808 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 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 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 CALL timeset(routinen, handle)
845
846 NULLIFY (mo_coeff_deriv, mo_coeff, mo_eigenvalues)
847 output_unit = cp_logger_get_default_io_unit()
848
849 DO ispin = 1, nspins
850 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
851 eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
852 IF (output_unit > 0) WRITE (output_unit, *) " "
853 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 IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
858
859 IF (scf_control%use_ot) THEN
860 ! Also rotate the OT derivs, since they are needed for force calculations
861 IF (ASSOCIATED(mo_derivs)) THEN
862 mo_coeff_deriv => mo_derivs(ispin)%matrix
863 ELSE
864 mo_coeff_deriv => null()
865 END IF
866
867 ! If we do ADMM, we add have to modify the Kohn-Sham matrix
868 IF (PRESENT(admm_env)) THEN
869 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 co_rotate_dbcsr=mo_coeff_deriv)
875
876 ! If we do ADMM, we restore the original Kohn-Sham matrix
877 IF (PRESENT(admm_env)) THEN
878 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
879 END IF
880 ELSE
881 IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo)
882 END IF
883 IF (.NOT. scf_control%diagonalization%mom) THEN
884 IF (PRESENT(hairy_probes) .EQV. .true.) THEN
885 scf_control%smear%do_smear = .false.
886 CALL set_mo_occupation(mo_set=mos(ispin), &
887 smear=scf_control%smear, &
888 probe=probe)
889 ELSE
890 CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear)
891 END IF
892 END IF
893 IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
894 "Fermi Energy [eV] :", mos(ispin)%mu*evolt
895 END DO
896
897 CALL timestop(handle)
898
899 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 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 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 CALL timeset(routinen, handle)
925
926 orbital_energies = 0.0_dp
927 max_offdiag = 0.0_dp
928 CALL cp_fm_get_info(matrix=orbitals, ncol_global=ncol_global, nrow_global=nrow_global)
929 IF (ncol_global > 0) THEN
930 NULLIFY (fm_struct_tmp)
931 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 context=orbitals%matrix_struct%context)
935 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
936 CALL cp_fm_struct_release(fm_struct_tmp)
937
938 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 orbitals, weighted_vectors, 0.0_dp, h_block)
941
942 ALLOCATE (fock(ncol_global, ncol_global))
943 CALL cp_fm_get_submatrix(h_block, fock)
944 DO i = 1, min(ncol_global, SIZE(orbital_energies))
945 orbital_energies(i) = fock(i, i)
946 END DO
947 DO i = 1, ncol_global
948 fock(i, i) = 0.0_dp
949 END DO
950 max_offdiag = maxval(abs(fock))
951 DEALLOCATE (fock)
952
953 CALL cp_fm_release(weighted_vectors)
954 CALL cp_fm_release(h_block)
955 END IF
956
957 CALL timestop(handle)
958
959 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 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 CALL timeset(routinen, handle)
998
999 rotated = .false.
1000 nmo = mo_set%nmo
1001 IF (mo_set%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mo_set%mo_coeff_b, mo_set%mo_coeff)
1002 ALLOCATE (energies(nmo))
1003 CALL calculate_orbital_fock_elements(mo_set%mo_coeff, ks_matrix, energies, max_offdiag)
1004 DEALLOCATE (energies)
1005
1006 IF (max_offdiag > eps_canonical) THEN
1007 DO i = 1, nmo
1008 IF (abs(mo_set%occupation_numbers(i)) > eps_occ .AND. &
1009 abs(mo_set%occupation_numbers(i) - mo_set%maxocc) > eps_occ) THEN
1010 cpabort("Canonicalization of orbitals with fractional occupations is not possible")
1011 END IF
1012 END DO
1013 do_derivs = .false.
1014 IF (PRESENT(mo_derivs)) do_derivs = ASSOCIATED(mo_derivs)
1015 ! occupied orbitals: those with nonzero occupation, they come first
1016 nocc = count(mo_set%occupation_numbers(1:nmo) > eps_occ)
1017 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, nrow_global=nrow_global)
1018 IF (do_derivs) THEN
1019 CALL cp_fm_create(derivs, mo_set%mo_coeff%matrix_struct, name="orbital derivatives")
1020 CALL copy_dbcsr_to_fm(mo_derivs, derivs)
1021 END IF
1022 DO iblock = 1, 2
1023 IF (iblock == 1) THEN
1024 first = 1
1025 n = nocc
1026 ELSE
1027 first = nocc + 1
1028 n = nmo - nocc
1029 END IF
1030 IF (n <= 0) cycle
1031 NULLIFY (fm_struct_tmp)
1032 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, ncol_global=n, &
1033 template_fmstruct=mo_set%mo_coeff%matrix_struct)
1034 CALL cp_fm_create(block_coeff, fm_struct_tmp, name="canonicalized orbitals")
1035 CALL cp_fm_to_fm(mo_set%mo_coeff, block_coeff, n, first, 1)
1036 IF (do_derivs) THEN
1037 CALL cp_fm_create(block_derivs, fm_struct_tmp, name="rotated orbital derivatives")
1038 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 do_rotation=.true., co_rotate=block_derivs)
1042 CALL cp_fm_to_fm(block_derivs, derivs, n, 1, first)
1043 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 do_rotation=.true.)
1048 END IF
1049 CALL cp_fm_struct_release(fm_struct_tmp)
1050 CALL cp_fm_to_fm(block_coeff, mo_set%mo_coeff, n, 1, first)
1051 CALL cp_fm_release(block_coeff)
1052 END DO
1053 IF (do_derivs) THEN
1054 CALL copy_fm_to_dbcsr(derivs, mo_derivs)
1055 CALL cp_fm_release(derivs)
1056 END IF
1057 IF (mo_set%use_mo_coeff_b) CALL copy_fm_to_dbcsr(mo_set%mo_coeff, mo_set%mo_coeff_b)
1058 rotated = .true.
1059 END IF
1060
1061 CALL timestop(handle)
1062
1063 END SUBROUTINE make_mo_set_canonical
1064
1065END MODULE qs_mo_methods
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:70
qs_mo_methods::make_basis_sv
Definition qs_mo_methods.F:75
qs_mo_occupation::set_mo_occupation
Definition qs_mo_occupation.F:50
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_utils
Contains methods used in the context of density fitting.
Definition admm_utils.F:15
admm_utils::admm_uncorrect_for_eigenvalues
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:127
admm_utils::admm_correct_for_eigenvalues
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:53
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_syevx
subroutine, public cp_dbcsr_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx, para_env, blacs_env)
compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack....
Definition cp_dbcsr_diag.F:126
cp_dbcsr_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:67
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:829
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:233
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:1173
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_syrk
subroutine, public cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
performs a rank-k update of a symmetric matrix_c matrix_c = beta * matrix_c + alpha * matrix_a * tran...
Definition cp_fm_basic_linalg.F:666
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_fm_basic_linalg.F:1414
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:62
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1144
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_diag::cp_fm_syevx
subroutine, public cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack....
Definition cp_fm_diag.F:814
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1044
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:931
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
hairy_probes
Definition hairy_probes.F:8
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_methods::make_basis_simple
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
Definition qs_mo_methods.F:358
qs_mo_methods::make_basis_lowdin
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
Definition qs_mo_methods.F:302
qs_mo_methods::make_basis_cholesky
subroutine, public make_basis_cholesky(vmatrix, ncol, ortho)
return a set of S orthonormal vectors (C^T S C == 1) where the cholesky decomposed form of S is passe...
Definition qs_mo_methods.F:251
qs_mo_methods::calculate_magnitude
subroutine, public calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
...
Definition qs_mo_methods.F:770
qs_mo_methods::make_mo_set_canonical
subroutine, public make_mo_set_canonical(mo_set, ks_matrix, eps_canonical, rotated, max_offdiag, mo_derivs)
Rotates the orbitals of an MO set into canonical orbitals (eigenvectors of the Fock matrix in the spa...
Definition qs_mo_methods.F:980
qs_mo_methods::calculate_orbital_fock_elements
subroutine, public calculate_orbital_fock_elements(orbitals, ks_matrix, orbital_energies, max_offdiag)
Orbital energies <phi_i|F|phi_i> of a set of orbitals and the largest off-diagonal Fock matrix elemen...
Definition qs_mo_methods.F:912
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition qs_mo_methods.F:91
qs_mo_methods::make_mo_eig
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env, hairy_probes, probe)
Calculate KS eigenvalues starting from OF MOS.
Definition qs_mo_methods.F:826
qs_mo_methods::calculate_orthonormality
subroutine, public calculate_orthonormality(orthonormality, mo_array, matrix_s)
...
Definition qs_mo_methods.F:691
qs_mo_occupation
Set occupation of molecular orbitals.
Definition qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::hairy_probes_type
Definition cp_control_types.F:49
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
scf_control_types::scf_control_type
Definition scf_control_types.F:135