Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, 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_reserve_diag_blocks, dbcsr_set, dbcsr_type, &
29 : dbcsr_type_no_symmetry
30 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : cp_dbcsr_plus_fm_fm_t,&
33 : dbcsr_allocate_matrix_set
34 : USE cp_fm_diag, ONLY: cp_fm_geeig
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_release,&
37 : cp_fm_struct_type
38 : USE cp_fm_types, ONLY: cp_fm_create,&
39 : cp_fm_release,&
40 : cp_fm_type
41 : USE input_constants, ONLY: mao_basis_ext,&
42 : mao_basis_orb,&
43 : mao_basis_prim
44 : USE iterate_matrix, ONLY: invert_Hotelling
45 : USE kinds, ONLY: dp
46 : USE kpoint_methods, ONLY: rskp_transform
47 : USE kpoint_types, ONLY: get_kpoint_info,&
48 : kpoint_type
49 : USE message_passing, ONLY: mp_comm_type,&
50 : mp_para_env_type
51 : USE particle_types, ONLY: particle_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
55 : USE qs_kind_types, ONLY: get_qs_kind,&
56 : qs_kind_type
57 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_methods'
64 :
65 : TYPE mblocks
66 : INTEGER :: n = -1, ma = -1
67 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat => NULL()
68 : REAL(KIND=dp), DIMENSION(:), POINTER :: eig => NULL()
69 : END TYPE mblocks
70 :
71 : PUBLIC :: mao_initialization, mao_function, mao_function_gradient, mao_orthogonalization, &
72 : mao_project_gradient, mao_scalar_product, mao_build_q, mao_basis_analysis, &
73 : mao_reference_basis, calculate_p_gamma
74 :
75 : ! **************************************************************************************************
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief ...
81 : !> \param mao_coef ...
82 : !> \param pmat ...
83 : !> \param smat ...
84 : !> \param eps1 ...
85 : !> \param iolevel ...
86 : !> \param iw ...
87 : ! **************************************************************************************************
88 16 : SUBROUTINE mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
89 : TYPE(dbcsr_type) :: mao_coef, pmat, smat
90 : REAL(KIND=dp), INTENT(IN) :: eps1
91 : INTEGER, INTENT(IN) :: iolevel, iw
92 :
93 : INTEGER :: group_handle, i, iatom, info, jatom, &
94 : lwork, m, n, nblk
95 16 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, row_blk, &
96 16 : row_blk_sizes
97 : LOGICAL :: found
98 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
99 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
100 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, pblock, sblock
101 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
102 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
103 16 : TYPE(mblocks), ALLOCATABLE, DIMENSION(:) :: mbl
104 : TYPE(mp_comm_type) :: group
105 :
106 16 : CALL dbcsr_get_info(mao_coef, nblkrows_total=nblk)
107 126 : ALLOCATE (mbl(nblk))
108 94 : DO i = 1, nblk
109 94 : NULLIFY (mbl(i)%mat, mbl(i)%eig)
110 : END DO
111 :
112 16 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
113 55 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
114 39 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
115 39 : CPASSERT(iatom == jatom)
116 39 : m = SIZE(cblock, 2)
117 39 : NULLIFY (pblock, sblock)
118 39 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pblock, found=found)
119 39 : CPASSERT(found)
120 39 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
121 39 : CPASSERT(found)
122 39 : n = SIZE(sblock, 1)
123 39 : lwork = MAX(n*n, 100)
124 390 : ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
125 20925 : amat(1:n, 1:n) = pblock(1:n, 1:n)
126 20925 : bmat(1:n, 1:n) = sblock(1:n, 1:n)
127 39 : info = 0
128 39 : CALL dsygv(1, "V", "U", n, amat, n, bmat, n, w, work, lwork, info)
129 39 : CPASSERT(info == 0)
130 234 : ALLOCATE (mbl(iatom)%mat(n, n), mbl(iatom)%eig(n))
131 39 : mbl(iatom)%n = n
132 39 : mbl(iatom)%ma = m
133 797 : DO i = 1, n
134 758 : mbl(iatom)%eig(i) = w(n - i + 1)
135 20925 : mbl(iatom)%mat(1:n, i) = amat(1:n, n - i + 1)
136 : END DO
137 2105 : cblock(1:n, 1:m) = amat(1:n, n:n - m + 1:-1)
138 172 : DEALLOCATE (amat, bmat, w, work)
139 : END DO
140 16 : CALL dbcsr_iterator_stop(dbcsr_iter)
141 :
142 16 : IF (eps1 < 10.0_dp) THEN
143 0 : CALL dbcsr_get_info(mao_coef, row_blk_size=row_blk_sizes, group=group_handle)
144 0 : CALL group%set_handle(group_handle)
145 0 : ALLOCATE (row_blk(nblk), mao_blk(nblk))
146 0 : mao_blk = 0
147 0 : row_blk = row_blk_sizes
148 0 : DO iatom = 1, nblk
149 0 : IF (ASSOCIATED(mbl(iatom)%mat)) THEN
150 0 : n = mbl(iatom)%n
151 0 : m = 0
152 0 : DO i = 1, n
153 0 : IF (mbl(iatom)%eig(i) < eps1) EXIT
154 0 : m = i
155 : END DO
156 0 : m = MAX(m, mbl(iatom)%ma)
157 0 : mbl(iatom)%ma = m
158 0 : mao_blk(iatom) = m
159 : END IF
160 : END DO
161 0 : CALL group%sum(mao_blk)
162 0 : CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist)
163 0 : CALL dbcsr_release(mao_coef)
164 : CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, &
165 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, &
166 0 : col_blk_size=mao_blk, nze=0)
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_handle)
186 12 : CALL group%set_handle(group_handle)
187 66 : DO iatom = 1, nblk
188 54 : n = row_blk_sizes(iatom)
189 54 : m = col_blk_sizes(iatom)
190 162 : ALLOCATE (w(n))
191 978 : w(1:n) = 0._dp
192 54 : IF (ASSOCIATED(mbl(iatom)%mat)) THEN
193 489 : w(1:n) = mbl(iatom)%eig(1:n)
194 : END IF
195 54 : CALL group%sum(w)
196 54 : IF (iw > 0) THEN
197 27 : WRITE (iw, '(A,i2,20F8.4)', ADVANCE="NO") " Spectrum/Gap ", iatom, w(1:m)
198 27 : WRITE (iw, '(A,F8.4)') " || ", w(m + 1)
199 : END IF
200 66 : DEALLOCATE (w)
201 : END DO
202 : END IF
203 :
204 16 : CALL mao_orthogonalization(mao_coef, smat)
205 :
206 94 : DO i = 1, nblk
207 78 : IF (ASSOCIATED(mbl(i)%mat)) THEN
208 39 : DEALLOCATE (mbl(i)%mat)
209 : END IF
210 94 : IF (ASSOCIATED(mbl(i)%eig)) THEN
211 39 : DEALLOCATE (mbl(i)%eig)
212 : END IF
213 : END DO
214 16 : DEALLOCATE (mbl)
215 :
216 48 : END SUBROUTINE mao_initialization
217 :
218 : ! **************************************************************************************************
219 : !> \brief ...
220 : !> \param mao_coef ...
221 : !> \param fval ...
222 : !> \param qmat ...
223 : !> \param smat ...
224 : !> \param binv ...
225 : !> \param reuse ...
226 : ! **************************************************************************************************
227 318 : SUBROUTINE mao_function(mao_coef, fval, qmat, smat, binv, reuse)
228 : TYPE(dbcsr_type) :: mao_coef
229 : REAL(KIND=dp), INTENT(OUT) :: fval
230 : TYPE(dbcsr_type) :: qmat, smat, binv
231 : LOGICAL, INTENT(IN) :: reuse
232 :
233 : REAL(KIND=dp) :: convergence, threshold
234 : TYPE(dbcsr_type) :: bmat, scmat, tmat
235 :
236 318 : threshold = 1.e-8_dp
237 318 : convergence = 1.e-6_dp
238 : ! temp matrices
239 318 : CALL dbcsr_create(scmat, template=mao_coef)
240 318 : CALL dbcsr_create(bmat, template=binv)
241 318 : CALL dbcsr_create(tmat, template=qmat)
242 : ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
243 318 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
244 318 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
245 : ! calculate inverse of B
246 : CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
247 318 : norm_convergence=convergence, silent=.TRUE.)
248 : ! calculate Binv*C and T=C(T)*Binv*C
249 318 : CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
250 318 : CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
251 : ! function = Tr(Q*T)
252 318 : CALL dbcsr_dot(qmat, tmat, fval)
253 : ! free temp matrices
254 318 : CALL dbcsr_release(scmat)
255 318 : CALL dbcsr_release(bmat)
256 318 : CALL dbcsr_release(tmat)
257 :
258 318 : END SUBROUTINE mao_function
259 :
260 : ! **************************************************************************************************
261 : !> \brief ...
262 : !> \param mao_coef ...
263 : !> \param fval ...
264 : !> \param mao_grad ...
265 : !> \param qmat ...
266 : !> \param smat ...
267 : !> \param binv ...
268 : !> \param reuse ...
269 : ! **************************************************************************************************
270 166 : SUBROUTINE mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
271 : TYPE(dbcsr_type) :: mao_coef
272 : REAL(KIND=dp), INTENT(OUT) :: fval
273 : TYPE(dbcsr_type) :: mao_grad, qmat, smat, binv
274 : LOGICAL, INTENT(IN) :: reuse
275 :
276 : REAL(KIND=dp) :: convergence, threshold
277 : TYPE(dbcsr_type) :: bmat, scmat, t2mat, tmat
278 :
279 166 : threshold = 1.e-8_dp
280 166 : convergence = 1.e-6_dp
281 : ! temp matrices
282 166 : CALL dbcsr_create(scmat, template=mao_coef)
283 166 : CALL dbcsr_create(bmat, template=binv)
284 166 : CALL dbcsr_create(tmat, template=qmat)
285 166 : CALL dbcsr_create(t2mat, template=scmat)
286 : ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
287 166 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
288 166 : CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
289 : ! calculate inverse of B
290 : CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
291 166 : norm_convergence=convergence, silent=.TRUE.)
292 : ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T)
293 166 : CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
294 166 : CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
295 : ! function = Tr(Q*T)
296 166 : CALL dbcsr_dot(qmat, tmat, fval)
297 : ! Gradient part 1: g = 2*Q*C*Binv = 2*Q*R
298 : CALL dbcsr_multiply("N", "N", 2.0_dp, qmat, scmat, 0.0_dp, mao_grad, &
299 166 : retain_sparsity=.TRUE.)
300 : ! Gradient part 2: g = -2*S*T*X; X = Q*R
301 166 : CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat)
302 166 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, t2mat, 0.0_dp, scmat)
303 : CALL dbcsr_multiply("N", "N", -2.0_dp, smat, scmat, 1.0_dp, mao_grad, &
304 166 : retain_sparsity=.TRUE.)
305 : ! free temp matrices
306 166 : CALL dbcsr_release(scmat)
307 166 : CALL dbcsr_release(bmat)
308 166 : CALL dbcsr_release(tmat)
309 166 : CALL dbcsr_release(t2mat)
310 :
311 166 : CALL mao_project_gradient(mao_coef, mao_grad, smat)
312 :
313 166 : END SUBROUTINE mao_function_gradient
314 :
315 : ! **************************************************************************************************
316 : !> \brief ...
317 : !> \param mao_coef ...
318 : !> \param smat ...
319 : ! **************************************************************************************************
320 484 : SUBROUTINE mao_orthogonalization(mao_coef, smat)
321 : TYPE(dbcsr_type) :: mao_coef, smat
322 :
323 : INTEGER :: i, iatom, info, jatom, lwork, m, n
324 : LOGICAL :: found
325 484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
326 484 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
327 484 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, sblock
328 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
329 :
330 484 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
331 1747 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
332 1263 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
333 1263 : CPASSERT(iatom == jatom)
334 1263 : m = SIZE(cblock, 2)
335 1263 : n = SIZE(cblock, 1)
336 1263 : NULLIFY (sblock)
337 1263 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
338 1263 : CPASSERT(found)
339 1263 : lwork = MAX(n*n, 100)
340 13893 : ALLOCATE (amat(n, m), bmat(m, m), w(m), work(lwork))
341 2667347 : amat(1:n, 1:m) = MATMUL(sblock(1:n, 1:n), cblock(1:n, 1:m))
342 278171 : bmat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), amat(1:n, 1:m))
343 1263 : info = 0
344 1263 : CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
345 1263 : CPASSERT(info == 0)
346 3789 : CPASSERT(ALL(w > 0.0_dp))
347 3789 : w = 1.0_dp/SQRT(w)
348 3789 : DO i = 1, m
349 11367 : amat(1:m, i) = bmat(1:m, i)*w(i)
350 : END DO
351 87147 : bmat(1:m, 1:m) = MATMUL(amat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
352 679283 : cblock(1:n, 1:m) = MATMUL(cblock(1:n, 1:m), bmat(1:m, 1:m))
353 3789 : DEALLOCATE (amat, bmat, w, work)
354 : END DO
355 484 : CALL dbcsr_iterator_stop(dbcsr_iter)
356 :
357 968 : END SUBROUTINE mao_orthogonalization
358 :
359 : ! **************************************************************************************************
360 : !> \brief ...
361 : !> \param mao_coef ...
362 : !> \param mao_grad ...
363 : !> \param smat ...
364 : ! **************************************************************************************************
365 166 : SUBROUTINE mao_project_gradient(mao_coef, mao_grad, smat)
366 : TYPE(dbcsr_type) :: mao_coef, mao_grad, smat
367 :
368 : INTEGER :: iatom, jatom, m, n
369 : LOGICAL :: found
370 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
371 166 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock
372 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
373 :
374 166 : CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
375 598 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
376 432 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
377 432 : CPASSERT(iatom == jatom)
378 432 : m = SIZE(cblock, 2)
379 432 : n = SIZE(cblock, 1)
380 432 : NULLIFY (sblock)
381 432 : CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
382 432 : CPASSERT(found)
383 432 : NULLIFY (gblock)
384 432 : CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found)
385 432 : CPASSERT(found)
386 1728 : ALLOCATE (amat(m, m))
387 1948800 : amat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), MATMUL(sblock(1:n, 1:n), gblock(1:n, 1:m)))
388 145440 : gblock(1:n, 1:m) = gblock(1:n, 1:m) - MATMUL(cblock(1:n, 1:m), amat(1:m, 1:m))
389 1728 : DEALLOCATE (amat)
390 : END DO
391 166 : CALL dbcsr_iterator_stop(dbcsr_iter)
392 :
393 332 : END SUBROUTINE mao_project_gradient
394 :
395 : ! **************************************************************************************************
396 : !> \brief ...
397 : !> \param fmat1 ...
398 : !> \param fmat2 ...
399 : !> \return ...
400 : ! **************************************************************************************************
401 332 : FUNCTION mao_scalar_product(fmat1, fmat2) RESULT(spro)
402 : TYPE(dbcsr_type) :: fmat1, fmat2
403 : REAL(KIND=dp) :: spro
404 :
405 : INTEGER :: group_handle, iatom, jatom, m, n
406 : LOGICAL :: found
407 166 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ablock, bblock
408 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
409 : TYPE(mp_comm_type) :: group
410 :
411 166 : spro = 0.0_dp
412 :
413 166 : CALL dbcsr_iterator_start(dbcsr_iter, fmat1)
414 598 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
415 432 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock)
416 432 : CPASSERT(iatom == jatom)
417 432 : m = SIZE(ablock, 2)
418 432 : n = SIZE(ablock, 1)
419 432 : CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found)
420 432 : CPASSERT(found)
421 26950 : spro = spro + SUM(ablock(1:n, 1:m)*bblock(1:n, 1:m))
422 : END DO
423 166 : CALL dbcsr_iterator_stop(dbcsr_iter)
424 :
425 166 : CALL dbcsr_get_info(fmat1, group=group_handle)
426 166 : CALL group%set_handle(group_handle)
427 166 : CALL group%sum(spro)
428 :
429 166 : END FUNCTION mao_scalar_product
430 :
431 : ! **************************************************************************************************
432 : !> \brief Calculate the density matrix at the Gamma point
433 : !> \param pmat ...
434 : !> \param ksmat ...
435 : !> \param smat ...
436 : !> \param kpoints Kpoint environment
437 : !> \param nmos Number of occupied orbitals
438 : !> \param occ Maximum occupation per orbital
439 : !> \par History
440 : !> 04.2016 created [JGH]
441 : ! **************************************************************************************************
442 0 : SUBROUTINE calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
443 :
444 : TYPE(dbcsr_type) :: pmat, ksmat, smat
445 : TYPE(kpoint_type), POINTER :: kpoints
446 : INTEGER, INTENT(IN) :: nmos
447 : REAL(KIND=dp), INTENT(IN) :: occ
448 :
449 : INTEGER :: norb
450 : REAL(KIND=dp) :: de
451 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
452 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
453 : TYPE(cp_fm_type) :: fmksmat, fmsmat, fmvec, fmwork
454 : TYPE(dbcsr_type) :: tempmat
455 :
456 : ! FM matrices
457 :
458 0 : CALL dbcsr_get_info(smat, nfullrows_total=norb)
459 : CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, &
460 0 : nrow_global=norb, ncol_global=norb)
461 0 : CALL cp_fm_create(fmksmat, matrix_struct)
462 0 : CALL cp_fm_create(fmsmat, matrix_struct)
463 0 : CALL cp_fm_create(fmvec, matrix_struct)
464 0 : CALL cp_fm_create(fmwork, matrix_struct)
465 0 : ALLOCATE (eigenvalues(norb))
466 :
467 : ! DBCSR matrix
468 0 : CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
469 :
470 : ! transfer to FM
471 0 : CALL dbcsr_desymmetrize(smat, tempmat)
472 0 : CALL copy_dbcsr_to_fm(tempmat, fmsmat)
473 0 : CALL dbcsr_desymmetrize(ksmat, tempmat)
474 0 : CALL copy_dbcsr_to_fm(tempmat, fmksmat)
475 :
476 : ! diagonalize
477 0 : CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
478 0 : de = eigenvalues(nmos + 1) - eigenvalues(nmos)
479 0 : IF (de < 0.001_dp) THEN
480 : CALL cp_warn(__LOCATION__, "MAO: No band gap at "// &
481 0 : "Gamma point. MAO analysis not reliable.")
482 : END IF
483 : ! density matrix
484 0 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ)
485 :
486 0 : DEALLOCATE (eigenvalues)
487 0 : CALL dbcsr_release(tempmat)
488 0 : CALL cp_fm_release(fmksmat)
489 0 : CALL cp_fm_release(fmsmat)
490 0 : CALL cp_fm_release(fmvec)
491 0 : CALL cp_fm_release(fmwork)
492 0 : CALL cp_fm_struct_release(matrix_struct)
493 :
494 0 : END SUBROUTINE calculate_p_gamma
495 :
496 : ! **************************************************************************************************
497 : !> \brief Define the MAO reference basis set
498 : !> \param qs_env ...
499 : !> \param mao_basis ...
500 : !> \param mao_basis_set_list ...
501 : !> \param orb_basis_set_list ...
502 : !> \param iunit ...
503 : !> \param print_basis ...
504 : !> \par History
505 : !> 07.2016 created [JGH]
506 : ! **************************************************************************************************
507 10 : SUBROUTINE mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, &
508 : iunit, print_basis)
509 :
510 : TYPE(qs_environment_type), POINTER :: qs_env
511 : INTEGER, INTENT(IN) :: mao_basis
512 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
513 : INTEGER, INTENT(IN), OPTIONAL :: iunit
514 : LOGICAL, INTENT(IN), OPTIONAL :: print_basis
515 :
516 : INTEGER :: ikind, nbas, nkind, unit_nr
517 : REAL(KIND=dp) :: eps_pgf_orb
518 : TYPE(dft_control_type), POINTER :: dft_control
519 : TYPE(gto_basis_set_type), POINTER :: basis_set, pbasis
520 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
521 : TYPE(qs_kind_type), POINTER :: qs_kind
522 :
523 : ! Reference basis set
524 0 : CPASSERT(.NOT. ASSOCIATED(mao_basis_set_list))
525 10 : CPASSERT(.NOT. ASSOCIATED(orb_basis_set_list))
526 :
527 : ! options
528 10 : IF (PRESENT(iunit)) THEN
529 10 : unit_nr = iunit
530 : ELSE
531 0 : unit_nr = -1
532 : END IF
533 :
534 10 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
535 10 : nkind = SIZE(qs_kind_set)
536 90 : ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
537 30 : DO ikind = 1, nkind
538 20 : NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
539 30 : NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
540 : END DO
541 : !
542 30 : DO ikind = 1, nkind
543 20 : qs_kind => qs_kind_set(ikind)
544 20 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
545 30 : IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set
546 : END DO
547 : !
548 12 : SELECT CASE (mao_basis)
549 : CASE (mao_basis_orb)
550 6 : DO ikind = 1, nkind
551 4 : qs_kind => qs_kind_set(ikind)
552 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
553 6 : IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set
554 : END DO
555 : CASE (mao_basis_prim)
556 6 : DO ikind = 1, nkind
557 4 : qs_kind => qs_kind_set(ikind)
558 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
559 4 : NULLIFY (pbasis)
560 6 : IF (ASSOCIATED(basis_set)) THEN
561 4 : CALL create_primitive_basis_set(basis_set, pbasis)
562 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
563 4 : eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
564 4 : CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb)
565 4 : pbasis%kind_radius = basis_set%kind_radius
566 4 : mao_basis_set_list(ikind)%gto_basis_set => pbasis
567 4 : CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO")
568 : END IF
569 : END DO
570 : CASE (mao_basis_ext)
571 18 : DO ikind = 1, nkind
572 12 : qs_kind => qs_kind_set(ikind)
573 12 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO")
574 18 : IF (ASSOCIATED(basis_set)) THEN
575 12 : basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
576 12 : mao_basis_set_list(ikind)%gto_basis_set => basis_set
577 : END IF
578 : END DO
579 : CASE DEFAULT
580 10 : CPABORT("Unknown option for MAO basis")
581 : END SELECT
582 10 : IF (unit_nr > 0) THEN
583 15 : DO ikind = 1, nkind
584 15 : IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
585 : WRITE (UNIT=unit_nr, FMT="(T2,A,I4)") &
586 0 : "WARNING: No MAO basis set associated with Kind ", ikind
587 : ELSE
588 10 : nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
589 : WRITE (UNIT=unit_nr, FMT="(T2,A,I4,T56,A,I10)") &
590 10 : "MAO basis set Kind ", ikind, " Number of BSF:", nbas
591 : END IF
592 : END DO
593 : END IF
594 :
595 10 : IF (PRESENT(print_basis)) THEN
596 10 : IF (print_basis) THEN
597 0 : DO ikind = 1, nkind
598 0 : basis_set => mao_basis_set_list(ikind)%gto_basis_set
599 0 : IF (ASSOCIATED(basis_set)) CALL write_gto_basis_set(basis_set, unit_nr, "MAO REFERENCE BASIS")
600 : END DO
601 : END IF
602 : END IF
603 :
604 10 : END SUBROUTINE mao_reference_basis
605 :
606 : ! **************************************************************************************************
607 : !> \brief Analyze the MAO basis, projection on angular functions
608 : !> \param mao_coef ...
609 : !> \param matrix_smm ...
610 : !> \param mao_basis_set_list ...
611 : !> \param particle_set ...
612 : !> \param qs_kind_set ...
613 : !> \param unit_nr ...
614 : !> \param para_env ...
615 : !> \par History
616 : !> 07.2016 created [JGH]
617 : ! **************************************************************************************************
618 10 : SUBROUTINE mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
619 : qs_kind_set, unit_nr, para_env)
620 :
621 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_smm
622 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list
623 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
624 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
625 : INTEGER, INTENT(IN) :: unit_nr
626 : TYPE(mp_para_env_type), POINTER :: para_env
627 :
628 : CHARACTER(len=2) :: element_symbol
629 : INTEGER :: ia, iab, iatom, ikind, iset, ishell, &
630 : ispin, l, lmax, lshell, m, ma, na, &
631 : natom, nspin
632 : LOGICAL :: found
633 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cmask, vec1, vec2
634 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weight
635 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, cmao
636 : TYPE(gto_basis_set_type), POINTER :: basis_set
637 :
638 : ! Analyze the MAO basis
639 10 : IF (unit_nr > 0) THEN
640 5 : WRITE (unit_nr, "(/,A)") " Analyze angular momentum character of MAOs "
641 : WRITE (unit_nr, "(T7,A,T15,A,T20,A,T40,A,T50,A,T60,A,T70,A,T80,A)") &
642 5 : "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G"
643 : END IF
644 10 : lmax = 4 ! analyze up to g-functions
645 10 : natom = SIZE(particle_set)
646 10 : nspin = SIZE(mao_coef)
647 58 : DO iatom = 1, natom
648 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
649 48 : element_symbol=element_symbol, kind_number=ikind)
650 48 : basis_set => mao_basis_set_list(ikind)%gto_basis_set
651 48 : CALL get_qs_kind(qs_kind_set(ikind), mao=na)
652 48 : CALL get_gto_basis_set(basis_set, nsgf=ma)
653 336 : ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na))
654 624 : weight = 0.0_dp
655 : CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, &
656 48 : block=block, found=found)
657 102 : DO ispin = 1, nspin
658 : CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, &
659 54 : block=cmao, found=found)
660 54 : IF (found) THEN
661 162 : DO l = 0, lmax
662 2445 : cmask = 0.0_dp
663 135 : iab = 0
664 675 : DO iset = 1, basis_set%nset
665 1635 : DO ishell = 1, basis_set%nshell(iset)
666 960 : lshell = basis_set%l(ishell, iset)
667 3810 : DO m = -lshell, lshell
668 2310 : iab = iab + 1
669 3270 : IF (l == lshell) cmask(iab) = 1.0_dp
670 : END DO
671 : END DO
672 : END DO
673 432 : DO ia = 1, na
674 6450 : vec1(1:ma) = cmask*cmao(1:ma, ia)
675 383430 : vec2(1:ma) = MATMUL(block, vec1)
676 6585 : weight(l, ia) = SUM(vec1(1:ma)*vec2(1:ma))
677 : END DO
678 : END DO
679 : END IF
680 54 : CALL para_env%sum(weight)
681 156 : IF (unit_nr > 0) THEN
682 81 : DO ia = 1, na
683 81 : IF (ispin == 1 .AND. ia == 1) THEN
684 : WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") &
685 24 : iatom, element_symbol, ispin, ia, weight(0:lmax, ia)
686 : ELSE
687 30 : WRITE (unit_nr, "(T17,i2,T20,i3,T31,5F10.4)") ispin, ia, weight(0:lmax, ia)
688 : END IF
689 : END DO
690 : END IF
691 : END DO
692 154 : DEALLOCATE (cmask, weight, vec1, vec2)
693 : END DO
694 20 : END SUBROUTINE mao_basis_analysis
695 :
696 : ! **************************************************************************************************
697 : !> \brief Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap
698 : !> \param matrix_q ...
699 : !> \param matrix_p ...
700 : !> \param matrix_s ...
701 : !> \param matrix_smm ...
702 : !> \param matrix_smo ...
703 : !> \param smm_list ...
704 : !> \param electra ...
705 : !> \param eps_filter ...
706 : !> \param nimages ...
707 : !> \param kpoints ...
708 : !> \param matrix_ks ...
709 : !> \param sab_orb ...
710 : !> \par History
711 : !> 08.2016 created [JGH]
712 : ! **************************************************************************************************
713 14 : SUBROUTINE mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, &
714 : electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
715 :
716 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q
717 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
718 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_smm, matrix_smo
719 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
720 : POINTER :: smm_list
721 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: electra
722 : REAL(KIND=dp), INTENT(IN) :: eps_filter
723 : INTEGER, INTENT(IN), OPTIONAL :: nimages
724 : TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
725 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
726 : POINTER :: matrix_ks
727 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
728 : OPTIONAL, POINTER :: sab_orb
729 :
730 : INTEGER :: im, ispin, nim, nocc, norb, nspin
731 14 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
732 : REAL(KIND=dp) :: elex, xkp(3)
733 : TYPE(dbcsr_type) :: ksmat, pmat, smat, tmat
734 :
735 14 : nim = 1
736 14 : IF (PRESENT(nimages)) nim = nimages
737 0 : IF (nim > 1) THEN
738 0 : CPASSERT(PRESENT(kpoints))
739 0 : CPASSERT(PRESENT(matrix_ks))
740 0 : CPASSERT(PRESENT(sab_orb))
741 : END IF
742 :
743 : ! Reference
744 14 : nspin = SIZE(matrix_p, 1)
745 30 : DO ispin = 1, nspin
746 16 : electra(ispin) = 0.0_dp
747 46 : DO im = 1, nim
748 16 : CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex)
749 32 : electra(ispin) = electra(ispin) + elex
750 : END DO
751 : END DO
752 :
753 : ! Q matrix
754 14 : NULLIFY (matrix_q)
755 14 : CALL dbcsr_allocate_matrix_set(matrix_q, nspin)
756 30 : DO ispin = 1, nspin
757 16 : ALLOCATE (matrix_q(ispin)%matrix)
758 16 : CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix)
759 30 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list)
760 : END DO
761 : ! temp matrix
762 14 : CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
763 : ! Q=APA(T)
764 30 : DO ispin = 1, nspin
765 30 : IF (nim == 1) THEN
766 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, &
767 16 : 0.0_dp, tmat, filter_eps=eps_filter)
768 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
769 16 : 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
770 : ELSE
771 : ! k-points
772 0 : CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix)
773 0 : CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix)
774 0 : CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix)
775 0 : CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb)
776 0 : CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb)
777 0 : CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb)
778 0 : NULLIFY (cell_to_index)
779 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
780 : ! calculate density matrix at gamma point
781 0 : xkp = 0.0_dp
782 : ! transform KS and S matrices to the gamma point
783 0 : CALL dbcsr_set(ksmat, 0.0_dp)
784 : CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, &
785 0 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
786 0 : CALL dbcsr_set(smat, 0.0_dp)
787 : CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, &
788 0 : xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
789 0 : norb = NINT(electra(ispin))
790 0 : nocc = MOD(2, nspin) + 1
791 0 : CALL calculate_p_gamma(pmat, ksmat, smat, kpoints, norb, REAL(nocc, KIND=dp))
792 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, pmat, &
793 0 : 0.0_dp, tmat, filter_eps=eps_filter)
794 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
795 0 : 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
796 0 : CALL dbcsr_release(pmat)
797 0 : CALL dbcsr_release(smat)
798 0 : CALL dbcsr_release(ksmat)
799 : END IF
800 : END DO
801 : ! free temp matrix
802 14 : CALL dbcsr_release(tmat)
803 :
804 14 : END SUBROUTINE mao_build_q
805 :
806 864 : END MODULE mao_methods
|
