(git:374b731)
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-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! **************************************************************************************************
15 USE admm_types, ONLY: admm_type
19 USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd,&
34 USE cp_fm_types, ONLY: cp_fm_create,&
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
51 USE physcon, ONLY: evolt
53 USE qs_mo_types, ONLY: get_mo_set,&
56#include "./base/base_uses.f90"
57
58 IMPLICIT NONE
59 PRIVATE
60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_methods'
61
65
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
76CONTAINS
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 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 IF (ncol .EQ. 0) RETURN
99
100 CALL timeset(routinen, handle)
101
102 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
103 IF (ncol .GT. ncol_global) cpabort("Wrong ncol value")
104
105 CALL cp_fm_create(svmatrix, vmatrix%matrix_struct, "SV")
106 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, svmatrix, ncol)
107
108 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 context=vmatrix%matrix_struct%context)
112 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
113 CALL cp_fm_struct_release(fm_struct_tmp)
114
115 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
116 CALL cp_fm_cholesky_decompose(overlap_vv)
117 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
118
119 CALL cp_fm_release(overlap_vv)
120 CALL cp_fm_release(svmatrix)
121
122 CALL timestop(handle)
123
124 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 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 IF (ncol .EQ. 0) RETURN
148
149 CALL timeset(routinen, handle)
150 NULLIFY (fm_struct_tmp)
151
152 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
153 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 context=vmatrix%matrix_struct%context)
158 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
159 CALL cp_fm_struct_release(fm_struct_tmp)
160
161 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, svmatrix, rzero, overlap_vv)
162 CALL cp_fm_cholesky_decompose(overlap_vv)
163 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
164 CALL cp_fm_triangular_multiply(overlap_vv, svmatrix, n_cols=ncol, side='R', invert_tr=.true.)
165
166 CALL cp_fm_release(overlap_vv)
167
168 CALL timestop(handle)
169
170 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 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 IF (ncol .EQ. 0) RETURN
196
197 CALL timeset(routinen, handle)
198
199 !CALL cp_fm_get_info(matrix=vmatrix,nrow_global=n,ncol_global=ncol_global)
200 CALL dbcsr_get_info(vmatrix, nfullrows_total=n, nfullcols_total=ncol_global)
201 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 ncol_global=ncol, para_env=para_env)
205 CALL cp_fm_create(overlap_vv, fm_struct_tmp, name="fm_overlap_vv")
206 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 ncol_global=ncol_global, para_env=para_env)
210 CALL cp_fm_create(fm_vmatrix, fm_struct_tmp, name="fm_vmatrix")
211 CALL cp_fm_create(fm_svmatrix, fm_struct_tmp, name="fm_svmatrix")
212 CALL cp_fm_struct_release(fm_struct_tmp)
213
214 CALL copy_dbcsr_to_fm(vmatrix, fm_vmatrix)
215 CALL copy_dbcsr_to_fm(svmatrix, fm_svmatrix)
216
217 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, fm_vmatrix, fm_svmatrix, rzero, overlap_vv)
218 CALL cp_fm_cholesky_decompose(overlap_vv)
219 CALL cp_fm_triangular_multiply(overlap_vv, fm_vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
220 CALL cp_fm_triangular_multiply(overlap_vv, fm_svmatrix, n_cols=ncol, side='R', invert_tr=.true.)
221
222 CALL copy_fm_to_dbcsr(fm_vmatrix, vmatrix)
223 CALL copy_fm_to_dbcsr(fm_svmatrix, svmatrix)
224
225 CALL cp_fm_release(overlap_vv)
226 CALL cp_fm_release(fm_vmatrix)
227 CALL cp_fm_release(fm_svmatrix)
228
229 CALL timestop(handle)
230
231 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 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 IF (ncol .EQ. 0) RETURN
260
261 CALL timeset(routinen, handle)
262 NULLIFY (fm_struct_tmp)
263
264 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
265 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 context=vmatrix%matrix_struct%context)
270 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
271 CALL cp_fm_struct_release(fm_struct_tmp)
272
273 CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol)
274 CALL cp_fm_syrk('U', 'T', n, rone, vmatrix, 1, 1, rzero, overlap_vv)
275 CALL cp_fm_cholesky_decompose(overlap_vv)
276 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
277 CALL cp_fm_triangular_multiply(ortho, vmatrix, n_cols=ncol, invert_tr=.true.)
278
279 CALL cp_fm_release(overlap_vv)
280
281 CALL timestop(handle)
282
283 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 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 IF (ncol .EQ. 0) RETURN
312
313 CALL timeset(routinen, handle)
314 NULLIFY (fm_struct_tmp)
315 threshold = 1.0e-7_dp
316 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
317 IF (ncol .GT. ncol_global) cpabort("Wrong ncol value")
318
319 CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
320 CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
321
322 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 context=vmatrix%matrix_struct%context)
326 CALL cp_fm_create(csc, fm_struct_tmp, "csc")
327 CALL cp_fm_create(work, fm_struct_tmp, "work")
328 CALL cp_fm_struct_release(fm_struct_tmp)
329
330 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
331 CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
332 CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
333 CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
334
335 CALL cp_fm_release(csc)
336 CALL cp_fm_release(sc)
337 CALL cp_fm_release(work)
338
339 CALL timestop(handle)
340
341 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 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 IF (ncol .EQ. 0) RETURN
364
365 CALL timeset(routinen, handle)
366
367 NULLIFY (fm_struct_tmp)
368
369 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
370 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 context=vmatrix%matrix_struct%context)
375 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
376 CALL cp_fm_struct_release(fm_struct_tmp)
377
378 CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, vmatrix, rzero, overlap_vv)
379 CALL cp_fm_cholesky_decompose(overlap_vv)
380 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
381
382 CALL cp_fm_release(overlap_vv)
383
384 CALL timestop(handle)
385
386 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 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 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 CALL timeset(routinen, handle)
427
428 do_rotation_local = .true.
429 IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
430
431 NULLIFY (fm_struct_tmp)
432 CALL cp_fm_get_info(matrix=orbitals, &
433 ncol_global=ncol_global, &
434 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 IF (ncol_global .GT. 0) THEN
446
447 ALLOCATE (evals(ncol_global))
448
449 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 context=orbitals%matrix_struct%context)
453 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
454 IF (compute_evecs) THEN
455 CALL cp_fm_create(e_vectors, fm_struct_tmp, name="e vectors")
456 END IF
457 CALL cp_fm_struct_release(fm_struct_tmp)
458
459 ! h subblock and diag
460 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 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 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 IF (do_rotation_local) THEN
474 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
475 orbitals, e_vectors, 0.0_dp, weighted_vectors)
476 CALL cp_fm_to_fm(weighted_vectors, orbitals)
477 IF (PRESENT(co_rotate)) THEN
478 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
479 co_rotate, e_vectors, 0.0_dp, weighted_vectors)
480 CALL cp_fm_to_fm(weighted_vectors, co_rotate)
481 END IF
482 IF (PRESENT(co_rotate_dbcsr)) THEN
483 IF (ASSOCIATED(co_rotate_dbcsr)) THEN
484 CALL cp_fm_create(weighted_vectors2, orbitals%matrix_struct, "weighted_vectors")
485 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 weighted_vectors2, e_vectors, 0.0_dp, weighted_vectors)
488 CALL copy_fm_to_dbcsr(weighted_vectors, co_rotate_dbcsr)
489 CALL cp_fm_release(weighted_vectors2)
490 END IF
491 END IF
492 END IF
493
494 ! give output
495 IF (PRESENT(evals_arg)) THEN
496 n = min(SIZE(evals_arg), SIZE(evals))
497 evals_arg(1:n) = evals(1:n)
498 END IF
499
500 IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
501 IF (.NOT. PRESENT(ionode)) cpabort("IONODE?")
502 IF (.NOT. PRESENT(scr)) cpabort("SCR?")
503 IF (ionode) THEN
504 DO i = 1, ncol_global, 4
505 j = min(3, ncol_global - i)
506 SELECT CASE (j)
507 CASE (3)
508 WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
509 CASE (2)
510 WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
511 CASE (1)
512 WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
513 CASE (0)
514 WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
515 END SELECT
516 END DO
517 END IF
518 END IF
519
520 CALL cp_fm_release(weighted_vectors)
521 CALL cp_fm_release(h_block)
522 IF (compute_evecs) THEN
523 CALL cp_fm_release(e_vectors)
524 END IF
525
526 DEALLOCATE (evals)
527
528 END IF
529
530 CALL timestop(handle)
531
532 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 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 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
563 TYPE(dbcsr_type), POINTER :: e_vectors, h_block, weighted_vectors
564
565 CALL timeset(routinen, handle)
566
567 do_rotation_local = .true.
568 IF (PRESENT(do_rotation)) do_rotation_local = do_rotation
569
570 NULLIFY (e_vectors, h_block, weighted_vectors)
571
572 CALL dbcsr_get_info(matrix=orbitals, &
573 nfullcols_total=ncol_global, &
574 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 IF (ncol_global .GT. 0) THEN
586
587 ALLOCATE (evals(ncol_global))
588
589 CALL dbcsr_init_p(weighted_vectors)
590 CALL dbcsr_copy(weighted_vectors, orbitals, name="weighted_vectors")
591
592 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 sym=dbcsr_type_no_symmetry)
595
596!!!!!!!!!!!!!!XXXXXXXXXXXXXXXXXXX!!!!!!!!!!!!!
597 IF (compute_evecs) THEN
598 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 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 0.0_dp, weighted_vectors)
606 !CALL cp_dbcsr_sm_fm_multiply(ks_matrix,orbitals,weighted_vectors, ncol_global)
607
608 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 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 IF (do_rotation_local) THEN
621 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 CALL dbcsr_copy(orbitals, weighted_vectors)
625 !CALL cp_fm_to_fm(weighted_vectors,orbitals)
626 IF (PRESENT(co_rotate)) THEN
627 IF (ASSOCIATED(co_rotate)) THEN
628 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 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 IF (PRESENT(evals_arg)) THEN
639 evals_arg(:) = evals(:)
640 END IF
641
642 IF (PRESENT(ionode) .OR. PRESENT(scr)) THEN
643 IF (.NOT. PRESENT(ionode)) cpabort("IONODE?")
644 IF (.NOT. PRESENT(scr)) cpabort("SCR?")
645 IF (ionode) THEN
646 DO i = 1, ncol_global, 4
647 j = min(3, ncol_global - i)
648 SELECT CASE (j)
649 CASE (3)
650 WRITE (scr, '(1X,4F16.8)') evals(i:i + j)
651 CASE (2)
652 WRITE (scr, '(1X,3F16.8)') evals(i:i + j)
653 CASE (1)
654 WRITE (scr, '(1X,2F16.8)') evals(i:i + j)
655 CASE (0)
656 WRITE (scr, '(1X,1F16.8)') evals(i:i + j)
657 END SELECT
658 END DO
659 END IF
660 END IF
661
662 CALL dbcsr_release_p(weighted_vectors)
663 CALL dbcsr_release_p(h_block)
664 IF (compute_evecs) THEN
665 CALL dbcsr_release_p(e_vectors)
666 END IF
667
668 DEALLOCATE (evals)
669
670 END IF
671
672 CALL timestop(handle)
673
674 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 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 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 NULLIFY (tmp_fm_struct)
699
700 CALL timeset(routinen, handle)
701
702 nspin = SIZE(mo_array)
703 max_alpha = 0.0_dp
704
705 DO ispin = 1, nspin
706 IF (PRESENT(matrix_s)) THEN
707 ! get S*C
708 CALL cp_fm_create(svec, mo_array(ispin)%mo_coeff%matrix_struct)
709 CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
710 nrow_global=n, ncol_global=k)
711 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_array(ispin)%mo_coeff, &
712 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
717 CALL cp_fm_create(overlap, tmp_fm_struct)
718 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 svec, 0.0_dp, overlap)
721 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 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
729 CALL cp_fm_create(overlap, tmp_fm_struct)
730 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 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 row_indices=row_indices, col_indices=col_indices)
736 DO i = 1, nrow_local
737 DO j = 1, ncol_local
738 alpha = overlap%local_data(i, j)
739 IF (row_indices(i) .EQ. col_indices(j)) alpha = alpha - 1.0_dp
740 max_alpha = max(max_alpha, abs(alpha))
741 END DO
742 END DO
743 CALL cp_fm_release(overlap)
744 END DO
745 CALL mo_array(1)%mo_coeff%matrix_struct%para_env%max(max_alpha)
746 orthonormality = max_alpha
747 ! write(*,*) "max deviation from orthonormalization ",orthonormality
748
749 CALL timestop(handle)
750
751 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 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 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 NULLIFY (tmp_fm_struct)
775
776 CALL timeset(routinen, handle)
777
778 nspin = SIZE(mo_array)
779 mo_mag_min = huge(0.0_dp)
780 mo_mag_max = -huge(0.0_dp)
781 DO ispin = 1, nspin
782 CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, &
783 nrow_global=n, ncol_global=k)
784 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 context=mo_array(ispin)%mo_coeff%matrix_struct%context)
788 CALL cp_fm_create(overlap, tmp_fm_struct)
789 CALL cp_fm_create(evecs, tmp_fm_struct)
790 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 mo_array(ispin)%mo_coeff, 0.0_dp, overlap)
793 CALL choose_eigv_solver(overlap, evecs, evals)
794 mo_mag_min = min(minval(evals), mo_mag_min)
795 mo_mag_max = max(maxval(evals), mo_mag_max)
796 CALL cp_fm_release(overlap)
797 CALL cp_fm_release(evecs)
798 DEALLOCATE (evals)
799 END DO
800 CALL timestop(handle)
801
802 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 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 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 CALL timeset(routinen, handle)
833
834 NULLIFY (mo_coeff_deriv, mo_coeff, mo_eigenvalues)
835 output_unit = cp_logger_get_default_io_unit()
836
837 DO ispin = 1, nspins
838 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
839 eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
840 IF (output_unit > 0) WRITE (output_unit, *) " "
841 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 IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
846
847 IF (scf_control%use_ot) THEN
848 ! Also rotate the OT derivs, since they are needed for force calculations
849 IF (ASSOCIATED(mo_derivs)) THEN
850 mo_coeff_deriv => mo_derivs(ispin)%matrix
851 ELSE
852 mo_coeff_deriv => null()
853 END IF
854
855 ! ** If we do ADMM, we add have to modify the kohn-sham matrix
856 IF (PRESENT(admm_env)) THEN
857 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 co_rotate_dbcsr=mo_coeff_deriv)
863
864 ! ** If we do ADMM, we restore the original kohn-sham matrix
865 IF (PRESENT(admm_env)) THEN
866 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
867 END IF
868 ELSE
869 IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(1:homo)
870 END IF
871 IF (.NOT. scf_control%diagonalization%mom) &
872 CALL set_mo_occupation(mo_set=mos(ispin), smear=scf_control%smear)
873 IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
874 "Fermi Energy [eV] :", mos(ispin)%mu*evolt
875 END DO
876
877 CALL timestop(handle)
878
879 END SUBROUTINE make_mo_eig
880
881END MODULE qs_mo_methods
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
Contains methods used in the context of density fitting.
Definition admm_utils.F:15
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:127
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:53
methods related to the blacs parallel environment
Interface to (sca)lapack for the Cholesky based procedures.
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....
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
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...
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...
various cholesky decomposition related routines
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,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:896
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:208
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:657
represent the structure of a full matrix
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
collects routines that perform operations directly related to MOs
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,...
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...
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env)
Calculate KS eigenvalues starting from OF MOS.
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...
subroutine, public calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
...
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
subroutine, public calculate_orthonormality(orthonormality, mo_array, matrix_s)
...
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
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)
Get the components of a MO set data structure.
parameters that control an scf iteration
stores some data used in wavefunction fitting
Definition admm_types.F:120
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
stores all the informations relevant to an mpi environment