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