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