Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate localized minimal basis
10 : !> \par History
11 : !> 12.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE minbas_methods
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_iterator_blocks_left, &
19 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
20 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks
22 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 : copy_fm_to_dbcsr,&
24 : dbcsr_allocate_matrix_set,&
25 : dbcsr_deallocate_matrix_set
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
27 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
28 : cp_fm_power
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_get_diag,&
34 : cp_fm_release,&
35 : cp_fm_to_fm_submat,&
36 : cp_fm_type
37 : USE kinds, ONLY: dp
38 : USE mao_basis, ONLY: mao_generate_basis
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE parallel_gemm_api, ONLY: parallel_gemm
41 : USE particle_methods, ONLY: get_particle_set
42 : USE particle_types, ONLY: particle_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kind_types, ONLY: qs_kind_type
46 : USE qs_ks_types, ONLY: get_ks_env,&
47 : qs_ks_env_type
48 : USE qs_mo_types, ONLY: get_mo_set,&
49 : mo_set_type
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_methods'
56 :
57 : PUBLIC :: minbas_calculation
58 :
59 : ! **************************************************************************************************
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief ...
65 : !> \param qs_env ...
66 : !> \param mos ...
67 : !> \param quambo ...
68 : !> \param mao ...
69 : !> \param iounit ...
70 : !> \param full_ortho ...
71 : !> \param eps_filter ...
72 : ! **************************************************************************************************
73 0 : SUBROUTINE minbas_calculation(qs_env, mos, quambo, mao, iounit, full_ortho, eps_filter)
74 : TYPE(qs_environment_type), POINTER :: qs_env
75 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quambo
77 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
78 : POINTER :: mao
79 : INTEGER, INTENT(IN), OPTIONAL :: iounit
80 : LOGICAL, INTENT(IN), OPTIONAL :: full_ortho
81 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
82 :
83 : CHARACTER(len=*), PARAMETER :: routineN = 'minbas_calculation'
84 :
85 : INTEGER :: handle, homo, i, iab, ispin, nao, natom, &
86 : ndep, nmao, nmo, nmx, np, np1, nspin, &
87 : nvirt, unit_nr
88 0 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
89 : LOGICAL :: do_minbas, my_full_ortho
90 : REAL(KIND=dp) :: my_eps_filter
91 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dval, dvalo, dvalv, eigval
92 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
93 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c, &
94 : fm_struct_d, fm_struct_e
95 : TYPE(cp_fm_type) :: fm1, fm2, fm3, fm4, fm5, fm6, fma, fmb, &
96 : fmwork
97 : TYPE(cp_fm_type), POINTER :: fm_mos
98 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
99 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
100 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
101 : TYPE(dbcsr_type) :: smao, sortho
102 : TYPE(dbcsr_type), POINTER :: smat
103 : TYPE(dft_control_type), POINTER :: dft_control
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107 : TYPE(qs_ks_env_type), POINTER :: ks_env
108 :
109 0 : CALL timeset(routineN, handle)
110 :
111 0 : IF (PRESENT(iounit)) THEN
112 0 : unit_nr = iounit
113 : ELSE
114 : unit_nr = -1
115 : END IF
116 :
117 0 : IF (PRESENT(full_ortho)) THEN
118 0 : my_full_ortho = full_ortho
119 : ELSE
120 : my_full_ortho = .FALSE.
121 : END IF
122 :
123 0 : IF (PRESENT(eps_filter)) THEN
124 0 : my_eps_filter = eps_filter
125 : ELSE
126 0 : my_eps_filter = 1.0e-10_dp
127 : END IF
128 :
129 0 : CALL get_qs_env(qs_env, dft_control=dft_control)
130 0 : nspin = dft_control%nspins
131 :
132 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
133 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
134 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
135 0 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
136 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
137 0 : CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
138 0 : nmao = SUM(col_blk_sizes)
139 : ! check if MAOs have been specified
140 0 : DO iab = 1, natom
141 0 : IF (col_blk_sizes(iab) < 0) THEN
142 0 : CPABORT("Number of MAOs has to be specified in KIND section for all elements")
143 : END IF
144 : END DO
145 0 : CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=nmo)
146 :
147 0 : IF (unit_nr > 0) THEN
148 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Atomic Basis Set Functions :', nao
149 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Minimal Basis Set Functions :', nmao
150 0 : IF (nspin == 1) THEN
151 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Molecular Orbitals available :', nmo
152 : ELSE
153 0 : DO ispin = 1, nspin
154 0 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
155 : WRITE (unit_nr, '(T2,A,i2,T71,I10)') &
156 0 : 'Total Number of Molecular Orbitals available for Spin ', ispin, nmx
157 : END DO
158 : END IF
159 : END IF
160 0 : CPASSERT(nmao <= nao)
161 0 : DO ispin = 1, nspin
162 0 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
163 0 : IF (nmx /= nmo) EXIT
164 : END DO
165 0 : do_minbas = .TRUE.
166 0 : IF (nmao > nmo) THEN
167 0 : IF (unit_nr > 0) THEN
168 0 : WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
169 : END IF
170 : do_minbas = .FALSE.
171 0 : ELSE IF (nmo /= nmx) THEN
172 0 : IF (unit_nr > 0) THEN
173 0 : WRITE (unit_nr, '(T2,A)') 'Different Number of Alpha and Beta MOs'
174 0 : WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
175 : END IF
176 : do_minbas = .FALSE.
177 : ELSE
178 0 : IF (nao > nmo) THEN
179 0 : IF (unit_nr > 0) THEN
180 0 : WRITE (unit_nr, '(T2,A)') 'WARNING: Only a subset of MOs is available: Analysis depends on MOs'
181 : END IF
182 : END IF
183 : END IF
184 :
185 0 : IF (do_minbas) THEN
186 : ! initialize QUAMBOs
187 0 : NULLIFY (quambo)
188 0 : CALL dbcsr_allocate_matrix_set(quambo, nspin)
189 0 : DO ispin = 1, nspin
190 : ! coeficients
191 0 : ALLOCATE (quambo(ispin)%matrix)
192 : CALL dbcsr_create(matrix=quambo(ispin)%matrix, &
193 : name="QUAMBO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
194 0 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
195 : END DO
196 :
197 : ! initialize MAOs
198 : ! optimize MAOs (mao_coef is allocated in the routine)
199 0 : CALL mao_generate_basis(qs_env, mao_coef)
200 :
201 : ! sortho (nmao x nmao)
202 : CALL dbcsr_create(sortho, name="SORTHO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
203 0 : row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes)
204 0 : CALL dbcsr_reserve_diag_blocks(matrix=sortho)
205 :
206 0 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
207 :
208 : ! temporary FM matrices
209 0 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
210 0 : NULLIFY (fm_struct_a, fm_struct_b)
211 : CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, &
212 0 : para_env=para_env, context=blacs_env)
213 : CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmo, ncol_global=nmao, &
214 0 : para_env=para_env, context=blacs_env)
215 0 : CALL cp_fm_create(fm1, fm_struct_a)
216 0 : CALL cp_fm_create(fm2, fm_struct_b)
217 0 : CALL cp_fm_create(fma, fm_struct_b)
218 0 : CALL cp_fm_create(fmb, fm_struct_b)
219 :
220 0 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
221 0 : smat => matrix_s(1, 1)%matrix
222 0 : DO ispin = 1, nspin
223 :
224 : ! SMAO = Overlap*MAOs
225 0 : CALL dbcsr_create(smao, name="S*MAO", template=mao_coef(1)%matrix)
226 0 : CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef(ispin)%matrix, 0.0_dp, smao)
227 : ! a(nj)* = C(vn)(T) * SMAO(vj)
228 0 : CALL copy_dbcsr_to_fm(smao, fm1)
229 0 : CALL get_mo_set(mos(ispin), mo_coeff=fm_mos)
230 0 : CALL parallel_gemm("T", "N", nmo, nmao, nao, 1.0_dp, fm_mos, fm1, 0.0_dp, fm2)
231 0 : CALL dbcsr_release(smao)
232 0 : CALL get_mo_set(mo_set=mos(ispin), homo=homo)
233 0 : IF (unit_nr > 0) THEN
234 0 : WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Occupied Valence Set', 'Spin ', ispin, homo
235 : END IF
236 0 : nvirt = nmo - homo
237 0 : NULLIFY (fm_struct_c)
238 : CALL cp_fm_struct_create(fm_struct_c, nrow_global=nvirt, ncol_global=nvirt, &
239 0 : para_env=para_env, context=blacs_env)
240 0 : CALL cp_fm_create(fm3, fm_struct_c)
241 0 : CALL cp_fm_create(fm4, fm_struct_c)
242 : ! B(vw) = a(vj)* * a(wj)*
243 : CALL parallel_gemm("N", "T", nvirt, nvirt, nmao, 1.0_dp, fm2, fm2, 0.0_dp, fm3, &
244 0 : a_first_row=homo + 1, b_first_row=homo + 1)
245 0 : ALLOCATE (eigval(nvirt))
246 0 : CALL choose_eigv_solver(fm3, fm4, eigval)
247 : ! SVD(B) -> select p largest eigenvalues and vectors
248 0 : np = nmao - homo
249 0 : np1 = nvirt - np + 1
250 0 : IF (unit_nr > 0) THEN
251 0 : WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Virtual Valence Set', 'Spin ', ispin, np
252 : END IF
253 : ! R(vw) = SUM_p T(vp)*T(wp)
254 : CALL parallel_gemm("N", "T", nvirt, nvirt, np, 1.0_dp, fm4, fm4, 0.0_dp, fm3, &
255 0 : a_first_col=np1, b_first_col=np1)
256 : !
257 0 : ALLOCATE (dval(nmao), dvalo(nmao), dvalv(nmao))
258 0 : NULLIFY (fm_struct_d)
259 : CALL cp_fm_struct_create(fm_struct_d, nrow_global=nvirt, ncol_global=nmao, &
260 0 : para_env=para_env, context=blacs_env)
261 0 : CALL cp_fm_create(fm5, fm_struct_d)
262 0 : NULLIFY (fm_struct_e)
263 : CALL cp_fm_struct_create(fm_struct_e, nrow_global=nmao, ncol_global=nmao, &
264 0 : para_env=para_env, context=blacs_env)
265 0 : CALL cp_fm_create(fm6, fm_struct_e)
266 : ! D(j) = SUM_n (a(nj)*)^2 + SUM_vw R(vw) * a(vj)* * a(wj)*
267 : CALL parallel_gemm("N", "N", nvirt, nmao, nvirt, 1.0_dp, fm3, fm2, 0.0_dp, fm5, &
268 0 : b_first_row=homo + 1)
269 : CALL parallel_gemm("T", "N", nmao, nmao, nvirt, 1.0_dp, fm2, fm5, 0.0_dp, fm6, &
270 0 : a_first_row=homo + 1)
271 0 : CALL cp_fm_get_diag(fm6, dvalv(1:nmao))
272 0 : CALL parallel_gemm("T", "N", nmao, nmao, homo, 1.0_dp, fm2, fm2, 0.0_dp, fm6)
273 0 : CALL cp_fm_get_diag(fm6, dvalo(1:nmao))
274 0 : DO i = 1, nmao
275 0 : dval(i) = 1.0_dp/SQRT(dvalo(i) + dvalv(i))
276 : END DO
277 : ! scale intermediate expansion
278 0 : CALL cp_fm_to_fm_submat(fm2, fma, homo, nmao, 1, 1, 1, 1)
279 0 : CALL cp_fm_to_fm_submat(fm5, fma, nvirt, nmao, 1, 1, homo + 1, 1)
280 0 : CALL cp_fm_column_scale(fma, dval)
281 : ! Orthogonalization
282 0 : CALL cp_fm_create(fmwork, fm_struct_e)
283 0 : CALL parallel_gemm("T", "N", nmao, nmao, nmo, 1.0_dp, fma, fma, 0.0_dp, fm6)
284 0 : IF (my_full_ortho) THEN
285 : ! full orthogonalization
286 0 : CALL cp_fm_power(fm6, fmwork, -0.5_dp, 1.0e-12_dp, ndep)
287 0 : IF (ndep > 0 .AND. unit_nr > 0) THEN
288 0 : WRITE (unit_nr, '(T2,A,T71,I10)') 'Warning: linear dependent basis ', ndep
289 : END IF
290 0 : CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
291 : ELSE
292 : ! orthogonalize on-atom blocks
293 0 : CALL copy_fm_to_dbcsr(fm6, sortho, keep_sparsity=.TRUE.)
294 0 : CALL diag_sqrt_invert(sortho)
295 0 : CALL copy_dbcsr_to_fm(sortho, fm6)
296 0 : CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
297 : END IF
298 : ! store as QUAMBO
299 0 : CALL parallel_gemm("N", "N", nao, nmao, nmo, 1.0_dp, fm_mos, fmb, 0.0_dp, fm1)
300 0 : CALL copy_fm_to_dbcsr(fm1, quambo(ispin)%matrix)
301 0 : CALL dbcsr_filter(quambo(ispin)%matrix, my_eps_filter)
302 : !
303 0 : DEALLOCATE (eigval, dval, dvalo, dvalv)
304 0 : CALL cp_fm_release(fm3)
305 0 : CALL cp_fm_release(fm4)
306 0 : CALL cp_fm_release(fm5)
307 0 : CALL cp_fm_release(fm6)
308 0 : CALL cp_fm_release(fmwork)
309 0 : CALL cp_fm_struct_release(fm_struct_c)
310 0 : CALL cp_fm_struct_release(fm_struct_d)
311 0 : CALL cp_fm_struct_release(fm_struct_e)
312 :
313 : END DO
314 :
315 : ! clean up
316 0 : CALL cp_fm_release(fm1)
317 0 : CALL cp_fm_release(fm2)
318 0 : CALL cp_fm_release(fma)
319 0 : CALL cp_fm_release(fmb)
320 0 : CALL cp_fm_struct_release(fm_struct_a)
321 0 : CALL cp_fm_struct_release(fm_struct_b)
322 0 : CALL dbcsr_release(sortho)
323 :
324 : ! return MAOs if requested
325 0 : IF (PRESENT(mao)) THEN
326 0 : mao => mao_coef
327 : ELSE
328 0 : CALL dbcsr_deallocate_matrix_set(mao_coef)
329 : END IF
330 :
331 : ELSE
332 0 : NULLIFY (quambo)
333 : END IF
334 :
335 0 : CALL timestop(handle)
336 :
337 0 : END SUBROUTINE minbas_calculation
338 :
339 : ! **************************************************************************************************
340 : !> \brief ...
341 : !> \param sortho ...
342 : ! **************************************************************************************************
343 0 : SUBROUTINE diag_sqrt_invert(sortho)
344 : TYPE(dbcsr_type) :: sortho
345 :
346 : INTEGER :: i, iatom, info, jatom, lwork, n
347 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
348 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat
349 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock
350 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
351 :
352 0 : CALL dbcsr_iterator_start(dbcsr_iter, sortho)
353 0 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
354 0 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock)
355 0 : CPASSERT(iatom == jatom)
356 0 : n = SIZE(sblock, 1)
357 0 : lwork = MAX(n*n, 100)
358 0 : ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
359 0 : amat(1:n, 1:n) = sblock(1:n, 1:n)
360 0 : info = 0
361 0 : CALL dsyev("V", "U", n, amat, n, w, work, lwork, info)
362 0 : CPASSERT(info == 0)
363 0 : w(1:n) = 1._dp/SQRT(w(1:n))
364 0 : DO i = 1, n
365 0 : bmat(1:n, i) = amat(1:n, i)*w(i)
366 : END DO
367 0 : sblock(1:n, 1:n) = MATMUL(amat, TRANSPOSE(bmat))
368 0 : DEALLOCATE (amat, bmat, w, work)
369 : END DO
370 0 : CALL dbcsr_iterator_stop(dbcsr_iter)
371 :
372 0 : END SUBROUTINE diag_sqrt_invert
373 :
374 : END MODULE minbas_methods
|
