(git:374b731)
Loading...
Searching...
No Matches
mao_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 Calculate MAO's and analyze wavefunctions
10!> \par History
11!> 03.2016 created [JGH]
12!> 12.2016 split into four modules [JGH]
13!> \author JGH
14! **************************************************************************************************
28 USE cp_fm_diag, ONLY: cp_fm_geeig
32 USE cp_fm_types, ONLY: cp_fm_create,&
35 USE dbcsr_api, ONLY: &
36 dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, dbcsr_get_block_p, &
37 dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
38 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
39 dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_set, dbcsr_type, &
40 dbcsr_type_no_symmetry
45 USE kinds, ONLY: dp
47 USE kpoint_types, ONLY: get_kpoint_info,&
49 USE lapack, ONLY: lapack_ssyev,&
51 USE message_passing, ONLY: mp_comm_type,&
57 USE qs_kind_types, ONLY: get_qs_kind,&
60#include "./base/base_uses.f90"
61
62 IMPLICIT NONE
63 PRIVATE
64
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_methods'
66
67 TYPE mblocks
68 INTEGER :: n, ma
69 REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat
70 REAL(KIND=dp), DIMENSION(:), POINTER :: eig
71 END TYPE mblocks
72
76
77! **************************************************************************************************
78
79CONTAINS
80
81! **************************************************************************************************
82!> \brief ...
83!> \param mao_coef ...
84!> \param pmat ...
85!> \param smat ...
86!> \param eps1 ...
87!> \param iolevel ...
88!> \param iw ...
89! **************************************************************************************************
90 SUBROUTINE mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
91 TYPE(dbcsr_type) :: mao_coef, pmat, smat
92 REAL(kind=dp), INTENT(IN) :: eps1
93 INTEGER, INTENT(IN) :: iolevel, iw
94
95 INTEGER :: group_handle, i, iatom, info, jatom, &
96 lwork, m, n, nblk
97 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, row_blk, &
98 row_blk_sizes
99 LOGICAL :: found
100 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
101 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
102 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, pblock, sblock
103 TYPE(dbcsr_distribution_type) :: dbcsr_dist
104 TYPE(dbcsr_iterator_type) :: dbcsr_iter
105 TYPE(mblocks), ALLOCATABLE, DIMENSION(:) :: mbl
106 TYPE(mp_comm_type) :: group
107
108 CALL dbcsr_get_info(mao_coef, nblkrows_total=nblk)
109 ALLOCATE (mbl(nblk))
110 DO i = 1, nblk
111 NULLIFY (mbl(i)%mat, mbl(i)%eig)
112 END DO
113
114 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
115 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
116 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
117 cpassert(iatom == jatom)
118 m = SIZE(cblock, 2)
119 NULLIFY (pblock, sblock)
120 CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pblock, found=found)
121 cpassert(found)
122 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
123 cpassert(found)
124 n = SIZE(sblock, 1)
125 lwork = max(n*n, 100)
126 ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
127 amat(1:n, 1:n) = pblock(1:n, 1:n)
128 bmat(1:n, 1:n) = sblock(1:n, 1:n)
129 info = 0
130 CALL lapack_ssygv(1, "V", "U", n, amat, n, bmat, n, w, work, lwork, info)
131 cpassert(info == 0)
132 ALLOCATE (mbl(iatom)%mat(n, n), mbl(iatom)%eig(n))
133 mbl(iatom)%n = n
134 mbl(iatom)%ma = m
135 DO i = 1, n
136 mbl(iatom)%eig(i) = w(n - i + 1)
137 mbl(iatom)%mat(1:n, i) = amat(1:n, n - i + 1)
138 END DO
139 cblock(1:n, 1:m) = amat(1:n, n:n - m + 1:-1)
140 DEALLOCATE (amat, bmat, w, work)
141 END DO
142 CALL dbcsr_iterator_stop(dbcsr_iter)
143
144 IF (eps1 < 10.0_dp) THEN
145 CALL dbcsr_get_info(mao_coef, row_blk_size=row_blk_sizes, group=group_handle)
146 CALL group%set_handle(group_handle)
147 ALLOCATE (row_blk(nblk), mao_blk(nblk))
148 mao_blk = 0
149 row_blk = row_blk_sizes
150 DO iatom = 1, nblk
151 IF (ASSOCIATED(mbl(iatom)%mat)) THEN
152 n = mbl(iatom)%n
153 m = 0
154 DO i = 1, n
155 IF (mbl(iatom)%eig(i) < eps1) EXIT
156 m = i
157 END DO
158 m = max(m, mbl(iatom)%ma)
159 mbl(iatom)%ma = m
160 mao_blk(iatom) = m
161 END IF
162 END DO
163 CALL group%sum(mao_blk)
164 CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist)
165 CALL dbcsr_release(mao_coef)
166 CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, &
167 matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, &
168 col_blk_size=mao_blk, nze=0)
169 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef)
170 DEALLOCATE (mao_blk, row_blk)
171 !
172 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
173 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
174 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
175 cpassert(iatom == jatom)
176 n = SIZE(cblock, 1)
177 m = SIZE(cblock, 2)
178 cpassert(n == mbl(iatom)%n .AND. m == mbl(iatom)%ma)
179 cblock(1:n, 1:m) = mbl(iatom)%mat(1:n, 1:m)
180 END DO
181 CALL dbcsr_iterator_stop(dbcsr_iter)
182 !
183 END IF
184
185 IF (iolevel > 2) THEN
186 CALL dbcsr_get_info(mao_coef, col_blk_size=col_blk_sizes, &
187 row_blk_size=row_blk_sizes, group=group_handle)
188 CALL group%set_handle(group_handle)
189 DO iatom = 1, nblk
190 n = row_blk_sizes(iatom)
191 m = col_blk_sizes(iatom)
192 ALLOCATE (w(n))
193 w(1:n) = 0._dp
194 IF (ASSOCIATED(mbl(iatom)%mat)) THEN
195 w(1:n) = mbl(iatom)%eig(1:n)
196 END IF
197 CALL group%sum(w)
198 IF (iw > 0) THEN
199 WRITE (iw, '(A,i2,20F8.4)', advance="NO") " Spectrum/Gap ", iatom, w(1:m)
200 WRITE (iw, '(A,F8.4)') " || ", w(m + 1)
201 END IF
202 DEALLOCATE (w)
203 END DO
204 END IF
205
206 CALL mao_orthogonalization(mao_coef, smat)
207
208 DO i = 1, nblk
209 IF (ASSOCIATED(mbl(i)%mat)) THEN
210 DEALLOCATE (mbl(i)%mat)
211 END IF
212 IF (ASSOCIATED(mbl(i)%eig)) THEN
213 DEALLOCATE (mbl(i)%eig)
214 END IF
215 END DO
216 DEALLOCATE (mbl)
217
218 END SUBROUTINE mao_initialization
219
220! **************************************************************************************************
221!> \brief ...
222!> \param mao_coef ...
223!> \param fval ...
224!> \param qmat ...
225!> \param smat ...
226!> \param binv ...
227!> \param reuse ...
228! **************************************************************************************************
229 SUBROUTINE mao_function(mao_coef, fval, qmat, smat, binv, reuse)
230 TYPE(dbcsr_type) :: mao_coef
231 REAL(kind=dp), INTENT(OUT) :: fval
232 TYPE(dbcsr_type) :: qmat, smat, binv
233 LOGICAL, INTENT(IN) :: reuse
234
235 REAL(kind=dp) :: convergence, threshold
236 TYPE(dbcsr_type) :: bmat, scmat, tmat
237
238 threshold = 1.e-8_dp
239 convergence = 1.e-6_dp
240 ! temp matrices
241 CALL dbcsr_create(scmat, template=mao_coef)
242 CALL dbcsr_create(bmat, template=binv)
243 CALL dbcsr_create(tmat, template=qmat)
244 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
245 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
246 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
247 ! calculate inverse of B
248 CALL invert_hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
249 norm_convergence=convergence, silent=.true.)
250 ! calculate Binv*C and T=C(T)*Binv*C
251 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
252 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
253 ! function = Tr(Q*T)
254 CALL dbcsr_dot(qmat, tmat, fval)
255 ! free temp matrices
256 CALL dbcsr_release(scmat)
257 CALL dbcsr_release(bmat)
258 CALL dbcsr_release(tmat)
259
260 END SUBROUTINE mao_function
261
262! **************************************************************************************************
263!> \brief ...
264!> \param mao_coef ...
265!> \param fval ...
266!> \param mao_grad ...
267!> \param qmat ...
268!> \param smat ...
269!> \param binv ...
270!> \param reuse ...
271! **************************************************************************************************
272 SUBROUTINE mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
273 TYPE(dbcsr_type) :: mao_coef
274 REAL(kind=dp), INTENT(OUT) :: fval
275 TYPE(dbcsr_type) :: mao_grad, qmat, smat, binv
276 LOGICAL, INTENT(IN) :: reuse
277
278 REAL(kind=dp) :: convergence, threshold
279 TYPE(dbcsr_type) :: bmat, scmat, t2mat, tmat
280
281 threshold = 1.e-8_dp
282 convergence = 1.e-6_dp
283 ! temp matrices
284 CALL dbcsr_create(scmat, template=mao_coef)
285 CALL dbcsr_create(bmat, template=binv)
286 CALL dbcsr_create(tmat, template=qmat)
287 CALL dbcsr_create(t2mat, template=scmat)
288 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
289 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
290 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
291 ! calculate inverse of B
292 CALL invert_hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
293 norm_convergence=convergence, silent=.true.)
294 ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T)
295 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
296 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
297 ! function = Tr(Q*T)
298 CALL dbcsr_dot(qmat, tmat, fval)
299 ! Gradient part 1: g = 2*Q*C*Binv = 2*Q*R
300 CALL dbcsr_multiply("N", "N", 2.0_dp, qmat, scmat, 0.0_dp, mao_grad, &
301 retain_sparsity=.true.)
302 ! Gradient part 2: g = -2*S*T*X; X = Q*R
303 CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat)
304 CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, t2mat, 0.0_dp, scmat)
305 CALL dbcsr_multiply("N", "N", -2.0_dp, smat, scmat, 1.0_dp, mao_grad, &
306 retain_sparsity=.true.)
307 ! free temp matrices
308 CALL dbcsr_release(scmat)
309 CALL dbcsr_release(bmat)
310 CALL dbcsr_release(tmat)
311 CALL dbcsr_release(t2mat)
312
313 CALL mao_project_gradient(mao_coef, mao_grad, smat)
314
315 END SUBROUTINE mao_function_gradient
316
317! **************************************************************************************************
318!> \brief ...
319!> \param mao_coef ...
320!> \param smat ...
321! **************************************************************************************************
322 SUBROUTINE mao_orthogonalization(mao_coef, smat)
323 TYPE(dbcsr_type) :: mao_coef, smat
324
325 INTEGER :: i, iatom, info, jatom, lwork, m, n
326 LOGICAL :: found
327 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
328 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
329 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, sblock
330 TYPE(dbcsr_iterator_type) :: dbcsr_iter
331
332 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
333 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
334 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
335 cpassert(iatom == jatom)
336 m = SIZE(cblock, 2)
337 n = SIZE(cblock, 1)
338 NULLIFY (sblock)
339 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
340 cpassert(found)
341 lwork = max(n*n, 100)
342 ALLOCATE (amat(n, m), bmat(m, m), w(m), work(lwork))
343 amat(1:n, 1:m) = matmul(sblock(1:n, 1:n), cblock(1:n, 1:m))
344 bmat(1:m, 1:m) = matmul(transpose(cblock(1:n, 1:m)), amat(1:n, 1:m))
345 info = 0
346 CALL lapack_ssyev("V", "U", m, bmat, m, w, work, lwork, info)
347 cpassert(info == 0)
348 cpassert(all(w > 0.0_dp))
349 w = 1.0_dp/sqrt(w)
350 DO i = 1, m
351 amat(1:m, i) = bmat(1:m, i)*w(i)
352 END DO
353 bmat(1:m, 1:m) = matmul(amat(1:m, 1:m), transpose(bmat(1:m, 1:m)))
354 cblock(1:n, 1:m) = matmul(cblock(1:n, 1:m), bmat(1:m, 1:m))
355 DEALLOCATE (amat, bmat, w, work)
356 END DO
357 CALL dbcsr_iterator_stop(dbcsr_iter)
358
359 END SUBROUTINE mao_orthogonalization
360
361! **************************************************************************************************
362!> \brief ...
363!> \param mao_coef ...
364!> \param mao_grad ...
365!> \param smat ...
366! **************************************************************************************************
367 SUBROUTINE mao_project_gradient(mao_coef, mao_grad, smat)
368 TYPE(dbcsr_type) :: mao_coef, mao_grad, smat
369
370 INTEGER :: iatom, jatom, m, n
371 LOGICAL :: found
372 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
373 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
374 TYPE(dbcsr_iterator_type) :: dbcsr_iter
375
376 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
377 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
378 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
379 cpassert(iatom == jatom)
380 m = SIZE(cblock, 2)
381 n = SIZE(cblock, 1)
382 NULLIFY (sblock)
383 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
384 cpassert(found)
385 NULLIFY (gblock)
386 CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found)
387 cpassert(found)
388 ALLOCATE (amat(m, m))
389 amat(1:m, 1:m) = matmul(transpose(cblock(1:n, 1:m)), matmul(sblock(1:n, 1:n), gblock(1:n, 1:m)))
390 gblock(1:n, 1:m) = gblock(1:n, 1:m) - matmul(cblock(1:n, 1:m), amat(1:m, 1:m))
391 DEALLOCATE (amat)
392 END DO
393 CALL dbcsr_iterator_stop(dbcsr_iter)
394
395 END SUBROUTINE mao_project_gradient
396
397! **************************************************************************************************
398!> \brief ...
399!> \param fmat1 ...
400!> \param fmat2 ...
401!> \return ...
402! **************************************************************************************************
403 FUNCTION mao_scalar_product(fmat1, fmat2) RESULT(spro)
404 TYPE(dbcsr_type) :: fmat1, fmat2
405 REAL(kind=dp) :: spro
406
407 INTEGER :: group_handle, iatom, jatom, m, n
408 LOGICAL :: found
409 REAL(kind=dp), DIMENSION(:, :), POINTER :: ablock, bblock
410 TYPE(dbcsr_iterator_type) :: dbcsr_iter
411 TYPE(mp_comm_type) :: group
412
413 spro = 0.0_dp
414
415 CALL dbcsr_iterator_start(dbcsr_iter, fmat1)
416 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
417 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock)
418 cpassert(iatom == jatom)
419 m = SIZE(ablock, 2)
420 n = SIZE(ablock, 1)
421 CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found)
422 cpassert(found)
423 spro = spro + sum(ablock(1:n, 1:m)*bblock(1:n, 1:m))
424 END DO
425 CALL dbcsr_iterator_stop(dbcsr_iter)
426
427 CALL dbcsr_get_info(fmat1, group=group_handle)
428 CALL group%set_handle(group_handle)
429 CALL group%sum(spro)
430
431 END FUNCTION mao_scalar_product
432
433! **************************************************************************************************
434!> \brief Calculate the density matrix at the Gamma point
435!> \param pmat ...
436!> \param ksmat ...
437!> \param smat ...
438!> \param kpoints Kpoint environment
439!> \param nmos Number of occupied orbitals
440!> \param occ Maximum occupation per orbital
441!> \par History
442!> 04.2016 created [JGH]
443! **************************************************************************************************
444 SUBROUTINE calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
445
446 TYPE(dbcsr_type) :: pmat, ksmat, smat
447 TYPE(kpoint_type), POINTER :: kpoints
448 INTEGER, INTENT(IN) :: nmos
449 REAL(kind=dp), INTENT(IN) :: occ
450
451 INTEGER :: norb
452 REAL(kind=dp) :: de
453 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
454 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
455 TYPE(cp_fm_type) :: fmksmat, fmsmat, fmvec, fmwork
456 TYPE(dbcsr_type) :: tempmat
457
458 ! FM matrices
459
460 CALL dbcsr_get_info(smat, nfullrows_total=norb)
461 CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, &
462 nrow_global=norb, ncol_global=norb)
463 CALL cp_fm_create(fmksmat, matrix_struct)
464 CALL cp_fm_create(fmsmat, matrix_struct)
465 CALL cp_fm_create(fmvec, matrix_struct)
466 CALL cp_fm_create(fmwork, matrix_struct)
467 ALLOCATE (eigenvalues(norb))
468
469 ! DBCSR matrix
470 CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
471
472 ! transfer to FM
473 CALL dbcsr_desymmetrize(smat, tempmat)
474 CALL copy_dbcsr_to_fm(tempmat, fmsmat)
475 CALL dbcsr_desymmetrize(ksmat, tempmat)
476 CALL copy_dbcsr_to_fm(tempmat, fmksmat)
477
478 ! diagonalize
479 CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
480 de = eigenvalues(nmos + 1) - eigenvalues(nmos)
481 IF (de < 0.001_dp) THEN
482 CALL cp_warn(__location__, "MAO: No band gap at "// &
483 "Gamma point. MAO analysis not reliable.")
484 END IF
485 ! density matrix
486 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ)
487
488 DEALLOCATE (eigenvalues)
489 CALL dbcsr_release(tempmat)
490 CALL cp_fm_release(fmksmat)
491 CALL cp_fm_release(fmsmat)
492 CALL cp_fm_release(fmvec)
493 CALL cp_fm_release(fmwork)
494 CALL cp_fm_struct_release(matrix_struct)
495
496 END SUBROUTINE calculate_p_gamma
497
498! **************************************************************************************************
499!> \brief Define the MAO reference basis set
500!> \param qs_env ...
501!> \param mao_basis ...
502!> \param mao_basis_set_list ...
503!> \param orb_basis_set_list ...
504!> \param iunit ...
505!> \param print_basis ...
506!> \par History
507!> 07.2016 created [JGH]
508! **************************************************************************************************
509 SUBROUTINE mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, &
510 iunit, print_basis)
511
512 TYPE(qs_environment_type), POINTER :: qs_env
513 INTEGER, INTENT(IN) :: mao_basis
514 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
515 INTEGER, INTENT(IN), OPTIONAL :: iunit
516 LOGICAL, INTENT(IN), OPTIONAL :: print_basis
517
518 INTEGER :: ikind, nbas, nkind, unit_nr
519 REAL(kind=dp) :: eps_pgf_orb
520 TYPE(dft_control_type), POINTER :: dft_control
521 TYPE(gto_basis_set_type), POINTER :: basis_set, pbasis
522 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
523 TYPE(qs_kind_type), POINTER :: qs_kind
524
525 ! Reference basis set
526 cpassert(.NOT. ASSOCIATED(mao_basis_set_list))
527 cpassert(.NOT. ASSOCIATED(orb_basis_set_list))
528
529 ! options
530 IF (PRESENT(iunit)) THEN
531 unit_nr = iunit
532 ELSE
533 unit_nr = -1
534 END IF
535
536 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
537 nkind = SIZE(qs_kind_set)
538 ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
539 DO ikind = 1, nkind
540 NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
541 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
542 END DO
543 !
544 DO ikind = 1, nkind
545 qs_kind => qs_kind_set(ikind)
546 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
547 IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set
548 END DO
549 !
550 SELECT CASE (mao_basis)
551 CASE (mao_basis_orb)
552 DO ikind = 1, nkind
553 qs_kind => qs_kind_set(ikind)
554 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
555 IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set
556 END DO
557 CASE (mao_basis_prim)
558 DO ikind = 1, nkind
559 qs_kind => qs_kind_set(ikind)
560 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
561 NULLIFY (pbasis)
562 IF (ASSOCIATED(basis_set)) THEN
563 CALL create_primitive_basis_set(basis_set, pbasis)
564 CALL get_qs_env(qs_env, dft_control=dft_control)
565 eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
566 CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb)
567 pbasis%kind_radius = basis_set%kind_radius
568 mao_basis_set_list(ikind)%gto_basis_set => pbasis
569 CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO")
570 END IF
571 END DO
572 CASE (mao_basis_ext)
573 DO ikind = 1, nkind
574 qs_kind => qs_kind_set(ikind)
575 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO")
576 IF (ASSOCIATED(basis_set)) THEN
577 basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
578 mao_basis_set_list(ikind)%gto_basis_set => basis_set
579 END IF
580 END DO
581 CASE DEFAULT
582 cpabort("Unknown option for MAO basis")
583 END SELECT
584 IF (unit_nr > 0) THEN
585 DO ikind = 1, nkind
586 IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
587 WRITE (unit=unit_nr, fmt="(T2,A,I4)") &
588 "WARNING: No MAO basis set associated with Kind ", ikind
589 ELSE
590 nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
591 WRITE (unit=unit_nr, fmt="(T2,A,I4,T56,A,I10)") &
592 "MAO basis set Kind ", ikind, " Number of BSF:", nbas
593 END IF
594 END DO
595 END IF
596
597 IF (PRESENT(print_basis)) THEN
598 IF (print_basis) THEN
599 DO ikind = 1, nkind
600 basis_set => mao_basis_set_list(ikind)%gto_basis_set
601 IF (ASSOCIATED(basis_set)) CALL write_gto_basis_set(basis_set, unit_nr, "MAO REFERENCE BASIS")
602 END DO
603 END IF
604 END IF
605
606 END SUBROUTINE mao_reference_basis
607
608! **************************************************************************************************
609!> \brief Analyze the MAO basis, projection on angular functions
610!> \param mao_coef ...
611!> \param matrix_smm ...
612!> \param mao_basis_set_list ...
613!> \param particle_set ...
614!> \param qs_kind_set ...
615!> \param unit_nr ...
616!> \param para_env ...
617!> \par History
618!> 07.2016 created [JGH]
619! **************************************************************************************************
620 SUBROUTINE mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
621 qs_kind_set, unit_nr, para_env)
622
623 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_smm
624 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list
625 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
626 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
627 INTEGER, INTENT(IN) :: unit_nr
628 TYPE(mp_para_env_type), POINTER :: para_env
629
630 CHARACTER(len=2) :: element_symbol
631 INTEGER :: ia, iab, iatom, ikind, iset, ishell, &
632 ispin, l, lmax, lshell, m, ma, na, &
633 natom, nspin
634 LOGICAL :: found
635 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cmask, vec1, vec2
636 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: weight
637 REAL(kind=dp), DIMENSION(:, :), POINTER :: block, cmao
638 TYPE(gto_basis_set_type), POINTER :: basis_set
639
640 ! Analyze the MAO basis
641 IF (unit_nr > 0) THEN
642 WRITE (unit_nr, "(/,A)") " Analyze angular momentum character of MAOs "
643 WRITE (unit_nr, "(T7,A,T15,A,T20,A,T40,A,T50,A,T60,A,T70,A,T80,A)") &
644 "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G"
645 END IF
646 lmax = 4 ! analyze up to g-functions
647 natom = SIZE(particle_set)
648 nspin = SIZE(mao_coef)
649 DO iatom = 1, natom
650 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
651 element_symbol=element_symbol, kind_number=ikind)
652 basis_set => mao_basis_set_list(ikind)%gto_basis_set
653 CALL get_qs_kind(qs_kind_set(ikind), mao=na)
654 CALL get_gto_basis_set(basis_set, nsgf=ma)
655 ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na))
656 weight = 0.0_dp
657 CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, &
658 block=block, found=found)
659 DO ispin = 1, nspin
660 CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, &
661 block=cmao, found=found)
662 IF (found) THEN
663 DO l = 0, lmax
664 cmask = 0.0_dp
665 iab = 0
666 DO iset = 1, basis_set%nset
667 DO ishell = 1, basis_set%nshell(iset)
668 lshell = basis_set%l(ishell, iset)
669 DO m = -lshell, lshell
670 iab = iab + 1
671 IF (l == lshell) cmask(iab) = 1.0_dp
672 END DO
673 END DO
674 END DO
675 DO ia = 1, na
676 vec1(1:ma) = cmask*cmao(1:ma, ia)
677 vec2(1:ma) = matmul(block, vec1)
678 weight(l, ia) = sum(vec1(1:ma)*vec2(1:ma))
679 END DO
680 END DO
681 END IF
682 CALL para_env%sum(weight)
683 IF (unit_nr > 0) THEN
684 DO ia = 1, na
685 IF (ispin == 1 .AND. ia == 1) THEN
686 WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") &
687 iatom, element_symbol, ispin, ia, weight(0:lmax, ia)
688 ELSE
689 WRITE (unit_nr, "(T17,i2,T20,i3,T31,5F10.4)") ispin, ia, weight(0:lmax, ia)
690 END IF
691 END DO
692 END IF
693 END DO
694 DEALLOCATE (cmask, weight, vec1, vec2)
695 END DO
696 END SUBROUTINE mao_basis_analysis
697
698! **************************************************************************************************
699!> \brief Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap
700!> \param matrix_q ...
701!> \param matrix_p ...
702!> \param matrix_s ...
703!> \param matrix_smm ...
704!> \param matrix_smo ...
705!> \param smm_list ...
706!> \param electra ...
707!> \param eps_filter ...
708!> \param nimages ...
709!> \param kpoints ...
710!> \param matrix_ks ...
711!> \param sab_orb ...
712!> \par History
713!> 08.2016 created [JGH]
714! **************************************************************************************************
715 SUBROUTINE mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, &
716 electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
717
718 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q
719 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
720 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_smm, matrix_smo
721 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
722 POINTER :: smm_list
723 REAL(kind=dp), DIMENSION(2), INTENT(OUT) :: electra
724 REAL(kind=dp), INTENT(IN) :: eps_filter
725 INTEGER, INTENT(IN), OPTIONAL :: nimages
726 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
727 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
728 POINTER :: matrix_ks
729 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
730 OPTIONAL, POINTER :: sab_orb
731
732 INTEGER :: im, ispin, nim, nocc, norb, nspin
733 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
734 REAL(kind=dp) :: elex, xkp(3)
735 TYPE(dbcsr_type) :: ksmat, pmat, smat, tmat
736
737 nim = 1
738 IF (PRESENT(nimages)) nim = nimages
739 IF (nim > 1) THEN
740 cpassert(PRESENT(kpoints))
741 cpassert(PRESENT(matrix_ks))
742 cpassert(PRESENT(sab_orb))
743 END IF
744
745 ! Reference
746 nspin = SIZE(matrix_p, 1)
747 DO ispin = 1, nspin
748 electra(ispin) = 0.0_dp
749 DO im = 1, nim
750 CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex)
751 electra(ispin) = electra(ispin) + elex
752 END DO
753 END DO
754
755 ! Q matrix
756 NULLIFY (matrix_q)
757 CALL dbcsr_allocate_matrix_set(matrix_q, nspin)
758 DO ispin = 1, nspin
759 ALLOCATE (matrix_q(ispin)%matrix)
760 CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix)
761 CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list)
762 END DO
763 ! temp matrix
764 CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
765 ! Q=APA(T)
766 DO ispin = 1, nspin
767 IF (nim == 1) THEN
768 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, &
769 0.0_dp, tmat, filter_eps=eps_filter)
770 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
771 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
772 ELSE
773 ! k-points
774 CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix)
775 CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix)
776 CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix)
777 CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb)
778 CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb)
779 CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb)
780 NULLIFY (cell_to_index)
781 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
782 ! calculate density matrix at gamma point
783 xkp = 0.0_dp
784 ! transform KS and S matrices to the gamma point
785 CALL dbcsr_set(ksmat, 0.0_dp)
786 CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, &
787 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
788 CALL dbcsr_set(smat, 0.0_dp)
789 CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, &
790 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
791 norb = nint(electra(ispin))
792 nocc = mod(2, nspin) + 1
793 CALL calculate_p_gamma(pmat, ksmat, smat, kpoints, norb, real(nocc, kind=dp))
794 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, pmat, &
795 0.0_dp, tmat, filter_eps=eps_filter)
796 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
797 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
798 CALL dbcsr_release(pmat)
799 CALL dbcsr_release(smat)
800 CALL dbcsr_release(ksmat)
801 END IF
802 END DO
803 ! free temp matrix
804 CALL dbcsr_release(tmat)
805
806 END SUBROUTINE mao_build_q
807
808END MODULE mao_methods
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
subroutine, public create_primitive_basis_set(basis_set, pbasis)
...
subroutine, public write_gto_basis_set(gto_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
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_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
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_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public mao_basis_orb
integer, parameter, public mao_basis_ext
integer, parameter, public mao_basis_prim
Routines useful for iterative matrix calculations.
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Routines needed for kpoint calculation.
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Interface to the LAPACK F77 library.
Definition lapack.F:17
Calculate MAO's and analyze wavefunctions.
Definition mao_basis.F:15
Calculate MAO's and analyze wavefunctions.
Definition mao_methods.F:15
subroutine, public mao_project_gradient(mao_coef, mao_grad, smat)
...
subroutine, public mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
...
subroutine, public mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, iunit, print_basis)
Define the MAO reference basis set.
subroutine, public mao_orthogonalization(mao_coef, smat)
...
subroutine, public calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
Calculate the density matrix at the Gamma point.
subroutine, public mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
...
Definition mao_methods.F:91
subroutine, public mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, qs_kind_set, unit_nr, para_env)
Analyze the MAO basis, projection on angular functions.
real(kind=dp) function, public mao_scalar_product(fmat1, fmat2)
...
subroutine, public mao_function(mao_coef, fval, qmat, smat, binv, reuse)
...
subroutine, public mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap.
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Define the neighbor list data types and the corresponding functionality.
keeps the information about the structure of a full matrix
represent a full matrix
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.