(git:1f9fd2c)
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-2026 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! **************************************************************************************************
15MODULE mao_methods
16 USE atomic_kind_types, ONLY: get_atomic_kind
17 USE basis_set_container_types, ONLY: add_basis_set_to_container
18 USE basis_set_types, ONLY: create_primitive_basis_set,&
19 get_gto_basis_set,&
20 gto_basis_set_p_type,&
21 gto_basis_set_type,&
22 write_gto_basis_set
23 USE cp_control_types, ONLY: dft_control_type
24 USE cp_dbcsr_api, ONLY: &
25 dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_get_block_p, &
26 dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
27 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
28 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
29 dbcsr_type_symmetric
30 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
31 dbcsr_reserve_diag_blocks
32 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 cp_dbcsr_plus_fm_fm_t,&
35 dbcsr_allocate_matrix_set
36 USE cp_fm_diag, ONLY: cp_fm_geeig
37 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
38 cp_fm_struct_release,&
39 cp_fm_struct_type
40 USE cp_fm_types, ONLY: cp_fm_create,&
41 cp_fm_release,&
42 cp_fm_type
43 USE input_constants, ONLY: mao_basis_ext,&
44 mao_basis_orb,&
45 mao_basis_prim
46 USE iterate_matrix, ONLY: invert_hotelling
47 USE kinds, ONLY: dp
48 USE kpoint_methods, ONLY: rskp_transform
49 USE kpoint_types, ONLY: get_kpoint_info,&
50 kpoint_type
51 USE mathlib, ONLY: invert_matrix
52 USE message_passing, ONLY: mp_comm_type,&
53 mp_para_env_type
54 USE particle_types, ONLY: particle_type
55 USE qs_environment_types, ONLY: get_qs_env,&
56 qs_environment_type
57 USE qs_interactions, ONLY: init_interaction_radii_orb_basis
58 USE qs_kind_types, ONLY: get_qs_kind,&
59 qs_kind_type
60 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
61#include "./base/base_uses.f90"
62
63 IMPLICIT NONE
64 PRIVATE
65
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_methods'
67
68 TYPE mblocks
69 INTEGER :: n = -1, ma = -1
70 REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat => null()
71 REAL(KIND=dp), DIMENSION(:), POINTER :: eig => null()
72 END TYPE mblocks
73
74 PUBLIC :: mao_initialization, mao_function, mao_function_gradient, mao_orthogonalization, &
75 mao_project_gradient, mao_scalar_product, mao_build_q, mao_basis_analysis, &
76 mao_reference_basis, calculate_p_gamma, mao_update_coef, mao_preconditioner
77
78! **************************************************************************************************
79
80CONTAINS
81
82! **************************************************************************************************
83!> \brief ...
84!> \param mao_coef ...
85!> \param pmat ...
86!> \param smat ...
87!> \param eps1 ...
88!> \param iolevel ...
89!> \param iw ...
90! **************************************************************************************************
91 SUBROUTINE mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
92 TYPE(dbcsr_type) :: mao_coef, pmat, smat
93 REAL(kind=dp), INTENT(IN) :: eps1
94 INTEGER, INTENT(IN) :: iolevel, iw
95
96 INTEGER :: i, iatom, info, jatom, 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 dsygv(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)
146 ALLOCATE (row_blk(nblk), mao_blk(nblk))
147 mao_blk = 0
148 row_blk = row_blk_sizes
149 DO iatom = 1, nblk
150 IF (ASSOCIATED(mbl(iatom)%mat)) THEN
151 n = mbl(iatom)%n
152 m = 0
153 DO i = 1, n
154 IF (mbl(iatom)%eig(i) < eps1) EXIT
155 m = i
156 END DO
157 m = max(m, mbl(iatom)%ma)
158 mbl(iatom)%ma = m
159 mao_blk(iatom) = m
160 END IF
161 END DO
162 CALL group%sum(mao_blk)
163 CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist)
164 CALL dbcsr_release(mao_coef)
165 CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, &
166 matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, col_blk_size=mao_blk)
167 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef)
168 DEALLOCATE (mao_blk, row_blk)
169 !
170 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
171 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
172 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
173 cpassert(iatom == jatom)
174 n = SIZE(cblock, 1)
175 m = SIZE(cblock, 2)
176 cpassert(n == mbl(iatom)%n .AND. m == mbl(iatom)%ma)
177 cblock(1:n, 1:m) = mbl(iatom)%mat(1:n, 1:m)
178 END DO
179 CALL dbcsr_iterator_stop(dbcsr_iter)
180 !
181 END IF
182
183 IF (iolevel > 2) THEN
184 CALL dbcsr_get_info(mao_coef, col_blk_size=col_blk_sizes, &
185 row_blk_size=row_blk_sizes, group=group)
186 DO iatom = 1, nblk
187 n = row_blk_sizes(iatom)
188 m = col_blk_sizes(iatom)
189 ALLOCATE (w(n))
190 w(1:n) = 0._dp
191 IF (ASSOCIATED(mbl(iatom)%mat)) THEN
192 w(1:n) = mbl(iatom)%eig(1:n)
193 END IF
194 CALL group%sum(w)
195 IF (iw > 0) THEN
196 WRITE (iw, '(A,i2,20F8.4)', advance="NO") " Spectrum/Gap ", iatom, w(1:m)
197 WRITE (iw, '(A,F8.4)') " || ", w(m + 1)
198 END IF
199 DEALLOCATE (w)
200 END DO
201 END IF
202
203 CALL mao_orthogonalization(mao_coef, smat)
204
205 DO i = 1, nblk
206 IF (ASSOCIATED(mbl(i)%mat)) THEN
207 DEALLOCATE (mbl(i)%mat)
208 END IF
209 IF (ASSOCIATED(mbl(i)%eig)) THEN
210 DEALLOCATE (mbl(i)%eig)
211 END IF
212 END DO
213 DEALLOCATE (mbl)
214
215 END SUBROUTINE mao_initialization
216
217! **************************************************************************************************
218!> \brief ...
219!> \param mao_coef ...
220!> \param fval ...
221!> \param qmat ...
222!> \param smat ...
223!> \param binv ...
224!> \param reuse ...
225! **************************************************************************************************
226 SUBROUTINE mao_function(mao_coef, fval, qmat, smat, binv, reuse)
227 TYPE(dbcsr_type) :: mao_coef
228 REAL(kind=dp), INTENT(OUT) :: fval
229 TYPE(dbcsr_type) :: qmat, smat, binv
230 LOGICAL, INTENT(IN) :: reuse
231
232 REAL(kind=dp) :: convergence, threshold
233 TYPE(dbcsr_type) :: bmat, scmat, tmat
234
235 threshold = 1.e-8_dp
236 convergence = 1.e-6_dp
237 ! temp matrices
238 CALL dbcsr_create(scmat, template=mao_coef)
239 CALL dbcsr_create(bmat, template=binv)
240 CALL dbcsr_create(tmat, template=qmat)
241 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
242 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
243 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
244 ! calculate inverse of B
245 CALL invert_hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
246 norm_convergence=convergence, silent=.true.)
247 ! calculate Binv*C and T=C(T)*Binv*C
248 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
249 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
250 ! function = Tr(Q*T)
251 CALL dbcsr_dot(qmat, tmat, fval)
252 ! free temp matrices
253 CALL dbcsr_release(scmat)
254 CALL dbcsr_release(bmat)
255 CALL dbcsr_release(tmat)
256
257 END SUBROUTINE mao_function
258
259! **************************************************************************************************
260!> \brief ...
261!> \param mao_coef ...
262!> \param fval ...
263!> \param mao_grad ...
264!> \param qmat ...
265!> \param smat ...
266!> \param binv ...
267!> \param reuse ...
268! **************************************************************************************************
269 SUBROUTINE mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
270 TYPE(dbcsr_type) :: mao_coef
271 REAL(kind=dp), INTENT(OUT) :: fval
272 TYPE(dbcsr_type) :: mao_grad, qmat, smat, binv
273 LOGICAL, INTENT(IN) :: reuse
274
275 REAL(kind=dp) :: convergence, threshold
276 TYPE(dbcsr_type) :: bmat, scmat, t2mat, tmat
277
278 threshold = 1.e-8_dp
279 convergence = 1.e-6_dp
280 ! temp matrices
281 CALL dbcsr_create(scmat, template=mao_coef)
282 CALL dbcsr_create(bmat, template=binv)
283 CALL dbcsr_create(tmat, template=qmat)
284 CALL dbcsr_create(t2mat, template=scmat)
285 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
286 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
287 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
288 ! calculate inverse of B
289 CALL invert_hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
290 norm_convergence=convergence, silent=.true.)
291 ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T)
292 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
293 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
294 ! function = Tr(Q*T)
295 CALL dbcsr_dot(qmat, tmat, fval)
296 ! Gradient part 1: g = 2*Q*C*Binv = 2*Q*R
297 CALL dbcsr_multiply("N", "N", 2.0_dp, qmat, scmat, 0.0_dp, mao_grad, &
298 retain_sparsity=.true.)
299 ! Gradient part 2: g = -2*S*T*X; X = Q*R
300 CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat)
301 CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, t2mat, 0.0_dp, scmat)
302 CALL dbcsr_multiply("N", "N", -2.0_dp, smat, scmat, 1.0_dp, mao_grad, &
303 retain_sparsity=.true.)
304 ! free temp matrices
305 CALL dbcsr_release(scmat)
306 CALL dbcsr_release(bmat)
307 CALL dbcsr_release(tmat)
308 CALL dbcsr_release(t2mat)
309
310 ! apply Lagrange multiplyer
311 CALL mao_grad_lagrange(mao_coef, mao_grad, smat)
312
313 END SUBROUTINE mao_function_gradient
314
315! **************************************************************************************************
316!> \brief Preconditioner for MAO optimization (seems not to work, but left for further tests)
317!> \param mao_coef ...
318!> \param mao_grad ...
319!> \param qmat ...
320!> \param smat ...
321! **************************************************************************************************
322 SUBROUTINE mao_preconditioner(mao_coef, mao_grad, qmat, smat)
323 TYPE(dbcsr_type) :: mao_coef, mao_grad, qmat, smat
324
325 INTEGER :: i, iatom, jatom, m, n
326 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes
327 LOGICAL :: found
328 REAL(kind=dp) :: eval_error, q_shift
329 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qinv, qsmat
330 REAL(kind=dp), DIMENSION(:, :), POINTER :: bblock, gblock, qblock
331 TYPE(dbcsr_distribution_type) :: dbcsr_dist
332 TYPE(dbcsr_iterator_type) :: dbcsr_iter
333 TYPE(dbcsr_type) :: bmat, scmat
334
335 ! temp matrices
336 CALL dbcsr_create(scmat, template=mao_coef)
337 CALL dbcsr_get_info(mao_coef, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
338 CALL dbcsr_create(bmat, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
339 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
340 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
341 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
342 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
343
344 q_shift = 0.001_dp
345
346 CALL dbcsr_iterator_start(dbcsr_iter, mao_grad)
347 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
348 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, gblock)
349 cpassert(iatom == jatom)
350 m = SIZE(gblock, 1)
351 n = SIZE(gblock, 2)
352 CALL dbcsr_get_block_p(matrix=bmat, row=iatom, col=jatom, block=bblock, found=found)
353 cpassert(found)
354 CALL dbcsr_get_block_p(matrix=qmat, row=iatom, col=jatom, block=qblock, found=found)
355 cpassert(found)
356 ALLOCATE (qinv(m, m), qsmat(m, m))
357 qsmat(1:m, 1:m) = qblock(1:m, 1:m)
358 DO i = 1, m
359 qsmat(i, i) = qsmat(i, i) + q_shift
360 END DO
361 CALL invert_matrix(qsmat, qinv, eval_error)
362 ! g = q * g * b
363 gblock(1:m, 1:n) = matmul(qinv(1:m, 1:m), matmul(gblock(1:m, 1:n), bblock(1:n, 1:n)))
364 DEALLOCATE (qinv, qsmat)
365 END DO
366 CALL dbcsr_iterator_stop(dbcsr_iter)
367
368 CALL mao_project_gradient(mao_coef, mao_grad, smat)
369
370 ! free temp matrices
371 CALL dbcsr_release(scmat)
372 CALL dbcsr_release(bmat)
373
374 END SUBROUTINE mao_preconditioner
375
376! **************************************************************************************************
377!> \brief ...
378!> \param mao_coef ...
379!> \param mao_grad ...
380!> \param smat ...
381! **************************************************************************************************
382 SUBROUTINE mao_grad_lagrange(mao_coef, mao_grad, smat)
383 TYPE(dbcsr_type) :: mao_coef, mao_grad, smat
384
385 INTEGER :: i, iatom, info, jatom, lwork, m, n
386 LOGICAL :: found
387 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
388 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, lmat, rmat
389 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
390 TYPE(dbcsr_iterator_type) :: dbcsr_iter
391
392 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
393 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
394 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
395 cpassert(iatom == jatom)
396 m = SIZE(cblock, 2)
397 n = SIZE(cblock, 1)
398 NULLIFY (sblock)
399 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
400 cpassert(found)
401 lwork = max(n*n, 100)
402 ALLOCATE (amat(n, m), bmat(m, m), rmat(m, m), lmat(n, n), w(m), work(lwork))
403 amat(1:n, 1:m) = matmul(sblock(1:n, 1:n), cblock(1:n, 1:m))
404 bmat(1:m, 1:m) = matmul(transpose(amat(1:n, 1:m)), amat(1:n, 1:m))
405 info = 0
406 CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
407 cpassert(info == 0)
408 cpassert(all(w > 0.0_dp))
409 w = 1.0_dp/w
410 DO i = 1, m
411 rmat(1:m, i) = bmat(1:m, i)*w(i)
412 END DO
413 bmat(1:m, 1:m) = matmul(rmat(1:m, 1:m), transpose(bmat(1:m, 1:m)))
414 lmat(1:n, 1:n) = matmul(matmul(amat(1:n, 1:m), bmat(1:m, 1:m)), transpose(amat(1:n, 1:m)))
415 !
416 CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=iatom, block=gblock, found=found)
417 cpassert(found)
418 gblock(1:n, 1:m) = gblock(1:n, 1:m) - matmul(lmat(1:n, 1:n), gblock(1:n, 1:m))
419 DEALLOCATE (amat, bmat, rmat, lmat, w, work)
420 END DO
421 CALL dbcsr_iterator_stop(dbcsr_iter)
422
423 END SUBROUTINE mao_grad_lagrange
424
425! **************************************************************************************************
426!> \brief ...
427!> \param mao_up ...
428!> \param mao_coef ...
429!> \param mao_grad ...
430!> \param smat ...
431!> \param alpha ...
432! **************************************************************************************************
433 SUBROUTINE mao_update_coef(mao_up, mao_coef, mao_grad, smat, alpha)
434 TYPE(dbcsr_type) :: mao_up, mao_coef, mao_grad, smat
435 REAL(kind=dp) :: alpha
436
437 INTEGER :: i, iatom, info, jatom, lwork, m, n
438 LOGICAL :: found
439 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
440 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, r1mat, r2mat
441 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock, ublock
442 TYPE(dbcsr_iterator_type) :: dbcsr_iter
443
444 CALL dbcsr_copy(mao_up, mao_coef)
445 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
446 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
447 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
448 cpassert(iatom == jatom)
449 CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=iatom, block=gblock, found=found)
450 cpassert(found)
451 m = SIZE(cblock, 2)
452 n = SIZE(cblock, 1)
453 NULLIFY (sblock)
454 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
455 cpassert(found)
456 CALL dbcsr_get_block_p(matrix=mao_up, row=iatom, col=iatom, block=ublock, found=found)
457 cpassert(found)
458 lwork = max(n*n, 100)
459 ALLOCATE (amat(n, m), bmat(m, m), r1mat(m, m), r2mat(m, m), w(m), work(lwork))
460 !
461 amat(1:n, 1:m) = alpha*gblock(1:n, 1:m)
462 bmat(1:m, 1:m) = matmul(transpose(amat(1:n, 1:m)), matmul(sblock(1:n, 1:n), amat(1:n, 1:m)))
463 info = 0
464 CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
465 cpassert(info == 0)
466 DO i = 1, m
467 IF (w(i) < 1.e-12_dp) THEN
468 w(i) = 0.0_dp
469 ELSE
470 w(i) = sqrt(w(i))
471 END IF
472 r1mat(1:m, i) = bmat(1:m, i)*cos(w(i))
473 IF (w(i) < 1.e-6_dp) THEN
474 r2mat(1:m, i) = bmat(1:m, i)*(1._dp - w(i)**2/6._dp)
475 ELSE
476 r2mat(1:m, i) = bmat(1:m, i)*sin(w(i))/w(i)
477 END IF
478 END DO
479 r1mat(1:m, 1:m) = matmul(r1mat(1:m, 1:m), transpose(bmat(1:m, 1:m)))
480 r2mat(1:m, 1:m) = matmul(r2mat(1:m, 1:m), transpose(bmat(1:m, 1:m)))
481 !
482 ublock(1:n, 1:m) = matmul(cblock(1:n, 1:m), r1mat(1:m, 1:m)) + &
483 matmul(amat(1:n, 1:m), r2mat(1:m, 1:m))
484 !
485 DEALLOCATE (amat, bmat, r1mat, r2mat, w, work)
486 END DO
487 CALL dbcsr_iterator_stop(dbcsr_iter)
488
489 END SUBROUTINE mao_update_coef
490
491! **************************************************************************************************
492!> \brief ...
493!> \param mao_coef ...
494!> \param smat ...
495! **************************************************************************************************
496 SUBROUTINE mao_orthogonalization(mao_coef, smat)
497 TYPE(dbcsr_type) :: mao_coef, smat
498
499 INTEGER :: i, iatom, info, jatom, lwork, m, n
500 LOGICAL :: found
501 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
502 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, rmat
503 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, sblock
504 TYPE(dbcsr_iterator_type) :: dbcsr_iter
505
506 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
507 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
508 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
509 cpassert(iatom == jatom)
510 m = SIZE(cblock, 2)
511 n = SIZE(cblock, 1)
512 NULLIFY (sblock)
513 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
514 cpassert(found)
515 lwork = max(n*n, 100)
516 ALLOCATE (amat(n, m), bmat(m, m), rmat(m, m), w(m), work(lwork))
517 amat(1:n, 1:m) = matmul(sblock(1:n, 1:n), cblock(1:n, 1:m))
518 bmat(1:m, 1:m) = matmul(transpose(cblock(1:n, 1:m)), amat(1:n, 1:m))
519 info = 0
520 CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
521 cpassert(info == 0)
522 cpassert(all(w > 0.0_dp))
523 w = 1.0_dp/sqrt(w)
524 DO i = 1, m
525 rmat(1:m, i) = bmat(1:m, i)*w(i)
526 END DO
527 bmat(1:m, 1:m) = matmul(rmat(1:m, 1:m), transpose(bmat(1:m, 1:m)))
528 cblock(1:n, 1:m) = matmul(cblock(1:n, 1:m), bmat(1:m, 1:m))
529 DEALLOCATE (amat, bmat, rmat, w, work)
530 END DO
531 CALL dbcsr_iterator_stop(dbcsr_iter)
532
533 END SUBROUTINE mao_orthogonalization
534
535! **************************************************************************************************
536!> \brief ...
537!> \param mao_coef ...
538!> \param mao_grad ...
539!> \param smat ...
540! **************************************************************************************************
541 SUBROUTINE mao_project_gradient(mao_coef, mao_grad, smat)
542 TYPE(dbcsr_type) :: mao_coef, mao_grad, smat
543
544 INTEGER :: iatom, jatom, m, n
545 LOGICAL :: found
546 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
547 REAL(kind=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
548 TYPE(dbcsr_iterator_type) :: dbcsr_iter
549
550 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
551 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
552 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
553 cpassert(iatom == jatom)
554 m = SIZE(cblock, 2)
555 n = SIZE(cblock, 1)
556 NULLIFY (sblock)
557 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
558 cpassert(found)
559 NULLIFY (gblock)
560 CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found)
561 cpassert(found)
562 ALLOCATE (amat(m, m))
563 amat(1:m, 1:m) = matmul(transpose(cblock(1:n, 1:m)), matmul(sblock(1:n, 1:n), gblock(1:n, 1:m)))
564 gblock(1:n, 1:m) = gblock(1:n, 1:m) - matmul(cblock(1:n, 1:m), amat(1:m, 1:m))
565 DEALLOCATE (amat)
566 END DO
567 CALL dbcsr_iterator_stop(dbcsr_iter)
568
569 END SUBROUTINE mao_project_gradient
570
571! **************************************************************************************************
572!> \brief ...
573!> \param fmat1 ...
574!> \param fmat2 ...
575!> \return ...
576! **************************************************************************************************
577 FUNCTION mao_scalar_product(fmat1, fmat2) RESULT(spro)
578 TYPE(dbcsr_type) :: fmat1, fmat2
579 REAL(kind=dp) :: spro
580
581 INTEGER :: iatom, jatom, m, n
582 LOGICAL :: found
583 REAL(kind=dp), DIMENSION(:, :), POINTER :: ablock, bblock
584 TYPE(dbcsr_iterator_type) :: dbcsr_iter
585 TYPE(mp_comm_type) :: group
586
587 spro = 0.0_dp
588
589 CALL dbcsr_iterator_start(dbcsr_iter, fmat1)
590 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
591 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock)
592 cpassert(iatom == jatom)
593 m = SIZE(ablock, 2)
594 n = SIZE(ablock, 1)
595 CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found)
596 cpassert(found)
597 spro = spro + sum(ablock(1:n, 1:m)*bblock(1:n, 1:m))
598 END DO
599 CALL dbcsr_iterator_stop(dbcsr_iter)
600
601 CALL dbcsr_get_info(fmat1, group=group)
602 CALL group%sum(spro)
603
604 END FUNCTION mao_scalar_product
605
606! **************************************************************************************************
607!> \brief Calculate the density matrix at the Gamma point
608!> \param pmat ...
609!> \param ksmat ...
610!> \param smat ...
611!> \param kpoints Kpoint environment
612!> \param nmos Number of occupied orbitals
613!> \param occ Maximum occupation per orbital
614!> \par History
615!> 04.2016 created [JGH]
616! **************************************************************************************************
617 SUBROUTINE calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
618
619 TYPE(dbcsr_type) :: pmat, ksmat, smat
620 TYPE(kpoint_type), POINTER :: kpoints
621 INTEGER, INTENT(IN) :: nmos
622 REAL(kind=dp), INTENT(IN) :: occ
623
624 INTEGER :: norb
625 REAL(kind=dp) :: de
626 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
627 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
628 TYPE(cp_fm_type) :: fmksmat, fmsmat, fmvec, fmwork
629 TYPE(dbcsr_type) :: tempmat
630
631 ! FM matrices
632
633 CALL dbcsr_get_info(smat, nfullrows_total=norb)
634 CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, &
635 nrow_global=norb, ncol_global=norb)
636 CALL cp_fm_create(fmksmat, matrix_struct)
637 CALL cp_fm_create(fmsmat, matrix_struct)
638 CALL cp_fm_create(fmvec, matrix_struct)
639 CALL cp_fm_create(fmwork, matrix_struct)
640 ALLOCATE (eigenvalues(norb))
641
642 ! DBCSR matrix
643 CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
644
645 ! transfer to FM
646 CALL dbcsr_desymmetrize(smat, tempmat)
647 CALL copy_dbcsr_to_fm(tempmat, fmsmat)
648 CALL dbcsr_desymmetrize(ksmat, tempmat)
649 CALL copy_dbcsr_to_fm(tempmat, fmksmat)
650
651 ! diagonalize
652 CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
653 de = eigenvalues(nmos + 1) - eigenvalues(nmos)
654 IF (de < 0.001_dp) THEN
655 CALL cp_warn(__location__, "MAO: No band gap at "// &
656 "Gamma point. MAO analysis not reliable.")
657 END IF
658 ! density matrix
659 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ)
660
661 DEALLOCATE (eigenvalues)
662 CALL dbcsr_release(tempmat)
663 CALL cp_fm_release(fmksmat)
664 CALL cp_fm_release(fmsmat)
665 CALL cp_fm_release(fmvec)
666 CALL cp_fm_release(fmwork)
667 CALL cp_fm_struct_release(matrix_struct)
668
669 END SUBROUTINE calculate_p_gamma
670
671! **************************************************************************************************
672!> \brief Define the MAO reference basis set
673!> \param qs_env ...
674!> \param mao_basis ...
675!> \param mao_basis_set_list ...
676!> \param orb_basis_set_list ...
677!> \param iunit ...
678!> \param print_basis ...
679!> \par History
680!> 07.2016 created [JGH]
681! **************************************************************************************************
682 SUBROUTINE mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, &
683 iunit, print_basis)
684
685 TYPE(qs_environment_type), POINTER :: qs_env
686 INTEGER, INTENT(IN) :: mao_basis
687 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
688 INTEGER, INTENT(IN), OPTIONAL :: iunit
689 LOGICAL, INTENT(IN), OPTIONAL :: print_basis
690
691 INTEGER :: ikind, nbas, nkind, unit_nr
692 REAL(kind=dp) :: eps_pgf_orb
693 TYPE(dft_control_type), POINTER :: dft_control
694 TYPE(gto_basis_set_type), POINTER :: basis_set, pbasis
695 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
696 TYPE(qs_kind_type), POINTER :: qs_kind
697
698 ! Reference basis set
699 cpassert(.NOT. ASSOCIATED(mao_basis_set_list))
700 cpassert(.NOT. ASSOCIATED(orb_basis_set_list))
701
702 ! options
703 IF (PRESENT(iunit)) THEN
704 unit_nr = iunit
705 ELSE
706 unit_nr = -1
707 END IF
708
709 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
710 nkind = SIZE(qs_kind_set)
711 ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
712 DO ikind = 1, nkind
713 NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
714 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
715 END DO
716 !
717 DO ikind = 1, nkind
718 qs_kind => qs_kind_set(ikind)
719 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
720 IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set
721 END DO
722 !
723 SELECT CASE (mao_basis)
724 CASE (mao_basis_orb)
725 DO ikind = 1, nkind
726 qs_kind => qs_kind_set(ikind)
727 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
728 IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set
729 END DO
730 CASE (mao_basis_prim)
731 DO ikind = 1, nkind
732 qs_kind => qs_kind_set(ikind)
733 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
734 NULLIFY (pbasis)
735 IF (ASSOCIATED(basis_set)) THEN
736 CALL create_primitive_basis_set(basis_set, pbasis)
737 CALL get_qs_env(qs_env, dft_control=dft_control)
738 eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
739 CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb)
740 pbasis%kind_radius = basis_set%kind_radius
741 mao_basis_set_list(ikind)%gto_basis_set => pbasis
742 CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO")
743 END IF
744 END DO
745 CASE (mao_basis_ext)
746 DO ikind = 1, nkind
747 qs_kind => qs_kind_set(ikind)
748 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO")
749 IF (ASSOCIATED(basis_set)) THEN
750 basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
751 mao_basis_set_list(ikind)%gto_basis_set => basis_set
752 END IF
753 END DO
754 CASE DEFAULT
755 cpabort("Unknown option for MAO basis")
756 END SELECT
757 IF (unit_nr > 0) THEN
758 DO ikind = 1, nkind
759 IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
760 WRITE (unit=unit_nr, fmt="(T2,A,I4)") &
761 "WARNING: No MAO basis set associated with Kind ", ikind
762 ELSE
763 nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
764 WRITE (unit=unit_nr, fmt="(T2,A,I4,T56,A,I10)") &
765 "MAO basis set Kind ", ikind, " Number of BSF:", nbas
766 END IF
767 END DO
768 END IF
769
770 IF (PRESENT(print_basis)) THEN
771 IF (print_basis) THEN
772 DO ikind = 1, nkind
773 basis_set => mao_basis_set_list(ikind)%gto_basis_set
774 IF (ASSOCIATED(basis_set)) CALL write_gto_basis_set(basis_set, unit_nr, "MAO REFERENCE BASIS")
775 END DO
776 END IF
777 END IF
778
779 END SUBROUTINE mao_reference_basis
780
781! **************************************************************************************************
782!> \brief Analyze the MAO basis, projection on angular functions
783!> \param mao_coef ...
784!> \param matrix_smm ...
785!> \param mao_basis_set_list ...
786!> \param particle_set ...
787!> \param qs_kind_set ...
788!> \param unit_nr ...
789!> \param para_env ...
790!> \par History
791!> 07.2016 created [JGH]
792! **************************************************************************************************
793 SUBROUTINE mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
794 qs_kind_set, unit_nr, para_env)
795
796 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_smm
797 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list
798 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
799 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
800 INTEGER, INTENT(IN) :: unit_nr
801 TYPE(mp_para_env_type), POINTER :: para_env
802
803 CHARACTER(len=2) :: element_symbol
804 INTEGER :: ia, iab, iatom, ikind, iset, ishell, &
805 ispin, l, lmax, lshell, m, ma, na, &
806 natom, nspin
807 LOGICAL :: found
808 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cmask, vec1, vec2
809 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: weight
810 REAL(kind=dp), DIMENSION(:, :), POINTER :: block, cmao
811 TYPE(gto_basis_set_type), POINTER :: basis_set
812
813 ! Analyze the MAO basis
814 IF (unit_nr > 0) THEN
815 WRITE (unit_nr, "(/,A)") " Analyze angular momentum character of MAOs "
816 WRITE (unit_nr, "(T7,A,T15,A,T20,A,T40,A,T50,A,T60,A,T70,A,T80,A)") &
817 "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G"
818 END IF
819 lmax = 4 ! analyze up to g-functions
820 natom = SIZE(particle_set)
821 nspin = SIZE(mao_coef)
822 DO iatom = 1, natom
823 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
824 element_symbol=element_symbol, kind_number=ikind)
825 basis_set => mao_basis_set_list(ikind)%gto_basis_set
826 CALL get_qs_kind(qs_kind_set(ikind), mao=na)
827 CALL get_gto_basis_set(basis_set, nsgf=ma)
828 ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na))
829 weight = 0.0_dp
830 CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, &
831 block=block, found=found)
832 DO ispin = 1, nspin
833 CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, &
834 block=cmao, found=found)
835 IF (found) THEN
836 DO l = 0, lmax
837 cmask = 0.0_dp
838 iab = 0
839 DO iset = 1, basis_set%nset
840 DO ishell = 1, basis_set%nshell(iset)
841 lshell = basis_set%l(ishell, iset)
842 DO m = -lshell, lshell
843 iab = iab + 1
844 IF (l == lshell) cmask(iab) = 1.0_dp
845 END DO
846 END DO
847 END DO
848 DO ia = 1, na
849 vec1(1:ma) = cmask*cmao(1:ma, ia)
850 vec2(1:ma) = matmul(block, vec1)
851 weight(l, ia) = sum(vec1(1:ma)*vec2(1:ma))
852 END DO
853 END DO
854 END IF
855 CALL para_env%sum(weight)
856 IF (unit_nr > 0) THEN
857 DO ia = 1, na
858 IF (ispin == 1 .AND. ia == 1) THEN
859 WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") &
860 iatom, element_symbol, ispin, ia, weight(0:lmax, ia)
861 ELSE
862 WRITE (unit_nr, "(T17,i2,T20,i3,T31,5F10.4)") ispin, ia, weight(0:lmax, ia)
863 END IF
864 END DO
865 END IF
866 END DO
867 DEALLOCATE (cmask, weight, vec1, vec2)
868 END DO
869 END SUBROUTINE mao_basis_analysis
870
871! **************************************************************************************************
872!> \brief Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap
873!> \param matrix_q ...
874!> \param matrix_p ...
875!> \param matrix_s ...
876!> \param matrix_smm ...
877!> \param matrix_smo ...
878!> \param smm_list ...
879!> \param electra ...
880!> \param eps_filter ...
881!> \param nimages ...
882!> \param kpoints ...
883!> \param matrix_ks ...
884!> \param sab_orb ...
885!> \par History
886!> 08.2016 created [JGH]
887! **************************************************************************************************
888 SUBROUTINE mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, &
889 electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
890
891 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q
892 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
893 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_smm, matrix_smo
894 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
895 POINTER :: smm_list
896 REAL(kind=dp), DIMENSION(2), INTENT(OUT) :: electra
897 REAL(kind=dp), INTENT(IN) :: eps_filter
898 INTEGER, INTENT(IN), OPTIONAL :: nimages
899 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
900 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
901 POINTER :: matrix_ks
902 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
903 OPTIONAL, POINTER :: sab_orb
904
905 INTEGER :: im, ispin, nim, nocc, norb, nspin
906 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
907 REAL(kind=dp) :: elex, xkp(3)
908 TYPE(dbcsr_type) :: ksmat, pmat, smat, tmat
909
910 nim = 1
911 IF (PRESENT(nimages)) nim = nimages
912 IF (nim > 1) THEN
913 cpassert(PRESENT(kpoints))
914 cpassert(PRESENT(matrix_ks))
915 cpassert(PRESENT(sab_orb))
916 END IF
917
918 ! Reference
919 nspin = SIZE(matrix_p, 1)
920 DO ispin = 1, nspin
921 electra(ispin) = 0.0_dp
922 DO im = 1, nim
923 CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex)
924 electra(ispin) = electra(ispin) + elex
925 END DO
926 END DO
927
928 ! Q matrix
929 NULLIFY (matrix_q)
930 CALL dbcsr_allocate_matrix_set(matrix_q, nspin)
931 DO ispin = 1, nspin
932 ALLOCATE (matrix_q(ispin)%matrix)
933 CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix)
934 CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list)
935 END DO
936 ! temp matrix
937 CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
938 ! Q=APA(T)
939 DO ispin = 1, nspin
940 IF (nim == 1) THEN
941 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, &
942 0.0_dp, tmat, filter_eps=eps_filter)
943 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
944 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
945 ELSE
946 ! k-points
947 CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix)
948 CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix)
949 CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix)
950 CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb)
951 CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb)
952 CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb)
953 NULLIFY (cell_to_index)
954 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
955 ! calculate density matrix at gamma point
956 xkp = 0.0_dp
957 ! transform KS and S matrices to the gamma point
958 CALL dbcsr_set(ksmat, 0.0_dp)
959 CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, &
960 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
961 CALL dbcsr_set(smat, 0.0_dp)
962 CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, &
963 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
964 norb = nint(electra(ispin))
965 nocc = mod(2, nspin) + 1
966 CALL calculate_p_gamma(pmat, ksmat, smat, kpoints, norb, real(nocc, kind=dp))
967 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, pmat, &
968 0.0_dp, tmat, filter_eps=eps_filter)
969 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
970 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
971 CALL dbcsr_release(pmat)
972 CALL dbcsr_release(smat)
973 CALL dbcsr_release(ksmat)
974 END IF
975 END DO
976 ! free temp matrix
977 CALL dbcsr_release(tmat)
978
979 END SUBROUTINE mao_build_q
980
981END MODULE mao_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
mathlib::invert_matrix
Definition mathlib.F:73
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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.
Definition atomic_kind_types.F:138
basis_set_container_types
Definition basis_set_container_types.F:13
basis_set_container_types::add_basis_set_to_container
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
Definition basis_set_container_types.F:161
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
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, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
basis_set_types::write_gto_basis_set
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.
Definition basis_set_types.F:1923
basis_set_types::create_primitive_basis_set
subroutine, public create_primitive_basis_set(basis_set, pbasis, lmax)
...
Definition basis_set_types.F:294
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_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1049
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:524
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:692
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:1095
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:807
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:970
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1011
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1203
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1141
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_reserve_diag_blocks
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
Definition cp_dbcsr_contrib.F:255
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
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,...
Definition cp_dbcsr_operations.F:691
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_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE. Use cuSOLVERMp directly when requested and large enough; otherwi...
Definition cp_fm_diag.F:1486
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:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::mao_basis_orb
integer, parameter, public mao_basis_orb
Definition input_constants.F:679
input_constants::mao_basis_ext
integer, parameter, public mao_basis_ext
Definition input_constants.F:679
input_constants::mao_basis_prim
integer, parameter, public mao_basis_prim
Definition input_constants.F:679
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::invert_hotelling
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 ...
Definition iterate_matrix.F:473
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::rskp_transform
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.
Definition kpoint_methods.F:1114
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:525
mao_basis
Calculate MAO's and analyze wavefunctions.
Definition mao_basis.F:15
mao_methods
Calculate MAO's and analyze wavefunctions.
Definition mao_methods.F:15
mao_methods::mao_project_gradient
subroutine, public mao_project_gradient(mao_coef, mao_grad, smat)
...
Definition mao_methods.F:542
mao_methods::mao_function_gradient
subroutine, public mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
...
Definition mao_methods.F:270
mao_methods::mao_reference_basis
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.
Definition mao_methods.F:684
mao_methods::mao_orthogonalization
subroutine, public mao_orthogonalization(mao_coef, smat)
...
Definition mao_methods.F:497
mao_methods::mao_update_coef
subroutine, public mao_update_coef(mao_up, mao_coef, mao_grad, smat, alpha)
...
Definition mao_methods.F:434
mao_methods::calculate_p_gamma
subroutine, public calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
Calculate the density matrix at the Gamma point.
Definition mao_methods.F:618
mao_methods::mao_initialization
subroutine, public mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
...
Definition mao_methods.F:92
mao_methods::mao_basis_analysis
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.
Definition mao_methods.F:795
mao_methods::mao_scalar_product
real(kind=dp) function, public mao_scalar_product(fmat1, fmat2)
...
Definition mao_methods.F:578
mao_methods::mao_function
subroutine, public mao_function(mao_coef, fval, qmat, smat, binv, reuse)
...
Definition mao_methods.F:227
mao_methods::mao_preconditioner
subroutine, public mao_preconditioner(mao_coef, mao_grad, qmat, smat)
Preconditioner for MAO optimization (seems not to work, but left for further tests)
Definition mao_methods.F:323
mao_methods::mao_build_q
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.
Definition mao_methods.F:890
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_interactions
Calculate the interaction radii for the operator matrix calculation.
Definition qs_interactions.F:17
qs_interactions::init_interaction_radii_orb_basis
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Definition qs_interactions.F:439
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:751
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
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
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:171
message_passing::mp_comm_type
Definition message_passing.F:155
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67