Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
80 : CONTAINS
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 16 : 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 16 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, row_blk, &
98 16 : row_blk_sizes
99 : LOGICAL :: found
100 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
101 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
102 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, pblock, sblock
103 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
104 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
105 16 : TYPE(mblocks), ALLOCATABLE, DIMENSION(:) :: mbl
106 : TYPE(mp_comm_type) :: group
107 :
108 16 : CALL dbcsr_get_info(mao_coef, nblkrows_total=nblk)
109 126 : ALLOCATE (mbl(nblk))
110 94 : DO i = 1, nblk
111 94 : NULLIFY (mbl(i)%mat, mbl(i)%eig)
112 : END DO
113 :
114 16 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
115 55 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
116 39 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
117 39 : CPASSERT(iatom == jatom)
118 39 : m = SIZE(cblock, 2)
119 39 : NULLIFY (pblock, sblock)
120 39 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pblock, found=found)
121 39 : CPASSERT(found)
122 39 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
123 39 : CPASSERT(found)
124 39 : n = SIZE(sblock, 1)
125 39 : lwork = MAX(n*n, 100)
126 390 : ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
127 20925 : amat(1:n, 1:n) = pblock(1:n, 1:n)
128 20925 : bmat(1:n, 1:n) = sblock(1:n, 1:n)
129 39 : info = 0
130 39 : CALL dsygv(1, "V", "U", n, amat, n, bmat, n, w, work, lwork, info)
131 39 : CPASSERT(info == 0)
132 234 : ALLOCATE (mbl(iatom)%mat(n, n), mbl(iatom)%eig(n))
133 39 : mbl(iatom)%n = n
134 39 : mbl(iatom)%ma = m
135 797 : DO i = 1, n
136 758 : mbl(iatom)%eig(i) = w(n - i + 1)
137 20925 : mbl(iatom)%mat(1:n, i) = amat(1:n, n - i + 1)
138 : END DO
139 2105 : cblock(1:n, 1:m) = amat(1:n, n:n - m + 1:-1)
140 172 : DEALLOCATE (amat, bmat, w, work)
141 : END DO
142 16 : CALL dbcsr_iterator_stop(dbcsr_iter)
143 :
144 16 : IF (eps1 < 10.0_dp) THEN
145 0 : CALL dbcsr_get_info(mao_coef, row_blk_size=row_blk_sizes, group=group)
146 0 : ALLOCATE (row_blk(nblk), mao_blk(nblk))
147 0 : mao_blk = 0
148 0 : row_blk = row_blk_sizes
149 0 : DO iatom = 1, nblk
150 0 : IF (ASSOCIATED(mbl(iatom)%mat)) THEN
151 0 : n = mbl(iatom)%n
152 0 : m = 0
153 0 : DO i = 1, n
154 0 : IF (mbl(iatom)%eig(i) < eps1) EXIT
155 0 : m = i
156 : END DO
157 0 : m = MAX(m, mbl(iatom)%ma)
158 0 : mbl(iatom)%ma = m
159 0 : mao_blk(iatom) = m
160 : END IF
161 : END DO
162 0 : CALL group%sum(mao_blk)
163 0 : CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist)
164 0 : CALL dbcsr_release(mao_coef)
165 : CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, &
166 0 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, col_blk_size=mao_blk)
167 0 : CALL dbcsr_reserve_diag_blocks(matrix=mao_coef)
168 0 : DEALLOCATE (mao_blk, row_blk)
169 : !
170 0 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
171 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
172 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
173 0 : CPASSERT(iatom == jatom)
174 0 : n = SIZE(cblock, 1)
175 0 : m = SIZE(cblock, 2)
176 0 : CPASSERT(n == mbl(iatom)%n .AND. m == mbl(iatom)%ma)
177 0 : cblock(1:n, 1:m) = mbl(iatom)%mat(1:n, 1:m)
178 : END DO
179 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
180 : !
181 : END IF
182 :
183 16 : IF (iolevel > 2) THEN
184 : CALL dbcsr_get_info(mao_coef, col_blk_size=col_blk_sizes, &
185 12 : row_blk_size=row_blk_sizes, group=group)
186 66 : DO iatom = 1, nblk
187 54 : n = row_blk_sizes(iatom)
188 54 : m = col_blk_sizes(iatom)
189 162 : ALLOCATE (w(n))
190 978 : w(1:n) = 0._dp
191 54 : IF (ASSOCIATED(mbl(iatom)%mat)) THEN
192 489 : w(1:n) = mbl(iatom)%eig(1:n)
193 : END IF
194 54 : CALL group%sum(w)
195 54 : IF (iw > 0) THEN
196 27 : WRITE (iw, '(A,i2,20F8.4)', ADVANCE="NO") " Spectrum/Gap ", iatom, w(1:m)
197 27 : WRITE (iw, '(A,F8.4)') " || ", w(m + 1)
198 : END IF
199 66 : DEALLOCATE (w)
200 : END DO
201 : END IF
202 :
203 16 : CALL mao_orthogonalization(mao_coef, smat)
204 :
205 94 : DO i = 1, nblk
206 78 : IF (ASSOCIATED(mbl(i)%mat)) THEN
207 39 : DEALLOCATE (mbl(i)%mat)
208 : END IF
209 94 : IF (ASSOCIATED(mbl(i)%eig)) THEN
210 39 : DEALLOCATE (mbl(i)%eig)
211 : END IF
212 : END DO
213 16 : DEALLOCATE (mbl)
214 :
215 48 : 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 804 : 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 402 : threshold = 1.e-8_dp
236 402 : convergence = 1.e-6_dp
237 : ! temp matrices
238 402 : CALL dbcsr_create(scmat, template=mao_coef)
239 402 : CALL dbcsr_create(bmat, template=binv)
240 402 : CALL dbcsr_create(tmat, template=qmat)
241 : ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
242 402 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
243 402 : 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 402 : norm_convergence=convergence, silent=.TRUE.)
247 : ! calculate Binv*C and T=C(T)*Binv*C
248 402 : CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
249 402 : CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
250 : ! function = Tr(Q*T)
251 402 : CALL dbcsr_dot(qmat, tmat, fval)
252 : ! free temp matrices
253 402 : CALL dbcsr_release(scmat)
254 402 : CALL dbcsr_release(bmat)
255 402 : CALL dbcsr_release(tmat)
256 :
257 402 : 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 416 : 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 208 : threshold = 1.e-8_dp
279 208 : convergence = 1.e-6_dp
280 : ! temp matrices
281 208 : CALL dbcsr_create(scmat, template=mao_coef)
282 208 : CALL dbcsr_create(bmat, template=binv)
283 208 : CALL dbcsr_create(tmat, template=qmat)
284 208 : CALL dbcsr_create(t2mat, template=scmat)
285 : ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
286 208 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
287 208 : 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 208 : norm_convergence=convergence, silent=.TRUE.)
291 : ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T)
292 208 : CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
293 208 : CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
294 : ! function = Tr(Q*T)
295 208 : 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 208 : retain_sparsity=.TRUE.)
299 : ! Gradient part 2: g = -2*S*T*X; X = Q*R
300 208 : CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat)
301 208 : 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 208 : retain_sparsity=.TRUE.)
304 : ! free temp matrices
305 208 : CALL dbcsr_release(scmat)
306 208 : CALL dbcsr_release(bmat)
307 208 : CALL dbcsr_release(tmat)
308 208 : CALL dbcsr_release(t2mat)
309 :
310 : ! apply Lagrange multiplyer
311 208 : CALL mao_grad_lagrange(mao_coef, mao_grad, smat)
312 :
313 208 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes
327 : LOGICAL :: found
328 : REAL(KIND=dp) :: eval_error, q_shift
329 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qinv, qsmat
330 0 : 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 0 : CALL dbcsr_create(scmat, template=mao_coef)
337 0 : 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 0 : 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 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
342 0 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
343 :
344 0 : q_shift = 0.001_dp
345 :
346 0 : CALL dbcsr_iterator_start(dbcsr_iter, mao_grad)
347 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
348 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, gblock)
349 0 : CPASSERT(iatom == jatom)
350 0 : m = SIZE(gblock, 1)
351 0 : n = SIZE(gblock, 2)
352 0 : CALL dbcsr_get_block_p(matrix=bmat, row=iatom, col=jatom, block=bblock, found=found)
353 0 : CPASSERT(found)
354 0 : CALL dbcsr_get_block_p(matrix=qmat, row=iatom, col=jatom, block=qblock, found=found)
355 0 : CPASSERT(found)
356 0 : ALLOCATE (qinv(m, m), qsmat(m, m))
357 0 : qsmat(1:m, 1:m) = qblock(1:m, 1:m)
358 0 : DO i = 1, m
359 0 : qsmat(i, i) = qsmat(i, i) + q_shift
360 : END DO
361 0 : CALL invert_matrix(qsmat, qinv, eval_error)
362 : ! g = q * g * b
363 0 : gblock(1:m, 1:n) = MATMUL(qinv(1:m, 1:m), MATMUL(gblock(1:m, 1:n), bblock(1:n, 1:n)))
364 0 : DEALLOCATE (qinv, qsmat)
365 : END DO
366 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
367 :
368 0 : CALL mao_project_gradient(mao_coef, mao_grad, smat)
369 :
370 : ! free temp matrices
371 0 : CALL dbcsr_release(scmat)
372 0 : CALL dbcsr_release(bmat)
373 :
374 0 : END SUBROUTINE mao_preconditioner
375 :
376 : ! **************************************************************************************************
377 : !> \brief ...
378 : !> \param mao_coef ...
379 : !> \param mao_grad ...
380 : !> \param smat ...
381 : ! **************************************************************************************************
382 208 : 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 208 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
388 208 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, lmat, rmat
389 208 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
390 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
391 :
392 208 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
393 766 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
394 558 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
395 558 : CPASSERT(iatom == jatom)
396 558 : m = SIZE(cblock, 2)
397 558 : n = SIZE(cblock, 1)
398 558 : NULLIFY (sblock)
399 558 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
400 558 : CPASSERT(found)
401 558 : lwork = MAX(n*n, 100)
402 8928 : ALLOCATE (amat(n, m), bmat(m, m), rmat(m, m), lmat(n, n), w(m), work(lwork))
403 1160634 : amat(1:n, 1:m) = MATMUL(sblock(1:n, 1:n), cblock(1:n, 1:m))
404 122610 : bmat(1:m, 1:m) = MATMUL(TRANSPOSE(amat(1:n, 1:m)), amat(1:n, 1:m))
405 558 : info = 0
406 558 : CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
407 558 : CPASSERT(info == 0)
408 1674 : CPASSERT(ALL(w > 0.0_dp))
409 1674 : w = 1.0_dp/w
410 1674 : DO i = 1, m
411 5022 : rmat(1:m, i) = bmat(1:m, i)*w(i)
412 : END DO
413 27342 : bmat(1:m, 1:m) = MATMUL(rmat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
414 2827626 : 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 558 : CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=iatom, block=gblock, found=found)
417 558 : CPASSERT(found)
418 1161192 : gblock(1:n, 1:m) = gblock(1:n, 1:m) - MATMUL(lmat(1:n, 1:n), gblock(1:n, 1:m))
419 2232 : DEALLOCATE (amat, bmat, rmat, lmat, w, work)
420 : END DO
421 208 : CALL dbcsr_iterator_stop(dbcsr_iter)
422 :
423 416 : 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 594 : 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 594 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
440 594 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, r1mat, r2mat
441 594 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock, ublock
442 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
443 :
444 594 : CALL dbcsr_copy(mao_up, mao_coef)
445 594 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
446 2196 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
447 1602 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
448 1602 : CPASSERT(iatom == jatom)
449 1602 : CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=iatom, block=gblock, found=found)
450 1602 : CPASSERT(found)
451 1602 : m = SIZE(cblock, 2)
452 1602 : n = SIZE(cblock, 1)
453 1602 : NULLIFY (sblock)
454 1602 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
455 1602 : CPASSERT(found)
456 1602 : CALL dbcsr_get_block_p(matrix=mao_up, row=iatom, col=iatom, block=ublock, found=found)
457 1602 : CPASSERT(found)
458 1602 : lwork = MAX(n*n, 100)
459 24030 : ALLOCATE (amat(n, m), bmat(m, m), r1mat(m, m), r2mat(m, m), w(m), work(lwork))
460 : !
461 98862 : amat(1:n, 1:m) = alpha*gblock(1:n, 1:m)
462 7085412 : bmat(1:m, 1:m) = MATMUL(TRANSPOSE(amat(1:n, 1:m)), MATMUL(sblock(1:n, 1:n), amat(1:n, 1:m)))
463 1602 : info = 0
464 1602 : CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
465 1602 : CPASSERT(info == 0)
466 4806 : DO i = 1, m
467 3204 : IF (w(i) < 1.E-12_dp) THEN
468 36 : w(i) = 0.0_dp
469 : ELSE
470 3168 : w(i) = SQRT(w(i))
471 : END IF
472 12816 : r1mat(1:m, i) = bmat(1:m, i)*COS(w(i))
473 4806 : IF (w(i) < 1.E-6_dp) THEN
474 180 : r2mat(1:m, i) = bmat(1:m, i)*(1._dp - w(i)**2/6._dp)
475 : ELSE
476 12636 : r2mat(1:m, i) = bmat(1:m, i)*SIN(w(i))/w(i)
477 : END IF
478 : END DO
479 75294 : r1mat(1:m, 1:m) = MATMUL(r1mat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
480 75294 : r2mat(1:m, 1:m) = MATMUL(r2mat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
481 : !
482 6408 : ublock(1:n, 1:m) = MATMUL(cblock(1:n, 1:m), r1mat(1:m, 1:m)) + &
483 777222 : MATMUL(amat(1:n, 1:m), r2mat(1:m, 1:m))
484 : !
485 8010 : DEALLOCATE (amat, bmat, r1mat, r2mat, w, work)
486 : END DO
487 594 : CALL dbcsr_iterator_stop(dbcsr_iter)
488 :
489 1188 : END SUBROUTINE mao_update_coef
490 :
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param mao_coef ...
494 : !> \param smat ...
495 : ! **************************************************************************************************
496 16 : 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 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
502 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, rmat
503 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, sblock
504 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
505 :
506 16 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
507 55 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
508 39 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
509 39 : CPASSERT(iatom == jatom)
510 39 : m = SIZE(cblock, 2)
511 39 : n = SIZE(cblock, 1)
512 39 : NULLIFY (sblock)
513 39 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
514 39 : CPASSERT(found)
515 39 : lwork = MAX(n*n, 100)
516 507 : ALLOCATE (amat(n, m), bmat(m, m), rmat(m, m), w(m), work(lwork))
517 66947 : amat(1:n, 1:m) = MATMUL(sblock(1:n, 1:n), cblock(1:n, 1:m))
518 7571 : bmat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), amat(1:n, 1:m))
519 39 : info = 0
520 39 : CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
521 39 : CPASSERT(info == 0)
522 117 : CPASSERT(ALL(w > 0.0_dp))
523 117 : w = 1.0_dp/SQRT(w)
524 117 : DO i = 1, m
525 351 : rmat(1:m, i) = bmat(1:m, i)*w(i)
526 : END DO
527 1833 : bmat(1:m, 1:m) = MATMUL(rmat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
528 11391 : cblock(1:n, 1:m) = MATMUL(cblock(1:n, 1:m), bmat(1:m, 1:m))
529 117 : DEALLOCATE (amat, bmat, rmat, w, work)
530 : END DO
531 16 : CALL dbcsr_iterator_stop(dbcsr_iter)
532 :
533 32 : END SUBROUTINE mao_orthogonalization
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param mao_coef ...
538 : !> \param mao_grad ...
539 : !> \param smat ...
540 : ! **************************************************************************************************
541 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
547 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
548 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
549 :
550 0 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
551 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
552 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
553 0 : CPASSERT(iatom == jatom)
554 0 : m = SIZE(cblock, 2)
555 0 : n = SIZE(cblock, 1)
556 0 : NULLIFY (sblock)
557 0 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
558 0 : CPASSERT(found)
559 0 : NULLIFY (gblock)
560 0 : CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found)
561 0 : CPASSERT(found)
562 0 : ALLOCATE (amat(m, m))
563 0 : amat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), MATMUL(sblock(1:n, 1:n), gblock(1:n, 1:m)))
564 0 : gblock(1:n, 1:m) = gblock(1:n, 1:m) - MATMUL(cblock(1:n, 1:m), amat(1:m, 1:m))
565 0 : DEALLOCATE (amat)
566 : END DO
567 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
568 :
569 0 : END SUBROUTINE mao_project_gradient
570 :
571 : ! **************************************************************************************************
572 : !> \brief ...
573 : !> \param fmat1 ...
574 : !> \param fmat2 ...
575 : !> \return ...
576 : ! **************************************************************************************************
577 416 : 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 208 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ablock, bblock
584 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
585 : TYPE(mp_comm_type) :: group
586 :
587 208 : spro = 0.0_dp
588 :
589 208 : CALL dbcsr_iterator_start(dbcsr_iter, fmat1)
590 766 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
591 558 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock)
592 558 : CPASSERT(iatom == jatom)
593 558 : m = SIZE(ablock, 2)
594 558 : n = SIZE(ablock, 1)
595 558 : CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found)
596 558 : CPASSERT(found)
597 34474 : spro = spro + SUM(ablock(1:n, 1:m)*bblock(1:n, 1:m))
598 : END DO
599 208 : CALL dbcsr_iterator_stop(dbcsr_iter)
600 :
601 208 : CALL dbcsr_get_info(fmat1, group=group)
602 208 : CALL group%sum(spro)
603 :
604 208 : 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 0 : 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 0 : CALL dbcsr_get_info(smat, nfullrows_total=norb)
634 : CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, &
635 0 : nrow_global=norb, ncol_global=norb)
636 0 : CALL cp_fm_create(fmksmat, matrix_struct)
637 0 : CALL cp_fm_create(fmsmat, matrix_struct)
638 0 : CALL cp_fm_create(fmvec, matrix_struct)
639 0 : CALL cp_fm_create(fmwork, matrix_struct)
640 0 : ALLOCATE (eigenvalues(norb))
641 :
642 : ! DBCSR matrix
643 0 : CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
644 :
645 : ! transfer to FM
646 0 : CALL dbcsr_desymmetrize(smat, tempmat)
647 0 : CALL copy_dbcsr_to_fm(tempmat, fmsmat)
648 0 : CALL dbcsr_desymmetrize(ksmat, tempmat)
649 0 : CALL copy_dbcsr_to_fm(tempmat, fmksmat)
650 :
651 : ! diagonalize
652 0 : CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
653 0 : de = eigenvalues(nmos + 1) - eigenvalues(nmos)
654 0 : IF (de < 0.001_dp) THEN
655 : CALL cp_warn(__LOCATION__, "MAO: No band gap at "// &
656 0 : "Gamma point. MAO analysis not reliable.")
657 : END IF
658 : ! density matrix
659 0 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ)
660 :
661 0 : DEALLOCATE (eigenvalues)
662 0 : CALL dbcsr_release(tempmat)
663 0 : CALL cp_fm_release(fmksmat)
664 0 : CALL cp_fm_release(fmsmat)
665 0 : CALL cp_fm_release(fmvec)
666 0 : CALL cp_fm_release(fmwork)
667 0 : CALL cp_fm_struct_release(matrix_struct)
668 :
669 0 : 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 10 : 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 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
696 : TYPE(qs_kind_type), POINTER :: qs_kind
697 :
698 : ! Reference basis set
699 0 : CPASSERT(.NOT. ASSOCIATED(mao_basis_set_list))
700 10 : CPASSERT(.NOT. ASSOCIATED(orb_basis_set_list))
701 :
702 : ! options
703 10 : IF (PRESENT(iunit)) THEN
704 10 : unit_nr = iunit
705 : ELSE
706 0 : unit_nr = -1
707 : END IF
708 :
709 10 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
710 10 : nkind = SIZE(qs_kind_set)
711 80 : ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
712 30 : DO ikind = 1, nkind
713 20 : NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
714 30 : NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
715 : END DO
716 : !
717 30 : DO ikind = 1, nkind
718 20 : qs_kind => qs_kind_set(ikind)
719 20 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
720 30 : IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set
721 : END DO
722 : !
723 12 : SELECT CASE (mao_basis)
724 : CASE (mao_basis_orb)
725 6 : DO ikind = 1, nkind
726 4 : qs_kind => qs_kind_set(ikind)
727 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
728 6 : IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set
729 : END DO
730 : CASE (mao_basis_prim)
731 6 : DO ikind = 1, nkind
732 4 : qs_kind => qs_kind_set(ikind)
733 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
734 4 : NULLIFY (pbasis)
735 6 : IF (ASSOCIATED(basis_set)) THEN
736 4 : CALL create_primitive_basis_set(basis_set, pbasis)
737 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
738 4 : eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
739 4 : CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb)
740 4 : pbasis%kind_radius = basis_set%kind_radius
741 4 : mao_basis_set_list(ikind)%gto_basis_set => pbasis
742 4 : CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO")
743 : END IF
744 : END DO
745 : CASE (mao_basis_ext)
746 18 : DO ikind = 1, nkind
747 12 : qs_kind => qs_kind_set(ikind)
748 12 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO")
749 18 : IF (ASSOCIATED(basis_set)) THEN
750 12 : basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
751 12 : mao_basis_set_list(ikind)%gto_basis_set => basis_set
752 : END IF
753 : END DO
754 : CASE DEFAULT
755 10 : CPABORT("Unknown option for MAO basis")
756 : END SELECT
757 10 : IF (unit_nr > 0) THEN
758 15 : DO ikind = 1, nkind
759 15 : IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
760 : WRITE (UNIT=unit_nr, FMT="(T2,A,I4)") &
761 0 : "WARNING: No MAO basis set associated with Kind ", ikind
762 : ELSE
763 10 : nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
764 : WRITE (UNIT=unit_nr, FMT="(T2,A,I4,T56,A,I10)") &
765 10 : "MAO basis set Kind ", ikind, " Number of BSF:", nbas
766 : END IF
767 : END DO
768 : END IF
769 :
770 10 : IF (PRESENT(print_basis)) THEN
771 10 : IF (print_basis) THEN
772 0 : DO ikind = 1, nkind
773 0 : basis_set => mao_basis_set_list(ikind)%gto_basis_set
774 0 : 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 10 : 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 10 : 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 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cmask, vec1, vec2
809 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weight
810 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, cmao
811 : TYPE(gto_basis_set_type), POINTER :: basis_set
812 :
813 : ! Analyze the MAO basis
814 10 : IF (unit_nr > 0) THEN
815 5 : 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 5 : "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G"
818 : END IF
819 10 : lmax = 4 ! analyze up to g-functions
820 10 : natom = SIZE(particle_set)
821 10 : nspin = SIZE(mao_coef)
822 58 : DO iatom = 1, natom
823 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
824 48 : element_symbol=element_symbol, kind_number=ikind)
825 48 : basis_set => mao_basis_set_list(ikind)%gto_basis_set
826 48 : CALL get_qs_kind(qs_kind_set(ikind), mao=na)
827 48 : CALL get_gto_basis_set(basis_set, nsgf=ma)
828 336 : ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na))
829 48 : weight = 0.0_dp
830 : CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, &
831 48 : block=block, found=found)
832 102 : DO ispin = 1, nspin
833 : CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, &
834 54 : block=cmao, found=found)
835 54 : IF (found) THEN
836 162 : DO l = 0, lmax
837 135 : cmask = 0.0_dp
838 135 : iab = 0
839 675 : DO iset = 1, basis_set%nset
840 1635 : DO ishell = 1, basis_set%nshell(iset)
841 960 : lshell = basis_set%l(ishell, iset)
842 3810 : DO m = -lshell, lshell
843 2310 : iab = iab + 1
844 3270 : IF (l == lshell) cmask(iab) = 1.0_dp
845 : END DO
846 : END DO
847 : END DO
848 432 : DO ia = 1, na
849 6450 : vec1(1:ma) = cmask*cmao(1:ma, ia)
850 383430 : vec2(1:ma) = MATMUL(block, vec1)
851 6585 : weight(l, ia) = SUM(vec1(1:ma)*vec2(1:ma))
852 : END DO
853 : END DO
854 : END IF
855 54 : CALL para_env%sum(weight)
856 156 : IF (unit_nr > 0) THEN
857 81 : DO ia = 1, na
858 81 : IF (ispin == 1 .AND. ia == 1) THEN
859 : WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") &
860 24 : iatom, element_symbol, ispin, ia, weight(0:lmax, ia)
861 : ELSE
862 30 : 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 154 : DEALLOCATE (cmask, weight, vec1, vec2)
868 : END DO
869 20 : 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 14 : 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 14 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
907 : REAL(KIND=dp) :: elex, xkp(3)
908 : TYPE(dbcsr_type) :: ksmat, pmat, smat, tmat
909 :
910 14 : nim = 1
911 14 : IF (PRESENT(nimages)) nim = nimages
912 0 : IF (nim > 1) THEN
913 0 : CPASSERT(PRESENT(kpoints))
914 0 : CPASSERT(PRESENT(matrix_ks))
915 0 : CPASSERT(PRESENT(sab_orb))
916 : END IF
917 :
918 : ! Reference
919 14 : nspin = SIZE(matrix_p, 1)
920 30 : DO ispin = 1, nspin
921 16 : electra(ispin) = 0.0_dp
922 46 : DO im = 1, nim
923 16 : CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex)
924 32 : electra(ispin) = electra(ispin) + elex
925 : END DO
926 : END DO
927 :
928 : ! Q matrix
929 14 : NULLIFY (matrix_q)
930 14 : CALL dbcsr_allocate_matrix_set(matrix_q, nspin)
931 30 : DO ispin = 1, nspin
932 16 : ALLOCATE (matrix_q(ispin)%matrix)
933 16 : CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix)
934 30 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list)
935 : END DO
936 : ! temp matrix
937 14 : CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
938 : ! Q=APA(T)
939 30 : DO ispin = 1, nspin
940 30 : IF (nim == 1) THEN
941 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, &
942 16 : 0.0_dp, tmat, filter_eps=eps_filter)
943 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
944 16 : 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
945 : ELSE
946 : ! k-points
947 0 : CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix)
948 0 : CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix)
949 0 : CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix)
950 0 : CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb)
951 0 : CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb)
952 0 : CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb)
953 0 : NULLIFY (cell_to_index)
954 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
955 : ! calculate density matrix at gamma point
956 0 : xkp = 0.0_dp
957 : ! transform KS and S matrices to the gamma point
958 0 : CALL dbcsr_set(ksmat, 0.0_dp)
959 : CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, &
960 0 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
961 0 : CALL dbcsr_set(smat, 0.0_dp)
962 : CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, &
963 0 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
964 0 : norb = NINT(electra(ispin))
965 0 : nocc = MOD(2, nspin) + 1
966 0 : 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.0_dp, tmat, filter_eps=eps_filter)
969 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
970 0 : 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
971 0 : CALL dbcsr_release(pmat)
972 0 : CALL dbcsr_release(smat)
973 0 : CALL dbcsr_release(ksmat)
974 : END IF
975 : END DO
976 : ! free temp matrix
977 14 : CALL dbcsr_release(tmat)
978 :
979 14 : END SUBROUTINE mao_build_q
980 :
981 4320 : END MODULE mao_methods
|
